C++ LinearSolver::solve方法代码示例

void testDoubleQR()
{
    if (os_) *os_ << "testDoubleQR()\n";

    LinearSolver<LinearSolverType_QR> solver;

    ublas::matrix<double> A(2,2);
    A(0,0) = 1.;
    A(0,1) = 2.;
    A(1,0) = 3.;
    A(1,1) = 4.;

    ublas::vector<double> y(2);
    y(0) = 5.;
    y(1) = 11.;

    ublas::vector<double> x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    if (os_) *os_ << x(0) << " - 1. = " << x(0) - 1. << endl;

    unit_assert_equal(x(0), 1., 1e-14);
    unit_assert_equal(x(1), 2., 1e-14);
}

void testBanded()
{
    if (os_) *os_ << "testBanded()\n";

    LinearSolver<> solver;

    ublas::banded_matrix<double> A(2,2,1,1);
    A(0,0) = 1.;
    A(0,1) = 2.;
    A(1,0) = 3.;
    A(1,1) = 4.;

    ublas::vector<double> y(2);
    y(0) = 5.;
    y(1) = 11.;

    ublas::vector<double> x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    unit_assert_equal(x(0), 1., 1e-14);
    unit_assert_equal(x(1), 2., 1e-14);
}

void testBandedComplex()
{
    if (os_) *os_ << "testBandedComplex()\n";

    LinearSolver<> solver;

    ublas::banded_matrix< complex<double> > A(2,2,1,1);
    A(0,0) = 1.;
    A(0,1) = 2.;
    A(1,0) = 3.;
    A(1,1) = 4.;

    ublas::vector< complex<double> > y(2);
    y(0) = 5.;
    y(1) = 11.;

    ublas::vector< complex<double> > x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    unit_assert(norm(x(0)-1.) < 1e-14);
    unit_assert(norm(x(1)-2.) < 1e-14);
}

