C++ Adapt::calc_err_exact方法代码示例

int main(int argc, char **args)
{
  // Test variable.
  int success_test = 1;

  if (argc < 2) error("Not enough parameters.");

  // Load the mesh.
	Mesh mesh;
  H3DReader mloader;
  if (!mloader.load(args[1], &mesh)) error("Loading mesh file '%s'.", args[1]);

	// Initialize the space.
#if defined NONLIN1
	Ord3 order(1, 1, 1);
#else
	Ord3 order(2, 2, 2);
#endif
	H1Space space(&mesh, bc_types, essential_bc_values, order);

#if defined NONLIN2
	// Do L2 projection of zero function.
	WeakForm proj_wf;
	proj_wf.add_matrix_form(biproj_form<double, scalar>, biproj_form<Ord, Ord>, HERMES_SYM);
	proj_wf.add_vector_form(liproj_form<double, scalar>, liproj_form<Ord, Ord>);

	bool is_linear = true;
	DiscreteProblem lp(&proj_wf, &space, is_linear);

  // 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_proj = create_linear_solver(matrix_solver, matrix, rhs);
  
  // Initialize the preconditioner in the case of SOLVER_AZTECOO.
  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver_proj)->set_solver(iterative_method);
    ((AztecOOSolver*) solver_proj)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }
  
  // Assemble the linear problem.
  info("Assembling (ndof: %d).", Space::get_num_dofs(&space));
  lp.assemble(matrix, rhs);
    
  // Solve the linear system.
  info("Solving.");
  if(!solver_proj->solve()) error ("Matrix solver failed.\n");

  delete matrix;
  delete rhs;
#endif

	// Initialize the weak formulation.
	WeakForm wf(1);
	wf.add_matrix_form(0, 0, jacobi_form<double, scalar>, jacobi_form<Ord, Ord>, HERMES_NONSYM);
	wf.add_vector_form(0, resid_form<double, scalar>, resid_form<Ord, Ord>);

	// Initialize the FE problem.
#if defined NONLIN2
  is_linear = false;
#else
  bool is_linear = false;
#endif
	DiscreteProblem dp(&wf, &space, is_linear);

	NoxSolver solver(&dp);

#if defined NONLIN2
solver.set_init_sln(solver_proj->get_solution());
delete solver_proj;
#endif

solver.set_conv_iters(10);
	info("Solving.");
	Solution sln(&mesh);
	if(solver.solve()) Solution::vector_to_solution(solver.get_solution(), &space, &sln);
  else error ("Matrix solver failed.\n");

		Solution ex_sln(&mesh);
#ifdef NONLIN1
		ex_sln.set_const(100.0);
#else
		ex_sln.set_exact(exact_solution);
#endif
		// Calculate exact error.
  info("Calculating exact error.");
  Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
  bool solutions_for_adapt = false;
  double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

  if (err_exact > EPS)
		// Calculated solution is not precise enough.
		success_test = 0;

  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
  }
//.........这里部分代码省略.........

