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

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

//.........这里部分代码省略.........
   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);
   U.Print("Resultado",outcm2);
   secondmesh->LoadSolution(U);
   secondmesh->Solution().Print("Mesh solution ",outcm2);
//   TPZElementMatrix ek,ef;
//   secondmesh->ElementVec()[0]->CalcStiff(ek,ef);
//	ek.fMat->Print();
//   ef.fMat->Print();
   delete secondmesh;
   delete firstmesh;
   return 0;
}