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

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 * 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 *MalhaMultifisicaOpt(TPZVec<TPZCompMesh *> meshvec, TPZGeoMesh *gmesh){
    
    //Creating computational mesh for multiphysic elements
    gmesh->ResetReference();
    TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
    
    //criando material
    int dim =2;
    
    TPZMatPoissonControl *material = new TPZMatPoissonControl(MatId,dim);
    
    
    //incluindo os dados do problema
    REAL k=1;
    REAL alpha=1;
    material-> SetParameters( k, alpha);
    
    //solucao exata
    TPZAutoPointer<TPZFunction<STATE> > solexata;
    solexata = new TPZDummyFunction<STATE>(StateAd, 5);
    material->SetForcingFunctionExact(solexata);
    
    //funcao do lado direito da equacao do problema
    TPZAutoPointer<TPZFunction<STATE> > force;
    TPZDummyFunction<STATE> *dum;
    dum = new TPZDummyFunction<STATE>(OptForcing, 5);
    dum->SetPolynomialOrder(20);
    force = dum;
    material->SetForcingFunction(force);
    
    //inserindo o material na malha computacional
    TPZMaterial *mat(material);
    mphysics->InsertMaterialObject(mat);
    mphysics->SetDimModel(dim);
    
    //Criando condicoes de contorno
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    
    
    TPZMaterial * BCond0 = material->CreateBC(mat, bc0, bcdirichlet, val1, val2);
    TPZMaterial * BCond1 = material->CreateBC(mat, bc1, bcdirichlet, val1, val2);
    TPZMaterial * BCond2 = material->CreateBC(mat, bc2, bcdirichlet, val1, val2);
    TPZMaterial * BCond3 = material->CreateBC(mat, bc3, bcdirichlet, val1, val2);
    
    ///Inserir condicoes de contorno
    mphysics->InsertMaterialObject(BCond0);
    mphysics->InsertMaterialObject(BCond1);
    mphysics->InsertMaterialObject(BCond2);
    mphysics->InsertMaterialObject(BCond3);
    
    mphysics->SetAllCreateFunctionsMultiphysicElem();
    
    //Fazendo auto build
    mphysics->AutoBuild();
    mphysics->AdjustBoundaryElements();
    mphysics->CleanUpUnconnectedNodes();
    
    TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
    TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
    TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
    
    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;
}

