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

int main2(int argc, char *argv[])
{
#ifdef LOG4CXX
    InitializePZLOG();
#endif
    
    REAL Lx = 1.;
    REAL Ly = 1.;
    
    {// Checando funcoes
        TPZVec<REAL> pto(2,0.);
        TPZVec<STATE> disp(1,0.),solp(1,0.);
        pto[0]=0.200935;
        pto[1]=0.598129;
        TPZFMatrix<STATE>flux(3,1,0.);
        SolExata(pto, solp,flux);
        Forcing(pto, disp);
        
        solp.Print(std::cout);
        flux.Print(std::cout);
        disp.Print(std::cout);
    }
    
   // ofstream saidaerro("../ErroPoissonHdivMalhaQuad.txt",ios::app);
    ofstream saidaerro("../ErroPoissonHdivMalhaTriang.txt",ios::app);

    for(int p = 1; p<2; p++)
    {
        int pq = p;
        int pp;
        if(ftriang==true){
            pp = pq-1;
        }else{
            pp = pq;
        }
        
        int ndiv;
        saidaerro<<"\n CALCULO DO ERRO, COM ORDEM POLINOMIAL pq = " << pq << " e pp = "<< pp <<endl;
        for (ndiv = 1; ndiv< 2; ndiv++)
        {
            
            //std::cout << "p order " << p << " number of divisions " << ndiv << std::endl;
            
            saidaerro<<"\n<<<<<< Numero de divisoes uniforme ndiv = " << ndiv <<" >>>>>>>>>>> "<<endl;
            
            TPZGeoMesh *gmesh = GMesh(ftriang, Lx, Ly);
            ofstream arg("gmesh1.txt");
            gmesh->Print(arg);
            
            UniformRefine(gmesh, ndiv);
            
            TPZCompMesh *cmesh1 = CMeshFlux(gmesh, pq);
            TPZCompMesh *cmesh2 = CMeshPressure(gmesh, pp);
            
//            ofstream arg1("cmeshflux.txt");
//            cmesh1->Print(arg1);
//            
//            ofstream arg2("cmeshpressure.txt");
//            cmesh2->Print(arg2);
//            
//            ofstream arg4("gmesh2.txt");
//            gmesh->Print(arg4);
            
            
            //malha multifisica
            TPZVec<TPZCompMesh *> meshvec(2);
            meshvec[0] = cmesh1;
            meshvec[1] = cmesh2;
            
            TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
            
            std::cout << "Number of equations " << mphysics->NEquations() << std::endl;
            int numthreads = 8;
            std::cout << "Number of threads " << numthreads << std::endl;

            TPZAnalysis an(mphysics);
            SolveSyst(an, mphysics,numthreads);
            
            
//            ofstream arg5("cmeshmultiphysics.txt");
//            mphysics->Print(arg5);

            //Calculo do erro
            TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
            TPZVec<REAL> erros;
    
            saidaerro << "Valor de epsilone " << EPSILON << std::endl;
            saidaerro << "Numero de threads " << numthreads << std::endl;
            saidaerro<<" \nErro da simulacao multifisica do fluxo (q)" <<endl;
            ErrorHDiv(cmesh1, saidaerro);
            
            saidaerro<<" Erro da simulacao multifisica da pressao (p)" <<endl;
            ErrorL2(cmesh2, saidaerro);
            
            std::cout << "Postprocessed\n";
            
            //Plot da solucao aproximada
//            string plotfile("Solution_mphysics.vtk");
//            PosProcessMultphysics(meshvec,  mphysics, an, plotfile);
        }
//.........这里部分代码省略.........

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

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

int main2()
{
	
#ifdef LOG4CXX
	{
		InitializePZLOG();
		std::stringstream sout;
		sout<< "Testando p adaptatividade"<<endl;
		LOGPZ_DEBUG(logger, sout.str().c_str());
	}
#endif
    
    //std::ofstream erro("TaxaArcTanQuadUni.txt");
   // std::ofstream erro("TaxaArcTanTriangNaoUni.txt");
    //std::ofstream erro("TaxaProbJuan.txt");
    //std::ofstream erro("TaxaCNMACQuadUni.txt");
    //    std::ofstream erro("TaxaCNMACQuadNaoUni.txt");
    std::ofstream erro("TaxaCNMACTriangUni.txt");
    
   //bool ftriang=false;
    bool prefine=false;
//    bool nouniform=false;//false para ref uniforme
//    REAL Lx=1;
//    REAL Ly=1;
    
    
    
    //TPZGeoMesh *gmesh = GMesh(ftriang, Lx,  Ly);
    //UniformRefine(gmesh, 1);
    //RefiningNearCircunference(2,gmesh,1,1);
    //std::ofstream filemesh2("MalhaGeoIni.vtk");
    //PrintGMeshVTK( gmesh, filemesh2);
   // TPZGeoMesh *gmesh =GMesh(ftriang, Lx,  Ly);//MalhaGeoT(h,hrefine);
    TPZGeoMesh *gmesh =MalhaGeo2(1);
    ofstream arg("gmeshInicial.txt");
    gmesh->Print(arg);
    
	TPZVec<REAL> calcErro;
	for (int porder=1; porder<2; porder++) {
        
//        TPZCompMesh *cmesh = CompMeshPAdap(*gmesh,porder,prefine);
//        ofstream arg2("CmeshInicial.txt");
//        cmesh->Print(arg2);
//        
//        SetDifferentOrderP(cmesh,porder);
//        ofstream arg22("CmeshDifOrder.txt");
//        cmesh->Print(arg22);

        
        erro<<"ordem "<<porder <<std::endl;
        for(int h=1;h<4;h++){
        //    erro << "NRef " << h << std::endl;
            
            erro<<std::endl;
            
            //1. Criacao da malha geom. e computacional
 //           bool hrefine=true;//true nao uniforme
           // TPZGeoMesh *gmesh =MalhaGeo2(h);//MalhaGeoT(h,nouniform);//MalhaGeo(h,nouniform);//GMesh(ftriang, Lx,  Ly);//MalhaGeo(h,nonuniform);
            

//            UniformRefine(gmesh, h);
//            ofstream arg3("gmeshRefinada.txt");
//            gmesh->Print(arg3);
//         NoUniformRefine(gmesh, h);
//            std::ofstream filemesh("MalhaGeoNaoUniT.vtk");
//            std::ofstream filemesh("MalhaGeoUniQCNMAC.vtk");
            //std::ofstream filemesh("MalhaGeoNaoUniQCNMAC.vtk");
           std::ofstream filemesh("MalhaGeo4ElCNMAC.vtk");
            PrintGMeshVTK( gmesh, filemesh);
            

//            RefiningNearCircunference(2,gmesh,h,1);
//            std::ofstream filemesh2("MalhaGeoQArcTanRefeineNearCirc.vtk");
//            PrintGMeshVTK( gmesh, filemesh2);
            
            TPZCompMesh *cmesh = CompMeshPAdap(*gmesh,porder,prefine);//CompMeshPAdapJuan(*gmesh,porder,prefine);
            ofstream arg2("CmeshInicial.txt");
            cmesh->Print(arg2);
            
            SetDifferentOrderPMesh4Elem(cmesh,porder);
            ofstream arg22("CmeshDifOrder.txt");
            cmesh->Print(arg22);
            
            
            int nDofs;
            nDofs=cmesh->NEquations();
            
            erro<< "\nNRefinamento "<<h<< "   NDofs "<<nDofs<<std::endl;
            
            std::cout << "Neq = " << nDofs << std::endl;
            
            //2. Resolve o problema
            
            TPZAnalysis analysis(cmesh);
            int numthreads = 1;
            SolveSyst(analysis, cmesh, numthreads);
            // SolveLU ( analysis );
            
            
            
//.........这里部分代码省略.........

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