C++ TPZStack::NElements方法代码示例

void TPZAnalysisError::ExpandConnected(TPZStack<TPZCompElSide> &singel){
	int64_t nelem = singel.NElements();
	TPZStack<TPZGeoElSide> gelstack;
	TPZStack<TPZCompElSide> celstack;
	int64_t iel;
	for(iel=0; iel<nelem; iel++) {
		TPZCompElSide celside = singel[iel];
		TPZGeoElSide gelside;
		gelside = celside.Reference();
		if(!gelside.Exists()) continue;
		gelstack.Resize(0);
		cout << "This part needs to be fixed\n";
		//		gelside.Element()->LowerDimensionSides(gelside.Side(),gelstack);
		while(gelstack.NElements()) {
			TPZGeoElSide gelsideloc;
			gelsideloc = gelstack.Pop();
			gelsideloc.EqualLevelCompElementList(celstack,1,0);
			while(celstack.NElements()) {
				TPZCompElSide celsideloc = celstack.Pop();
				if(! celsideloc.Exists()) continue;
				int64_t nelsing = singel.NElements();
				int64_t smel;
				for(smel=0; smel<nelsing; smel++) if(singel[smel].Element() == celsideloc.Element()) break;
				if(smel != nelsing) singel.Push(celsideloc);
			}
		}
	}
}

void TPZMultiphysicsElement::CreateInterfaces()
{
    //nao verifica-se caso o elemento de contorno
    //eh maior em tamanho que o interface associado
    //caso AdjustBoundaryElement nao for aplicado
    //a malha eh criada consistentemente
    TPZGeoEl *ref = Reference();
    int nsides = ref->NSides();
    int InterfaceDimension = Mesh()->Dimension() - 1;
    int side;
    nsides--;//last face
    for(side=nsides;side>=0;side--) {
        if(ref->SideDimension(side) != InterfaceDimension) continue;
        TPZCompElSide thisside(this,side);
        if(this->ExistsInterface(side)) {
            //      cout << "TPZCompElDisc::CreateInterface inconsistent: interface already exists\n";
            continue;
        }
        TPZStack<TPZCompElSide> highlist;
        thisside.HigherLevelElementList(highlist,0,1);
        //a interface se cria uma vez so quando existem ambos
        //elementos esquerdo e direito (compu tacionais)
        if(!highlist.NElements()) {
            this->CreateInterface(side);//s�tem iguais ou grande => pode criar a interface
        } else {
            int64_t ns = highlist.NElements();
            int64_t is;
            for(is=0; is<ns; is++) {//existem pequenos ligados ao lado atual
                const int higheldim = highlist[is].Reference().Dimension();
                if(higheldim != InterfaceDimension) continue;
                // 	TPZCompElDisc *del = dynamic_cast<TPZCompElDisc *> (highlist[is].Element());
                // 	if(!del) continue;
                
                TPZCompEl *del = highlist[is].Element();
                if(!del) continue;
                
                TPZCompElSide delside( del, highlist[is].Side() );
                TPZMultiphysicsElement * delSp = dynamic_cast<TPZMultiphysicsElement *>(del);
                if (!delSp){
                    PZError << "\nERROR AT " << __PRETTY_FUNCTION__ <<  " - CASE NOT AVAILABLE\n";
                    return;
                }
                if ( delSp->ExistsInterface(highlist[is].Side()) ) {
                    //          cout << "TPZCompElDisc::CreateInterface inconsistent: interface already exists\n";
                }
                else {
                    delSp->CreateInterface(highlist[is].Side());
                }
            }
        }
    }
}

int TPZErrorIndicator::GetRefSide(TPZCompEl *cel, int sidedim, int sidestate, TPZMatrix *errormat){
  int e,s,nsides = cel->Reference()->NSides();
  REAL sum = 0.;
  int side = -1;
  for (s=0;s<nsides;s++){
    TPZCompElSide celside (cel,s);
    if (celside.Reference().Dimension() != sidedim) continue;
    TPZStack<TPZCompElSide> elsidevec;
    celside.EqualLevelElementList(elsidevec,0,0);
    REAL auxsum = 0;
    int neigbyside = elsidevec.NElements();
    for (e=0;e<neigbyside;e++){
      int index = elsidevec[e].Element()->Index();
      double val = 0.;
      val = errormat->Get(e,sidestate);
      auxsum += val;
    }
    auxsum /= ((REAL) neigbyside);
    if(auxsum > sum){
      sum = auxsum;
      side = s;
    }
  }
  return side;
}