int main(int argc, char **args)
{
    // Test variable.
    int success_test = 1;

    // Load the initial mesh.
    Mesh mesh;
    H3DReader mesh_loader;
    mesh_loader.load("../hexahedron.mesh3d", &mesh);

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

    // Create H1 space with default shapeset.
    H1Space space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));

    // Construct initial solution and set it to zero.
    Solution sln_prev(&mesh);
    sln_prev.set_zero();

    // Initialize weak formulation.
    WeakForm wf;
    wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
    wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>, HERMES_ANY_INT, &sln_prev);

    // Initialize discrete problem.
    bool is_linear = true;
    DiscreteProblem dp(&wf, &space, is_linear);

    // 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);

    // Initialize the preconditioner in the case of SOLVER_AZTECOO.
    if (matrix_solver == SOLVER_AZTECOO)
    {
        ((AztecOOSolver*) solver)->set_solver(iterative_method);
        ((AztecOOSolver*) solver)->set_precond(preconditioner);
        // Using default iteration parameters (see solver/aztecoo.h).
    }

    // Exact error for testing purposes.
    double err_exact;

    // Time stepping.
    int nsteps = (int) (FINAL_TIME/TAU + 0.5);
    for (int ts = 0; ts < nsteps;  ts++)
    {
        info("---- Time step %d, time %3.5f.", ts, TIME);

        // Assemble the linear problem.
        info("Assembling the linear problem (ndof: %d).", Space::get_num_dofs(&space));

        if (ts == 0) dp.assemble(matrix, rhs);
        else dp.assemble(NULL, rhs);

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

        // Output solution.
        if (solution_output)
            out_fn_vtk(&sln, "sln", ts);


        // Calculate exact error.
        ExactSolution esln(&mesh, fndd);

        info("Calculating exact error.");
        Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
        bool solutions_for_adapt = false;
        err_exact = adaptivity->calc_err_exact(&sln, &esln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS) * 100;
        info("Err. exact: %g%%.", err_exact);

        // Next time step.
        sln_prev = sln;
        TIME += TAU;

        // Cleanup.
        delete adaptivity;
    }

    if(err_exact > 3.00)
        success_test = 0;

    // Clean up.
    delete matrix;
    delete rhs;
    delete solver;

    if (success_test) {
        info("Success!");
        return ERR_SUCCESS;
    }
    else {
        info("Failure!");
        return ERR_FAILURE;
//.........这里部分代码省略.........

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 **args) 
{
  // Test variable.
  int success_test = 1;

	if (argc < 2) error("Not enough parameters.");

  // Load the mesh.
	Mesh mesh;
  H3DReader mloader;
  if (!mloader.load(args[1], &mesh)) error("Loading mesh file '%s'.", args[1]);

  // Initialize the space 1.
	Ord3 o1(2, 2, 2);
	H1Space space1(&mesh, bc_types, NULL, o1);

	// Initialize the space 2.
	Ord3 o2(4, 4, 4);
	H1Space space2(&mesh, bc_types, NULL, o2);

	WeakForm wf(2);
	wf.add_matrix_form(0, 0, bilinear_form_1_1<double, scalar>, bilinear_form_1_1<Ord, Ord>, HERMES_SYM);
	wf.add_matrix_form(0, 1, bilinear_form_1_2<double, scalar>, bilinear_form_1_2<Ord, Ord>, HERMES_SYM);
	wf.add_vector_form(0, linear_form_1<double, scalar>, linear_form_1<Ord, Ord>);
	wf.add_matrix_form(1, 1, bilinear_form_2_2<double, scalar>, bilinear_form_2_2<Ord, Ord>, HERMES_SYM);
	wf.add_vector_form(1, linear_form_2<double, scalar>, linear_form_2<Ord, Ord>);

  // Initialize the FE problem.
  bool is_linear = true;
  DiscreteProblem dp(&wf, Tuple<Space *>(&space1, &space2), is_linear);

  // Initialize the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  initialize_solution_environment(matrix_solver, argc, args);

  // 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);

  // Initialize the preconditioner in the case of SOLVER_AZTECOO.
  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver)->set_solver(iterative_method);
    ((AztecOOSolver*) solver)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }

  // Assemble the linear problem.
  info("Assembling (ndof: %d).", Space::get_num_dofs(Tuple<Space *>(&space1, &space2)));
  dp.assemble(matrix, rhs);

  // Solve the linear system. If successful, obtain the solution.
  info("Solving.");
	Solution sln1(&mesh);
	Solution sln2(&mesh);
  if(solver->solve()) Solution::vector_to_solutions(solver->get_solution(), Tuple<Space *>(&space1, &space2), Tuple<Solution *>(&sln1, &sln2));
  else error ("Matrix solver failed.\n");

  ExactSolution ex_sln1(&mesh, exact_sln_fn_1);
  ExactSolution ex_sln2(&mesh, exact_sln_fn_2);

  // Calculate exact error.
  info("Calculating exact error.");
  Adapt *adaptivity = new Adapt(Tuple<Space *>(&space1, &space2), Tuple<ProjNormType>(HERMES_H1_NORM, HERMES_H1_NORM));
  bool solutions_for_adapt = false;
  double err_exact = adaptivity->calc_err_exact(Tuple<Solution *>(&sln1, &sln2), Tuple<Solution *>(&ex_sln1, &ex_sln2), solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

  if (err_exact > EPS)
		// Calculated solution is not precise enough.
		success_test = 0;

  // Clean up.
  delete matrix;
  delete rhs;
  delete solver;
  delete adaptivity;

  // Properly terminate the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  finalize_solution_environment(matrix_solver);
  
  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
	}
	else {
    info("Failure!");
    return ERR_FAILURE;
	}
}

