C++ Linearizer::plot_solution方法代码示例

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);
  mesh->set_bc_right_dirichlet(0, Val_dir_right);
  mesh->assign_dofs();

  // Create discrete problem on coarse mesh
  DiscreteProblem *dp = new DiscreteProblem();
  dp->add_matrix_form(0, 0, jacobian);
  dp->add_vector_form(0, residual);

  // 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("max FTR error", "k", "--");

  // Main adaptivity loop
  int adapt_iterations = 1;
  double ftr_errors[MAX_ELEM_NUM];      // This array decides what 
                                         // elements will be refined.
  ElemPtr2 ref_ftr_pairs[MAX_ELEM_NUM]; // To store element pairs from the 
                                         // FTR solution. Decides how 
                                         // elements will be hp-refined. 
  for (int i=0; i < MAX_ELEM_NUM; i++) {
    ref_ftr_pairs[i][0] = new Element();
    ref_ftr_pairs[i][1] = new Element();
  }
  while(1) {
    printf("============ Adaptivity step %d ============\n", adapt_iterations); 

    printf("N_dof = %d\n", mesh->get_n_dof());
 
    // Newton's loop on coarse mesh
    newton(dp, mesh, NULL, NEWTON_TOL_COARSE, NEWTON_MAXITER);

    // For every element perform its fast trial refinement (FTR),
    // calculate the norm of the difference between the FTR
    // solution and the coarse mesh solution, and store the
    // error in the ftr_errors[] array.
    int n_elem = mesh->get_n_active_elem();
    for (int i=0; i < n_elem; i++) {

      printf("=== Starting FTR of Elem [%d]\n", i);

      // Replicate coarse mesh including solution.
      Mesh *mesh_ref_local = mesh->replicate();

      // Perform FTR of element 'i'
      mesh_ref_local->reference_refinement(i, 1);
      printf("Elem [%d]: fine mesh created (%d DOF).\n", 
             i, mesh_ref_local->assign_dofs());

      // Newton's loop on the FTR mesh
      newton(dp, mesh_ref_local, NULL, NEWTON_TOL_REF, NEWTON_MAXITER);

      // Print FTR solution (enumerated) 
      Linearizer *lxx = new Linearizer(mesh_ref_local);
      char out_filename[255];
      sprintf(out_filename, "solution_ref_%d.gp", i);
      lxx->plot_solution(out_filename);
      delete lxx;

      // Calculate norm of the difference between the coarse mesh 
      // and FTR solutions.
      // NOTE: later we want to look at the difference in some quantity 
      // of interest rather than error in global norm.
      double err_est_array[MAX_ELEM_NUM];
      ftr_errors[i] = calc_error_estimate(NORM, mesh, mesh_ref_local, 
                      err_est_array);
      //printf("Elem [%d]: absolute error (est) = %g\n", i, ftr_errors[i]);

      // Copy the reference element pair for element 'i'
      // into the ref_ftr_pairs[i][] array
      Iterator *I = new Iterator(mesh);
      Iterator *I_ref = new Iterator(mesh_ref_local);
      Element *e, *e_ref;
      while (1) {
        e = I->next_active_element();
        e_ref = I_ref->next_active_element();
        if (e->id == i) {
  	  e_ref->copy_into(ref_ftr_pairs[e->id][0]);
          // coarse element 'e' was split in space
          if (e->level != e_ref->level) {
            e_ref = I_ref->next_active_element();
            e_ref->copy_into(ref_ftr_pairs[e->id][1]);
          }
          break;
        }
      }

      delete I;
      delete I_ref;
      delete mesh_ref_local;
    }  

    // If exact solution available, also calculate exact error
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if(res_l2_norm < NEWTON_TOL_REF && 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_ref_local);

        it++;
      }
      
      // Cleanup.
      delete matrix;
      delete rhs;
      delete solver;
      delete dp;
      delete [] coeff_vec;

      // Print FTR solution (enumerated). 
      Linearizer *lxx = new Linearizer(space_ref_local);
      char out_filename[255];
      sprintf(out_filename, "solution_ref_%d.gp", i);
      lxx->plot_solution(out_filename);
      delete lxx;

      // Calculate norm of the difference between the coarse space 
      // and FTR solutions.
      // NOTE: later we want to look at the difference in some quantity 
      // of interest rather than error in global norm.
      double err_est_array[MAX_ELEM_NUM];
      elem_errors[i] = calc_err_est(NORM, space, space_ref_local, 
                       err_est_array) * 100;
      info("Elem [%d]: absolute error (est) = %g%%", i, elem_errors[i]);

      // Copy the reference element pair for element 'i'.
      // into the ref_elem_pairs[i][] array
      Iterator *I = new Iterator(space);
      Iterator *I_ref = new Iterator(space_ref_local);
      Element *e, *e_ref;
      while (1) {
        e = I->next_active_element();
        e_ref = I_ref->next_active_element();
        if (e->id == i) {
  	  e_ref->copy_into(ref_elem_pairs[e->id][0]);
          // coarse element 'e' was split in space.
          if (e->level != e_ref->level) {
            e_ref = I_ref->next_active_element();
            e_ref->copy_into(ref_elem_pairs[e->id][1]);
          }
          break;
        }
      }

      delete I;
      delete I_ref;