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

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

  void setUp()
  {
    this->build_mesh();

    // libMesh *should* renumber now, or a ParallelMesh might not have
    // contiguous ids, which is a requirement to write xda files.
    _mesh->allow_renumbering(true);

    _es = new EquationSystems(*_mesh);
    _sys = &_es->add_system<System> ("SimpleSystem");
    _sys->add_variable("u", FIRST);

    _es->init();
    SlitFunc slitfunc;
    _sys->project_solution(&slitfunc);

#ifdef LIBMESH_ENABLE_AMR
    MeshRefinement(*_mesh).uniformly_refine(1);
    _es->reinit();
    MeshRefinement(*_mesh).uniformly_refine(1);
    _es->reinit();
#endif
  }

int main (int argc, char** argv)
{
  LibMeshInit init(argc, argv);

  if (argc < 4)
    libMesh::out << "Usage: ./prog -d DIM filename" << std::endl;

  // Variables to get us started
  const unsigned int dim = atoi(argv[2]);

  std::string meshname  (argv[3]);

  // declare a mesh...
  Mesh mesh(init.comm(), dim);

  // Read a mesh
  mesh.read(meshname);

  GMVIO(mesh).write ("out_0.gmv");

  mesh.elem(0)->set_refinement_flag (Elem::REFINE);

  MeshRefinement mesh_refinement (mesh);

  mesh_refinement.refine_and_coarsen_elements ();
  mesh_refinement.uniformly_refine (2);

  mesh.print_info();


  // Set up the equation system(s)
  EquationSystems es (mesh);

  LinearImplicitSystem& primary =
    es.add_system<LinearImplicitSystem>("primary");

  primary.add_variable ("U", FIRST);
  primary.add_variable ("V", FIRST);

  primary.get_dof_map()._dof_coupling->resize(2);
  (*primary.get_dof_map()._dof_coupling)(0,0) = 1;
  (*primary.get_dof_map()._dof_coupling)(1,1) = 1;

  primary.attach_assemble_function(assemble);

  es.init ();

  es.print_info ();
  primary.get_dof_map().print_dof_constraints ();

  // call the solver.
  primary.solve ();

  GMVIO(mesh).write_equation_systems ("out_1.gmv",
                                      es);



  // Refine uniformly
  mesh_refinement.uniformly_refine (1);
  es.reinit ();

  // Write out the projected solution
  GMVIO(mesh).write_equation_systems ("out_2.gmv",
                                      es);

  // Solve again. Output the refined solution
  primary.solve ();
  GMVIO(mesh).write_equation_systems ("out_3.gmv",
                                      es);

  return 0;
}

void assemble_and_solve(MeshBase & mesh,
                        EquationSystems & equation_systems)
{
  mesh.print_info();

  LinearImplicitSystem & system =
    equation_systems.add_system<LinearImplicitSystem> ("Poisson");

  unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);

  system.attach_assemble_function (assemble_poisson);

  // the cube has boundaries IDs 0, 1, 2, 3, 4 and 5
  std::set<boundary_id_type> boundary_ids;
  for (int j = 0; j<6; ++j)
    boundary_ids.insert(j);

  // Create a vector storing the variable numbers which the BC applies to
  std::vector<unsigned int> variables(1);
  variables[0] = u_var;

  ZeroFunction<> zf;
  DirichletBoundary dirichlet_bc(boundary_ids,
                                 variables,
                                 &zf);
  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);

  equation_systems.init();
  equation_systems.print_info();

#ifdef LIBMESH_ENABLE_AMR
  MeshRefinement mesh_refinement(mesh);

  mesh_refinement.refine_fraction()  = 0.7;
  mesh_refinement.coarsen_fraction() = 0.3;
  mesh_refinement.max_h_level()      = 5;

  const unsigned int max_r_steps = 2;

  for (unsigned int r_step=0; r_step<=max_r_steps; r_step++)
    {
      system.solve();
      if (r_step != max_r_steps)
        {
          ErrorVector error;
          KellyErrorEstimator error_estimator;

          error_estimator.estimate_error(system, error);

          libMesh::out << "Error estimate\nl2 norm = "
                       << error.l2_norm()
                       << "\nmaximum = "
                       << error.maximum()
                       << std::endl;

          mesh_refinement.flag_elements_by_error_fraction (error);

          mesh_refinement.refine_and_coarsen_elements();

          equation_systems.reinit();
        }
    }
#else
  system.solve();
#endif
}