当前位置: 首页>>代码示例>>C++>>正文


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

本文整理汇总了C++中LinearSolver::get_sln_vector方法的典型用法代码示例。如果您正苦于以下问题:C++ LinearSolver::get_sln_vector方法的具体用法?C++ LinearSolver::get_sln_vector怎么用?C++ LinearSolver::get_sln_vector使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在LinearSolver的用法示例。


在下文中一共展示了LinearSolver::get_sln_vector方法的10个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: main

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;
}
开发者ID:fauzisd,项目名称:hermes-examples,代码行数:54,代码来源:main.cpp

示例2: main


//.........这里部分代码省略.........
        if(Space<double>::get_num_dofs(*ref_spaces) == ndofs_prev)
          selector.set_error_weights(2.0 * selector.get_error_weight_h(), 1.0, 1.0);
        else
          selector.set_error_weights(1.0, 1.0, 1.0);

      ndofs_prev = Space<double>::get_num_dofs(*ref_spaces);

      // 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), 
        Hermes::vector<Solution<double>*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), matrix_solver_type, Hermes::vector<Hermes::Hermes2D::ProjNormType>(), iteration > 1);

      // 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<double>::get_num_dofs(*ref_spaces));

      // Assemble the reference problem.
      info("Solving on reference mesh.");
      DiscreteProblem<double> dp(&wf, *ref_spaces);

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

    wf.set_time_step(time_step);

    dp.assemble(matrix, rhs);
    
    // Solve the matrix problem.
    info("Solving the matrix problem.");
    if(solver->solve())
      if(!SHOCK_CAPTURING)
          Solution<double>::vector_to_solutions(solver->get_sln_vector(), *ref_spaces, 
          Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e));
      else
        {      
          FluxLimiter flux_limiter(FluxLimiter::Kuzmin, solver->get_sln_vector(), *ref_spaces, true);
          
          flux_limiter.limit_second_orders_according_to_detector(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
            &space_rho_v_y, &space_e));

          flux_limiter.limit_according_to_detector(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
            &space_rho_v_y, &space_e));

          flux_limiter.get_limited_solutions(Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e));
        }
    else
      error ("Matrix solver failed.\n");

      // Project the fine mesh solution onto the coarse mesh.
      info("Projecting reference solution on coarse mesh.");
      OGProjection<double>::project_global(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e), Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), 
        Hermes::vector<Solution<double>*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e), matrix_solver_type, 
        Hermes::vector<ProjNormType>(HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM)); 

      // Calculate element errors and total error estimate.
      info("Calculating error estimate.");
      Adapt<double>* adaptivity = new Adapt<double>(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e), Hermes::vector<ProjNormType>(HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM));
      double err_est_rel_total = adaptivity->calc_err_est(Hermes::vector<Solution<double>*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e),
        Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e)) * 100;

      CFL.calculate_semi_implicit(Hermes::vector<Solution<double> *>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), (*ref_spaces)[0]->get_mesh(), time_step);
开发者ID:JordanBlocher,项目名称:hermes-examples,代码行数:66,代码来源:main.cpp

示例3: main

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

//.........这里部分代码省略.........
开发者ID:fauzisd,项目名称:hermes-examples,代码行数:101,代码来源:main.cpp

示例4: main

