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

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

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

            std::cout << "Refined mesh to "
                      << mesh.n_active_elem()
                      << " active elements and "
                      << equation_systems.n_active_dofs()
                      << " active dofs." << std::endl;
        } // End loop over adaptive steps

//.........这里部分代码省略.........
      // 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
}

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

          std::cout << "Refined mesh to "
                    << mesh.n_active_elem()
                    << " active elements and "
                    << equation_systems.n_active_dofs()
                    << " active dofs." << std::endl;
        }

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

          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();
              equation_systems.reinit();
              a_step = max_adaptivesteps;
              break;
            }
        }

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

      libMesh::out << "Refined mesh to "
                   << mesh.n_active_elem()
                   << " active elements and "
                   << equation_systems.n_active_dofs()

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

        std::cout << "Refined mesh to "
                  << mesh.n_active_elem()
                  << " active elements and "
                  << equation_systems.n_active_dofs()
                  << " active dofs." << std::endl;
      }