C++ EquationSystems::get_mesh方法代码示例

void assemble_poisson(EquationSystems & es,
                      const std::string & system_name)
{
  libmesh_assert_equal_to (system_name, "Poisson");

  const MeshBase & mesh = es.get_mesh();
  const unsigned int dim = mesh.mesh_dimension();
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");

  const DofMap & dof_map = system.get_dof_map();

  FEType fe_type = dof_map.variable_type(0);
  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
  QGauss qrule (dim, FIFTH);
  fe->attach_quadrature_rule (&qrule);
  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  QGauss qface(dim-1, FIFTH);
  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();

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

  std::vector<dof_id_type> dof_indices;

  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);

      fe->reinit (elem);

      Ke.resize (dof_indices.size(),
                 dof_indices.size());

      Fe.resize (dof_indices.size());

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

      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

      system.matrix->add_matrix (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);
    }
}

void transform_mesh_and_plot(EquationSystems & es,
                             Real curvature,
                             const std::string & filename)
{
  // Loop over the mesh nodes and move them!
  MeshBase & mesh = es.get_mesh();

  MeshBase::node_iterator       node_it  = mesh.nodes_begin();
  const MeshBase::node_iterator node_end = mesh.nodes_end();

  for ( ; node_it != node_end; node_it++)
    {
      Node * node = *node_it;

      Real x = (*node)(0);
      Real z = (*node)(2);

      (*node)(0) = -1./curvature + (1./curvature + x)*cos(curvature*z);
      (*node)(2) = (1./curvature + x)*sin(curvature*z);
    }

#ifdef LIBMESH_HAVE_EXODUS_API
  ExodusII_IO(mesh).write_equation_systems(filename, es);
#endif
}

void write_output_solvedata(EquationSystems & es,
                            unsigned int a_step,
                            unsigned int newton_steps,
                            unsigned int krylov_steps,
                            unsigned int tv_sec,
                            unsigned int tv_usec)
{
    MeshBase & mesh = es.get_mesh();
    unsigned int n_active_elem = mesh.n_active_elem();
    unsigned int n_active_dofs = es.n_active_dofs();

    if (mesh.processor_id() == 0)
    {
        // Write out the number of elements/dofs used
        std::ofstream activemesh ("out_activemesh.m",
                                  std::ios_base::app | std::ios_base::out);
        activemesh.precision(17);
        activemesh << (a_step + 1) << ' '
                   << n_active_elem << ' '
                   << n_active_dofs << std::endl;

        // Write out the number of solver steps used
        std::ofstream solvesteps ("out_solvesteps.m",
                                  std::ios_base::app | std::ios_base::out);
        solvesteps.precision(17);
        solvesteps << newton_steps << ' '
                   << krylov_steps << std::endl;

        // Write out the clock time
        std::ofstream clocktime ("out_clocktime.m",
                                 std::ios_base::app | std::ios_base::out);
        clocktime.precision(17);
        clocktime << tv_sec << '.' << tv_usec << std::endl;
    }
}

