C++ Space::plot方法代码示例

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

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

int main() {

  // Create space.
  // Transform input data to the format used by the "Space" constructor.
  SpaceData *md = new SpaceData();
  Space* space = new Space(md->N_macroel, md->interfaces, md->poly_orders, md->material_markers, md->subdivisions, N_GRP, N_SLN);  
  delete md;
  
  // Enumerate basis functions, info for user.
  info("N_dof = %d", Space::get_num_dofs(space));
  // Plot the space.
  space->plot("space.gp");

  for (int g = 0; g < N_GRP; g++)  {
  	space->set_bc_right_dirichlet(g, flux_right_surf[g]);
	}
  
  // Initialize the weak formulation.
  WeakForm wf(2);
  wf.add_matrix_form(0, 0, jacobian_fuel_0_0, NULL, fuel);
  wf.add_matrix_form(0, 1, jacobian_fuel_0_1, NULL, fuel);
  wf.add_matrix_form(1, 0, jacobian_fuel_1_0, NULL, fuel);    
  wf.add_matrix_form(1, 1, jacobian_fuel_1_1, NULL, fuel);
    
  wf.add_vector_form(0, residual_fuel_0, NULL, fuel);  
  wf.add_vector_form(1, residual_fuel_1, NULL, fuel); 

  wf.add_vector_form_surf(0, residual_surf_left_0, BOUNDARY_LEFT);
  wf.add_vector_form_surf(1, residual_surf_left_1, BOUNDARY_LEFT);

  // 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)];
  solution_to_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(matrix, rhs);

    // Calculate the l2-norm of residual vector.
    double res_norm = 0;
    for(int i=0; i<ndof; i++) res_norm += rhs->get(i)*rhs->get(i);
    res_norm = sqrt(res_norm);

    // Info for user.
    info("---- Newton iter %d, residual norm: %.15f", it, res_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_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.
    vector_to_solution(coeff_vec, space);

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

	// Calculate flux integral for comparison with the reference value.
	double I = calc_integrated_flux(space, 1, 60., 80.);
	double Iref = 134.9238787715397;
	info("I = %.13f, err = %.13f%%", I, 100.*(I - Iref)/Iref );
	
  info("Done.");
  return 1;
}

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);
    wf.add_vector_form(residual);

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

    // Allocate coefficient vector.
    double *coeff_vec = new double[ndof];
    memset(coeff_vec, 0, ndof*sizeof(double));

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

    // Time stepping loop.
    double current_time = 0.0;
    while (current_time < T_FINAL)
    {
        // Newton's loop.
        // Fill vector coeff_vec using dof and coeffs arrays in elements.
        get_coeff_vector(space, coeff_vec);

        int it = 1;
        while (true)
        {
            // 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);
        char filename[100];
        sprintf(filename, "solution_%g.gp", current_time);
        l.plot_solution(filename);
        info("Solution %s written.", filename);

        current_time += TAU;
    }

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

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

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