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

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

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

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 *MeshL2(TPZGeoMesh *gmesh, int pOrder, int dim)
{
    /// 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
//    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>(OptForcing);
//    dum->SetPolynomialOrder(20);
//    force = dum;
//    material->SetForcingFunction(force);
    
    
    cmesh->SetDefaultOrder(pOrder);
    cmesh->SetDimModel(dim);
    
    
    cmesh->SetAllCreateFunctionsDiscontinuous();
    
    //cmesh->SetAllCreateFunctionsContinuous();
    
    //Ajuste da estrutura de dados computacional
    cmesh->AutoBuild();
    
    
    return cmesh;
    
}

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

TPZCompMesh * CMeshFooting2D(TPZGeoMesh * gmesh, int p_order){
    
    unsigned int dim  = gmesh->Dimension();
    const std::string name("ElastoPlastic Footing Problem ");

    // Setting up attributes
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    cmesh->SetName(name);
    cmesh->SetDefaultOrder(p_order);
    cmesh->SetDimModel(dim);
    
    // Mohr Coulomb data
    REAL mc_cohesion    = 10.0;
    REAL mc_phi         = (20.0*M_PI/180);
    REAL mc_psi         = mc_phi;
    
    /// ElastoPlastic Material using Mohr Coulomb
    // Elastic predictor
    TPZElasticResponse ER;
    REAL G = 400*mc_cohesion;
    REAL nu = 0.3;
    REAL E = 2.0*G*(1+nu);
    

    TPZPlasticStepPV<TPZYCMohrCoulombPV, TPZElasticResponse> LEMC;
    ER.SetEngineeringData(E, nu);
    LEMC.SetElasticResponse(ER);
    LEMC.fYC.SetUp(mc_phi, mc_psi, mc_cohesion, ER);
    int PlaneStrain = 1;
    
//    TPZElasticCriterion MatEla;
//    MatEla.SetElasticResponse(ER);
//    TPZMatElastoPlastic2D < TPZElasticCriterion, TPZElastoPlasticMem > * material = new TPZMatElastoPlastic2D < TPZElasticCriterion, TPZElastoPlasticMem >(ERock,PlaneStrain);
//    material->SetPlasticityModel(MatEla);
//    cmesh->InsertMaterialObject(material);
    
    TPZMatElastoPlastic2D < TPZPlasticStepPV<TPZYCMohrCoulombPV, TPZElasticResponse>, TPZElastoPlasticMem > * material = new TPZMatElastoPlastic2D < TPZPlasticStepPV<TPZYCMohrCoulombPV, TPZElasticResponse>, TPZElastoPlasticMem >(ERock,PlaneStrain);
    material->SetPlasticityModel(LEMC);
    cmesh->InsertMaterialObject(material);
    
    
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    
    val2(0,0) = 0;
    val2(1,0) = 1;
    TPZBndCond * bc_bottom = material->CreateBC(material, EBottomBC, Eu_null, val1, val2);
    
    val2(0,0) = 1;
    val2(1,0) = 0;
    TPZBndCond * bc_lateral = material->CreateBC(material, ELateralBC, Eu_null, val1, val2);
    
//    val2(0,0) = 0;
//    val2(1,0) = 0;
//    val1(0,0) = 0;
//    val1(1,1) = 1;
    val2(0,0) = 0;
    val2(1,0) = 0;
    TPZBndCond * bc_top = material->CreateBC(material, ETopBC, ETn, val1, val2);
    
    val2(0,0) = 0;
    val2(1,0) = 0;
    TPZBndCond * bc_top_null = material->CreateBC(material, ETopNullBC, ETn, val1, val2);
    
    cmesh->InsertMaterialObject(bc_bottom);
    cmesh->InsertMaterialObject(bc_lateral);
    cmesh->InsertMaterialObject(bc_top);
//    cmesh->InsertMaterialObject(bc_top_null);
    
    cmesh->SetAllCreateFunctionsContinuousWithMem();
    cmesh->AutoBuild();
    
#ifdef PZDEBUG
    std::ofstream out("cmesh.txt");
    cmesh->Print(out);
#endif
    
    return cmesh;
}