void MeshOutput<MT>::write_equation_systems (const std::string & fname,
                                             const EquationSystems & es,
                                             const std::set<std::string> * system_names)
{
  START_LOG("write_equation_systems()", "MeshOutput");

  // We may need to gather and/or renumber a ParallelMesh to output
  // it, making that const qualifier in our constructor a dirty lie
  MT & my_mesh = const_cast<MT &>(*_obj);

  // If we're asked to write data that's associated with a different
  // mesh, output files full of garbage are the result.
  libmesh_assert_equal_to(&es.get_mesh(), _obj);

  // A non-renumbered mesh may not have a contiguous numbering, and
  // that needs to be fixed before we can build a solution vector.
  if (my_mesh.max_elem_id() != my_mesh.n_elem() ||
      my_mesh.max_node_id() != my_mesh.n_nodes())
    {
      // If we were allowed to renumber then we should have already
      // been properly renumbered...
      libmesh_assert(!my_mesh.allow_renumbering());

      libmesh_do_once(libMesh::out <<
                      "Warning:  This MeshOutput subclass only supports meshes which are contiguously renumbered!"
                      << std::endl;);

void scale_mesh_and_plot(EquationSystems & es,
                         const RBParameters & mu,
                         const std::string & filename)
{
  // Loop over the mesh nodes and move them!
  MeshBase & mesh = es.get_mesh();

  MeshBase::node_iterator       node_it  = mesh.nodes_begin();
  const MeshBase::node_iterator node_end = mesh.nodes_end();

  for ( ; node_it != node_end; node_it++)
    {
      Node * node = *node_it;

      (*node)(0) *= mu.get_value("x_scaling");
    }

  // Post-process the solution to compute the stresses
  compute_stresses(es);

#ifdef LIBMESH_HAVE_EXODUS_API
  ExodusII_IO (mesh).write_equation_systems (filename, es);
#endif

  // Loop over the mesh nodes and move them!
  node_it = mesh.nodes_begin();

  for ( ; node_it != node_end; node_it++)
    {
      Node * node = *node_it;

      (*node)(0) /= mu.get_value("x_scaling");
    }
}

void run_timestepping(EquationSystems& systems, GetPot& args)
{
    TransientExplicitSystem& aux_system = systems.get_system<TransientExplicitSystem>("auxiliary");

    SolidSystem& solid_system = systems.get_system<SolidSystem>("solid");

    AutoPtr<VTKIO> io = AutoPtr<VTKIO>(new VTKIO(systems.get_mesh()));

    Real duration = args("duration", 1.0);

    for (unsigned int t_step = 0; t_step < duration/solid_system.deltat; t_step++) {
        // Progress in current phase [0..1]
        Real progress = t_step * solid_system.deltat / duration;
        systems.parameters.set<Real>("progress") = progress;
        systems.parameters.set<unsigned int>("step") = t_step;

        // Update message

        out << "===== Time Step " << std::setw(4) << t_step;
        out << " (" << std::fixed << std::setprecision(2) << std::setw(6) << (progress*100.) << "%)";
        out << ", time = " << std::setw(7) << solid_system.time;
        out << " =====" << std::endl;

        // Advance timestep in auxiliary system
        aux_system.current_local_solution->close();
        aux_system.old_local_solution->close();
        *aux_system.older_local_solution = *aux_system.old_local_solution;
        *aux_system.old_local_solution = *aux_system.current_local_solution;

        out << "Solving Solid" << std::endl;
        solid_system.solve();
        aux_system.reinit();

        // Carry out the adaptive mesh refinement/coarsening
        out << "Doing a reinit of the equation systems" << std::endl;
        systems.reinit();

        if (t_step % args("output/frequency", 1) == 0) {
            std::stringstream file_name;
            file_name << args("results_directory", "./") << "fem_";
            file_name << std::setw(6) << std::setfill('0') << t_step;
            file_name << ".pvtu";

            io->write_equation_systems(file_name.str(), systems);
        }
        // Advance to the next timestep in a transient problem
        out << "Advancing to next step" << std::endl;
        solid_system.time_solver->advance_timestep();
    }
}

/*
 * FIXME: This is known to write nonsense on AMR meshes
 * and it strips the imaginary parts of complex Numbers
 */
void VTKIO::system_vectors_to_vtk(const EquationSystems& es, vtkUnstructuredGrid*& grid)
{
  if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
    {
      std::map<std::string, std::vector<Number> > vecs;
      for (unsigned int i=0; i<es.n_systems(); ++i)
        {
          const System& sys = es.get_system(i);
          System::const_vectors_iterator v_end = sys.vectors_end();
          System::const_vectors_iterator it = sys.vectors_begin();
          for (; it!= v_end; ++it)
            {
              // for all vectors on this system
              std::vector<Number> values;
              // libMesh::out<<"it "<<it->first<<std::endl;

              it->second->localize_to_one(values, 0);
              // libMesh::out<<"finish localize"<<std::endl;
              vecs[it->first] = values;
            }
        }

      std::map<std::string, std::vector<Number> >::iterator it = vecs.begin();

      for (; it!=vecs.end(); ++it)
        {
          vtkDoubleArray *data = vtkDoubleArray::New();
          data->SetName(it->first.c_str());
          libmesh_assert_equal_to (it->second.size(), es.get_mesh().n_nodes());
          data->SetNumberOfValues(it->second.size());

          for (unsigned int i=0; i<it->second.size(); ++i)
            {
#ifdef LIBMESH_USE_COMPLEX_NUMBERS
              libmesh_do_once (libMesh::err << "Only writing the real part for complex numbers!\n"
                               << "if you need this support contact " << LIBMESH_PACKAGE_BUGREPORT
                               << std::endl);
              data->SetValue(i, it->second[i].real());
#else
              data->SetValue(i, it->second[i]);
#endif

            }
          grid->GetPointData()->AddArray(data);
          data->Delete();
        }
    }
}

void write_output(EquationSystems &es,		                   
		  unsigned int a_step, // The adaptive step count
		  std::string solution_type) // primal or adjoint solve
{
  MeshBase &mesh = es.get_mesh();

#ifdef LIBMESH_HAVE_GMV
  
  OStringStream file_name_gmv;
  file_name_gmv << solution_type << ".out.gmv.";      
  OSSRealzeroright(file_name_gmv,2,0,a_step);
  
  GMVIO(mesh).write_equation_systems
    (file_name_gmv.str(), es);
    
#endif    
}

void write_output(EquationSystems &es,
                  unsigned int a_step,       // The adaptive step count
                  std::string solution_type) // primal or adjoint solve
{
#ifdef LIBMESH_HAVE_GMV
  MeshBase &mesh = es.get_mesh();

  std::ostringstream file_name_gmv;
  file_name_gmv << solution_type
                << ".out.gmv."
                << std::setw(2)
                << std::setfill('0')
                << std::right
                << a_step;

  GMVIO(mesh).write_equation_systems
    (file_name_gmv.str(), es);
#endif
}

void assemble_stokes (EquationSystems& es,
                      const std::string& system_name)
{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Stokes");

  // 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();

  // Get a reference to the Convection-Diffusion system object.
  LinearImplicitSystem & system =
    es.get_system<LinearImplicitSystem> ("Stokes");

  // Numeric ids corresponding to each variable in the system
  const unsigned int u_var = system.variable_number ("u");
  const unsigned int v_var = system.variable_number ("v");
  const unsigned int p_var = system.variable_number ("p");

  // Get the Finite Element type for "u".  Note this will be
  // the same as the type for "v".
  FEType fe_vel_type = system.variable_type(u_var);

  // Get the Finite Element type for "p".
  FEType fe_pres_type = system.variable_type(p_var);

  // Build a Finite Element object of the specified type for
  // the velocity variables.
  UniquePtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));

  // Build a Finite Element object of the specified type for
  // the pressure variables.
  UniquePtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));

  // A Gauss quadrature rule for numerical integration.
  // Let the \p FEType object decide what order rule is appropriate.
  QGauss qrule (dim, fe_vel_type.default_quadrature_order());

  // Tell the finite element objects to use our quadrature rule.
  fe_vel->attach_quadrature_rule (&qrule);
  fe_pres->attach_quadrature_rule (&qrule);

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

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

  // The element shape functions for the pressure variable
  // evaluated at the quadrature points.
  const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();

  // 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();

  // Define data structures to contain the element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke" and "Fe".
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  DenseSubMatrix<Number>
    Kuu(Ke), Kuv(Ke), Kup(Ke),
    Kvu(Ke), Kvv(Ke), Kvp(Ke),
    Kpu(Ke), Kpv(Ke), Kpp(Ke);

  DenseSubVector<Number>
    Fu(Fe),
    Fv(Fe),
    Fp(Fe);

  // 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<dof_id_type> dof_indices;
  std::vector<dof_id_type> dof_indices_u;
  std::vector<dof_id_type> dof_indices_v;
  std::vector<dof_id_type> dof_indices_p;

  // Now we will loop over all the elements in the mesh that
  // live on the local processor. We will compute the element
  // matrix and right-hand-side contribution.  In case users later
  // modify this program to include refinement, we will be safe and
  // will only consider the active elements; hence we use a variant of
  // the \p active_elem_iterator.

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

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

