C++ DiscreteProblem::assemble方法代码示例

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

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_neumann(Hermes::vector<int>(BDY_LEFT_RIGHT, BDY_TOP_BOTTOM));

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, P_INIT);
  int ndof = Space::get_num_dofs(&space);
  info("ndof = %d", ndof);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form_vol));
  wf.add_vector_form(callback(linear_form_vol));
  wf.add_vector_form_surf(callback(linear_form_surf_right), 2);
  wf.add_vector_form_surf(callback(linear_form_surf_left), 2);

  Solution sln; 

  // NON-ADAPTIVE VERSION
  
  // Initialize the linear problem.
  bool is_linear = true;
  DiscreteProblem* dp = new DiscreteProblem(&wf, &space, is_linear);

  // Select matrix solver.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Assemble stiffness matrix and rhs.
  dp->assemble(matrix, rhs);
 
  // Solve the linear system of the reference problem. If successful, obtain the solutions.
  if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
  else error ("Matrix solver failed.\n");

  // Visualize the solution.
  ScalarView view("Solution", new WinGeom(0, 0, 440, 350));
  view.show(&sln);

  // Calculate error wrt. exact solution.
  Solution sln_exact;
  sln_exact.set_exact(&mesh, exact);
  double err = calc_abs_error(&sln, &sln_exact, HERMES_H1_NORM);
  printf("err = %g, err_squared = %g\n\n", err, err*err);
 
  // Wait for all views to be closed.
  View::wait();
  return 0;
}

void OGProjection::project_internal(Hermes::Tuple<Space *> spaces, WeakForm* wf, 
                                    scalar* target_vec, MatrixSolverType matrix_solver)
{
  _F_
  unsigned int n = spaces.size();

  // sanity checks
  if (n <= 0 || n > 10) error("Wrong number of projected functions in project_internal().");
  for (unsigned int i = 0; i < n; i++) if(spaces[i] == NULL) error("this->spaces[%d] == NULL in project_internal().", i);
  if (spaces.size() != n) error("Number of spaces must match number of projected functions in project_internal().");

  // this is needed since spaces may have their DOFs enumerated only locally.
  int ndof = Space::assign_dofs(spaces);

  // Initialize DiscreteProblem.
  bool is_linear = true;
  DiscreteProblem* dp = new DiscreteProblem(wf, spaces, is_linear);

  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  dp->assemble(matrix, rhs, false);

  // Calculate the coefficient vector.
  bool solved = solver->solve();
  scalar* coeffs;
  if (solved) 
    coeffs = solver->get_solution();

  if (target_vec != NULL) 
    for (int i=0; i<ndof; i++) target_vec[i] = coeffs[i];
    
  delete solver;
  delete matrix;
  delete rhs;
  delete dp;
  delete wf;
}

