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

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

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

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

  // Previous time level solution (initialized by the initial condition).
  Solution u_prev_time(&mesh, init_cond);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(jac), HERMES_NONSYM, HERMES_ANY);
  wf.add_vector_form(callback(res), HERMES_ANY, &u_prev_time);

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

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

  // Initialize views.
  ScalarView sview("Solution", new WinGeom(0, 0, 500, 400));
  OrderView oview("Mesh", new WinGeom(510, 0, 460, 400));
  oview.show(&space);

  // Time stepping loop:
  double current_time = 0.0; int ts = 1;
  do 
  {
    info("---- Time step %d, t = %g s.", ts, current_time); ts++;

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

    // Update previous time level solution.
    Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);

    // Update time.
    current_time += TAU;

    // Show the new time level solution.
    char title[100];
    sprintf(title, "Solution, t = %g", current_time);
    sview.set_title(title);
    sview.show(&u_prev_time);
    oview.show(&space);
  } 
  while (current_time < T_FINAL);

  // Cleanup.
  delete [] coeff_vec;
  delete matrix;
  delete rhs;
  delete solver;

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

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

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

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

    // Enter boundary markers.
    BCTypes bc_types;
    bc_types.add_bc_dirichlet(Hermes::Tuple<int>(BDY_1, BDY_6));
    bc_types.add_bc_newton(Hermes::Tuple<int>(BDY_2, BDY_3, BDY_4, BDY_5));

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

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

    // Initialize the weak formulation.
    WeakForm wf;
    wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
    wf.add_matrix_form_surf(callback(bilinear_form_surf));
    wf.add_vector_form_surf(linear_form_surf, linear_form_surf_ord);

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

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

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

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

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

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

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

        // Time measurement.
        cpu_time.tick();

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

        // Time measurement.
        cpu_time.tick();

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

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

        // Calculate exact error,
        bool solutions_for_adapt = false;
        double err_exact_rel = adaptivity->calc_err_exact(&sln, &sln_exact, solutions_for_adapt) * 100;

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

        // Time measurement.
        cpu_time.tick();

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