void
MultiAppProjectionTransfer::assembleL2From(EquationSystems & es, const std::string & system_name)
{
  unsigned int n_apps = _multi_app->numGlobalApps();
  std::vector<NumericVector<Number> *> from_slns(n_apps, NULL);
  std::vector<MeshFunction *> from_fns(n_apps, NULL);
  std::vector<MeshTools::BoundingBox *> from_bbs(n_apps, NULL);

  // get bounding box, mesh function and solution for each subapp
  for (unsigned int i = 0; i < n_apps; i++)
  {
    if (!_multi_app->hasLocalApp(i))
      continue;

    MPI_Comm swapped = Moose::swapLibMeshComm(_multi_app->comm());

    FEProblem & from_problem = *_multi_app->appProblem(i);
    EquationSystems & from_es = from_problem.es();
    MeshBase & from_mesh = from_es.get_mesh();
    MeshTools::BoundingBox * app_box = new MeshTools::BoundingBox(MeshTools::processor_bounding_box(from_mesh, from_mesh.processor_id()));
    from_bbs[i] = app_box;

    MooseVariable & from_var = from_problem.getVariable(0, _from_var_name);
    System & from_sys = from_var.sys().system();
    unsigned int from_var_num = from_sys.variable_number(from_var.name());

    NumericVector<Number> * serialized_from_solution = NumericVector<Number>::build(from_sys.comm()).release();
    serialized_from_solution->init(from_sys.n_dofs(), false, SERIAL);
    // Need to pull down a full copy of this vector on every processor so we can get values in parallel
    from_sys.solution->localize(*serialized_from_solution);
    from_slns[i] = serialized_from_solution;

    MeshFunction * from_func = new MeshFunction(from_es, *serialized_from_solution, from_sys.get_dof_map(), from_var_num);
    from_func->init(Trees::ELEMENTS);
    from_func->enable_out_of_mesh_mode(NOTFOUND);
    from_fns[i] = from_func;

    Moose::swapLibMeshComm(swapped);
  }


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

  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);

  FEType fe_type = system.variable_type(0);
  AutoPtr<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<Point> & xyz = fe->get_xyz();

  const DofMap& dof_map = system.get_dof_map();
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
  std::vector<dof_id_type> dof_indices;

  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;

    fe->reinit (elem);

    dof_map.dof_indices (elem, dof_indices);
    Ke.resize (dof_indices.size(), dof_indices.size());
    Fe.resize (dof_indices.size());

    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
    {
      Point qpt = xyz[qp];
      Real f = 0.;
      for (unsigned int app = 0; app < n_apps; app++)
      {
        Point pt = qpt - _multi_app->position(app);
        if (from_bbs[app] != NULL && from_bbs[app]->contains_point(pt))
        {
          MPI_Comm swapped = Moose::swapLibMeshComm(_multi_app->comm());
          f = (*from_fns[app])(pt);
          Moose::swapLibMeshComm(swapped);
          break;
        }
      }

      // Now compute the element matrix and RHS contributions.
      for (unsigned int i=0; i<phi.size(); i++)
      {
        // RHS
        Fe(i) += JxW[qp] * (f * phi[i][qp]);

        if (_compute_matrix)
          for (unsigned int j = 0; j < phi.size(); j++)
          {
            // The matrix contribution
            Ke(i,j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
          }
      }
//.........这里部分代码省略.........

void
MultiAppProjectionTransfer::assembleL2To(EquationSystems & es, const std::string & system_name)
{
  unsigned int app = es.parameters.get<unsigned int>("app");

  FEProblem & from_problem = *_multi_app->problem();
  EquationSystems & from_es = from_problem.es();

  MooseVariable & from_var = from_problem.getVariable(0, _from_var_name);
  System & from_sys = from_var.sys().system();
  unsigned int from_var_num = from_sys.variable_number(from_var.name());

  NumericVector<Number> * serialized_from_solution = NumericVector<Number>::build(from_sys.comm()).release();
  serialized_from_solution->init(from_sys.n_dofs(), false, SERIAL);
  // Need to pull down a full copy of this vector on every processor so we can get values in parallel
  from_sys.solution->localize(*serialized_from_solution);

  MeshFunction from_func(from_es, *serialized_from_solution, from_sys.get_dof_map(), from_var_num);
  from_func.init(Trees::ELEMENTS);
  from_func.enable_out_of_mesh_mode(0.);


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

  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);

  FEType fe_type = system.variable_type(0);
  AutoPtr<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<Point> & xyz = fe->get_xyz();

  const DofMap& dof_map = system.get_dof_map();
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
  std::vector<dof_id_type> dof_indices;

  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;

    fe->reinit (elem);

    dof_map.dof_indices (elem, dof_indices);
    Ke.resize (dof_indices.size(), dof_indices.size());
    Fe.resize (dof_indices.size());

    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
    {
      Point qpt = xyz[qp];
      Point pt = qpt + _multi_app->position(app);
      Real f = from_func(pt);

      // Now compute the element matrix and RHS contributions.
      for (unsigned int i=0; i<phi.size(); i++)
      {
        // RHS
        Fe(i) += JxW[qp] * (f * phi[i][qp]);

        if (_compute_matrix)
          for (unsigned int j = 0; j < phi.size(); j++)
          {
            // The matrix contribution
            Ke(i,j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
          }
      }
      dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);

      if (_compute_matrix)
        system.matrix->add_matrix(Ke, dof_indices);
      system.rhs->add_vector(Fe, dof_indices);
    }
  }
}

