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 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 TPZAdaptMesh::CreateClones(){
	// asserting references of the original meshes
    fReferenceCompMesh->Reference()->ResetReference();
    fReferenceCompMesh->LoadReferences();
    TPZGeoMesh *geomesh = fReferenceCompMesh->Reference();
    
    TPZStack<TPZGeoEl*> patch;
    
    int clid,elid;
    for (clid=0; clid<fPatchIndex.NElements()-1;clid++) {
		// making clone of the original geometric mesh, only to construct computational clone
        TPZGeoCloneMesh *geoclone = new TPZGeoCloneMesh(geomesh);
        TPZStack<TPZGeoEl*> patch;
        for (elid=fPatchIndex[clid];elid<fPatchIndex[clid+1];elid++){
            patch.Push(fPatch[elid]);
        }
        geoclone->SetElements(patch,fGeoRef[clid]);
        TPZVec<TPZGeoEl *> sub;
        //    int ngcel = geoclone->ElementVec().NElements();
        int printing = 0;
        if(printing) {
            ofstream out("testAdaptMesh.txt",ios::app);
            geoclone->Print(out);
        }
        
        TPZCompCloneMesh *clonecompmesh = new TPZCompCloneMesh(geoclone,fReferenceCompMesh);
        clonecompmesh->AutoBuild(/*fMaxP*/);
		// Computational mesh clone is stored
        fCloneMeshes.Push(clonecompmesh);    
    }
}

	void TPZPrism::LowerDimensionSides(int side,TPZStack<int> &smallsides)
	{
		smallsides.Resize(0);
		int nsidecon = NContainedSides(side);
		int is;
		for(is=0; is<nsidecon-1; is++)
			smallsides.Push(ContainedSideLocId(side,is));
	}

void TPZGeoCloneMesh::Write(TPZStream &buf, int withclassid) const
{
    TPZGeoMesh::Write(buf,withclassid);
	try
	{

        if(!fGeoReference) {
            std::cout << "Cloned geo mesh without geometric mesh from which this mesh is cloned." << std::endl;
            DebugStop();
        }
        TPZPersistenceManager::WritePointer(fGeoReference, &buf);
        
        buf.Write(fMapNodes);
        
        // Maps elements with original elements
        std::map<int64_t,int64_t> MappingElements;
        TPZGeoEl *gel;
        int64_t indexorig, indexcloned;
        std::map<TPZGeoEl* , TPZGeoEl* >::const_iterator it;
        for(it=fMapElements.begin();it!=fMapElements.end();it++) {
            gel = it->first;
            indexorig = gel->Index();
            gel = it->second;
            indexcloned = gel->Index();
            MappingElements.insert(std::make_pair(indexorig,indexcloned));
        }
        buf.Write(MappingElements);
        
        // Writing index of the elements in fReferenceElement
        TPZStack<int64_t> RefElements;
        int64_t sz = fReferenceElement.size();
        for(int64_t ii=0;ii<sz;ii++) {
            RefElements.push_back(fReferenceElement[ii]->Index());
        }
        buf.Write(RefElements);
        // Have to save a compact structure of the vector??

        // Writing index of the elements in fPatchElement
        std::set<int> PatchElements;
        std::set<TPZGeoEl* >::iterator itpatch = fPatchElements.begin();
        while(itpatch != fPatchElements.end()) {
            PatchElements.insert((*itpatch)->Index());
            itpatch++;
        }
        buf.Write(PatchElements);
        
        int64_t rootindex = fGeoRoot->Index();
        buf.Write(&rootindex,1);
	}
	catch(const exception& e)
	{
		cout << "Exception catched! " << e.what() << std::endl;
		cout.flush();
		DebugStop();
	}
}//method

	void TPZRefQuad::GetSubElements(TPZGeoEl *father,int side, TPZStack<TPZGeoElSide> &subel){
		
		subel.Resize(0);
		if(side<0 || side>TPZShapeQuad::NSides || !father->HasSubElement()){
			PZError << "TPZRefQuad::GetSubelements called with error arguments\n";
			return;
		}
		int nsub = NSideSubElements(side);
		for(int i=0;i<nsub;i++)
			subel.Push(TPZGeoElSide(father->SubElement(subeldata[side][i][0]),subeldata[side][i][1]));
	}

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