int main(int argc, char* argv[])
{
  // Load mesh.
  load_mesh(mesh, "domain.xml", INIT_REF_NUM);

  // Create an H1 space with default shapeset.
  SpaceSharedPtr<double> space(new H1Space<double>(mesh, &essential_bcs, P_INIT));
  solver.set_weak_formulation(&weak_formulation);
  solver.set_space(space);

#pragma region Time stepping loop.
  /*
    solver.set_time_step(time_step_length);
    do
    {
    std::cout << "Time step: " << time_step_number << std::endl;

    #pragma region Spatial adaptivity loop.
    adaptivity.set_space(space);
    int adaptivity_step = 1;
    do
    {
    std::cout << "Adaptivity step: " << adaptivity_step << std::endl;
    #pragma region Construct globally refined reference mesh and setup reference space.
    MeshSharedPtr ref_mesh = ref_mesh_creator.create_ref_mesh();
    SpaceSharedPtr<double> ref_space = ref_space_creator.create_ref_space(space, ref_mesh);
    solver.set_space(ref_space);
    #pragma endregion
    try
    {
    // Solving.
    solver.solve(get_initial_Newton_guess(adaptivity_step, &weak_formulation, space, ref_space, sln_time_prev));
    Solution<double>::vector_to_solution(solver.get_sln_vector(), ref_space, sln_time_new);
    }
    catch(Exceptions::Exception& e)
    {
    std::cout << e.info();
    }
    catch(std::exception& e)
    {
    std::cout << e.what();
    }

    // Project the fine mesh solution onto the coarse mesh.
    OGProjection<double>::project_global(space, sln_time_new, sln_time_new_coarse);

    // Calculate element errors and error estimate.
    errorCalculator.calculate_errors(sln_time_new_coarse, sln_time_new);
    double error_estimate = errorCalculator.get_total_error_squared() * 100;
    std::cout << "Error estimate: " << error_estimate << "%" << std::endl;

    // Visualize the solution and mesh.
    display(sln_time_new, ref_space);

    // If err_est too large, adapt the mesh.
    if (error_estimate < ERR_STOP)
    break;
    else
    adaptivity.adapt(&refinement_selector);

    adaptivity_step++;
    }
    while(true);
    #pragma endregion

    #pragma region No adaptivity in space.
    try
    {
    // Solving.
    solver.solve(sln_time_prev);

    // Get the solution for visualization etc. from the coefficient vector.
    Solution<double>::vector_to_solution(solver.get_sln_vector(), space, sln_time_new);

    // Visualize the solution and mesh.
    display(sln_time_new, space);
    }
    catch(Exceptions::Exception& e)
    {
    std::cout << e.info();
    }
    catch(std::exception& e)
    {
    std::cout << e.what();
    }
    #pragma endregion

    sln_time_prev->copy(sln_time_new);

    // Increase current time and counter of time steps.
    current_time += time_step_length;
    time_step_number++;
    }
    while (current_time < T_FINAL);
    */
#pragma endregion

#pragma region No time stepping (= stationary problem).
  try
  {
//.........这里部分代码省略.........
开发者ID:HPeX,项目名称:hermes,代码行数:101,代码来源:main.cpp

示例5: main


//.........这里部分代码省略.........
  // 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.
      info("Solving nonlinear problem:");
      Hermes::Hermes2D::NewtonSolver<double> newton(&dp, matrix_solver_type);
      try
      {
        newton.solve(coeff_vec, NEWTON_TOL, NEWTON_MAX_ITER);
      }
      catch(Hermes::Exceptions::Exception e)
      {
        e.printMsg();
        error("Newton's iteration failed.");
      };

      // Update previous time level solutions.
      Solution<double>::vector_to_solutions(coeff_vec, Hermes::vector<Space<double>*>(&xvel_space, &yvel_space, &p_space),
                                    Hermes::vector<Solution<double>*>(&xvel_prev_time, &yvel_prev_time, &p_prev_time));
    }
    else 
    {
      // Linear solve.
      info("Assembling and solving linear problem.");
      dp.assemble(matrix, rhs, false);
      if(solver->solve())
        Solution<double>::vector_to_solutions(solver->get_sln_vector(),
                  Hermes::vector<Space<double>*>(&xvel_space, &yvel_space, &p_space),
                  Hermes::vector<Solution<double>*>(&xvel_prev_time, &yvel_prev_time, &p_prev_time));
      else
        error ("Matrix solver failed.\n");
    }

    // 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, HERMES_EPS_LOW);
    sprintf(title, "Pressure, time %g", current_time);
    pview.set_title(title);
    pview.show(&p_prev_time);
  }

  delete [] coeff_vec;
  delete matrix;
  delete rhs;
  delete solver;

  // Wait for all views to be closed.
  View::wait();
  return 0;
}
开发者ID:fauzisd,项目名称:hermes-examples,代码行数:101,代码来源:main.cpp

