C++ MeshReaderH2D类代码示例

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable bt(butcher_table_type);
  if (bt.is_explicit()) Hermes::Mixins::Loggable::Static::info("Using a %d-stage explicit R-K method.", bt.get_size());
  if (bt.is_diagonally_implicit()) Hermes::Mixins::Loggable::Static::info("Using a %d-stage diagonally implicit R-K method.", bt.get_size());
  if (bt.is_fully_implicit()) Hermes::Mixins::Loggable::Static::info("Using a %d-stage fully implicit R-K method.", bt.get_size());

  // Turn off adaptive time stepping if R-K method is not embedded.
  if (bt.is_embedded() == false && ADAPTIVE_TIME_STEP_ON == true) {
    throw Hermes::Exceptions::Exception("R-K method not embedded, turning off adaptive time stepping.");
    ADAPTIVE_TIME_STEP_ON = false;
  }

  // Load the mesh.
  Mesh mesh, basemesh;
  MeshReaderH2D mloader;
  mloader.load("square.mesh", &basemesh);
  mesh.copy(&basemesh);

  // Initial mesh refinements.
  for(int i = 0; i < INIT_GLOB_REF_NUM; i++) mesh.refine_all_elements();
  mesh.refine_towards_boundary("Bdy", INIT_BDY_REF_NUM);

  // Initialize boundary conditions.
  EssentialBCNonConst bc_essential("Bdy");
  EssentialBCs<double> bcs(&bc_essential);

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

  // Convert initial condition into a Solution.
  CustomInitialCondition sln_time_prev(&mesh);

  // Initialize the weak formulation
  CustomNonlinearity lambda(alpha);
  Hermes2DFunction<double> f(heat_src);
  WeakFormsH1::DefaultWeakFormPoisson<double> wf(HERMES_ANY, &lambda, &f);

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

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

  // Visualize initial condition.
  char title[100];
  ScalarView sln_view("Initial condition", new WinGeom(0, 0, 440, 350));
  sln_view.show_mesh(false);
  OrderView ordview("Initial mesh", new WinGeom(445, 0, 440, 350));
  ScalarView time_error_view("Temporal error", new WinGeom(0, 400, 440, 350));
  time_error_view.fix_scale_width(60);
  ScalarView space_error_view("Spatial error", new WinGeom(445, 400, 440, 350));
  space_error_view.fix_scale_width(60);
  sln_view.show(&sln_time_prev);
  ordview.show(&space);

  // Graph for time step history.
  SimpleGraph time_step_graph;
  if (ADAPTIVE_TIME_STEP_ON) Hermes::Mixins::Loggable::Static::info("Time step history will be saved to file time_step_history.dat.");
  
  // Time stepping loop.
  double current_time = 0.0; int ts = 1;
  do 
  {
    Hermes::Mixins::Loggable::Static::info("Begin time step %d.", ts);
    // Periodic global derefinement.
    if (ts > 1 && ts % UNREF_FREQ == 0) 
    {
      Hermes::Mixins::Loggable::Static::info("Global mesh derefinement.");
      switch (UNREF_METHOD) {
        case 1: mesh.copy(&basemesh);
                space.set_uniform_order(P_INIT);
                break;
        case 2: mesh.unrefine_all_elements();
                space.set_uniform_order(P_INIT);
                break;
        case 3: mesh.unrefine_all_elements();
                space.adjust_element_order(-1, -1, P_INIT, P_INIT);
                break;
      }

      ndof = Space<double>::get_num_dofs(&space);
    }
    Hermes::Mixins::Loggable::Static::info("ndof: %d", ndof);

    // Spatial adaptivity loop. Note: sln_time_prev must not be 
    // changed during spatial adaptivity. 
    Solution<double> ref_sln;
    Solution<double>* time_error_fn;
    if (bt.is_embedded() == true) time_error_fn = new Solution<double>(&mesh);
    else time_error_fn = NULL;
    bool done = false; int as = 1;
    double err_est;
    do {
      // Construct globally refined reference mesh and setup reference space.
      Space<double>* ref_space = Space<double>::construct_refined_space(&space);

      RungeKutta<double> runge_kutta(&wf, ref_space, &bt);
//.........这里部分代码省略.........

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

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

  // Set exact solution.
  CustomExactSolution exact_sln(&mesh, alpha);

  // Define right-hand side.
  CustomRightHandSide f(alpha);

  // Initialize weak formulation.
  CustomWeakForm wf(&f);

  // Initialize boundary conditions
  DefaultEssentialBCNonConst<double> bc_essential("Bdy", &exact_sln);
  EssentialBCs<double> bcs(&bc_essential);

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

  // Initialize approximate solution.
  Solution<double> 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, 440, 350));
  sview.show_mesh(false);
  sview.fix_scale_width(50);
  Views::OrderView  oview("Polynomial orders", new Views::WinGeom(450, 0, 420, 350));

  // DOF and CPU convergence graphs.
  SimpleGraph graph_dof_est, graph_cpu_est, graph_dof_exact, graph_cpu_exact;

  // Time measurement.
  Hermes::TimePeriod cpu_time;

  // Adaptivity loop:
  int as = 1; bool done = false;
  do
  {
    cpu_time.tick();

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

    info("---- Adaptivity step %d (%d DOF):", as, ndof_ref);
    cpu_time.tick();
    
    info("Solving on reference mesh.");
    
    // Assemble the discrete problem.    
    DiscreteProblem<double> dp(&wf, ref_space);
    
    NewtonSolver<double> newton(&dp, matrix_solver);
    newton.set_verbose_output(false);
    
    Solution<double> ref_sln;
    try
    {
      newton.solve();
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.printMsg();
      error("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);
    
    cpu_time.tick();
    verbose("Solution: %g s", cpu_time.last());
    
    // Project the fine mesh solution onto the coarse mesh.
    info("Calculating error estimate and exact error.");
    OGProjection<double>::project_global(&space, &ref_sln, &sln, matrix_solver);

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

    // Calculate exact error.
    double err_exact_rel = Global<double>::calc_rel_error(&sln, &exact_sln, HERMES_H1_NORM) * 100;

    cpu_time.tick();
    verbose("Error calculation: %g s", cpu_time.last());
    
    // Report results.
    info("ndof_coarse: %d, ndof_fine: %d", space.get_num_dofs(), ref_space->get_num_dofs());
    info("err_est_rel: %g%%, err_exact_rel: %g%%", err_est_rel, err_exact_rel);

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

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable bt(butcher_table_type);
  if (bt.is_explicit()) info("Using a %d-stage explicit R-K method.", bt.get_size());
  if (bt.is_diagonally_implicit()) info("Using a %d-stage diagonally implicit R-K method.", bt.get_size());
  if (bt.is_fully_implicit()) info("Using a %d-stage fully implicit R-K method.", bt.get_size());

  // 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 solutions.
  CustomInitialConditionWave E_sln(&mesh);
  ZeroSolution B_sln(&mesh);
  Hermes::vector<Solution<double>*> slns(&E_sln, &B_sln);

  // Initialize the weak formulation.
  CustomWeakFormWave wf(C_SQUARED);
  
  // Initialize boundary conditions
  DefaultEssentialBCConst<double> bc_essential("Perfect conductor", 0.0);
  EssentialBCs<double> bcs_E(&bc_essential);
  EssentialBCs<double> bcs_B;

  // Create x- and y- displacement space using the default H1 shapeset.
  HcurlSpace<double> E_space(&mesh, &bcs_E, P_INIT);
  H1Space<double> B_space(&mesh, &bcs_B, P_INIT);
  Hermes::vector<const Space<double> *> spaces(&E_space, &B_space);
  Hermes::vector<Space<double> *> spaces_mutable(&E_space, &B_space);
  info("ndof = %d.", Space<double>::get_num_dofs(spaces));

  // 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 B_view("Solution B", new WinGeom(0, 405, 400, 350));
  B_view.fix_scale_width(50);

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

  // Time stepping loop.
  double current_time = time_step; int ts = 1;
  do
  {
    // Perform one Runge-Kutta time step according to the selected Butcher's table.
    info("Runge-Kutta time step (t = %g s, time_step = %g s, stages: %d).", 
         current_time, time_step, bt.get_size());
    bool jacobian_changed = false;
    bool verbose = true;
    
    try
    {
      runge_kutta.rk_time_step_newton(current_time, time_step, slns, slns, jacobian_changed, verbose);
    }
    catch(Exceptions::Exception& e)
    {
      e.printMsg();
      error("Runge-Kutta time step failed");
    }

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

    // Update time.
    current_time += time_step;
  
  } while (current_time < T_FINAL);

  // Wait for the view to be closed.
  View::wait();

  return 0;
}

int main(int argc, char* argv[])
{
  // Define nonlinear thermal conductivity lambda(u) via a cubic spline.
  // Step 1: Fill the x values and use lambda_macro(u) = 1 + u^4 for the y values.
#define lambda_macro(x) (1 + Hermes::pow(x, 4))
  Hermes::vector<double> lambda_pts(-2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0);
  Hermes::vector<double> lambda_val;
  for (unsigned int i = 0; i < lambda_pts.size(); i++) 
    lambda_val.push_back(lambda_macro(lambda_pts[i]));
  // Step 2: Create the cubic spline (and plot it for visual control). 
  double bc_left = 0.0;
  double bc_right = 0.0;
  bool first_der_left = false;
  bool first_der_right = false;
  bool extrapolate_der_left = true;
  bool extrapolate_der_right = true;
  CubicSpline lambda(lambda_pts, lambda_val, bc_left, bc_right, first_der_left, first_der_right,
                     extrapolate_der_left, extrapolate_der_right);
  Hermes::Mixins::Loggable::Static::info("Saving cubic spline into a Pylab file spline.dat.");
  // The interval of definition of the spline will be 
  // extended by "interval_extension" on both sides.
  double interval_extension = 3.0;
  lambda.calculate_coeffs();
  lambda.plot("spline.dat", interval_extension);

  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2D mloader;
  mloader.load("square.mesh", mesh);

  // Perform initial mesh refinements.
  for(int i = 0; i < INIT_GLOB_REF_NUM; i++) mesh->refine_all_elements();
  mesh->refine_towards_boundary("Bdy", INIT_BDY_REF_NUM);

  // Initialize boundary conditions.
  CustomEssentialBCNonConst bc_essential("Bdy");
  EssentialBCs<double> bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  SpaceSharedPtr<double> space(new H1Space<double>(mesh, &bcs, P_INIT));
  int ndof = space->get_num_dofs();
  Hermes::Mixins::Loggable::Static::info("ndof: %d", ndof);

  // Initialize the weak formulation.
  Hermes2DFunction<double> src(-heat_src);
  DefaultWeakFormPoisson<double> wf(HERMES_ANY, &lambda, &src);

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

  // 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];
  MeshFunctionSharedPtr<double> init_sln(new CustomInitialCondition(mesh));
  OGProjection<double>::project_global(space, init_sln, coeff_vec); 

  // Initialize Newton solver.
  NewtonSolver<double> newton(&dp);
  newton.set_max_allowed_iterations(NEWTON_MAX_ITER);
  newton.set_tolerance(NEWTON_TOL, Hermes::Solvers::ResidualNormAbsolute);
    
  // Perform Newton's iteration.
  try
  {
    newton.solve(coeff_vec);
  }
  catch(std::exception& e)
  {
    std::cout << e.what();
  }

  // Translate the resulting coefficient vector into a Solution.
  MeshFunctionSharedPtr<double> sln(new Solution<double>);
  Solution<double>::vector_to_solution(newton.get_sln_vector(), space, sln);

  // Get info about time spent during assembling in its respective parts.
  //dp.get_all_profiling_output(std::cout);

  // Clean up.
  delete [] coeff_vec;

  // Visualise the solution and mesh.
  ScalarView s_view("Solution", new WinGeom(0, 0, 440, 350));
  s_view.show_mesh(false);
  s_view.show(sln);
  OrderView o_view("Mesh", new WinGeom(450, 0, 400, 350));
  o_view.show(space);

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

int main(int argc, char* argv[])
{
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("square.mesh", &mesh);
  
  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, P_INIT);

  // Initialize the weak formulation.
  WeakForm<double> wf_dummy;

  // Initialize coarse and reference mesh solution.
  Solution<double> sln;
  ExactSolutionCustom* ref_sln = NULL;

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

  // Initialize views.
  ScalarView sview("Scalar potential Phi", new WinGeom(0, 0, 610, 300));
  sview.fix_scale_width(40);
  sview.show_mesh(false);
  OrderView  oview("Mesh", new WinGeom(620, 0, 600, 300));

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

  // Adapt<double>ivity loop:
  int as = 1; bool done = false;
  do
  {
    info("---- Adapt<double>ivity step %d:", as);

    // Construct globally refined reference mesh and setup reference space.
    Space<double>* ref_space = Space<double>::construct_refined_space(&space);

    // Assign the function f() to the fine mesh.
    info("Assigning f() to the fine mesh.");
    if(ref_sln != NULL) delete ref_sln;
    ref_sln = new ExactSolutionCustom(ref_space->get_mesh());

    // Time measurement.
    cpu_time.tick();
    
    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection<double>::project_global(&space, ref_sln, &sln, matrix_solver); 
   
    // View the coarse mesh solution and polynomial orders.
    sview.show(&sln);
    oview.show(&space);

    // Calculate element errors and total error estimate.
    info("Calculating exact error."); 
    Adapt<double>* adaptivity = new Adapt<double>(&space);
    // Note: the error estimate is now equal to the exact error.
    double err_exact_rel = adaptivity->calc_err_est(&sln, ref_sln) * 100;

    // Report results.
    info("ndof_coarse: %d, ndof_fine: %d, err_exact_rel: %g%%", 
      Space<double>::get_num_dofs(&space), Space<double>::get_num_dofs(ref_space), err_exact_rel);

    // Time measurement.
    cpu_time.tick();

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

    // If err_exact_rel too large, adapt the mesh.
    if (err_exact_rel < ERR_STOP) done = true;
    else 
    {
      info("Adapting coarse mesh.");
      done = adaptivity->adapt(&selector, THRESHOLD, STRATEGY, MESH_REGULARITY);
      
      // Increase the counter of performed adaptivity steps.
      if (done == false)  as++;
    }
    if (Space<double>::get_num_dofs(&space) >= NDOF_STOP) done = true;

    // Clean up.
    delete adaptivity;
    if (done == false)
      delete ref_space->get_mesh();
    delete ref_space;
  }
  while (done == false);
  
  verbose("Total running time: %g s", cpu_time.accumulated());

  // Show the reference solution - the final result.
//.........这里部分代码省略.........

int main(int argc, char* args[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new 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.
  SpaceSharedPtr<double> fine_space(new L2Space<double>(mesh, USE_TAYLOR_SHAPESET ? std::max(P_INIT, 2) : P_INIT, (USE_TAYLOR_SHAPESET ? (Shapeset*)(new L2ShapesetTaylor) : (Shapeset*)(new L2ShapesetLegendre))));

  // Initialize refinement selector.
  L2ProjBasedSelector<double> selector(CAND_LIST);
  selector.set_error_weights(1., 1., 1.);

  MeshFunctionSharedPtr<double> sln(new Solution<double>);
  MeshFunctionSharedPtr<double> refsln(new Solution<double>);

  // Initialize the weak formulation.
  WeakFormSharedPtr<double> wf(new CustomWeakForm("Bdy_bottom_left", mesh));
  ScalarView view1("Solution", new WinGeom(900, 0, 450, 350));
  view1.fix_scale_width(60);

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

  adaptivity.set_space(fine_space);

  int as = 1; bool done = false;
  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<double>::ReferenceSpaceCreator refspace_creator(fine_space, ref_mesh, 0);
    SpaceSharedPtr<double> refspace = refspace_creator.create_ref_space();

    try
    {
      linear_solver.set_space(refspace);
      linear_solver.solve();

      if (USE_TAYLOR_SHAPESET)
      {
        PostProcessing::VertexBasedLimiter limiter(refspace, linear_solver.get_sln_vector(), P_INIT);
        refsln = limiter.get_solution();
      }
      else
      {
        Solution<double>::vector_to_solution(linear_solver.get_sln_vector(), refspace, refsln);
      }

      view1.show(refsln);
      OGProjection<double>::project_global(fine_space, refsln, sln, HERMES_L2_NORM);
    }
    catch (Exceptions::Exception& e)
    {
      std::cout << e.info();
    }
    catch (std::exception& e)
    {
      std::cout << e.what();
    }

    // Calculate element errors and total error estimate.
    errorCalculator.calculate_errors(sln, refsln);
    double err_est_rel = errorCalculator.get_total_error_squared() * 100;

    std::cout << "Error: " << err_est_rel << "%." << std::endl;

    // If err_est_rel too large, adapt the mesh.
    if (err_est_rel < ERR_STOP)
      done = true;
    else
      done = adaptivity.adapt(&selector);
    as++;
  } while (done == false);

  // Wait for keyboard or mouse input.
  View::wait();
  return 0;
}

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

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

    mesh->refine_by_criterion(criterion, INIT_REF_CRITERION);

    ScalarView view1("Solution - Discontinuous Galerkin FEM", new WinGeom(900, 0, 450, 350));
    ScalarView view2("Solution - Standard continuous FEM", new WinGeom(900, 400, 450, 350));

    if(WANT_DG)
    {
        // Create an L2 space.
        SpaceSharedPtr<double> space_l2(new L2Space<double>(mesh, P_INIT));

        // Initialize the solution.
        MeshFunctionSharedPtr<double> sln_l2(new Solution<double>);

        // Initialize the weak formulation.
        CustomWeakForm wf_l2(BDY_BOTTOM_LEFT);

        // Initialize the FE problem.
        DiscreteProblem<double> dp_l2(&wf_l2, space_l2);
        dp_l2.set_linear();

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

        Hermes::Mixins::Loggable::Static::info("Assembling Discontinuous Galerkin (nelem: %d, ndof: %d).", mesh->get_num_active_elements(), space_l2->get_num_dofs());

        // Solve the linear system. If successful, obtain the solution.
        Hermes::Mixins::Loggable::Static::info("Solving Discontinuous Galerkin.");
        try
        {
            linear_solver.solve();
            if(DG_SHOCK_CAPTURING)
            {
                FluxLimiter flux_limiter(FluxLimiter::Kuzmin, linear_solver.get_sln_vector(), space_l2, true);

                flux_limiter.limit_second_orders_according_to_detector();

                flux_limiter.limit_according_to_detector();

                flux_limiter.get_limited_solution(sln_l2);

                view1.set_title("Solution - limited Discontinuous Galerkin FEM");
            }
            else
                Solution<double>::vector_to_solution(linear_solver.get_sln_vector(), space_l2, sln_l2);

            // View the solution.
            view1.show(sln_l2);
        }
        catch(std::exception& e)
        {
            std::cout << e.what();

        }
    }
    if(WANT_FEM)
    {
        // Create an H1 space.
        SpaceSharedPtr<double> space_h1(new H1Space<double>(mesh, P_INIT));

        // Initialize the solution.
        MeshFunctionSharedPtr<double> sln_h1(new Solution<double>);

        // Initialize the weak formulation.
        CustomWeakForm wf_h1(BDY_BOTTOM_LEFT, false);

        // Initialize the FE problem.
        DiscreteProblem<double> dp_h1(&wf_h1, space_h1);
        dp_h1.set_linear();

        Hermes::Mixins::Loggable::Static::info("Assembling Continuous FEM (nelem: %d, ndof: %d).", mesh->get_num_active_elements(), space_h1->get_num_dofs());

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

        // Solve the linear system. If successful, obtain the solution.
        Hermes::Mixins::Loggable::Static::info("Solving Continuous FEM.");
        try
        {
            linear_solver.solve();
            Solution<double>::vector_to_solution(linear_solver.get_sln_vector(), space_h1, sln_h1);

            // View the solution.
            view2.show(sln_h1);
        }
        catch(std::exception& e)
        {
            std::cout << e.what();
        }
    }
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable bt(butcher_table_type);
  if (bt.is_explicit()) info("Using a %d-stage explicit R-K method.", bt.get_size());
  if (bt.is_diagonally_implicit()) info("Using a %d-stage diagonally implicit R-K method.", bt.get_size());
  if (bt.is_fully_implicit()) info("Using a %d-stage fully implicit R-K method.", bt.get_size());

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

  // Perform initial mesh refinements.
  for(int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
  mesh.refine_towards_boundary(BDY_RIGHT, 2);
  mesh.refine_towards_boundary(BDY_FIRE, INIT_REF_NUM_BDY);

  // Initialize essential boundary conditions (none).
  EssentialBCs<double> bcs;

  // Initialize an H1 space with default shapeset.
  H1Space<double> space(&mesh, &bcs, P_INIT);
  int ndof = Space<double>::get_num_dofs(&space);
  info("ndof = %d.", ndof);
 
  // Previous and next time level solutions.
  ConstantSolution<double> sln_time_prev(&mesh, TEMP_INIT);
  ZeroSolution sln_time_new(&mesh);
  ConstantSolution<double> time_error_fn(&mesh, 0.0);

  // Initialize the weak formulation.
  double current_time = 0;
  CustomWeakFormHeatRK wf(BDY_FIRE, BDY_AIR, ALPHA_FIRE, ALPHA_AIR,
                          RHO, HEATCAP, TEMP_EXT_AIR, TEMP_INIT, &current_time);

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

  // Initialize views.
  ScalarView Tview("Temperature", new WinGeom(0, 0, 1500, 400));
  Tview.fix_scale_width(40);
  ScalarView eview("Temporal error", new WinGeom(0, 450, 1500, 400));
  eview.fix_scale_width(40);

  // Graph for time step history.
  SimpleGraph time_step_graph;
  info("Time step history will be saved to file time_step_history.dat.");

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

  // Time stepping loop:
  int ts = 1;
  do 
  {
    // Perform one Runge-Kutta time step according to the selected Butcher's table.
    info("Runge-Kutta time step (t = %g s, tau = %g s, stages: %d).", 
         current_time, time_step, bt.get_size());
    bool jacobian_changed = false;
    bool verbose = true;

    try
    {
      runge_kutta.rk_time_step_newton(current_time, time_step, &sln_time_prev, &sln_time_new, 
                                  &time_error_fn, !jacobian_changed, false, verbose);
    }
    catch(Exceptions::Exception& e)
    {
      e.printMsg();
      error("Runge-Kutta time step failed");
    }

    // Plot error function.
    char title[100];
    sprintf(title, "Temporal error, t = %g", current_time);
    eview.set_title(title);
    AbsFilter abs_tef(&time_error_fn);
    eview.show(&abs_tef, HERMES_EPS_VERYHIGH);

    // Calculate relative time stepping error and decide whether the 
    // time step can be accepted. If not, then the time step size is 
    // reduced and the entire time step repeated. If yes, then another
    // check is run, and if the relative error is very low, time step 
    // is increased.
    double rel_err_time = Global<double>::calc_norm(&time_error_fn, HERMES_H1_NORM) 
                          / Global<double>::calc_norm(&sln_time_new, HERMES_H1_NORM) * 100;
    info("rel_err_time = %g%%", rel_err_time);
    if (rel_err_time > TIME_TOL_UPPER) {
      info("rel_err_time above upper limit %g%% -> decreasing time step from %g to %g and restarting time step.", 
           TIME_TOL_UPPER, time_step, time_step * TIME_STEP_DEC_RATIO);
      time_step *= TIME_STEP_DEC_RATIO;
      continue;
    }
    if (rel_err_time < TIME_TOL_LOWER) {
      info("rel_err_time = below lower limit %g%% -> increasing time step from %g to %g", 
           TIME_TOL_UPPER, time_step, time_step * TIME_STEP_INC_RATIO);
      time_step *= TIME_STEP_INC_RATIO;
    }

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