int main() {
  // Three macroelements are defined above via the interfaces[] array.
  // poly_orders[]... initial poly degrees of macroelements.
  // material_markers[]... material markers of macroelements.
  // subdivisions[]... equidistant subdivision of macroelements.
  int poly_orders[N_MAT] = {P_init_inner, P_init_outer, P_init_reflector };
  int material_markers[N_MAT] = {Marker_inner, Marker_outer, Marker_reflector };
  int subdivisions[N_MAT] = {N_subdiv_inner, N_subdiv_outer, N_subdiv_reflector };

  // Create space.
  Space* space = new Space(N_MAT, interfaces, poly_orders, material_markers, subdivisions, N_GRP, N_SLN);
  // Enumerate basis functions, info for user.
  info("N_dof = %d", Space::get_num_dofs(space));

  // Initial approximation: u = 1.
  double K_EFF_old;
  double init_val = 1.0;
  set_vertex_dofs_constant(space, init_val, 0);
  
  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(jacobian_vol_inner, NULL, Marker_inner);
  wf.add_matrix_form(jacobian_vol_outer, NULL, Marker_outer);
  wf.add_matrix_form(jacobian_vol_reflector, NULL, Marker_reflector);
  wf.add_vector_form(residual_vol_inner, NULL, Marker_inner);
  wf.add_vector_form(residual_vol_outer, NULL, Marker_outer);
  wf.add_vector_form(residual_vol_reflector, NULL, Marker_reflector);
  wf.add_vector_form_surf(residual_surf_left, BOUNDARY_LEFT);
  wf.add_matrix_form_surf(jacobian_surf_right, BOUNDARY_RIGHT);
  wf.add_vector_form_surf(residual_surf_right, BOUNDARY_RIGHT);

  // Initialize the FE problem.
  bool is_linear = false;
  DiscreteProblem *dp = new DiscreteProblem(&wf, space, is_linear);

  // Source iteration (power method).
  for (int i = 0; i < Max_SI; i++)
  {	
    // Obtain fission source.
    int current_solution = 0, previous_solution = 1;
    copy_dofs(current_solution, previous_solution, space);
		
    // Obtain the number of degrees of freedom.
    int ndof = Space::get_num_dofs(space);

    // Fill vector coeff_vec using dof and coeffs arrays in elements.
  double *coeff_vec = new double[Space::get_num_dofs(space)];
    solution_to_vector(space, coeff_vec);
  
    // Set up the solver, matrix, and rhs according to the solver selection.
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
  
    int it = 1;
  while (1) {
    // Obtain the number of degrees of freedom.
    int ndof = Space::get_num_dofs(space);

      // Assemble the Jacobian matrix and residual vector.
      dp->assemble(matrix, rhs);

      // Calculate the l2-norm of residual vector.
      double res_norm = 0;
      for(int i=0; i<ndof; i++) res_norm += rhs->get(i)*rhs->get(i);
      res_norm = sqrt(res_norm);

      // Info for user.
      info("---- Newton iter %d, residual norm: %.15f", it, res_norm);

      // If l2 norm of the residual vector is within tolerance, then quit.
      // NOTE: at least one full iteration forced
      //       here because sometimes the initial
      //       residual on fine mesh is too small.
      if(res_norm < NEWTON_TOL && it > 1) break;

      // Multiply the residual vector with -1 since the matrix 
      // equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
      for(int i=0; i<ndof; i++) rhs->set(i, -rhs->get(i));

      // Solve the linear system.
      if(!solver->solve())
        error ("Matrix solver failed.\n");

      // Add \deltaY^{n+1} to Y^n.
      for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];

      // If the maximum number of iteration has been reached, then quit.
      if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
      
      // Copy coefficients from vector y to elements.
      vector_to_solution(coeff_vec, space);

      it++;
    }
    
    // Cleanup.
    delete matrix;
    delete rhs;
    delete solver;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    // Instantiate a class with global functions.
    Hermes2D hermes2d;

    // Load the mesh.
    Mesh mesh;
    H2DReader 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(&mesh, EPSILON);

    // Define right-hand side.
    CustomRightHandSide rhs(EPSILON);

    // Initialize the weak formulation.
    CustomWeakForm wf(&rhs);

    // Initialize boundary conditions
    DefaultEssentialBCNonConst bc_essential(BDY_DIRICHLET, &exact);
    EssentialBCs bcs(&bc_essential);

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

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

    // Initialize views.
    ScalarView sview("Solution", new WinGeom(0, 0, 460, 350));
    sview.show_mesh(false);
    sview.fix_scale_width(70);
    OrderView  oview("Polynomial orders", new WinGeom(470, 0, 410, 350));

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

    // Time measurement.
    TimePeriod cpu_time;
    cpu_time.tick();

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

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

        // Set up the solver, matrix, and rhs according to the solver selection.
        SparseMatrix* matrix = create_matrix(matrix_solver);
        Vector* rhs = create_vector(matrix_solver);
        Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

        // Assemble the reference problem.
        info("Solving on reference mesh.");
        bool is_linear = true;
        DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
        dp->assemble(matrix, rhs);

        // Time measurement.
        cpu_time.tick();

        // Solve the linear system of the reference problem. If successful, obtain the solution.
        Solution ref_sln;
        if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
        else error ("Matrix solver failed.\n");

        // Time measurement.
        cpu_time.tick();

        // Project the fine mesh solution onto the coarse mesh.
        Solution sln;
        info("Projecting reference solution on coarse mesh.");
        OGProjection::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 error estimate and exact error.");
        Adapt* adaptivity = new Adapt(&space);
        double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln) * 100;

        // Calculate exact error for each solution component.
        double err_exact_rel = hermes2d.calc_rel_error(&sln, &exact, HERMES_H1_NORM) * 100;

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

        // Time measurement.
        cpu_time.tick();
