C++ TPZCompMesh::CleanUpUnconnectedNodes方法代码示例

TPZCompMesh *MalhaCompDois(TPZGeoMesh * gmesh, int pOrder, bool isdiscontinuous)
{
    /// criar materiais
    int dim = 2;
    TPZMatPoisson3d *material;
    material = new TPZMatPoisson3d(matId,dim);
    TPZMaterial * mat(material);
    
    material->SetNoPenalty();
    material->SetNonSymmetric();
    
    REAL diff = -1.;
    REAL conv = 0.;
    TPZVec<REAL> convdir(3,0.);
    REAL flux = 0.;
    
    material->SetParameters(diff, conv, convdir);
    material->SetInternalFlux(flux);
    material->NStateVariables();
    
    TPZCompEl::SetgOrder(pOrder);
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    cmesh->SetDimModel(dim);
    //cmesh->SetAllCreateFunctionsContinuous();
    cmesh->InsertMaterialObject(mat);
    
    TPZAutoPointer<TPZFunction<STATE> > forcef = new TPZDummyFunction<STATE>(ForcingF, 5);
    material->SetForcingFunction(forcef);
    
    ///Inserir condicao de contorno
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    TPZMaterial * BCond0 = material->CreateBC(mat, bc0,dirichlet, val1, val2);
    TPZMaterial * BCond2 = material->CreateBC(mat, bc2,dirichlet, val1, val2);
    TPZMaterial * BCond1 = material->CreateBC(mat, bc1,dirichlet, val1, val2);
    TPZMaterial * BCond3 = material->CreateBC(mat, bc3,dirichlet, val1, val2);
    
    cmesh->InsertMaterialObject(BCond0);
    cmesh->InsertMaterialObject(BCond1);
    cmesh->InsertMaterialObject(BCond2);
    cmesh->InsertMaterialObject(BCond3);
    
    //Ajuste da estrutura de dados computacional
    if (isdiscontinuous==true) {
        //Set discontinuous functions
        cmesh->SetAllCreateFunctionsDiscontinuous();
        cmesh->AutoBuild();
        cmesh->ExpandSolution();
        cmesh->AdjustBoundaryElements();
        cmesh->CleanUpUnconnectedNodes();
    }
    else{
        cmesh->SetAllCreateFunctionsContinuous();
        cmesh->AutoBuild();
        cmesh->ExpandSolution();
        cmesh->AdjustBoundaryElements();
        cmesh->CleanUpUnconnectedNodes();
    }
    
    return cmesh;
}

/**
 * @brief transform in low order Raviar Tomas
 */
void TPZCreateApproximationSpace::UndoMakeRaviartTomas(TPZCompMesh &cmesh)
{
    int numcell = cmesh.NElements();
    int el;
    for (el = 0; el<numcell ; el++) {
        TPZCompEl *cel = cmesh.ElementVec()[el];
        TPZInterpolatedElement *intel = dynamic_cast<TPZInterpolatedElement *>(cel);
        if (!intel) {
            continue;
        }
        TPZGeoEl *gel = intel->Reference();
        int geldim = gel->Dimension();
        int is;
        int nsides = gel->NSides();
        for (is=0; is<nsides; is++) {
            if (gel->SideDimension(is) != geldim-1) {
                continue;
            }
            int nsconnects = intel->NSideConnects(is);
            // only interested in HDiv elements
            if (nsconnects != 1) {
                continue;
            }
//            int cindex = intel->SideConnectIndex(0, is);
            TPZConnect &c = intel->Connect(intel->SideConnectLocId(0,is));
            if (c.HasDependency()) {
                c.RemoveDepend();
            }
        }
    }
    cmesh.ExpandSolution();
    cmesh.CleanUpUnconnectedNodes();
}

TPZCompMesh * ComputationalMesh(TPZGeoMesh * gmesh, int p)
{
    int matid = 1;
    int dim = 2;
    REAL wavespeed = 1.0;

    ///Computational Mesh
    TPZCompEl::SetgOrder(p);
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    cmesh->SetDimModel(dim);    
    cmesh->SetAllCreateFunctionsContinuous();
    
    TPZMaterial * Air = new TPZLinearWave(matid,dim);
    cmesh->InsertMaterialObject(Air);
    
    {
        //Boundary Conditions
        TPZFMatrix<STATE> k1(dim,dim,0.), k2(dim,dim,0.);
        TPZMaterial * BCD = Air->CreateBC(Air, 2, 0, k1, k2);
        cmesh->InsertMaterialObject(BCD);
        
        TPZMaterial * BCN = Air->CreateBC(Air, 3, 1, k1, k2);
        cmesh->InsertMaterialObject(BCN);
    }   
        
    cmesh->AutoBuild();
    cmesh->AdjustBoundaryElements();
    cmesh->CleanUpUnconnectedNodes();    
    
    return cmesh;
}

