C++ NumericVector::close方法代码示例

//.........这里部分代码省略.........
      for (unsigned jdim = 0; jdim < dim; jdim++) {
        x[jdim][i] = (*msh->_coordinate->_Sol[jdim])(xDof);      // global extraction and local storage for the element coordinates
      }
    }

    if (level == levelMax || !el->GetRefinedElementIndex(kel)) {      // do not care about this if now (it is used for the AMR)
      // start a new recording of all the operations involving adept::adouble variables
      s.new_recording();

      // *** Gauss point loop ***
      for (unsigned ig = 0; ig < msh->_finiteElement[kelGeom][soluType]->GetGaussPointNumber(); ig++) {
        // *** get gauss point weight, test function and test function partial derivatives ***
        msh->_finiteElement[kelGeom][soluType]->Jacobian(x, ig, weight, phi, phi_x, phi_xx);

        // evaluate the solution, the solution derivatives and the coordinates in the gauss point
        adept::adouble soluGauss = 0;
        vector < adept::adouble > soluGauss_x(dim, 0.);

        adept::adouble solvGauss = 0;
        vector < adept::adouble > solvGauss_x(dim, 0.);

        vector < double > xGauss(dim, 0.);

        for (unsigned i = 0; i < nDofs; i++) {
          soluGauss += phi[i] * solu[i];
          solvGauss += phi[i] * solv[i];

          for (unsigned jdim = 0; jdim < dim; jdim++) {
            soluGauss_x[jdim] += phi_x[i * dim + jdim] * solu[i];
            solvGauss_x[jdim] += phi_x[i * dim + jdim] * solv[i];
            xGauss[jdim] += x[jdim][i] * phi[i];
          }
        }

        // *** phi_i loop ***
        for (unsigned i = 0; i < nDofs; i++) {

          adept::adouble Laplace_u = 0.;
          adept::adouble Laplace_v = 0.;

          for (unsigned jdim = 0; jdim < dim; jdim++) {
            Laplace_u   +=  - phi_x[i * dim + jdim] * soluGauss_x[jdim];
            Laplace_v   +=  - phi_x[i * dim + jdim] * solvGauss_x[jdim];
          }

          double exactSolValue = GetExactSolutionValue(xGauss);

          double pi = acos(-1.);

          aResv[i] += (solvGauss * phi[i] -  Laplace_u) * weight;
          aResu[i] += (4.*pi * pi * pi * pi * exactSolValue * phi[i] -  Laplace_v) * weight;

        } // end phi_i loop
      } // end gauss point loop
    } // endif single element not refined or fine grid loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store
    for (int i = 0; i < nDofs; i++) {
      Res[i]       = aResu[i].value();
      Res[nDofs + i] = aResv[i].value();
    }

    RES->add_vector_blocked(Res, KKDof);

    if (assembleMatrix) {

      // define the dependent variables
      s.dependent(&aResu[0], nDofs);
      s.dependent(&aResv[0], nDofs);

      // define the independent variables
      s.independent(&solu[0], nDofs);
      s.independent(&solv[0], nDofs);
      // get the jacobian matrix (ordered by column)
      s.jacobian(&Jac[0]);

      // get the jacobian matrix (ordered by raw, i.e. Jact=Jac^t)
      for (int inode = 0; inode < 2 * nDofs; inode++) {
        for (int jnode = 0; jnode < 2 * nDofs; jnode++) {
          Jact[inode * 2 * nDofs + jnode] = -Jac[jnode * 2 * nDofs + inode];
        }
      }

      //store Jact in the global matrix KK
      KK->add_matrix_blocked(Jact, KKDof, KKDof);

      s.clear_independents();
      s.clear_dependents();
    }
  } //end element loop for each process

  RES->close();

  if (assembleMatrix) KK->close();

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
            
