C++ ErrorVector::l2_norm方法代码示例

int main(int argc, char** argv)
{
  // Initialize the library.  This is necessary because the library
  // may depend on a number of other libraries (i.e. MPI and PETSc)
  // that require initialization before use.  When the LibMeshInit
  // object goes out of scope, other libraries and resources are
  // finalized.
  LibMeshInit init (argc, argv);

  // Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
  libmesh_example_requires(false, "--enable-amr");
#else

  // Create a mesh, with dimension to be overridden later, on the
  // default MPI communicator.
  Mesh mesh(init.comm());

  GetPot command_line (argc, argv);

  int n = 4;
  if ( command_line.search(1, "-n") )
    n = command_line.next(n);

  // Build a 1D mesh with 4 elements from x=0 to x=1, using
  // EDGE3 (i.e. quadratic) 1D elements. They are called EDGE3 elements
  // because a quadratic element contains 3 nodes.
  MeshTools::Generation::build_line(mesh,n,0.,1.,EDGE3);

  // Define the equation systems object and the system we are going
  // to solve. See Introduction Example 2 for more details.
  EquationSystems equation_systems(mesh);
  LinearImplicitSystem& system = equation_systems.add_system
    <LinearImplicitSystem>("1D");

  // Add a variable "u" to the system, using second-order approximation
  system.add_variable("u",SECOND);

  // Give the system a pointer to the matrix assembly function. This
  // will be called when needed by the library.
  system.attach_assemble_function(assemble_1D);

  // Define the mesh refinement object that takes care of adaptively
  // refining the mesh.
  MeshRefinement mesh_refinement(mesh);

  // These parameters determine the proportion of elements that will
  // be refined and coarsened. Any element within 30% of the maximum
  // error on any element will be refined, and any element within 30%
  // of the minimum error on any element might be coarsened
  mesh_refinement.refine_fraction()  = 0.7;
  mesh_refinement.coarsen_fraction() = 0.3;
  // We won't refine any element more than 5 times in total
  mesh_refinement.max_h_level()      = 5;

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Refinement parameters
  const unsigned int max_r_steps = 5; // Refine the mesh 5 times

  // Define the refinement loop
  for(unsigned int r_step=0; r_step<=max_r_steps; r_step++)
    {
      // Solve the equation system
      equation_systems.get_system("1D").solve();

      // We need to ensure that the mesh is not refined on the last iteration
      // of this loop, since we do not want to refine the mesh unless we are
      // going to solve the equation system for that refined mesh.
      if(r_step != max_r_steps)
        {
          // Error estimation objects, see Adaptivity Example 2 for details
          ErrorVector error;
          KellyErrorEstimator error_estimator;

          // Compute the error for each active element
          error_estimator.estimate_error(system, error);

          // Output error estimate magnitude
          libMesh::out << "Error estimate\nl2 norm = " << error.l2_norm() <<
            "\nmaximum = " << error.maximum() << std::endl;

          // Flag elements to be refined and coarsened
          mesh_refinement.flag_elements_by_error_fraction (error);

          // Perform refinement and coarsening
          mesh_refinement.refine_and_coarsen_elements();

          // Reinitialize the equation_systems object for the newly refined
          // mesh. One of the steps in this is project the solution onto the
          // new mesh
          equation_systems.reinit();
        }
    }

  // Construct gnuplot plotting object, pass in mesh, title of plot
  // and boolean to indicate use of grid in plot. The grid is used to
  // show the edges of each element in the mesh.
  GnuPlotIO plot(mesh,"Adaptivity Example 1", GnuPlotIO::GRID_ON);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                UniformRefinementEstimator *u =
                    new UniformRefinementEstimator;

                // The lid-driven cavity problem isn't in H1, so
                // lets estimate L2 error
                u->error_norm = L2;

                error_estimator.reset(u);
            }
            else
            {
                // If we aren't adapting to a tolerance we need a
                // target mesh size
                libmesh_assert_greater (nelem_target, 0);

                // Kelly is a lousy estimator to use for a problem
                // not in H1 - if we were doing more than a few
                // timesteps we'd need to turn off or limit the
                // maximum level of our adaptivity eventually
                error_estimator.reset(new KellyErrorEstimator);
            }

            // Calculate error
            std::vector<Real> weights(9,1.0);  // based on u, v, p, c, their adjoints, and source parameter

            // Keep the same default norm type.
            std::vector<FEMNormType>
            norms(1, error_estimator->error_norm.type(0));
            error_estimator->error_norm = SystemNorm(norms, weights);

            error_estimator->estimate_error(system, error);

            // Print out status at each adaptive step.
            Real global_error = error.l2_norm();
            std::cout << "Adaptive step " << a_step << ": " << std::endl;
            if (global_tolerance != 0.)
                std::cout << "Global_error = " << global_error
                          << std::endl;
            if (global_tolerance != 0.)
                std::cout << "Worst element error = " << error.maximum()
                          << ", mean = " << error.mean() << std::endl;

            if (global_tolerance != 0.)
            {
                // If we've reached our desired tolerance, we
                // don't need any more adaptive steps
                if (global_error < global_tolerance)
                    break;
                mesh_refinement.flag_elements_by_error_tolerance(error);
            }
            else
            {
                // If flag_elements_by_nelem_target returns true, this
                // should be our last adaptive step.
                if (mesh_refinement.flag_elements_by_nelem_target(error))
                {
                    mesh_refinement.refine_and_coarsen_elements();
                    equation_systems.reinit();
                    a_step = max_adaptivesteps;
                    break;
                }
            }

            // Carry out the adaptive mesh refinement/coarsening
            mesh_refinement.refine_and_coarsen_elements();
            equation_systems.reinit();

