C++ OGProjection::project_global方法代码示例

int main(int argc, char* argv[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2D mloader;
  mloader.load("square.mesh", mesh);

  // Perform uniform mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh->refine_all_elements();

  // Create an L2 space with default shapeset.
  SpaceSharedPtr<double> space(new L2Space<double>(mesh, P_INIT));

  // View basis functions.
  BaseView<double> bview("BaseView", new WinGeom(0, 0, 600, 500));
  bview.show(space);
  // View::wait(H2DV_WAIT_KEYPRESS);

  // Initialize the exact and projected solution.
  MeshFunctionSharedPtr<double> sln(new Solution<double>);
  MeshFunctionSharedPtr<double> sln_exact(new CustomExactSolution(mesh));

  // Project the exact function on the FE space.
  OGProjection<double> ogProjection;
  ogProjection.project_global(space, sln_exact, sln);

  // Visualize the projection.
  ScalarView view1("Projection", new WinGeom(610, 0, 600, 500));
  view1.show(sln);

  // Wait for all views to be closed.
  View::wait();
  return 0;
}

int main(int argc, char* argv[])
{
  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("domain.mesh", &mesh);

  // Perform initial mesh refinemets.
  for (int i=0; i < INIT_REF_NUM; i++)  mesh.refine_all_elements();

  // Initialize boundary conditions.
  DefaultEssentialBCConst<double> left_t("Left", 1.0);
  EssentialBCs<double> bcs_t(&left_t);

  DefaultEssentialBCConst<double> left_c("Left", 0.0);
  EssentialBCs<double> bcs_c(&left_c);

  // Create H1 spaces with default shapesets.
  H1Space<double>* t_space = new H1Space<double>(&mesh, &bcs_t, P_INIT);
  H1Space<double>* c_space = new H1Space<double>(&mesh, &bcs_c, P_INIT);
  int ndof = Space<double>::get_num_dofs(Hermes::vector<const Space<double>*>(t_space, c_space));
  Hermes::Mixins::Loggable::Static::info("ndof = %d.", ndof);

  // Define initial conditions.
  InitialSolutionTemperature t_prev_time_1(&mesh, x1);
  InitialSolutionConcentration c_prev_time_1(&mesh, x1, Le);
  InitialSolutionTemperature t_prev_time_2(&mesh, x1);
  InitialSolutionConcentration c_prev_time_2(&mesh, x1, Le);
  Solution<double> t_prev_newton;
  Solution<double> c_prev_newton;

  // Filters for the reaction rate omega and its derivatives.
  CustomFilter omega(Hermes::vector<Solution<double>*>(&t_prev_time_1, &c_prev_time_1), Le, alpha, beta, kappa, x1, TAU);
  CustomFilterDt omega_dt(Hermes::vector<Solution<double>*>(&t_prev_time_1, &c_prev_time_1), Le, alpha, beta, kappa, x1, TAU);
  CustomFilterDc omega_dc(Hermes::vector<Solution<double>*>(&t_prev_time_1, &c_prev_time_1), Le, alpha, beta, kappa, x1, TAU);

  // Initialize visualization.
  ScalarView rview("Reaction rate", new WinGeom(0, 0, 800, 230));

  // Initialize weak formulation.
  CustomWeakForm wf(Le, alpha, beta, kappa, x1, TAU, TRILINOS_JFNK, PRECOND, &omega, &omega_dt, 
                    &omega_dc, &t_prev_time_1, &c_prev_time_1, &t_prev_time_2, &c_prev_time_2);

  // Project the functions "t_prev_time_1" and "c_prev_time_1" on the FE space 
  // in order to obtain initial vector for NOX. 
  Hermes::Mixins::Loggable::Static::info("Projecting initial solutions on the FE meshes.");
  double* coeff_vec = new double[ndof];
  OGProjection<double> ogProjection; ogProjection.project_global(Hermes::vector<const Space<double> *>(t_space, c_space), 
                                       Hermes::vector<MeshFunction<double>*>(&t_prev_time_1, &c_prev_time_1),
                                       coeff_vec);

  // Measure the projection time.
  double proj_time = cpu_time.tick().last();

  // Initialize finite element problem.
  DiscreteProblem<double> dp(&wf, Hermes::vector<const Space<double>*>(t_space, c_space));

  // Initialize NOX solver and preconditioner.
  NewtonSolverNOX<double> solver(&dp);
  MlPrecond<double> pc("sa");
  if (PRECOND)
  {
    if (TRILINOS_JFNK) 
      solver.set_precond(pc);
    else 
      solver.set_precond("New Ifpack");
  }
  if (TRILINOS_OUTPUT)
    solver.set_output_flags(NOX::Utils::Error | NOX::Utils::OuterIteration |
                            NOX::Utils::OuterIterationStatusTest |
                            NOX::Utils::LinearSolverDetails);

  // Time stepping loop:
  double total_time = 0.0;
  cpu_time.tick_reset();
  for (int ts = 1; total_time <= T_FINAL; ts++)
  {
    Hermes::Mixins::Loggable::Static::info("---- Time step %d, t = %g s", ts, total_time + TAU);

    cpu_time.tick();
    try
    {
      solver.solve(coeff_vec);
    }
    catch(std::exception& e)
    {
      std::cout << e.what();
      
    }

    Solution<double>::vector_to_solutions(solver.get_sln_vector(), Hermes::vector<const Space<double> *>(t_space, c_space), 
              Hermes::vector<Solution<double> *>(&t_prev_newton, &c_prev_newton));

    cpu_time.tick();
    Hermes::Mixins::Loggable::Static::info("Number of nonlin iterations: %d (norm of residual: %g)",
        solver.get_num_iters(), solver.get_residual());
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2D mloader;
  mloader.load("domain.mesh", mesh);

  // Initial mesh refinements.
  mesh->refine_towards_boundary(BDY_OBSTACLE, 2, false);
  mesh->refine_towards_boundary(BDY_TOP, 2, true);     // '4' is the number of levels,
  mesh->refine_towards_boundary(BDY_BOTTOM, 2, true);  // 'true' stands for anisotropic refinements.
  mesh->refine_all_elements();

  // Initialize boundary conditions.
  EssentialBCNonConst bc_left_vel_x(BDY_LEFT, VEL_INLET, H, STARTUP_TIME);
  Hermes::Hermes2D::DefaultEssentialBCConst<double> bc_other_vel_x(Hermes::vector<std::string>(BDY_BOTTOM, BDY_TOP, BDY_OBSTACLE), 0.0);
  Hermes::Hermes2D::EssentialBCs<double> bcs_vel_x(Hermes::vector<EssentialBoundaryCondition<double> *>(&bc_left_vel_x, &bc_other_vel_x));
  Hermes::Hermes2D::DefaultEssentialBCConst<double> bc_vel_y(Hermes::vector<std::string>(BDY_LEFT, BDY_BOTTOM, BDY_TOP, BDY_OBSTACLE), 0.0);
  Hermes::Hermes2D::EssentialBCs<double> bcs_vel_y(&bc_vel_y);
  Hermes::Hermes2D::EssentialBCs<double> bcs_pressure;

  // Spaces for velocity components and pressure.
  SpaceSharedPtr<double> xvel_space(new H1Space<double>(mesh, &bcs_vel_x, P_INIT_VEL));
  SpaceSharedPtr<double> yvel_space(new H1Space<double>(mesh, &bcs_vel_y, P_INIT_VEL));
#ifdef PRESSURE_IN_L2
  SpaceSharedPtr<double> p_space(new L2Space<double> (mesh, P_INIT_PRESSURE));
#else
  SpaceSharedPtr<double> p_space(new H1Space<double> (mesh, &bcs_pressure, P_INIT_PRESSURE));
#endif

  // Calculate and report the number of degrees of freedom.
  int ndof = Space<double>::get_num_dofs(Hermes::vector<SpaceSharedPtr<double> >(xvel_space, yvel_space, p_space));

  // Define projection norms.
  NormType vel_proj_norm = HERMES_H1_NORM;
#ifdef PRESSURE_IN_L2
  NormType p_proj_norm = HERMES_L2_NORM;
#else
  NormType p_proj_norm = HERMES_H1_NORM;
#endif

  // Solutions for the Newton's iteration and time stepping.
  MeshFunctionSharedPtr<double> xvel_prev_time(new ConstantSolution<double> (mesh, 0.0));
  MeshFunctionSharedPtr<double> yvel_prev_time(new ConstantSolution<double> (mesh, 0.0));
  MeshFunctionSharedPtr<double> p_prev_time(new ConstantSolution<double> (mesh, 0.0));

  // Initialize weak formulation.
  WeakFormNSNewton wf(STOKES, RE, TAU, xvel_prev_time, yvel_prev_time);
  UExtFunctionSharedPtr<double> fn_0(new CustomUExtFunction(0));
  UExtFunctionSharedPtr<double> fn_1(new CustomUExtFunction(1));
  wf.set_ext(Hermes::vector<MeshFunctionSharedPtr<double> >(xvel_prev_time, yvel_prev_time));
  wf.set_u_ext_fn(Hermes::vector<UExtFunctionSharedPtr<double> >(fn_0, fn_1));

  // Initialize the Newton solver.
  Hermes::Hermes2D::NewtonSolver<double> newton;
  newton.set_weak_formulation(&wf);
  Hermes::vector<SpaceSharedPtr<double> > spaces(xvel_space, yvel_space, p_space);
  newton.set_spaces(spaces);

  // Initialize views.
  Views::VectorView vview("velocity[m/s]", new Views::WinGeom(0, 0, 750, 240));
  Views::ScalarView pview("pressure[Pa]", new Views::WinGeom(0, 290, 750, 240));
  vview.set_min_max_range(0, 1.6);
  vview.fix_scale_width(80);
  //pview.set_min_max_range(-0.9, 1.0);
  pview.fix_scale_width(80);
  if(HERMES_VISUALIZATION)
    pview.show_mesh(true);

  // Project the initial condition on the FE space to obtain initial
  // coefficient vector for the Newton's method.
  double* coeff_vec = new double[Space<double>::get_num_dofs(Hermes::vector<SpaceSharedPtr<double> >(xvel_space, yvel_space, p_space))];
  OGProjection<double> ogProjection;

  ogProjection.project_global(Hermes::vector<SpaceSharedPtr<double> >(xvel_space, yvel_space, p_space),
    Hermes::vector<MeshFunctionSharedPtr<double> >(xvel_prev_time, yvel_prev_time, p_prev_time),
    coeff_vec, Hermes::vector<NormType>(vel_proj_norm, vel_proj_norm, p_proj_norm));

  newton.set_max_allowed_iterations(max_allowed_iterations);
  newton.set_tolerance(NEWTON_TOL, Hermes::Solvers::ResidualNormAbsolute);
  newton.set_sufficient_improvement_factor_jacobian(1e-2);
  //newton.set_jacobian_constant();

  // Time-stepping loop:
  char title[100];
  int num_time_steps = T_FINAL / TAU;
  for (int ts = 1; ts <= num_time_steps; ts++)
  {
    current_time += TAU;
    Hermes::Mixins::Loggable::Static::info("Time step %i, time %f.", ts, current_time);

    // Update time-dependent essential BCs.
    newton.set_time(current_time);

    // Perform Newton's iteration and translate the resulting coefficient vector into previous time level solutions.
    try
    {
      newton.solve(coeff_vec);
    }
    catch(Hermes::Exceptions::Exception& e)
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  Hermes::Mixins::TimeMeasurable m;
  m.tick();

  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("domain.mesh", &mesh);

  // Perform initial mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

  // Initialize boundary conditions.
  Hermes::Hermes2D::DefaultEssentialBCConst<std::complex<double> > bc_essential("Dirichlet", std::complex<double>(0.0, 0.0));
  EssentialBCs<std::complex<double> > bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  H1Space<std::complex<double> > space(&mesh, &bcs, P_INIT);
  int ndof = space.get_num_dofs();

  // Initialize the weak formulation.
  CustomWeakForm wf("Air", MU_0, "Iron", MU_IRON, GAMMA_IRON,
    "Wire", MU_0, std::complex<double>(J_EXT, 0.0), OMEGA);

  // Initialize coarse and reference mesh solution.
  Solution<std::complex<double> > sln, ref_sln;

  // Initialize refinement selector.
  H1ProjBasedSelector<std::complex<double> > selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Initialize views.
  Views::VectorView sview("Solution", new Views::WinGeom(0, 0, 600, 350));
  Views::OrderView oview("Polynomial orders", new Views::WinGeom(610, 0, 520, 350));

  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;

  DiscreteProblem<std::complex<double> > dp(&wf, &space);

  // Perform Newton's iteration and translate the resulting coefficient vector into a Solution.
  Hermes::Hermes2D::NewtonSolver<std::complex<double> > newton(&dp);
    
  Views::MeshView m1, m2;
  Views::OrderView o1, o2;
  // Adaptivity loop:
  int as = 1; bool done = false;
  do
  {
    // Construct globally refined reference mesh and setup reference space.
    Space<std::complex<double> >::ReferenceSpaceCreator ref_space_creator(&space, &mesh);
    Space<std::complex<double> >* ref_space = ref_space_creator.create_ref_space();
    
    newton.set_space(ref_space);

    int ndof_ref = ref_space->get_num_dofs();

    // Initialize reference problem.

    // Initial coefficient vector for the Newton's method.
    std::complex<double>* coeff_vec = new std::complex<double>[ndof_ref];
    memset(coeff_vec, 0, ndof_ref * sizeof(std::complex<double>));

    // Perform Newton's iteration and translate the resulting coefficient vector into a Solution.
    // For iterative solver.
    if(matrix_solver_type == SOLVER_AZTECOO)
    {
      newton.set_iterative_method(iterative_method);
      newton.set_preconditioner(preconditioner);
    }
    try
    {
      newton.solve_keep_jacobian(coeff_vec);
    }
    catch(Hermes::Exceptions::Exception& e)
    {
      e.print_msg();
    }

    Hermes::Hermes2D::Solution<std::complex<double> >::vector_to_solution(newton.get_sln_vector(), ref_space, &ref_sln);

    // Project the fine mesh solution onto the coarse mesh.
    OGProjection<std::complex<double> > ogProjection;
    ogProjection.project_global(&space, &ref_sln, &sln);

    // View the coarse mesh solution and polynomial orders.
    RealFilter real_filter(&sln);
    sview.show(&real_filter, &real_filter);

    oview.show(&space);

    // Calculate element errors and total error estimate.
    Adapt<std::complex<double> >* adaptivity = new Adapt<std::complex<double> >(&space);
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln) * 100;
    std::cout << (std::string)"Relative error: " << err_est_rel << std::endl;

    // Add entry to DOF and CPU convergence graphs.
    graph_dof.add_values(space.get_num_dofs(), err_est_rel);
    graph_dof.save("conv_dof_est.dat");

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

//.........这里部分代码省略.........
    MeshSharedPtr ref_v_mesh = refMeshCreatorV.create_ref_mesh();

    Space<double>::ReferenceSpaceCreator refSpaceCreatorV(v_space, ref_v_mesh);
    SpaceSharedPtr<double> ref_v_space = refSpaceCreatorV.create_ref_space();

    Hermes::vector<SpaceSharedPtr<double> > ref_spaces(ref_u_space, ref_v_space);

    int ndof_ref = Space<double>::get_num_dofs(ref_spaces);

    Hermes::Mixins::Loggable::Static::info("---- Adaptivity step %d (%d DOF):", as, ndof_ref);
    cpu_time.tick();

    Hermes::Mixins::Loggable::Static::info("Solving on reference mesh.");

    // Assemble the discrete problem.    
    DiscreteProblem<double> dp(&wf, ref_spaces);

    NewtonSolver<double> newton(&dp);
    
    try
    {
      newton.solve();
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.print_msg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    };

    // Translate the resulting coefficient vector into the instance of Solution.
    Solution<double>::vector_to_solutions(newton.get_sln_vector(), ref_spaces, ref_slns);

    cpu_time.tick();
    Hermes::Mixins::Loggable::Static::info("Solution: %g s", cpu_time.last());

    // Project the fine mesh solution onto the coarse mesh.
    Hermes::Mixins::Loggable::Static::info("Calculating error estimate and exact error.");
    OGProjection<double> ogProjection; ogProjection.project_global(spaces, ref_slns, slns);

    // Calculate element errors and total error estimate.
    DefaultErrorCalculator<double, HERMES_H1_NORM> error_calculator(errorType, 2);
    error_calculator.calculate_errors(slns, exact_slns);
    double err_exact_rel_total = error_calculator.get_total_error_squared() * 100.0;
    error_calculator.calculate_errors(slns, ref_slns);
    double err_est_rel_total = error_calculator.get_total_error_squared() * 100.0;

    Adapt<double> adaptivity(spaces, &error_calculator);
    adaptivity.set_strategy(&stoppingCriterion);

    cpu_time.tick();
    Hermes::Mixins::Loggable::Static::info("Error calculation: %g s", cpu_time.last());

    // Time measurement.
    cpu_time.tick();
    double accum_time = cpu_time.accumulated();

    // View the coarse mesh solution and polynomial orders.
    s_view_u.show(u_sln); 
    o_view_u.show(u_space);
    s_view_v.show(v_sln); 
    o_view_v.show(v_space);
    MeshFunctionSharedPtr<double> stress(new VonMisesFilter(Hermes::vector<MeshFunctionSharedPtr<double> >(u_sln, v_sln), lambda, mu));
    mises_view.show(stress, H2D_FN_VAL_0, u_sln, v_sln, 0.03);

    // Add entry to DOF and CPU convergence graphs.
    graph_dof_est.add_values(Space<double>::get_num_dofs(Hermes::vector<SpaceSharedPtr<double> >(u_space, v_space)), err_est_rel_total);
    graph_dof_est.save("conv_dof_est.dat");
    graph_cpu_est.add_values(cpu_time.accumulated(), err_est_rel_total);
    graph_cpu_est.save("conv_cpu_est.dat");
    graph_dof_exact.add_values(Space<double>::get_num_dofs(Hermes::vector<SpaceSharedPtr<double> >(u_space, v_space)), err_exact_rel_total);
    graph_dof_exact.save("conv_dof_exact.dat");
    graph_cpu_exact.add_values(cpu_time.accumulated(), err_exact_rel_total);
    graph_cpu_exact.save("conv_cpu_exact.dat");

    cpu_time.tick(Hermes::Mixins::TimeMeasurable::HERMES_SKIP);

    // If err_est too large, adapt the mesh.
    if (err_est_rel_total < ERR_STOP) 
      done = true;
    else 
    {
      Hermes::Mixins::Loggable::Static::info("Adapting coarse mesh.");
      done = adaptivity.adapt(Hermes::vector<RefinementSelectors::Selector<double> *>(&selector, &selector));
    }

    cpu_time.tick();
    Hermes::Mixins::Loggable::Static::info("Adaptation: %g s", cpu_time.last());

    // Increase the counter of adaptivity steps.
    if (done == false)  
      as++;
  }
  while (done == false);

  Hermes::Mixins::Loggable::Static::info("Total running time: %g s", cpu_time.accumulated());

  // Wait for all views to be closed.
  Views::View::wait();
  return 0;
}

//.........这里部分代码省略.........
    SpaceSharedPtr<double> v_ref_space = v_ref_space_creator.create_ref_space();

    Hermes::vector<SpaceSharedPtr<double> > ref_spaces_const(u_ref_space, v_ref_space);

    newton.set_spaces(ref_spaces_const);

    int ndof_ref = Space<double>::get_num_dofs(ref_spaces_const);

    // Initialize reference problem.
    Hermes::Mixins::Loggable::Static::info("Solving on reference mesh.");

    // Time measurement.
    cpu_time.tick();

    // Perform Newton's iteration.
    try
    {
      newton.solve();
    }
    catch (Hermes::Exceptions::Exception& e)
    {
      std::cout << e.info();
    }
    catch (std::exception& e)
    {
      std::cout << e.what();
    }

    // Translate the resulting coefficient vector into the instance of Solution.
    Solution<double>::vector_to_solutions(newton.get_sln_vector(), ref_spaces_const, Hermes::vector<MeshFunctionSharedPtr<double> >(u_ref_sln, v_ref_sln));

    // Project the fine mesh solution onto the coarse mesh.
    Hermes::Mixins::Loggable::Static::info("Projecting reference solution on coarse mesh.");
    OGProjection<double> ogProjection; ogProjection.project_global(Hermes::vector<SpaceSharedPtr<double> >(u_space, v_space), ref_slns, slns);

    cpu_time.tick();

    // View the coarse mesh solution and polynomial orders.
    s_view_0.show(u_sln);
    o_view_0.show(u_space);
    s_view_1.show(v_sln);
    o_view_1.show(v_space);

    // Calculate element errors.
    Hermes::Mixins::Loggable::Static::info("Calculating error estimate and exact error.");
    errorCalculator.calculate_errors(slns, exact_slns, false);
    double err_exact_rel_total = errorCalculator.get_total_error_squared() * 100;
    Hermes::vector<double> err_exact_rel;
    err_exact_rel.push_back(errorCalculator.get_error_squared(0) * 100);
    err_exact_rel.push_back(errorCalculator.get_error_squared(1) * 100);

    errorCalculator.calculate_errors(slns, ref_slns, true);
    double err_est_rel_total = errorCalculator.get_total_error_squared() * 100;
    Hermes::vector<double> err_est_rel;
    err_est_rel.push_back(errorCalculator.get_error_squared(0) * 100);
    err_est_rel.push_back(errorCalculator.get_error_squared(1) * 100);

    adaptivity.set_spaces(Hermes::vector<SpaceSharedPtr<double> >(u_space, v_space));

    // Time measurement.
    cpu_time.tick();

    // Report results.
    Hermes::Mixins::Loggable::Static::info("ndof_coarse[0]: %d, ndof_fine[0]: %d",
      u_space->get_num_dofs(), u_ref_space->get_num_dofs());
    Hermes::Mixins::Loggable::Static::info("err_est_rel[0]: %g%%, err_exact_rel[0]: %g%%", err_est_rel[0], err_exact_rel[0]);

//.........这里部分代码省略.........
      case 2: space->unrefine_all_mesh_elements();
        space->set_uniform_order(P_INIT);
        break;
      case 3: space->unrefine_all_mesh_elements();
        //space->adjust_element_order(-1, P_INIT);
        space->adjust_element_order(-1, -1, P_INIT, P_INIT);
        break;
      default: throw Hermes::Exceptions::Exception("Wrong global derefinement method.");
      }

      space->assign_dofs();
      ndof = Space<double>::get_num_dofs(space);
    }

    // Spatial adaptivity loop. Note: sln_prev_time must not be 
    // changed during spatial adaptivity. 
    MeshFunctionSharedPtr<double> ref_sln(new Solution<double>());
    MeshFunctionSharedPtr<double> time_error_fn(new Solution<double>(mesh));
    bool done = false; int as = 1;
    double err_est;
    do {
      // Construct globally refined reference mesh and setup reference space.
      Mesh::ReferenceMeshCreator refMeshCreator(mesh);
      MeshSharedPtr ref_mesh = refMeshCreator.create_ref_mesh();

      Space<double>::ReferenceSpaceCreator refSpaceCreator(space, ref_mesh);
      SpaceSharedPtr<double> ref_space = refSpaceCreator.create_ref_space();

      // Initialize Runge-Kutta time stepping on the reference mesh.
      RungeKutta<double> runge_kutta(&wf, ref_space, &bt);

      try
      {
        ogProjection.project_global(ref_space, sln_prev_time, 
          sln_prev_time);
      }
      catch(Exceptions::Exception& e)
      {
        std::cout << e.what() << std::endl;
        Hermes::Mixins::Loggable::Static::error("Projection failed.");

        return -1;
      }

      // Runge-Kutta step on the fine mesh->
      Hermes::Mixins::Loggable::Static::info("Runge-Kutta time step on fine mesh (t = %g s, tau = %g s, stages: %d).", 
        current_time, time_step, bt.get_size());
      bool verbose = true;
      bool jacobian_changed = false;

      try
      {
        runge_kutta.set_time(current_time);
        runge_kutta.set_time_step(time_step);
        runge_kutta.set_max_allowed_iterations(NEWTON_MAX_ITER);
        runge_kutta.set_tolerance(NEWTON_TOL_FINE);
        runge_kutta.rk_time_step_newton(sln_prev_time, ref_sln, bt.is_embedded() ? time_error_fn : NULL);
      }
      catch(Exceptions::Exception& e)
      {
        std::cout << e.what() << std::endl;
        Hermes::Mixins::Loggable::Static::error("Runge-Kutta time step failed");

        return -1;
      }

int main(int argc, char* argv[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2D mloader;
  mloader.load("domain.mesh", mesh);

  // Perform initial mesh refinemets.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh->refine_all_elements();

  // Initialize solutions.
  MeshFunctionSharedPtr<double> E_sln(new CustomInitialConditionWave(mesh));
  MeshFunctionSharedPtr<double>  F_sln(new ZeroSolutionVector<double>(mesh));
  Hermes::vector<MeshFunctionSharedPtr<double> > slns(E_sln, F_sln);

  // Initialize the weak formulation.
  CustomWeakFormWaveIE wf(time_step, C_SQUARED, E_sln, F_sln);

  // Initialize boundary conditions
  DefaultEssentialBCConst<double> bc_essential("Perfect conductor", 0.0);
  EssentialBCs<double> bcs(&bc_essential);

  SpaceSharedPtr<double> E_space(new HcurlSpace<double>(mesh, &bcs, P_INIT));
  SpaceSharedPtr<double> F_space(new HcurlSpace<double>(mesh, &bcs, P_INIT));

  Hermes::vector<SpaceSharedPtr<double> > spaces = Hermes::vector<SpaceSharedPtr<double> >(E_space, F_space);
  int ndof = HcurlSpace<double>::get_num_dofs(spaces);
  Hermes::Mixins::Loggable::Static::info("ndof = %d.", ndof);

  // Initialize the FE problem.
  DiscreteProblem<double> dp(&wf, spaces);

  // Project the initial condition on the FE space to obtain initial 
  // coefficient vector for the Newton's method.
  // NOTE: If you want to start from the zero vector, just define 
  // coeff_vec to be a vector of ndof zeros (no projection is needed).
  Hermes::Mixins::Loggable::Static::info("Projecting to obtain initial vector for the Newton's method.");
  double* coeff_vec = new double[ndof];
  OGProjection<double> ogProjection; ogProjection.project_global(spaces, slns, coeff_vec); 

  // Initialize Newton solver.
  NewtonSolver<double> newton(&dp);

  // Initialize views.
  ScalarView E1_view("Solution E1", new WinGeom(0, 0, 400, 350));
  E1_view.fix_scale_width(50);
  ScalarView E2_view("Solution E2", new WinGeom(410, 0, 400, 350));
  E2_view.fix_scale_width(50);
  ScalarView F1_view("Solution F1", new WinGeom(0, 410, 400, 350));
  F1_view.fix_scale_width(50);
  ScalarView F2_view("Solution F2", new WinGeom(410, 410, 400, 350));
  F2_view.fix_scale_width(50);

  // Time stepping loop.
  double current_time = 0; int ts = 1;
  do
  {
    // Perform one implicit Euler time step.
    Hermes::Mixins::Loggable::Static::info("Implicit Euler time step (t = %g s, time_step = %g s).", current_time, time_step);

    // Perform Newton's iteration.
    try
    {
      newton.set_max_allowed_iterations(NEWTON_MAX_ITER);
      newton.set_tolerance(NEWTON_TOL, Hermes::Solvers::ResidualNormAbsolute);
      newton.solve(coeff_vec);
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.print_msg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    }

    // Translate the resulting coefficient vector into Solutions.
    Solution<double>::vector_to_solutions(newton.get_sln_vector(), spaces, slns);

    // Visualize the solutions.
    char title[100];
    sprintf(title, "E1, t = %g", current_time + time_step);
    E1_view.set_title(title);
    E1_view.show(E_sln, H2D_FN_VAL_0);
    sprintf(title, "E2, t = %g", current_time + time_step);
    E2_view.set_title(title);
    E2_view.show(E_sln, H2D_FN_VAL_1);

    sprintf(title, "F1, t = %g", current_time + time_step);
    F1_view.set_title(title);
    F1_view.show(F_sln, H2D_FN_VAL_0);
    sprintf(title, "F2, t = %g", current_time + time_step);
    F2_view.set_title(title);
    F2_view.show(F_sln, H2D_FN_VAL_1);

    //View::wait();

    // Update time.
    current_time += time_step;

  } while (current_time < T_FINAL);

  // Wait for the view to be closed.
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  if (ALIGN_MESH) 
    mloader.load("oven_load_circle.mesh", &mesh);
  else 
    mloader.load("oven_load_square.mesh", &mesh);

  // Perform initial mesh refinemets.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

  // Initialize boundary conditions
  DefaultEssentialBCConst<std::complex<double> > bc_essential(BDY_PERFECT_CONDUCTOR, std::complex<double>(0.0, 0.0));

  EssentialBCs<std::complex<double> > bcs(&bc_essential);

  // Create an Hcurl space with default shapeset.
  HcurlSpace<std::complex<double> > space(&mesh, &bcs, P_INIT);
  int ndof = space.get_num_dofs();
  Hermes::Mixins::Loggable::Static::info("ndof = %d", ndof);

  // Initialize the weak formulation.
  CustomWeakForm wf(e_0, mu_0, mu_r, kappa, omega, J, ALIGN_MESH, &mesh, BDY_CURRENT);

  // Initialize coarse and reference mesh solution.
  Solution<std::complex<double> > sln, ref_sln;

  // Initialize refinements selector.
  HcurlProjBasedSelector<std::complex<double> > selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Initialize views.
  ScalarView eview("Electric field", new WinGeom(0, 0, 580, 400));
  OrderView  oview("Polynomial orders", new WinGeom(590, 0, 550, 400));
  
  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;
  
  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

  // Adaptivity loop:
  int as = 1; bool done = false;
  do
  {
    Hermes::Mixins::Loggable::Static::info("---- Adaptivity step %d:", as);

    // Construct globally refined reference mesh and setup reference space.
    Mesh::ReferenceMeshCreator refMeshCreator(&mesh);
    Mesh* ref_mesh = refMeshCreator.create_ref_mesh();

    Space<std::complex<double> >::ReferenceSpaceCreator refSpaceCreator(&space, ref_mesh);
    Space<std::complex<double> >* ref_space = refSpaceCreator.create_ref_space();
    int ndof_ref = Space<std::complex<double> >::get_num_dofs(ref_space);

    // Initialize reference problem.
    Hermes::Mixins::Loggable::Static::info("Solving on reference mesh.");
    DiscreteProblem<std::complex<double> > dp(&wf, ref_space);

    // Time measurement.
    cpu_time.tick();

    // Perform Newton's iteration.
    Hermes::Hermes2D::NewtonSolver<std::complex<double> > newton(&dp);
    try
    {
      newton.set_newton_max_iter(NEWTON_MAX_ITER);
      newton.set_newton_tol(NEWTON_TOL);
      newton.solve();
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.print_msg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    };
    // Translate the resulting coefficient vector into the Solution<std::complex<double> > sln.
    Hermes::Hermes2D::Solution<std::complex<double> >::vector_to_solution(newton.get_sln_vector(), ref_space, &ref_sln);
  
    // Project the fine mesh solution onto the coarse mesh.
    Hermes::Mixins::Loggable::Static::info("Projecting reference solution on coarse mesh.");
    OGProjection<std::complex<double> > ogProjection; ogProjection.project_global(&space, &ref_sln, &sln); 
   
    // View the coarse mesh solution and polynomial orders.
    RealFilter real(&sln);
    MagFilter<double> magn(&real);
    ValFilter limited_magn(&magn, 0.0, 4e3);
    char title[100];
    sprintf(title, "Electric field, adaptivity step %d", as);
    eview.set_title(title);
    //eview.set_min_max_range(0.0, 4e3);
    eview.show(&limited_magn);
    sprintf(title, "Polynomial orders, adaptivity step %d", as);
    oview.set_title(title);
    oview.show(&space);

    // Calculate element errors and total error estimate.
    Hermes::Mixins::Loggable::Static::info("Calculating error estimate."); 
    Adapt<std::complex<double> >* adaptivity = new Adapt<std::complex<double> >(&space);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("domain.mesh", &mesh);

  // Initialize the weak formulation.
  CustomWeakFormPoisson wf("Motor", EPS_MOTOR, "Air", EPS_AIR);
  
  // Initialize boundary conditions
  DefaultEssentialBCConst<double> bc_essential_out("Outer", 0.0);
  DefaultEssentialBCConst<double> bc_essential_stator("Stator", VOLTAGE);
  EssentialBCs<double> bcs(Hermes::vector<EssentialBoundaryCondition<double> *>(&bc_essential_out, &bc_essential_stator));

  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, &bcs, P_INIT);

  // Initialize coarse and fine mesh solution.
  Solution<double> sln, ref_sln;

  // Initialize refinement selector.
  H1ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Initialize views.
  Views::ScalarView sview("Solution", new Views::WinGeom(0, 0, 410, 600));
  sview.fix_scale_width(50);
  sview.show_mesh(false);
  Views::OrderView  oview("Polynomial orders", new Views::WinGeom(420, 0, 400, 600));

  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;

  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;

  DiscreteProblem<double> dp(&wf, &space);
  NewtonSolver<double> newton(&dp);
  newton.set_verbose_output(true);

  // Adaptivity loop:
  int as = 1; bool done = false;
  do
  {
    Hermes::Mixins::Loggable::Static::info("---- Adaptivity step %d:", as);
    
    // Time measurement.
    cpu_time.tick();

    // Construct globally refined mesh and setup fine mesh space.
    Mesh::ReferenceMeshCreator ref_mesh_creator(&mesh);
    Mesh* ref_mesh = ref_mesh_creator.create_ref_mesh();
    Space<double>::ReferenceSpaceCreator ref_space_creator(&space, ref_mesh);
    Space<double>* ref_space = ref_space_creator.create_ref_space();
    int ndof_ref = ref_space->get_num_dofs();

    // Initialize fine mesh problem.
    Hermes::Mixins::Loggable::Static::info("Solving on fine mesh.");
    
    newton.set_space(ref_space);


    // Perform Newton's iteration.
    try
    {
      newton.solve();
    }
    catch(std::exception& e)
    {
      std::cout << e.what();
      
    }

    // Translate the resulting coefficient vector into the instance of Solution.
    Solution<double>::vector_to_solution(newton.get_sln_vector(), ref_space, &ref_sln);
    
    // Project the fine mesh solution onto the coarse mesh.
    Hermes::Mixins::Loggable::Static::info("Projecting fine mesh solution on coarse mesh.");
    OGProjection<double> ogProjection; ogProjection.project_global(&space, &ref_sln, &sln);

    // Time measurement.
    cpu_time.tick();

    // VTK output.
    if (VTK_VISUALIZATION) 
    {
      // Output solution in VTK format.
      Views::Linearizer lin;
      char* title = new char[100];
      sprintf(title, "sln-%d.vtk", as);
      lin.save_solution_vtk(&sln, title, "Potential", false);
      Hermes::Mixins::Loggable::Static::info("Solution in VTK format saved to file %s.", title);

      // Output mesh and element orders in VTK format.
      Views::Orderizer ord;
      sprintf(title, "ord-%d.vtk", as);
      ord.save_orders_vtk(&space, title);
      Hermes::Mixins::Loggable::Static::info("Element orders in VTK format saved to file %s.", title);
    }

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

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("domain.mesh", &mesh);
  
  // Perform initial uniform mesh refinement.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

  // Set essential boundary conditions.
  DefaultEssentialBCConst<double> bc_essential(Hermes::vector<std::string>("right", "top"), 0.0);
  EssentialBCs<double> bcs(&bc_essential);
  
  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, &bcs, P_INIT);

  // Associate element markers (corresponding to physical regions) 
  // with material properties (diffusion coefficient, absorption 
  // cross-section, external sources).
  Hermes::vector<std::string> regions("e1", "e2", "e3", "e4", "e5");
  Hermes::vector<double> D_map(D_1, D_2, D_3, D_4, D_5);
  Hermes::vector<double> Sigma_a_map(SIGMA_A_1, SIGMA_A_2, SIGMA_A_3, SIGMA_A_4, SIGMA_A_5);
  Hermes::vector<double> Sources_map(Q_EXT_1, 0.0, Q_EXT_3, 0.0, 0.0);
  
  // Initialize the weak formulation.
  WeakFormsNeutronics::Monoenergetic::Diffusion::DefaultWeakFormFixedSource<double>
    wf(regions, D_map, Sigma_a_map, Sources_map);

  // Initialize coarse and reference mesh solution.
  Solution<double> sln, ref_sln;
  
  // Initialize refinement selector.
  H1ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Initialize views.
  ScalarView sview("Solution", new WinGeom(0, 0, 440, 350));
  sview.fix_scale_width(50);
  sview.show_mesh(false);
  OrderView  oview("Polynomial orders", new WinGeom(450, 0, 400, 350));
  
  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;
  
  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

  // Adaptivity loop:
  int as = 1; bool done = false;
  do
  {
    Hermes::Mixins::Loggable::Static::info("---- Adaptivity step %d:", as);
    
    // Time measurement.
    cpu_time.tick();

    // Construct globally refined mesh and setup fine mesh space.
    Space<double>* ref_space = Space<double>::construct_refined_space(&space);
    int ndof_ref = ref_space->get_num_dofs();

    // Initialize fine mesh problem.
    Hermes::Mixins::Loggable::Static::info("Solving on fine mesh.");
    DiscreteProblem<double> dp(&wf, ref_space);
    
    NewtonSolver<double> newton(&dp);
    newton.set_verbose_output(false);

    // Perform Newton's iteration.
    try
    {
      newton.solve();
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.printMsg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    }

    // Translate the resulting coefficient vector into the instance of Solution.
    Solution<double>::vector_to_solution(newton.get_sln_vector(), ref_space, &ref_sln);
    
    // Project the fine mesh solution onto the coarse mesh.
    Hermes::Mixins::Loggable::Static::info("Projecting fine mesh solution on coarse mesh.");
    OGProjection<double> ogProjection; ogProjection.project_global(&space, &ref_sln, &sln);

    // Time measurement.
    cpu_time.tick();

    // Visualize the solution and mesh.
    sview.show(&sln);
    oview.show(&space);

    // Skip visualization time.
    cpu_time.tick(Hermes::Mixins::TimeMeasurable::HERMES_SKIP);

    // Calculate element errors and total error estimate.
    Hermes::Mixins::Loggable::Static::info("Calculating error estimate.");
    Adapt<double> adaptivity(&space);
    bool solutions_for_adapt = true;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      // Construct globally refined reference mesh and setup reference space.
      // FIXME: This should be increase in the x-direction only.
      int order_increase = 1;          
      // FIXME: This should be '2' but that leads to a segfault.
      int refinement_type = 0;         
      Mesh::ReferenceMeshCreator refMeshCreator(&mesh);
      Mesh* ref_mesh = refMeshCreator.create_ref_mesh();

      Space<double>::ReferenceSpaceCreator refSpaceCreator(&space, ref_mesh);
      Space<double>* ref_space = refSpaceCreator.create_ref_space();
      int ndof_ref = ref_space->get_num_dofs();

      // Initialize Runge-Kutta time stepping.
      RungeKutta<double> runge_kutta(&wf, ref_space, &bt);

      // Perform one Runge-Kutta time step according to the selected Butcher's table.
      Hermes::Mixins::Loggable::Static::info("Runge-Kutta time step (t = %g s, tau = %g s, stages: %d).",
          current_time, time_step, bt.get_size());
      bool freeze_jacobian = true;
      bool block_diagonal_jacobian = false;
      bool verbose = true;
      
      try
      {
        runge_kutta.set_time(current_time);
        runge_kutta.set_time_step(time_step);
        runge_kutta.set_newton_max_iter(NEWTON_MAX_ITER);
        runge_kutta.set_newton_tol(NEWTON_TOL);
        runge_kutta.rk_time_step_newton(&sln_time_prev, &sln_time_new);
      }
      catch(Exceptions::Exception& e)
      {
        e.print_msg();
        throw Hermes::Exceptions::Exception("Runge-Kutta time step failed");
      }

      // Project the fine mesh solution onto the coarse mesh.
      Solution<double> sln_coarse;
      Hermes::Mixins::Loggable::Static::info("Projecting fine mesh solution on coarse mesh for error estimation.");
      OGProjection<double> ogProjection; ogProjection.project_global(&space, &sln_time_new, &sln_coarse); 

      // Calculate element errors and total error estimate.
      Hermes::Mixins::Loggable::Static::info("Calculating error estimate.");
      Adapt<double>* adaptivity = new Adapt<double>(&space);
      double err_est_rel_total = adaptivity->calc_err_est(&sln_coarse, &sln_time_new) * 100;

      // Report results.
      Hermes::Mixins::Loggable::Static::info("ndof_coarse: %d, ndof_ref: %d, err_est_rel: %g%%", 
           space.get_num_dofs(), ref_space->get_num_dofs(), err_est_rel_total);

      // If err_est too large, adapt the mesh.
      if (err_est_rel_total < ERR_STOP) done = true;
      else 
      {
        Hermes::Mixins::Loggable::Static::info("Adapting the coarse mesh.");
        done = adaptivity->adapt(&selector, THRESHOLD, STRATEGY, MESH_REGULARITY);

        if (space.get_num_dofs() >= NDOF_STOP) 
          done = true;
        else
          // Increase the counter of performed adaptivity steps.
          as++;
      }
 
      // Clean up.
      delete adaptivity;
      if(!done)
      {
        delete sln_time_new.get_mesh();
        delete ref_space;
      }
    }
    while (done == false);

    // Visualize the solution and mesh.
    char title[100];
    sprintf(title, "Solution, time %g", current_time);
    sview.set_title(title);
    sview.show_mesh(false);
    sview.show(&sln_time_new);
    sprintf(title, "Mesh, time %g", current_time);
    oview.set_title(title);
    oview.show(&space);

    // Copy last reference solution into sln_time_prev.
    sln_time_prev.copy(&sln_time_new);

    dof_history_graph.add_values(current_time, space.get_num_dofs());
    dof_history_graph.save("dof_history.dat");

    // Increase current time and counter of time steps.
    current_time += time_step;
    ts++;
  }
  while (current_time < T_FINAL);

  // Wait for all views to be closed.
  Views::View::wait();
  return 0;
}

int main(int argc, char* argv[])
{
    // Load the mesh.
    MeshSharedPtr mesh(new Mesh);
    MeshReaderH2D mloader;
    mloader.load("domain.mesh", mesh);

    // Perform initial mesh refinements.
    for (int i = 0; i < INIT_REF_NUM; i++)
        mesh->refine_all_elements();

    // Initialize boundary conditions.
    Hermes::Hermes2D::DefaultEssentialBCConst<complex> bc_essential("Dirichlet", complex(0.0, 0.0));
    EssentialBCs<complex> bcs(&bc_essential);

    // Create an H1 space with default shapeset.
    SpaceSharedPtr<complex> space(new H1Space<complex>(mesh, &bcs, P_INIT));

    // Initialize the weak formulation.
    CustomWeakForm wf("Air", MU_0, "Iron", MU_IRON, GAMMA_IRON,
                      "Wire", MU_0, complex(J_EXT, 0.0), OMEGA);

    // Initialize coarse and reference mesh solution.
    MeshFunctionSharedPtr<complex> sln(new Hermes::Hermes2D::Solution<complex>());
    MeshFunctionSharedPtr<complex> ref_sln(new Hermes::Hermes2D::Solution<complex>());

    // Initialize refinement selector.
    H1ProjBasedSelector<complex> selector(CAND_LIST);

    // DOF and CPU convergence graphs initialization.
    SimpleGraph graph_dof, graph_cpu;

    DiscreteProblem<complex> dp(&wf, space);

    // Perform Newton's iteration and translate the resulting coefficient vector into a Solution.
    Hermes::Hermes2D::NewtonSolver<complex> newton(&dp);

    // Adaptivity loop:
    int as = 1;
    bool done = false;
    adaptivity.set_space(space);
    do
    {
        // Construct globally refined reference mesh and setup reference space->
        Mesh::ReferenceMeshCreator ref_mesh_creator(mesh);
        MeshSharedPtr ref_mesh = ref_mesh_creator.create_ref_mesh();
        Space<complex>::ReferenceSpaceCreator ref_space_creator(space, ref_mesh);
        SpaceSharedPtr<complex> ref_space = ref_space_creator.create_ref_space();

        newton.set_space(ref_space);

        int ndof_ref = ref_space->get_num_dofs();

        // Initialize reference problem.

        // Initial coefficient vector for the Newton's method.
        complex* coeff_vec = new complex[ndof_ref];
        memset(coeff_vec, 0, ndof_ref * sizeof(complex));

        // Perform Newton's iteration and translate the resulting coefficient vector into a Solution.
        try
        {
            newton.solve(coeff_vec);
        }
        catch(Hermes::Exceptions::Exception& e)
        {
            e.print_msg();
        }

        Hermes::Hermes2D::Solution<complex>::vector_to_solution(newton.get_sln_vector(), ref_space, ref_sln);

        // Project the fine mesh solution onto the coarse mesh.
        OGProjection<complex> ogProjection;
        ogProjection.project_global(space, ref_sln, sln);

        // Calculate element errors and total error estimate.
        errorCalculator.calculate_errors(sln, ref_sln);

        // If err_est too large, adapt the mesh->
        if(errorCalculator.get_total_error_squared()  * 100. < ERR_STOP)
            done = true;
        else
        {
            adaptivity.adapt(&selector);
        }

        // Clean up.
        delete [] coeff_vec;

        // Increase counter.
        as++;
    }
    while (done == false);

    complex sum = 0;
    for (int i = 0; i < space->get_num_dofs(); i++)
        sum += newton.get_sln_vector()[i];
    printf("coefficient sum = %f\n", sum);

    complex expected_sum;
//.........这里部分代码省略.........

int main(int argc, char* args[])
{
  // Load the mesh->
  HermesSharedPtr mesh;
  MeshReaderH2D mloader;
  mloader.load("../square.mesh", mesh);

  // Perform initial mesh refinement.
  for (int i=0; i<INIT_REF; i++) 
    mesh->refine_all_elements();

  // Create an L2 space->
  L2Space<double> space(mesh, P_INIT);

  // Initialize refinement selector.
  L2ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Disable weighting of refinement candidates.
  selector.set_error_weights(1, 1, 1);

  Solution<double> sln;
  Solution<double> ref_sln;

  // Initialize the weak formulation.
  CustomWeakForm wf("Bdy_bottom_left", mesh);

  // Initialize the FE problem.
  DiscreteProblemLinear<double> dp(&wf, space);

  // Initialize linear solver.
  Hermes::Hermes2D::LinearSolver<double> linear_solver(&dp);

  int as = 1; bool done = false;
  do
  {
    // Construct globally refined reference mesh
    // and setup reference space->
    Mesh::ReferenceMeshCreator ref_mesh_creator(mesh);
    Mesh* ref_mesh = ref_mesh_creator.create_ref_mesh();
    Space<double>::ReferenceSpaceCreator ref_space_creator(space, ref_mesh);
    Space<double>* ref_space = ref_space_creator.create_ref_space();

    dp.set_space(ref_space);

    // Solve the linear system. If successful, obtain the solution.
    try
    {
      linear_solver.solve();
      Solution<double>::vector_to_solution(linear_solver.get_sln_vector(), ref_space, &ref_sln);
    }
    catch(std::exception& e)
    {
      std::cout << e.what();
    }
    // Project the fine mesh solution onto the coarse mesh->
    OGProjection<double> ogProjection;
    ogProjection.project_global(space, &ref_sln, &sln, HERMES_L2_NORM);

    // Calculate element errors and total error estimate.
    Adapt<double>* adaptivity = new Adapt<double>(space);
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln) * 100;

    // If err_est_rel too large, adapt the mesh->
    if(err_est_rel < ERR_STOP) done = true;
    else
    {
      done = adaptivity->adapt(&selector, THRESHOLD, STRATEGY, MESH_REGULARITY);

      if(Space<double>::get_num_dofs(space) >= NDOF_STOP)
      {
        done = true;
        break;
      }
    }

    // Clean up.
    delete adaptivity;
    if(done == false)
      delete ref_space->get_mesh();
    delete ref_space;

    as++;
  }
  while (done == false);

  if(done)
  {
    printf("Success!\n");
    return 0;
  }
  else
  {
    printf("Failure!\n");
    return -1;
  }
}

//.........这里部分代码省略.........
      ndof_coarse = Space<double>::get_num_dofs(space);
    }

    // Spatial adaptivity loop. Note: h_time_prev must not be changed 
    // during spatial adaptivity. 
    bool done = false; int as = 1;
    double err_est;
    do {
      Hermes::Mixins::Loggable::Static::info("Time step %d, adaptivity step %d:", ts, as);

      // Construct globally refined reference mesh and setup reference space.
      Mesh::ReferenceMeshCreator refMeshCreator(mesh);
      MeshSharedPtr ref_mesh = refMeshCreator.create_ref_mesh();

      Space<double>::ReferenceSpaceCreator refSpaceCreator(space, ref_mesh);
      SpaceSharedPtr<double> ref_space = refSpaceCreator.create_ref_space();
      int ndof_ref = Space<double>::get_num_dofs(ref_space);

      // Time measurement.
      cpu_time.tick();

      // Initialize Runge-Kutta time stepping.
      RungeKutta<double> runge_kutta(&wf, ref_space, &bt);

      // Perform one Runge-Kutta time step according to the selected Butcher's table.
      Hermes::Mixins::Loggable::Static::info("Runge-Kutta time step (t = %g s, tau = %g s, stages: %d).",
           current_time, time_step, bt.get_size());
      try
      {
        runge_kutta.set_time(current_time);
        runge_kutta.set_time_step(time_step);
        runge_kutta.set_max_allowed_iterations(NEWTON_MAX_ITER);
        runge_kutta.set_tolerance(NEWTON_TOL);
        runge_kutta.rk_time_step_newton(h_time_prev, h_time_new);
      }
      catch(Exceptions::Exception& e)
      {
        e.print_msg();
        throw Hermes::Exceptions::Exception("Runge-Kutta time step failed");
      }

      // Project the fine mesh solution onto the coarse mesh.
      MeshFunctionSharedPtr<double> sln_coarse(new Solution<double>);
      Hermes::Mixins::Loggable::Static::info("Projecting fine mesh solution on coarse mesh for error estimation.");
      OGProjection<double> ogProjection; ogProjection.project_global(space, h_time_new, sln_coarse); 

      // Calculate element errors and total error estimate.
      Hermes::Mixins::Loggable::Static::info("Calculating error estimate.");
      errorCalculator.calculate_errors(sln_coarse, h_time_new, true);
      double err_est_rel_total = errorCalculator.get_total_error_squared() * 100;

      // Report results.
      Hermes::Mixins::Loggable::Static::info("ndof_coarse: %d, ndof_ref: %d, err_est_rel: %g%%", 
           Space<double>::get_num_dofs(space), Space<double>::get_num_dofs(ref_space), err_est_rel_total);

      // Time measurement.
      cpu_time.tick();

      // If err_est too large, adapt the mesh.
      if (err_est_rel_total < ERR_STOP) done = true;
      else 
      {
        Hermes::Mixins::Loggable::Static::info("Adapting the coarse mesh.");
        done = adaptivity.adapt(&selector);

        // Increase the counter of performed adaptivity steps.
        as++;
      }
    }
    while (done == false);

    // Add entry to DOF and CPU convergence graphs.
    graph_dof.add_values(current_time, Space<double>::get_num_dofs(space));
    graph_dof.save("conv_dof_est.dat");
    graph_cpu.add_values(current_time, cpu_time.accumulated());
    graph_cpu.save("conv_cpu_est.dat");

    // Visualize the solution and mesh->
    char title[100];
    sprintf(title, "Solution, time %g", current_time);
    view.set_title(title);
    view.show_mesh(false);
    view.show(h_time_new);
    sprintf(title, "Mesh, time %g", current_time);
    ordview.set_title(title);
    ordview.show(space);

    // Copy last reference solution into h_time_prev.
    h_time_prev->copy(h_time_new);

    // Increase current time and counter of time steps.
    current_time += time_step;
    ts++;
  }
  while (current_time < T_FINAL);

  // Wait for all views to be closed.
  View::wait();
  return 0;
}