C++ MeshReaderH2DXML类代码示例

 void Continuity<Scalar>::Record::load_mesh(Mesh* mesh)
 {
   MeshReaderH2DXML reader;
   std::stringstream filename;
   filename << Continuity<Scalar>::meshFileName << 0 << '_' << (std::string)"t = " << this->time << (std::string)"n = " << this->number << (std::string)".h2d";
   reader.load(filename.str().c_str(), mesh);
 }

void load_mesh(MeshSharedPtr& mesh, const char* filename, int num_initial_refinements)
{
  MeshReaderH2DXML mloader;
  mloader.load(filename, mesh);

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

 void CalculationContinuity<Scalar>::Record::save_mesh(Mesh* mesh)
 {
   MeshReaderH2DXML reader;
   std::stringstream filename;
   filename << CalculationContinuity<Scalar>::mesh_file_name << 0 << '_' << (std::string)"t = " << this->time << (std::string)"n = " << this->number << (std::string)".h2d";
   try
   {
     reader.save(filename.str().c_str(), mesh);
   }
   catch(std::exception& e)
   {
     throw IOCalculationContinuityException(CalculationContinuityException::meshes, IOCalculationContinuityException::output, filename.str().c_str(), e.what());
   }
 }

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh_whole_domain, mesh_bottom_left_corner, mesh_complement;
  Hermes::vector<Mesh*> meshes (&mesh_whole_domain, &mesh_bottom_left_corner, &mesh_complement);
  MeshReaderH2DXML mloader;
  mloader.load("subdomains.xml", meshes);

  // Perform initial mesh refinements (optional).
  for(int i = 0; i < INIT_REF_NUM; i++)
    for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
      meshes[meshes_i]->refine_all_elements();

  mloader.save("subdomains2.xml", meshes);
  mloader.load("subdomains2.xml", meshes);

  // Initialize essential boundary conditions.
  DefaultEssentialBCConst<double> bc_essential_whole_domain(Hermes::vector<std::string>("Bottom Left", "Bottom Right", "Top Left", "Top Right"), 0.0);
  EssentialBCs<double> bcs_whole_domain(&bc_essential_whole_domain);

  DefaultEssentialBCConst<double> bc_essential_bottom_left_corner(Hermes::vector<std::string>("Bottom Left", "Horizontal Left"), 0.0);
  EssentialBCs<double> bcs_bottom_left_corner(&bc_essential_bottom_left_corner);

  DefaultEssentialBCConst<double> bc_essential_complement(Hermes::vector<std::string>("Bottom Right", "Top Right", "Top Left", "Horizontal Left", "Vertical Bottom"), 0.0);
  EssentialBCs<double> bcs_complement(&bc_essential_complement);

  // Create H1 spaces with default shapeset.
  H1Space<double> space_whole_domain(&mesh_whole_domain, &bcs_whole_domain, P_INIT);
  int ndof_whole_domain = space_whole_domain.get_num_dofs();

  H1Space<double> space_bottom_left_corner(&mesh_bottom_left_corner, &bcs_bottom_left_corner, P_INIT);
  int ndof_bottom_left_corner = space_bottom_left_corner.get_num_dofs();

  H1Space<double> space_complement(&mesh_complement, &bcs_complement, P_INIT);
  int ndof_complement = space_complement.get_num_dofs();

  if(ndof_whole_domain == 225 && ndof_bottom_left_corner == 56 && ndof_complement == 161)
  {
    info("Success!");
    return TEST_SUCCESS;
  }
  else
  {
    info("Failure!");
    return TEST_FAILURE;
  }

  return 0;
}