TPZCompMesh *MalhaCompMultifisica(TPZGeoMesh * gmesh,TPZVec<TPZCompMesh *> meshvec, TPZMatUncoupledPoissonDisc* &mymaterial){
    
    
    // Creating computational mesh for multiphysic elements
    gmesh->ResetReference();
    TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
    
    mphysics->SetAllCreateFunctionsMultiphysicElem();
    
    int dim = 2;
    mphysics->SetDimModel(dim);
    mymaterial = new TPZMatUncoupledPoissonDisc(matId, mphysics->Dimension());
    
    mymaterial->SetParameters(1., 1.);
    mymaterial->SetInternalFlux(-8.,0.);
    //mymaterial->SetInternalFlux(0.,0.);
    
    mymaterial->SetNonSymmetricOne();
    mymaterial->SetNonSymmetricTwo();
    mymaterial->SetPenaltyConstant(0., 0.);
    
    TPZMaterial * mat(mymaterial);
    mphysics->InsertMaterialObject(mat);
    
    TPZAutoPointer<TPZFunction<STATE> > forcef = new TPZDummyFunction<STATE>(ForcingF, 5);
    
    //
    //    TPZAutoPointer<TPZFunction<STATE> > forcef = new TPZDummyFunction<STATE>(ForcingF);
    mymaterial->SetForcingFunction(forcef);
    
    
    ///Inserir condicao de contorno
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    
    TPZMaterial * BCond0 = mymaterial->CreateBC(mat, bc0,neumann_dirichlet, val1, val2);
    TPZMaterial * BCond2 = mymaterial->CreateBC(mat, bc2,neumann_dirichlet, val1, val2);
    TPZMaterial * BCond1 = mymaterial->CreateBC(mat, bc1,dirichlet, val1, val2);
    TPZMaterial * BCond3 = mymaterial->CreateBC(mat, bc3,dirichlet, val1, val2);
    
    //    TPZMaterial * BCond0 = mymaterial->CreateBC(mat, bc0,dirichlet, val1, val2);
    //    TPZMaterial * BCond2 = mymaterial->CreateBC(mat, bc2,dirichlet, val1, val2);
    //    TPZMaterial * BCond1 = mymaterial->CreateBC(mat, bc1,dirichlet, val1, val2);
    //    TPZMaterial * BCond3 = mymaterial->CreateBC(mat, bc3,dirichlet, val1, val2);
    
    
    mphysics->InsertMaterialObject(BCond0);
    mphysics->InsertMaterialObject(BCond1);
    mphysics->InsertMaterialObject(BCond2);
    mphysics->InsertMaterialObject(BCond3);
    
    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);
    
    mphysics->Reference()->ResetReference();
    mphysics->LoadReferences();
    
    if (disc_functions==true){
        //criar elementos de interface
        int nel = mphysics->ElementVec().NElements();
        for(int el = 0; el < nel; el++)
        {
            TPZCompEl * compEl = mphysics->ElementVec()[el];
            if(!compEl) continue;
            int index = compEl ->Index();
            if(compEl->Dimension() == mphysics->Dimension())
            {
                TPZMultiphysicsElement * InterpEl = dynamic_cast<TPZMultiphysicsElement *>(mphysics->ElementVec()[index]);
                if(!InterpEl) continue;
                InterpEl->CreateInterfaces();
                
            }
        }
    }
    
    return mphysics;
    
}

