C++ NonLinearImplicitSystem::GetLevelToAssemble方法代码示例

void AssemblePoisson_AD(MultiLevelProblem& ml_prob) {
  //  ml_prob is the global object from/to where get/set all the data
  //  level is the level of the PDE system to be assembled
  //  levelMax is the Maximum level of the MultiLevelProblem
  //  assembleMatrix is a flag that tells if only the residual or also the matrix should be assembled

  // call the adept stack object
  adept::Stack& s = FemusInit::_adeptStack;

  //  extract pointers to the several objects that we are going to use
  NonLinearImplicitSystem* mlPdeSys   = &ml_prob.get_system<NonLinearImplicitSystem> ("Poisson");   // pointer to the linear implicit system named "Poisson"
  const unsigned level = mlPdeSys->GetLevelToAssemble();

  Mesh*           msh         = ml_prob._ml_msh->GetLevel(level);    // pointer to the mesh (level) object
  elem*           el          = msh->el;  // pointer to the elem object in msh (level)

  MultiLevelSolution*   mlSol = ml_prob._ml_sol;  // pointer to the multilevel solution object
  Solution*   sol             = ml_prob._ml_sol->GetSolutionLevel(level);    // pointer to the solution (level) object


  LinearEquationSolver* pdeSys  = mlPdeSys->_LinSolver[level];  // pointer to the equation (level) object
  SparseMatrix*   KK          = pdeSys->_KK;  // pointer to the global stifness matrix object in pdeSys (level)
  NumericVector*  RES         = pdeSys->_RES; // pointer to the global residual vector object in pdeSys (level)

  const unsigned  dim = msh->GetDimension(); // get the domain dimension of the problem
  unsigned dim2 = (3 * (dim - 1) + !(dim - 1));        // dim2 is the number of second order partial derivatives (1,3,6 depending on the dimension)
  unsigned    iproc = msh->processor_id(); // get the process_id (for parallel computation)

  // reserve memory for the local standar vectors
  const unsigned maxSize = static_cast< unsigned >(ceil(pow(3, dim)));          // conservative: based on line3, quad9, hex27

  //solution variable
  unsigned solUIndex;
  solUIndex = mlSol->GetIndex("U");    // get the position of "U" in the ml_sol object = 0
  unsigned solUType = mlSol->GetSolutionType(solUIndex);    // get the finite element type for "T"



  unsigned solUPdeIndex;
  solUPdeIndex = mlPdeSys->GetSolPdeIndex("U");    // get the position of "U" in the pdeSys object = 0

  std::cout << solUIndex << " " << solUPdeIndex << std::endl;


  vector < adept::adouble >  solU; // local solution
  vector< adept::adouble > aResU; // local redidual vector

  vector < vector < double > > crdX(dim);    // local coordinates
  unsigned crdXType = 2; // get the finite element type for "x", it is always 2 (LAGRANGE QUADRATIC)

  solU.reserve(maxSize);
  aResU.reserve(maxSize);

  for (unsigned  k = 0; k < dim; k++) {
    crdX[k].reserve(maxSize);
  }

  vector <double> phi;  // local test function
  vector <double> phi_x; // local test function first order partial derivatives
  vector <double> phi_xx; // local test function second order partial derivatives

  phi.reserve(maxSize);
  phi_x.reserve(maxSize * dim);
  phi_xx.reserve(maxSize * dim2);

  double weight; // gauss point weight

  vector< int > sysDof; // local to global pdeSys dofs
  sysDof.reserve(maxSize);

  vector< double > ResU; // local redidual vector
  ResU.reserve(maxSize);

  vector < double > Jac;
  Jac.reserve(maxSize * maxSize);

  KK->zero(); // Set to zero all the entries of the Global Matrix

  // element loop: each process loops only on the elements that owns
  for (int iel = msh->_elementOffset[iproc]; iel < msh->_elementOffset[iproc + 1]; iel++) {

    unsigned kel = iel; //msh->IS_Mts2Gmt_elem[iel]; // mapping between paralell dof and mesh dof
    short unsigned kelGeom = el->GetElementType(kel);      // element geometry type

    unsigned nDofsU = el->GetElementDofNumber(kel, solUType);      // number of solution element dofs
    unsigned nDofsX = el->GetElementDofNumber(kel, crdXType);      // number of solution element dofs

    // resize local arrays
    sysDof.resize(nDofsU);

    solU.resize(nDofsU);

    for (unsigned  k = 0; k < dim; k++) {
      crdX[k].resize(nDofsX);
    }

    aResU.assign(nDofsU, 0);

    // local storage of global mapping and solution
    for (unsigned i = 0; i < nDofsU; i++) {
//.........这里部分代码省略.........

void AssembleBoussinesqAppoximation(MultiLevelProblem& ml_prob)
{
  //  ml_prob is the global object from/to where get/set all the data
  //  level is the level of the PDE system to be assembled
  //  levelMax is the Maximum level of the MultiLevelProblem
  //  assembleMatrix is a flag that tells if only the residual or also the matrix should be assembled

  //  extract pointers to the several objects that we are going to use

  double Mu;
  
  if(counter < c0 ) Mu = 2. * Miu;
  else if ( counter <= cn ) {
    Mu = 2*Miu*(cn-counter)/(cn-c0) + Miu*(counter-c0)/(cn-c0);
  }
  else{
    Mu = Miu;
  }
  std::cout << counter << " " << Mu <<std::endl;
  counter++;

  NonLinearImplicitSystem* mlPdeSys   = &ml_prob.get_system<NonLinearImplicitSystem> ("NS");   // pointer to the linear implicit system named "Poisson"
  const unsigned level = mlPdeSys->GetLevelToAssemble();

  Mesh*           msh         = ml_prob._ml_msh->GetLevel(level);    // pointer to the mesh (level) object
  elem*           el          = msh->el;  // pointer to the elem object in msh (level)

  MultiLevelSolution*   mlSol         = ml_prob._ml_sol;  // pointer to the multilevel solution object
  Solution*   sol         = ml_prob._ml_sol->GetSolutionLevel(level);    // pointer to the solution (level) object


  LinearEquationSolver* pdeSys        = mlPdeSys->_LinSolver[level];  // pointer to the equation (level) object

  bool assembleMatrix = mlPdeSys->GetAssembleMatrix();
  // call the adept stack object

  SparseMatrix*   KK          = pdeSys->_KK;  // pointer to the global stifness matrix object in pdeSys (level)
  NumericVector*  RES         = pdeSys->_RES; // pointer to the global residual vector object in pdeSys (level)

  const unsigned  dim = msh->GetDimension(); // get the domain dimension of the problem
  unsigned dim2 = (3 * (dim - 1) + !(dim - 1));        // dim2 is the number of second order partial derivatives (1,3,6 depending on the dimension)
  unsigned    iproc = msh->processor_id(); // get the process_id (for parallel computation)

  // reserve memory for the local standar vectors
  const unsigned maxSize = static_cast< unsigned >(ceil(pow(3, dim)));          // conservative: based on line3, quad9, hex27

  //solution variable
  vector < unsigned > solVIndex(dim);
  solVIndex[0] = mlSol->GetIndex("U");    // get the position of "U" in the ml_sol object
  solVIndex[1] = mlSol->GetIndex("V");    // get the position of "V" in the ml_sol object
  if (dim == 3) solVIndex[2] = mlSol->GetIndex("W");      // get the position of "V" in the ml_sol object

  unsigned solVType = mlSol->GetSolutionType(solVIndex[0]);    // get the finite element type for "u"

  unsigned solPIndex;
  solPIndex = mlSol->GetIndex("P");    // get the position of "P" in the ml_sol object
  unsigned solPType = mlSol->GetSolutionType(solPIndex);    // get the finite element type for "u"

  vector < unsigned > solVPdeIndex(dim);
  solVPdeIndex[0] = mlPdeSys->GetSolPdeIndex("U");    // get the position of "U" in the pdeSys object
  solVPdeIndex[1] = mlPdeSys->GetSolPdeIndex("V");    // get the position of "V" in the pdeSys object
  if (dim == 3) solVPdeIndex[2] = mlPdeSys->GetSolPdeIndex("W");

  unsigned solPPdeIndex;
  solPPdeIndex = mlPdeSys->GetSolPdeIndex("P");    // get the position of "P" in the pdeSys object

  vector < vector < double > >  solV(dim);    // local solution
  vector < double >  solP; // local solution

  vector < vector < double > > coordX(dim);    // local coordinates
  unsigned coordXType = 2; // get the finite element type for "x", it is always 2 (LAGRANGE QUADRATIC)

  for (unsigned  k = 0; k < dim; k++) {
    solV[k].reserve(maxSize);
    coordX[k].reserve(maxSize);
  }

  solP.reserve(maxSize);

  vector <double> phiV;  // local test function
  vector <double> phiV_x; // local test function first order partial derivatives
  vector <double> phiV_xx; // local test function second order partial derivatives

  phiV.reserve(maxSize);
  phiV_x.reserve(maxSize * dim);
  phiV_xx.reserve(maxSize * dim2);

  double* phiP;
  double weight; // gauss point weight

  vector< int > sysDof; // local to global pdeSys dofs
  sysDof.reserve((dim + 2) *maxSize);

  vector< double > Res; // local redidual vector
  Res.reserve((dim + 2) *maxSize);

  vector < double > Jac;
  Jac.reserve((dim + 2) *maxSize * (dim + 2) *maxSize);

  if (assembleMatrix)
//.........这里部分代码省略.........

std::pair <vector<double>, vector <double> > GetVaribleValues(MultiLevelProblem& ml_prob, const unsigned &elem, const std::vector<double>&xi) {

  NonLinearImplicitSystem* mlPdeSys   = &ml_prob.get_system<NonLinearImplicitSystem> ("NS");   // pointer to the linear implicit system named "Poisson"
  const unsigned level = mlPdeSys->GetLevelToAssemble();

  Mesh* msh = ml_prob._ml_msh->GetLevel(level);

  MultiLevelSolution* mlSol = ml_prob._ml_sol;  // pointer to the multilevel solution object
  Solution* sol = ml_prob._ml_sol->GetSolutionLevel(level);    // pointer to the solution (level) object
  LinearEquationSolver* pdeSys = mlPdeSys->_LinSolver[level];  // pointer to the equation (level) object
  const unsigned  dim = msh->GetDimension(); // get the domain dimension of the problem
  
  //solution variable
  unsigned solTIndex;
  solTIndex = mlSol->GetIndex("T");    // get the position of "T" in the ml_sol object
  unsigned solTType = mlSol->GetSolutionType(solTIndex);    // get the finite element type for "T"
  vector < double >  solT; // local solution
  
  vector <unsigned> solVIndex(dim);
  solVIndex[0] = mlSol->GetIndex("U");
  solVIndex[1] = mlSol->GetIndex("V");
  if (dim==3) solVIndex[2] = mlSol->GetIndex("W");
  unsigned solVType = mlSol->GetSolutionType(solVIndex[0]);
  vector < vector <double> > solV(dim);
  
  unsigned iproc = msh->processor_id(); // get the process_id (for parallel computation)
  MyVector <double> solTXi(1, 0.);
  solTXi.stack();
  MyVector <double> solUXi(1, 0.);
  solUXi.stack();
  MyVector <double> solVXi(1, 0.);
  solVXi.stack();
  /*
  if (dim==3) {
  	MyVector <double> solWXi(1, 0.);
    solWXi.stack();
  }
  */
  //BEGIN local dof number extraction
  if(elem >= msh->_elementOffset[iproc] && elem  < msh->_elementOffset[iproc + 1]) {
	unsigned nDofsT = msh->GetElementDofNumber(elem, solTType);  //temperature
  	solT.reserve(nDofsT);
  
    unsigned nDofsV = msh->GetElementDofNumber(elem, solVType); //velocity
    for(unsigned  k = 0; k < dim; k++) solV[k].reserve(nDofsV);

  	//BEGIN global to local extraction
  	for(unsigned i = 0; i < nDofsT; i++) { //temperature
      unsigned solTDof = msh->GetSolutionDof(i, elem, solTType);  //local to global solution dof
      solT[i] = (*sol->_Sol[solTIndex])(solTDof);  //global to local solution value
  	}

    for(unsigned i = 0; i < nDofsV; i++){ //velocity
      unsigned solVDof = msh->GetSolutionDof(i, elem, solVType);
      for(unsigned  k = 0; k < dim; k++) {
      solV[k][i] = (*sol->_Sol[solVIndex[k]])(solVDof); // global extraction and local storage for the solution
     }
   }
  
   short unsigned ielGeom = msh->GetElementType(elem);
   for(unsigned i = 0; i < nDofsT; i++) {
      basis *base = msh->_finiteElement[ielGeom][solTType]->GetBasis();
      double phiT = base->eval_phi(base->GetIND(i), &xi[0]);
      solTXi[iproc] += phiT * solT[i];
   }
   
   for(unsigned i = 0; i < nDofsV; i++) {
    basis *base = msh->_finiteElement[ielGeom][solVType]->GetBasis();
    double phiV = base->eval_phi(base->GetIND(i), &xi[0]);
    solUXi[iproc] += phiV * solV[0][i];
    solVXi[iproc] += phiV * solV[1][i];
 //   if(dim==3) solWXi[iproc] += phiV * solV[2][i];
   }
  }
  
  std::pair <vector <double>, vector <double> > out_value;
  unsigned mproc = msh->IsdomBisectionSearch(elem , 3);
  solUXi.broadcast(mproc);
  solVXi.broadcast(mproc);
  solTXi.broadcast(mproc);
  
  vector <double> solV_pt(dim);
  solV_pt[0]= solUXi[mproc];
  solV_pt[1]= solVXi[mproc];
  out_value.first = solV_pt;
  
  vector <double> solPT_pt(2);
  solPT_pt[0]= 0.;
  solPT_pt[1]=solTXi[mproc];
  out_value.second = solPT_pt;
  
  solUXi.clearBroadcast();
  solVXi.clearBroadcast();
  solTXi.clearBroadcast();
  return out_value;
}