C++ TPZStack类代码示例

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