//             for (unsigned j = 0; j < nDofsP; j++) {
//               M[i] += phiP[i] * phiP[j] * weight;
//             }
            
          }

        }
      }

      //END phiP_i loop

    }
    //END Gauss point loop

    
    std::vector<std::vector<double> > BtMinv (dim * nDofsV);
    for(unsigned i=0; i< dim * nDofsV; i++){
      BtMinv[i].resize(nDofsP);
    }
    
    std::vector<std::vector<double> > B(nDofsP);
    for(unsigned i=0; i< nDofsP; i++){
      B[i].resize(dim * nDofsV);
    }
    
    for(unsigned i = 0; i < dim * nDofsV; i++){
      unsigned irow =  i * nDofsVP;
      for(unsigned j = 0; j < nDofsP; j++){
	unsigned jcol = (dim * nDofsV) + j;
	BtMinv[i][j] = .1 * Jac[ irow +  jcol] / M[j];
      }
    }
    
    for(unsigned i = 0; i < nDofsP; i++){
      unsigned irow = ( (dim * nDofsV) + i) * nDofsVP;
      for(unsigned j = 0; j < dim * nDofsV; j++){
	B[i][j] = Jac[irow + j];
      }
    }
    
    std::vector<std::vector<double> > Jg (dim * nDofsV);
    for(unsigned i=0; i< dim * nDofsV; i++){
      Jg[i].resize(dim * nDofsV);
    }
        
    for(unsigned i = 0; i < dim * nDofsV; i++){
      for(unsigned j = 0; j < dim * nDofsV; j++){
	Jg[i][j] = 0.;
	for(unsigned k = 0; k < nDofsP; k++){
	  Jg[i][j] += BtMinv[i][k] * B[k][j];
	}
      }
    }
   
    std::vector<double> fg (dim * nDofsV);
    for(unsigned i = 0; i < dim * nDofsV; i++){
      fg[i] = 0;
      for(unsigned j = 0; j < nDofsP; j++){
	fg[i] += BtMinv[i][j] * Res[dim * nDofsV + j];
      }
    }
    
    
     for(unsigned i = 0; i < dim * nDofsV; i++){
	unsigned irow =  i * nDofsVP;
	for(unsigned j = 0; j < dim * nDofsV; j++){
	  Jac[irow + j] += Jg[i][j];
	  //std::cout<< Jg[i][j]<<" ";
	  //std::cout<< Jac[irow + j]<<" ";
	}
	//std::cout<<"\n";
	Res[i] += fg[i];
     }
    //std::cout<<"\n";
    
    
    
    
    
    //BEGIN local to global Matrix/Vector assembly
    RES->add_vector_blocked(Res, sysDof);

    if (assembleMatrix) {
      KK->add_matrix_blocked(Jac, sysDof, sysDof);
    }

    //END local to global Matrix/Vector assembly

  }

  //END element loop

  RES->close();

  if (assembleMatrix) {
    KK->close();
  }
  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
    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++) {
      unsigned solUDof = msh->GetSolutionDof(i, iel, solUType);    // local to global mapping of the solution U
      solU[i] = (*sol->_Sol[solUIndex])(solUDof);      // value of the solution U in the dofs
      sysDof[i] = pdeSys->GetSystemDof(solUIndex, solUPdeIndex, i, iel);    // local to global mapping between solution U and system
    }

    // local storage of coordinates
    for (unsigned i = 0; i < nDofsX; i++) {
      unsigned coordXDof  = msh->GetSolutionDof(i, iel, crdXType);   // local to global mapping of the coordinate X[dim]

      for (unsigned k = 0; k < dim; k++) {
        crdX[k][i] = (*msh->_topology->_Sol[k])(coordXDof);      // value of the solution X[dim]
      }
    }

    s.new_recording();

    // *** Gauss point loop ***
    for (unsigned ig = 0; ig < msh->_finiteElement[kelGeom][solUType]->GetGaussPointNumber(); ig++) {
      // *** get gauss point weight, test function and test function partial derivatives ***
      msh->_finiteElement[kelGeom][solUType]->Jacobian(crdX, ig, weight, phi, phi_x, phi_xx);

      adept::adouble solUig = 0; // solution U in the gauss point
      vector < adept::adouble > gradSolUig(dim, 0.); // graadient of solution U in the gauss point

      for (unsigned i = 0; i < nDofsU; i++) {
        solUig += phi[i] * solU[i];

        for (unsigned j = 0; j < dim; j++) {
          gradSolUig[j] += phi_x[i * dim + j] * solU[i];
        }
      }

      double nu = 1.;

      // *** phiU_i loop ***
      for (unsigned i = 0; i < nDofsU; i++) {
        adept::adouble LaplaceU = 0.;

        for (unsigned j = 0; j < dim; j++) {
          LaplaceU +=  nu * phi_x[i * dim + j] * gradSolUig[j];
        }

        aResU[i] += (phi[i] - LaplaceU) * weight;
      } // end phiU_i loop
    } // end gauss point loop

    // } // endif single element not refined or fine grid loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store them in RES
    ResU.resize(nDofsU);    //resize

    for (int i = 0; i < nDofsU; i++) {
      ResU[i] = -aResU[i].value();
    }

    RES->add_vector_blocked(ResU, sysDof);

    //Extarct and store the Jacobian

    Jac.resize(nDofsU * nDofsU);
    // define the dependent variables
    s.dependent(&aResU[0], nDofsU);

    // define the independent variables
    s.independent(&solU[0], nDofsU);

    // get the and store jacobian matrix (row-major)
    s.jacobian(&Jac[0] , true);
    KK->add_matrix_blocked(Jac, sysDof, sysDof);

    s.clear_independents();
    s.clear_dependents();

  } //end element loop for each process

  RES->close();

  KK->close();

  // ***************** END ASSEMBLY *******************
}