//.........这里部分代码省略.........
   TPZBndCond *bc;
   REAL big = 1.e12;
   TPZFMatrix val1(6,6,0.),val2(6,1,0.);

   // engastes nos lados 4 e 7 do elemento 0
   TPZGeoElBC(elg[0],4,-2,*firstmesh);
   TPZGeoElBC(elg[0],7,-2,*firstmesh);

   // engaste no lado 4 do elemento 1
   TPZGeoElBC(elg[1],4,-2,*firstmesh);

   // engaste no lado 4 do elemento 2
   TPZGeoElBC(elg[2],4,-2,*firstmesh);

   // engaste no lado 4 do elemento 3
   TPZGeoElBC(elg[3],4,-2,*firstmesh);

   // imposicao do valor zero associado a condicao -2 (engaste)
   bc = pl->CreateBC(-2,0,val1,val2);
   secondmesh->InsertMaterialObject(bc);

   // imposicao da condicao de simetria no lado 5 do elemento 4
   val1(0,0)=big;
   val1(1,1)=0.;
   val1(2,2)=0.;
   val1(3,3)=0.;
   val1(4,4)=big;
   val1(5,5)=big;
   TPZGeoElBC(elg[3],5,-3,*firstmesh);
   bc = pl->CreateBC(-3,2,val1,val2);
   secondmesh->InsertMaterialObject(bc);

   //ordem de interpolacao
   int ord;
   cout << "Entre ordem 1,2,3,4,5 : ";
   cin >> ord;
//   TPZCompEl::gOrder = ord;
   firstmesh.SetDefaultOrder(order);
   //construção malha computacional
   TPZVec<int> csub(0);
   TPZManVector<TPZGeoEl *> pv(4);
   int n1=1,level=0;
   cout << "\nDividir ate nivel ? ";
   int resp;
   cin >> resp;
   int nelc = firstmesh->ElementVec().NElements();
   int el;
   TPZGeoEl *cpel;
   for(el=0;el<firstmesh->ElementVec().NElements();el++) {
     cpel = firstmesh->ElementVec()[el];
     if(cpel && cpel->Level() < resp)
		cpel->Divide(pv);

   }
   cout << "\nDividir o elemento esquerdo superior quantas vezes? ";
   cin >> resp;
   cpel = firstmesh->ElementVec()[0];
   for(el=0; el<resp; el++) {
		cpel->Divide(pv);
		cpel = pv[3];
   }
   //analysis
   secondmesh->AutoBuild();
   firstmesh->Print(outgm1);
   outgm1.flush();
   secondmesh->AdjustBoundaryElements();
   secondmesh->InitializeBlock();
   secondmesh->Print(outcm1);
   TPZAnalysis an(secondmesh,outcm1);
   int numeq = secondmesh->NEquations();
   secondmesh->Print(outcm1);
   outcm1.flush();
   TPZVec<int> skyline;
   secondmesh->Skyline(skyline);
   TPZSkylMatrix *stiff = new TPZSkylMatrix(numeq,skyline);
   an.SetMatrix(stiff);
   an.Solver().SetDirect(ECholesky);
   secondmesh->SetName("Malha Computacional :  Connects e Elementos");
   // Posprocessamento
   an.Run(outcm2);
   TPZVec<char *> scalnames(5);
   scalnames[0] = "Mn1";
   scalnames[1] = "Mn2";
   scalnames[2] = "Vn1";
   scalnames[3] = "Vn2";
   scalnames[4] = "Deslocz";
   TPZVec<char *> vecnames(0);
   char plotfile[] =  "placaPos.pos";
   char pltfile[] =  "placaView.plt";
   an.DefineGraphMesh(2, scalnames, vecnames, plotfile);
   an.Print("FEM SOLUTION ",outcm1);
   an.PostProcess(3);
   an.DefineGraphMesh(2, scalnames, vecnames, pltfile);
   an.PostProcess(2);
   firstmesh->Print(outgm1);
   outgm1.flush();
   delete secondmesh;
   delete firstmesh;
   return 0;
}

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