//.........这里部分代码省略.........

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_BOTTOM);
  bc_types.add_bc_neumann(Hermes::Tuple<int>(BDY_RIGHT, BDY_TOP, BDY_LEFT));

  // Enter Dirichlet boudnary values.
  BCValues bc_values;
  bc_values.add_zero(BDY_BOTTOM);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, &bc_values, P_INIT);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
  wf.add_vector_form(callback(linear_form));
  wf.add_vector_form_surf(callback(linear_form_surf), BDY_TOP);

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

  // Set exact solution.
  ExactSolution exact(&mesh, fndd);

  // Initialize views.
  ScalarView sview("Solution", new WinGeom(0, 0, 440, 350));
  sview.show_mesh(false);
  OrderView  oview("Polynomial orders", new WinGeom(450, 0, 400, 350));

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

  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

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

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

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system of the reference problem. If successful, obtain the solution.
    Solution ref_sln;
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");

    // Time measurement.
    cpu_time.tick();

    // Project the fine mesh solution onto the coarse mesh.
    Solution sln;
    info("Projecting reference solution on coarse mesh.");
    OGProjection::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 error estimate and exact error.");
    Adapt* adaptivity = new Adapt(&space, HERMES_H1_NORM);
    bool solutions_for_adapt = true;
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL) * 100;

    // Calculate exact error.   
    solutions_for_adapt = false;
    double err_exact_rel = adaptivity->calc_err_exact(&sln, &exact, solutions_for_adapt, HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL) * 100;

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

    // Time measurement.
    cpu_time.tick();

    // Add entry to DOF and CPU convergence graphs.
//.........这里部分代码省略.........

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

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

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, bc_types, essential_bc_values, P_INIT);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
  wf.add_vector_form(callback(linear_form));

  // Initialize coarse and reference mesh solution.
  Solution sln, ref_sln;
  
  // Initialize refinement selector.
  H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
  
  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;
  
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

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

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

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();
    
    // Solve the linear system of the reference problem. 
    // If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver); 

    // Calculate element errors and total error estimate.
    info("Calculating error estimate."); 
    Adapt* adaptivity = new Adapt(&space, HERMES_H1_NORM);
    bool solutions_for_adapt = true;
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL) * 100;

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

    // Time measurement.
    cpu_time.tick();

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

    // If err_est too large, adapt the mesh.
    if (err_est_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::get_num_dofs(&space) >= NDOF_STOP) done = true;

    // Clean up.
    delete solver;
    delete matrix;
    delete rhs;
    delete adaptivity;
    if(done == false) delete ref_space->get_mesh();