TPZCompMesh *CompMesh1D(TPZGeoMesh *gmesh,int p, TPZMaterial *material,TPZVec<int> &bc,TPZVec<int> &bcType) {
	if(!material || bc.NElements()<2 || bcType.NElements() != bc.NElements()) return NULL;
	int dim = 1;
	
	
	TPZAutoPointer<TPZMaterial> mat(material);
	
	// related to interpolation space
	TPZCompEl::SetgOrder(p);
	TPZCompMesh *cmesh = new TPZCompMesh(gmesh);
	cmesh->SetDimModel(dim);
	cmesh->SetAllCreateFunctionsContinuous();
	cmesh->InsertMaterialObject(mat);
	
	// Related to boundary conditions
	//	REAL uN=1-cosh(1.)/sinh(1.);
	TPZFMatrix<STATE> val1(1,1,0.), val2(1,1,0.);
	if(!bcType[0])  // dirichlet
		val2.PutVal(0,0,0.0);
	TPZAutoPointer<TPZMaterial> BCond1 = material->CreateBC(mat, bc[0],bcType[0], val1, val2);
	cmesh->InsertMaterialObject(BCond1);
	
	if(!bcType[1])  // dirichlet
		val2.PutVal(0,0,0.0);
	TPZAutoPointer<TPZMaterial> BCond2 = material->CreateBC(mat, bc[1],bcType[1], val1, val2);
	cmesh->InsertMaterialObject(BCond2);
	
	//Adjusting data
	cmesh->AutoBuild();
	cmesh->AdjustBoundaryElements(); 
	cmesh->CleanUpUnconnectedNodes();
	
	return cmesh;
}

TPZCompMesh *MalhaMultifisicaOpt(TPZVec<TPZCompMesh *> meshvec, TPZGeoMesh *gmesh){
    
    //Creating computational mesh for multiphysic elements
    gmesh->ResetReference();
    TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
    
    //criando material
    int dim =2;
    
    TPZMatPoissonControl *material = new TPZMatPoissonControl(MatId,dim);
    
    
    //incluindo os dados do problema
    REAL k=1;
    REAL alpha=1;
    material-> SetParameters( k, alpha);
    
    //solucao exata
    TPZAutoPointer<TPZFunction<STATE> > solexata;
    solexata = new TPZDummyFunction<STATE>(StateAd, 5);
    material->SetForcingFunctionExact(solexata);
    
    //funcao do lado direito da equacao do problema
    TPZAutoPointer<TPZFunction<STATE> > force;
    TPZDummyFunction<STATE> *dum;
    dum = new TPZDummyFunction<STATE>(OptForcing, 5);
    dum->SetPolynomialOrder(20);
    force = dum;
    material->SetForcingFunction(force);
    
    //inserindo o material na malha computacional
    TPZMaterial *mat(material);
    mphysics->InsertMaterialObject(mat);
    mphysics->SetDimModel(dim);
    
    //Criando condicoes de contorno
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    
    
    TPZMaterial * BCond0 = material->CreateBC(mat, bc0, bcdirichlet, val1, val2);
    TPZMaterial * BCond1 = material->CreateBC(mat, bc1, bcdirichlet, val1, val2);
    TPZMaterial * BCond2 = material->CreateBC(mat, bc2, bcdirichlet, val1, val2);
    TPZMaterial * BCond3 = material->CreateBC(mat, bc3, bcdirichlet, val1, val2);
    
    ///Inserir condicoes de contorno
    mphysics->InsertMaterialObject(BCond0);
    mphysics->InsertMaterialObject(BCond1);
    mphysics->InsertMaterialObject(BCond2);
    mphysics->InsertMaterialObject(BCond3);
    
    mphysics->SetAllCreateFunctionsMultiphysicElem();
    
    //Fazendo auto build
    mphysics->AutoBuild();
    mphysics->AdjustBoundaryElements();
    mphysics->CleanUpUnconnectedNodes();
    
    TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
    TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
    TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
    
    return mphysics;
    
}