示例6: main

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

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

    mesh.refine_by_criterion(criterion, INIT_REF_CRITERION);

    MeshView m;
    m.show(&mesh);

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

    ScalarView view1("Solution - Discontinuous Galerkin FEM", new WinGeom(900, 0, 450, 350));
    ScalarView view2("Solution - Standard continuous FEM", new WinGeom(900, 400, 450, 350));

    if(WANT_DG)
    {
        // Create an L2 space.
        L2Space<double> space_l2(&mesh, P_INIT);

        // Initialize the solution.
        Solution<double> sln_l2;

        // Initialize the weak formulation.
        CustomWeakForm wf_l2(BDY_BOTTOM_LEFT);


        // Initialize the FE problem.
        DiscreteProblem<double> dp_l2(&wf_l2, &space_l2);

        info("Assembling Discontinuous Galerkin (nelem: %d, ndof: %d).", mesh.get_num_active_elements(), space_l2.get_num_dofs());
        dp_l2.assemble(matrix, rhs);

        // Solve the linear system. If successful, obtain the solution.
        info("Solving Discontinuous Galerkin.");
        if(solver->solve())
            if(DG_SHOCK_CAPTURING)
            {
                FluxLimiter flux_limiter(FluxLimiter::Kuzmin, solver->get_sln_vector(), &space_l2, true);

                flux_limiter.limit_second_orders_according_to_detector();

                flux_limiter.limit_according_to_detector();

                flux_limiter.get_limited_solution(&sln_l2);

                view1.set_title("Solution - limited Discontinuous Galerkin FEM");
            }
            else
                Solution<double>::vector_to_solution(solver->get_sln_vector(), &space_l2, &sln_l2);
        else
            error ("Matrix solver failed.\n");

        // View the solution.
        view1.show(&sln_l2);
    }
    if(WANT_FEM)
    {
        // Create an H1 space.
        H1Space<double> space_h1(&mesh, P_INIT);

        // Initialize the solution.
        Solution<double> sln_h1;

        // Initialize the weak formulation.
        CustomWeakForm wf_h1(BDY_BOTTOM_LEFT, false);


        // Initialize the FE problem.
        DiscreteProblem<double> dp_h1(&wf_h1, &space_h1);

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

        info("Assembling Continuous FEM (nelem: %d, ndof: %d).", mesh.get_num_active_elements(), space_h1.get_num_dofs());
        dp_h1.assemble(matrix, rhs);

        // Solve the linear system. If successful, obtain the solution.
        info("Solving Continuous FEM.");
        if(solver->solve())
            Solution<double>::vector_to_solution(solver->get_sln_vector(), &space_h1, &sln_h1);
        else
            error ("Matrix solver failed.\n");

        // View the solution.
        view2.show(&sln_h1);
    }

    // Clean up.
//.........这里部分代码省略.........
开发者ID:Zhonghua,项目名称:hermes-tutorial,代码行数:101,代码来源:main.cpp

示例7: main