//.........这里部分代码省略.........

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(Hermes::Tuple<int>(OUTER_BDY, STATOR_BDY));

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_const(STATOR_BDY, VOLTAGE);
  bc_values.add_const(OUTER_BDY, 0.0);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, &bc_values, P_INIT);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(biform1), HERMES_SYM, MATERIAL_1);
  wf.add_matrix_form(callback(biform2), HERMES_SYM, MATERIAL_2);

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

  // Initialize views.
  ScalarView sview("Solution", new WinGeom(0, 0, 410, 600));
  sview.fix_scale_width(50);
  sview.show_mesh(false);
  OrderView  oview("Polynomial orders", new WinGeom(420, 0, 400, 600));
  
  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;
  
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

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

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

    // Initialize matrix solver.
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

    // Assemble reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();
    
    // Solve the linear system of the reference problem. 
    // If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver); 

    // Time measurement.
    cpu_time.tick();
   
    // View the coarse mesh solution and polynomial orders.
    sview.show(&sln);
    oview.show(&space);

    // Skip visualization time.
    cpu_time.tick(HERMES_SKIP);

    // Calculate element errors and total error estimate.
    info("Calculating error estimate."); 
    Adapt* adaptivity = new Adapt(&space);
    bool solutions_for_adapt = true;
    // In the following function, the Boolean parameter "solutions_for_adapt" determines whether 
    // the calculated errors are intended for use with adaptivity (this may not be the case, for example,
    // when error wrt. an exact solution is calculated). The default value is solutions_for_adapt = true, 
    // The last parameter "error_flags" determine whether the total and element errors are treated as 
    // absolute or relative. Its default value is error_flags = HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL.
    // In subsequent examples and benchmarks, these two parameters will be often used with 
    // their default values, and thus they will not be present in the code explicitly.
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt, 
                         HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL) * 100;

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

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

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("domain2.mesh", &mesh);

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

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, bc_types, essential_bc_values, P_INIT);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form_iron), HERMES_SYM, 3);
  wf.add_matrix_form(callback(bilinear_form_wire), HERMES_SYM, 2);
  wf.add_matrix_form(callback(bilinear_form_air), HERMES_SYM, 1);
  wf.add_vector_form(callback(linear_form_wire), 2);

  // Initialize coarse and reference mesh solution.
  Solution sln, ref_sln;

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

  // Initialize views.
  ScalarView sview("Solution", new WinGeom(0, 0, 600, 350));
  sview.show_mesh(false);
  OrderView  oview("Polynomial orders", new WinGeom(610, 0, 520, 350));
  
  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;

  initialize_solution_environment(matrix_solver, argc, argv);
  
  // Adaptivity loop:
  int as = 1; 
  bool done = false;
  do
  {
    info("---- Adaptivity step %d:", as);

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

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
    
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    
    if (matrix_solver == SOLVER_AZTECOO) {
      ((AztecOOSolver*) solver)->set_solver(iterative_method);
      ((AztecOOSolver*) solver)->set_precond(preconditioner);
      // Using default iteration parameters (see solver/aztecoo.h).
    }
    
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();
    
    // Solve the linear system of the reference problem. If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");
  
    // Time measurement.
    cpu_time.tick();

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection::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 error estimate."); 
    Adapt* adaptivity = new Adapt(&space, HERMES_H1_NORM);
    bool solutions_for_adapt = true;
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL) * 100;

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

    // Time measurement.
    cpu_time.tick();

    // Add entry to DOF and CPU convergence graphs.
    graph_dof.add_values(Space::get_num_dofs(&space), err_est_rel);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Instantiate a class with global functions.
  Hermes2D hermes2d;

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;

  switch (PARAM) {
    case 0: mloader.load("../geom0.mesh", &mesh); break;
    case 1: mloader.load("../geom1.mesh", &mesh); break;
    case 2: mloader.load("../geom2.mesh", &mesh); break;
    case 3: mloader.load("../geom3.mesh", &mesh); break;
    default: error("Admissible values of PARAM are 0, 1, 2, 3.");
  }

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

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

  // Initialize the weak formulation.
  DefaultWeakFormLaplace wf;
  
  // Initialize boundary conditions
  DefaultEssentialBCNonConst bc_essential(BDY_DIRICHLET, &exact);
  EssentialBCs bcs(&bc_essential);

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

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

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

  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

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

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

    // Set up the solver, matrix, and rhs according to the solver selection.
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system of the reference problem. If successful, obtain the solution.
    Solution ref_sln;
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");

    // Time measurement.
    cpu_time.tick();

    // Project the fine mesh solution onto the coarse mesh.
    Solution sln;
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver);

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

    // Calculate exact error.
    double err_exact_rel = hermes2d.calc_rel_error(&sln, &exact, HERMES_H1_NORM) * 100;

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

    // Time measurement.
    cpu_time.tick();

    // Add entry to DOF and CPU convergence graphs.
    graph_dof.add_values(Space::get_num_dofs(&space), err_est_rel);
    graph_dof.save("conv_dof_est.dat");
    graph_cpu.add_values(cpu_time.accumulated(), err_est_rel);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  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, graph_dof_exact, graph_cpu_exact;
  
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
  
  // Adaptivity loop:
  int as = 1; 
  bool done = false;
  Space* actual_sln_space;
  do
  {
    info("---- Adaptivity step %d:", as);

    if (STRATEGY == -1)
      actual_sln_space = space;
    else
      // Construct globally refined reference mesh and setup reference space.
      actual_sln_space = Space::construct_refined_space(space, ORDER_INCREASE);

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, actual_sln_space, is_linear);
    dp->assemble(matrix, rhs);
    
    // Solve the linear system of the reference problem. 
    // If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), actual_sln_space, &ref_sln);
    else error ("Matrix solver failed.\n");
    
    // Instantiate adaptivity and error calculation driver. Space is used only for adaptivity, it is ignored when 
    // STRATEGY == -1 and only the exact error is calculated by this object.
    Adapt* adaptivity = new Adapt(space, norm);
    
    // Calculate exact error.
    bool solutions_for_adapt = false;
    double err_exact_rel = calc_rel_error(&sln, &exact, HERMES_H1_NORM) * 100;
    info("ndof_fine: %d, err_exact_rel: %g%%", Space::get_num_dofs(actual_sln_space), err_exact_rel);

    // Time measurement.
    cpu_time.tick();
    
    // View the fine mesh solution and polynomial orders.
    sview.show(&ref_sln);
    oview.show(actual_sln_space);
    
    // Skip visualization time.
    cpu_time.tick(HERMES_SKIP);

    if (STRATEGY == -1) done = true;  // Do not adapt.
    else
    {  
      // Project the fine mesh solution onto the coarse mesh.
      info("Projecting reference solution on coarse mesh.");
      OGProjection::project_global(space, &ref_sln, &sln, matrix_solver, norm); 

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh basemesh, T_mesh, M_mesh;
  H2DReader mloader;
  mloader.load("domain2.mesh", &basemesh);

  // Create temperature and moisture meshes.
  // This also initializes the multimesh hp-FEM.
  T_mesh.copy(&basemesh);
  M_mesh.copy(&basemesh);

  // Create H1 spaces with default shapesets.
  H1Space T_space(&T_mesh, temp_bc_type, essential_bc_values_T, P_INIT);
  H1Space M_space(MULTI ? &M_mesh : &T_mesh, moist_bc_type, NULL, P_INIT);

  // Define constant initial conditions.
  info("Setting initial conditions.");
  Solution T_prev, M_prev;
  T_prev.set_const(&T_mesh, TEMP_INITIAL);
  M_prev.set_const(&M_mesh, MOIST_INITIAL);

  // Initialize the weak formulation.
  WeakForm wf(2);
  wf.add_matrix_form(0, 0, callback(bilinear_form_sym_0_0));
  wf.add_matrix_form(0, 1, callback(bilinear_form_sym_0_1));
  wf.add_matrix_form(1, 1, callback(bilinear_form_sym_1_1));
  wf.add_matrix_form(1, 0, callback(bilinear_form_sym_1_0));
  wf.add_vector_form(0, callback(linear_form_0), HERMES_ANY, &T_prev);
  wf.add_vector_form(1, callback(linear_form_1), HERMES_ANY, &M_prev);
  wf.add_matrix_form_surf(0, 0, callback(bilinear_form_surf_0_0_ext), MARKER_EXTERIOR_WALL);
  wf.add_matrix_form_surf(1, 1, callback(bilinear_form_surf_1_1_ext), MARKER_EXTERIOR_WALL);
  wf.add_vector_form_surf(0, callback(linear_form_surf_0_ext), MARKER_EXTERIOR_WALL);
  wf.add_vector_form_surf(1, callback(linear_form_surf_1_ext), MARKER_EXTERIOR_WALL);

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

  // Solutions.
  Solution T_coarse, M_coarse, T_fine, M_fine;

  // Geometry and position of visualization windows.
  WinGeom* T_sln_win_geom = new WinGeom(0, 0, 300, 450);
  WinGeom* M_sln_win_geom = new WinGeom(310, 0, 300, 450);
  WinGeom* T_mesh_win_geom = new WinGeom(620, 0, 280, 450);
  WinGeom* M_mesh_win_geom = new WinGeom(910, 0, 280, 450);

  // Initialize views.
  ScalarView T_sln_view("Temperature", T_sln_win_geom);
  ScalarView M_sln_view("Moisture", M_sln_win_geom);
  OrderView T_order_view("Temperature mesh", T_mesh_win_geom);
  OrderView M_order_view("Moisture mesh", M_mesh_win_geom);

  // Show initial conditions.
  T_sln_view.show(&T_prev);
  M_sln_view.show(&M_prev);
  T_order_view.show(&T_space);
  M_order_view.show(&M_space);

  // Time stepping loop:
  bool verbose = true;  // Print info during adaptivity.
  double comp_time = 0.0;
  static int ts = 1;
  while (CURRENT_TIME < SIMULATION_TIME)
  {
    info("Simulation time = %g s (%d h, %d d, %d y)",
        (CURRENT_TIME + TAU), (int) (CURRENT_TIME + TAU) / 3600,
        (int) (CURRENT_TIME + TAU) / (3600*24), (int) (CURRENT_TIME + TAU) / (3600*24*364));

    // Uniform mesh derefinement.
    if (ts > 1) {
      info("Global mesh derefinement.");
      T_mesh.copy(&basemesh);
      M_mesh.copy(&basemesh);
      T_space.set_uniform_order(P_INIT);
      M_space.set_uniform_order(P_INIT);
    }

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

      // Construct globally refined reference mesh and setup reference space.
      Tuple<Space *>* ref_spaces = construct_refined_spaces(Tuple<Space *>(&T_space, &M_space));

      // Assemble the reference problem.
      info("Solving on reference mesh.");
      bool is_linear = true;
      DiscreteProblem* dp = new DiscreteProblem(&wf, *ref_spaces, is_linear);
      SparseMatrix* matrix = create_matrix(matrix_solver);
      Vector* rhs = create_vector(matrix_solver);
      Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
      dp->assemble(matrix, rhs);

      // Now we can deallocate the previous fine meshes.
      if(as > 1){ delete T_fine.get_mesh(); delete M_fine.get_mesh(); }

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

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

    // Enter boundary markers.
    BCTypes bc_types;
    bc_types.add_bc_dirichlet(BDY_BOTTOM);
    bc_types.add_bc_neumann(Hermes::Tuple<int>(BDY_TOP_NE, BDY_TOP_NW));
    bc_types.add_bc_neumann(Hermes::Tuple<std::string>(BDY_VERTICAL_SE, BDY_VERTICAL_NE, BDY_VERTICAL_NW, BDY_VERTICAL_SW));

    // Enter Dirichlet boundary values.
    BCValues bc_values;
    bc_values.add_zero(BDY_BOTTOM);

    // Create an H1 space with default shapeset.
    H1Space space(&mesh, &bc_types, &bc_values, P_INIT);

    // Initialize the weak formulation.
    WeakForm wf;
    wf.add_matrix_form(callback(bilinear_form_vol_SE), HERMES_UNSYM, SOUTH_EAST);
    wf.add_matrix_form(callback(bilinear_form_vol_NE), HERMES_UNSYM, NORTH_EAST);
    wf.add_matrix_form(callback(bilinear_form_vol_NW), HERMES_UNSYM, NORTH_WEST);
    wf.add_matrix_form(callback(bilinear_form_vol_SW), HERMES_UNSYM, SOUTH_WEST);

    wf.add_vector_form(callback(linear_form_vol));

    wf.add_vector_form_surf(callback(linear_form_surf_VERTICAL_SE), BDY_VERTICAL_SE);
    wf.add_vector_form_surf(callback(linear_form_surf_VERTICAL_NE), BDY_VERTICAL_NE);
    wf.add_vector_form_surf(callback(linear_form_surf_VERTICAL_NW), BDY_VERTICAL_NW);
    wf.add_vector_form_surf(callback(linear_form_surf_VERTICAL_SW), BDY_VERTICAL_SW);
    wf.add_vector_form_surf(callback(linear_form_surf_TOP_NE), BDY_TOP_NE);
    wf.add_vector_form_surf(callback(linear_form_surf_TOP_NW), BDY_TOP_NW);

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

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

    // Time measurement.
    TimePeriod cpu_time;
    cpu_time.tick();

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

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

        // Assemble the reference problem.
        info("Solving on reference mesh.");
        bool is_linear = true;
        DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
        SparseMatrix* matrix = create_matrix(matrix_solver);
        Vector* rhs = create_vector(matrix_solver);
        Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
        dp->assemble(matrix, rhs);

        // Time measurement.
        cpu_time.tick();

        // Solve the linear system of the reference problem. If successful, obtain the solution.
        Solution ref_sln;
        if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
        else error ("Matrix solver failed.\n");

        // Time measurement.
        cpu_time.tick();

        // Project the fine mesh solution onto the coarse mesh.
        Solution sln;
        info("Projecting reference solution on the coarse mesh.");
        OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver);

        // Calculate element errors and total error estimate.
        info("Calculating error estimate.");
        Adapt* adaptivity = new Adapt(&space, HERMES_H1_NORM);
        bool solutions_for_adapt = true;
        double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL) * 100;

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

        // Time measurement.
        cpu_time.tick();

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

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

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

  // Load the mesh.
  Mesh u_mesh, v_mesh;
  H2DReader mloader;
  mloader.load("elasticity.mesh", &u_mesh);

  // Create initial mesh (master mesh).
  v_mesh.copy(&u_mesh);

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_DIRICHLET);

  // Enter Dirichlet boundary values.
  BCValues bc_values_u, bc_values_v;
  bc_values_u.add_function(BDY_DIRICHLET, essential_bc_values_u);
  bc_values_v.add_function(BDY_DIRICHLET, essential_bc_values_v);

  // Create H1 spaces with default shapeset for both displacement components.
  H1Space u_space(&u_mesh, &bc_types, &bc_values_u, P_INIT_U);
  H1Space v_space(&v_mesh, &bc_types, &bc_values_v, P_INIT_V);

  // Initialize the weak formulation.
  WeakForm wf(2);
  wf.add_matrix_form(0, 0, callback(bilinear_form_0_0), HERMES_SYM);
  wf.add_matrix_form(0, 1, callback(bilinear_form_0_1), HERMES_SYM);
  wf.add_matrix_form(1, 1, callback(bilinear_form_1_1), HERMES_SYM);
  //wf.add_vector_form(0, linear_form_u, linear_form_0_ord, HERMES_SYM);
  //wf.add_vector_form(1, linear_form_v, linear_form_1_ord, HERMES_SYM);

  Solution u_sln, v_sln, u_ref_sln, v_ref_sln;

  // Initialize exact solutions.
  ExactSolution u_exact(&u_mesh, u_fndd);
  ExactSolution v_exact(&v_mesh, v_fndd);

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

  // Initialize views.
  ScalarView s_view_u("Solution for u", new WinGeom(0, 0, 440, 350));
  s_view_u.show_mesh(false);
  OrderView  o_view_u("Mesh for u", new WinGeom(450, 0, 420, 350));

  ScalarView s_view_v("Solution for v", new WinGeom(880, 0, 440, 350));
  s_view_v.show_mesh(false);
  OrderView  o_view_v("Mesh for v", new WinGeom(1330, 0, 420, 350));

  /*  
  ScalarView sview_u_exact("", new WinGeom(550, 0, 500, 400));
  sview_u_exact.fix_scale_width(50);
  ScalarView sview_v_exact("", new WinGeom(550, 500, 500, 400));
  sview_v_exact.fix_scale_width(50);
  char title[100];
  */

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

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

    // Construct globally refined reference mesh and setup reference space.
    Hermes::vector<Space *>* ref_spaces = construct_refined_spaces(Hermes::vector<Space *>(&u_space, &v_space));

    // Initialize matrix solver.
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    solver->set_factorization_scheme(HERMES_REUSE_MATRIX_REORDERING);

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, *ref_spaces, is_linear);
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system of the reference problem. If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solutions(solver->get_solution(), *ref_spaces,
                                                      Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln));
    else error ("Matrix solver failed.\n");

    // Time measurement.
    cpu_time.tick();