//.........这里部分代码省略.........
    int in;
    nodeindex.Resize(numnos[nel]);
    for(in=0; in<numnos[nel]; in++) {
      nodeindex[in] = nodind[nel][in];
    }
    int index;
    switch(nel) {
    case 0:
      //      elvec[el] = gmesh->CreateGeoElement(ECube,nodeindex,1,index);
//      gelvec[nel]=new TPZGeoElC3d(nodeindex,matid,*gmesh);
      break;
    case 1:
      gelvec[nel] = gmesh->CreateGeoElement(EPiramide,nodeindex,matid,index);
      //       gelvec[nel]=new TPZGeoElPi3d(nodeindex,matid,*gmesh);
      break;
    case 2:
      gelvec[nel] = gmesh->CreateGeoElement(ETetraedro,nodeindex,matid,index);
    //       gelvec[nel]=new TPZGeoElT3d(nodeindex,matid,*gmesh);
      break;
    case 3:
//       gelvec[nel]=new TPZGeoElPr3d(nodeindex,matid,*gmesh);
//      gelvec[nel] = gmesh->CreateGeoElement(EPrisma,nodeindex,matid,index);
      break;
    case 4:
      //      gelvec[nel]=new TPZGeoEl1d(nodeindex,2,*gmesh);
      break;
    case 5:
      //      gelvec[nel]=new TPZGeoElQ2d(nodeindex,3,*gmesh);
      break;
    case 6:
      //      gelvec[nel]=new TPZGeoElT2d(nodeindex,3,*gmesh);
      break;
    default:
      break;
    }
  }
  gmesh->BuildConnectivity2();

  //TPZVec<TPZGeoEl *> sub;
  //elvec[0]->Divide(sub);
  //   	elvec[1]->Divide(sub);
  //   	elvec[2]->Divide(sub);

//  TPZGeoElBC gbc;

  // bc -1 -> Dirichlet
//  TPZGeoElBC gbc1(gelvec[0],20,-1,*gmesh);
  TPZGeoElBC gbc11(gelvec[1],14,-1,*gmesh);
//  TPZGeoElBC gbc12(gelvec[3],15,-1,*gmesh);



  // bc -2 -> Neumann at the right x==1
//  TPZGeoElBC gbc2(gelvec[0],25,-2,*gmesh);
//  TPZGeoElBC gbc21(gelvec[3],19,-2,*gmesh);
  TPZGeoElBC gbc22(gelvec[2],10,-2,*gmesh);

  TPZCompMesh *cmesh = new TPZCompMesh(gmesh);

  TPZAutoPointer<TPZMaterial> mat;
  //  if(nstate == 3) {
    mat = new TPZMaterialTest3D(1);
    TPZFMatrix mp (3,1,0.);
    TPZMaterialTest3D * mataux = dynamic_cast<TPZMaterialTest3D *> (mat.operator ->());
    TPZMaterialTest3D::geq3=1;
    mataux->SetMaterial(mp);
    /*  } else {
    TPZMat2dLin *mat2d = new TPZMat2dLin(1);
    int ist,jst;
    TPZFMatrix xk(nstate,nstate,1.),xc(nstate,nstate,0.),xf(nstate,1,0.);
    for(ist=0; ist<nstate; ist++) {
      if(nstate != 1) xf(ist,0) = 1.;
      for(jst=0; jst<nstate; jst++) {
	if(ist != jst) xk(ist,jst) = 0.;
      }
    }
    mat2d->SetMaterial(xk,xc,xf);
    mat = mat2d;
    }*/

  TPZFMatrix val1(3,3,0.),val2(3,1,0.);
  TPZAutoPointer<TPZMaterial> bc[2];

  bc[0] = mat->CreateBC(mat,-1,0,val1,val2);
  val2(0,0) = 1.;
  bc[1] = mat->CreateBC(mat,-2,1,val1,val2);
  cmesh->InsertMaterialObject(mat);

  int i;
  for(i=0; i<2; i++) cmesh->InsertMaterialObject(bc[i]);

  gmesh->Print(cout);

  cmesh->AutoBuild();
  cmesh->AdjustBoundaryElements();
  cmesh->CleanUpUnconnectedNodes();

  gmesh->Print(cout);
  return cmesh;
}