//.........这里部分代码省略.........
    bool is_linear = true;
    DiscreteProblem dp(&wf, ref_space, is_linear);

    // 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);
    
    // Initialize the preconditioner in the case of SOLVER_AZTECOO.
    if (matrix_solver == SOLVER_AZTECOO) 
    {
      ((AztecOOSolver*) solver)->set_solver(iterative_method);
      ((AztecOOSolver*) solver)->set_precond(preconditioner);
      // Using default iteration parameters (see solver/aztecoo.h).
    }
  
    // Assemble the reference problem.
    info("Assembling on reference mesh (ndof: %d).", Space::get_num_dofs(ref_space));
    dp.assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system on reference mesh. If successful, obtain the solution.
    info("Solving on reference mesh.");
    Solution ref_sln(ref_space->get_mesh());
    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 reference solution on the coarse mesh.
    Solution sln(space.get_mesh());
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver);

    // Time measurement.
    cpu_time.tick();

    // Output solution and mesh with polynomial orders.
    if (solution_output) 
    {
      out_fn_vtk(&sln, "sln", as);
      out_orders_vtk(&space, "order", as);
    }

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

    // 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) * 100;

    // Calculate exact error.
    solutions_for_adapt = false;
    double err_exact_rel = adaptivity->calc_err_exact(&sln, &exact, solutions_for_adapt) * 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);

    // Add entry to DOF and CPU convergence graphs.
    graph_dof_est.add_values(Space::get_num_dofs(&space), err_est_rel);
    graph_dof_est.save("conv_dof_est.dat");
    graph_cpu_est.add_values(cpu_time.accumulated(), err_est_rel);
    graph_cpu_est.save("conv_cpu_est.dat");
    graph_dof_exact.add_values(Space::get_num_dofs(&space), err_exact_rel);
    graph_dof_exact.save("conv_dof_exact.dat");
    graph_cpu_exact.add_values(cpu_time.accumulated(), err_exact_rel);
    graph_cpu_exact.save("conv_cpu_exact.dat");

    // If err_est_rel is too large, adapt the mesh. 
    if (err_est_rel < ERR_STOP) done = true;
    else 
    {
      info("Adapting coarse mesh.");
      adaptivity->adapt(THRESHOLD);
    }
    if (Space::get_num_dofs(&space) >= NDOF_STOP) done = true;

    // Clean up.
    delete ref_space->get_mesh();
    delete ref_space;
    delete matrix;
    delete rhs;
    delete solver;
    delete adaptivity;

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

  // Properly terminate the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  finalize_solution_environment(matrix_solver);

  return 1;
}

//.........这里部分代码省略.........
        OGProjection::project_global(spaces, Hermes::vector<MeshFunction*>((MeshFunction*)&T_fine, (MeshFunction*)&phi_fine), 
                        coeff_vec_coarse, matrix_solver_coarse);
        
        // Initialize the FE problem.
        DiscreteProblem dp_coarse(&wf, spaces, is_linear);

        // Perform Newton's iteration.
        info("Newton's solve on coarse meshes.");
        bool verbose = false;
        if (!solve_newton(coeff_vec_coarse, &dp_coarse, solver_coarse, matrix_coarse, rhs_coarse, NEWTON_TOL_COARSE, NEWTON_MAX_ITER, verbose)) 
          error("Newton's iteration failed.");
        
        // Translate the resulting coefficient vector into the actual solutions. 
        Solution::vector_to_solutions(coeff_vec_coarse, spaces, coarse_mesh_solutions);
        delete [] coeff_vec_coarse;
      } 
      else {
        // Projection onto the new coarse meshes.
        info("Projecting fine mesh solutions back onto coarse meshes.");
        OGProjection::project_global(spaces, fine_mesh_solutions, coarse_mesh_solutions, matrix_solver_coarse); 
      }

      // Calculate element errors.
      info("Calculating error estimate and exact error."); 
      Adapt* adaptivity = new Adapt(spaces);

      // Calculate error estimate for each solution component and the total error estimate.
      Hermes::vector<double> err_est_rel;
      double err_est_rel_total = adaptivity->calc_err_est(coarse_mesh_solutions, fine_mesh_solutions, &err_est_rel) * 100;

      // Calculate exact error for each solution component and the total exact error.
      bool solutions_for_adapt = false;
      Hermes::vector<double> err_exact_rel;
      double err_exact_rel_total = adaptivity->calc_err_exact(coarse_mesh_solutions, 
                                   Hermes::vector<Solution *>(&T_exact_solution, &phi_exact_solution), 
                                   &err_exact_rel, solutions_for_adapt) * 100;

      info("T: ndof_coarse: %d, ndof_fine: %d, err_est: %g %%, err_exact: %g %%", 
            space_T.get_num_dofs(), (*ref_spaces)[0]->get_num_dofs(), err_est_rel[0]*100, err_exact_rel[0]*100);
      info("phi: ndof_coarse: %d, ndof_fine: %d, err_est: %g %%, err_exact: %g %%", 
            space_phi.get_num_dofs(), (*ref_spaces)[1]->get_num_dofs(), err_est_rel[1]*100, err_exact_rel[1]*100);
      
      // If err_est too large, adapt the mesh.
      if (err_est_rel_total < ERR_STOP) done = true;
      else {
        info("Adapting the coarse meshes.");
        done = adaptivity->adapt(Hermes::vector<RefinementSelectors::Selector*> (&selector, &selector), THRESHOLD, STRATEGY, MESH_REGULARITY);
        if (Space::get_num_dofs(spaces) >= NDOF_STOP) done = true; 
      }
      
      delete adaptivity;
      delete ref_spaces;
    }
    while (!done);

    // Make the fine mesh solution at current time level the previous time level solution in the following time step.
    T_prev_time.copy(&T_fine);
    phi_prev_time.copy(&phi_fine);
  }

  info("Coordinate ( 0.0,  0.0) T value = %lf", T_prev_time.get_pt_value(0.0, 0.0));
  info("Coordinate ( 25.0, 25.0) T value = %lf", T_prev_time.get_pt_value(25.0, 25.0));
  info("Coordinate ( 25.0, 75.0) T value = %lf", T_prev_time.get_pt_value(25.0, 75.0));
  info("Coordinate ( 75.0, 25.0) T value = %lf", T_prev_time.get_pt_value(75.0, 25.0));
  info("Coordinate ( 75.0, 75.0) T value = %lf", T_prev_time.get_pt_value(75.0, 75.0));

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

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("lshape3q.mesh", &mesh);    // quadrilaterals
  //mloader.load("lshape3t.mesh", &mesh);  // triangles

  // Perform initial mesh refinemets.
  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::vector<int>(BDY_1, BDY_6));
  bc_types.add_bc_newton(Hermes::vector<int>(BDY_2, BDY_3, BDY_4, BDY_5));

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_zero(Hermes::vector<int>(BDY_1, BDY_6));

  // Create an Hcurl space with default shapeset.
  HcurlSpace 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_matrix_form_surf(callback(bilinear_form_surf));
  wf.add_vector_form_surf(linear_form_surf, linear_form_surf_ord);

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

  // Initialize exact solution.
  ExactSolution sln_exact(&mesh, exact);

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

  // Initialize views.
  VectorView v_view("Solution (magnitude)", new WinGeom(0, 0, 460, 350));
  v_view.set_min_max_range(0, 1.5);
  OrderView  o_view("Polynomial orders", new WinGeom(470, 0, 400, 350));
  
  // 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.
    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 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.
    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.
    v_view.show(&sln);
    o_view.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,
    bool solutions_for_adapt = false;
    double err_exact_rel = adaptivity->calc_err_exact(&sln, &sln_exact, solutions_for_adapt) * 100;

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

