C++ DiscreteProblem::add_vector_form方法代码示例

int main() {
 // Create coarse mesh
  MeshData *md = new MeshData();		// transform input data to the format used by the "Mesh" constructor
  Mesh *mesh = new Mesh(md->N_macroel, md->interfaces, md->poly_orders, md->material_markers, md->subdivisions, N_GRP, N_SLN);  
  delete md;
  
  printf("N_dof = %d\n", mesh->assign_dofs());
  mesh->plot("mesh.gp");

  for (int g = 0; g < N_GRP; g++)  {
  	mesh->set_bc_left_dirichlet(g, flux_left_surf[g]);
  	mesh->set_bc_right_dirichlet(g, flux_right_surf[g]);
	}
  
  // Register weak forms
  DiscreteProblem *dp = new DiscreteProblem();
  
  dp->add_matrix_form(0, 0, jacobian_mat1_0_0, mat1);
  dp->add_matrix_form(0, 0, jacobian_mat2_0_0, mat2);
  dp->add_matrix_form(0, 0, jacobian_mat3_0_0, mat3);
  
  dp->add_matrix_form(0, 1, jacobian_mat1_0_1, mat1);
  dp->add_matrix_form(0, 1, jacobian_mat2_0_1, mat2);
  dp->add_matrix_form(0, 1, jacobian_mat3_0_1, mat3);
  
  dp->add_matrix_form(1, 0, jacobian_mat1_1_0, mat1);    
  dp->add_matrix_form(1, 0, jacobian_mat2_1_0, mat2);
  dp->add_matrix_form(1, 0, jacobian_mat3_1_0, mat3);
    
  dp->add_matrix_form(1, 1, jacobian_mat1_1_1, mat1);
	dp->add_matrix_form(1, 1, jacobian_mat2_1_1, mat2);
	dp->add_matrix_form(1, 1, jacobian_mat3_1_1, mat3);
  
  dp->add_vector_form(0, residual_mat1_0, mat1);  
	dp->add_vector_form(0, residual_mat2_0, mat2);  
	dp->add_vector_form(0, residual_mat3_0, mat3);
	    
  dp->add_vector_form(1, residual_mat1_1, mat1);
  dp->add_vector_form(1, residual_mat2_1, mat2); 
  dp->add_vector_form(1, residual_mat3_1, mat3);  
	  	
  // Newton's loop		
  newton(dp, mesh, NULL, NEWTON_TOL, NEWTON_MAXITER, verbose);
	 
  // Plot the resulting neutron flux
  Linearizer l(mesh);
  l.plot_solution("solution.gp");
	
  printf("Done.\n");
  return 1;
}

int main() {
    // create space
    Space Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, P_INIT, NEQ);
    info("N_dof = %d", Space::get_num_dofs(&space));

    // Initialize the weak formulation.
    WeakForm wf;

    // Initialize the FE problem.
    DiscreteProblem *dp = new DiscreteProblem();
    dp->add_matrix_form(0, 0, jacobian_1_1);
    dp->add_matrix_form(0, 2, jacobian_1_3);
    dp->add_matrix_form(0, 3, jacobian_1_4);
    dp->add_matrix_form(1, 1, jacobian_2_2);
    dp->add_matrix_form(1, 2, jacobian_2_3);
    dp->add_matrix_form(1, 3, jacobian_2_4);
    dp->add_matrix_form(2, 0, jacobian_3_1);
    dp->add_matrix_form(2, 1, jacobian_3_2);
    dp->add_matrix_form(2, 2, jacobian_3_3);
    dp->add_matrix_form(3, 0, jacobian_4_1);
    dp->add_matrix_form(3, 1, jacobian_4_2);
    dp->add_matrix_form(3, 3, jacobian_4_4);
    dp->add_vector_form(0, residual_1);
    dp->add_vector_form(1, residual_2);
    dp->add_vector_form(2, residual_3);
    dp->add_vector_form(3, residual_4);
    dp->add_matrix_form_surf(0, 0, jacobian_surf_right_U_Re, BOUNDARY_RIGHT);
    dp->add_matrix_form_surf(0, 2, jacobian_surf_right_U_Im, BOUNDARY_RIGHT);
    dp->add_matrix_form_surf(1, 1, jacobian_surf_right_I_Re, BOUNDARY_RIGHT);
    dp->add_matrix_form_surf(1, 3, jacobian_surf_right_I_Im, BOUNDARY_RIGHT);

    // Newton's loop
    newton(dp, space, MATRIX_SOLVER, MATRIX_SOLVER_TOL, MATRIX_SOLVER_MAXITER,
           NEWTON_TOL, NEWTON_MAXITER);

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

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

