C++ H2DReader类代码示例

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

  // initialize the shapeset and the cache
  H1Shapeset shapeset;
  PrecalcShapeset pss(&shapeset);

  // create the x displacement space
  H1Space xdisp(&mesh, &shapeset);
  xdisp.set_bc_types(bc_types);
  xdisp.set_bc_values(bc_values);

  // create the y displacement space
  H1Space ydisp(&mesh, &shapeset);
  ydisp.set_bc_types(bc_types);
  ydisp.set_bc_values(bc_values);

  // initialize the weak formulation
  WeakForm wf(2);
  wf.add_biform(0, 0, callback(bilinear_form_0_0), SYM);  // Note that only one symmetric part is
  wf.add_biform(0, 1, callback(bilinear_form_0_1), SYM);  // added in the case of symmetric bilinear
  wf.add_biform(1, 1, callback(bilinear_form_1_1), SYM);  // forms.
  wf.add_liform_surf(0, callback(linear_form_surf_0), 3);
  wf.add_liform_surf(1, callback(linear_form_surf_1), 3);

  // initialize the linear system and solver
  UmfpackSolver umfpack;
  LinSystem sys(&wf, &umfpack);
  sys.set_spaces(2, &xdisp, &ydisp);
  sys.set_pss(1, &pss);

  // testing n_dof and correctness of solution vector
  // for p_init = 1, 2, ..., 10
  int success = 1;
  for (int p_init = 1; p_init <= 10; p_init++) {
    printf("********* p_init = %d *********\n", p_init);
    xdisp.set_uniform_order(p_init);
    int ndofs = xdisp.assign_dofs(0);
    ydisp.set_uniform_order(p_init);
    ndofs += ydisp.assign_dofs(ndofs);

    // assemble the stiffness matrix and solve the system
    Solution xsln, ysln;
    sys.assemble();
    sys.solve(2, &xsln, &ysln);

    scalar *sol_vector;
    int n_dof;
    sys.get_solution_vector(sol_vector, n_dof);
    printf("n_dof = %d\n", n_dof);
    double sum = 0;
    for (int i=0; i < n_dof; i++) sum += sol_vector[i];
    printf("coefficient sum = %g\n", sum);

    // Actual test. The values of 'sum' depend on the
    // current shapeset. If you change the shapeset,
    // you need to correct these numbers.
    if (p_init == 1 && fabs(sum - 3.50185e-06) > 1e-3) success = 0;
    if (p_init == 2 && fabs(sum - 4.34916e-06) > 1e-3) success = 0;
    if (p_init == 3 && fabs(sum - 4.60553e-06) > 1e-3) success = 0;
    if (p_init == 4 && fabs(sum - 4.65616e-06) > 1e-3) success = 0;
    if (p_init == 5 && fabs(sum - 4.62893e-06) > 1e-3) success = 0;
    if (p_init == 6 && fabs(sum - 4.64336e-06) > 1e-3) success = 0;
    if (p_init == 7 && fabs(sum - 4.63724e-06) > 1e-3) success = 0;
    if (p_init == 8 && fabs(sum - 4.64491e-06) > 1e-3) success = 0;
    if (p_init == 9 && fabs(sum - 4.64582e-06) > 1e-3) success = 0;
    if (p_init == 10 && fabs(sum - 4.65028e-06) > 1e-3) success = 0;
  }

#define ERROR_SUCCESS                               0
#define ERROR_FAILURE                               -1
  if (success == 1) {
    printf("Success!\n");
    return ERROR_SUCCESS;
  }
  else {
    printf("Failure!\n");
    return ERROR_FAILURE;
  }
}

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

  // Perform initial mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++) basemesh.refine_all_elements();
  basemesh.refine_by_criterion(criterion, 4);
  mesh.copy(&basemesh);

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_neumann(Hermes::Tuple<int>(BDY_SOLID_WALL, BDY_INLET_OUTLET));

  // Create L2 spaces with default shapesets.
  L2Space space_rho(&mesh, &bc_types, P_INIT);
  L2Space space_rho_v_x(&mesh, &bc_types, P_INIT);
  L2Space space_rho_v_y(&mesh, &bc_types, P_INIT);
  L2Space space_e(&mesh, &bc_types, P_INIT);

  // Initialize solutions, set initial conditions.
  Solution sln_rho, sln_rho_v_x, sln_rho_v_y, sln_e, prev_rho, prev_rho_v_x, prev_rho_v_y, prev_e;
  Solution rsln_rho, rsln_rho_v_x, rsln_rho_v_y, rsln_e;
  sln_rho.set_exact(&mesh, ic_density);
  sln_rho_v_x.set_exact(&mesh, ic_density_vel_x);
  sln_rho_v_y.set_exact(&mesh, ic_density_vel_y);
  sln_e.set_exact(&mesh, ic_energy);
  prev_rho.set_exact(&mesh, ic_density);
  prev_rho_v_x.set_exact(&mesh, ic_density_vel_x);
  prev_rho_v_y.set_exact(&mesh, ic_density_vel_y);
  prev_e.set_exact(&mesh, ic_energy);

  // Initialize weak formulation.
  WeakForm wf(4);

  // Bilinear forms coming from time discretization by explicit Euler's method.
  wf.add_matrix_form(0,0,callback(bilinear_form_0_0_time));
  wf.add_matrix_form(1,1,callback(bilinear_form_1_1_time));
  wf.add_matrix_form(2,2,callback(bilinear_form_2_2_time));
  wf.add_matrix_form(3,3,callback(bilinear_form_3_3_time));

  // Volumetric linear forms.
  // Linear forms coming from the linearization by taking the Eulerian fluxes' Jacobian matrices from the previous time step.
  // First flux.
  /*
  wf.add_vector_form(0,callback(linear_form_0_1), HERMES_ANY, Hermes::Tuple<MeshFunction*>(&prev_rho_v_x));
  wf.add_vector_form(1, callback(linear_form_1_0_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_1_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_2_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_3_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(2, callback(linear_form_2_0_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(2, callback(linear_form_2_1_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(2, callback(linear_form_2_2_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(2, callback(linear_form_2_3_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_0_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_1_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_2_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_3_first_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  // Second flux.
  
  wf.add_vector_form(0,callback(linear_form_0_2),HERMES_ANY, Hermes::Tuple<MeshFunction*>(&prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_0_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_1_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_2_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(1, callback(linear_form_1_3_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(2, callback(linear_form_2_0_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(2, callback(linear_form_2_1_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(2, callback(linear_form_2_2_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y));
  wf.add_vector_form(2, callback(linear_form_2_3_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_0_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_1_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_2_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  wf.add_vector_form(3, callback(linear_form_3_3_second_flux), HERMES_ANY, 
                     Hermes::Tuple<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
  */
