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

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

void read_output(EquationSystems & es,
                 unsigned int t_step,
                 unsigned int a_step,
                 std::string solution_type,
                 FEMParameters & param)
{
    MeshBase & mesh = es.get_mesh();

    std::string file_name_mesh, file_name_soln;
    // Look for ASCII files first
    if (param.output_xda)
    {
        file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param);
        file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xda", param);
    }
    else if (param.output_xdr)
    {
        file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param);
        file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param);
    }

    // Read in the mesh
    mesh.read(file_name_mesh);

    // And the stored solution
    es.read(file_name_soln, READ,
            EquationSystems::READ_HEADER |
            EquationSystems::READ_DATA |
            EquationSystems::READ_ADDITIONAL_DATA);

    // Put systems in a consistent state
    for (unsigned int i = 0; i != es.n_systems(); ++i)
        es.get_system<FEMSystem>(i).update();

    // Figure out the current time
    Real current_time = 0., current_timestep = 0.;

    if (param.timesolver_tolerance)
    {
        std::ifstream times ("out_time.m");
        std::ifstream timesteps ("out_timesteps.m");
        if (times.is_open() && timesteps.is_open())
        {
            // Read headers
            const unsigned int headersize = 25;
            char header[headersize];
            timesteps.getline (header, headersize);
            if (strcmp(header, "vector_timesteps = [") != 0)
                libmesh_error_msg("Bad header in out_timesteps.m:\n" << header);

            times.getline (header, headersize);
            if (strcmp(header, "vector_time = [") != 0)
                libmesh_error_msg("Bad header in out_time.m:\n" << header);

            // Read each timestep
            for (unsigned int i = 0; i != t_step; ++i)
            {
                if (!times.good())
                    libmesh_error_msg("Error: File out_time.m is in non-good state.");
                times >> current_time;
                timesteps >> current_timestep;
            }
            // Remember to increment the last timestep; out_times.m
            // lists each *start* time
            current_time += current_timestep;
        }
        else