int main() {
// Create coarse mesh
    MeshData *md = new MeshData();		// transform input data to the format used by the "Mesh" constructor
    Mesh *mesh = new Mesh(md->N_macroel, md->interfaces, md->poly_orders, md->material_markers, md->subdivisions, N_GRP, N_SLN);
    delete md;

    printf("N_dof = %d\n", mesh->assign_dofs());
    mesh->plot("mesh.gp");

    for (int g = 0; g < N_GRP; g++)  {
        mesh->set_bc_right_dirichlet(g, flux_right_surf[g]);
    }

    // Register weak forms
    DiscreteProblem *dp = new DiscreteProblem();

    dp->add_matrix_form(0, 0, jacobian_fuel_0_0, fuel);
    dp->add_matrix_form(0, 1, jacobian_fuel_0_1, fuel);
    dp->add_matrix_form(1, 0, jacobian_fuel_1_0, fuel);
    dp->add_matrix_form(1, 1, jacobian_fuel_1_1, fuel);

    dp->add_vector_form(0, residual_fuel_0, fuel);
    dp->add_vector_form(1, residual_fuel_1, fuel);

    dp->add_vector_form_surf(0, residual_surf_left_0, BOUNDARY_LEFT);
    dp->add_vector_form_surf(1, residual_surf_left_1, BOUNDARY_LEFT);

    // Newton's loop
    newton(dp, mesh, NULL, NEWTON_TOL, NEWTON_MAXITER, verbose);

    // Plot the resulting neutron flux
    Linearizer l(mesh);
    l.plot_solution("solution.gp");

    // Calculate flux integral for comparison with the reference value
    double I = calc_integrated_flux(mesh, 1, 60., 80.);
    double Iref = 134.9238787715397;
    printf("I = %.13f, err = %.13f%%\n", I, 100.*(I - Iref)/Iref );

    printf("Done.\n");
    return 1;
}

// ********************************************************************
int main() {
  // Create coarse mesh, set Dirichlet BC, enumerate basis functions
  Mesh *mesh = new Mesh(A, B, N_elem, P_init, N_eq);
  mesh->set_bc_left_dirichlet(0, Val_dir_left_0);
  mesh->set_bc_left_dirichlet(1, Val_dir_left_1);
  printf("N_dof = %d\n", mesh->assign_dofs());

  // Create discrete problem on coarse mesh
  DiscreteProblem *dp = new DiscreteProblem();
  dp->add_matrix_form(0, 0, jacobian_0_0);
  dp->add_matrix_form(0, 1, jacobian_0_1);
  dp->add_matrix_form(1, 0, jacobian_1_0);
  dp->add_matrix_form(1, 1, jacobian_1_1);
  dp->add_vector_form(0, residual_0);
  dp->add_vector_form(1, residual_1);

  // scipy umfpack solver
  CommonSolverSciPyUmfpack solver;

  // Initial Newton's loop on coarse mesh
  newton(dp, mesh, &solver, NEWTON_TOL_COARSE, NEWTON_MAXITER);

  // Replicate coarse mesh including dof arrays
  Mesh *mesh_ref = mesh->replicate();

  // Refine entire mesh_ref uniformly in 'h' and 'p'
  int start_elem_id = 0; 
  int num_to_ref = mesh_ref->get_n_active_elem();
  mesh_ref->reference_refinement(start_elem_id, num_to_ref);
  printf("Fine mesh created (%d DOF).\n", mesh_ref->get_n_dof());

  // Convergence graph wrt. the number of degrees of freedom
  GnuplotGraph graph;
  graph.set_log_y();
  graph.set_captions("Convergence History", "Degrees of Freedom", "Error [%]");
  graph.add_row("exact error", "k", "-", "o");
  graph.add_row("error estimate", "k", "--");

  // Main adaptivity loop
  int adapt_iterations = 1;
  while(1) {
    printf("============ Adaptivity step %d ============\n", adapt_iterations); 

    // Newton's loop on fine mesh
    newton(dp, mesh_ref, &solver, NEWTON_TOL_REF, NEWTON_MAXITER);

    // Starting with second adaptivity step, obtain new coarse 
    // mesh solution via Newton's method. Initial condition is 
    // the last coarse mesh solution.
    if (adapt_iterations > 1) {

      // Newton's loop on coarse mesh
      newton(dp, mesh, &solver, NEWTON_TOL_COARSE, NEWTON_MAXITER);
    }

    // In the next step, estimate element errors based on 
    // the difference between the fine mesh and coarse mesh solutions. 
    double err_est_array[MAX_ELEM_NUM]; 
    double err_est_total = calc_error_estimate(NORM, mesh, mesh_ref, 
                           err_est_array);

    // Calculate the norm of the fine mesh solution
    double ref_sol_norm = calc_solution_norm(NORM, mesh_ref);

    // Calculate an estimate of the global relative error
    double err_est_rel = err_est_total/ref_sol_norm;
    printf("Relative error (est) = %g %%\n", 100.*err_est_rel);

    // If exact solution available, also calculate exact error
    double err_exact_rel;  
    if (EXACT_SOL_PROVIDED) {
      // Calculate element errors wrt. exact solution
      double err_exact_total = calc_error_exact(NORM, mesh, exact_sol);
     
      // Calculate the norm of the exact solution
      // (using a fine subdivision and high-order quadrature)
      int subdivision = 500; // heuristic parameter
      int order = 20;        // heuristic parameter
      double exact_sol_norm = calc_solution_norm(NORM, exact_sol, N_eq, A, B,
                                                  subdivision, order);
      // Calculate an estimate of the global relative error
      err_exact_rel = err_exact_total/exact_sol_norm;
      printf("Relative error (exact) = %g %%\n", 100.*err_exact_rel);
      graph.add_values(0, mesh->get_n_dof(), 100 * err_exact_rel);
    }

    // add entry to DOF convergence graph
    graph.add_values(1, mesh->get_n_dof(), 100 * err_est_rel);

    // Decide whether the relative error is sufficiently small
    if(err_est_rel*100 < TOL_ERR_REL) break;

    // debug
    // (adapt_iterations == 2) break;
 
    // Returns updated coarse and fine meshes, with the last 
    // coarse and fine mesh solutions on them, respectively. 
    // The coefficient vectors and numbers of degrees of freedom 
    // on both meshes are also updated. 
//.........这里部分代码省略.........