//malha multifisica para o metodo da dupla projecao
TPZCompMesh *MalhaMDP(TPZVec<TPZCompMesh *> meshvec,TPZGeoMesh * gmesh){
    
    //Creating computational mesh for multiphysic elements
	gmesh->ResetReference();
	TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
    
    //criando material
    int dim =2;
    TPZMDPMaterial *material = new TPZMDPMaterial(1,dim);
    
    //incluindo os dados do problema
    REAL coefk = 1.;
    material->SetParameters(coefk, 0.);
    
    //solucao exata
    TPZAutoPointer<TPZFunction<STATE> > solexata;
    solexata = new TPZDummyFunction<STATE>(SolSuave);
    material->SetForcingFunctionExact(solexata);
    
    //funcao do lado direito da equacao do problema
    TPZAutoPointer<TPZFunction<STATE> > force;
    TPZDummyFunction<STATE> *dum;
    dum = new TPZDummyFunction<STATE>(ForceSuave);
    dum->SetPolynomialOrder(20);
    force = dum;
    material->SetForcingFunction(force);
    
    //inserindo o material na malha computacional
    TPZMaterial *mat(material);
    mphysics->InsertMaterialObject(mat);
    
    //Criando condicoes de contorno
    TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
    int boundcond = dirichlet;
    //BC -1
    TPZMaterial * BCondD1 = material->CreateBC(mat, bc1,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet1 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD1->SetForcingFunction(bcmatDirichlet1);
    mphysics->InsertMaterialObject(BCondD1);
    
    //BC -2
	TPZMaterial * BCondD2 = material->CreateBC(mat, bc2,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet2 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD2->SetForcingFunction(bcmatDirichlet2);
    mphysics->InsertMaterialObject(BCondD2);
    
    //BC -3
	TPZMaterial * BCondD3 = material->CreateBC(mat, bc3,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet3 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD3->SetForcingFunction(bcmatDirichlet3);
    mphysics->InsertMaterialObject(BCondD3);
    
    //BC -4
	TPZMaterial * BCondD4 = material->CreateBC(mat, bc4,boundcond, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet4 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD4->SetForcingFunction(bcmatDirichlet4);
    mphysics->InsertMaterialObject(BCondD4);

    
    mphysics->InsertMaterialObject(BCondD1);
    mphysics->InsertMaterialObject(BCondD2);
    mphysics->InsertMaterialObject(BCondD3);
    mphysics->InsertMaterialObject(BCondD4);
    
    
    //set multiphysics element
    mphysics->SetDimModel(dim);
    mphysics->SetAllCreateFunctionsMultiphysicElem();
    
    //Fazendo auto build
    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);
    
    return mphysics;
}

TPZCompMesh * CompMesh(TPZGeoMesh *gmesh, int porder)
{
    /// criar materiais
	int dim = gmesh->Dimension();
    
    TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
    
    TPZMatLaplacian *material = new TPZMatLaplacian(1,dim);
    
//    TPZAutoPointer<TPZFunction<REAL> > forcef = new TPZDummyFunction<REAL>(ForceSuave);
//    material->SetForcingFunction(forcef);
    
    TPZAutoPointer<TPZFunction<STATE> > force;
    TPZDummyFunction<STATE> *dum;
    dum = new TPZDummyFunction<STATE>(ForceSuave);
    dum->SetPolynomialOrder(20);
    force = dum;
    material->SetForcingFunction(force);
    
    
    TPZAutoPointer<TPZFunction<STATE> > solExata= new TPZDummyFunction<STATE>(SolSuave);
    material->SetForcingFunctionExact(solExata);
    
    TPZMaterial * mat(material);
    cmesh->InsertMaterialObject(mat);
    
    cmesh->SetDimModel(dim);
    cmesh->SetDefaultOrder(porder);

	///Inserir condicao de contorno
	TPZFMatrix<STATE> val1(2,2,1.), val2(2,1,0.);
	
    //BC -1
    TPZMaterial * BCondD1 = material->CreateBC(mat, bc1,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet1 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD1->SetForcingFunction(bcmatDirichlet1);
    cmesh->InsertMaterialObject(BCondD1);
    
    //BC -2
	TPZMaterial * BCondD2 = material->CreateBC(mat, bc2,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet2 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD2->SetForcingFunction(bcmatDirichlet2);
    cmesh->InsertMaterialObject(BCondD2);
    
    //BC -3
	TPZMaterial * BCondD3 = material->CreateBC(mat, bc3,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet3 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD3->SetForcingFunction(bcmatDirichlet3);
    cmesh->InsertMaterialObject(BCondD3);
    
    //BC -4
	TPZMaterial * BCondD4 = material->CreateBC(mat, bc4,dirichlet, val1, val2);
    TPZAutoPointer<TPZFunction<REAL> > bcmatDirichlet4 = new TPZDummyFunction<REAL>(DirichletSuave);
    BCondD4->SetForcingFunction(bcmatDirichlet4);
    cmesh->InsertMaterialObject(BCondD4);
    
    
    //Fazendo auto build
    cmesh->SetAllCreateFunctionsContinuous();
    cmesh->AutoBuild();
	cmesh->AdjustBoundaryElements();
	cmesh->CleanUpUnconnectedNodes();
    
    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;
}

TPZCompMesh *TPZAdaptMesh::CreateCompMesh (TPZCompMesh *mesh,                                          //malha a refinar
                                           TPZVec<TPZGeoEl *> &gelstack,   //
                                           TPZVec<int> &porders) {
    
    //Cria um ponteiro para a malha geom�trica de mesh
    TPZGeoMesh *gmesh = mesh->Reference();
    if(!gmesh) {
        cout << "TPZAdaptMesh::CreateCompMesh encountered no geometric mesh\n";
        return 0;
    }
    
    //Reseta as refer�ncias do ponteiro para a malha geom�trica criada
    //e cria uma nova malha computacional baseada nesta malha geom�trica
    gmesh->ResetReference();
    TPZCompMesh *cmesh = new TPZCompMesh(gmesh);
//    int nmat = mesh->MaterialVec().size();
  //  int m;
    
    //Cria um clone do vetor de materiais da malha mesh
    mesh->CopyMaterials(*cmesh);
    /*  for(m=0; m<nmat; m++) {
     TPZMaterial * mat = mesh->MaterialVec()[m];
     if(!mat) continue;
     mat->Clone(cmesh->MaterialVec());
     }
     */
    //Idenifica o vetor de elementos computacionais de mesh
    //  TPZAdmChunkVector<TPZCompEl *> &elementvec = mesh->ElementVec();
    
    int el,nelem = gelstack.NElements();
    //  cmesh->SetName("Antes PRefine");
    //  cmesh->Print(cout);
    for(el=0; el<nelem; el++) {
        
        //identifica os elementos geom�tricos passados em gelstack
        TPZGeoEl *gel = gelstack[el];
        if(!gel) {
            cout << "TPZAdaptMesh::CreateCompMesh encountered an null element\n";
            continue;
        }
        long celindex;
        
        //Cria um TPZIntel baseado no gel identificado
        TPZInterpolatedElement *csint;
        csint = dynamic_cast<TPZInterpolatedElement *> (cmesh->CreateCompEl(gel,celindex));
        if(!csint) continue;
        
        //Refina em p o elemento criado
        //	cmesh->SetName("depois criar elemento");
        //	cmesh->Print(cout);
        
        csint->PRefine(porders[el]);
        //	cmesh->SetName("depois prefine no elemento");
        //	cmesh->Print(cout);
    }
#ifndef CLONEBCTOO
    nelem = gmesh->NElements();
    for (el=0; el<nelem; el++) {
        TPZGeoEl *gel = gmesh->ElementVec()[el];
        if (!gel || gel->Reference()) {
            continue;
        }
        int matid = gel->MaterialId();
        if (matid < 0) {
            TPZStack<TPZCompElSide> celstack;
            int ns = gel->NSides();
            TPZGeoElSide gelside(gel,ns-1);
            gelside.HigherLevelCompElementList2(celstack, 1, 1);
            if (celstack.size()) {
                TPZStack<TPZGeoEl *> subels;
                gel->Divide(subels);
            }
        }
    }
    nelem = gmesh->NElements();
    for (el=0; el<nelem; el++) {
        TPZGeoEl *gel = gmesh->ElementVec()[el];
        if (!gel || gel->Reference()) {
            continue;
        }
        int matid = gel->MaterialId();
        if (matid < 0) {
            TPZStack<TPZCompElSide> celstack;
            int ns = gel->NSides();
            TPZGeoElSide gelside(gel,ns-1);
            gelside.EqualLevelCompElementList(celstack, 1, 0);
            if (celstack.size()) {
                long index;
                cmesh->CreateCompEl(gel, index);
            }
        }
    }
#endif
    //Mais einh!!
    //	cmesh->SetName("Antes Adjust");
    //	cmesh->Print(cout);
    cmesh->AdjustBoundaryElements();
    //  cmesh->SetName("Depois");
    //  cmesh->Print(cout);
    return cmesh;
//.........这里部分代码省略.........