int main(int argc, char* argv[]) 
{
  // This is a hack I used to run the code a dozen of times when plotting convergence graphs.
  if (argc > 1) {
    if (argv[1][0] == 'e') method = IE;
    else if (argv[1][0] == 's') method = SDIRK;
    else error("what are you doing?");
  }
  if (argc > 2) {
    TAU = std::atof(argv[2]);
  }

  // This is important to make sure we compare solution at exact same point in time when studying convergence.
  int N_STEP = std::ceil(T_FINAL / TAU);
  if (fabs(T_FINAL - N_STEP * TAU) > 1e-10) {
    error("bad choice of TAU");
  }

  info("t_final = %g, tau = %g, n = %i", T_FINAL, TAU, N_STEP);

  // 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(BDY_DIRICHLET, INIT_BDY_REF_NUM);

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_DIRICHLET);
  
  // Enter Dirichlet boudnary values.
  BCValues bc_values;
  bc_values.add_zero(BDY_DIRICHLET);

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

  // Previous time level solution (initialized by the initial condition).
  Solution u_prev_time(&mesh, exact_solution);

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

  // 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.
  ScalarView sview("Solution", new WinGeom(0, 0, 500, 400));
  OrderView oview("Mesh", new WinGeom(520, 0, 450, 400));
  oview.show(&space);

  if (method == IE) {
    info("IMPLICIT EULER METHOD");
    // Initialize the weak formulation.
    WeakForm wf;
    wf.add_matrix_form(callback(jac), HERMES_NONSYM, HERMES_ANY);
    wf.add_vector_form(callback(res), HERMES_ANY, &u_prev_time);

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

    // Time stepping loop:
    int ts = 0;
    do {
      info("---- Time step %d, t = %g s.", ++ts, TIME);

      info("We are computing solution at next time step TIME+TAU = %g s.", TIME+TAU);

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

      // Update previous time level solution.
      Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);

      // Update time.
      TIME = TIME + TAU;

      // Compute exact error.
      Solution exact_sln(&mesh, exact_solution);
      double exact_l2_error = calc_abs_error(&u_prev_time, &exact_sln, HERMES_L2_NORM);
      info("TIME: %g s.", TIME);
      info("Exact error in l2-norm: %g.", exact_l2_error);

      // Show the new time level solution.
      char title[100];
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable* bt = new ButcherTable(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());

  // Turn off adaptive time stepping if R-K method is not embedded.
  if (bt->is_embedded() == false && ADAPTIVE_TIME_STEP_ON == true) {
    warn("R-K method not embedded, turning off adaptive time stepping.");
    ADAPTIVE_TIME_STEP_ON = false;
  }

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

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

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

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

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

  // Convert initial condition into a Solution.
  Solution* sln_prev_time = new Solution(&mesh, init_cond);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(stac_jacobian), HERMES_NONSYM, HERMES_ANY, sln_prev_time);
  wf.add_vector_form(callback(stac_residual), HERMES_ANY, sln_prev_time);

  // Initialize the discrete problem.
  bool is_linear = false;
  DiscreteProblem dp_coarse(&wf, &space, is_linear);

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

  // Graph for time step history.
  SimpleGraph time_step_graph;
  if (ADAPTIVE_TIME_STEP_ON) info("Time step history will be saved to file time_step_history.dat.");
  
  // Time stepping loop.
  double current_time = time_step; int ts = 1;
  do 
  {
    info("Begin time step %d.", ts);
    // Periodic global derefinement.
    if (ts > 1 && ts % UNREF_FREQ == 0) 
    {
      info("Global mesh derefinement.");
      if (UNREF_LEVEL == 1) mesh.unrefine_all_elements();
      else mesh.copy(&basemesh);
      space.set_uniform_order(P_INIT);
      ndof = Space::get_num_dofs(&space);
    }

    // Spatial adaptivity loop. Note: sln_prev_time must not be 
    // changed during spatial adaptivity. 
    Solution ref_sln;
    Solution* time_error_fn;
    if (bt->is_embedded() == true) time_error_fn = new Solution(&mesh);
    else time_error_fn = NULL;
    bool done = false; int as = 1;
    double err_est;
    do {
      // Construct globally refined reference mesh and setup reference space.
      Space* ref_space = construct_refined_space(&space);

      // Initialize discrete problem on reference mesh.
      DiscreteProblem* ref_dp = new DiscreteProblem(&wf, ref_space);

      // Runge-Kutta step on the fine mesh.
      info("Runge-Kutta time step on fine mesh (t = %g s, tau = %g s, stages: %d).", 
         current_time, time_step, bt->get_size());
      bool verbose = true;
      bool is_linear = false;
      if (!rk_time_step(current_time, time_step, bt, sln_prev_time, &ref_sln, time_error_fn,
                        ref_dp, matrix_solver, verbose, is_linear, NEWTON_TOL_FINE, NEWTON_MAX_ITER)) {
        error("Runge-Kutta time step failed, try to decrease time step size.");
      }

      /* If ADAPTIVE_TIME_STEP_ON == true, estimate temporal error. 
         If too large or too small, then adjust it and restart the time step. */

      double rel_err_time;
      if (bt->is_embedded() == true) {
        info("Calculating temporal error estimate.");
//.........这里部分代码省略.........

int main(int argc, char **argv) {
	int res = ERR_SUCCESS;

#ifdef WITH_PETSC
	PetscInitialize(&argc, &argv, (char *) PETSC_NULL, PETSC_NULL);
#endif
	set_verbose(false);

	if (argc < 3) error("Not enough parameters");

	HcurlShapesetLobattoHex shapeset;

	printf("* Loading mesh '%s'\n", argv[1]);
	Mesh mesh;
	Mesh3DReader mesh_loader;
	if (!mesh_loader.load(argv[1], &mesh)) error("Loading mesh file '%s'\n", argv[1]);

	printf("* Setting the space up\n");
	HcurlSpace space(&mesh, &shapeset);
	space.set_bc_types(bc_types);

	int order;
	sscanf(argv[2], "%d", &order);
	int dir_x = order, dir_y = order, dir_z = order;
	order3_t o(dir_x, dir_y, dir_z);
	printf("  - Setting uniform order to (%d, %d, %d)\n", o.x, o.y ,o.z);
	space.set_uniform_order(o);

	int ndofs = space.assign_dofs();
	printf("  - Number of DOFs: %d\n", ndofs);

	printf("* Calculating a solution\n");

#if defined WITH_UMFPACK
	UMFPackMatrix mat;
	UMFPackVector rhs;
	UMFPackLinearSolver solver(&mat, &rhs);
#elif defined WITH_PARDISO
	PardisoMatrix mat;
	PardisoVector rhs;
	PardisoSolver solver(&mat, &rhs);
#elif defined WITH_PETSC
	PetscMatrix mat;
	PetscVector rhs;
	PetscLinearSolver solver(&mat, &rhs);
#elif defined WITH_MUMPS
	MumpsMatrix mat;
	MumpsVector rhs;
	MumpsSolver solver(&mat, &rhs);
#endif

	WeakForm wf;
	wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<ord_t, ord_t>, SYM);
	wf.add_matrix_form_surf(bilinear_form_surf<double, scalar>, bilinear_form_surf<ord_t, ord_t>);
	wf.add_vector_form(linear_form<double, scalar>, linear_form<ord_t, ord_t>);
	wf.add_vector_form_surf(linear_form_surf<double, scalar>, linear_form_surf<ord_t, ord_t>);

	LinearProblem lp(&wf, &space);

	// assemble stiffness matrix
	Timer assemble_timer("Assembling stiffness matrix");
	assemble_timer.start();
	lp.assemble(&mat, &rhs);
	assemble_timer.stop();

	// solve the stiffness matrix
	Timer solve_timer("Solving stiffness matrix");
	solve_timer.start();
	bool solved = solver.solve();
	solve_timer.stop();

//#ifdef OUTPUT_DIR
	mat.dump(stdout, "a");
	rhs.dump(stdout, "b");
//#endif

	if (solved) {
		scalar *s = solver.get_solution();

		Solution sln(&mesh);
		sln.set_coeff_vector(&space, s);

		printf("* Solution:\n");
		for (int i = 1; i <= ndofs; i++) {
			printf(" x[% 3d] = " SCALAR_FMT "\n", i, SCALAR(s[i]));
		}

		// output the measured values
		printf("%s: %s (%lf secs)\n", assemble_timer.get_name(), assemble_timer.get_human_time(), assemble_timer.get_seconds());
		printf("%s: %s (%lf secs)\n", solve_timer.get_name(), solve_timer.get_human_time(), solve_timer.get_seconds());

		// norm
		ExactSolution ex_sln(&mesh, exact_solution);
		double hcurl_sln_norm = hcurl_norm(&sln);
		double hcurl_err_norm = hcurl_error(&sln, &ex_sln);
		printf(" - Hcurl solution norm: % le\n", hcurl_sln_norm);
		printf(" - Hcurl error norm:    % le\n", hcurl_err_norm);

		double l2_sln_norm = l2_norm_hcurl(&sln);
		double l2_err_norm = l2_error_hcurl(&sln, &ex_sln);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  ExodusIIReader mloader;
  if (!mloader.load("iron-water.e", &mesh)) error("ExodusII mesh load failed.");

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(Hermes::vector<int>(WATER_2, IRON));
  bc_types.add_bc_neumann(WATER_1); 

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

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

  // Initialize the weak formulation
  WeakForm wf;
  wf.add_matrix_form(bilinear_form_water, bilinear_form_ord, HERMES_SYM, WATER_1);
  wf.add_matrix_form(bilinear_form_water, bilinear_form_ord, HERMES_SYM, WATER_2);
  wf.add_matrix_form(bilinear_form_iron, bilinear_form_ord, HERMES_SYM, IRON);
  wf.add_vector_form(linear_form_source, linear_form_ord, WATER_1);

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

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

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

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

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

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

    // Calculate element errors and total error estimate.
    info("Calculating error estimate."); 
    Adapt* adaptivity = new Adapt(&space);
    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++;
//.........这里部分代码省略.........

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

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

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

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

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

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

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

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

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

    // Exact error for testing purposes.
    double err_exact;

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

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

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

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

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


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

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

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

        // Cleanup.
        delete adaptivity;
    }

    if(err_exact > 3.00)
        success_test = 0;

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

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

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

  // Check the number of command-line parameters.
  if (argc < 2) {
    info("Use x, y, z, xy, xz, yz, or xyz as a command-line parameter.");
    error("Not enough command-line parameters.");
  }

  // Determine anisotropy type from the command-line parameter.
  ANISO_TYPE = parse_aniso_type(args[1]);

  // Load the mesh.
  Mesh mesh;
  H3DReader mloader;
  mloader.load("hex-0-1.mesh3d", &mesh);

  // Assign the lowest possible directional polynomial degrees so that the problem's NDOF >= 1.
  assign_poly_degrees();

  // Create an H1 space with default shapeset.
  info("Setting directional polynomial degrees %d, %d, %d.", P_INIT_X, P_INIT_Y, P_INIT_Z);
  H1Space space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));

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

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

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

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

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

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

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

    // Set up the solver, matrix, and rhs according to the solver selection.
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    
    // Initialize the preconditioner in the case of SOLVER_AZTECOO.
    if (matrix_solver == SOLVER_AZTECOO) 
    {
      ((AztecOOSolver*) solver)->set_solver(iterative_method);
      ((AztecOOSolver*) solver)->set_precond(preconditioner);
      // Using default iteration parameters (see solver/aztecoo.h).
    }
  
    // Assemble the reference problem.
    info("Assembling on reference mesh (ndof: %d).", Space::get_num_dofs(ref_space));
    dp.assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

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

    // Time measurement.
    cpu_time.tick();

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

    // Time measurement.
    cpu_time.tick();

    // Output solution and mesh with polynomial orders.
    if (solution_output) 
    {
      out_fn_vtk(&sln, "sln", as);
      out_orders_vtk(&space, "order", as);
//.........这里部分代码省略.........

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

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

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

  // Initialize the space according to the
  // command-line parameters passed.
	int o;
	sscanf(args[2], "%d", &o);
	Ord3 order(o, o, o);
	HcurlSpace space(&mesh, bc_types, NULL, order);

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

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

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

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

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

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

  ExactSolution ex_sln(&mesh, exact_solution);

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

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

  // Clean up.
  delete matrix;
  delete rhs;
  delete solver;
  delete adaptivity;
  
  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
	}
	else {
    info("Failure!");
    return ERR_FAILURE;
	}
}

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable* bt = new ButcherTable(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());

  // Turn off adaptive time stepping if R-K method is not embedded.
  if (bt->is_embedded() == false && ADAPTIVE_TIME_STEP_ON == true) {
    warn("R-K method not embedded, turning off adaptive time stepping.");
    ADAPTIVE_TIME_STEP_ON = false;
  }

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

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

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

  // Initialize an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, NULL, P_INIT);
  int ndof = Space::get_num_dofs(&space);
  info("ndof = %d.", ndof);
 
  // Convert initial condition into a Solution.
  Solution* sln_prev_time = new Solution(&mesh, TEMP_INIT);

  // Initialize weak formulation.
  WeakForm wf;
  wf.add_matrix_form(stac_jacobian_vol, stac_jacobian_vol_ord, HERMES_NONSYM, HERMES_ANY, sln_prev_time);
  wf.add_vector_form(stac_residual_vol, stac_residual_vol_ord, HERMES_ANY, sln_prev_time);
  wf.add_matrix_form_surf(stac_jacobian_bottom, stac_jacobian_bottom_ord, BDY_BOTTOM, sln_prev_time);
  wf.add_vector_form_surf(stac_residual_bottom, stac_residual_bottom_ord, BDY_BOTTOM, sln_prev_time);
  wf.add_matrix_form_surf(stac_jacobian_top, stac_jacobian_top_ord, BDY_TOP, sln_prev_time);
  wf.add_vector_form_surf(stac_residual_top, stac_residual_top_ord, BDY_TOP, sln_prev_time);

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

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

  // Visualize initial condition.
  char title[100];
  ScalarView sln_view("Initial condition", new WinGeom(0, 0, 1500, 360));
  OrderView ordview("Initial mesh", new WinGeom(0, 410, 1500, 360));
  ScalarView time_error_view("Temporal error", new WinGeom(0, 800, 1500, 360));
  time_error_view.fix_scale_width(40);
  ScalarView space_error_view("Spatial error", new WinGeom(0, 1220, 1500, 360));
  space_error_view.fix_scale_width(40);
  sln_view.show(sln_prev_time, HERMES_EPS_VERYHIGH);
  ordview.show(&space);

  // Graph for time step history.
  SimpleGraph time_step_graph;
  if (ADAPTIVE_TIME_STEP_ON) info("Time step history will be saved to file time_step_history.dat.");

  // Time stepping loop:
  double current_time = 0; int ts = 1;
  do 
  {
    info("Begin time step %d.", ts);
    // Periodic global derefinement.
    if (ts > 1 && ts % UNREF_FREQ == 0) 
    {
      info("Global mesh derefinement.");
      if (UNREF_LEVEL == 1) mesh.unrefine_all_elements();
      else mesh.copy(&basemesh);
      space.set_uniform_order(P_INIT);
      ndof = Space::get_num_dofs(&space);
    }

    // Spatial adaptivity loop. Note: sln_prev_time must not be 
    // changed during spatial adaptivity. 
    Solution ref_sln;
    Solution* time_error_fn;
    if (bt->is_embedded() == true) time_error_fn = new Solution(&mesh);
    else time_error_fn = NULL;
    bool done = false; int as = 1;
    double err_est;
    do {
      // Construct globally refined reference mesh and setup reference space.
      Space* ref_space = construct_refined_space(&space);

      // Initialize discrete problem on reference mesh.
      DiscreteProblem* ref_dp = new DiscreteProblem(&wf, ref_space);

      // Runge-Kutta step on the fine mesh.
      info("Runge-Kutta time step on fine mesh (t = %g s, tau = %g s, stages: %d).", 
//.........这里部分代码省略.........

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

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(Hermes::vector<std::string>(BDY_GROUND));
  bc_types.add_bc_newton(BDY_AIR);

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_const(BDY_GROUND, TEMP_INIT);

  // Initialize 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);
 
  // Previous time level solution (initialized by the external temperature).
  Solution u_prev_time(&mesh, TEMP_INIT);

  // Initialize weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(stac_jacobian_vol));
  wf.add_vector_form(callback(stac_residual_vol));
  wf.add_matrix_form_surf(callback(stac_jacobian_surf), BDY_AIR);
  wf.add_vector_form_surf(callback(stac_residual_surf), BDY_AIR);

  // Project the initial condition on the FE space to obtain initial solution coefficient vector.
  info("Projecting initial condition to translate initial condition into a vector.");
  scalar* coeff_vec = new scalar[ndof];
  OGProjection::project_global(&space, &u_prev_time, coeff_vec, matrix_solver);

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

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

  // Time stepping loop:
  double current_time = 0.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, tau = %g, stages: %d).", 
         current_time, time_step, bt.get_size());
    bool verbose = true;
    bool is_linear = true;
    if (!rk_time_step(current_time, time_step, &bt, coeff_vec, &dp, matrix_solver,
		      verbose, is_linear)) {
      error("Runge-Kutta time step failed, try to decrease time step size.");
    }

    // Convert coeff_vec into a new time level solution.
    Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);

    // Update time.
    current_time += time_step;

    // Show the new time level solution.
    char title[100];
    sprintf(title, "Time %3.2f, exterior temperature %3.5f", current_time, temp_ext(current_time));
    Tview.set_title(title);
    Tview.show(&u_prev_time);

    // Increase counter of time steps.
    ts++;
  } 
  while (current_time < T_FINAL);

  // Cleanup.
  delete [] coeff_vec;

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