//.........这里部分代码省略.........

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

  // a-priori mesh refinements
  mesh.refine_all_elements();
  mesh.refine_towards_boundary(5, 4, false);
  mesh.refine_towards_boundary(1, 4);
  mesh.refine_towards_boundary(3, 4);

  // initialize shapesets and the cache
  H1ShapesetBeuchler h1_shapeset;
  PrecalcShapeset h1_pss(&h1_shapeset);
#ifdef PRESSURE_IN_L2
  L2Shapeset l2_shapeset;
  PrecalcShapeset l2_pss(&l2_shapeset);
#endif

  // spaces for velocities and pressure
  H1Space xvel_space(&mesh, &h1_shapeset);
  H1Space yvel_space(&mesh, &h1_shapeset);
#ifdef PRESSURE_IN_L2
  L2Space p_space(&mesh, &l2_shapeset);
#else
  H1Space p_space(&mesh, &h1_shapeset);
#endif

  // initialize boundary conditions
  xvel_space.set_bc_types(xvel_bc_type);
  xvel_space.set_bc_values(xvel_bc_value);
  yvel_space.set_bc_types(yvel_bc_type);
  p_space.set_bc_types(p_bc_type);

  // set velocity and pressure polynomial degrees
  xvel_space.set_uniform_order(P_INIT_VEL);
  yvel_space.set_uniform_order(P_INIT_VEL);
  p_space.set_uniform_order(P_INIT_PRESSURE);

  // assign degrees of freedom
  int ndofs = 0;
  ndofs += xvel_space.assign_dofs(ndofs);
  ndofs += yvel_space.assign_dofs(ndofs);
  ndofs += p_space.assign_dofs(ndofs);

  // solutions for the Newton's iteration and time stepping
  Solution xvel_prev_time, yvel_prev_time, xvel_prev_newton, yvel_prev_newton, p_prev;
  xvel_prev_time.set_zero(&mesh);
  yvel_prev_time.set_zero(&mesh);
  xvel_prev_newton.set_zero(&mesh);
  yvel_prev_newton.set_zero(&mesh);
  p_prev.set_zero(&mesh);

  // set up weak formulation
  WeakForm wf(3);
  if (NEWTON) {
    wf.add_biform(0, 0, callback(bilinear_form_sym_0_0_1_1), SYM);
    wf.add_biform(0, 0, callback(newton_bilinear_form_unsym_0_0), UNSYM, ANY, 2, &xvel_prev_newton, &yvel_prev_newton);
    wf.add_biform(0, 1, callback(newton_bilinear_form_unsym_0_1), UNSYM, ANY, 1, &xvel_prev_newton);
    wf.add_biform(0, 2, callback(bilinear_form_unsym_0_2), ANTISYM);
    wf.add_biform(1, 0, callback(newton_bilinear_form_unsym_1_0), UNSYM, ANY, 1, &yvel_prev_newton);
    wf.add_biform(1, 1, callback(bilinear_form_sym_0_0_1_1), SYM);
    wf.add_biform(1, 1, callback(newton_bilinear_form_unsym_1_1), UNSYM, ANY, 2, &xvel_prev_newton, &yvel_prev_newton);
    wf.add_biform(1, 2, callback(bilinear_form_unsym_1_2), ANTISYM);
    wf.add_liform(0, callback(newton_F_0), ANY, 5, &xvel_prev_time, &yvel_prev_time, &xvel_prev_newton, &yvel_prev_newton, &p_prev);
    wf.add_liform(1, callback(newton_F_1), ANY, 5, &xvel_prev_time, &yvel_prev_time, &xvel_prev_newton, &yvel_prev_newton, &p_prev);
    wf.add_liform(2, callback(newton_F_2), ANY, 2, &xvel_prev_newton, &yvel_prev_newton);
  }
  else {
    wf.add_biform(0, 0, callback(bilinear_form_sym_0_0_1_1), SYM);
    wf.add_biform(0, 0, callback(simple_bilinear_form_unsym_0_0_1_1), UNSYM, ANY, 2, &xvel_prev_time, &yvel_prev_time);
    wf.add_biform(1, 1, callback(bilinear_form_sym_0_0_1_1), SYM);
    wf.add_biform(1, 1, callback(simple_bilinear_form_unsym_0_0_1_1), UNSYM, ANY, 2, &xvel_prev_time, &yvel_prev_time);
    wf.add_biform(0, 2, callback(bilinear_form_unsym_0_2), ANTISYM);
    wf.add_biform(1, 2, callback(bilinear_form_unsym_1_2), ANTISYM);
    wf.add_liform(0, callback(simple_linear_form), ANY, 1, &xvel_prev_time);
    wf.add_liform(1, callback(simple_linear_form), ANY, 1, &yvel_prev_time);
  }

  // visualization
  VectorView vview("velocity [m/s]", 0, 0, 1500, 470);
  ScalarView pview("pressure [Pa]", 0, 530, 1500, 470);
  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);

  // matrix solver
  UmfpackSolver umfpack;

  // linear system
  LinSystem ls(&wf, &umfpack);

  // nonlinear system
  NonlinSystem nls(&wf, &umfpack);

  if (NEWTON) {
//.........这里部分代码省略.........

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

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

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

  // Initialize solutions.
  CustomInitialConditionWave E_sln(&mesh);
  Solution B_sln(&mesh, 0.0);
  Hermes::vector<Solution*> slns(&E_sln, &B_sln);

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

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

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

  // Initialize Runge-Kutta time stepping.
  RungeKutta runge_kutta(&dp, &bt, matrix_solver);

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

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

  double coord_x[4] = {0.3, 0.6, 0.9, 1.4};
  double coord_y[4] = {0, 0.3, 0.5, 0.7};

  info("Coordinate (0.3, 0.0) value = %lf", B_sln.get_pt_value(coord_x[0], coord_y[0]));
  info("Coordinate (0.6, 0.3) value = %lf", B_sln.get_pt_value(coord_x[1], coord_y[1]));
  info("Coordinate (0.9, 0.5) value = %lf", B_sln.get_pt_value(coord_x[2], coord_y[2]));
  info("Coordinate (1.4, 0.7) value = %lf", B_sln.get_pt_value(coord_x[3], coord_y[3]));

  double t_value[4] = {0.0, -0.065100, -0.146515, -0.247677};
  bool success = true;

  for (int i = 0; i < 4; i++)
  {
    if (fabs(t_value[i] - B_sln.get_pt_value(coord_x[i], coord_y[i])) > 1E-6) success = false;
  }

  if (success) {  
    printf("Success!\n");
    return ERR_SUCCESS;
  }
  else {
    printf("Failure!\n");
    return ERR_FAILURE;
  }
}

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

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

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

  // Initialize the weak formulation
  CustomNonlinearity lambda(alpha);
  HermesFunction src(-heat_src);
  WeakFormsH1::DefaultWeakFormPoisson wf(HERMES_ANY, &lambda, &src);

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

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

  // Project the initial condition on the FE space to obtain initial 
  // coefficient vector for the Newton's method.
  // NOTE: If you want to start from the zero vector, just define 
  // coeff_vec to be a vector of ndof zeros (no projection is needed).
  info("Projecting to obtain initial vector for the Newton's method.");
  scalar* coeff_vec = new scalar[Space::get_num_dofs(&space)] ;
  CustomInitialCondition init_sln(&mesh);
  OGProjection::project_global(&space, &init_sln, coeff_vec, matrix_solver); 

  // Perform Newton's iteration.
  bool verbose = true;
  bool jacobian_changed = true;
  if (!hermes2d.solve_newton(coeff_vec, &dp, solver, matrix, rhs, jacobian_changed,
      NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");

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

  // Clean up.
  delete [] coeff_vec;
  delete matrix;
  delete rhs;
  delete solver;

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

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

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

    // Perform initial mesh refinements.
    for (int i = 0; i < INIT_REF_NUM; i++)
        mesh.refine_all_elements(0, true);
    mesh.refine_towards_boundary(BDY_SOLID_WALL_BOTTOM, INIT_REF_NUM_BOUNDARY_ANISO, true, false, true);
    mesh.refine_towards_boundary(BDY_SOLID_WALL_BOTTOM, INIT_REF_NUM_BOUNDARY_ISO, false, false, true);

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

    // Initialize solutions, set initial conditions.
    InitialSolutionEulerDensity sln_rho(&mesh, RHO_EXT);
    InitialSolutionEulerDensityVelX sln_rho_v_x(&mesh, RHO_EXT * V1_EXT);
    InitialSolutionEulerDensityVelY sln_rho_v_y(&mesh, RHO_EXT * V2_EXT);
    InitialSolutionEulerDensityEnergy sln_e(&mesh, QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));

    InitialSolutionEulerDensity prev_rho(&mesh, RHO_EXT);
    InitialSolutionEulerDensityVelX prev_rho_v_x(&mesh, RHO_EXT * V1_EXT);
    InitialSolutionEulerDensityVelY prev_rho_v_y(&mesh, RHO_EXT * V2_EXT);
    InitialSolutionEulerDensityEnergy prev_e(&mesh, QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));

    Solution rsln_rho, rsln_rho_v_x, rsln_rho_v_y, rsln_e;

    // Numerical flux.
    OsherSolomonNumericalFlux num_flux(KAPPA);

    // Initialize weak formulation.
    EulerEquationsWeakFormImplicitMultiComponent wf(&num_flux, KAPPA, RHO_EXT, V1_EXT, V2_EXT, P_EXT, BDY_SOLID_WALL_BOTTOM, BDY_SOLID_WALL_TOP,
            BDY_INLET, BDY_OUTLET, &prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e, PRECONDITIONING);

    wf.set_time_step(time_step);

    // Filters for visualization of Mach number, pressure and entropy.
    MachNumberFilter Mach_number(Hermes::vector<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA);
    PressureFilter pressure(Hermes::vector<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA);
    EntropyFilter entropy(Hermes::vector<MeshFunction*>(&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("1", new WinGeom(0, 0, 600, 300));
    ScalarView s2("2", new WinGeom(700, 0, 600, 300));
    ScalarView s3("3", new WinGeom(0, 400, 600, 300));
    ScalarView s4("4", new WinGeom(700, 400, 600, 300));
    */

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

    // Select preconditioner.
    RCP<Precond> pc = rcp(new IfpackPrecond("point-relax"));

    int iteration = 0;
    double t = 0;
    for(t = 0.0; t < 3.0; t += time_step) {
        info("---- Time step %d, time %3.5f.", iteration++, t);

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

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

            // Construct globally refined reference mesh and setup reference space.
            // Global polynomial order increase;
            int order_increase = 1;
            Hermes::vector<Space *>* ref_spaces = Space::construct_refined_spaces(Hermes::vector<Space *>(&space_rho, &space_rho_v_x,
                                                  &space_rho_v_y, &space_e), order_increase);

            // Report NDOFs.
            info("ndof_coarse: %d, ndof_fine: %d.",
                 Space::get_num_dofs(Hermes::vector<Space *>(&space_rho, &space_rho_v_x,
                                     &space_rho_v_y, &space_e)), Space::get_num_dofs(*ref_spaces));

            // Project the previous time level solution onto the new fine mesh
            // in order to obtain initial vector for NOX.
//.........这里部分代码省略.........

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

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

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

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

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_zero(Hermes::vector<int>(BDY_BOTTOM, BDY_RIGHT, BDY_TOP, BDY_LEFT));

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

  info("Assembling by DiscreteProblem, solving by Umfpack:");

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

  // Initialize weak formulation,
  WeakForm wf1;
  wf1.add_matrix_form(callback(jacobian_form_hermes), HERMES_NONSYM, HERMES_ANY);
  wf1.add_vector_form(callback(residual_form_hermes), HERMES_ANY);

  // Initialize the FE problem.
  bool is_linear = false;
  DiscreteProblem dp1(&wf1, &space, is_linear);

  // Set up the solver, matrix, and rhs for the coarse mesh 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 solution.
  Solution sln_hermes;

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

  // Project the initial condition on the FE space to obtain initial
  // coefficient vector for the Newton's method.
  info("Projecting to obtain initial vector for the Newton's method.");
  scalar* coeff_vec = new scalar[Space::get_num_dofs(&space)] ;
  Solution* sln_tmp = new Solution(&mesh, init_cond);
  OGProjection::project_global(&space, sln_tmp, coeff_vec, matrix_solver);
  delete sln_tmp;

  // Perform Newton's iteration.
  int it = 1;
  while (1)
  {
    // Obtain the number of degrees of freedom.
    int ndof = Space::get_num_dofs(&space);

    // Assemble the Jacobian matrix and residual vector.
    dp1.assemble(coeff_vec, matrix, rhs, false);

    // Multiply the residual vector with -1 since the matrix 
    // equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
    rhs->change_sign();
    
    // Calculate the l2-norm of residual vector.
    double res_l2_norm = get_l2_norm(rhs);

    // Info for user.
    info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(&space), res_l2_norm);

    // If l2 norm of the residual vector is within tolerance, or the maximum number 
    // of iteration has been reached, then quit.
    if (res_l2_norm < NEWTON_TOL || it > NEWTON_MAX_ITER) break;

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

    // Add \deltaY^{n+1} to Y^n.
    for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];
    
    if (it >= NEWTON_MAX_ITER)
      error ("Newton method did not converge.");

    it++;
  }
//.........这里部分代码省略.........

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

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

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

  // Initialize solutions.
  CustomInitialConditionWave E_sln(&mesh);
  Solution F_sln(&mesh, 0.0, 0.0);
  Hermes::vector<Solution*> slns(&E_sln, &F_sln);

  // Initialize the weak formulation.
  CustomWeakFormWave wf(C_SQUARED);
  
  // Initialize boundary conditions
  DefaultEssentialBCConst bc_essential(BDY, 0.0);
  EssentialBCs bcs(&bc_essential);

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

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

  // Initialize views.
  ScalarView E1_view("Solution E1", new WinGeom(0, 0, 400, 350));
  E1_view.fix_scale_width(50);
  ScalarView E2_view("Solution E2", new WinGeom(410, 0, 400, 350));
  E2_view.fix_scale_width(50);

  // Initialize Runge-Kutta time stepping.
  RungeKutta runge_kutta(&dp, &bt, matrix_solver);

  // Time stepping loop.
  double current_time = 0; int ts = 1;
  do
  {
    // Perform one Runge-Kutta time step according to the selected Butcher's table.
    info("Runge-Kutta time step (t = %g s, time_step = %g s, stages: %d).", 
         current_time, time_step, bt.get_size());
    bool verbose = true;
    bool jacobian_changed = true;
    if (!runge_kutta.rk_time_step(current_time, time_step, slns, slns, jacobian_changed, verbose))
      error("Runge-Kutta time step failed, try to decrease time step size.");

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

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

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

  return 0;
}

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

  // Initialize boundary conditions.
  DefaultEssentialBCConst essential_bc(BDY_BOTTOM, 0.0);
  EssentialBCs bcs(&essential_bc);

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

  // Initialize the weak formulation.
  CustomWeakForm wf(A_SE, A_NE, A_SW, A_NW, RHS);
  
  // Initialize refinement selector.
  H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

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

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

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

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

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

    // Time measurement.
    cpu_time.tick();

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

    // Time measurement.
    cpu_time.tick();

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

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

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

    // Time measurement.
    cpu_time.tick();

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

    // If err_est too large, adapt the mesh.
    if (err_est_rel < ERR_STOP) done = true;
    else
    {
      info("Adapting coarse mesh.");
      done = adaptivity->adapt(&selector, THRESHOLD, STRATEGY, MESH_REGULARITY);

      // Increase the counter of performed adaptivity steps.
      if (done == false)  as++;
    }
    if (Space::get_num_dofs(&space) >= NDOF_STOP) done = true;

    // Clean up.
    delete solver;
    delete matrix;
    delete rhs;
    delete adaptivity;
    if(done == false) delete ref_space->get_mesh();
    delete ref_space;
    delete dp;

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

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

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

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

  // Solution for the previous time level.
  Solution u_prev_time;

  // Initialize the weak formulation.
  WeakForm wf;
  if (TIME_INTEGRATION == 1) {
    wf.add_matrix_form(jac_euler, jac_ord, HERMES_UNSYM, HERMES_ANY, &u_prev_time);
    wf.add_vector_form(res_euler, res_ord, HERMES_ANY, &u_prev_time);
  }
  else {
    wf.add_matrix_form(jac_cranic, jac_ord, HERMES_UNSYM, HERMES_ANY, &u_prev_time);
    wf.add_vector_form(res_cranic, res_ord, HERMES_ANY, &u_prev_time);
  }

  // Project the function init_cond() on the FE space
  // to obtain initial coefficient vector for the Newton's method.
  scalar* coeff_vec = new scalar[Space::get_num_dofs(&space)];
  info("Projecting initial condition to obtain initial vector for the Newton's method.");
  OGProjection::project_global(&space, init_cond, coeff_vec, matrix_solver);
  Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);
  
  // Time stepping loop:
  double current_time = 0.0;
  int t_step = 1;
  do {
    info("---- Time step %d, t = %g s.", t_step, current_time); t_step++;

    // Initialize the FE problem.
    bool is_linear = false;
    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);

    // Perform Newton's iteration.
    int it = 1;
    while (1)
    {
      // Obtain the number of degrees of freedom.
      int ndof = Space::get_num_dofs(&space);

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

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

      // Info for user.
      info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(&space), res_l2_norm);

      // If l2 norm of the residual vector is within tolerance, or the maximum number 
      // of iteration has been reached, then quit.
      if (res_l2_norm < NEWTON_TOL || it > NEWTON_MAX_ITER) break;

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

      // Add \deltaY^{n+1} to Y^n.
      for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];
      
      if (it >= NEWTON_MAX_ITER)
        error ("Newton method did not converge.");

      it++;
    }

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

    // Cleanup.
    delete matrix;
    delete rhs;
    delete solver;

    // Update time.
    current_time += TAU;

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

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

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

  // Initial mesh refinements.
  for (int i=0; i < 4; i++) mesh.refine_all_elements();
  mesh.refine_towards_boundary(HERMES_ANY, 4);

  // Initialize boundary conditions.
  DefaultEssentialBCConst zero_vel_bc_x_brl(Hermes::vector<std::string>("Bottom", "Right", "Left"), 0.0);
  DefaultEssentialBCConst vel_bc_x_top(Hermes::vector<std::string>("Top"), XVEL_TOP);
  EssentialBCs bcs_vel_x(Hermes::vector<EssentialBoundaryCondition*>(&vel_bc_x_top, &zero_vel_bc_x_brl));
  DefaultEssentialBCConst zero_vel_bc_y(Hermes::vector<std::string>("Bottom", "Right", "Top", "Left"), 0.0);
  EssentialBCs bcs_vel_y(&zero_vel_bc_y);
  EssentialBCs bcs_pressure;

  // Spaces for velocity components and pressure.
  H1Space xvel_space(&mesh, &bcs_vel_x, P_INIT_VEL);
  H1Space yvel_space(&mesh, &bcs_vel_y, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
  L2Space p_space(&mesh, &bcs_pressure, P_INIT_PRESSURE);
#else
  H1Space p_space(&mesh, &bcs_pressure, P_INIT_PRESSURE);
#endif
  Hermes::vector<Space *> spaces = Hermes::vector<Space *>(&xvel_space, &yvel_space, &p_space);

  // Calculate and report the number of degrees of freedom.
  int ndof = Space::get_num_dofs(Hermes::vector<Space *>(&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
  ProjNormType t_proj_norm = HERMES_H1_NORM;

  // Solutions for the Newton's iteration and time stepping.
  info("Setting initial conditions.");
  Solution xvel_prev_time, yvel_prev_time, p_prev_time; 
  xvel_prev_time.set_zero(&mesh);
  yvel_prev_time.set_zero(&mesh);
  p_prev_time.set_zero(&mesh);
  Hermes::vector<Solution*> slns = Hermes::vector<Solution*>(&xvel_prev_time, &yvel_prev_time, 
                                                             &p_prev_time);

  // Initialize weak formulation.
  WeakForm* wf = new WeakFormDrivenCavity(Re, "Top", time_step, 
                                          &xvel_prev_time, &yvel_prev_time);

  // Initialize the FE problem.
  DiscreteProblem dp(wf, spaces);

  // 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 views.
  VectorView vview("velocity", new WinGeom(0, 0, 400, 400));
  ScalarView pview("pressure", new WinGeom(410, 0, 400, 400));
  //vview.set_min_max_range(0, 1.6);
  vview.fix_scale_width(80);
  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.
  scalar* coeff_vec = new scalar[Space::get_num_dofs(spaces)];
  info("Projecting initial condition to obtain initial vector for the Newton's method.");
  OGProjection::project_global(spaces, slns, coeff_vec, matrix_solver, 
                               Hermes::vector<ProjNormType>(vel_proj_norm, vel_proj_norm, p_proj_norm, t_proj_norm));

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

    // Perform Newton's iteration.
    info("Solving nonlinear problem:");
    bool verbose = true;
    if (!hermes_2D.solve_newton(coeff_vec, &dp, solver, matrix, rhs, 
        NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");

    // Update previous time level solutions.
    Solution::vector_to_solutions(coeff_vec, spaces, slns);
    // Show the solution at the end of time step.
    sprintf(title, "Velocity, time %g", current_time);
    vview.set_title(title);
    vview.show(&xvel_prev_time, &yvel_prev_time);
//.........这里部分代码省略.........

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

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

  // Perform initial mesh refinements.
  for (int i=0; i<INIT_REF; i++) mesh.refine_all_elements();
  
  // Create an L2 space with default shapeset.
  L2Space space(&mesh, bc_types, NULL, Ord2(P_H, P_V));
  int ndof = Space::get_num_dofs(&space);
  info("ndof = %d", ndof);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form));
  wf.add_vector_form(callback(linear_form));
  wf.add_matrix_form_surf(callback(bilinear_form_boundary), H2D_DG_BOUNDARY_EDGE);
  wf.add_vector_form_surf(callback(linear_form_boundary), H2D_DG_BOUNDARY_EDGE);
  wf.add_matrix_form_surf(callback(bilinear_form_interface), H2D_DG_INNER_EDGE);

  // 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).
  }
    
  // Initialize the solution.
  Solution sln;
  
  // Assemble the stiffness matrix and right-hand side vector.
  info("Assembling the stiffness matrix and right-hand side vector.");
  dp.assemble(matrix, rhs);
  
  // Solve the linear system and if successful, obtain the solution.
  info("Solving the matrix problem.");
  if(solver->solve())
    Solution::vector_to_solution(solver->get_solution(), &space, &sln);
  else
    error ("Matrix solver failed.\n");
  
  // Time measurement.
  cpu_time.tick();
  
  // Clean up.
  delete solver;
  delete matrix;
  delete rhs;
  
  info("ndof = %d", ndof);
  info("Coordinate ( 0.1, 0.1) value = %lf", sln.get_pt_value(0.1, 0.1));
  info("Coordinate ( 0.3, 0.3) value = %lf", sln.get_pt_value(0.3, 0.3));
  info("Coordinate ( 0.5, 0.5) value = %lf", sln.get_pt_value(0.5, 0.5));
  info("Coordinate ( 0.7, 0.7) value = %lf", sln.get_pt_value(0.7, 0.7));

  double coor_xy[4] = {0.1, 0.3, 0.5, 0.7};
  double value[4] = {0.999885, 0.844340, 0.000000, 0.000000};
  for (int i = 0; i < 4; i++)
  {
    if ((value[i] - sln.get_pt_value(coor_xy[i], coor_xy[i])) < 1E-6)
    {
      printf("Success!\n");
    }
    else
    {
      printf("Failure!\n");
      return ERR_FAILURE;
    }
  }
  return ERR_SUCCESS;
}

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

  // Perform initial mesh refinements.
  mesh.refine_towards_vertex(3, CORNER_REF_LEVEL);

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

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

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

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form));
  wf.add_vector_form(callback(linear_form));
  wf.add_vector_form_surf(callback(linear_form_surf_bottom), BDY_BOTTOM);
  wf.add_vector_form_surf(callback(linear_form_surf_outer), BDY_OUTER);
  wf.add_vector_form_surf(callback(linear_form_surf_left), BDY_LEFT);

  // 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 solution.
  Solution sln;

  // Assemble the stiffness matrix and right-hand side vector.
  info("Assembling the stiffness matrix and right-hand side vector.");
  dp.assemble(matrix, rhs);

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

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

  // Compute and show gradient magnitude.
  // (Note that the gradient at the re-entrant
  // corner needs to be truncated for visualization purposes.)
  ScalarView gradview("Gradient", new WinGeom(450, 0, 400, 350));
  MagFilter grad(Hermes::Tuple<MeshFunction *>(&sln, &sln), Hermes::Tuple<int>(H2D_FN_DX, H2D_FN_DY));
  gradview.show(&grad);

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

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

  return 0;
}

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("../square.mesh", &mesh);

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

  // Set exact solution.
  CustomExactSolution exact(&mesh);
  
  // Initialize the weak formulation.
  WeakFormPoisson wf1;

  // Initialize boundary conditions
  DefaultEssentialBCNonConst bc_essential(BDY_DIRICHLET, &exact);
  EssentialBCs bcs(&bc_essential);
 
  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bcs, P_INIT);
  int ndof = Space::get_num_dofs(&space);
  info("ndof: %d", ndof);

  info("---- Assembling by DiscreteProblem, solving by %s:", MatrixSolverNames[matrix_solver].c_str());

  // Initialize the solution.
  Solution sln1;
  
  // Initialize the linear discrete problem.
  bool is_linear = true;
  DiscreteProblem dp1(&wf1, &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);
  
  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver)->set_solver(iterative_method);
    ((AztecOOSolver*) solver)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }
  
  // Begin time measurement of assembly.
  cpu_time.tick(HERMES_SKIP);

  // Assemble the stiffness matrix and right-hand side vector.
  info("Assembling the stiffness matrix and right-hand side vector.");
  dp1.assemble(matrix, rhs);
  
  // Record assembly time.
  double time1 = cpu_time.tick().last();
  cpu_time.reset();

  // Solve the linear system and if successful, obtain the solution.
  info("Solving the matrix problem by %s.", MatrixSolverNames[matrix_solver].c_str());
  if(solver->solve())
    Solution::vector_to_solution(solver->get_solution(), &space, &sln1);
  else
    error ("Matrix solver failed.\n");

  // CPU time needed by UMFpack to solve the matrix problem.
  double time2 = cpu_time.tick().last();

  // Calculate errors.
  double rel_err_1 = hermes2d.calc_rel_error(&sln1, &exact, HERMES_H1_NORM) * 100;
  info("Assembly time: %g s, matrix solver time: %g s.", time1, time2);
  info("Xxact H1 error: %g%%.", rel_err_1);

  delete(matrix);
  delete(rhs);
  delete(solver);
    
  // View the solution and mesh.
  //ScalarView sview("Solution", new WinGeom(0, 0, 440, 350));
  //sview.show(&sln1);
  //OrderView  oview("Polynomial orders", new WinGeom(450, 0, 400, 350));
  //oview.show(&space);
  
  // TRILINOS PART:

  info("---- Assembling by DiscreteProblem, solving by NOX:");

  // Initialize the weak formulation for Trilinos.
  bool is_matrix_free = JFNK;
  WeakFormPoissonNox wf2(is_matrix_free);
  
  // Initialize DiscreteProblem.
  is_linear = false;
  DiscreteProblem dp2(&wf2, &space, is_linear);
