C++ TPZCompMesh类代码示例

TPZIntelGen<TSHAPE>::TPZIntelGen(TPZCompMesh &mesh, TPZGeoEl *gel, int &index) :
TPZInterpolatedElement(mesh,gel,index) {
	int i;
	fPreferredOrder = mesh.GetDefaultOrder();
	for(i=0;i<TSHAPE::NSides-TSHAPE::NCornerNodes;i++) {
		//    fSideOrder[i] = gOrder;
	}
	for(i=0; i<TSHAPE::NSides; i++) fConnectIndexes[i]=-1;
	//  RemoveSideRestraintsII(EInsert);
	gel->SetReference(this);
	for(i=0;i<TSHAPE::NCornerNodes;i++) {
		fConnectIndexes[i] = CreateMidSideConnect(i);
		mesh.ConnectVec()[fConnectIndexes[i]].IncrementElConnected();
	}
	for(;i<TSHAPE::NSides;i++) {
		fConnectIndexes[i] = CreateMidSideConnect(i);
		mesh.ConnectVec()[fConnectIndexes[i]].IncrementElConnected();
		IdentifySideOrder(i);
	}
	
	int sideorder = SideOrder(TSHAPE::NSides-1);
	sideorder = 2*sideorder;
	if (sideorder > fIntRule.GetMaxOrder()) sideorder = fIntRule.GetMaxOrder();
	//  TPZManVector<int,3> order(3,2*sideorder+2);
	TPZManVector<int,3> order(3,sideorder);
	//TPZManVector<int,3> order(3,20);
	fIntRule.SetOrder(order);
	
}

///Funcao principal do programa
int main(int argc, char *argv[])
{
    int dim = 2;//dimensao do problema
    int uNDiv=3,
    vNDiv=4;//numero de divisoes feitas no dominio
	int nel = uNDiv*vNDiv; //numero de elementos a serem utilizados
    int pOrder = 4; //ordem polinomial de aproximacao
	TPZGeoMesh *gmesh = CreateGMesh(nel, uNDiv, vNDiv); //funcao para criar a malha geometrica
	
	TPZCompMesh *cmesh = CMesh(gmesh, pOrder); //funcao para criar a malha computacional

	// Resolvendo o Sistema
    bool optimizeBandwidth = false; //impede a renumeracao das equacoes do problema(para obter o mesmo resultado do Oden)
	TPZAnalysis an(cmesh, optimizeBandwidth); //cria objeto de analise que gerenciaria a analise do problema
	an.Run();//assembla a matriz de rigidez (e o vetor de carga) global e inverte o sistema de equacoes
    
	TPZFMatrix<STATE> solucao=cmesh->Solution();//Pegando o vetor de solucao, alphaj
	solucao.Print("Sol",cout,EMathematicaInput);//imprime na formatacao do Mathematica
	
    //fazendo pos processamento para paraview
    TPZStack<string> scalnames, vecnames;
    scalnames.Push("State");//setando para imprimir u
    string plotfile = "ModelProblemSol.vtk";//arquivo de saida que estara na pasta debug
    an.DefineGraphMesh(dim, scalnames, vecnames, plotfile);//define malha grafica
	int postProcessResolution = 3;//define resolucao do pos processamento
    an.PostProcess(postProcessResolution);//realiza pos processamento
    
	std::cout << "FINISHED!" << std::endl;
	
	return 0;
}

