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

void AssembleOptimization::assemble_A_and_F()
{
  A_matrix->zero();
  F_vector->zero();

  const MeshBase & mesh = _sys.get_mesh();

  const unsigned int dim = mesh.mesh_dimension();
  const unsigned int u_var = _sys.variable_number ("u");

  const DofMap & dof_map = _sys.get_dof_map();
  FEType fe_type = dof_map.variable_type(u_var);
  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
  QGauss qrule (dim, fe_type.default_quadrature_order());
  fe->attach_quadrature_rule (&qrule);

  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<Real> > & phi = fe->get_phi();
  const std::vector<std::vector<RealGradient> > & dphi = fe->get_dphi();

  std::vector<dof_id_type> dof_indices;

  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      const Elem * elem = *el;

      dof_map.dof_indices (elem, dof_indices);

      const unsigned int n_dofs = dof_indices.size();

      fe->reinit (elem);

      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);

      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          for (unsigned int dof_i=0; dof_i<n_dofs; dof_i++)
            {
              for (unsigned int dof_j=0; dof_j<n_dofs; dof_j++)
                {
                  Ke(dof_i, dof_j) += JxW[qp] * (dphi[dof_j][qp]* dphi[dof_i][qp]);
                }
              Fe(dof_i) += JxW[qp] * phi[dof_i][qp];
            }
        }

      A_matrix->add_matrix (Ke, dof_indices);
      F_vector->add_vector (Fe, dof_indices);
    }

  A_matrix->close();
  F_vector->close();
}

void AssembleOptimization::gradient (const NumericVector<Number> & soln,
                                     NumericVector<Number> & grad_f,
                                     OptimizationSystem & /*sys*/)
{
  grad_f.zero();

  A_matrix->vector_mult(grad_f, soln);
  grad_f.add(-1, *F_vector);
}

void SumShellMatrix<T>::get_diagonal (NumericVector<T>& dest) const
{
  AutoPtr<NumericVector<T> > a = dest.clone();
  dest.zero();
  for(numeric_index_type i=matrices.size(); i-->0; )
    {
      matrices[i]->get_diagonal(*a);
      dest += *a;
    }
}

void AssembleOptimization::inequality_constraints (const NumericVector<Number> & X,
                                                   NumericVector<Number> & C_ineq,
                                                   OptimizationSystem & /*sys*/)
{
  C_ineq.zero();

  std::unique_ptr<NumericVector<Number>> X_localized =
    NumericVector<Number>::build(X.comm());
  X_localized->init(X.size(), false, SERIAL);
  X.localize(*X_localized);

  std::vector<Number> constraint_values(1);
  constraint_values[0] = (*X_localized)(200)*(*X_localized)(200) + (*X_localized)(201) - 5.;

  for (std::size_t i=0; i<constraint_values.size(); i++)
    if ((C_ineq.first_local_index() <= i) && (i < C_ineq.last_local_index()))
      C_ineq.set(i, constraint_values[i]);
}

void AssembleOptimization::equality_constraints (const NumericVector<Number> & X,
                                                 NumericVector<Number> & C_eq,
                                                 OptimizationSystem & /*sys*/)
{
  C_eq.zero();

  UniquePtr<NumericVector<Number> > X_localized =
    NumericVector<Number>::build(X.comm());
  X_localized->init(X.size(), false, SERIAL);
  X.localize(*X_localized);

  std::vector<Number> constraint_values(3);
  constraint_values[0] = (*X_localized)(17);
  constraint_values[1] = (*X_localized)(23);
  constraint_values[2] = (*X_localized)(98) + (*X_localized)(185);

  for (unsigned int i=0; i<constraint_values.size(); i++)
    if ((C_eq.first_local_index() <= i) &&
        (i < C_eq.last_local_index()))
      C_eq.set(i, constraint_values[i]);
}