//.........这里部分代码省略.........

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

    // Perform initial mesh refinements.
    for (int i=0; i < INIT_REF_NUM; i++) basemesh.refine_all_elements();
    basemesh.refine_towards_boundary(bdy_obstacle, INIT_REF_NUM_BDY, false); // 'true' stands for anisotropic refinements,
    basemesh.refine_towards_boundary(bdy_top, INIT_REF_NUM_BDY, true);       // 'false' for isotropic.
    basemesh.refine_towards_boundary(bdy_bottom, INIT_REF_NUM_BDY, true);
    mesh.copy(&basemesh);

    // Create spaces with default shapesets.
    H1Space* xvel_space = new H1Space(&mesh, xvel_bc_type, essential_bc_values_xvel, P_INIT_VEL);
    H1Space* yvel_space = new H1Space(&mesh, yvel_bc_type, essential_bc_values_yvel, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
    L2Space* p_space = new L2Space(&mesh, P_INIT_PRESSURE);
#else
    H1Space* p_space = new H1Space(&mesh, p_bc_type, NULL, P_INIT_PRESSURE);
#endif

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

    // Define projection norms.
    int vel_proj_norm = H2D_H1_NORM;
#ifdef PRESSURE_IN_L2
    int p_proj_norm = H2D_L2_NORM;
#else
    int p_proj_norm = H2D_H1_NORM;
#endif

    // Solutions for the Newton's iteration and time stepping.
    info("Setting initial conditions.");
//  Solution xvel_fine, yvel_fine, p_fine;
    Solution xvel_sln, yvel_sln, p_sln;
    Solution xvel_ref_sln, yvel_ref_sln, p_ref_sln;
    Solution xvel_prev_time, yvel_prev_time, p_prev_time;

    // Define initial conditions on the coarse mesh.
    xvel_prev_time.set_zero(&mesh);
    yvel_prev_time.set_zero(&mesh);
    p_prev_time.set_zero(&mesh);

    // Initialize the weak formulation.
    WeakForm wf(3);
    wf.add_matrix_form(0, 0, callback(bilinear_form_sym_0_0_1_1), H2D_SYM);
    wf.add_matrix_form(0, 0, callback(newton_bilinear_form_unsym_0_0), H2D_UNSYM, H2D_ANY);
    wf.add_matrix_form(0, 1, callback(newton_bilinear_form_unsym_0_1), H2D_UNSYM, H2D_ANY);
    wf.add_matrix_form(0, 2, callback(bilinear_form_unsym_0_2), H2D_ANTISYM);
    wf.add_matrix_form(1, 0, callback(newton_bilinear_form_unsym_1_0), H2D_UNSYM, H2D_ANY);
    wf.add_matrix_form(1, 1, callback(bilinear_form_sym_0_0_1_1), H2D_SYM);
    wf.add_matrix_form(1, 1, callback(newton_bilinear_form_unsym_1_1), H2D_UNSYM, H2D_ANY);
    wf.add_matrix_form(1, 2, callback(bilinear_form_unsym_1_2), H2D_ANTISYM);
    wf.add_vector_form(0, callback(newton_F_0), H2D_ANY, Tuple<MeshFunction*>(&xvel_prev_time, &yvel_prev_time));
    wf.add_vector_form(1, callback(newton_F_1), H2D_ANY, Tuple<MeshFunction*>(&xvel_prev_time, &yvel_prev_time));
    wf.add_vector_form(2, callback(newton_F_2), H2D_ANY);

    // Initialize adaptivity parameters.
    AdaptivityParamType apt(ERR_STOP, NDOF_STOP, THRESHOLD, STRATEGY, MESH_REGULARITY);
    // Create a selector which will select optimal candidate.
    H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
    // Assign initial condition to mesh.
    Vector *coeff_vec = new AVector(ndof);

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

        // Periodic global derefinements.
        if (ts > 1 && ts % UNREF_FREQ == 0) {
            info("Global mesh derefinement.");
            mesh.copy(&basemesh);
            xvel_space->set_uniform_order(P_INIT_VEL);
            yvel_space->set_uniform_order(P_INIT_VEL);
            p_space->set_uniform_order(P_INIT_PRESSURE);
        }

        // Update the coefficient vector and u_prev_time.
        info("Projecting to obtain coefficient vector on coarse mesh.");
        project_global(Tuple<Space *>(xvel_space, yvel_space, p_space),
                       Tuple<int>(vel_proj_norm, vel_proj_norm, p_proj_norm),
                       Tuple<MeshFunction*>(&xvel_prev_time, &yvel_prev_time, &p_prev_time),
                       Tuple<Solution*>(&xvel_prev_time, &yvel_prev_time, &p_prev_time),
                       coeff_vec);

        // Adaptivity loop (in space):
        bool verbose = true;     // Print info during adaptivity.
        info("Projecting coarse mesh solution to obtain initial vector on new fine mesh.");
        // The NULL pointers mean that we are not interested in visualization during the Newton's loop.
        solve_newton_adapt(Tuple<Space *>(xvel_space, yvel_space, p_space), &wf, coeff_vec, matrix_solver,
                           Tuple<int>(vel_proj_norm, vel_proj_norm, p_proj_norm),
                           Tuple<Solution *>(&xvel_sln, &yvel_sln, &p_sln),
                           Tuple<Solution *>(&xvel_ref_sln, &yvel_ref_sln, &p_ref_sln),
//.........这里部分代码省略.........