TPZCompMesh *TetraMesh(){
  REAL Coord [8][3] = {
    {0.,0.,0.},{1.,0.,0.},{1.,1.,0.},{0.,1.,0.},
    {0.,0.,1.},{1.,0.,1.},{1.,1.,1.},{0.,1.,1.}
  };

  int Connects [5][4] = {
    {0,1,3,4},
    {1,2,3,6},
    {5,6,4,1},
    {7,6,4,3},
    {1,3,4,6}
  };

  int i,j;
  TPZGeoMesh *gmesh = new TPZGeoMesh();
  TPZGeoEl * elvec[5];
  TPZVec <REAL> coord (3,0.);
  int index;
  //Nodes initialization
  for(i = 0; i < 8; i++){
    for(j=0;j<3;j++){
      coord[j] = Coord[i][j];
    }
    index = gmesh->NodeVec().AllocateNewElement();
    gmesh->NodeVec()[index] = TPZGeoNode(i,coord,*gmesh);
  }

  TPZVec<TPZRefPattern *> refinement_Patterns(6,0);
//  refinement_Patterns.Resize(6);
  refinement_Patterns[0] = new TPZRefPattern("/home/pos/cesar/RefPattern/Tetra_Rib_Side_4.rpt");
  refinement_Patterns[1] = new TPZRefPattern("/home/pos/cesar/RefPattern/Tetra_Rib_Side_5.rpt");
  refinement_Patterns[2] = new TPZRefPattern("/home/pos/cesar/RefPattern/Tetra_Rib_Side_6.rpt");
  refinement_Patterns[3] = new TPZRefPattern("/home/pos/cesar/RefPattern/Tetra_Rib_Side_7.rpt");
  refinement_Patterns[4] = new TPZRefPattern("/home/pos/cesar/RefPattern/Tetra_Rib_Side_8.rpt");
  refinement_Patterns[5] = new TPZRefPattern("/home/pos/cesar/RefPattern/Tetra_Rib_Side_9.rpt");

  for (i=0;i<6;i++) gmesh->InsertRefPattern(refinement_Patterns[i]);

  for (i=0;i<5;i++){
    int ncon = 4;
    TPZVec <int> connect(ncon,0);
    for(j=0; j<ncon;j++){
      connect[j] = Connects[i][j];
    }
    elvec[i] = GeoElementRefPattern(gmesh,7,connect,1,i,refinement_Patterns);
  }
  //Generate neighborhod information
//  gmesh->Print(cout);
  gmesh->BuildConnectivity();
//  gmesh->Print(cout);
  //Create computational mesh
  TPZCompMesh *cmesh = new TPZCompMesh(gmesh);
  TPZMaterial *mat;
  mat = new TPZMaterialTest3D (1);
  cmesh->InsertMaterialObject(mat);

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

/**
 * @brief Undo the encapsulate elements
 */
void TPZCreateApproximationSpace::UndoCondenseLocalEquations(TPZCompMesh &cmesh)
{
    int nel = cmesh.NElements();
    int iel;
    for (iel=0; iel<nel; iel++) {
        TPZCompEl *cel = cmesh.ElementVec()[iel];
        TPZCondensedCompEl *condel = dynamic_cast<TPZCondensedCompEl *>(cel);
        if(!condel) {
            continue;
        }
        condel->Unwrap();
    }    
}

/**
 * @brief Encapsulate the elements in condensed computational elements
 */
void TPZCreateApproximationSpace::CondenseLocalEquations(TPZCompMesh &cmesh)
{
    int nel = cmesh.NElements();
    int iel;
    for (iel=0; iel<nel; iel++) {
        TPZCompEl *cel = cmesh.ElementVec()[iel];
        if(!cel) {
            continue;
        }
        new TPZCondensedCompEl(cel);
    }

}

TPZCompMesh *MeshH1(TPZGeoMesh *gmesh, int pOrder, int dim,bool hasbc)
{
    /// criar materiais
    dim = 2;
    TPZMatPoisson3d *material = new TPZMatPoisson3d( MatId,  dim);
    material->NStateVariables();
	
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    cmesh->SetDimModel(dim);
    TPZMaterial * mat(material);
    cmesh->InsertMaterialObject(mat);
	
    ///Inserir condicao de contorno
    if(hasbc){
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
	
        TPZMaterial * BCond0 = material->CreateBC(mat, bc0,dirichlet, val1, val2);
        TPZMaterial * BCond1 = material->CreateBC(mat, bc1,dirichlet, val1, val2);
        TPZMaterial * BCond2 = material->CreateBC(mat, bc2,dirichlet, val1, val2);
        TPZMaterial * BCond3 = material->CreateBC(mat, bc3,dirichlet, val1, val2);
	
        cmesh->InsertMaterialObject(BCond0);
        cmesh->InsertMaterialObject(BCond1);
        cmesh->InsertMaterialObject(BCond2);
        cmesh->InsertMaterialObject(BCond3);
    }
    
    //solucao exata
//    TPZAutoPointer<TPZFunction<STATE> > solexata;
//    solexata = new TPZDummyFunction<STATE>(EstadoAd);
//    material->SetForcingFunctionExact(solexata);
//
//    //funcao do lado direito da equacao do problema
//    TPZAutoPointer<TPZFunction<STATE> > force;
//    TPZDummyFunction<STATE> *dum;
//    
//    dum = new TPZDummyFunction<STATE>(ForcingOpt);
//    dum->SetPolynomialOrder(20);
//    force = dum;
//    material->SetForcingFunction(force);
	
    
    cmesh->SetDefaultOrder(pOrder);
    cmesh->SetDimModel(dim);
	
	
    cmesh->SetAllCreateFunctionsContinuous();
    
    //Ajuste da estrutura de dados computacional
    cmesh->AutoBuild();
    
    
    return cmesh;
    
}

TPZCompElHDivFull<TSHAPE>::TPZCompElHDivFull(TPZCompMesh &mesh, TPZGeoEl *gel, int64_t &index) :
TPZRegisterClassId(&TPZCompElHDivFull::ClassId),
TPZCompElHDiv<TSHAPE>(mesh,gel,index) {
	
	int i;
	int nconflux= TPZCompElHDiv<TSHAPE>::NConnects();
    this->fConnectIndexes.Resize(nconflux);
	gel->SetReference(this);
	
	for(i=0;i< nconflux;i++)
	{
        int sideaux= i + TSHAPE::NCornerNodes;
		this->fConnectIndexes[i] = this->CreateMidSideConnect(sideaux);
#ifdef LOG4CXX
        if (logger->isDebugEnabled())
        {
            std::stringstream sout;
            sout << "After creating last flux connect " << i << std::endl;
            //	this->Print(sout);
            LOGPZ_DEBUG(logger,sout.str())
        }
#endif
        
		mesh.ConnectVec()[this->fConnectIndexes[i]].IncrementElConnected();
		this->IdentifySideOrder(sideaux);
    }	

void TPZBndCond::Read(TPZStream &buf, void *context)
{
    TPZMaterial::Read(buf, context);
    buf.Read(&fType, 1);
    fBCVal1.Read(buf, 0);
    fBCVal2.Read(buf, 0);
    int MatId;
    buf.Read(&MatId,1);
    TPZCompMesh * pCM = (TPZCompMesh * )/*dynamic_cast<TPZCompMesh *>*/(context);
    fMaterial = pCM->FindMaterial(MatId);
    if(!fMaterial)
    {
        std::cout << " reading a boundary condition without material object!!\n";
#ifdef LOG4CXX
        LOGPZ_FATAL(logger,"reading a boundary condition without material object!!");
#endif
    }
}

TPZCompMesh *CMesh(TPZGeoMesh *gmesh, int pOrder)
{
	const int dim = 2; //dimensao do problema
	const int matId = 1, bc0 = -1, bc1 = -2, bc2=-3, bc3=-4; //MESMOS ids da malha geometrica
    const int dirichlet = 0, neumann = 1;
//    const int mixed = 2; //tipo da condicao de contorno do problema ->default dirichlet na esquerda e na direita
	
    
	///criar malha computacional
	TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
	cmesh->SetDefaultOrder(pOrder);//seta ordem polimonial de aproximacao
	cmesh->SetDimModel(dim);//seta dimensao do modelo
	
    // Criando material
    TPZMatExSimples2D *material = new TPZMatExSimples2D(matId);//criando material que implementa a formulacao fraca do problema modelo
    
	// Inserindo material na malha
	cmesh->InsertMaterialObject(material);
		
	///Inserir condicao de contorno esquerda
	TPZFMatrix<STATE> val1(1,1,0.), val2(1,1,0.);
	TPZMaterial * BCond0 = material->CreateBC(material, bc0, neumann, val1, val2);//cria material que implementa a condicao de contorno da esquerda
	
    cmesh->InsertMaterialObject(BCond0);//insere material na malha
    
	// Condicao de contorno da direita
	TPZMaterial * BCond1 = material->CreateBC(material, bc1, neumann, val1, val2);//cria material que implementa a condicao de contorno da direita
    
    cmesh->InsertMaterialObject(BCond1);//insere material na malha
    
    val2(0,0) = 1.0;//potencial na placa inferior
    // Condicao de contorno da placa inferior
    TPZMaterial * BCond2 = material->CreateBC(material, bc2, dirichlet, val1, val2);//cria material que implementa a condicao de contorno da placa inferior
    
    cmesh->InsertMaterialObject(BCond2);//insere material na malha
    
    val2(0,0) = 1.5;//potencial na placa superior
    // Condicao de contorno da placa superior
    TPZMaterial * BCond3 = material->CreateBC(material, bc3, dirichlet, val1, val2);//cria material que implementa a condicao de contorno da placa superior
	
    cmesh->InsertMaterialObject(BCond3);//insere material na malha
	
	//Cria elementos computacionais que gerenciarao o espaco de aproximacao da malha
	cmesh->AutoBuild();
	
	return cmesh;
	
}

TPZCompMesh  *TPZMGAnalysis::UniformlyRefineMesh(TPZCompMesh *mesh, bool withP) {
	
	TPZGeoMesh *gmesh = mesh->Reference();
	if(!gmesh) {
		cout << "TPZMGAnalysis::UniformlyRefineMesh encountered no geometric mesh\n";
		return 0;
	}
	gmesh->ResetReference();
	TPZCompMesh *cmesh = new TPZCompMesh(gmesh);
	mesh->CopyMaterials(*cmesh);
	TPZAdmChunkVector<TPZCompEl *> &elementvec = mesh->ElementVec();
	int el,nelem = elementvec.NElements();
	for(el=0; el<nelem; el++) {
		TPZCompEl *cel = elementvec[el];
		if(!cel) continue;
		TPZInterpolatedElement *cint = dynamic_cast<TPZInterpolatedElement *> (cel);
		if(!cint) {
			cout << "TPZMGAnalysis::UniformlyRefineMesh encountered a non interpolated element\n";
			continue;
		}
		int ncon = cint->NConnects();
		int porder = cint->PreferredSideOrder(ncon-1);
		
		TPZGeoEl *gel = cint->Reference();
		if(!gel) {
			cout << "TPZMGAnalysis::UniformlyRefineMesh encountered an element without geometric reference\n";
			continue;
		}
		TPZVec<TPZGeoEl *> sub;
		gel->Divide(sub);
		int nsub = sub.NElements();
		int isub;
		int celindex;
		for(isub=0; isub<nsub; isub++) {
			TPZInterpolatedElement *csint;
			csint = (TPZInterpolatedElement *) cmesh->CreateCompEl(sub[isub],celindex);
			if(withP) csint->PRefine(porder+1);
			else csint->PRefine(porder);
		}
	}
	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();
}

TPZCompMesh *L2ProjectionP(TPZGeoMesh *gmesh, int pOrder, TPZVec<STATE> &solini)
{
    /// criar materiais
    int dim = 2;
    TPZL2Projection *material;
    material = new TPZL2Projection(1, dim, 1, solini, pOrder);
    
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    cmesh->SetDimModel(dim);
    TPZMaterial * mat(material);
    cmesh->InsertMaterialObject(mat);
    TPZAutoPointer<TPZFunction<STATE> > forcef = new TPZDummyFunction<STATE>(InitialPressure);
    material->SetForcingFunction(forcef);
    cmesh->SetAllCreateFunctionsContinuous();
    cmesh->SetDefaultOrder(pOrder);
    cmesh->SetDimModel(dim);
    cmesh->AutoBuild();
    
    
    return cmesh;
    
}