void testComplex()
{
    if (os_) *os_ << "testComplex()\n";

    LinearSolver<> solver;

    ublas::matrix< complex<double> > A(2,2);
    A(0,0) = 1;
    A(0,1) = 2;
    A(1,0) = 3;
    A(1,1) = 4;

    ublas::vector< complex<double> > y(2);
    y(0) = 5;
    y(1) = 11;

    ublas::vector< complex<double> > x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    unit_assert(x(0) == 1.);
    unit_assert(x(1) == 2.);
}

    boost::numeric::ublas::vector<T> solve(
        const boost::numeric::ublas::matrix<T>& A,
        const boost::numeric::ublas::vector<T>& x)
    {
        LinearSolver<LinearSolverType_QR> solver;

        boost::numeric::ublas::vector<T> y = solver.solve(A, x);

        return y;
    }

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

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

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

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

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

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

  Solution<double> sln; 

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

  // Select matrix solver.
  SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
  Vector<double>* rhs = create_vector<double>(matrix_solver_type);
  LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, matrix, rhs);

  // Assemble stiffness matrix and rhs.
  dp.assemble(matrix, rhs);
 
  // Solve the linear system of the reference problem. If successful, obtain the solutions.
  if(solver->solve()) Solution<double>::vector_to_solution(solver->get_sln_vector(), &space, &sln);
  else error ("Matrix solver failed.\n");

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

  // Calculate error wrt. exact solution.
  CustomExactSolution sln_exact(&mesh, K);
  double err = Global<double>::calc_abs_error(&sln, &sln_exact, HERMES_H1_NORM);
  printf("err = %g, err_squared = %g\n\n", err, err*err);
 
  // Wait for all views to be closed.
  View::wait();
  return 0;
}

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("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);
  ZeroSolutionVector F_sln(&mesh);
  Hermes::vector<Solution<double>*> slns(&E_sln, &F_sln);

  // Initialize the weak formulation.
  CustomWeakFormWave wf(time_step, C_SQUARED, &E_sln, &F_sln);
  
  // Initialize boundary conditions
  DefaultEssentialBCConst<double> bc_essential("Perfect conductor", 0.0);
  EssentialBCs<double> bcs(&bc_essential);

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

  info("ndof = %d.", Space<double>::get_num_dofs(spaces));

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

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

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

  // Time stepping loop.
  double current_time = 0; int ts = 1;
  do
  {
    // Perform one implicit Euler time step.
    info("Implicit Euler time step (t = %g s, time_step = %g s).", current_time, time_step);

    // First time assemble both the stiffness matrix and right-hand side vector,
    // then just the right-hand side vector.
    if (ts == 1) {
      info("Assembling the stiffness matrix and right-hand side vector.");
      dp.assemble(matrix, rhs);
      static char file_name[1024];
      sprintf(file_name, "matrix.m");
      FILE *f = fopen(file_name, "w");
      matrix->dump(f, "A");
      fclose(f);
    }
    else {
      info("Assembling the right-hand side vector (only).");
      dp.assemble(rhs);
    }

    // Solve the linear system and if successful, obtain the solution.
    info("Solving the matrix problem.");
    if(solver->solve()) Solution<double>::vector_to_solutions(solver->get_sln_vector(), spaces, slns);
    else error ("Matrix solver failed.\n");

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

//.........这里部分代码省略.........
        space_rho.unrefine_all_mesh_elements();
        if(CAND_LIST_FLOW == H2D_HP_ANISO)
          space_rho.adjust_element_order(-1, P_INIT_FLOW);
        space_rho_v_x.copy_orders(&space_rho);
        space_rho_v_y.copy_orders(&space_rho);
        space_e.copy_orders(&space_rho);
      }
      if(REFINEMENT_COUNT_CONCENTRATION > 0) {
        REFINEMENT_COUNT_CONCENTRATION = 0;
        space_c.unrefine_all_mesh_elements();
        space_c.adjust_element_order(-1, P_INIT_CONCENTRATION);
      }
    }

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

      // Construct globally refined reference mesh and setup reference space.
      int order_increase = 0;
      if(CAND_LIST_FLOW == H2D_HP_ANISO)
        order_increase = 1;
      Hermes::vector<Space<double> *>* ref_spaces = Space<double>::construct_refined_spaces(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e, &space_c), order_increase);
      if(CAND_LIST_FLOW != H2D_HP_ANISO)
        (*ref_spaces)[4]->adjust_element_order(+1, P_INIT_CONCENTRATION);

      // Project the previous time level solution onto the new fine mesh.
      info("Projecting the previous time level solution onto the new fine mesh.");
      OGProjection<double>::project_global(*ref_spaces, Hermes::vector<Solution<double>*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e, &prev_c), 
        Hermes::vector<Solution<double>*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e, &prev_c), matrix_solver_type);

      if(iteration == 1) {
        if(CAND_LIST_FLOW == H2D_HP_ANISO)
        {
          prev_rho.set_const((*ref_spaces)[4]->get_mesh(), RHO_EXT);
          prev_rho_v_x.set_const((*ref_spaces)[4]->get_mesh(), RHO_EXT * V1_EXT);
          prev_rho_v_y.set_const((*ref_spaces)[4]->get_mesh(), RHO_EXT * V2_EXT);
          prev_e.set_const((*ref_spaces)[4]->get_mesh(), QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));
        }
        prev_c.set_const((*ref_spaces)[4]->get_mesh(), 0.0);
      }

      if(as > 1) {
        delete rsln_rho.get_mesh();
        delete rsln_rho_v_x.get_mesh();
        delete rsln_rho_v_y.get_mesh();
        delete rsln_e.get_mesh();
        delete rsln_c.get_mesh();
      }

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

      // Very imporant, set the meshes for the flow as the same.
      (*ref_spaces)[1]->get_mesh()->set_seq((*ref_spaces)[0]->get_mesh()->get_seq());
      (*ref_spaces)[2]->get_mesh()->set_seq((*ref_spaces)[0]->get_mesh()->get_seq());
      (*ref_spaces)[3]->get_mesh()->set_seq((*ref_spaces)[0]->get_mesh()->get_seq());

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

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

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

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

  //1/2Exception test

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    // Construct globally refined reference mesh and setup reference space.
    Hermes::vector<Space<double> *>* ref_spaces = Space<double>::construct_refined_spaces(Hermes::vector<Space<double> *>(&xdisp, &ydisp, &temp));

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

    // Time measurement.
    cpu_time.tick();
    
    // Solve the linear system of the reference problem. If successful, obtain the solutions.
    if(solver->solve()) Solution::vector_to_solutions(solver->get_solution(), *ref_spaces, 
                                            Hermes::vector<Solution *>(&ref_xdisp_sln, &ref_ydisp_sln, &ref_temp_sln));
    else error ("Matrix solver failed.\n");
  
    // Time measurement.
    cpu_time.tick();

    // Project the fine mesh solution onto the coarse mesh.
    info("Projecting reference solution on coarse mesh.");
