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

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::MakeRaviartTomas(TPZCompMesh &cmesh)
{
    int numcell = cmesh.NElements();
    int el;
    for (el = 0; el<numcell ; el++) {
        TPZCompEl *cel = cmesh.ElementVec()[el];
        TPZInterpolationSpace *intel = dynamic_cast<TPZInterpolationSpace *>(cel);
        if (!intel) {
            continue;
        }
        intel->SetPreferredOrder(1);
    }
    cmesh.ExpandSolution();
    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()) {
                continue;
            }
            int nshape = 1;
            int nstate = 1;
            int order = 0;
            int cindex2 = cmesh.AllocateNewConnect(nshape, nstate, order);
//            TPZConnect &c2 = cmesh.ConnectVec()[cindex];
            TPZFNMatrix<2> depmat(2,1,1.);
            c.AddDependency(cindex, cindex2, depmat, 0, 0, 2, 1);
        }
    }
    cmesh.ExpandSolution();
}

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

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