void assemble_and_solve(MeshBase & mesh,
                        EquationSystems & equation_systems)
{
  mesh.print_info();

  LinearImplicitSystem & system =
    equation_systems.add_system<LinearImplicitSystem> ("Poisson");

  unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);

  system.attach_assemble_function (assemble_poisson);

  // the cube has boundaries IDs 0, 1, 2, 3, 4 and 5
  std::set<boundary_id_type> boundary_ids;
  for (int j = 0; j<6; ++j)
    boundary_ids.insert(j);

  // Create a vector storing the variable numbers which the BC applies to
  std::vector<unsigned int> variables(1);
  variables[0] = u_var;

  ZeroFunction<> zf;
  DirichletBoundary dirichlet_bc(boundary_ids,
                                 variables,
                                 &zf);
  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);

  equation_systems.init();
  equation_systems.print_info();

#ifdef LIBMESH_ENABLE_AMR
  MeshRefinement mesh_refinement(mesh);

  mesh_refinement.refine_fraction()  = 0.7;
  mesh_refinement.coarsen_fraction() = 0.3;
  mesh_refinement.max_h_level()      = 5;

  const unsigned int max_r_steps = 2;

  for (unsigned int r_step=0; r_step<=max_r_steps; r_step++)
    {
      system.solve();
      if (r_step != max_r_steps)
        {
          ErrorVector error;
          KellyErrorEstimator error_estimator;

          error_estimator.estimate_error(system, error);

          libMesh::out << "Error estimate\nl2 norm = "
                       << error.l2_norm()
                       << "\nmaximum = "
                       << error.maximum()
                       << std::endl;

          mesh_refinement.flag_elements_by_error_fraction (error);

          mesh_refinement.refine_and_coarsen_elements();

          equation_systems.reinit();
        }
    }
#else
  system.solve();
#endif
}

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

              // The lid-driven cavity problem isn't in H1, so
              // lets estimate L2 error
              u->error_norm = L2;

              error_estimator.reset(u);
            }
          else
            {
              // If we aren't adapting to a tolerance we need a
              // target mesh size
              libmesh_assert_greater (nelem_target, 0);

              // Kelly is a lousy estimator to use for a problem
              // not in H1 - if we were doing more than a few
              // timesteps we'd need to turn off or limit the
              // maximum level of our adaptivity eventually
              error_estimator.reset(new KellyErrorEstimator);
            }

          // Calculate error based on u and v (and w?) but not p
	  std::vector<Real> weights(2,1.0);  // u, v
          if (dim == 3)
            weights.push_back(1.0);          // w
          weights.push_back(0.0);            // p
	  // Keep the same default norm type.
	  std::vector<FEMNormType>
	    norms(1, error_estimator->error_norm.type(0));
	  error_estimator->error_norm = SystemNorm(norms, weights);

          error_estimator->estimate_error(system, error);

          // Print out status at each adaptive step.
          Real global_error = error.l2_norm();
          std::cout << "Adaptive step " << a_step << ": " << std::endl;
          if (global_tolerance != 0.)
            std::cout << "Global_error = " << global_error 
                      << std::endl;
          if (global_tolerance != 0.)
            std::cout << "Worst element error = " << error.maximum()
                      << ", mean = " << error.mean() << std::endl;

          if (global_tolerance != 0.)
            {
              // If we've reached our desired tolerance, we
              // don't need any more adaptive steps
              if (global_error < global_tolerance)
                break;
              mesh_refinement.flag_elements_by_error_tolerance(error);
            }
          else
            {
              // If flag_elements_by_nelem_target returns true, this
              // should be our last adaptive step.
              if (mesh_refinement.flag_elements_by_nelem_target(error))
                {
                  mesh_refinement.refine_and_coarsen_elements();
                  equation_systems.reinit();
                  a_step = max_adaptivesteps;
                  break;
                }
            }

          // Carry out the adaptive mesh refinement/coarsening
          mesh_refinement.refine_and_coarsen_elements();
          equation_systems.reinit();

