C++ TPZAutoPointer::NPoints方法代碼示例

REAL TPZMGAnalysis::ElementError(TPZInterpolatedElement *fine, TPZInterpolatedElement *coarse, TPZTransform &tr,
								 void (*f)(TPZVec<REAL> &loc, TPZVec<REAL> &val, TPZFMatrix<REAL> &deriv),REAL &truerror){
	// accumulates the transfer coefficients between the current element and the
	// coarse element into the transfer matrix, using the transformation t
	int locnod = fine->NConnects();
	int cornod = coarse->NConnects();
	int locmatsize = fine->NShapeF();
	int cormatsize = coarse->NShapeF();
	REAL error = 0.;
	truerror = 0.;
	
	REAL loclocmatstore[500] = {0.},loccormatstore[500] = {0.};
	TPZFMatrix<REAL> loclocmat(locmatsize,locmatsize,loclocmatstore,500);
	TPZFMatrix<REAL> loccormat(locmatsize,cormatsize,loccormatstore,500);
	
	TPZAutoPointer<TPZIntPoints> intrule = fine->GetIntegrationRule().Clone();
	int dimension = fine->Dimension();
	int numdof = fine->Material()->NStateVariables();
	TPZBlock<REAL> &locblock = fine->Mesh()->Block();
	TPZFMatrix<REAL> &locsolmesh = fine->Mesh()->Solution();
	
	TPZBlock<REAL> &corblock = coarse->Mesh()->Block();
	TPZFMatrix<REAL> &corsolmesh = coarse->Mesh()->Solution();
	
	TPZVec<REAL> locsol(numdof);
	TPZFMatrix<REAL> locdsol(dimension,numdof);
	
	TPZVec<REAL> corsol(numdof);
	TPZFMatrix<REAL> cordsol(dimension,numdof);
	
	TPZManVector<int> prevorder(dimension),
    order(dimension);
	intrule->GetOrder(prevorder);
	
	TPZManVector<int> interpolation(dimension);
	fine->GetInterpolationOrder(interpolation);
	
	// compute the interpolation order of the shapefunctions squared
	int dim;
	int maxorder = interpolation[0];
	for(dim=0; dim<interpolation.NElements(); dim++) {
		maxorder = interpolation[dim] < maxorder ? maxorder : interpolation[dim];
	}
	for(dim=0; dim<dimension; dim++) {
		order[dim] = 20;
	}
	intrule->SetOrder(order);
	
	
	REAL locphistore[50]={0.},locdphistore[150]={0.};
	TPZFMatrix<REAL> locphi(locmatsize,1,locphistore,50);
	TPZFMatrix<REAL> locdphi(dimension,locmatsize,locdphistore,150),locdphix(dimension,locmatsize);
	// derivative of the shape function
	// in the master domain
	
	REAL corphistore[50]={0.},cordphistore[150]={0.};
	TPZFMatrix<REAL> corphi(cormatsize,1,corphistore,50);
	TPZFMatrix<REAL> cordphi(dimension,cormatsize,cordphistore,150), cordphix(dimension,cormatsize);
	// derivative of the shape function
	// in the master domain
	
	REAL jacobianstore[9],
    axesstore[9];
	TPZManVector<REAL> int_point(dimension),
    coarse_int_point(dimension);
	TPZFMatrix<REAL> jacfine(dimension,dimension,jacobianstore,9),jacinvfine(dimension,dimension);
	TPZFMatrix<REAL> axesfine(3,3,axesstore,9);
	TPZManVector<REAL> xfine(3);
	TPZFMatrix<REAL> jaccoarse(dimension,dimension),jacinvcoarse(dimension,dimension);
	TPZFMatrix<REAL> axescoarse(3,3), axesinner(dimension,dimension);
	TPZManVector<REAL> xcoarse(3);
	
	REAL jacdetcoarse;
	int numintpoints = intrule->NPoints();
	REAL weight;
	int i,j,k;
	
	TPZVec<REAL> truesol(numdof);
	TPZFMatrix<REAL> truedsol(dimension,numdof);
	for(int int_ind = 0; int_ind < numintpoints; ++int_ind) {
		
		intrule->Point(int_ind,int_point,weight);
		REAL jacdetfine;
		fine->Reference()->Jacobian( int_point, jacfine , axesfine, jacdetfine, jacinvfine);
		fine->Reference()->X(int_point, xfine);
		if(f) f(xfine,truesol,truedsol);
		fine->Shape(int_point,locphi,locdphi);
		tr.Apply(int_point,coarse_int_point);
		coarse->Shape(coarse_int_point,corphi,cordphi);
		coarse->Reference()->Jacobian( coarse_int_point,jaccoarse,axescoarse, jacdetcoarse, jacinvcoarse);
		coarse->Reference()->X(coarse_int_point,xcoarse);
		REAL dist = (xfine[0]-xcoarse[0])*(xfine[0]-xcoarse[0])+(xfine[1]-xcoarse[1])*(xfine[1]-xcoarse[1])+(xfine[2]-xcoarse[2])*(xfine[2]-xcoarse[2]);
		if(dist > 1.e-6) cout << "TPZMGAnalysis::ElementError transformation between fine and coarse is wrong\n";
		for(i=0; i<dimension; i++) {
			for(j=0; j<dimension; j++) {
				axesinner(i,j) = 0.;
				for(k=0; k<3; k++) axesinner(i,j) += axesfine(i,k)*axescoarse(j,k);
			}
		}
		if(fabs(axesinner(0,0)-1.) > 1.e-6 || fabs(axesinner(1,1)-1.) > 1.e-6 || fabs(axesinner(0,1)) > 1.e-6 || fabs(axesinner(1,0)) > 1.e-6) {
//.........這裏部分代碼省略.........

REAL TPZAdaptMesh::UseTrueError(TPZInterpolatedElement *coarse, 
                                void (*f)(const TPZVec<REAL> &loc, TPZVec<REAL> &val, TPZFMatrix<REAL> &deriv)){
    if (coarse->Material()->Id() < 0) return 0.0;
    
    REAL error = 0.;
    
    //  REAL loclocmatstore[500] = {0.},loccormatstore[500] = {0.};
    // TPZFMatrix loccormat(locmatsize,cormatsize,loccormatstore,500);
    
    //Cesar 25/06/03 - Uso a ordem m�xima???
    TPZAutoPointer<TPZIntPoints> intrule = coarse->GetIntegrationRule().Clone();
    
    int dimension = coarse->Dimension();
    int numdof = coarse->Material()->NStateVariables();
        
    //TPZSolVec corsol;
    //TPZGradSolVec cordsol;
    TPZGradSolVec cordsolxy;
//    TPZVec<REAL> corsol(numdof);
//    TPZFNMatrix<9,REAL> cordsol(dimension,numdof),cordsolxy(dimension,numdof);
    
    TPZManVector<int> order(dimension,20);
    intrule->SetOrder(order);
    
    // derivative of the shape function
    // in the master domain
    TPZManVector<REAL,3> coarse_int_point(dimension);
//    TPZFNMatrix<9,REAL> jaccoarse(dimension,dimension),jacinvcoarse(dimension,dimension);
//    TPZFNMatrix<9,REAL> axescoarse(3,3);
//    TPZManVector<REAL,3> xcoarse(3);
    TPZFNMatrix<9,REAL> axesinner(3,3);
    
    
    TPZMaterialData datacoarse;
    coarse->InitMaterialData(datacoarse);
    
//    REAL jacdetcoarse;
    int numintpoints = intrule->NPoints();
    REAL weight;
    
    TPZVec<REAL> truesol(numdof);
    TPZFMatrix<REAL> truedsol(dimension,numdof);
    for(int int_ind = 0; int_ind < numintpoints; ++int_ind) {
        intrule->Point(int_ind,coarse_int_point,weight);
        //coarse->Reference()->X(coarse_int_point, xcoarse);
        coarse->Reference()->X(coarse_int_point, datacoarse.x);
        //if(f) f(xcoarse,truesol,truedsol);
        if(f) f(datacoarse.x,truesol,truedsol);

//        coarse->Reference()->Jacobian(coarse_int_point, jaccoarse, axescoarse, jacdetcoarse, jacinvcoarse);
        coarse->Reference()->Jacobian(coarse_int_point, datacoarse.jacobian, datacoarse.axes, datacoarse.detjac, datacoarse.jacinv);
        //weight *= fabs(jacdetcoarse);
        weight *= fabs(datacoarse.detjac);
//Er        int iv=0;
//        corsol[0].Fill(0.);
//        cordsol[0].Zero();

        //coarse->ComputeSolution(coarse_int_point, corsol, cordsol, axescoarse);
        coarse->ComputeShape(coarse_int_point,datacoarse);
        coarse->ComputeSolution(coarse_int_point,datacoarse);
        
        //int nc = cordsol[0].Cols();
        int nc = datacoarse.dsol[0].Cols();
        for (int col=0; col<nc; col++)
        {
            for (int d=0; d<dimension; d++) {
                REAL deriv = 0.;
                for (int d2=0; d2<dimension; d2++) {
                    deriv += datacoarse.dsol[0](d2,col)*datacoarse.axes(d2,d);
                }
               // cordsolxy[0](d,col) = deriv;
            }
        }
        int jn;
        for(jn=0; jn<numdof; jn++) {
            for(int d=0; d<dimension; d++) {
                error += (datacoarse.dsol[0](d,jn)-truedsol(d,jn))*(datacoarse.dsol[0](d,jn)-truedsol(d,jn))*weight;
            }
        }
    }
    return error;
}