//.........这里部分代码省略.........

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  // Initialize boundary conditions.
  DefaultEssentialBCConst<std::complex<double> > bc(Hermes::vector<std::string>(BDY_BOTTOM, BDY_RIGHT, BDY_TOP, BDY_LEFT), std::complex<double>(0.0, 0.0));

  EssentialBCs<std::complex<double> > bcs(&bc);

  // Create an H1 space.
  H1Space<std::complex<double> >* phi_space = new H1Space<std::complex<double> >(&mesh, &bcs, P_INIT);
  H1Space<std::complex<double> >* psi_space = new H1Space<std::complex<double> >(&mesh, &bcs, P_INIT);

  int ndof = Space<std::complex<double> >::get_num_dofs(Hermes::vector<Space<std::complex<double> > *>(phi_space, psi_space));
  info("ndof = %d.", ndof);

  // Initialize previous time level solutions.
  ConstantSolution<std::complex<double> > phi_prev_time(&mesh, init_cond_phi);
  ConstantSolution<std::complex<double> > psi_prev_time(&mesh, init_cond_psi);

  // Initialize the weak formulation.
  WeakForm<std::complex<double> > wf(2);
  wf.add_matrix_form(0, 0, callback(biform_euler_0_0));
  wf.add_matrix_form(0, 1, callback(biform_euler_0_1));
  wf.add_matrix_form(1, 0, callback(biform_euler_1_0));
  wf.add_matrix_form(1, 1, callback(biform_euler_1_1));
  wf.add_vector_form(0, callback(liform_euler_0), HERMES_ANY, &phi_prev_time);
  wf.add_vector_form(1, callback(liform_euler_1), HERMES_ANY, &psi_prev_time);

  // Initialize views.
  ScalarView view("Psi", new WinGeom(0, 0, 600, 500));
  view.fix_scale_width(80);

  // Time stepping loop:
  int nstep = (int)(T_FINAL/TAU + 0.5);
  for(int ts = 1; ts <= nstep; ts++)
  {

    info("Time step %d:", ts);

    info("Solving linear system.");
    // Initialize the FE problem.
    bool is_linear = true;
    DiscreteProblem<std::complex<double> > dp(&wf, Hermes::vector<Space<double>* *>(phi_space, psi_space), is_linear);
   
    SparseMatrix<std::complex<double> >* matrix = create_matrix<std::complex<double> >(matrix_solver_type);
    Vector<std::complex<double> >* rhs = create_vector<std::complex<double> >(matrix_solver_type);
    LinearSolver<std::complex<double> >* solver = create_linear_solver<std::complex<double> >(matrix_solver_type, matrix, rhs);

    // 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<std::complex<double> >::vector_to_solutions(solver->get_sln_vector(), Hermes::vector<Space<std::complex<double> >*>(phi_space, psi_space), Hermes::vector<Solution<std::complex<double> > *>(&phi_prev_time, &psi_prev_time));
    else
      error ("Matrix solver failed.\n");

    // Show the new time level solution.
    char title[100];
    sprintf(title, "Time step %d", ts);
    view.set_title(title);
    view.show(&psi_prev_time);
  }

  // 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 u1_mesh, u2_mesh;
  MeshReaderH2D mloader;
  mloader.load("bracket.mesh", &u1_mesh);

  // Initial mesh refinements.
  u1_mesh.refine_element_id(1);
  u1_mesh.refine_element_id(4);

  // Create initial mesh for the vertical displacement component.
  // This also initializes the multimesh hp-FEM.
  u2_mesh.copy(&u1_mesh);

  // Initialize boundary conditions.
  DefaultEssentialBCConst<double> zero_disp(BDY_RIGHT, 0.0);
  EssentialBCs<double> bcs(&zero_disp);

  // Create x- and y- displacement space using the default H1 shapeset.
  H1Space<double> u1_space(&u1_mesh, &bcs, P_INIT);
  H1Space<double> u2_space(&u2_mesh, &bcs, P_INIT);
  info("ndof = %d.", Space<double>::get_num_dofs(Hermes::vector<Space<double> *>(&u1_space, &u2_space)));

  // Initialize the weak formulation.
  // NOTE; These weak forms are identical to those in example P01-linear/08-system.
  CustomWeakForm wf(E, nu, rho*g1, BDY_TOP, f0, f1);

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

  // Initialize coarse and reference mesh solutions.
  Solution<double> u1_sln, u2_sln, u1_ref_sln, u2_ref_sln;

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

  // Initialize views.
  ScalarView s_view_0("Solution (x-displacement)", new WinGeom(0, 0, 400, 350));
  s_view_0.show_mesh(false);
  ScalarView s_view_1("Solution (y-displacement)", new WinGeom(760, 0, 400, 350));
  s_view_1.show_mesh(false);
  OrderView  o_view_0("Mesh (x-displacement)", new WinGeom(410, 0, 340, 350));
  OrderView  o_view_1("Mesh (y-displacement)", new WinGeom(1170, 0, 340, 350));
  ScalarView mises_view("Von Mises stress [Pa]", new WinGeom(0, 405, 400, 350));

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

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

    // Construct globally refined reference mesh and setup reference space.
    Hermes::vector<Space<double> *>* ref_spaces = Space<double>::construct_refined_spaces(Hermes::vector<Space<double> *>(&u1_space, &u2_space));

    // Initialize matrix solver.
    SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
    Vector<double>* rhs = create_vector<double>(matrix_solver_type);
    LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, matrix, rhs);

    // Assemble the reference problem.
    info("Solving on reference mesh.");
    DiscreteProblem<double> dp(&wf, *ref_spaces);
    dp.assemble(matrix, rhs);
    // Time measurement.
    cpu_time.tick();
    
    // Solve the linear system of the reference problem. If successful, obtain the solutions.
    if(solver->solve()) Solution<double>::vector_to_solutions(solver->get_sln_vector(), *ref_spaces, 
                                            Hermes::vector<Solution *>(&u1_ref_sln, &u2_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<double>::project_global(Hermes::vector<Space<double> *>(&u1_space, &u2_space), 
                                 Hermes::vector<Solution<double> *>(&u1_ref_sln, &u2_ref_sln), 
                                 Hermes::vector<Solution<double> *>(&u1_sln, &u2_sln), matrix_solver_type); 
   
    // View the coarse mesh solution and polynomial orders.
    s_view_0.show(&u1_sln); 
    o_view_0.show(&u1_space);
    s_view_1.show(&u2_sln); 
    o_view_1.show(&u2_space);
    // For von Mises stress Filter.
    double lambda = (E * nu) / ((1 + nu) * (1 - 2*nu));
    double mu = E / (2*(1 + nu));
    VonMisesFilter stress(Hermes::vector<MeshFunction<double> *>(&u1_sln, &u2_sln), lambda, mu);
    mises_view.show(&stress, HERMES_EPS_HIGH, H2D_FN_VAL_0, &u1_sln, &u2_sln, 1e4);

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

bool DuffingFloquet()
{
  int np = MPIComm::world().getNProc();
  TEUCHOS_TEST_FOR_EXCEPT(np != 1);

  const double pi = 4.0*atan(1.0);

  /* We will do our linear algebra using Epetra */
  VectorType<double> vecType = new EpetraVectorType();

  /* Create a periodic mesh */
  int nx = 128;

  MeshType meshType = new PeriodicMeshType1D();
  MeshSource mesher = new PeriodicLineMesher(0.0, 2.0*pi, nx, meshType);
  Mesh mesh = mesher.getMesh();

  /* Create a cell filter that will identify the maximal cells
   * in the interior of the domain */
  CellFilter interior = new MaximalCellFilter();
  CellFilter pts = new DimensionalCellFilter(0);
      
  CellFilter left = pts.subset(new CoordinateValueCellPredicate(0,0.0));
  CellFilter right = pts.subset(new CoordinateValueCellPredicate(0,2.0*pi));
      
  /* Create unknown and test functions, discretized using first-order
   * Lagrange interpolants */
  Expr u1 = new UnknownFunction(new Lagrange(1), "u1");
  Expr u2 = new UnknownFunction(new Lagrange(1), "u2");
  Expr v1 = new TestFunction(new Lagrange(1), "v1");
  Expr v2 = new TestFunction(new Lagrange(1), "v2");

  /* Create differential operator and coordinate function */
  Expr dx = new Derivative(0);
  Expr x = new CoordExpr(0);

  /* We need a quadrature rule for doing the integrations */
  QuadratureFamily quad = new GaussianQuadrature(4);

  double F0 = 0.5;
  double gamma = 2.0/3.0;
  double a0 = 1.0;
  double w0 = 1.0;
  double eps = 0.5;

  Expr u1Guess = -0.75*cos(x) + 0.237*sin(x);
  Expr u2Guess = 0.237*cos(x) + 0.75*sin(x);

  DiscreteSpace discSpace(mesh, 
    List(new Lagrange(1), new Lagrange(1)),
    vecType);
  L2Projector proj(discSpace, List(u1Guess, u2Guess));
  Expr u0 = proj.project();


  Expr rhs1 = u2;
  Expr rhs2 = -w0*w0*u1 - gamma*u2 - eps*w0*w0*pow(u1,3.0)/a0/a0 
    + F0*w0*w0*sin(x);

  /* Define the weak form */
  Expr eqn = Integral(interior, 
    v1*(dx*u1 - rhs1) + v2*(dx*u2 - rhs2),
    quad);
  Expr dummyBC ; 

  NonlinearProblem prob(mesh, eqn, dummyBC, List(v1,v2), List(u1,u2), 
    u0, vecType);


  ParameterXMLFileReader reader("nox.xml");
  ParameterList solverParams = reader.getParameters();

  Out::root() << "finding periodic solution" << endl;
  NOXSolver solver(solverParams);
  prob.solve(solver);

  /* unfold the solution onto a non-periodic mesh */
      
  Expr uP = unfoldPeriodicDiscreteFunction(u0, "u_p");
  Out::root() << "uP=" << uP << endl;
      
  Mesh unfoldedMesh = DiscreteFunction::discFunc(uP)->mesh();
  DiscreteSpace unfDiscSpace = DiscreteFunction::discFunc(uP)->discreteSpace();

  FieldWriter writer = new MatlabWriter("Floquet.dat");
  writer.addMesh(unfoldedMesh);
  writer.addField("u_p[0]", new ExprFieldWrapper(uP[0]));
  writer.addField("u_p[1]", new ExprFieldWrapper(uP[1]));

  Array<Expr> a(2);
  a[0] = new Sundance::Parameter(0.0, "a1");
  a[1] = new Sundance::Parameter(0.0, "a2");


  Expr bc = EssentialBC(left, v1*(u1-uP[0]-a[0]) + v2*(u2-uP[1]-a[1]), quad);

  NonlinearProblem unfProb(unfoldedMesh, eqn, bc, 
    List(v1,v2), List(u1,u2), uP, vecType);

  unfProb.setEvalPoint(uP);
//.........这里部分代码省略.........

// Test code.
void solve(LinearSolver<double> &solver, int n) {
  if (!solver.solve())
    printf("Unable to solve.\n");
}