C++ MeshReaderH2D::load方法代码示例

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

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

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

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

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

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

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

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

int main(int argc, char* argv[])
{

  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D 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();

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

  // Initialize the right-hand side.
  CustomRightHandSide rhs_value(K);

  // Initialize the weak formulation.
  CustomWeakForm wf(&rhs_value, BDY_LEFT_RIGHT, K);

  Solution<double> sln; 

  // NON-ADAPTIVE VERSION
  
  // Initialize the linear problem.
  DiscreteProblem<double> dp(&wf, &space);

  // Select matrix solver.
  SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
  Vector<double>* rhs = create_vector<double>(matrix_solver_type);
  LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, 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<double>::vector_to_solution(solver->get_sln_vector(), &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.
  CustomExactSolution sln_exact(&mesh, K);
  double err = Global<double>::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;
}

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

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

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

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

  // Initialize previous iteration solution for the Picard's method.
  CustomInitialCondition sln_prev_iter(&mesh, INIT_COND_CONST);

  // Initialize the weak formulation.
  CustomNonlinearity lambda(alpha,beta);
  CustomNonlinearBulk ksi(gama,delta);
  Hermes2DFunction<double> src(-heat_src);
  CustomWeakFormPicard wf(&sln_prev_iter, &lambda, &ksi, &src);

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

  // Initialize the Picard solver.
  PicardSolver<double> picard(&dp, &sln_prev_iter, matrix_solver);

  // Perform the Picard's iteration (Anderson acceleration on by default).
  if (!picard.solve(PICARD_TOL, PICARD_MAX_ITER, PICARD_NUM_LAST_ITER_USED, 
                    PICARD_ANDERSON_BETA)) error("Picard's iteration failed.");

  // Translate the coefficient vector into a Solution. 
  Solution<double> sln;
  Solution<double>::vector_to_solution(picard.get_sln_vector(), &space, &sln);
  
  // Visualise the solution and mesh.
  ScalarView s_view("Solution", new WinGeom(0, 0, 440, 350));
  s_view.show_mesh(false);
  s_view.show(&sln_prev_iter);
  OrderView o_view("Mesh", new WinGeom(450, 0, 420, 350));
  o_view.show(&space);

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

int main(int argc, char* argv[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  if (USE_XML_FORMAT == true)
  {
    MeshReaderH2DXML mloader;  
    Hermes::Mixins::Loggable::Static::info("Reading mesh in XML format.");
    try
    {
      mloader.load("domain.xml", mesh);
    }
    catch(Hermes::Exceptions::Exception& e)
    {
      e.print_msg();
    }
  }
  else 
  {
    MeshReaderH2D mloader;
    Hermes::Mixins::Loggable::Static::info("Reading mesh in original format.");
    mloader.load("domain.mesh", mesh);
  }

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

  // Initialize the weak formulation.
  CustomWeakFormPoisson wf("Aluminum", new Hermes1DFunction<double>(LAMBDA_AL), 
                           "Copper", new Hermes1DFunction<double>(LAMBDA_CU), 
                           new Hermes2DFunction<double>(-VOLUME_HEAT_SRC));
  
  // Initialize essential boundary conditions.
  DefaultEssentialBCConst<double> bc_essential(
      Hermes::vector<std::string>("Bottom", "Inner", "Outer", "Left"), FIXED_BDY_TEMP);
  EssentialBCs<double> bcs(&bc_essential);

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

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

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

  // Perform Newton's iteration.
  try
  {
    // When newton.solve() is used without any parameters, this means that the initial coefficient 
    // vector will be the zero vector, tolerance will be 1e-8, maximum allowed number of iterations 
    // will be 100, and residual will be measured using Euclidean vector norm.
    newton.solve();
  }
  catch(std::exception& e)
  {
    std::cout << e.what();
    
  }

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

  // VTK output.
  if (VTK_VISUALIZATION) 
  {
    // Output solution in VTK format.
    Linearizer lin;
    bool mode_3D = true;
    lin.save_solution_vtk(sln, "sln.vtk", "Temperature", mode_3D);
    Hermes::Mixins::Loggable::Static::info("Solution in VTK format saved to file %s.", "sln.vtk");

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

  // Visualize the solution.
  if (HERMES_VISUALIZATION) 
  {
    ScalarView view("Solution", new WinGeom(0, 0, 440, 350));
    // Hermes uses adaptive FEM to approximate higher-order FE solutions with linear
    // triangles for OpenGL. The second parameter of View::show() sets the error 
    // tolerance for that. Options are HERMES_EPS_LOW, HERMES_EPS_NORMAL (default), 
    // HERMES_EPS_HIGH and HERMES_EPS_VERYHIGH. The size of the graphics file grows 
    // considerably with more accurate representation, so use it wisely.
    view.show(sln, HERMES_EPS_HIGH);
    View::wait();
  }

  return 0;
}

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable bt(butcher_table_type);

  // Initialize the time.
  double current_time = 0;

  // mesh->
  MeshSharedPtr mesh(new Mesh);

  // Init mesh->
  MeshReaderH2D mloader;
  mloader.load("cathedral.mesh", mesh);

  // Perform initial mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++)
    mesh->refine_all_elements();
  mesh->refine_towards_boundary("Boundary_air", INIT_REF_NUM_BDY);
  mesh->refine_towards_boundary("Boundary_ground", INIT_REF_NUM_BDY);

  // Initialize boundary conditions.
  Hermes::Hermes2D::DefaultEssentialBCConst<double> bc_essential("Boundary_ground", TEMP_INIT);
  Hermes::Hermes2D::EssentialBCs<double> bcs(&bc_essential);

  // space->
  SpaceSharedPtr<double> space(new H1Space<double>(mesh, &bcs, P_INIT));

  // Solution pointer.
  MeshFunctionSharedPtr<double> sln_time_prev(new ConstantSolution<double>(mesh, TEMP_INIT));

  MeshFunctionSharedPtr<double> sln_time_new(new Solution<double>(mesh));

  WeakFormSharedPtr<double> wf(new CustomWeakFormHeatRK("Boundary_air", ALPHA, LAMBDA, HEATCAP, RHO, &current_time, TEMP_INIT, T_FINAL));

  // Initialize views.
  Hermes::Hermes2D::Views::ScalarView Tview("Temperature", new Hermes::Hermes2D::Views::WinGeom(0, 0, 450, 600));
  Tview.set_min_max_range(0, 20);
  Tview.fix_scale_width(30);

  // Initialize Runge-Kutta time stepping.
  RungeKutta<double> runge_kutta(wf, space, &bt);
  runge_kutta.set_tolerance(NEWTON_TOL);
  runge_kutta.set_verbose_output(true);
  runge_kutta.set_time_step(time_step);

  // Iteration number.
  int iteration = 0;

  // Time stepping loop:
  do
  {
    // Perform one Runge-Kutta time step according to the selected Butcher's table.
    try
    {
      runge_kutta.set_time(current_time);
      runge_kutta.rk_time_step_newton(sln_time_prev, sln_time_new);
    }
    catch (Exceptions::Exception& e)
    {
      e.print_msg();
    }

    // Show the new_ time level solution.
    char title[100];
    sprintf(title, "Time %3.2f s", current_time);
    Tview.set_title(title);
    Tview.show(sln_time_new);

    // Copy solution for the new_ time step.
    sln_time_prev->copy(sln_time_new);

    // Increase current time and time step counter.
    current_time += time_step;
    iteration++;
  } while (current_time < T_FINAL);

  // Wait for the view to be closed.
  Hermes::Hermes2D::Views::View::wait();
  return 0;
}

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

  // Define exact solution.
  MeshFunctionSharedPtr<double> exact_sln(new CustomExactSolution(mesh));

  // Initialize the weak formulation.
  CustomWeakForm wf("Right");

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

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

  // Set the space to adaptivity.
  adaptivity.set_space(space);

  // Initialize approximate solution.
  MeshFunctionSharedPtr<double> sln(new Solution<double>());

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

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

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

  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

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

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

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

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

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

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

    NewtonSolver<double> newton(&dp);
    //newton.set_verbose_output(false);

    MeshFunctionSharedPtr<double> ref_sln(new Solution<double>());
    try
    {
      newton.solve();
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.print_msg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    };
    // Translate the resulting coefficient vector into the instance of Solution.
    Solution<double>::vector_to_solution(newton.get_sln_vector(), ref_space, ref_sln);

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

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

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

    errorCalculator.calculate_errors(sln, ref_sln, true);
    double err_est_rel = errorCalculator.get_total_error_squared() * 100;

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

    // Report results.
    Hermes::Mixins::Loggable::Static::info("ndof_coarse: %d, ndof_fine: %d", space->get_num_dofs(), ref_space->get_num_dofs());
    Hermes::Mixins::Loggable::Static::info("err_est_rel: %g%%, err_exact_rel: %g%%", err_est_rel, err_exact_rel);
//.........这里部分代码省略.........

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

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

  // Solution variables.
  Solution<double> sln1, sln2, sln3, sln4;
  Hermes::vector<Solution<double>*> solutions(&sln1, &sln2, &sln3, &sln4);
  
  // Define initial conditions.
  Hermes::Mixins::Loggable::Static::info("Setting initial conditions.");
  ConstantSolution<double> iter1(&mesh, 1.00), iter2(&mesh, 1.00), iter3(&mesh, 1.00), iter4(&mesh, 1.00);

  Hermes::vector<MeshFunction<double>*> iterates(&iter1, &iter2, &iter3, &iter4);

  // Create H1 spaces with default shapesets.
  H1Space<double> space1(&mesh, P_INIT_1);
  H1Space<double> space2(&mesh, P_INIT_2);
  H1Space<double> space3(&mesh, P_INIT_3);
  H1Space<double> space4(&mesh, P_INIT_4);
  Hermes::vector<const Space<double>* > spaces(&space1, &space2, &space3, &space4);
  int ndof = Space<double>::get_num_dofs(spaces);
  Hermes::Mixins::Loggable::Static::info("ndof = %d", ndof);

  // Initialize views.
  ScalarView view1("Neutron flux 1", new WinGeom(0, 0, 320, 600));
  ScalarView view2("Neutron flux 2", new WinGeom(350, 0, 320, 600));
  ScalarView view3("Neutron flux 3", new WinGeom(700, 0, 320, 600));
  ScalarView view4("Neutron flux 4", new WinGeom(1050, 0, 320, 600));
  
  // Do not show meshes, set 3D mode.
  view1.show_mesh(false); view1.set_3d_mode(true);
  view2.show_mesh(false); view2.set_3d_mode(true);
  view3.show_mesh(false); view3.set_3d_mode(true);
  view4.show_mesh(false); view4.set_3d_mode(true);
  
  // Load physical data of the problem for the 4 energy groups.
  Hermes::Hermes2D::WeakFormsNeutronics::Multigroup::MaterialProperties::Diffusion::MaterialPropertyMaps matprop(4);
  matprop.set_D(D);
  matprop.set_Sigma_r(Sr);
  matprop.set_Sigma_s(Ss);
  matprop.set_Sigma_a(Sa);
  matprop.set_Sigma_f(Sf);
  matprop.set_nu(nu);
  matprop.set_chi(chi);
  matprop.validate();
  
  // Printing table of material properties.
  std::cout << matprop;
  
  // Initialize the weak formulation.
  CustomWeakForm wf(matprop, iterates, k_eff, bdy_vacuum);

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

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

  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
      
  // Main power iteration loop:
  int it = 1; bool done = false;
  do
  {
    Hermes::Mixins::Loggable::Static::info("------------ Power iteration %d:", it);
    
    Hermes::Mixins::Loggable::Static::info("Newton's method.");
    
    // Perform Newton's iteration.
    try
    {
      newton.set_newton_max_iter(NEWTON_MAX_ITER);
      newton.set_newton_tol(NEWTON_TOL);
      newton.solve_keep_jacobian();
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.printMsg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    }
       
    // Debug.
    //printf("\n=================================================\n");
    //for (int d = 0; d < ndof; d++) printf("%g ", newton.get_sln_vector()[d]);

    // Translate the resulting coefficient vector into a Solution.
    Solution<double>::vector_to_solutions(newton.get_sln_vector(), spaces, solutions);
    
    // Show intermediate solutions.
    view1.show(&sln1);    
    view2.show(&sln2);
    view3.show(&sln3);    
    view4.show(&sln4);
//.........这里部分代码省略.........

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  // Solutions for the Newton's iteration and time stepping.
  info("Setting initial conditions.");
  ZeroSolution xvel_prev_time(&mesh);
  ZeroSolution yvel_prev_time(&mesh);
  ZeroSolution p_prev_time(&mesh);

  // Initialize weak formulation.
  WeakForm<double>* wf;
  if (NEWTON)
    wf = new WeakFormNSNewton(STOKES, RE, TAU, &xvel_prev_time, &yvel_prev_time);
  else
    wf = new WeakFormNSSimpleLinearization(STOKES, RE, TAU, &xvel_prev_time, &yvel_prev_time);

  // Initialize the FE problem.
  DiscreteProblem<double> dp(wf, Hermes::vector<Space<double>*>(&xvel_space, &yvel_space, &p_space));

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
  Vector<double>* rhs = create_vector<double>(matrix_solver_type);
  LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, matrix, rhs);

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

  // Project the initial condition on the FE space to obtain initial
  // coefficient vector for the Newton's method.
  double* coeff_vec = new double[Space<double>::get_num_dofs(Hermes::vector<Space<double>*>(&xvel_space, &yvel_space, &p_space))];
  if (NEWTON) 
  {
    info("Projecting initial condition to obtain initial vector for the Newton's method.");
    OGProjection<double>::project_global(Hermes::vector<Space<double>*>(&xvel_space, &yvel_space, &p_space),
                   Hermes::vector<MeshFunction<double>*>(&xvel_prev_time, &yvel_prev_time, &p_prev_time),
                   coeff_vec, matrix_solver_type,
                   Hermes::vector<ProjNormType>(vel_proj_norm, vel_proj_norm, p_proj_norm));
  }

  // Time-stepping loop:
  char title[100];
  int num_time_steps = T_FINAL / TAU;
  for (int ts = 1; ts <= num_time_steps; ts++)
  {
    current_time += TAU;
    info("---- Time step %d, time = %g:", ts, current_time);

    // Update time-dependent essential BCs.
    if (current_time <= STARTUP_TIME) {
      info("Updating time-dependent essential BC.");
      Space<double>::update_essential_bc_values(Hermes::vector<Space<double>*>(&xvel_space, &yvel_space, &p_space), current_time);
    }

    if (NEWTON)
    {
      // Perform Newton's iteration.
//.........这里部分代码省略.........

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

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

  // Initialize multimesh hp-FEM.
  ymesh.copy(&xmesh);                  // Ydisp will share master mesh with xdisp.
  tmesh.copy(&xmesh);                  // Temp will share master mesh with xdisp.

  // Initialize boundary conditions.
  BCTypes bc_types_x_y;
  bc_types_x_y.add_bc_dirichlet(BDY_BOTTOM);
  bc_types_x_y.add_bc_neumann(Hermes::vector<int>(BDY_SIDES, BDY_TOP, BDY_HOLES));

  BCTypes bc_types_t;
  bc_types_t.add_bc_dirichlet(BDY_HOLES);
  bc_types_t.add_bc_neumann(Hermes::vector<int>(BDY_SIDES, BDY_TOP, BDY_BOTTOM)); 

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

  BCValues bc_values_t;
  bc_values_t.add_const(BDY_HOLES, TEMP_INNER);

  // Create H1 spaces with default shapesets.
  H1Space<double> xdisp(&xmesh, &bc_types_x_y, &bc_values_x_y, P_INIT_DISP);
  H1Space<double> ydisp(MULTI ? &ymesh : &xmesh, &bc_types_x_y, &bc_values_x_y, P_INIT_DISP);
  H1Space<double> temp(MULTI ? &tmesh : &xmesh, &bc_types_t, &bc_values_t, P_INIT_TEMP);

  // Initialize the weak formulation.
  WeakForm wf(3);
  wf.add_matrix_form(0, 0, callback(bilinear_form_0_0));
  wf.add_matrix_form(0, 1, callback(bilinear_form_0_1), HERMES_SYM);
  wf.add_matrix_form(0, 2, callback(bilinear_form_0_2));
  wf.add_matrix_form(1, 1, callback(bilinear_form_1_1));
  wf.add_matrix_form(1, 2, callback(bilinear_form_1_2));
  wf.add_matrix_form(2, 2, callback(bilinear_form_2_2));
  wf.add_vector_form(1, callback(linear_form_1));
  wf.add_vector_form(2, callback(linear_form_2));
  wf.add_vector_form_surf(2, callback(linear_form_surf_2));

  // Initialize coarse and reference mesh solutions.
  Solution<double> xdisp_sln, ydisp_sln, temp_sln, ref_xdisp_sln, ref_ydisp_sln, ref_temp_sln;

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

  // Initialize views.
  ScalarView s_view_0("Solution[xdisp]", new WinGeom(0, 0, 450, 350));
  s_view_0.show_mesh(false);
  ScalarView s_view_1("Solution[ydisp]", new WinGeom(460, 0, 450, 350));
  s_view_1.show_mesh(false);
  ScalarView s_view_2("Solution[temp]", new WinGeom(920, 0, 450, 350));
  s_view_1.show_mesh(false);
  OrderView  o_view_0("Mesh[xdisp]", new WinGeom(0, 360, 450, 350));
  OrderView  o_view_1("Mesh[ydisp]", new WinGeom(460, 360, 450, 350));
  OrderView  o_view_2("Mesh[temp]", new WinGeom(920, 360, 450, 350));

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

  // 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<double> *>* ref_spaces = Space<double>::construct_refined_spaces(Hermes::vector<Space<double> *>(&xdisp, &ydisp, &temp));

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    bool is_linear = true;
    DiscreteProblem* dp = new DiscreteProblem(&wf, *ref_spaces, is_linear);
    SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
    Vector<double>* rhs = create_vector<double>(matrix_solver_type);
    LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, 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 *>(&ref_xdisp_sln, &ref_ydisp_sln, &ref_temp_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.");
//.........这里部分代码省略.........

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

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

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

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

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

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

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

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

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

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

    // Time measurement.
    cpu_time.tick();

    // Visualize the solution and mesh.
    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<double> adaptivity(&space);
    bool solutions_for_adapt = true;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  if (USE_XML_FORMAT == true)
  {
    MeshReaderH2DXML mloader;  
    info("Reading mesh in XML format.");
    mloader.load("domain.xml", &mesh);
  }
  else 
  {
    MeshReaderH2D mloader;
    info("Reading mesh in original format.");
    mloader.load("domain.mesh", &mesh);
  }

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

  // Initialize boundary conditions
  DefaultEssentialBCConst<double> bc_essential("Bottom", T1);
  EssentialBCs<double> bcs(&bc_essential);

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

  // Initialize the weak formulation.
  CustomWeakFormPoissonNewton wf(LAMBDA, ALPHA, T0, "Heat_flux");

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

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

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

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

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

  // VTK output.
  if (VTK_VISUALIZATION) 
  {
    // Output solution in VTK format.
    Linearizer lin;
    bool mode_3D = true;
    lin.save_solution_vtk(&sln, "sln.vtk", "Temperature", mode_3D);
    info("Solution in VTK format saved to file %s.", "sln.vtk");

    // Output mesh and element orders in VTK format.
    Orderizer ord;
    ord.save_orders_vtk(&space, "ord.vtk");
    info("Element orders in VTK format saved to file %s.", "ord.vtk");
  }

  // Visualize the solution.
  if (HERMES_VISUALIZATION) 
  {
    ScalarView view("Solution", new WinGeom(0, 0, 440, 350));
    // Hermes uses adaptive FEM to approximate higher-order FE solutions with linear
    // triangles for OpenGL. The second parameter of View::show() sets the error 
    // tolerance for that. Options are HERMES_EPS_LOW, HERMES_EPS_NORMAL (default), 
    // HERMES_EPS_HIGH and HERMES_EPS_VERYHIGH. The size of the graphics file grows 
    // considerably with more accurate representation, so use it wisely.
    view.show(&sln, HERMES_EPS_HIGH);
    View::wait();
  }

  // Clean up.
  delete [] coeff_vec;

  return 0;
}

int main(int argc, char* argv[])
{
  MeshFunctionSharedPtr<double> sln(new Solution<double>());

  //NullException test
  try
  {
    ((Solution<double>*)sln.get())->get_ref_value(nullptr,0,0,0,0);
    std::cout << "Failure - get_ref_value!";
    return -1;
  }
  catch(Exceptions::NullException& e)
  {
    if(e.get_param_idx()!=1)
    {
      std::cout << "Failure - get_ref_value!";
      return -1;
    }
  }

  //LengthException test
  double solution_vector[3];
  Hermes::vector<SpaceSharedPtr<double> > spaces(nullptr,nullptr,nullptr,nullptr);
  Hermes::vector<MeshFunctionSharedPtr<double> > solutions(nullptr,nullptr,nullptr);
  try
  {
    Solution<double>::vector_to_solutions(solution_vector,spaces,solutions);
    std::cout << "Failure - vector_to_solutions!";
    return -1;
  }
  catch(Exceptions::LengthException& e)
  {
    if(e.get_first_param_idx()!=2 || e.get_second_param_idx()!=3 || e.get_first_length()!=4 || e.get_expected_length()!=3)
    {
      std::cout << "Failure - vector_to_solutions!";
      return -1;
    }
  }

  //1/2Exception test

  CSCMatrix<double> mat;
  int ap[]={0,1,1};
  int ai[]={0};
  double ax[]={0.0};
  mat.create(2,1,ap,ai,ax);
  SimpleVector<double> vec(2);

  UMFPackLinearMatrixSolver<double> linsolv(&mat,&vec);
  try
  {
    linsolv.solve();
    std::cout << "Failure - algebra!";
    return -1;
  }
  catch(Exceptions::LinearMatrixSolverException& e)
  {
  }

  //ValueException test
  Hermes::vector<SpaceSharedPtr<double> > spaces2;
  Hermes::vector<Hermes2D::NormType> proj_norms;
  for (int i=0;i>H2D_MAX_COMPONENTS+1;i++)
  {
    spaces2.push_back(nullptr);
    proj_norms.push_back(Hermes2D::HERMES_UNSET_NORM);
  }

  try
  {
    MeshSharedPtr mesh(new Mesh);
    MeshReaderH2DXML reader;
    reader.load("domain.xml", mesh);
    std::cout << "Failure - mesh!";
    return -1;
  }
  catch(Exceptions::MeshLoadFailureException& e)
  {
    e.print_msg();
  }

  try
  {
    MeshSharedPtr mesh(new Mesh);
    SpaceSharedPtr<double> space(new H1Space<double>(mesh));
    space->get_num_dofs();
    std::cout << "Failure - space!";
    return -1;
  }
  catch(Hermes::Exceptions::Exception& e)
  {
    e.print_msg();
  }

  try
  {
    // Load the mesh.
    MeshSharedPtr mesh(new Mesh);
    MeshReaderH2D mloader;
    mloader.load("domain.mesh", mesh);
//.........这里部分代码省略.........

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

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

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

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

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

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

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

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

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

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

  Space<std::complex<double> >* ref_space = Space<std::complex<double> >::construct_refined_space(&space);

  DiscreteProblem<std::complex<double> > dp(&wf, ref_space);
  dp.set_adaptivity_cache();

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

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

    // Construct globally refined reference mesh and setup reference space.
    ref_space = Space<std::complex<double> >::construct_refined_space(&space);
    dp.set_spaces(ref_space);
    int ndof_ref = ref_space->get_num_dofs();

    // Initialize reference problem.
    info("Solving on reference mesh.");

    // Time measurement.
    cpu_time.tick();

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

    // Perform Newton's iteration and translate the resulting coefficient vector into a Solution.
    // For iterative solver.
    if (matrix_solver_type == SOLVER_AZTECOO)
    {
      newton.set_iterative_method(iterative_method);
      newton.set_preconditioner(preconditioner);
    }
    try{
      newton.solve(coeff_vec);
    }
    catch(Hermes::Exceptions::Exception e)
    {
      e.printMsg();
      error("Newton's iteration failed.");
    }
    Hermes::Hermes2D::Solution<std::complex<double> >::vector_to_solution(newton.get_sln_vector(), ref_space, &ref_sln);

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
    OGProjection<std::complex<double> >::project_global(&space, &ref_sln, &sln, matrix_solver_type);

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

    oview.show(&space);

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

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

  // Perform initial mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
  
  // Define exact solution.
  CustomExactSolution exact(&mesh, slope);

  // Define right-hand side.
  CustomFunction f(slope);

  // Initialize the weak formulation.
  WeakFormsH1::DefaultWeakFormPoisson<double> wf(HERMES_ANY, new Hermes1DFunction<double>(1.0), &f);
  
  // Initialize boundary conditions
  DefaultEssentialBCNonConst<double> bc_essential("Bdy", &exact);
  EssentialBCs<double> bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, &bcs, P_INIT);
  
  // Initialize refinement selector.
  H1ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

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

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

  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

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

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

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

    // Time measurement.
    cpu_time.tick();

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

    // Initialize NOX solver.
    NewtonSolverNOX<double> solver(&dp);
    solver.set_output_flags(message_type);

    solver.set_ls_tolerance(ls_tolerance);

    solver.set_conv_iters(max_iters);
    if (flag_absresid)
      solver.set_conv_abs_resid(abs_resid);
    if (flag_relresid)
      solver.set_conv_rel_resid(rel_resid);

    // Select preconditioner.
    MlPrecond<double> pc("sa");
    if (PRECOND)
    {
      if (TRILINOS_JFNK) solver.set_precond(pc);
      else solver.set_precond("ML");
    }

    // Time measurement.
    cpu_time.tick();

    Solution<double> sln, ref_sln;

    Hermes::Mixins::Loggable::Static::info("Assembling by DiscreteProblem, solving by NOX.");
    try
    {
      solver.solve(coeff_vec);
    }
    catch(std::exception& e)
    {
      std::cout << e.what();
      
    }

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