//.........这里部分代码省略.........
        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);
      Vector<double>* rhs = create_vector<double>(matrix_solver_type);
      LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, matrix, rhs);

      // Initialize the FE problem.
      bool is_linear = true;
      DiscreteProblem<double> dp(wf, *ref_spaces);
      if(SEMI_IMPLICIT)
        static_cast<EulerEquationsWeakFormSemiImplicitCoupled*>(wf)->set_time_step(time_step);
      else
        static_cast<EulerEquationsWeakFormExplicitCoupled*>(wf)->set_time_step(time_step);

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

      // Solve the matrix problem.
      info("Solving the matrix problem.");
      if (solver->solve())
        Solution<double>::vector_to_solutions(solver->get_sln_vector(), *ref_spaces, 
        Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e, &rsln_c));
      else
        error ("Matrix solver failed.\n");

      Hermes::vector<Space<double>*> flow_spaces((*ref_spaces)[0], (*ref_spaces)[1], (*ref_spaces)[2], (*ref_spaces)[3]);

      double* flow_solution_vector = new double[Space<double>::get_num_dofs(flow_spaces)];

      OGProjection<double>::project_global(flow_spaces, Hermes::vector<MeshFunction<double> *>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), flow_solution_vector);

			FluxLimiter flux_limiter(FluxLimiter::Krivodonova, flow_solution_vector, flow_spaces);

			flux_limiter.limit_according_to_detector();

			flux_limiter.get_limited_solutions(Hermes::vector<Solution<double> *>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
			
      if(SHOCK_CAPTURING)
        flux_limiter.get_limited_solutions(Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e));
      
      // Project the fine mesh solution onto the coarse mesh.
      info("Projecting reference solution on coarse mesh.");
      OGProjection<double>::project_global(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x,
        &space_rho_v_y, &space_e, &space_c), Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e, &rsln_c),
        Hermes::vector<Solution<double>*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e, &sln_c), matrix_solver_type,
        Hermes::vector<ProjNormType>(HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM));

      util_time_step = time_step;
      if(SEMI_IMPLICIT)
        CFL.calculate_semi_implicit(Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), &mesh_flow, util_time_step);
      else
        CFL.calculate(Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), &mesh_flow, util_time_step);
开发者ID:davidquantum,项目名称:hermes-examples,代码行数:66,代码来源:main.cpp

示例8: main


//.........这里部分代码省略.........
          rsln_rho_v_x.own_mesh = false;
          delete rsln_rho_v_y.get_mesh();
          delete rsln_rho_v_y.get_space();
          rsln_rho_v_y.own_mesh = false;
          delete rsln_e.get_mesh();
          delete rsln_e.get_space();
          rsln_e.own_mesh = false;
        }
      }

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

      // Assemble the reference problem.
      info("Solving on reference mesh.");
      DiscreteProblem<double> dp(&wf, ref_spaces_const);

      SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver);
      Vector<double>* rhs = create_vector<double>(matrix_solver);
      LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver, matrix, rhs);

      wf.set_time_step(time_step);

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

      // Solve the matrix problem.
      info("Solving the matrix problem.");
      if(solver->solve())
        if(!SHOCK_CAPTURING)
          Solution<double>::vector_to_solutions(solver->get_sln_vector(), ref_spaces_const, 
          Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e));
        else
        {      
          FluxLimiter flux_limiter(FluxLimiter::Kuzmin, solver->get_sln_vector(), ref_spaces_const, true);
          
          flux_limiter.limit_second_orders_according_to_detector(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
            &space_rho_v_y, &space_e));
          
          flux_limiter.limit_according_to_detector(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
            &space_rho_v_y, &space_e));

          flux_limiter.get_limited_solutions(Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e));
        }
      else
        error ("Matrix solver failed.\n");
      
      // Project the fine mesh solution onto the coarse mesh.
      info("Projecting reference solution on coarse mesh.");
      OGProjection<double>::project_global(Hermes::vector<const Space<double> *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e), Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), 
        Hermes::vector<Solution<double>*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e), matrix_solver, 
        Hermes::vector<ProjNormType>(HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM)); 

      // Calculate element errors and total error estimate.
      info("Calculating error estimate.");
      Adapt<double>* adaptivity = new Adapt<double>(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e), Hermes::vector<ProjNormType>(HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM, HERMES_L2_NORM));
      double err_est_rel_total = adaptivity->calc_err_est(Hermes::vector<Solution<double>*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e),
        Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e)) * 100;

      CFL.calculate_semi_implicit(Hermes::vector<Solution<double> *>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), (*ref_spaces)[0]->get_mesh(), time_step);
开发者ID:certik,项目名称:hermes-examples,代码行数:66,代码来源:main.cpp

示例9: main

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;
}
开发者ID:JordanBlocher,项目名称:hermes-examples,代码行数:90,代码来源:main.cpp

示例10: main

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;
}
开发者ID:fauzisd,项目名称:hermes-examples,代码行数:77,代码来源:main.cpp


注:本文中的LinearSolver::get_sln_vector方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。