int TPZCheckGeom::CheckRefinement(TPZGeoEl *gel){
	
	int check = 0;
	if(!gel || !gel->HasSubElement()) return check;
	int nsides = gel->NSides();
	int is;
	for(is=0; is<nsides; is++) {
		TPZStack<TPZGeoElSide> subel;
		gel->GetSubElements2(is,subel);
		int nsub = subel.NElements();
		int isub;
		for(isub=0; isub<nsub; isub++) {
			TPZGeoElSide fath = subel[isub].Father2();
			int son = subel[isub].Element()->WhichSubel();
			if(fath.Side() != is) {
				PZError << "TPZCheckGeom::CheckRefinement non corresponding subelement/sides son "
				<< son << " sonside " << subel[isub].Side() << " fathside " << is <<
				" fath2side " << fath.Side() << endl;
				gel->Print();
				check = 1;
			}
		}
	}
	int nsub = gel->NSubElements();
	for(is=0; is<nsub; is++) {
		TPZGeoEl *sub = gel->SubElement(is);
		int nsubsides = sub->NSides();
		int iss;
		for(iss=0; iss<nsubsides; iss++) {
			check = (CheckSubFatherTransform(sub,iss) || check);
		}
	}
	
	return check;
}

void TPZGeoCloneMesh::SetElements(TPZStack <TPZGeoEl *> &patch, TPZGeoEl *ref){
    
    int64_t i;
    if (!ref){
        cout << "TPZGeoCloneMesh::Error\n Reference element must not be null\n";
        DebugStop();
    }
    //fGeoRoot = ref;
    CloneElement(ref);
    fGeoRoot = fMapElements[ref];
    
    //  cout << "\n\n\nTeste\n\n\n";
    //  Print(cout);
    
    int64_t nel = patch.NElements();
    //  fReferenceElement.Resize(nel);
    for (i=0; i<nel; i++){
        if(patch[i]) {
            TPZGeoEl *gel = patch[i];
            TPZGeoEl *father = gel->Father();
            while(father) {
                gel = father;
                father = gel->Father();
            }
            //      cout << "\nElemento a ser clonado:\n";
            //      gel->Print(cout);
            CloneElement(gel);
            // verificar se neighbour.Element ja esta no map --->>>> j� � feito no CloneElement
            TPZGeoEl *localpatch = fMapElements[patch[i]];
#ifdef PZDEBUG 
			if (localpatch == 0) {
				DebugStop();
			}
#endif
            fPatchElements.insert(localpatch);
            AddBoundaryConditionElements(patch[i]);

//			cout << "Printing fPatchElements("<< i << "_NEl_" << fPatchElements.size() << ") : " << endl;
//			std::set<TPZGeoEl *>::iterator it;
//			for(it=fPatchElements.begin();it!=fPatchElements.end();it++) {
//				(*it)->Print(cout);
//			}
		}
    }
}

void TPZFrontNonSym::Compress(){
	//	PrintGlobal("Before COmpress");
	//	Print("Before Compress", cout);
	TPZStack <int> from;
	int nfound;
	int i, j;
	for(i = 0; i < fFront; i++){
		if(fGlobal[i] != -1) from.Push(i);
	}
	
	/**
	 *First fLocal initialization
	 *Any needed updates is done on next loop
	 */
	nfound = from.NElements();
	for(i=0;i<nfound;i++) {
		fGlobal[i]=fGlobal[from[i]];
		//fGlobal[from[i]] = -1;
		fLocal[fGlobal[i]] = i;
	}
	for(;i<fGlobal.NElements();i++) fGlobal[i] = -1;
	
	if(nfound+fFree.NElements()!=fFront) cout << "TPZFront.Compress inconsistent data structure\n";
	
	fFront = nfound;
	fFree.Resize(0);	
	fGlobal.Resize(fFront);
	if(fData.NElements()==0) return;
	
	for(j = 0; j < nfound; j++){
		for(i = 0; i < nfound; i++){
			Element(i,j) = Element(from[i], from[j]);
			if(from[i]!=i || from[j]!=j) Element(from[i],from[j])=0.;
		}
		//		fGlobal[i] = fGlobal[from[i]];
		//		fLocal[fGlobal[i]] = i;
	}
	
	//	Print("After Compress", cout);
	//	PrintGlobal("After Compress",output);
}