void GetSolutionNorm(MultiLevelSolution& mlSol, const unsigned & group, std::vector <double> &data)
{

  int  iproc, nprocs;
  MPI_Comm_rank(MPI_COMM_WORLD, &iproc);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

  NumericVector* p2;
  NumericVector* v2;
  NumericVector* vol;
  NumericVector* vol0;
  p2 = NumericVector::build().release();
  v2 = NumericVector::build().release();
  vol = NumericVector::build().release();
  vol0 = NumericVector::build().release();

  if(nprocs == 1) {
    p2->init(nprocs, 1, false, SERIAL);
    v2->init(nprocs, 1, false, SERIAL);
    vol->init(nprocs, 1, false, SERIAL);
    vol0->init(nprocs, 1, false, SERIAL);
  }
  else {
    p2->init(nprocs, 1, false, PARALLEL);
    v2->init(nprocs, 1, false, PARALLEL);
    vol->init(nprocs, 1, false, PARALLEL);
    vol0->init(nprocs, 1, false, PARALLEL);
  }

  p2->zero();
  v2->zero();
  vol->zero();
  vol0->zero();

  unsigned level = mlSol._mlMesh->GetNumberOfLevels() - 1;

  Solution* solution  = mlSol.GetSolutionLevel(level);
  Mesh* msh = mlSol._mlMesh->GetLevel(level);


  const unsigned dim = msh->GetDimension();


  const unsigned max_size = static_cast< unsigned >(ceil(pow(3, dim)));

  vector< double > solP;
  vector< vector < double> >  solV(dim);
  vector< vector < double> > x0(dim);
  vector< vector < double> > x(dim);

  solP.reserve(max_size);
  for(unsigned d = 0; d < dim; d++) {
    solV[d].reserve(max_size);
    x0[d].reserve(max_size);
    x[d].reserve(max_size);
  }
  double weight;
  double weight0;

  vector <double> phiV;
  vector <double> gradphiV;
  vector <double> nablaphiV;

  double *phiP;

  phiV.reserve(max_size);
  gradphiV.reserve(max_size * dim);
  nablaphiV.reserve(max_size * (3 * (dim - 1) + !(dim - 1)));

  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"

  vector < unsigned > solDIndex(dim);
  solDIndex[0] = mlSol.GetIndex("DX");    // get the position of "U" in the ml_sol object
  solDIndex[1] = mlSol.GetIndex("DY");    // get the position of "V" in the ml_sol object
  if(dim == 3) solDIndex[2] = mlSol.GetIndex("DZ");       // get the position of "V" in the ml_sol object

  unsigned solDType = mlSol.GetSolutionType(solDIndex[0]);

  unsigned solPIndex;
  solPIndex = mlSol.GetIndex("PS");
  unsigned solPType = mlSol.GetSolutionType(solPIndex);

  for(int iel = msh->_elementOffset[iproc]; iel < msh->_elementOffset[iproc + 1]; iel++) {
    if(msh->GetElementGroup(iel) == group) {
      short unsigned ielt = msh->GetElementType(iel);
      unsigned ndofV = msh->GetElementDofNumber(iel, solVType);
      unsigned ndofP = msh->GetElementDofNumber(iel, solPType);
      unsigned ndofD = msh->GetElementDofNumber(iel, solDType);
      // resize

      phiV.resize(ndofV);
      gradphiV.resize(ndofV * dim);
      nablaphiV.resize(ndofV * (3 * (dim - 1) + !(dim - 1)));

      solP.resize(ndofP);
//.........这里部分代码省略.........

void ETD(MultiLevelProblem& ml_prob)
{

  const unsigned& NLayers = NumberOfLayers;

  adept::Stack& s = FemusInit::_adeptStack;

  LinearImplicitSystem* mlPdeSys  = &ml_prob.get_system<LinearImplicitSystem> ("SW");   // pointer to the linear implicit system named "Poisson"

  unsigned level = ml_prob._ml_msh->GetNumberOfLevels() - 1u;

  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)
  NumericVector* EPS = pdeSys->_EPS; // pointer to the global residual vector object in pdeSys (level)

  const unsigned  dim = msh->GetDimension(); // get the domain dimension of the problem

  unsigned    iproc = msh->processor_id(); // get the process_id (for parallel computation)
  unsigned    nprocs = msh->n_processors(); // get the process_id (for parallel computation)


  //solution variable
  std::vector < unsigned > solIndexh(NLayers);
  std::vector < unsigned > solPdeIndexh(NLayers);

  std::vector < unsigned > solIndexv(NLayers);
  std::vector < unsigned > solPdeIndexv(NLayers);
  
  std::vector < unsigned > solIndexHT(NLayers);
  std::vector < unsigned > solPdeIndexHT(NLayers);
  
  
   std::vector < unsigned > solIndexT(NLayers);

  vector< int > l2GMap; // local to global mapping
  
  for(unsigned i = 0; i < NLayers; i++) {
    char name[10];
    sprintf(name, "h%d", i);
    solIndexh[i] = mlSol->GetIndex(name); // get the position of "hi" in the sol object
    solPdeIndexh[i] = mlPdeSys->GetSolPdeIndex(name); // get the position of "hi" in the pdeSys object

    sprintf(name, "v%d", i);
    solIndexv[i] = mlSol->GetIndex(name); // get the position of "vi" in the sol object
    solPdeIndexv[i] = mlPdeSys->GetSolPdeIndex(name); // get the position of "vi" in the pdeSys object
    
    sprintf(name, "HT%d", i);
    solIndexHT[i] = mlSol->GetIndex(name); // get the position of "Ti" in the sol object
    solPdeIndexHT[i] = mlPdeSys->GetSolPdeIndex(name); // get the position of "Ti" in the pdeSys object
    
    sprintf(name, "T%d", i);
    solIndexT[i] = mlSol->GetIndex(name); // get the position of "Ti" in the sol object
    
  }

  unsigned solTypeh = mlSol->GetSolutionType(solIndexh[0]);    // get the finite element type for "hi"
  unsigned solTypev = mlSol->GetSolutionType(solIndexv[0]);    // get the finite element type for "vi"
  unsigned solTypeHT = mlSol->GetSolutionType(solIndexHT[0]);    // get the finite element type for "Ti"
  
  vector < double > x;    // local coordinates
  vector< vector < adept::adouble > > solh(NLayers);    // local coordinates
  vector< vector < adept::adouble > > solv(NLayers);    // local coordinates
  vector< vector < adept::adouble > > solHT(NLayers);    // local coordinates
  vector< vector < double > > solT(NLayers);    // local coordinates
  
  vector< vector < bool > > bdch(NLayers);    // local coordinates
  vector< vector < bool > > bdcv(NLayers);    // local coordinates
  vector< vector < bool > > bdcHT(NLayers);    // local coordinates
  
  unsigned xType = 2; // get the finite element type for "x", it is always 2 (LAGRANGE QUADRATIC)

  vector < vector< adept::adouble > > aResh(NLayers);
  vector < vector< adept::adouble > > aResv(NLayers);
  vector < vector< adept::adouble > > aResHT(NLayers);
  
  KK->zero();
  RES->zero();

  double maxWaveSpeed = 0.;
  
  double dx;
  
  for(unsigned k=0; k<NumberOfLayers; k++){
    for(unsigned i =  msh->_dofOffset[solTypeHT][iproc]; i <  msh->_dofOffset[solTypeHT][iproc + 1]; i++){
      double valueT = (*sol->_Sol[solIndexT[k]])(i);
      double valueH = (*sol->_Sol[solIndexh[k]])(i);
            
      double valueHT = valueT * valueH;
    
      sol->_Sol[solIndexHT[k]]->set(i, valueHT);
    }
    sol->_Sol[solIndexHT[k]]->close();
//.........这里部分代码省略.........

void SparseMatrix<T>::vector_mult (NumericVector<T>& dest,
				   const NumericVector<T>& arg) const
{
  dest.zero();
  this->vector_mult_add(dest,arg);
}

//.........这里部分代码省略.........
  //----------- quantities (at dof objects) ------------------------------
  vector < vector < double > >     sol_eldofs(n_unknowns);
  for(int k=0; k<n_unknowns; k++)  sol_eldofs[k].reserve(max_size);
  
//************** act flag (at dof objects) **************************** 
   std::string act_flag_name = "act_flag";
   unsigned int solIndex_act_flag = ml_sol->GetIndex(act_flag_name.c_str());
   unsigned int solFEType_act_flag = ml_sol->GetSolutionType(solIndex_act_flag); 
      if(sol->GetSolutionTimeOrder(solIndex_act_flag) == 2) {
        *(sol->_SolOld[solIndex_act_flag]) = *(sol->_Sol[solIndex_act_flag]);
      }

  //********* variables for ineq constraints (at dof objects) *****************
  const int ineq_flag = INEQ_FLAG;
  const double c_compl = C_COMPL;
  vector < double/*int*/ >  sol_actflag;   sol_actflag.reserve(max_size); //flag for active set
  vector < double >          ctrl_lower;    ctrl_lower.reserve(max_size);
  vector < double >          ctrl_upper;    ctrl_upper.reserve(max_size);
      
  //------------ quantities (at quadrature points) ---------------------
            vector<double>        sol_qp(n_unknowns);
    vector< vector<double> > sol_grad_qp(n_unknowns);
    
      std::fill(sol_qp.begin(), sol_qp.end(), 0.);
    for (unsigned  k = 0; k < n_unknowns; k++) {
        sol_grad_qp[k].resize(dim);
        std::fill(sol_grad_qp[k].begin(), sol_grad_qp[k].end(), 0.);
    }

      
  //******* EQUATION RELATED STUFF ********************************************  
  
 int m_b_f[n_unknowns][n_unknowns];
     m_b_f[pos_state][pos_state] = 1; //nonzero
     m_b_f[pos_state][pos_ctrl]  = 0; //THIS IS ZERO IN NON-LIFTING APPROACHES (there are also pieces that go on top on blocks that are already meant to be different from zero)
     m_b_f[pos_state][pos_adj]   = 1; //nonzero
     m_b_f[pos_state][pos_mu]    = 0;  //this is zero without state constraints
     
     m_b_f[pos_ctrl][pos_state]  = 0;//THIS IS ZERO IN NON-LIFTING APPROACHES
     m_b_f[pos_ctrl][pos_ctrl]   = 1;//nonzero
     m_b_f[pos_ctrl][pos_adj]    = 1;//nonzero
     m_b_f[pos_ctrl][pos_mu]     = 1;//nonzero
     
     m_b_f[pos_adj][pos_state]  = 1; //nonzero
     m_b_f[pos_adj][pos_ctrl]   = 1; //nonzero
     m_b_f[pos_adj][pos_adj]    = 0; //this must always be zero 
     m_b_f[pos_adj][pos_mu]     = 0; //this must always be zero 
     
     m_b_f[pos_mu][pos_state]  = 0; //this is zero without state constraints
     m_b_f[pos_mu][pos_ctrl]   = 1; //nonzero
     m_b_f[pos_mu][pos_adj]    = 0; //this must always be zero 
     m_b_f[pos_mu][pos_mu]     = 1; //nonzero
 
  //***************************************************  
  vector < vector < int > > L2G_dofmap(n_unknowns);     for(int i = 0; i < n_unknowns; i++) { L2G_dofmap[i].reserve(max_size); }
            vector< int >   L2G_dofmap_AllVars; L2G_dofmap_AllVars.reserve( n_unknowns*max_size );
            vector< double >         Res;                      Res.reserve( n_unknowns*max_size );
            vector< double >         Jac;                      Jac.reserve( n_unknowns*max_size * n_unknowns*max_size);
  //***************************************************  
    
    
  //********************* DATA ************************ 
  double alpha = ALPHA_CTRL_VOL;
  double beta  = BETA_CTRL_VOL;
  double penalty_strong = 10e+14;
 //***************************************************  

// Here we compute the residual R(x) = K(x)*x - f. The current solution
// x is passed in the soln vector
void compute_residual (const NumericVector<Number>& soln,
                       NumericVector<Number>& residual,
                       NonlinearImplicitSystem& sys)
{
    EquationSystems &es = *_equation_system;

    // Get a constant reference to the mesh object.
    const MeshBase& mesh = es.get_mesh();

    // The dimension that we are running
    const unsigned int dim = mesh.mesh_dimension();
    libmesh_assert_equal_to (dim, 2);

    // Get a reference to the NonlinearImplicitSystem we are solving
    NonlinearImplicitSystem& system =
        es.get_system<NonlinearImplicitSystem>("Laplace-Young");

    // A reference to the \p DofMap object for this system.  The \p DofMap
    // object handles the index translation from node and element numbers
    // to degree of freedom numbers.  We will talk more about the \p DofMap
    // in future examples.
    const DofMap& dof_map = system.get_dof_map();

    // Get a constant reference to the Finite Element type
    // for the first (and only) variable in the system.
    FEType fe_type = dof_map.variable_type(0);

    // Build a Finite Element object of the specified type.  Since the
    // \p FEBase::build() member dynamically creates memory we will
    // store the object as an \p AutoPtr<FEBase>.  This can be thought
    // of as a pointer that will clean up after itself.
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));

    // A 5th order Gauss quadrature rule for numerical integration.
    QGauss qrule (dim, FIFTH);

    // Tell the finite element object to use our quadrature rule.
    fe->attach_quadrature_rule (&qrule);

    // Declare a special finite element object for
    // boundary integration.
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));

    // Boundary integration requires one quadraure rule,
    // with dimensionality one less than the dimensionality
    // of the element.
    QGauss qface(dim-1, FIFTH);

    // Tell the finte element object to use our
    // quadrature rule.
    fe_face->attach_quadrature_rule (&qface);

    // Here we define some references to cell-specific data that
    // will be used to assemble the linear system.
    // We begin with the element Jacobian * quadrature weight at each
    // integration point.
    const std::vector<Real>& JxW = fe->get_JxW();

    // The element shape functions evaluated at the quadrature points.
    const std::vector<std::vector<Real> >& phi = fe->get_phi();

    // The element shape function gradients evaluated at the quadrature
    // points.
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();

    // Define data structures to contain the resdual contributions
    DenseVector<Number> Re;

    // This vector will hold the degree of freedom indices for
    // the element.  These define where in the global system
    // the element degrees of freedom get mapped.
    std::vector<unsigned int> dof_indices;

    // Now we will loop over all the active elements in the mesh which
    // are local to this processor.
    // We will compute the element residual.
    residual.zero();

    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

    for ( ; el != end_el; ++el)
    {
        // Store a pointer to the element we are currently
        // working on.  This allows for nicer syntax later.
        const Elem* elem = *el;

        // Get the degree of freedom indices for the
        // current element.  These define where in the global
        // matrix and right-hand-side this element will
        // contribute to.
        dof_map.dof_indices (elem, dof_indices);

        // Compute the element-specific data for the current
        // element.  This involves computing the location of the
        // quadrature points (q_point) and the shape functions
        // (phi, dphi) for the current element.
        fe->reinit (elem);
//.........这里部分代码省略.........

void AssembleOptimization::assemble_A_and_F()
{
  A_matrix->zero();
  F_vector->zero();

  const MeshBase & mesh = _sys.get_mesh();

  const unsigned int dim = mesh.mesh_dimension();
  const unsigned int u_var = _sys.variable_number ("u");

  const DofMap & dof_map = _sys.get_dof_map();
  FEType fe_type = dof_map.variable_type(u_var);
  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
  QGauss qrule (dim, fe_type.default_quadrature_order());
  fe->attach_quadrature_rule (&qrule);

  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<Real> > & phi = fe->get_phi();
  const std::vector<std::vector<RealGradient> > & dphi = fe->get_dphi();

  std::vector<dof_id_type> dof_indices;

  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      const Elem * elem = *el;

      dof_map.dof_indices (elem, dof_indices);

      const unsigned int n_dofs = dof_indices.size();

      fe->reinit (elem);

      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);

      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          for (unsigned int dof_i=0; dof_i<n_dofs; dof_i++)
            {
              for (unsigned int dof_j=0; dof_j<n_dofs; dof_j++)
                {
                  Ke(dof_i, dof_j) += JxW[qp] * (dphi[dof_j][qp]* dphi[dof_i][qp]);
                }
              Fe(dof_i) += JxW[qp] * phi[dof_i][qp];
            }
        }

      // This will zero off-diagonal entries of Ke corresponding to
      // Dirichlet dofs.
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

      // We want the diagonal of constrained dofs to be zero too
      for (unsigned int local_dof_index=0; local_dof_index<n_dofs; local_dof_index++)
        {
          dof_id_type global_dof_index = dof_indices[local_dof_index];
          if (dof_map.is_constrained_dof(global_dof_index))
            {
              Ke(local_dof_index, local_dof_index) = 0.;
            }
        }

      A_matrix->add_matrix (Ke, dof_indices);
      F_vector->add_vector (Fe, dof_indices);
    }

  A_matrix->close();
  F_vector->close();
}

  /**
   * Evaluate the residual of the nonlinear system.
   */
  virtual void residual (const NumericVector<Number> & soln,
                         NumericVector<Number> & residual,
                         NonlinearImplicitSystem & /*sys*/)
  {
    const Real young_modulus = es.parameters.get<Real>("young_modulus");
    const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
    const Real forcing_magnitude = es.parameters.get<Real>("forcing_magnitude");

    const MeshBase & mesh = es.get_mesh();
    const unsigned int dim = mesh.mesh_dimension();

    NonlinearImplicitSystem & system =
      es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");

    const unsigned int u_var = system.variable_number ("u");

    const DofMap & dof_map = system.get_dof_map();

    FEType fe_type = dof_map.variable_type(u_var);
    std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);

    std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
    QGauss qface (dim-1, fe_type.default_quadrature_order());
    fe_face->attach_quadrature_rule (&qface);

    const std::vector<Real> & JxW = fe->get_JxW();
    const std::vector<std::vector<Real>> & phi = fe->get_phi();
    const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

    DenseVector<Number> Re;

    DenseSubVector<Number> Re_var[3] =
      {DenseSubVector<Number>(Re),
       DenseSubVector<Number>(Re),
       DenseSubVector<Number>(Re)};

    std::vector<dof_id_type> dof_indices;
    std::vector<std::vector<dof_id_type>> dof_indices_var(3);

    residual.zero();

    for (const auto & elem : mesh.active_local_element_ptr_range())
      {
        dof_map.dof_indices (elem, dof_indices);
        for (unsigned int var=0; var<3; var++)
          dof_map.dof_indices (elem, dof_indices_var[var], var);

        const unsigned int n_dofs = dof_indices.size();
        const unsigned int n_var_dofs = dof_indices_var[0].size();

        fe->reinit (elem);

        Re.resize (n_dofs);
        for (unsigned int var=0; var<3; var++)
          Re_var[var].reposition (var*n_var_dofs, n_var_dofs);

        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            DenseVector<Number> u_vec(3);
            DenseMatrix<Number> grad_u(3, 3);
            for (unsigned int var_i=0; var_i<3; var_i++)
              {
                for (unsigned int j=0; j<n_var_dofs; j++)
                  u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);

                // Row is variable u, v, or w column is x, y, or z
                for (unsigned int var_j=0; var_j<3; var_j++)
                  for (unsigned int j=0; j<n_var_dofs; j++)
                    grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
              }

            DenseMatrix<Number> strain_tensor(3, 3);
            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                {
                  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));

                  for (unsigned int k=0; k<3; k++)
                    strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
                }

            // Define the deformation gradient
            DenseMatrix<Number> F(3, 3);
            F = grad_u;
            for (unsigned int var=0; var<3; var++)
              F(var, var) += 1.;

            DenseMatrix<Number> stress_tensor(3, 3);

            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                for (unsigned int k=0; k<3; k++)
                  for (unsigned int l=0; l<3; l++)
                    stress_tensor(i,j) +=
                      elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k,l);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  std::fill(sol_qp.begin(), sol_qp.end(), 0.);
  std::fill(sol_old_qp.begin(), sol_old_qp.end(), 0.);
  for (unsigned  k = 0; k < n_unknowns; k++) {
        sol_grad_qp[k].resize(dim);
        std::fill(sol_grad_qp[k].begin(), sol_grad_qp[k].end(), 0.);
        sol_old_grad_qp[k].resize(dim);
        std::fill(sol_old_grad_qp[k].begin(), sol_old_grad_qp[k].end(), 0.);
    }

  //----------- quantities (at dof objects) ------------------------------
  vector < vector < double > >     sol_eldofs(n_unknowns);
  vector < vector < double > >     sol_old_eldofs(n_unknowns);
  for(int k=0; k<n_unknowns; k++) { sol_eldofs[k].reserve(max_size);
                                sol_old_eldofs[k].reserve(max_size);
   }
  
  //******** linear system *******************************************  
  vector < vector < int > > L2G_dofmap(n_unknowns);     for(int i = 0; i < n_unknowns; i++) { L2G_dofmap[i].reserve(max_size); }
            vector< int >   L2G_dofmap_AllVars; L2G_dofmap_AllVars.reserve( n_unknowns*max_size );
          
  vector< vector< double > > Res_el(n_unknowns);
  vector< vector< vector< double > > > Jac_el(n_unknowns);

   for(int i = 0; i < n_unknowns; i++) Res_el[i].reserve(max_size);

   for(int i = 0; i < n_unknowns; i++) {
     Jac_el[i].resize(n_unknowns);
     for(int j = 0; j < n_unknowns; j++) {
        Jac_el[i][j].reserve(max_size*max_size);
     }
   }
  
  // Set to zero all the entries of the matrix
  RES->zero();
  JAC->zero();
  
  
  const double deltat_term = 1.;
  const double lapl_term = 1.;
  const double delta_g_term = 1.;

  
  // *** element loop ***
  for (int iel = msh->_elementOffset[iproc]; iel < msh->_elementOffset[iproc+1]; iel++) {

    short unsigned ielGeom = msh->GetElementType(iel);    // element geometry type
    
    //******************** GEOMETRY ********************* 
    unsigned nDofx = msh->GetElementDofNumber(iel, coords_fe_type);       // number of coordinate element dofs
    for (int i = 0; i < dim; i++)  coords_at_dofs[i].resize(nDofx);
    // local storage of coordinates
    for (unsigned i = 0; i < nDofx; i++) {
      unsigned xDof  = msh->GetSolutionDof(i, iel, coords_fe_type);       // global to global mapping between coordinates node and coordinate dof

      for (unsigned jdim = 0; jdim < dim; jdim++) {
        coords_at_dofs[jdim][i] = (*msh->_topology->_Sol[jdim])(xDof);    // global extraction and local storage for the element coordinates
      }
    }
  //***************************************************  
    
  //all vars###################################################################  
  for (unsigned  k = 0; k < n_unknowns; k++) {
    unsigned  ndofs_unk = msh->GetElementDofNumber(iel, SolFEType[k]);
	   Sol_n_el_dofs[k] = ndofs_unk;
          sol_eldofs[k].resize(ndofs_unk);
      sol_old_eldofs[k].resize(ndofs_unk);