int main(int argc, char* argv[])
{
  // Load the mesh.
  MeshSharedPtr mesh_whole_domain(new Mesh), mesh_bottom_left_corner(new Mesh), mesh_complement(new Mesh);
  Hermes::vector<MeshSharedPtr> meshes (mesh_bottom_left_corner, mesh_whole_domain, mesh_complement);
  MeshReaderH2DXML mloader;
  mloader.set_validation(false);
  mloader.load("subdomains.xml", meshes);

  // Perform initial mesh refinements (optional).
  for(int i = 0; i < INIT_REF_NUM; i++)
    for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
      meshes[meshes_i]->refine_all_elements();

  mloader.save("subdomains2.xml", meshes);
  mloader.load("subdomains2.xml", meshes);

  // Initialize essential boundary conditions.
  DefaultEssentialBCConst<double> bc_essential_whole_domain(Hermes::vector<std::string>("Bottom Left", "Bottom Right", "Top Left", "Top Right"), 0.0);
  EssentialBCs<double> bcs_whole_domain(&bc_essential_whole_domain);

  DefaultEssentialBCConst<double> bc_essential_bottom_left_corner(Hermes::vector<std::string>("Bottom Left", "Horizontal Left"), 0.0);
  EssentialBCs<double> bcs_bottom_left_corner(&bc_essential_bottom_left_corner);

  DefaultEssentialBCConst<double> bc_essential_complement(Hermes::vector<std::string>("Bottom Right", "Top Right", "Top Left", "Horizontal Left", "Vertical Bottom"), 0.0);
  EssentialBCs<double> bcs_complement(&bc_essential_complement);

  // Create H1 spaces with default shapeset.
  SpaceSharedPtr<double> space_whole_domain(new H1Space<double>(mesh_whole_domain, &bcs_whole_domain, P_INIT));
  int ndof_whole_domain = space_whole_domain->get_num_dofs();

  SpaceSharedPtr<double> space_bottom_left_corner(new H1Space<double>(mesh_bottom_left_corner, &bcs_bottom_left_corner, P_INIT));
  int ndof_bottom_left_corner = space_bottom_left_corner->get_num_dofs();

  SpaceSharedPtr<double> space_complement(new H1Space<double>(mesh_complement, &bcs_complement, P_INIT));
  int ndof_complement = space_complement->get_num_dofs();

  if(ndof_whole_domain == 225 && ndof_bottom_left_corner == 56 && ndof_complement == 161)
  {
    return 0;
  }
  else
  {
    return -1;
  }

  return 0;
}

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh_whole_domain, mesh_without_hole;
  Hermes::vector<Mesh*> meshes (&mesh_whole_domain, &mesh_without_hole);
  MeshReaderH2DXML mloader;
  mloader.load("subdomains.xml", meshes);

  // Perform initial mesh refinements (optional).
  for(int i = 0; i < INIT_REF_NUM; i++)
    for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
      meshes[meshes_i]->refine_all_elements();

  // Perform refinement towards the hole.
  for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
    meshes[meshes_i]->refine_towards_boundary("Inner", INIT_REF_NUM_HOLE);

  // Initialize boundary conditions.
  // Flow.
  EssentialBCNonConst bc_inlet_vel_x("Inlet", VEL_INLET, H, STARTUP_TIME);
  DefaultEssentialBCConst<double> bc_other_vel_x(Hermes::vector<std::string>("Outer", "Inner"), 0.0);
  EssentialBCs<double> bcs_vel_x(Hermes::vector<EssentialBoundaryCondition<double> *>(&bc_inlet_vel_x, &bc_other_vel_x));
  DefaultEssentialBCConst<double> bc_vel_y(Hermes::vector<std::string>("Inlet", "Outer", "Inner"), 0.0);
  EssentialBCs<double> bcs_vel_y(&bc_vel_y);
  EssentialBCs<double> bcs_pressure;

  // Temperature.
  DefaultEssentialBCConst<double> bc_temperature(Hermes::vector<std::string>("Inlet", "Outer"), 20.0);
  EssentialBCs<double> bcs_temperature(&bc_temperature);

  // Spaces for velocity components and pressure.
  H1Space<double> xvel_space(&mesh_without_hole, &bcs_vel_x, P_INIT_VEL);
  H1Space<double> yvel_space(&mesh_without_hole, &bcs_vel_y, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
  L2Space<double> p_space(&mesh_without_hole, P_INIT_PRESSURE);
#else
  H1Space<double> p_space(&mesh_without_hole, &bcs_pressure, P_INIT_PRESSURE);
#endif
  // Space<double> for temperature.
  H1Space<double> temperature_space(&mesh_whole_domain, &bcs_temperature, P_INIT_TEMP);

  // 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, &temperature_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
  ProjNormType temperature_proj_norm = HERMES_H1_NORM;

  // Solutions for the Newton's iteration and time stepping.
  info("Setting initial conditions.");
  ZeroSolution xvel_prev_time(&mesh_without_hole), yvel_prev_time(&mesh_without_hole), 
                                      p_prev_time(&mesh_without_hole);
  ConstantSolution<double>  temperature_prev_time(&mesh_whole_domain, TEMP_INIT); 

  // Calculate Reynolds number.
  double reynolds_number = VEL_INLET * OBSTACLE_DIAMETER / KINEMATIC_VISCOSITY_WATER;
  info("RE = %g", reynolds_number);

  // Initialize weak formulation.
  CustomWeakFormHeatAndFlow wf(STOKES, reynolds_number, time_step, &xvel_prev_time, &yvel_prev_time, &temperature_prev_time, 
      HEAT_SOURCE_GRAPHITE, SPECIFIC_HEAT_GRAPHITE, SPECIFIC_HEAT_WATER, RHO_GRAPHITE, RHO_WATER, 
      THERMAL_CONDUCTIVITY_GRAPHITE, THERMAL_CONDUCTIVITY_WATER);
  
  // Initialize the FE problem.
  DiscreteProblem<double> dp(&wf, Hermes::vector<Space<double> *>(&xvel_space, &yvel_space, &p_space, &temperature_space));

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

  // Initialize views.
  Views::VectorView vview("velocity [m/s]", new Views::WinGeom(0, 0, 500, 300));
  Views::ScalarView pview("pressure [Pa]", new Views::WinGeom(0, 310, 500, 300));
  Views::ScalarView tempview("temperature [C]", new Views::WinGeom(510, 0, 500, 300));
  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[ndof];
  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, &temperature_space), 
    Hermes::vector<MeshFunction<double> *>(&xvel_prev_time, &yvel_prev_time, &p_prev_time, &temperature_prev_time), 
    coeff_vec, matrix_solver_type, 
    Hermes::vector<ProjNormType>(vel_proj_norm, vel_proj_norm, p_proj_norm, temperature_proj_norm));

  // Time-stepping loop:
  char title[100];
  int num_time_steps = T_FINAL / time_step;
  double current_time = 0.0;
  for (int ts = 1; ts <= num_time_steps; ts++)
  {
//.........这里部分代码省略.........

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

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

   // Show mesh.
   MeshView mv("Mesh", new WinGeom(0, 0, 580, 400));
   mv.show(mesh);

   // Exact lambda
   MeshFunctionSharedPtr<double> exact_lambda(new ExactSolutionLambda(mesh,E,nu));
   ScalarView viewLam("lambda [Pa]", new WinGeom(0, 460, 530, 350));
   viewLam.show_mesh(false);
   viewLam.show(exact_lambda);

   // Exact lambda
   MeshFunctionSharedPtr<double> exact_mu(new ExactSolutionMu(mesh,E,nu));
   ScalarView viewMu("mu [Pa]", new WinGeom(550, 460, 530, 350));
   viewMu.show_mesh(false);
   viewMu.show(exact_mu);

   // Initialize boundary conditions.
   DefaultEssentialBCConst<double> disp_bot_top_x(Hermes::vector<std::string>("Bottom","Top"), 0.0);
   DefaultEssentialBCConst<double> disp_bot_y("Bottom", 0.0);
   DefaultEssentialBCConst<double> disp_top_y("Top", 0.1);
   EssentialBCs<double> bcs_x(&disp_bot_top_x);
   EssentialBCs<double> bcs_y(Hermes::vector<EssentialBoundaryCondition<double> *>(&disp_bot_y, &disp_top_y));

   // Create x- and y- displacement space using the default H1 shapeset.
   SpaceSharedPtr<double> u1_space(new H1Space<double>(mesh, &bcs_x, P_INIT));
   SpaceSharedPtr<double> u2_space(new H1Space<double>(mesh, &bcs_y, P_INIT));
   Hermes::vector<SpaceSharedPtr<double> > spaces(u1_space, u2_space);
   int ndof = Space<double>::get_num_dofs(spaces);
   Hermes::Mixins::Loggable::Static::info("ndof = %d", ndof);

   // Initialize the weak formulation.
   CustomWeakFormLinearElasticity wf(E, nu, &lambdaFun, &muFun, rho*g1, "Top", f0, f1);

   // Initialize Newton solver.
   NewtonSolver<double> newton(&wf, spaces);
   newton.set_verbose_output(true);

   // Perform Newton's iteration.
   try
   {
	   newton.solve();
   }
   catch(std::exception& e)
   {
     std::cout << e.what();
	 Hermes::Mixins::Loggable::Static::info("Newton's iteration failed.");
   }

   // Translate the resulting coefficient vector into the Solution sln.
   MeshFunctionSharedPtr<double> u1_sln(new Solution<double>), u2_sln(new Solution<double>);
   Hermes::vector<MeshFunctionSharedPtr<double> > solutions(u1_sln, u2_sln);

   Solution<double>::vector_to_solutions(newton.get_sln_vector(), spaces, solutions);

   // Visualize the solution.
   ScalarView view("Von Mises stress [Pa]", new WinGeom(590, 0, 700, 400));
   // First Lame constant.

   double lambda = (E * nu) / ((1 + nu) * (1 - 2*nu));
   // Second Lame constant.
   double mu = E / (2*(1 + nu));

   MeshFunctionSharedPtr<double> stress(new VonMisesFilter(solutions, lambda, mu));
   MeshFunctionSharedPtr<double> S11(new CustomFilterS11(solutions, &muFun, &lambdaFun));
   MeshFunctionSharedPtr<double> S12(new CustomFilterS12(solutions, mu));
   MeshFunctionSharedPtr<double> S22(new CustomFilterS22(solutions, mu, lambda));

   view.show_mesh(false);
   view.show(stress, HERMES_EPS_HIGH, H2D_FN_VAL_0, u1_sln, u2_sln, 1.0);

   ScalarView viewS11("S11 [Pa]", new WinGeom(0, 260, 530, 350));
   viewS11.show_mesh(false);
   viewS11.show(S11, HERMES_EPS_HIGH, H2D_FN_VAL_0, u1_sln, u2_sln, 1.0);

   ScalarView viewS12("S12 [Pa]", new WinGeom(540, 260, 530, 350));
   viewS12.show_mesh(false);
   viewS12.show(S12, HERMES_EPS_HIGH, H2D_FN_VAL_0, u1_sln, u2_sln, 1.0);

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

int main(int argc, char* argv[])
{
  // Check number of command-line parameters.
  if (argc < 2)
    error("Not enough parameters.");

  // Load the mesh.
  Mesh mesh;
  MeshReaderH2DXML mloader;
  if (strcasecmp(argv[1], "1") == 0)
    mloader.load(mesh_file_1, &mesh);
  if (strcasecmp(argv[1], "2") == 0)
    mloader.load(mesh_file_2, &mesh);
  if (strcasecmp(argv[1], "3") == 0)
    mloader.load(mesh_file_3, &mesh);

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

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

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

  // Initialize the weak formulation.
  WeakFormsH1::DefaultWeakFormPoisson<double> wf(new Hermes1DFunction<double>(1.0), new Hermes2DFunction<double>(-const_f));

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

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

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

  // Perform Newton's iteration.
  Hermes::Hermes2D::Solution<double> sln;
  Hermes::Hermes2D::NewtonSolver<double> newton(&dp, matrix_solver);
  try{
    newton.solve(coeff_vec);
  }
  catch(Hermes::Exceptions::Exception e)
  {
    e.printMsg();
    error("Newton's iteration failed.");
  }
  Hermes::Hermes2D::Solution<double>::vector_to_solution(newton.get_sln_vector(), &space, &sln);

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

  info("Coordinate ( 0.3, 0.5) value = %lf", sln.get_pt_value(0.3, 0.5));
  info("Coordinate ( 0.7, 0.5) value = %lf", sln.get_pt_value(0.7, 0.5));
  info("Coordinate ( 1.3, 0.5) value = %lf", sln.get_pt_value(1.3, 0.5));
  info("Coordinate ( 1.7, 0.5) value = %lf", sln.get_pt_value(1.7, 0.5));

  double coor_x[4] = {0.3, 0.7, 1.3, 1.7};
  double coor_y = 0.5;

  double value[4] = {0.102569, 0.167907, 0.174203, 0.109630};
  if (strcasecmp(argv[1], "2") == 0)
  {
    value[0] = 0.062896;
    value[1] = 0.096658;
    value[2] = 0.114445;
    value[3] = 0.081221;
  }

  if (strcasecmp(argv[1], "3") == 0)
  {
    value[0] = 0.048752;
    value[1] = 0.028585;
    value[2] = 0.028585;
    value[3] = 0.048752;
  }

  bool success = true;
  for (int i = 0; i < 4; i++)
  {
    if (Hermes::abs(value[i] - sln.get_pt_value(coor_x[i], coor_y)) > 1E-6)
      success = false;
  }
  if (success)
  {
    printf("Success!\n");
    return TEST_SUCCESS;
  }
  else
  {
    printf("Failure!\n");
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh_whole_domain, mesh_with_hole;
  Hermes::vector<Mesh*> meshes (&mesh_whole_domain, &mesh_with_hole);
  MeshReaderH2DXML mloader;
  mloader.load("domain.xml", meshes);

  // Temperature mesh: Initial uniform mesh refinements in graphite.
  meshes[0]->refine_by_criterion(element_in_graphite, INIT_REF_NUM_TEMPERATURE_GRAPHITE);

  // Temperature mesh: Initial uniform mesh refinements in fluid.
  meshes[0]->refine_by_criterion(element_in_fluid, INIT_REF_NUM_TEMPERATURE_FLUID);

  // Fluid mesh: Initial uniform mesh refinements.
  for(int i = 0; i < INIT_REF_NUM_FLUID; i++)
    meshes[1]->refine_all_elements();

  // Initial refinements towards boundary of graphite.
  for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
    meshes[meshes_i]->refine_towards_boundary("Inner Wall", INIT_REF_NUM_BDY_GRAPHITE);

  // Initial refinements towards the top and bottom edges.
  for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
    meshes[meshes_i]->refine_towards_boundary("Outer Wall", INIT_REF_NUM_BDY_WALL);

  /* View both meshes. */
  MeshView m1("Mesh for temperature"), m2("Mesh for fluid");
  m1.show(&mesh_whole_domain);
  m2.show(&mesh_with_hole);

  // Initialize boundary conditions.
  EssentialBCNonConst bc_inlet_vel_x("Inlet", VEL_INLET, H, STARTUP_TIME);
  DefaultEssentialBCConst<double> bc_other_vel_x(Hermes::vector<std::string>("Outer Wall", "Inner Wall"), 0.0);
  EssentialBCs<double> bcs_vel_x(Hermes::vector<EssentialBoundaryCondition<double> *>(&bc_inlet_vel_x, &bc_other_vel_x));
  DefaultEssentialBCConst<double> bc_vel_y(Hermes::vector<std::string>("Inlet", "Outer Wall", "Inner Wall"), 0.0);
  EssentialBCs<double> bcs_vel_y(&bc_vel_y);
  EssentialBCs<double> bcs_pressure;
  DefaultEssentialBCConst<double> bc_temperature(Hermes::vector<std::string>("Outer Wall", "Inlet"), 20.0);
  EssentialBCs<double> bcs_temperature(&bc_temperature);

  // Spaces for velocity components, pressure and temperature.
  H1Space<double> xvel_space(&mesh_with_hole, &bcs_vel_x, P_INIT_VEL);
  H1Space<double> yvel_space(&mesh_with_hole, &bcs_vel_y, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
  L2Space<double> p_space(&mesh_with_hole, P_INIT_PRESSURE);
#else
  H1Space<double> p_space(&mesh_with_hole, &bcs_pressure, P_INIT_PRESSURE);
#endif
  H1Space<double> temperature_space(&mesh_whole_domain, &bcs_temperature, P_INIT_TEMPERATURE);
  Hermes::vector<Space<double> *> all_spaces(&xvel_space, 
      &yvel_space, &p_space, &temperature_space);
  Hermes::vector<const Space<double> *> all_spaces_const(&xvel_space, 
      &yvel_space, &p_space, &temperature_space);

  // Calculate and report the number of degrees of freedom.
  int ndof = Space<double>::get_num_dofs(Hermes::vector<const Space<double> *>(&xvel_space, 
      &yvel_space, &p_space, &temperature_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
  ProjNormType temperature_proj_norm = HERMES_H1_NORM;
  Hermes::vector<ProjNormType> all_proj_norms = Hermes::vector<ProjNormType>(vel_proj_norm, 
      vel_proj_norm, p_proj_norm, temperature_proj_norm);

  // Initial conditions and such.
  info("Setting initial conditions.");
  ZeroSolution xvel_prev_time(&mesh_with_hole), yvel_prev_time(&mesh_with_hole), p_prev_time(&mesh_with_hole);
  CustomInitialConditionTemperature temperature_init_cond(&mesh_whole_domain, HOLE_MID_X, HOLE_MID_Y, 
      0.5*OBSTACLE_DIAMETER, TEMPERATURE_INIT_FLUID, TEMPERATURE_INIT_GRAPHITE); 
  Solution<double> temperature_prev_time;
  Hermes::vector<Solution<double> *> all_solutions = Hermes::vector<Solution<double> *>(&xvel_prev_time, 
      &yvel_prev_time, &p_prev_time, &temperature_prev_time);
  Hermes::vector<MeshFunction<double> *> all_meshfns = Hermes::vector<MeshFunction<double> *>(&xvel_prev_time, 
      &yvel_prev_time, &p_prev_time, &temperature_init_cond);

  // Project all initial conditions on their FE spaces to obtain aninitial
  // coefficient vector for the Newton's method. We use local projection
  // to avoid oscillations in temperature on the graphite-fluid interface
  // FIXME - currently the LocalProjection only does the lowest-order part (linear
  // interpolation) at the moment. Higher-order part needs to be added.
  double* coeff_vec = new double[ndof];
  info("Projecting initial condition to obtain initial vector for the Newton's method.");
  //OGProjection<double>::project_global(all_spaces, all_meshfns, coeff_vec, matrix_solver, all_proj_norms);
  LocalProjection<double>::project_local(all_spaces_const, all_meshfns, coeff_vec, matrix_solver, all_proj_norms);

  // Translate the solution vector back to Solutions. This is needed to replace
  // the discontinuous initial condition for temperature_prev_time with its projection.
  Solution<double>::vector_to_solutions(coeff_vec, all_spaces_const, all_solutions);

  // Calculate Reynolds number.
  double reynolds_number = VEL_INLET * OBSTACLE_DIAMETER / KINEMATIC_VISCOSITY_FLUID;
  info("RE = %g", reynolds_number);
  if (reynolds_number < 1e-8) error("Re == 0 will not work - the equations use 1/Re.");
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Check number of command-line parameters.
  if (argc < 2)
    throw Hermes::Exceptions::Exception("Not enough parameters.");

  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2DXML mloader;
  mloader.set_validation(false);
  if (strcasecmp(argv[1], "1") == 0)
    mloader.load(mesh_file_1, mesh);
  if (strcasecmp(argv[1], "2") == 0)
    mloader.load(mesh_file_2, mesh);
  if (strcasecmp(argv[1], "3") == 0)
    mloader.load(mesh_file_3, mesh);

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

  // Initialize boundary conditions.
  DefaultEssentialBCConst<double> bc_essential("Bdy", 0.0);
  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<double>::get_num_dofs(space);

  // Initialize the weak formulation.
  WeakFormsH1::DefaultWeakFormPoisson<double> wf(HERMES_ANY, new Hermes1DFunction<double>(1.0), new Hermes2DFunction<double>(-const_f));

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

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

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

  // Perform Newton's iteration.
  MeshFunctionSharedPtr<double> sln(new Solution<double>());
  NewtonSolver<double> newton(&dp);
  try{
    newton.solve(coeff_vec);
  }
  catch (Hermes::Exceptions::Exception& e)
  {
    e.print_msg();
  }
  Solution<double>::vector_to_solution(newton.get_sln_vector(), space, sln);

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

  double coor_x[4] = { 0.3, 0.7, 1.3, 1.7 };
  double coor_y = 0.5;

  double value[4] = { 0.102569, 0.167907, 0.174203, 0.109630 };
  if (strcasecmp(argv[1], "2") == 0)
  {
    value[0] = 0.062896;
    value[1] = 0.096658;
    value[2] = 0.114445;
    value[3] = 0.081221;
  }

  if (strcasecmp(argv[1], "3") == 0)
  {
    value[0] = 0.048752;
    value[1] = 0.028585;
    value[2] = 0.028585;
    value[3] = 0.048752;
  }

  bool success = true;
  for (int i = 0; i < 4; i++)
    success = Testing::test_value(sln->get_pt_value(coor_x[i], coor_y)->val[0], value[i], "value") && success;
  if (success)
  {
    printf("Success!\n");
    return 0;
  }
  else
  {
    printf("Failure!\n");
    return -1;
  }
}

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.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

  // Load the mesh.
  Mesh mesh;
  if (USE_XML_FORMAT == true)
  {
    MeshReaderH2DXML mloader;  
    Hermes::Mixins::Loggable::Static::info("Reading mesh in XML format.");
    mloader.load("domain.xml", &mesh);
  }
  else 
  {
    MeshReaderH2D mloader;
    Hermes::Mixins::Loggable::Static::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
  CustomEssentialBCNonConst bc_essential("Horizontal");
  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();
  Hermes::Mixins::Loggable::Static::info("ndof = %d", ndof);

  // Initialize the weak formulation.
  CustomWeakFormGeneral wf("Horizontal");

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

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

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

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

  // Time measurement.
  cpu_time.tick();

  // View the solution and mesh.
  ScalarView sview("Solution", new WinGeom(0, 0, 440, 350));
  sview.show(&sln);
  OrderView oview("Polynomial orders", new WinGeom(450, 0, 405, 350));
  oview.show(&space);

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

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

  return 0;
}

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[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2DXML mloader;
  mloader.load("domain-arcs.xml", mesh);
  
  mesh->refine_towards_boundary(BDY_SOLID_WALL_PROFILE, INIT_REF_NUM_BOUNDARY_ANISO);
  mesh->refine_towards_vertex(0, INIT_REF_NUM_VERTEX);

SpaceSharedPtr<double> space_rho(mesh, P_INIT);
  L2Space<double> space_rho_v_x(mesh, P_INIT);
  L2Space<double> space_rho_v_y(mesh, P_INIT);
  L2Space<double> space_e(mesh, P_INIT);
  L2Space<double> space_stabilization(new // Initialize boundary condition types and spaces with default shapesets.
  L2Space<double>(mesh, 0));
  int ndof = Space<double>::get_num_dofs({&space_rho, &space_rho_v_x, &space_rho_v_y, &space_e});
  Hermes::Mixins::Loggable::Static::info("ndof: %d", ndof);
        
  // Initialize solutions, set initial conditions.
  ConstantSolution<double> prev_rho(mesh, RHO_EXT);
  ConstantSolution<double> prev_rho_v_x(mesh, RHO_EXT * V1_EXT);
  ConstantSolution<double> prev_rho_v_y(mesh, RHO_EXT * V2_EXT);
  ConstantSolution<double> prev_e(mesh, QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));

  // Filters for visualization of Mach number, pressure and entropy.
  MachNumberFilter Mach_number({&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e}, KAPPA);
  PressureFilter pressure({&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e}, KAPPA);
  EntropyFilter entropy({&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e}, KAPPA, RHO_EXT, P_EXT);

  ScalarView pressure_view("Pressure", new WinGeom(0, 0, 600, 300));
  ScalarView Mach_number_view("Mach number", new WinGeom(700, 0, 600, 300));
  ScalarView entropy_production_view("Entropy estimate", new WinGeom(0, 400, 600, 300));
  ScalarView s1("prev_rho", new WinGeom(0, 0, 600, 300));
  ScalarView s2("prev_rho_v_x", new WinGeom(700, 0, 600, 300));
  ScalarView s3("prev_rho_v_y", new WinGeom(0, 400, 600, 300));
  ScalarView s4("prev_e", new WinGeom(700, 400, 600, 300));

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix<double>* matrix = create_matrix<double>();
  Vector<double>* rhs = create_vector<double>();
  Vector<double>* rhs_stabilization = create_vector<double>();
  Hermes::Solvers::LinearMatrixSolver<double>* solver = create_linear_solver<double>( matrix, rhs);

  // Set up CFL calculation class.
  CFLCalculation CFL(CFL_NUMBER, KAPPA);

  // Look for a saved solution on the disk.
  CalculationContinuity<double> continuity(CalculationContinuity<double>::onlyTime);
  int iteration = 0; double t = 0;

  if(REUSE_SOLUTION && continuity.have_record_available())
  {
    continuity.get_last_record()->load_mesh(mesh);
SpaceSharedPtr<double> *> spaceVector = continuity.get_last_record()->load_spaces(new   std::vector<Space<double>({mesh, mesh, mesh, mesh}));
    space_rho.copy(spaceVector[0], mesh);
    space_rho_v_x.copy(spaceVector[1], mesh);
    space_rho_v_y.copy(spaceVector[2], mesh);
    space_e.copy(spaceVector[3], mesh);
    continuity.get_last_record()->load_solutions({&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e},{&space_rho, &space_rho_v_x, 
      &space_rho_v_y, &space_e, &space_rho, &space_rho_v_x, &space_rho_v_y, &space_e});
    continuity.get_last_record()->load_time_step_length(time_step_n);
    t = continuity.get_last_record()->get_time();
    iteration = continuity.get_num();
  }

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

  mesh.refine_towards_boundary(BDY_SOLID_WALL_PROFILE, INIT_REF_NUM_BOUNDARY_ANISO, true, true);
  mesh.refine_towards_vertex(0, INIT_REF_NUM_VERTEX, true);

  MeshView m;
  m.show(&mesh);
  m.wait_for_close();

  // Initialize boundary condition types and spaces with default shapesets.
  L2Space<double>space_rho(&mesh, P_INIT);
  L2Space<double>space_rho_v_x(&mesh, P_INIT);
  L2Space<double>space_rho_v_y(&mesh, P_INIT);
  L2Space<double>space_e(&mesh, P_INIT);
  int ndof = Space<double>::get_num_dofs(Hermes::vector<const Space<double>*>(&space_rho, &space_rho_v_x, &space_rho_v_y, &space_e));
  info("Initial coarse ndof: %d", ndof);

  // Initialize solutions, set initial conditions.
  ConstantSolution<double> sln_rho(&mesh, RHO_EXT);
  ConstantSolution<double> sln_rho_v_x(&mesh, RHO_EXT * V1_EXT);
  ConstantSolution<double> sln_rho_v_y(&mesh, RHO_EXT * V2_EXT);
  ConstantSolution<double> sln_e(&mesh, QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));

  ConstantSolution<double> prev_rho(&mesh, RHO_EXT);
  ConstantSolution<double> prev_rho_v_x(&mesh, RHO_EXT * V1_EXT);
  ConstantSolution<double> prev_rho_v_y(&mesh, RHO_EXT * V2_EXT);
  ConstantSolution<double> prev_e(&mesh, QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));

  Solution<double> rsln_rho, rsln_rho_v_x, rsln_rho_v_y, rsln_e;

  // Numerical flux.
  VijayasundaramNumericalFlux num_flux(KAPPA);
  
  // Initialize weak formulation.
  EulerEquationsWeakFormSemiImplicitMultiComponent wf(&num_flux, KAPPA, RHO_EXT, V1_EXT, V2_EXT, P_EXT, BDY_SOLID_WALL, BDY_SOLID_WALL_PROFILE, 
    BDY_INLET, BDY_OUTLET, &prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e);

  // Filters for visualization of Mach number, pressure and entropy.
  MachNumberFilter Mach_number(Hermes::vector<MeshFunction<double>*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA);
  PressureFilter pressure(Hermes::vector<MeshFunction<double>*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA);
  EntropyFilter entropy(Hermes::vector<MeshFunction<double>*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA, RHO_EXT, P_EXT);

  ScalarView pressure_view("Pressure", new WinGeom(0, 0, 600, 300));
  ScalarView Mach_number_view("Mach number", new WinGeom(700, 0, 600, 300));
  ScalarView entropy_production_view("Entropy estimate", new WinGeom(0, 400, 600, 300));
  OrderView space_view("Space", new WinGeom(700, 400, 600, 300));
  
  // Initialize refinement selector.
  L2ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, MAX_P_ORDER);
  selector.set_error_weights(1.0, 1.0, 1.0);

  // Set up CFL calculation class.
  CFLCalculation CFL(CFL_NUMBER, KAPPA);

  // Look for a saved solution on the disk.
  Continuity<double> continuity(Continuity<double>::onlyTime);
  int iteration = 0; double t = 0;
  bool loaded_now = false;

  if(REUSE_SOLUTION && continuity.have_record_available())
  {
    continuity.get_last_record()->load_mesh(&mesh);
    continuity.get_last_record()->load_spaces(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
      &space_rho_v_y, &space_e), Hermes::vector<SpaceType>(HERMES_L2_SPACE, HERMES_L2_SPACE, HERMES_L2_SPACE, HERMES_L2_SPACE), Hermes::vector<Mesh *>(&mesh, &mesh, 
      &mesh, &mesh));
    continuity.get_last_record()->load_time_step_length(time_step);
    t = continuity.get_last_record()->get_time() + time_step;
    iteration = continuity.get_num() * EVERY_NTH_STEP + 1;
    loaded_now = true;
  }

  // Time stepping loop.
  for(; t < 5.0; t += time_step)
  {
    CFL.set_number(CFL_NUMBER + (t/5.0) * 10.0);
    info("---- Time step %d, time %3.5f.", iteration++, t);

    // Periodic global derefinements.
    if (iteration > 1 && iteration % UNREF_FREQ == 0 && REFINEMENT_COUNT > 0) 
    {
      info("Global mesh derefinement.");
      REFINEMENT_COUNT = 0;
      
      space_rho.unrefine_all_mesh_elements(true);
      
      space_rho.adjust_element_order(-1, P_INIT);
      space_rho_v_x.copy_orders(&space_rho);
      space_rho_v_y.copy_orders(&space_rho);
      space_e.copy_orders(&space_rho);
    }

    // Adaptivity loop:
    int as = 1; 
    int ndofs_prev = 0;
    bool done = false;
//.........这里部分代码省略.........