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

void MeshRefinement::flag_elements_by_mean_stddev (const ErrorVector & error_per_cell,
                                                   const Real refine_frac,
                                                   const Real coarsen_frac,
                                                   const unsigned int max_l)
{
  // The function arguments are currently just there for
  // backwards_compatibility
  if (!_use_member_parameters)
    {
      // If the user used non-default parameters, lets warn
      // that they're deprecated
      if (refine_frac != 0.3 ||
          coarsen_frac != 0.0 ||
          max_l != libMesh::invalid_uint)
        libmesh_deprecated();

      _refine_fraction = refine_frac;
      _coarsen_fraction = coarsen_frac;
      _max_h_level = max_l;
    }

  // Get the mean value from the error vector
  const Real mean = error_per_cell.mean();

  // Get the standard deviation.  This equals the
  // square-root of the variance
  const Real stddev = std::sqrt (error_per_cell.variance());

  // Check for valid fractions
  libmesh_assert_greater_equal (_refine_fraction, 0);
  libmesh_assert_less_equal (_refine_fraction, 1);
  libmesh_assert_greater_equal (_coarsen_fraction, 0);
  libmesh_assert_less_equal (_coarsen_fraction, 1);

  // The refine and coarsen cutoff
  const Real refine_cutoff  =  mean + _refine_fraction  * stddev;
  const Real coarsen_cutoff =  std::max(mean - _coarsen_fraction * stddev, 0.);

  // Loop over the elements and flag them for coarsening or
  // refinement based on the element error
  for (auto & elem : _mesh.active_element_ptr_range())
    {
      const dof_id_type id  = elem->id();

      libmesh_assert_less (id, error_per_cell.size());

      const ErrorVectorReal elem_error = error_per_cell[id];

      // Possibly flag the element for coarsening ...
      if (elem_error <= coarsen_cutoff)
        elem->set_refinement_flag(Elem::COARSEN);

      // ... or refinement
      if ((elem_error >= refine_cutoff) && (elem->level() < _max_h_level))
        elem->set_refinement_flag(Elem::REFINE);
    }
}

//.........这里部分代码省略.........
      // Solve the system "Biharmonic", just like example 2.
      system.solve();

      std::cout << "Linear solver converged at step: "
                << system.n_linear_iterations()
                << ", final residual: "
                << system.final_linear_residual()
                << std::endl;

      // Compute the error.
      exact_sol.compute_error("Biharmonic", "u");
      // Compute the norm.
      zero_sol.compute_error("Biharmonic", "u");

      // Print out the error values
      std::cout << "L2-Norm is: "
                << zero_sol.l2_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H1-Norm is: "
                << zero_sol.h1_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H2-Norm is: "
                << zero_sol.h2_error("Biharmonic", "u")
                << std::endl
                << std::endl;
      std::cout << "L2-Error is: "
                << exact_sol.l2_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H1-Error is: "
                << exact_sol.h1_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H2-Error is: "
                << exact_sol.h2_error("Biharmonic", "u")
                << std::endl
                << std::endl;

      // Print to output file
      out << equation_systems.n_active_dofs() << " "
          << exact_sol.l2_error("Biharmonic", "u") << " "
          << exact_sol.h1_error("Biharmonic", "u") << " "
          << exact_sol.h2_error("Biharmonic", "u") << std::endl;

      // Possibly refine the mesh
      if (r_step+1 != max_r_steps)
        {
          std::cout << "  Refining the mesh..." << std::endl;

          if (uniform_refine == 0)
            {
              ErrorVector error;
              LaplacianErrorEstimator error_estimator;

              error_estimator.estimate_error(system, error);
              mesh_refinement.flag_elements_by_elem_fraction (error);

              std::cout << "Mean Error: " << error.mean() <<
                              std::endl;
              std::cout << "Error Variance: " << error.variance() <<
                              std::endl;

              mesh_refinement.refine_and_coarsen_elements();
            }
          else
            {
              mesh_refinement.uniformly_refine(1);
            }
              
          // This call reinitializes the \p EquationSystems object for
          // the newly refined mesh.  One of the steps in the
          // reinitialization is projecting the \p solution,
          // \p old_solution, etc... vectors from the old mesh to
          // the current one.
          equation_systems.reinit ();
        }
    }            
  
#ifdef LIBMESH_HAVE_EXODUS_API
  // Write out the solution
  // After solving the system write the solution
  // to a ExodusII-formatted plot file.
  ExodusII_IO (mesh).write_equation_systems (exd_file,
                                       equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API

  // Close up the output file.
  out << "];\n"
      << "hold on\n"
      << "loglog(e(:,1), e(:,2), 'bo-');\n"
      << "loglog(e(:,1), e(:,3), 'ro-');\n"
      << "loglog(e(:,1), e(:,4), 'go-');\n"
      << "xlabel('log(dofs)');\n"
      << "ylabel('log(error)');\n"
      << "title('C1 " << approx_type << " elements');\n"
      << "legend('L2-error', 'H1-error', 'H2-error');\n";
  
  // All done.  
  return 0;
#endif // #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
#endif // #ifndef LIBMESH_ENABLE_AMR
}