// We now define the matrix assembly function for the
// Poisson system.  We need to first compute element
// matrices and right-hand sides, and then take into
// account the boundary conditions.
void assemble_poisson(EquationSystems& es,
                      const std::string& system_name)
{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Poisson");


  // Declare a performance log.  Give it a descriptive
  // string to identify what part of the code we are
  // logging, since there may be many PerfLogs in an
  // application.
  PerfLog perf_log ("Matrix Assembly");

  // 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();

  // Get a reference to the LinearImplicitSystem we are solving
  LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson");

  // 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 UniquePtr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.
  UniquePtr<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.
  UniquePtr<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 physical XY locations of the quadrature points on the element.
  // These might be useful for evaluating spatially varying material
  // properties at the quadrature points.
  const std::vector<Point>& q_point = fe->get_xyz();

  // 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 element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke" and "Fe". More detail is in example 3.
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  // 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<dof_id_type> dof_indices;

  // Now we will loop over all the elements in the mesh.
  // We will compute the element matrix and right-hand-side
  // contribution.  See example 3 for a discussion of the
  // element iterators.
  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)
//.........这里部分代码省略.........

void
MultiAppProjectionTransfer::assembleL2(EquationSystems & es, const std::string & system_name)
{
  // Get the system and mesh from the input arguments.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);
  MeshBase & to_mesh = es.get_mesh();

  // Get the meshfunction evaluations and the map that was stashed in the es.
  std::vector<Real> & final_evals = * es.parameters.get<std::vector<Real>*>("final_evals");
  std::map<unsigned int, unsigned int> & element_map = * es.parameters.get<std::map<unsigned int, unsigned int>*>("element_map");

  // Setup system vectors and matrices.
  FEType fe_type = system.variable_type(0);
  std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
  QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
  fe->attach_quadrature_rule(&qrule);
  const DofMap& dof_map = system.get_dof_map();
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
  std::vector<dof_id_type> dof_indices;
  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<Real> > & phi = fe->get_phi();

  const MeshBase::const_element_iterator end_el = to_mesh.active_local_elements_end();
  for (MeshBase::const_element_iterator el = to_mesh.active_local_elements_begin(); el != end_el; ++el)
  {
    const Elem* elem = *el;
    fe->reinit (elem);

    dof_map.dof_indices (elem, dof_indices);
    Ke.resize (dof_indices.size(), dof_indices.size());
    Fe.resize (dof_indices.size());

    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
    {
      Real meshfun_eval = 0.;
      if (element_map.find(elem->id()) != element_map.end())
      {
        // We have evaluations for this element.
        meshfun_eval = final_evals[element_map[elem->id()] + qp];
      }

      // Now compute the element matrix and RHS contributions.
      for (unsigned int i=0; i<phi.size(); i++)
      {
        // RHS
        Fe(i) += JxW[qp] * (meshfun_eval * phi[i][qp]);

        if (_compute_matrix)
          for (unsigned int j = 0; j < phi.size(); j++)
          {
            // The matrix contribution
            Ke(i,j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
          }
      }
      dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);

      if (_compute_matrix)
        system.matrix->add_matrix(Ke, dof_indices);
      system.rhs->add_vector(Fe, dof_indices);
    }
  }
}