int main() {

   //malha geometrica
   TPZGeoMesh *firstmesh = new TPZGeoMesh;
   firstmesh->NodeVec().Resize(3);
   TPZVec<REAL> coord(2);
   coord[0] = 0.;
   coord[1] = 0.;
   //nos geometricos
   firstmesh->NodeVec()[0].Initialize(coord,*firstmesh);
   coord[0] = 1.0;
   firstmesh->NodeVec()[1].Initialize(coord,*firstmesh);
   coord[1] = 1.0;
   firstmesh->NodeVec()[2].Initialize(coord,*firstmesh);
//   coord[0] = 0.0;
//   firstmesh->NodeVec()[3].Initialize(coord,*firstmesh);
   TPZVec<int> nodeindexes(3);//triangulo
   nodeindexes[0] = 0;//local[i] = global[i] , i=0,1,2,3
   nodeindexes[1] = 1;
   nodeindexes[2] = 2;
   //elementos geometricos
   TPZGeoElT2d *elq1 = new TPZGeoElT2d(nodeindexes,1,*firstmesh);
 //orientacao local de um segundo elemento superposto
   int i,sen;;
   cout<<"Sentido local antihorario/horario : 0/1 ?  ";
   cin>>sen;
   cout<<"Entre primeiro no = 0,1,2 : ";
   cin>>i;
   if(sen==0) {//direito
        nodeindexes[0] = (0+i)%3;//local[i] = global[j] , i,j em {0,1,2}
        nodeindexes[1] = (1+i)%3;
        nodeindexes[2] = (2+i)%3;
	} else {//inverso
        nodeindexes[0] = (0+i)%3;//local[i] = global[j] , i,j em {0,1,2}
        nodeindexes[1] = (2+i)%3;
        nodeindexes[2] = (1+i)%3;
   }
/*   nodeindexes[0] = 1;//local[i] = global[i] , i=0,1,2,3
   nodeindexes[1] = 2;
   nodeindexes[2] = 3;*/
   TPZGeoElT2d *elq2 = new TPZGeoElT2d(nodeindexes,1,*firstmesh);//segundo elemento superposto ao primeiro
/*   coord[1] = 0.0;
   coord[0] = 2.0;
   firstmesh->NodeVec()[4].Initialize(coord,*firstmesh);
   coord[1] = 1.0;
   firstmesh->NodeVec()[5].Initialize(coord,*firstmesh);
   nodeindexes[0] = 1;//local[i] = global[i] , i=0,1,2,3
   nodeindexes[1] = 4;
   nodeindexes[2] = 5;
   nodeindexes[3] = 2;
   TPZGeoElT2d *elq2 = new TPZGeoElT2d(nodeindexes,1,*firstmesh);    */
   //Arquivos de saida
	ofstream outgm1("outgm1.dat");
   ofstream outcm1("outcm1.dat");
	ofstream outcm2("outcm2.dat");
   //montagem de conectividades entre elementos
   firstmesh->BuildConnectivity();
 	firstmesh->Print(outgm1);
   outgm1.flush();
  	//teste de divisao geometrica : 1 elemento
   TPZVec<TPZGeoEl *> vecsub,vecsub1;
   elq1->Divide(vecsub);//divide 0
   elq2->Divide(vecsub);//divide 1
/*   vecsub[2]->Divide(vecsub1);//
   vecsub1[3]->Divide(vecsub1);
	vecsub[0]->Divide(vecsub1);//divide 1
   vecsub1[2]->Divide(vecsub1); */
 	firstmesh->Print(outgm1);
   outgm1.flush();
   //malha computacional
   TPZCompMesh *secondmesh = new TPZCompMesh(firstmesh);
   //material
   int matindex = secondmesh->MaterialVec().AllocateNewElement();
   TPZFMatrix k(1,1,1.),f(1,1,0.),c(1,2,1.);
   TPZMat2dLin * mat = new TPZMat2dLin(1);
   mat->SetMaterial(k,c,f);
   //mat->SetForcingFunction(force);
   mat->SetForcingFunction(derivforce);
   secondmesh->MaterialVec()[matindex] = mat;
   //CC : condicao de contorno
   //ordem de interpolacao
//   TPZCompEl::gOrder = 3;
   cmesh.SetDefaultOrder(3);
   //constroe a malha computacional
   secondmesh->AutoBuild();
   secondmesh->InitializeBlock();
   secondmesh->ComputeConnectSequence();
   secondmesh->Print(outcm1);
   outcm1.flush();
	//Resolucao do sistema
   TPZFMatrix Rhs(secondmesh->NEquations(),1),Stiff(secondmesh->NEquations(),secondmesh->NEquations()),U;
   Stiff.Zero();
   Rhs.Zero();
   secondmesh->Assemble(Stiff,Rhs);
   Rhs.Print("Rhs teste",outcm2);
   Stiff.Print("Bloco teste",outcm2);
	Rhs.Print("Computational Mesh -> fBlock",outcm2);
   TPZMatrixSolver solver(&Stiff);
   solver.SetDirect(ELU);
   solver.Solve(Rhs,U);
//.........这里部分代码省略.........

TPZCompMesh *ComputationalElasticityMesh2D(TPZAutoPointer<TPZGeoMesh>  gmesh,int pOrder)
{
    
    // remove some connectivities 3, 5
    TPZGeoEl *gel = gmesh->Element(0);
    TPZGeoElSide gelside(gel,3);
    gelside.RemoveConnectivity();
    gelside.SetSide(5);
    gelside.RemoveConnectivity();
    gelside.SetSide(4);
    TPZGeoElSide neighbour = gelside.NNeighbours();
    int matid = neighbour.Element()->MaterialId();
    gel->SetMaterialId(matid);
    neighbour.Element()->RemoveConnectivities();
    int64_t index = neighbour.Element()->Index();
    delete neighbour.Element();
    gmesh->ElementVec()[index] = 0;
    
    
    // Plane strain assumption
    int planestress = 0;
    
    // Getting mesh dimension
    int dim = 2;
    
    TPZMatElasticity2D *materialConcrete;
    materialConcrete = new TPZMatElasticity2D(EMatConcrete);
    
    TPZMatElasticity2D *materialSteel;
    materialSteel = new TPZMatElasticity2D(EMatSteel);
    
    
    // Setting up paremeters
    materialConcrete->SetfPlaneProblem(planestress);
    materialConcrete->SetElasticity(25.e6, 0.25);
    materialSteel->SetElasticity(205.e6, 0.25);
    
    //material->SetBiotAlpha(Alpha);cade o metodo?
    
    
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    cmesh->SetDefaultOrder(pOrder);
    cmesh->SetDimModel(dim);
    
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    
    val2(0,0) = 0.0;
    val2(1,0) = 0.0;
    val1(1,1) = 1.e12;
    TPZMaterial * BCond2 = materialConcrete->CreateBC(materialConcrete,EBottom,3, val1, val2);
    
    val2(0,0) = 0.0;
    val2(1,0) = 0.0;
    val1.Zero();
    val1(0,0) = 1.e12;
    TPZMaterial * BCond3 = materialConcrete->CreateBC(materialConcrete,ELateral,3, val1, val2);
    
    val2(0,0) = 0.0;
    val2(1,0) = -1000.0;
    val1.Zero();
    TPZMaterial * BCond4 = materialSteel->CreateBC(materialSteel,EBeam,1, val1, val2);
    
    cmesh->SetAllCreateFunctionsContinuous();
    cmesh->InsertMaterialObject(materialConcrete);
    cmesh->InsertMaterialObject(materialSteel);
    cmesh->InsertMaterialObject(BCond2);
    cmesh->InsertMaterialObject(BCond3);
    cmesh->InsertMaterialObject(BCond4);
    cmesh->AutoBuild();
    return cmesh;
    
}

//.........这里部分代码省略.........
   REAL big = 1.e12;
   TPZFMatrix val1(6,6,0.),val2(6,1,0.);

   val1(0,0)=big;
   val1(1,1)=big;
   val1(2,2)=big;
   val1(3,3)=0.;
   val1(4,4)=0.;
   val1(5,5)=0.;

   TPZGeoElBC(elg0,5,-2,*firstmesh);
   bc = pl->CreateBC(-2,2,val1,val2);
   secondmesh->InsertMaterialObject(bc);


   TPZGeoElBC(elg0,6,-3,*firstmesh);
   bc = pl->CreateBC(-3,2,val1,val2);
   secondmesh->InsertMaterialObject(bc);

   val1(0,0)=0.;
   val1(1,1)=big;
   val1(2,2)=0.;
   val1(3,3)=big;
   val1(4,4)=0.;
   val1(5,5)=0.;

	TPZGeoElBC(elg0,4,-1,*firstmesh);
   bc = pl->CreateBC(-1,2,val1,val2);
   secondmesh->InsertMaterialObject(bc);

   val1(0,0)=big;
   val1(1,1)=0.;
   val1(2,2)=0.;
   val1(3,3)=0.;
   val1(4,4)=big;
   val1(5,5)=0.;

	TPZGeoElBC(elg0,7,-4,*firstmesh);
   bc = pl->CreateBC(-4,2,val1,val2);
   secondmesh->InsertMaterialObject(bc);

   //ordem de interpolacao
   int ord;
   cout << "Entre ordem 1,2,3,4,5 : ";
   cin >> ord;
//   TPZCompEl::gOrder = ord;
   cmesh.SetDefaultOrder(ord);
   //construção malha computacional
   TPZVec<int> csub(0);
   TPZManVector<TPZGeoEl *> pv(4);
   int n1=1,level=0;
   cout << "\nDividir ate nivel ? ";
   int resp;
   cin >> resp;
   int nelc = firstmesh->ElementVec().NElements();
   int el;
   TPZGeoEl *cpel;
   for(el=0;el<firstmesh->ElementVec().NElements();el++) {
     cpel = firstmesh->ElementVec()[el];
     if(cpel && cpel->Level() < resp)
		cpel->Divide(pv);

   }
   //analysis
   secondmesh->AutoBuild();
   secondmesh->AdjustBoundaryElements();
   secondmesh->InitializeBlock();
   secondmesh->Print(outcm1);
   TPZAnalysis an(secondmesh,outcm1);
   int numeq = secondmesh->NEquations();
   secondmesh->Print(outcm1);
   outcm1.flush();
   TPZVec<int> skyline;
   secondmesh->Skyline(skyline);
   TPZSkylMatrix *stiff = new TPZSkylMatrix(numeq,skyline);
   an.SetMatrix(stiff);
   an.Solver().SetDirect(ECholesky);
   secondmesh->SetName("Malha Computacional :  Connects e Elementos");
   // Posprocessamento
   an.Run(outcm2);
   TPZVec<char *> scalnames(5);
   scalnames[0] = "Mn1";
   scalnames[1] = "Mn2";
   scalnames[2] = "Sign1";
   scalnames[3] = "Sign2";
   scalnames[4] = "Deslocz";
   TPZVec<char *> vecnames(0);
   char plotfile[] =  "placaPos.pos";
   char pltfile[] =  "placaView.plt";
   an.DefineGraphMesh(2, scalnames, vecnames, plotfile);
   an.Print("FEM SOLUTION ",outcm1);
   an.PostProcess(2);
   an.DefineGraphMesh(2, scalnames, vecnames, pltfile);
   an.PostProcess(2);
   firstmesh->Print(outgm1);
   outgm1.flush();
   delete secondmesh;
   delete firstmesh;
   return 0;
}

/** @brief Creates the computational elements, and the degree of freedom nodes */
void TPZCreateApproximationSpace::AutoBuild(TPZCompMesh &cmesh){
    cmesh.AutoBuild();
}