//.........这里部分代码省略.........
      // To solve to a tolerance in this problem we
      // need a better estimator than Kelly
      if (global_tolerance != 0.)
        {
          // We can't adapt to both a tolerance and a mesh
          // size at once
          libmesh_assert_equal_to (nelem_target, 0);

          UniformRefinementEstimator * u =
            new UniformRefinementEstimator;

          // The lid-driven cavity problem isn't in H1, so
          // lets estimate L2 error
          u->error_norm = L2;

          error_estimator.reset(u);
        }
      else
        {
          // If we aren't adapting to a tolerance we need a
          // target mesh size
          libmesh_assert_greater (nelem_target, 0);

          // Kelly is a lousy estimator to use for a problem
          // not in H1 - if we were doing more than a few
          // timesteps we'd need to turn off or limit the
          // maximum level of our adaptivity eventually
          error_estimator.reset(new KellyErrorEstimator);
        }

      error_estimator->estimate_error(system, error);

      // Print out status at each adaptive step.
      Real global_error = error.l2_norm();
      libMesh::out << "Adaptive step "
                   << a_step
                   << ": "
                   << std::endl;

      if (global_tolerance != 0.)
        libMesh::out << "Global_error = "
                     << global_error
                     << std::endl;

      if (global_tolerance != 0.)
        libMesh::out << "Worst element error = "
                     << error.maximum()
                     << ", mean = "
                     << error.mean()
                     << std::endl;

      if (global_tolerance != 0.)
        {
          // If we've reached our desired tolerance, we
          // don't need any more adaptive steps
          if (global_error < global_tolerance)
            break;
          mesh_refinement.flag_elements_by_error_tolerance(error);
        }
      else
        {
          // If flag_elements_by_nelem_target returns true, this
          // should be our last adaptive step.
          if (mesh_refinement.flag_elements_by_nelem_target(error))
            {
              mesh_refinement.refine_and_coarsen_elements();

//.........这里部分代码省略.........
            UniformRefinementEstimator *u =
              new UniformRefinementEstimator;

            // The lid-driven cavity problem isn't in H1, so
            // lets estimate L2 error
            u->error_norm = L2;

            error_estimator.reset(u);
          }
        else
          {
            // If we aren't adapting to a tolerance we need a
            // target mesh size
            libmesh_assert_greater (nelem_target, 0);

            // Kelly is a lousy estimator to use for a problem
            // not in H1 - if we were doing more than a few
            // timesteps we'd need to turn off or limit the
            // maximum level of our adaptivity eventually
            error_estimator.reset(new KellyErrorEstimator);
          }

        // Calculate error
        std::vector<Real> weights(3,1.0); 

        // Keep the same default norm type.
        std::vector<FEMNormType>
          norms(1, error_estimator->error_norm.type(0));
        error_estimator->error_norm = SystemNorm(norms, weights);

        error_estimator->estimate_error(system, error);

        // Print out status at each adaptive step.
        Real global_error = error.l2_norm();
        std::cout << "Adaptive step " << a_step << ": " << std::endl;
        if (global_tolerance != 0.)
          std::cout << "Global_error = " << global_error
                    << std::endl;
        if (global_tolerance != 0.)
          std::cout << "Worst element error = " << error.maximum()
                    << ", mean = " << error.mean() << std::endl;

        if (global_tolerance != 0.)
          {
            // If we've reached our desired tolerance, we
            // don't need any more adaptive steps
            if (global_error < global_tolerance)
              break;
            mesh_refinement.flag_elements_by_error_tolerance(error);
          }
        else
          {
            // If flag_elements_by_nelem_target returns true, this
            // should be our last adaptive step.
            if (mesh_refinement.flag_elements_by_nelem_target(error))
              {
                mesh_refinement.refine_and_coarsen_elements();
                equation_systems.reinit();
                a_step = max_adaptivesteps;
                break;
              }
          }

        // Carry out the adaptive mesh refinement/coarsening
        mesh_refinement.refine_and_coarsen_elements();
        equation_systems.reinit();