void assemble_matrices(EquationSystems & es,
                       const std::string & libmesh_dbg_var(system_name))
{

  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Eigensystem");

#ifdef LIBMESH_HAVE_SLEPC

  // 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();

  // Get a reference to our system.
  EigenSystem & eigen_system = es.get_system<EigenSystem> ("Eigensystem");

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

  // A reference to the two system matrices
  SparseMatrix<Number> & matrix_A = *eigen_system.matrix_A;
  SparseMatrix<Number> & matrix_B = *eigen_system.matrix_B;

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

  // A  Gauss quadrature rule for numerical integration.
  // Use the default quadrature order.
  QGauss qrule (dim, fe_type.default_quadrature_order());

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

  // 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();

  // A reference to the DofMap object for this system.  The DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.
  const DofMap & dof_map = eigen_system.get_dof_map();

  // The element mass and stiffness matrices.
  DenseMatrix<Number> Me;
  DenseMatrix<Number> Ke;

  // 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<dof_id_type> dof_indices;


  // Now we will loop over all the elements in the mesh that
  // live on the local processor. We will compute the element
  // matrix and right-hand-side contribution.  In case users
  // later modify this program to include refinement, we will
  // be safe and will only consider the active elements;
  // hence we use a variant of the active_elem_iterator.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // 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);

      // Zero the element matrices before
      // summing them.  We use the resize member here because
      // the number of degrees of freedom might have changed from
      // the last element.  Note that this will be the case if the
      // element type is different (i.e. the last element was a
      // triangle, now we are on a quadrilateral).
      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());
      Ke.resize (n_dofs, n_dofs);
      Me.resize (n_dofs, n_dofs);

      // Now loop over the quadrature points.  This handles
      // the numeric integration.
      //
      // We will build the element matrix.  This involves
//.........这里部分代码省略.........