C++ EquationSystems类代码示例

/**
 * everything that is identical for the systems, and
 * should _not_ go into EquationSystems::compare(),
 * can go in this do_compare().
 */
bool do_compare (EquationSystems & les,
                 EquationSystems & res,
                 double threshold,
                 bool verbose)
{

    if (verbose)
    {
        libMesh::out        << "*********   LEFT SYSTEM    *********" << std::endl;
        les.print_info  ();
        libMesh::out << "*********   RIGHT SYSTEM   *********" << std::endl;
        res.print_info ();
        libMesh::out << "********* COMPARISON PHASE *********" << std::endl
                     << std::endl;
    }

    /**
     * start comparing
     */
    bool result = les.compare(res, threshold, verbose);
    if (verbose)
    {
        libMesh::out        << "*********     FINISHED     *********" << std::endl;
    }
    return result;
}

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

bool EquationSystems::compare (const EquationSystems& other_es,
                               const Real threshold,
                               const bool verbose) const
{
  // safety check, whether we handle at least the same number
  // of systems
  std::vector<bool> os_result;

  if (this->n_systems() != other_es.n_systems())
    {
      if (verbose)
        {
          libMesh::out << "  Fatal difference. This system handles "
                       << this->n_systems() << " systems," << std::endl
                       << "  while the other system handles "
                       << other_es.n_systems()
                       << " systems." << std::endl
                       << "  Aborting comparison." << std::endl;
        }
      return false;
    }
  else
    {
      // start comparing each system
      const_system_iterator       pos = _systems.begin();
      const const_system_iterator end = _systems.end();

      for (; pos != end; ++pos)
        {
          const std::string& sys_name = pos->first;
          const System&  system        = *(pos->second);

          // get the other system
          const System& other_system   = other_es.get_system (sys_name);

          os_result.push_back (system.compare (other_system, threshold, verbose));

        }

    }


  // sum up the results
  if (os_result.size()==0)
    return true;
  else
    {
      bool os_identical;
      unsigned int n = 0;
      do
        {
          os_identical = os_result[n];
          n++;
        }
      while (os_identical && n<os_result.size());
      return os_identical;
    }
}

void ExodusII_IO::write_discontinuous_exodusII (const std::string& name,
				     const EquationSystems& es)
{
  std::vector<std::string> solution_names;
  std::vector<Number>      v;

  es.build_variable_names  (solution_names);
  es.build_discontinuous_solution_vector (v);

  this->write_nodal_data_discontinuous(name, v, solution_names);
}

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

CondensedEigenSystem::CondensedEigenSystem (EquationSystems& es,
                                            const std::string& name,
                                            const unsigned int number)
  : Parent(es, name, number),
    condensed_matrix_A(SparseMatrix<Number>::build(es.comm())),
    condensed_matrix_B(SparseMatrix<Number>::build(es.comm())),
    condensed_dofs_initialized(false)
{
}

/**
 * FIXME: This is a default implementation - derived classes should
 * reimplement it for efficiency.
 */
void ErrorEstimator::estimate_errors(const EquationSystems & equation_systems,
                                     ErrorMap & errors_per_cell,
                                     const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
                                     bool estimate_parent_error)
{
  SystemNorm old_error_norm = this->error_norm;

  // Find the requested error values from each system
  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
    {
      const System & sys = equation_systems.get_system(s);

      unsigned int n_vars = sys.n_vars();

      for (unsigned int v = 0; v != n_vars; ++v)
        {
          // Only fill in ErrorVectors the user asks for
          if (errors_per_cell.find(std::make_pair(&sys, v)) ==
              errors_per_cell.end())
            continue;

          // Calculate error in only one variable
          std::vector<Real> weights(n_vars, 0.0);
          weights[v] = 1.0;
          this->error_norm =
            SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
                       weights);

          const NumericVector<Number> * solution_vector = nullptr;
          if (solution_vectors &&
              solution_vectors->find(&sys) != solution_vectors->end())
            solution_vector = solution_vectors->find(&sys)->second;

          this->estimate_error
            (sys, *errors_per_cell[std::make_pair(&sys, v)],
             solution_vector, estimate_parent_error);
        }
    }

  // Restore our old state before returning
  this->error_norm = old_error_norm;
}

void ExodusII_IO::write_element_data (const EquationSystems & es)
{
  // The first step is to collect the element data onto this processor.
  // We want the constant monomial data.

  std::vector<Number> soln;
  std::vector<std::string> names;

  // If _output_variables is populated we need to filter the monomials we output
  if (_output_variables.size())
  {
    std::vector<std::string> monomials;
    const FEType type(CONSTANT, MONOMIAL);
    es.build_variable_names(monomials, &type);

    for (std::vector<std::string>::iterator it = monomials.begin(); it != monomials.end(); ++it)
      if (std::find(_output_variables.begin(), _output_variables.end(), *it) != _output_variables.end())
        names.push_back(*it);
  }

  // If we pass in a list of names to "get_solution" it'll filter the variables comming back
  es.get_solution( soln, names );

  // The data must ultimately be written block by block.  This means that this data
  // must be sorted appropriately.

  if (!exio_helper->created())
    {
      libMesh::err << "ERROR, ExodusII file must be initialized "
                   << "before outputting element variables.\n"
                   << std::endl;
      libmesh_error();
    }

  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();

  exio_helper->initialize_element_variables( mesh, names );
  exio_helper->write_element_values(mesh,soln,_timestep);
}

// ------------------------------------------------------------
// LinearImplicitSystem implementation
LinearImplicitSystem::LinearImplicitSystem (EquationSystems & es,
                                            const std::string & name_in,
                                            const unsigned int number_in) :

  Parent                 (es, name_in, number_in),
  linear_solver          (LinearSolver<Number>::build(es.comm())),
  _n_linear_iterations   (0),
  _final_linear_residual (1.e20),
  _shell_matrix(nullptr),
  _subset(nullptr),
  _subset_solve_mode(SUBSET_ZERO)
{
}

  void testRestart()
  {
    SlitFunc slitfunc;

    _mesh->write("slit_mesh.xda");
    _es->write("slit_solution.xda",
               EquationSystems::WRITE_DATA |
               EquationSystems::WRITE_SERIAL_FILES);

    Mesh mesh2(*TestCommWorld);
    mesh2.read("slit_mesh.xda");
    EquationSystems es2(mesh2);
    es2.read("slit_solution.xda");

    System & sys2 = es2.get_system<System> ("SimpleSystem");

    unsigned int dim = 2;

    CPPUNIT_ASSERT_EQUAL( sys2.n_vars(), 1u );

    FEMContext context(sys2);
    FEBase * fe = NULL;
    context.get_element_fe( 0, fe, dim );
    const std::vector<Point> & xyz = fe->get_xyz();
    fe->get_phi();

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

    for (; el != end_el; ++el)
      {
        const Elem * elem = *el;
        context.pre_fe_reinit(sys2, elem);
        context.elem_fe_reinit();

        const unsigned int n_qp = xyz.size();

        for (unsigned int qp=0; qp != n_qp; ++qp)
          {
            const Number exact_val = slitfunc(context, xyz[qp]);

            const Number discrete_val = context.interior_value(0, qp);

            CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(exact_val),
                                         libmesh_real(discrete_val),
                                         TOLERANCE*TOLERANCE);
          }
      }
  }