void TPZGenPartialGrid::SetBC(TPZGeoMesh &g, int side, int bc) {
	if(side == 0 && fRangey[0] != 0) return;
	if(side == 1 && fRangex[1] != fNx[0]) return;
	if(side == 2 && fRangey[1] != fNx[1]) return;
	if(side == 3 && fRangex[0] != 0) return;
	TPZStack<TPZGeoEl*> ElementVec;
	TPZStack<int> Sides;
	TPZVec<int64_t> cornernodes(4);
	cornernodes[0] = NodeIndex(fRangex[0],fRangey[0]);
	cornernodes[1] = NodeIndex(fRangex[1],fRangey[0]);
	cornernodes[2] = NodeIndex(fRangex[1],fRangey[1]);
	cornernodes[3] = NodeIndex(fRangex[0],fRangey[1]);
	g.GetBoundaryElements(cornernodes[side],cornernodes[(side+1)%4],ElementVec, Sides);
	int64_t numel = ElementVec.NElements();
	for(int64_t el=0; el<numel; el++) {
		TPZGeoEl *gel = (TPZGeoEl *) ElementVec[el];
		if(gel) {
			TPZGeoElBC(gel,Sides[el],bc);
		}
	}
}

int TPZCheckGeom::CheckElement(TPZGeoEl *gel) {
	
	int check = 0;
	int nsides = gel->NSides();
	int geldim = gel->Dimension();
	int is;
	for(is=nsides-1; is>=0; is--) {
		TPZStack<TPZGeoElSide> highdim;
		int dim;
		int sidedim = gel->SideDimension(is);
		for(dim = sidedim+1; dim<= geldim; dim++) {
			gel->AllHigherDimensionSides(is,dim,highdim);
		}
		int nhighdim = highdim.NElements();
		int idim;
		for(idim=0; idim<nhighdim; idim++) {
			check = (CheckSideTransform(gel,is,highdim[idim].Side()) || check);
		}
	}
	return check;
}