void RefinamentoSingular(TPZAutoPointer<TPZGeoMesh> gmesh,int nref)
{
    int64_t nnodes = gmesh->NNodes();
    int64_t in;
    for (in=0; in<nnodes; in++) {
        TPZGeoNode *gno = &gmesh->NodeVec()[in];
        if (abs(gno->Coord(0))< 1.e-6 && abs(gno->Coord(1)) < 1.e-6) {
            break;
        }
    }
    if (in == nnodes) {
        DebugStop();
    }
    TPZGeoElSide gelside;
    int64_t nelem = gmesh->NElements();
    for (int64_t el = 0; el<nelem; el++) {
        TPZGeoEl *gel = gmesh->ElementVec()[el];
        int ncorner = gel->NCornerNodes();
        for (int ic=0; ic<ncorner; ic++) {
            int64_t nodeindex = gel->NodeIndex(ic);
            if (nodeindex == in) {
                gelside = TPZGeoElSide(gel, ic);
                break;
            }
        }
        if (gelside.Element()) {
            break;
        }
    }
    if (!gelside.Element()) {
        DebugStop();
    }
    for (int iref = 0; iref <nref; iref++) {
        TPZStack<TPZGeoElSide> gelstack;
        gelstack.Push(gelside);
        TPZGeoElSide neighbour = gelside.Neighbour();
        while (neighbour != gelside) {
            gelstack.Push(neighbour);
            neighbour = neighbour.Neighbour();
        }
        int64_t nstack = gelstack.size();
        for (int64_t ist=0; ist < nstack; ist++) {
            if (!gelstack[ist].Element()->HasSubElement()) {
                TPZVec<TPZGeoEl *> subel;
                gelstack[ist].Element()->Divide(subel);
            }
        }
    }
}

void TPZGeoCloneMesh::Read(TPZStream &buf, void *context)
{
    TPZGeoMesh::Read(buf,context);
	try
	{
		
        fGeoReference = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::GetInstance(&buf));
        
        buf.Read<int64_t>(fMapNodes);
        
        std::map<int64_t,int64_t> MappingElements;
        
        TPZGeoEl *gelorig, *gelcloned;
        buf.Read<int64_t>(MappingElements);
        std::map<int64_t,int64_t>::iterator it = MappingElements.begin();
        while(it != MappingElements.end()) {
            gelorig = fGeoReference->ElementVec()[it->first];
            gelcloned = ElementVec()[it->second];
            fMapElements.insert(std::make_pair(gelorig,gelcloned));
            it++;
        }
        
        // Writing index of the elements in fReferenceElement
        TPZStack<int64_t> RefElements;
        buf.Read(RefElements);
        for(int64_t ii=0;ii<RefElements.size();ii++) {
            fReferenceElement.push_back(fGeoReference->ElementVec()[RefElements[ii]]);
        }

        // Reading index of the elements in fPatchElement
        std::set<int> PatchElements;
        buf.Read(PatchElements);
        std::set<int>::iterator itpatch = PatchElements.begin();
        while(itpatch != PatchElements.end()) {
            fPatchElements.insert(fGeoReference->ElementVec()[*itpatch]);
            itpatch++;
        }
        
        int64_t indexroot;
        buf.Read(&indexroot);
        fGeoRoot = ElementVec()[indexroot];

	}
	catch(const exception& e)
	{
		cout << "Exception catched! " << e.what() << std::endl;
		cout.flush();
	}
}

	void TPZPrism::HigherDimensionSides(int side, TPZStack<int> &high)
	{
		if(side <0 || side >= NSides) {
			PZError << "TPZPrism::HigherDimensionSides side "<< side << endl;
		}
		int is;
		for(is=0; is<nhighdimsides[side]; is++) high.Push(highsides[side][is]);
		
	}

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

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

//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);
//.........这里部分代码省略.........

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();
//.........这里部分代码省略.........