// Main function.
int main(int argc, char* argv[])
{
  
  // Either use exact constitutive relations (slow) (method 0) or precalculate 
  // their linear approximations (faster) (method 1) or
  // precalculate their quintic polynomial approximations (method 2) -- managed by following loop "Initializing polynomial approximation".
  if (CONSTITUTIVE_TABLE_METHOD == 1)
    CONSTITUTIVE_TABLES_READY = get_constitutive_tables(ITERATIVE_METHOD);
  // Points to be used for polynomial approximation of K(h).
  double* points = new double[NUM_OF_INSIDE_PTS];

  // The van Genuchten + Mualem K(h) function is approximated by polynomials close to zero in case of CONSTITUTIVE_TABLE_METHOD=1.
  // In case of CONSTITUTIVE_TABLE_METHOD=2, all constitutive functions are approximated by polynomials.
  info("Initializing polynomial approximations.");
  for (int i=0; i < MATERIAL_COUNT; i++) {
    info("Processing layer %d", i);
    init_polynomials(6 + NUM_OF_INSIDE_PTS, LOW_LIMIT, points, NUM_OF_INSIDE_PTS, i);
  }
  POLYNOMIALS_READY = true;
  if (CONSTITUTIVE_TABLE_METHOD == 2) {
    CONSTITUTIVE_TABLES_READY = true ;
    //Assign table limit to global definition.
    TABLE_LIMIT = INTERVALS_4_APPROX[NUM_OF_INTERVALS-1];
  }
  

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

  // Load the mesh.
  Mesh mesh, basemesh;
  H2DReader mloader;
  mloader.load(mesh_file, &basemesh);
  
  // Initial refinements.
  //basemesh.refine_towards_boundary(1, 1);
  //basemesh.refine_towards_boundary(3, 1);
  //basemesh.refine_towards_boundary(4, 1);

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_1);
  bc_types.add_bc_neumann(Hermes::Tuple<int>(BDY_2, BDY_3, BDY_4));

  // Enter Dirichlet boundary values.
  BCValues bc_values(&TIME);
  bc_values.add_timedep_function(BDY_1, essential_bc_values);

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

  // Create a selector which will select optimal candidate.
  H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Solutions for the time stepping and the Newton's method.
  Solution sln, ref_sln, sln_prev_time, sln_prev_iter;
  
  // Assign the function f() to the fine mesh.
  sln_prev_time.set_exact(&mesh, init_cond);
  sln_prev_iter.set_exact(&mesh, init_cond);

  // Initialize the weak formulation.
  WeakForm wf;
  if (ITERATIVE_METHOD == 1) {
    if (TIME_INTEGRATION == 1) {
      info("Registering forms for the Newton's method (implicit Euler in time).");
      wf.add_matrix_form(jac_form_vol_euler, jac_form_vol_ord, HERMES_NONSYM, HERMES_ANY, 
	                 &sln_prev_time);
      wf.add_vector_form(res_form_vol_euler, res_form_vol_ord, HERMES_ANY, 
			 &sln_prev_time);
    }
    else {
      info("Registering forms for the Newton's method (Crank-Nicolson in time).");
      wf.add_matrix_form(jac_form_vol_cranic, jac_form_vol_ord, HERMES_NONSYM, HERMES_ANY, 
      		         &sln_prev_time);
      wf.add_vector_form(res_form_vol_cranic, res_form_vol_ord, HERMES_ANY, 
			 &sln_prev_time);
    }
  }
  else {
    if (TIME_INTEGRATION == 1) {
      info("Registering forms for the Picard's method (implicit Euler in time).");
      wf.add_matrix_form(bilinear_form_picard_euler, bilinear_form_picard_euler_ord, HERMES_NONSYM, HERMES_ANY, 
	                 &sln_prev_iter);
      wf.add_vector_form(linear_form_picard_euler, linear_form_picard_euler_ord, HERMES_ANY, 
			 Hermes::Tuple<MeshFunction*>(&sln_prev_iter, &sln_prev_time));
    }
    else {
      info("Registering forms for the Picard's method (Crank-Nicolson in time).");
      error("Not implemented yet.");
      wf.add_matrix_form(bilinear_form_picard_euler, bilinear_form_picard_euler_ord, HERMES_NONSYM, HERMES_ANY, 
	                 &sln_prev_iter);
//.........这里部分代码省略.........

int main()
{
    // Create space, set Dirichlet BC, enumerate basis functions.
    Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, P_INIT, NEQ);
    int ndof = Space::get_num_dofs(space);
    info("ndof: %d", ndof);

    // Initialize the weak formulation.
    WeakForm wf;
    wf.add_matrix_form(jacobian_vol);
    wf.add_vector_form(residual_vol);
    wf.add_vector_form_surf(0, residual_surf_right, BOUNDARY_RIGHT);

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

    // Newton's loop.
    // Fill vector coeff_vec using dof and coeffs arrays in elements.
    double *coeff_vec = new double[Space::get_num_dofs(space)];
    get_coeff_vector(space, coeff_vec);

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

    int it = 1;
    while (1) {
        // Obtain the number of degrees of freedom.
        int ndof = Space::get_num_dofs(space);

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

        // 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, then quit.
        // NOTE: at least one full iteration forced
        //       here because sometimes the initial
        //       residual on fine mesh is too small.
        if(res_l2_norm < NEWTON_TOL && it > 1) break;

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

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

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

        // If the maximum number of iteration has been reached, then quit.
        if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");

        // Copy coefficients from vector y to elements.
        set_coeff_vector(coeff_vec, space);

        it++;
    }

    // Plot the solution.
    Linearizer l(space);
    l.plot_solution("solution.gp");

    // Plot the resulting space.
    space->plot("space.gp");

    // Cleaning
    delete dp;
    delete rhs;
    delete solver;
    delete[] coeff_vec;
    delete space;
    delete bc;
    delete matrix;

    info("Done.");
    return 0;
}