TPZMultiphysicsInterfaceElement * TPZMultiphysicsElement::CreateInterface(int side)
{
	//  LOGPZ_INFO(logger, "Entering CreateInterface");
	TPZMultiphysicsInterfaceElement * newcreatedinterface = NULL;
	
	TPZGeoEl *ref = Reference();
	if(!ref) {
		LOGPZ_WARN(logger, "Exiting CreateInterface Null reference reached - NULL interface returned");
		return newcreatedinterface;
	}
	
	TPZCompElSide thisside(this,side);
	TPZStack<TPZCompElSide> list;
	list.Resize(0);
	thisside.EqualLevelElementList(list,0,0);//retorna distinto ao atual ou nulo
	int64_t size = list.NElements();
	//espera-se ter os elementos computacionais esquerdo e direito
	//ja criados antes de criar o elemento interface
    // try to create an interface element between myself and an equal sized neighbour
    for (int64_t is=0; is<size; is++)
    {
		//Interface has the same material of the neighbour with lesser dimension.
		//It makes the interface have the same material of boundary conditions (TPZCompElDisc with interface dimension)
		int matid;
		int thisdim = this->Dimension();
		int neighbourdim = list[is].Element()->Dimension();
        
        if(thisdim != neighbourdim)
        {
            if (thisdim < neighbourdim)
            {
                // return the material id of boundary condition IF the associated material is derived from bndcond
                TPZMaterial *mat = this->Material();
                TPZBndCond *bnd = dynamic_cast<TPZBndCond *>(mat);
                if(bnd)
                {
                    matid = this->Material()->Id();
                }
                else
                {
                    matid = this->Mesh()->Reference()->InterfaceMaterial(this->Reference()->MaterialId(),list[0].Element()->Reference()->MaterialId());
                    continue;
                }
            }
			else
            {
                TPZMaterial *mat = list[is].Element()->Material();
                TPZBndCond *bnd = dynamic_cast<TPZBndCond *>(mat);
                if (bnd) {
                    matid = bnd->Id();
                }
                else {
                    matid = this->Mesh()->Reference()->InterfaceMaterial(this->Reference()->MaterialId(), list[0].Element()->Reference()->MaterialId());
                    continue;
                }
            }
        }else
        {
            matid = this->Mesh()->Reference()->InterfaceMaterial(this->Reference()->MaterialId(), list[0].Element()->Reference()->MaterialId());
            if(matid == GMESHNOMATERIAL)
            {
                continue;
            }
        }
		
		int64_t index;
		
		
		TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid); //isto acertou as vizinhanas da interface geometrica com o atual
		if(!gel){
			DebugStop();
#ifdef LOG4CXX
            if (logger->isDebugEnabled())
			{
				std::stringstream sout;
				sout << "CreateBCGeoEl devolveu [email protected]@@@";
				LOGPZ_DEBUG(logger,sout.str());
			}
#endif
		}
        bool withmem = fMesh->ApproxSpace().NeedsMemory();
		if(Dimension() > list[is].Reference().Dimension()) {
			//o de volume eh o direito caso um deles seja BC
			//a normal aponta para fora do contorno
			TPZCompElSide thiscompelside(this, thisside.Side());
			TPZCompElSide neighcompelside(list[is]);
            if (!withmem) {
                newcreatedinterface = new TPZMultiphysicsInterfaceElement(*fMesh,gel,index,thiscompelside,neighcompelside);
            }
            else
            {
                newcreatedinterface = new TPZCompElWithMem<TPZMultiphysicsInterfaceElement>(*fMesh,gel,index,thiscompelside,neighcompelside);
            }
		} else {
			//caso contrario ou caso ambos sejam de volume
			TPZCompElSide thiscompelside(this, thisside.Side());
			TPZCompElSide neighcompelside(list[is]);
            if (!withmem) {
                newcreatedinterface = new TPZMultiphysicsInterfaceElement(*fMesh,gel,index,neighcompelside,thiscompelside);
            }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	comp->AutoBuild();
	comp->Print(output);
  	output.flush();
    
    /**********************************************************************
     * Cria��o de uma malha computacional clone
     * ********************************************************************/
 	comp->GetRefPatches(patch);
	
	geomesh.ResetReference();
	TPZStack <int64_t> patchel;
	TPZStack <TPZGeoEl *> toclonegel;
	TPZStack <int64_t> patchindex;
	TPZVec<int64_t> n2elgraph;
	TPZVec<int64_t> n2elgraphid;
	TPZStack<int64_t> elgraph;
	TPZVec<int64_t> elgraphindex;
	int64_t k;
	TPZCompMesh *clonecmesh = new TPZCompMesh(&geomesh);
	cout << "Check 1: number of reference elements for patch before createcompel: " << patch.size() << endl;
    std::set<TPZGeoEl *>::iterator it;
    for (it=patch.begin(); it!=patch.end(); it++)
    {
		//patch[i]->Print(cout);
        int64_t index;
        TPZGeoEl *gel = *it;
        clonecmesh->CreateCompEl(gel, index);
        //		patch[i]->CreateCompEl(*clonecmesh,i);
	}
    //	cout << "Check 2: number of reference elements for patch after createcompel: " << patch.NElements() << endl;
	clonecmesh->CleanUpUnconnectedNodes();
    //	clonecmesh->Print(cout);
	clonecmesh->GetNodeToElGraph(n2elgraph,n2elgraphid,elgraph,elgraphindex);
	int64_t clnel = clonecmesh->NElements();
    //	cout << "Number of elements in clonemessh: " << clnel << endl;
	//o primeiro patch come�a em zero
	patchindex.Push(0);
	for (i=0; i<clnel; i++){
		//cout << endl << endl << "Evaluating patch for element: " << i << endl;
		clonecmesh->GetElementPatch(n2elgraph,n2elgraphid,elgraph,elgraphindex,i,patchel);
		cout << "Patch elements: " << patchel.NElements() << endl;
		/*for (k=0;k<patchel.NElements();k++){
         clonecmesh->ElementVec()[patchel[k]]->Reference()->Print();
         cout << endl;
         }*/
		for (j=0; j<patchel.NElements(); j++){
			//obten��o do elemento geom�trico do patch
			//cout << "Creating geometric clone elements for computational element :" << j << endl;
			TPZGeoEl *gel = clonecmesh->ElementVec()[patchel[j]]->Reference();
			//gel->Print(cout);
			//inserir todos os pais do elemento geom�trico do patch
			int64_t count = 0;
			//cout << "Inserting father element:" << "\t"; 
			while(gel){	
				TPZGeoEl *father = gel->Father();
				if (father){
					//father->Print(cout);
					gel = father;
					continue;
				}
				else toclonegel.Push(gel);
				gel = father;
				//cout <<  count << "\t";
				count ++;
			}
			//cout << endl;

//SingularElements(..)
void TPZAnalysisError::ZoomInSingularity(REAL csi, TPZCompElSide elside, REAL singularity_strength) {

	REAL hn = 1./pow(csi,1./singularity_strength);
	REAL Q=2.;
	REAL NcReal = log( 1.+(1./hn - 1.)*(Q - 1.) )/log(Q);
	int Nc = 0;
	while(REAL(Nc) < (NcReal+0.5)) Nc++;
	int minporder = 2;
	
	TPZStack<TPZCompElSide> ElToRefine;
	TPZStack<int> POrder;
	TPZStack<TPZGeoElSide> subelements;
	TPZStack<int64_t> csubindex;
	ElToRefine.Push(elside);
	POrder.Push(Nc);
	while(ElToRefine.NElements()) {
		/** Take the next element and its interpolation order from the stack*/
		TPZCompElSide curelside = ElToRefine.Pop();
		int curporder = POrder.Pop();
		if(!curelside.Exists()) continue;
		int64_t cindex = curelside.Element()->Index();
		if(cindex < 0) continue;
		
		/** Cast the element to an interpolated element if possible*/
		TPZCompEl *cel = curelside.Element();
		TPZInterpolatedElement *cintel = 0;
		cintel = dynamic_cast<TPZInterpolatedElement *> (cel);
		/** If the element is not interpolated, nothing to do */
		if(!cintel) continue;
		/** Set the interpolation order of the current element to curporder*/
		if(curporder == minporder) {
			cintel->PRefine(Nc);
			fSingular.Push(curelside);
		} else {
			cintel->PRefine(curporder);
			cintel->Divide(cindex,csubindex,1);
			/** Identify the subelements along the side and push them on the stack*/
		}
		TPZGeoElSide gelside = curelside.Reference();
		if(!gelside.Exists()) continue;
		gelside.GetSubElements2(subelements);
		int64_t ns = subelements.NElements();
		curporder--;
		int64_t is;
		for(is=0; is<ns; is++) {
			TPZGeoElSide sub = subelements[is];
			TPZCompElSide csub = sub.Reference();
			if(csub.Exists()) {
				ElToRefine.Push(csub);
				POrder.Push(curporder);
			}
		}
	}
	ExpandConnected(fSingular);
	
	/*
	 REAL H1_error,L2_error,estimate;
	 TPZBlock *flux=0;
	 int64_t nel = fElIndexes.NElements();
	 for(int64_t elloc=0;elloc<nel;elloc++) {
	 int64_t el = fElIndexes[elloc];
	 estimate = fElErrors[elloc];
	 REAL csi = estimate / fAdmissibleError;
	 REAL h = h_Parameter(intellist[el]);
	 REAL hn = h/pow(csi,1./.9);
	 REAL Nc = log( 1.+(h/hn - 1.)*(Q - 1.) )/log(Q);
	 if(hn > 1.3*h) hn = 2.0*h*hn / (h + hn);
	 REAL hsub = h;//100.0;//pode ser = h ; Cedric
	 TPZCompEl *locel = intellist[el];
	 //obter um subelemento que contem o ponto singular e tem tamanho <= hn
	 TPZAdmChunkVector<TPZCompEl *> sublist;
	 while(hsub > hn) {
	 TPZVec<int64_t> indexsubs;
	 int64_t index = locel->Index();
	 locel->Divide(index,indexsubs,1);
	 int64_t nsub = indexsubs.NElements();
	 TPZAdmChunkVector<TPZCompEl *> listsub(0);
	 for(int64_t k=0;k<nsub;k++) {
	 index = listsub.AllocateNewElement();
	 listsub[index] = Mesh()->ElementVec()[indexsubs[k]];
	 }
	 //existe um unico filho que contem o ponto singular
	 SingularElement(point,listsub,sublist);
	 hsub = h_Parameter(sublist[0]);
	 }
	 TPZInterpolatedElement *intel = (TPZInterpolatedElement *) locel;
	 intel->PRefine(Nc+1);
	 indexlist.Push(intel->Index());
	 //os elemento viz devem ter ordens menores a cel quanto mais longe de point
	 TPZInterpolatedElement *neighkeep,*neigh;
	 //feito s�para o caso 1d , extender para o caso geral
	 int dim = intel->Dimension();
	 if(dim != 1) {
	 cout << "TPZAnalysisError::Step3 not dimension implemented , dimension = " << intellist[el]->Dimension() << endl;
	 return ;//exit(1);
	 }
	 for(int side=0;side<2;side++) {
	 int ly = 1;
	 TPZGeoElSide neighside = intel->Reference()->Neighbour(side);
//.........这里部分代码省略.........

void TPZAdaptMesh::BuildReferencePatch() {
    
    // the fGeoRef elements are a partition of the computational domain (should be)
    // create a computational element based on each reference element
    TPZGeoMesh *gmesh = fReferenceCompMesh->Reference();
    gmesh->ResetReference();
    TPZCompMesh *tmpcmesh = new TPZCompMesh (gmesh);
    int i,j;
    for (i=0;i<fGeoRef.NElements();i++){
        long index;
        tmpcmesh->CreateCompEl(fGeoRef[i],index);
    } 
    tmpcmesh->CleanUpUnconnectedNodes();
	tmpcmesh->ExpandSolution();
    TPZStack <long> patchelindex;
    TPZStack <TPZGeoEl *> toclonegel;
    TPZStack<long> elgraph;
    TPZVec<long> n2elgraph;
    TPZVec<long> n2elgraphid;
    TPZVec<long> elgraphindex;

    tmpcmesh->GetNodeToElGraph(n2elgraph,n2elgraphid,elgraph,elgraphindex);
    // we use the  node to elgraph structure to decide which elements will be included
    int clnel = tmpcmesh->NElements();
    // clnel corresponds to the number of patches
    // fPatch and fPatchIndex form a compacted list which form the patches.
    // Boundary elements will be added to each patch.
    fPatchIndex.Push(0);
    for (int ipatch=0; ipatch<clnel; ipatch++){
        tmpcmesh->GetElementPatch(n2elgraph,n2elgraphid,elgraph,elgraphindex,ipatch,patchelindex);
        for (j=0; j<patchelindex.NElements(); j++){
            TPZGeoEl *gel = tmpcmesh->ElementVec()[patchelindex[j]]->Reference();
            //      int count = 0;
            if(gel) fPatch.Push(gel);
        }
        int sum = fPatch.NElements();
        fPatchIndex.Push(sum);
    }
	
#ifdef DEBUG2 
	// CAJU TOOL
	{
		std::string filename("cMeshVtk.");
		{
			std::stringstream finalname;
			finalname << filename << 0 << ".vtk";
			ofstream file(finalname.str().c_str());
			/** @brief Generate an output of all geometric elements that have a computational counterpart to VTK */
			//static void PrintCMeshVTK(TPZGeoMesh *gmesh, std::ofstream &file, bool matColor = false);
			TPZVTKGeoMesh::PrintCMeshVTK(gmesh,file,true);
		}
		for (int ip=0; ip<clnel; ip++) {
			int firstindex = fPatchIndex[ip];
			int lastindex = fPatchIndex[ip+1];
			gmesh->ResetReference();
			tmpcmesh->LoadReferences();
			std::set<TPZGeoEl *> loaded;
			for (int ind=firstindex; ind<lastindex; ind++) {
				TPZGeoEl *gel = fPatch[ind];
				loaded.insert(gel);
			}
			int ngel = gmesh->NElements();
			for (int el=0; el<ngel; el++) {
				TPZGeoEl *gel = gmesh->ElementVec()[el];
				if (!gel) {
					continue;
				}
				if (gel->Reference() && loaded.find(gel) == loaded.end()) {
					gel->ResetReference();
				}
			}
			std::stringstream finalname;
			finalname << filename << ip+1 << ".vtk";
			ofstream file(finalname.str().c_str());
			/** @brief Generate an output of all geometric elements that have a computational counterpart to VTK */
			//static void PrintCMeshVTK(TPZGeoMesh *gmesh, std::ofstream &file, bool matColor = false);
			TPZVTKGeoMesh::PrintCMeshVTK(gmesh,file,true);
			
		}
	}
#endif
	// cleaning reference to computational elements into temporary cmesh
    gmesh->ResetReference();
    delete tmpcmesh;
	// loading references between geometric and computational meshes (originals)
    fReferenceCompMesh->LoadReferences();
}