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

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

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

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

// The matrix assembly function to be called at each time step to
// prepare for the linear solve.
void assemble_solid (EquationSystems& es,
                      const std::string& system_name)
{

//es.print_info();

#if LOG_ASSEMBLE_PERFORMANCE
  PerfLog perf_log("Assemble");
  perf_log.push("assemble stiffness");
#endif

    // Get a reference to the auxiliary system
  //TransientExplicitSystem& aux_system = es.get_system<TransientExplicitSystem>("Newton-update");

  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert (system_name == "Newton-update");
  
  // 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 Stokes system object.
  TransientLinearImplicitSystem & newton_update =
   es.get_system<TransientLinearImplicitSystem> ("Newton-update");

  // New
   TransientLinearImplicitSystem & last_non_linear_soln =
    es.get_system<TransientLinearImplicitSystem> ("Last-non-linear-soln");

 TransientLinearImplicitSystem & fluid_system_vel =
    es.get_system<TransientLinearImplicitSystem> ("fluid-system-vel");

#if VELOCITY
TransientLinearImplicitSystem&  velocity = es.get_system<TransientLinearImplicitSystem>("velocity-system");
#endif

#if UN_MINUS_ONE
TransientLinearImplicitSystem & unm1 =
    es.get_system<TransientLinearImplicitSystem> ("unm1-system");
#endif
test(62);
const System & ref_sys = es.get_system("Reference-Configuration"); 
test(63);
  
  // Numeric ids corresponding to each variable in the system
  const unsigned int u_var = last_non_linear_soln .variable_number ("u");
  const unsigned int v_var = last_non_linear_soln .variable_number ("v");
  const unsigned int w_var = last_non_linear_soln .variable_number ("w");
#if INCOMPRESSIBLE
  const unsigned int p_var = last_non_linear_soln .variable_number ("p");
#endif 



#if FLUID_P_CONST
    const unsigned int m_var = fluid_system_vel.variable_number ("fluid_M");
  std::vector<unsigned int> dof_indices_m;
#endif


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


test(64);

#if INCOMPRESSIBLE
  // Get the Finite Element type for "p".
  FEType fe_pres_type = last_non_linear_soln .variable_type(p_var);
#endif 

  // Build a Finite Element object of the specified type for
  // the velocity variables.
  AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));
  
#if INCOMPRESSIBLE 
  // Build a Finite Element object of the specified type for
  // the pressure variables.
  AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
#endif 
  // 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);
test(65);
//        AutoPtr<QBase> qrule2(fe_vel_type.default_quadrature_rule(dim));
// fe_vel->attach_quadrature_rule (qrule2.get());

#if INCOMPRESSIBLE 
  fe_pres->attach_quadrature_rule (&qrule);
#endif
  // The element Jacobian * quadrature weight at each integration point.   
//.........这里部分代码省略.........