//.........这里部分代码省略.........

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

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_DIRICHLET);

  // Enter Dirichlet boudnary values.
  BCValues bc_values;
  bc_values.add_function(BDY_DIRICHLET, essential_bc_values);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, &bc_values, P_INIT);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
  wf.add_vector_form(callback(linear_form));
  
  // Initialize approximate solution.
  Solution sln;

  // Set exact solution.
  ExactSolution exact(&mesh, fndd);

  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();
  
  // Initialize the matrix solver.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Adaptivity loop:
  int as = 1;
  bool done = false;
  double err_exact_rel;
  do
  {
    info("---- Adaptivity step %d:", as);

    // Assemble the discrete problem.
    info("Solving.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, &space, is_linear);
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system. If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
    else error ("Matrix solver failed.\n");

    // Calculate element errors and total error estimate.
    info("Calculating error estimate and exact error.");
    BasicKellyAdapt* adaptivity = new BasicKellyAdapt(&space);
    
    if (USE_RESIDUAL_ESTIMATOR)
      adaptivity->add_error_estimator_vol(callback(residual_estimator));
    
    double err_est_rel = adaptivity->calc_err_est(&sln) * 100;  
    err_exact_rel = adaptivity->calc_err_exact(&sln, &exact) * 100;   
    
    // Time measurement.
    cpu_time.tick();
    
    // Report results.
    info("ndof_coarse: %d", Space::get_num_dofs(&space));
    info("err_est_rel: %1.15g%%, err_exact_rel: %1.15g%%", err_est_rel, err_exact_rel);
    
    // This is to ensure that the two possible approaches to interface error estimators accumulation give
    // same results.
    KellyTypeAdapt* adaptivity2 = new KellyTypeAdapt(&space, Hermes::vector<ProjNormType>(), false);
    adaptivity2->disable_aposteriori_interface_scaling();
    adaptivity2->add_error_estimator_surf(callback(interface_estimator));
    
    double err_est_rel2 = adaptivity2->calc_err_est(&sln) * 100;  
    double err_exact_rel2 = adaptivity2->calc_err_exact(&sln, &exact) * 100;
    
    info("err_est_rel_2: %1.15g%%, err_exact_rel_2: %1.15g%%", err_est_rel2, err_exact_rel2);
    
    if (fabs(err_est_rel2 - err_est_rel) >= 1e-13 || fabs(err_exact_rel2 - err_exact_rel) >= 1e-13)
    {
      info("err_est_rel diff: %g, err_exact_rel diff: %g", 
           std::abs(err_est_rel2 - err_est_rel), std::abs(err_exact_rel2 - err_exact_rel));
      err_est_rel = 0;      // to immediately exit the adaptivity loop
      err_exact_rel = 1e20; // to fail the test
    }
     
    // Time measurement.
//.........这里部分代码省略.........

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

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("domain2.mesh", &mesh);

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(Hermes::Tuple<int>(BDY_RIGHT, BDY_TOP, BDY_LEFT));
  bc_types.add_bc_neumann(BDY_BUTTOM);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, essential_bc_values, P_INIT);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form_iron), HERMES_SYM, 3);
  wf.add_matrix_form(callback(bilinear_form_wire), HERMES_SYM, 2);
  wf.add_matrix_form(callback(bilinear_form_air), HERMES_SYM, 1);
  wf.add_vector_form(callback(linear_form_wire), 2);

  // Initialize coarse and reference mesh solution.
  Solution sln, ref_sln;

  // Initialize refinement selector.
  H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
  
  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;
  
  // Adaptivity loop:
  int as = 1; 
  bool done = false;
  do
  {
    info("---- Adaptivity step %d:", as);

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

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
    
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    
    if (matrix_solver == SOLVER_AZTECOO) {
      ((AztecOOSolver*) solver)->set_solver(iterative_method);
      ((AztecOOSolver*) solver)->set_precond(preconditioner);
      // Using default iteration parameters (see solver/aztecoo.h).
    }
    
    dp->assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();
    
    // Solve the linear system of the reference problem. If successful, obtain the solution.
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");
  
    // Time measurement.
    cpu_time.tick();

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver); 

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

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

    // Time measurement.
    cpu_time.tick();

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

    // If err_est too large, adapt the mesh.
    if (err_est_rel < ERR_STOP) done = true;
    else 
//.........这里部分代码省略.........