//.........这里部分代码省略.........
                                 Hermes::vector<Solution *>(&u_sln, &v_sln), matrix_solver); 
   
    // 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);

    /*
    // Exact solution for comparison with computational results.
    u_exact.update(&u_mesh, u_fndd);
    v_exact.update(&v_mesh, v_fndd);

    // Show exact solution.
    sview_u_exact.show(&u_exact);
    sprintf(title, "Exact solution for u.");
    sview_u_exact.set_title(title);
    
    sview_v_exact.show(&v_exact);
    sprintf(title, "Exact solution for v.");
    sview_v_exact.set_title(title);
    */

    // Calculate element errors.
    info("Calculating error estimate and exact error."); 
    Adapt* adaptivity = new Adapt(Hermes::vector<Space *>(&u_space, &v_space));
    
    // Calculate error estimate for each solution component and the total error estimate.
    Hermes::vector<double> err_est_rel;
    double err_est_rel_total = adaptivity->calc_err_est(Hermes::vector<Solution *>(&u_sln, &v_sln), 
                               Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln), &err_est_rel) * 100;

    // Calculate exact error for each solution component and the total exact error.
    Hermes::vector<double> err_exact_rel;
    bool solutions_for_adapt = false;
    double err_exact_rel_total = adaptivity->calc_err_exact(Hermes::vector<Solution *>(&u_sln, &v_sln), 
                                 Hermes::vector<Solution *>(&u_exact, &v_exact), 
                                 &err_exact_rel, solutions_for_adapt) * 100;

    // Time measurement.
    cpu_time.tick();

    // Report results.
    info("ndof_coarse[u]: %d, ndof_fine[u]: %d",
         u_space.Space::get_num_dofs(), Space::get_num_dofs((*ref_spaces)[0]));

    info("err_est_rel[u]: %g%%, err_exact_rel[u]: %g%%", err_est_rel[0]*100, err_exact_rel[0]*100);

    info("ndof_coarse[v]: %d, ndof_fine[v]: %d",
         v_space.Space::get_num_dofs(), Space::get_num_dofs((*ref_spaces)[1]));

    info("err_est_rel[v]: %g%%, err_exact_rel[v]: %g%%", err_est_rel[1]*100, err_exact_rel[1]*100);

    info("ndof_coarse_total: %d, ndof_fine_total: %d",
         Space::get_num_dofs(Hermes::vector<Space *>(&u_space, &v_space)), Space::get_num_dofs(*ref_spaces));

    info("err_est_rel_total: %g%%, err_est_exact_total: %g%%", err_est_rel_total, err_exact_rel_total);

    // Add entry to DOF and CPU convergence graphs.
    graph_dof_est.add_values(Space::get_num_dofs(Hermes::vector<Space *>(&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::get_num_dofs(Hermes::vector<Space *>(&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");

    // If err_est too large, adapt the mesh.
    if (err_est_rel_total < ERR_STOP) 
      done = true;
    else 
    {
      info("Adapting coarse mesh.");
      done = adaptivity->adapt(Hermes::vector<RefinementSelectors::Selector *>(&selector, &selector), 
                               THRESHOLD, STRATEGY, MESH_REGULARITY);
    }
    if (Space::get_num_dofs(Hermes::vector<Space *>(&u_space, &v_space)) >= NDOF_STOP) done = true;

    // Clean up.
    delete solver;
    delete matrix;
    delete rhs;
    delete adaptivity;
    if(done == false)
      for(unsigned int i = 0; i < ref_spaces->size(); i++)
        delete (*ref_spaces)[i]->get_mesh();
    delete ref_spaces;
    
    // Increase counter.
    as++;
  }
  while (done == false);

  verbose("Total running time: %g s", cpu_time.accumulated());

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

//.........这里部分代码省略.........
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

    // Initialize the preconditioner in the case of SOLVER_AZTECOO.
    if (matrix_solver == SOLVER_AZTECOO) 
    {
      ((AztecOOSolver*) solver)->set_solver(iterative_method);
      ((AztecOOSolver*) solver)->set_precond(preconditioner);
      // Using default iteration parameters (see solver/aztecoo.h).
    }

    // Assemble the reference problem.
    info("Assembling on reference mesh (ndof: %d).", Space::get_num_dofs(ref_space));
    lp.assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system on reference mesh. If successful, obtain the solution.
    info("Solving on reference mesh.");
    Solution ref_sln(ref_space->get_mesh());
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else {
		  printf("Matrix solver failed.\n");
		  success_test = 0;
	  }
    
    // Time measurement.
    cpu_time.tick();

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

    // Time measurement.
    cpu_time.tick();

	  // 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) * 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_est_rel);

	  // If err_est_rel is too large, adapt the mesh. 
    if (err_est_rel < ERR_STOP) {
		  done = true;
      ExactSolution ex_sln(&mesh, exact_solution);
		  
      // Calculate exact error.
      info("Calculating exact error.");
      Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
      bool solutions_for_adapt = false;
      double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

      if (err_exact > EPS)
		    // Calculated solution is not precise enough.
		    success_test = 0;
	 
      break;
	  }	
    else {
      info("Adapting coarse mesh.");
      adaptivity->adapt(THRESHOLD);
    }

    // If we reached the maximum allowed number of degrees of freedom, set the return flag to failure.
    if (Space::get_num_dofs(&space) >= NDOF_STOP)
    {
      done = true;
      success_test = 0;
    }

	  // Clean up.
    delete ref_space->get_mesh();
    delete ref_space;
    delete matrix;
    delete rhs;
    delete solver;
    delete adaptivity;

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

  // Properly terminate the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  finalize_solution_environment(matrix_solver);
  
  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
  }
  else {
    info("Failure!");
    return ERR_FAILURE;
  }
}

//.........这里部分代码省略.........
    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 solutions.
    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();

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(Hermes::vector<Space *>(&u_space, &v_space), Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln), 
                   Hermes::vector<Solution *>(&u_sln, &v_sln), matrix_solver); 
   
    // Calculate element errors.
    info("Calculating error estimate and exact error."); 
    Adapt* adaptivity = new Adapt(Hermes::vector<Space *>(&u_space, &v_space));
    
    // Calculate error estimate for each solution component and the total error estimate.
    Hermes::vector<double> err_est_rel;
    double err_est_rel_total = adaptivity->calc_err_est(Hermes::vector<Solution *>(&u_sln, &v_sln), 
                               Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln), &err_est_rel) * 100;

    // Calculate exact error for each solution component and the total exact error.
    Hermes::vector<double> err_exact_rel;
    bool solutions_for_adapt = false;
    double err_exact_rel_total = adaptivity->calc_err_exact(Hermes::vector<Solution *>(&u_sln, &v_sln), 
							    Hermes::vector<Solution *>(&u_exact, &v_exact), 
                                                            &err_exact_rel, solutions_for_adapt) * 100;

    // Time measurement.
    cpu_time.tick();

    // Report results.
    info("ndof_coarse[0]: %d, ndof_fine[0]: %d",
         u_space.Space::get_num_dofs(), (*ref_spaces)[0]->Space::get_num_dofs());
    info("err_est_rel[0]: %g%%, err_exact_rel[0]: %g%%", err_est_rel[0]*100, err_exact_rel[0]*100);
    info("ndof_coarse[1]: %d, ndof_fine[1]: %d",
         v_space.Space::get_num_dofs(), (*ref_spaces)[1]->Space::get_num_dofs());
    info("err_est_rel[1]: %g%%, err_exact_rel[1]: %g%%", err_est_rel[1]*100, err_exact_rel[1]*100);
    info("ndof_coarse_total: %d, ndof_fine_total: %d",
         Space::get_num_dofs(Hermes::vector<Space *>(&u_space, &v_space)), Space::get_num_dofs(*ref_spaces));
    info("err_est_rel_total: %g%%, err_est_exact_total: %g%%", err_est_rel_total, err_exact_rel_total);

    // Add entry to DOF and CPU convergence graphs.
    graph_dof_est.add_values(Space::get_num_dofs(Hermes::vector<Space *>(&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::get_num_dofs(Hermes::vector<Space *>(&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");

    // If err_est too large, adapt the mesh.
    if (err_est_rel_total < ERR_STOP) 
      done = true;
    else 
    {

int main(int argc, char **args) 
{
  // Test variable.
  int success_test = 1;

	for (int i = 0; i < 48; i++) {
		for (int j = 0; j < 48; j++) {
			info("Config: %d, %d ", i, j);

			Mesh mesh;

			for (unsigned int k = 0; k < countof(vtcs); k++)
				mesh.add_vertex(vtcs[k].x, vtcs[k].y, vtcs[k].z);
			unsigned int h1[] = {
					hexs[0][i][0] + 1, hexs[0][i][1] + 1, hexs[0][i][2] + 1, hexs[0][i][3] + 1,
					hexs[0][i][4] + 1, hexs[0][i][5] + 1, hexs[0][i][6] + 1, hexs[0][i][7] + 1 };
			mesh.add_hex(h1);
			unsigned int h2[] = {
					hexs[1][j][0] + 1, hexs[1][j][1] + 1, hexs[1][j][2] + 1, hexs[1][j][3] + 1,
					hexs[1][j][4] + 1, hexs[1][j][5] + 1, hexs[1][j][6] + 1, hexs[1][j][7] + 1 };
			mesh.add_hex(h2);
			// bc
			for (unsigned int k = 0; k < countof(bnd); k++) {
				unsigned int facet_idxs[Quad::NUM_VERTICES] = { bnd[k][0] + 1, bnd[k][1] + 1, bnd[k][2] + 1, bnd[k][3] + 1 };
				mesh.add_quad_boundary(facet_idxs, bnd[k][4]);
			}

			mesh.ugh();

      // Initialize the space.
			H1Space space(&mesh, bc_types, essential_bc_values);
			
#ifdef XM_YN_ZO
			Ord3 ord(4, 4, 4);
#elif defined XM_YN_ZO_2
			Ord3 ord(4, 4, 4);
#elif defined X2_Y2_Z2
			Ord3 ord(2, 2, 2);
#endif
			space.set_uniform_order(ord);

      // Initialize the weak formulation.
      WeakForm wf;
#ifdef DIRICHLET
      wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
      wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
#elif defined NEWTON
      wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
      wf.add_matrix_form_surf(bilinear_form_surf<double, scalar>, bilinear_form_surf<Ord, Ord>);
      wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
      wf.add_vector_form_surf(linear_form_surf<double, scalar>, linear_form_surf<Ord, Ord>);
#endif

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

      // 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);
      
      // Initialize the preconditioner in the case of SOLVER_AZTECOO.
      if (matrix_solver == SOLVER_AZTECOO) 
      {
        ((AztecOOSolver*) solver)->set_solver(iterative_method);
        ((AztecOOSolver*) solver)->set_precond(preconditioner);
        // Using default iteration parameters (see solver/aztecoo.h).
      }

      // Assemble the linear problem.
      info("Assembling (ndof: %d).", Space::get_num_dofs(&space));
      dp.assemble(matrix, rhs);
        
      // Solve the linear system. If successful, obtain the solution.
      info("Solving.");
      Solution sln(space.get_mesh());
      if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
      else error ("Matrix solver failed.\n");


      ExactSolution ex_sln(&mesh, exact_solution);

      // Calculate exact error.
      info("Calculating exact error.");
      Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
      bool solutions_for_adapt = false;
      double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

      if (err_exact > EPS)
      {
        // Calculated solution is not precise enough.
	      success_test = 0;
        info("failed, error:%g", err_exact);
      }
      else
        info("passed");

      // Clean up.
      delete matrix;
//.........这里部分代码省略.........

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

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

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("../lshape3q.mesh", &mesh);    // quadrilaterals
  //mloader.load("../lshape3t.mesh", &mesh);  // triangles

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

  // Initialize boundary conditions.
  DefaultEssentialBCConst bc_essential(Hermes::vector<std::string>("Corner horizontal",
                                                                   "Corner vertical"), 0);
  EssentialBCs bcs(&bc_essential);

  // Create an Hcurl space with default shapeset.
  HcurlSpace space(&mesh, &bcs, P_INIT);
  int ndof = space.get_num_dofs();
  info("ndof = %d", ndof);

  // Initialize the weak formulation.
  CustomWeakForm wf(MU_R, KAPPA);

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

  // Initialize exact solution.
  CustomExactSolution sln_exact(&mesh);

  // Initialize refinement selector.
  HcurlProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
  
  // 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.
    Space* ref_space = Space::construct_refined_space(&space);
    int ndof_ref = Space::get_num_dofs(ref_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);

    // Initialize reference problem.
    info("Solving on reference mesh.");
    DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space);

    // Time measurement.
    cpu_time.tick();

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

    // Perform Newton's iteration.
    if (!hermes2d.solve_newton(coeff_vec, dp, solver, matrix, rhs)) error("Newton's iteration failed.");

    // Translate the resulting coefficient vector into the Solution sln.
    Solution::vector_to_solution(coeff_vec, ref_space, &ref_sln);

    // 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 and exact error."); 
    Adapt* adaptivity = new Adapt(&space);
    double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln) * 100;

    // Calculate exact error,
    bool solutions_for_adapt = false;
    double err_exact_rel = adaptivity->calc_err_exact(&sln, &sln_exact, solutions_for_adapt) * 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 **args) 
{
  // Test variable.
  int success_test = 1;

  if (argc < 2) error("Not enough parameters.");

  // Load the mesh.
	Mesh mesh;
  H3DReader mloader;
  if (!mloader.load(args[1], &mesh)) error("Loading mesh file '%s'.", args[1]);

	// Initialize the space.
	int mx = 2;
	Ord3 order(mx, mx, mx);
	H1Space space(&mesh, bc_types, NULL, order);
#if defined LIN_DIRICHLET || defined NLN_DIRICHLET
	space.set_essential_bc_values(essential_bc_values);
#endif
	// Initialize the weak formulation.
	WeakForm wf;
	wf.add_vector_form(form_0<double, scalar>, form_0<Ord, Ord>);
#if defined LIN_NEUMANN || defined LIN_NEWTON
	wf.add_vector_form_surf(form_0_surf<double, scalar>, form_0_surf<Ord, Ord>);
#endif
#if defined LIN_DIRICHLET || defined NLN_DIRICHLET
	// preconditioner
	wf.add_matrix_form(precond_0_0<double, scalar>, precond_0_0<Ord, Ord>, HERMES_SYM);
#endif

	// Initialize the FE problem.
	DiscreteProblem fep(&wf, &space);

#if defined LIN_DIRICHLET || defined NLN_DIRICHLET
	// use ML preconditioner to speed-up things
	MlPrecond pc("sa");
	pc.set_param("max levels", 6);
	pc.set_param("increasing or decreasing", "decreasing");
	pc.set_param("aggregation: type", "MIS");
	pc.set_param("coarse: type", "Amesos-KLU");
#endif

	NoxSolver solver(&fep);
#if defined LIN_DIRICHLET || defined NLN_DIRICHLET
//	solver.set_precond(&pc);
#endif

	info("Solving.");
	Solution sln(&mesh);
	if(solver.solve()) Solution::vector_to_solution(solver.get_solution(), &space, &sln);
  else error ("Matrix solver failed.\n");
	

		Solution ex_sln(&mesh);
		ex_sln.set_exact(exact_solution);

		// Calculate exact error.
  info("Calculating exact error.");
  Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
  bool solutions_for_adapt = false;
  double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

  if (err_exact > EPS)
		// Calculated solution is not precise enough.
		success_test = 0;

  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
  }
  else {
    info("Failure!");
    return ERR_FAILURE;
  }
}

//.........这里部分代码省略.........
                                            Hermes::vector<Solution<double> *>(&u_ref_sln, &v_ref_sln));

    // Translate the resulting coefficient vector into the Solution<double> sln.
    Solution<double>::vector_to_solutions(coeff_vec, *ref_spaces, Hermes::vector<Solution<double> *>(&u_ref_sln, &v_ref_sln));

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection<double>::project_global(Hermes::vector<Space<double> *>(&u_space, &v_space), 
                                 Hermes::vector<Solution<double> *>(&u_ref_sln, &v_ref_sln), 
                                 Hermes::vector<Solution<double> *>(&u_sln, &v_sln), matrix_solver_type); 
   
    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.
    info("Calculating error estimate and exact error."); 
    Adapt<double>* adaptivity = new Adapt<double>(Hermes::vector<Space<double> *>(&u_space, &v_space));
    
    // Calculate error estimate for each solution component and the total error estimate.
    Hermes::vector<double> err_est_rel;
    double err_est_rel_total = adaptivity->calc_err_est(Hermes::vector<Solution<double> *>(&u_sln, &v_sln), 
                                                        Hermes::vector<Solution<double> *>(&u_ref_sln, &v_ref_sln), 
                                                        &err_est_rel) * 100;

#ifdef WITH_EXACT_SOLUTION
    // Calculate exact error for each solution component and the total exact error.
    Hermes::vector<double> err_exact_rel;
    bool solutions_for_adapt = false;
    double err_exact_rel_total = adaptivity->calc_err_exact(Hermes::vector<Solution<double> *>(&u_sln, &v_sln), 
                                                            Hermes::vector<Solution<double> *>(&exact_u, &exact_v), 
                                                            &err_exact_rel, solutions_for_adapt) * 100;
#endif

    // Time measurement.
    cpu_time.tick();

    // Report results.
#ifdef WITH_EXACT_SOLUTION
    info("ndof_coarse[0]: %d, ndof_fine[0]: %d",
         u_space.get_num_dofs(), u_ref_space->get_num_dofs());
    info("err_est_rel[0]: %g%%, err_exact_rel[0]: %g%%", err_est_rel[0]*100, err_exact_rel[0]*100);
    info("ndof_coarse[1]: %d, ndof_fine[1]: %d",
         v_space.get_num_dofs(), v_ref_space->get_num_dofs());
    info("err_est_rel[1]: %g%%, err_exact_rel[1]: %g%%", err_est_rel[1]*100, err_exact_rel[1]*100);
    info("ndof_coarse_total: %d, ndof_fine_total: %d",
         Space<double>::get_num_dofs(Hermes::vector<Space<double> *>(&u_space, &v_space)), 
         Space<double>::get_num_dofs(*ref_spaces));
    info("err_est_rel_total: %g%%, err_est_exact_total: %g%%", err_est_rel_total, err_exact_rel_total);
#else
    info("ndof_coarse[0]: %d, ndof_fine[0]: %d", u_space.get_num_dofs(), u_ref_space->get_num_dofs());
    info("err_est_rel[0]: %g%%", err_est_rel[0]*100);
    info("ndof_coarse[1]: %d, ndof_fine[1]: %d", v_space.get_num_dofs(), v_ref_space->get_num_dofs());
    info("err_est_rel[1]: %g%%", err_est_rel[1]*100);
    info("ndof_coarse_total: %d, ndof_fine_total: %d",
         Space<double>::get_num_dofs(Hermes::vector<Space<double> *>(&u_space, &v_space)), 
         Space<double>::get_num_dofs(*ref_spaces));
    info("err_est_rel_total: %g%%", err_est_rel_total);
#endif

    // Add entry to DOF and CPU convergence graphs.
    graph_dof_est.add_values(Space<double>::get_num_dofs(Hermes::vector<Space<double> *>(&u_space, &v_space)), 

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

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system on reference mesh. If successful, obtain the solution.
    info("Solving on reference mesh.");
    Solution ref_sln(ref_space->get_mesh());
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else {
		  error ("Matrix solver failed.\n");
		  success_test = 0;
	  }

    // Time measurement.
    cpu_time.tick();

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

    // Time measurement.
    cpu_time.tick();

    // Output solution and mesh with polynomial orders.
    if (solution_output) 
    {
      out_fn_vtk(&sln, "sln", as);
      out_orders_vtk(&space, "order", as);
    }

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

    // 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) * 100;

    // Calculate exact error.
    solutions_for_adapt = false;
    double err_exact_rel = adaptivity->calc_err_exact(&sln, &exact, solutions_for_adapt) * 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);

    // If err_est_rel is too large, adapt the mesh. 
    if (err_est_rel < ERR_STOP) done = true;
    else 
    {
      info("Adapting coarse mesh.");
      adaptivity->adapt(THRESHOLD);
    }
    if (Space::get_num_dofs(&space) >= NDOF_STOP) done = true;

    // Clean up.
    delete ref_space->get_mesh();
    delete ref_space;
    delete matrix;
    delete rhs;
    delete solver;
    delete adaptivity;

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

  // This is the actual test.
#define ERROR_SUCCESS                               0
#define ERROR_FAILURE                               -1
  int ndof_allowed;
  switch (ANISO_TYPE) {
  case ANISO_X: ndof_allowed = 28; break;
    case ANISO_Y: ndof_allowed = 28; break;
    case ANISO_Z: ndof_allowed = 28; break;
    case ANISO_X | ANISO_Y: ndof_allowed = 98; break;
    case ANISO_X | ANISO_Z: ndof_allowed = 98; break;
    case ANISO_Y | ANISO_Z: ndof_allowed = 98; break;
  case ANISO_X | ANISO_Y | ANISO_Z: ndof_allowed = 343; break; 
    default: error("Admissible command-line options are x, y, x, xy, xz, yz, xyz.");
  }

  int ndof = Space::get_num_dofs(&space);

  info("ndof_actual = %d", ndof);
  info("ndof_allowed = %d", ndof_allowed); 
  if (ndof > ndof_allowed)
    success_test = 0;
  
  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
  }
  else {
    info("Failure!");
    return ERR_FAILURE;
  }
}