//malha multifisica para o metodo da dupla projecao
TPZCompMesh *MalhaMDP(TPZVec<TPZCompMesh *> meshvec,TPZGeoMesh * gmesh){
    
    //Creating computational mesh for multiphysic elements
	gmesh->ResetReference();
	TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
    
    //criando material
    int dim =2;
    TPZMDPMaterial *material = new TPZMDPMaterial(1,dim);
    
    //incluindo os dados do problema
    REAL coefk = 1.;
    material->SetParameters(coefk, 0.);
    
    //solucao exata
    TPZAutoPointer<TPZFunction<STATE> > solexata;
    solexata = new TPZDummyFunction<STATE>(SolSuave);
    material->SetForcingFunctionExact(solexata);
    
    //funcao do lado direito da equacao do problema
    TPZAutoPointer<TPZFunction<STATE> > force;
    TPZDummyFunction<STATE> *dum;
    dum = new TPZDummyFunction<STATE>(ForceSuave);
    dum->SetPolynomialOrder(20);
    force = dum;
    material->SetForcingFunction(force);
    
    //inserindo o material na malha computacional
    TPZMaterial *mat(material);
    mphysics->InsertMaterialObject(mat);
    
    //Criando condicoes de contorno
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    int boundcond = dirichlet;
    //BC -1
    TPZMaterial * BCondD1 = material->CreateBC(mat, bc1,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet1 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD1->SetForcingFunction(bcmatDirichlet1);
    mphysics->InsertMaterialObject(BCondD1);
    
    //BC -2
	TPZMaterial * BCondD2 = material->CreateBC(mat, bc2,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet2 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD2->SetForcingFunction(bcmatDirichlet2);
    mphysics->InsertMaterialObject(BCondD2);
    
    //BC -3
	TPZMaterial * BCondD3 = material->CreateBC(mat, bc3,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet3 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD3->SetForcingFunction(bcmatDirichlet3);
    mphysics->InsertMaterialObject(BCondD3);
    
    //BC -4
	TPZMaterial * BCondD4 = material->CreateBC(mat, bc4,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet4 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD4->SetForcingFunction(bcmatDirichlet4);
    mphysics->InsertMaterialObject(BCondD4);

    
    mphysics->InsertMaterialObject(BCondD1);
    mphysics->InsertMaterialObject(BCondD2);
    mphysics->InsertMaterialObject(BCondD3);
    mphysics->InsertMaterialObject(BCondD4);
    
    
    //set multiphysics element
    mphysics->SetDimModel(dim);
    mphysics->SetAllCreateFunctionsMultiphysicElem();
    
    //Fazendo auto build
    mphysics->AutoBuild();
	mphysics->AdjustBoundaryElements();
	mphysics->CleanUpUnconnectedNodes();
    
    // Creating multiphysic elements into mphysics computational mesh
	TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
	TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
	TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
    
    return mphysics;
}

TPZCompMesh * CompMesh(TPZGeoMesh *gmesh, int porder)
{
    /// criar materiais
	int dim = gmesh->Dimension();
    
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    
    TPZMatLaplacian *material = new TPZMatLaplacian(1,dim);
    
//    TPZAutoPointer<TPZFunction<REAL> > forcef = new TPZDummyFunction<REAL>(ForceSuave);
//    material->SetForcingFunction(forcef);
    
    TPZAutoPointer<TPZFunction<STATE> > force;
    TPZDummyFunction<STATE> *dum;
    dum = new TPZDummyFunction<STATE>(ForceSuave);
    dum->SetPolynomialOrder(20);
    force = dum;
    material->SetForcingFunction(force);
    
    
    TPZAutoPointer<TPZFunction<STATE> > solExata= new TPZDummyFunction<STATE>(SolSuave);
    material->SetForcingFunctionExact(solExata);
    
    TPZMaterial * mat(material);
    cmesh->InsertMaterialObject(mat);
    
    cmesh->SetDimModel(dim);
    cmesh->SetDefaultOrder(porder);

	///Inserir condicao de contorno
	TPZFMatrix<STATE> val1(2,2,1.), val2(2,1,0.);
	
    //BC -1
    TPZMaterial * BCondD1 = material->CreateBC(mat, bc1,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet1 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD1->SetForcingFunction(bcmatDirichlet1);
    cmesh->InsertMaterialObject(BCondD1);
    
    //BC -2
	TPZMaterial * BCondD2 = material->CreateBC(mat, bc2,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet2 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD2->SetForcingFunction(bcmatDirichlet2);
    cmesh->InsertMaterialObject(BCondD2);
    
    //BC -3
	TPZMaterial * BCondD3 = material->CreateBC(mat, bc3,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet3 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD3->SetForcingFunction(bcmatDirichlet3);
    cmesh->InsertMaterialObject(BCondD3);
    
    //BC -4
	TPZMaterial * BCondD4 = material->CreateBC(mat, bc4,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet4 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD4->SetForcingFunction(bcmatDirichlet4);
    cmesh->InsertMaterialObject(BCondD4);
    
    
    //Fazendo auto build
    cmesh->SetAllCreateFunctionsContinuous();
    cmesh->AutoBuild();
	cmesh->AdjustBoundaryElements();
	cmesh->CleanUpUnconnectedNodes();
    
    return cmesh;
}

int TPZGeoCloneMesh::main(){
	cout << "**************************************" << endl;
  	cout << "****** Getting Patchs!************" << endl;
	cout << "**************************************" << endl;
    
    /*******************************************************
     * Constru��o da malha
     * *****************************************************/
  	//malha quadrada de nr x nc
	const int numrel = 3;
  	const int numcel = 3;
//  	int numel = numrel*numcel;
  	TPZVec<REAL> coord(2,0.);
    
  	// criar um objeto tipo malha geometrica
  	TPZGeoMesh geomesh;
    
  	// criar nos
  	int i,j;
  	for(i=0; i<(numrel+1); i++) {
        for (j=0; j<(numcel+1); j++) {
            int64_t nodind = geomesh.NodeVec().AllocateNewElement();
            TPZVec<REAL> coord(2);
            coord[0] = j;//co[nod][0];
            coord[1] = i;//co[nod][1];
            geomesh.NodeVec()[nodind] = TPZGeoNode(i*(numrel+1)+j,coord,geomesh);
        }
  	}
  	// cria��o dos elementos
  	int elc, elr;
  	TPZGeoEl *gel[numrel*numcel];
  	TPZVec<int64_t> indices(4);
  	for(elr=0; elr<numrel; elr++) {  
        for(elc=0; elc<numcel; elc++) {
            indices[0] = (numrel+1)*elr+elc;
            indices[1] = indices[0]+1;
            indices[3] = indices[0]+numrel+1;
            indices[2] = indices[1]+numrel+1;
            // O proprio construtor vai inserir o elemento na malha
			int64_t index;
			gel[elr*numrel+elc] = geomesh.CreateGeoElement(EQuadrilateral,indices,1,index);
            //gel[elr*numrel+elc] = new TPZGeoElQ2d(elr*numrel+elc,indices,1,geomesh);
        }
  	}
	//Divis�o dos elementos
  	TPZVec<TPZGeoEl *> sub;
  	gel[0]->Divide(sub);
    //  	gel[1]->Divide(sub);
    //  	gel[3]->Divide(sub); 
  	ofstream output("patches.dat");
  	geomesh.Print(output);
    //  	TPZGeoElBC t3(gel[0],4,-1,geomesh); 
    //  	TPZGeoElBC t4(gel[numel-1],6,-2,geomesh); 
  	geomesh.Print(output);
	geomesh.BuildConnectivity();
    std::set <TPZGeoEl *> patch;
    
	
  	TPZCompMesh *comp = new TPZCompMesh(&geomesh);
 	// inserir os materiais
  	TPZMaterial *meumat = new TPZElasticityMaterial(1,1.e5,0.2,0,0);
  	comp->InsertMaterialObject(meumat);
  	// inserir a condicao de contorno
    //  	TPZFMatrix val1(3,3,0.),val2(3,1,0.);
    //  	TPZMaterial *bnd = meumat->CreateBC (-1,0,val1,val2);
    //  	comp->InsertMaterialObject(bnd);
    //  	TPZFMatrix val3(3,3,1);
    // 	bnd = meumat->CreateBC (-2,1,val3,val2);
    //  	comp->InsertMaterialObject(bnd);
	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();
//.........这里部分代码省略.........

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();
}