Real NewtonSolver::line_search(Real tol,
                               Real last_residual,
                               Real &current_residual,
                               NumericVector<Number> &newton_iterate,
                               const NumericVector<Number> &linear_solution)
{
  // Take a full step if we got a residual reduction or if we
  // aren't substepping
  if ((current_residual < last_residual) ||
      (!require_residual_reduction &&
       (!require_finite_residual || !libmesh_isnan(current_residual))))
    return 1.;

  // The residual vector
  NumericVector<Number> &rhs = *(_system.rhs);

  Real ax = 0.;  // First abscissa, don't take negative steps
  Real cx = 1.;  // Second abscissa, don't extrapolate steps

  // Find bx, a step length that gives lower residual than ax or cx
  Real bx = 1.;

  while (libmesh_isnan(current_residual) ||
         (current_residual > last_residual &&
          require_residual_reduction))
    {
      // Reduce step size to 1/2, 1/4, etc.
      Real substepdivision;
      if (brent_line_search && !libmesh_isnan(current_residual))
        {
          substepdivision = std::min(0.5, last_residual/current_residual);
          substepdivision = std::max(substepdivision, tol*2.);
        }
      else
        substepdivision = 0.5;

      newton_iterate.add (bx * (1.-substepdivision),
                          linear_solution);
      newton_iterate.close();
      bx *= substepdivision;
      if (verbose)
        libMesh::out << "  Shrinking Newton step to "
                  << bx << std::endl;

      // Check residual with fractional Newton step
      _system.assembly (true, false);

      rhs.close();
      current_residual = rhs.l2_norm();
      if (verbose)
        libMesh::out << "  Current Residual: "
                  << current_residual << std::endl;

      if (bx/2. < minsteplength &&
          (libmesh_isnan(current_residual) ||
           (current_residual > last_residual)))
        {
          libMesh::out << "Inexact Newton step FAILED at step "
                    << _outer_iterations << std::endl;

          if (!continue_after_backtrack_failure)
            {
              libmesh_convergence_failure();
            }
          else
            {
              libMesh::out << "Continuing anyway ..." << std::endl;
              _solve_result = DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
              return bx;
            }
        }
    } // end while (current_residual > last_residual)

  // Now return that reduced-residual step, or  use Brent's method to
  // find a more optimal step.

  if (!brent_line_search)
    return bx;

  // Brent's method adapted from Numerical Recipes in C, ch. 10.2
  Real e = 0.;

  Real x = bx, w = bx, v = bx;

  // Residuals at bx
  Real fx = current_residual,
       fw = current_residual,
       fv = current_residual;

  // Max iterations for Brent's method loop
  const unsigned int max_i = 20;

  // for golden ratio steps
  const Real golden_ratio = 1.-(std::sqrt(5.)-1.)/2.;

  for (unsigned int i=1; i <= max_i; i++)
    {
      Real xm = (ax+cx)*0.5;
      Real tol1 = tol * std::abs(x) + tol*tol;
      Real tol2 = 2.0 * tol1;
//.........这里部分代码省略.........

//.........这里部分代码省略.........


      // *** phiV_i loop ***
      for(unsigned i = 0; i < nDofsV; i++) {
        vector < adept::adouble > NSV(dim, 0.);
        vector < double > NSVold(dim, 0.);

        for(unsigned j = 0; j < dim; j++) {
          for(unsigned  k = 0; k < dim; k++) {
            NSV[k]   +=  sqrt(Pr / Ra) * phiV_x[i * dim + j] * (gradSolV_gss[k][j] + gradSolV_gss[j][k]);
            NSV[k]   +=  phiV[i] * (solV_gss[j] * gradSolV_gss[k][j]);

            NSVold[k]   +=  sqrt(Pr / Ra) * phiV_x[i * dim + j] * (gradSolVold_gss[k][j] + gradSolVold_gss[j][k]);
            NSVold[k]   +=  phiV[i] * (solVold_gss[j] * gradSolVold_gss[k][j]);

          }
        }

        for(unsigned  k = 0; k < dim; k++) {
          NSV[k] += -solP_gss * phiV_x[i * dim + k];
          NSVold[k] += -solPold_gss * phiV_x[i * dim + k];
        }

        NSV[1] += -beta * solT_gss * phiV[i];
        NSVold[1] += -beta * solTold_gss * phiV[i];

        for(unsigned  k = 0; k < dim; k++) {
          aResV[k][i] += (- (solV_gss[k] - solVold_gss[k]) * phiV[i] / dt - 0.5 * (NSV[k] + NSVold[k])) * weight;
        }
      } // end phiV_i loop

      // *** phiP_i loop ***
      for(unsigned i = 0; i < nDofsP; i++) {
        for(int k = 0; k < dim; k++) {
          aResP[i] += - (gradSolV_gss[k][k]) * phiP[i]  * weight;
        }
      } // end phiP_i loop

    } // end gauss point loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store them in RES
    Res.resize(nDofsTVP);    //resize

    for(int i = 0; i < nDofsT; i++) {
      Res[i] = -aResT[i].value();
    }

    for(int i = 0; i < nDofsV; i++) {
      for(unsigned  k = 0; k < dim; k++) {
        Res[ i + nDofsT + k * nDofsV ] = -aResV[k][i].value();
      }
    }

    for(int i = 0; i < nDofsP; i++) {
      Res[ i + nDofsT + dim * nDofsV ] = -aResP[i].value();
    }

    RES->add_vector_blocked(Res, sysDof);

    //Extarct and store the Jacobian
    if(assembleMatrix) {
      Jac.resize(nDofsTVP * nDofsTVP);
      // define the dependent variables
      s.dependent(&aResT[0], nDofsT);

      for(unsigned  k = 0; k < dim; k++) {
        s.dependent(&aResV[k][0], nDofsV);
      }

      s.dependent(&aResP[0], nDofsP);

      // define the independent variables
      s.independent(&solT[0], nDofsT);

      for(unsigned  k = 0; k < dim; k++) {
        s.independent(&solV[k][0], nDofsV);
      }

      s.independent(&solP[0], nDofsP);

      // get the and store jacobian matrix (row-major)
      s.jacobian(&Jac[0] , true);
      KK->add_matrix_blocked(Jac, sysDof, sysDof);

      s.clear_independents();
      s.clear_dependents();
    }
  } //end element loop for each process

  RES->close();

  if(assembleMatrix) {
    KK->close();
  }

  // ***************** END ASSEMBLY *******************
}