libMesh::SmoothnessEstimator::EstimateSmoothness Class Reference

Class to compute the error contribution for a range of elements. More...

Public Member Functions
	EstimateSmoothness (const System &sys, const SmoothnessEstimator &ee, ErrorVector &epc)
void	operator() (const ConstElemRange &range) const

Private Attributes
const System &	system
const SmoothnessEstimator &	smoothness_estimator
ErrorVector &	smoothness_per_cell

Detailed Description

Class to compute the error contribution for a range of elements.

May be executed in parallel on separate threads.

Definition at line 118 of file smoothness_estimator.h.

Constructor & Destructor Documentation

EstimateSmoothness()

libMesh::SmoothnessEstimator::EstimateSmoothness::EstimateSmoothness ( const System &  sys, const SmoothnessEstimator &  ee, ErrorVector &  epc  )

inline

Definition at line 121 of file smoothness_estimator.h.

                                      :
      system(sys),
      smoothness_estimator(ee),
      smoothness_per_cell(epc)
    {}

Member Function Documentation

operator()()

void libMesh::SmoothnessEstimator::EstimateSmoothness::operator()

(

const ConstElemRange &

range

)

const

Definition at line 197 of file smoothness_estimator.C.

{
  // The current mesh
  const MeshBase & mesh = system.get_mesh();
 
  // The dimensionality of the mesh
  const unsigned int dim = mesh.mesh_dimension();
 
  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();
 
  // The DofMap for this system
  const DofMap & dof_map = system.get_dof_map();
 
  //------------------------------------------------------------
  // Iterate over all the elements in the range.
for (const auto & elem : range)
{
  const dof_id_type e_id = elem->id();
  // Loop over each variable
  for (unsigned int var=0; var<n_vars; var++)
  {
    const FEType & fe_type = dof_map.variable_type(var);
    const Order element_order = fe_type.order + elem->p_level();
 
    std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
    std::unique_ptr<QBase> qrule(fe_type.default_quadrature_rule(dim, smoothness_estimator._extra_order));
    fe->attach_quadrature_rule(qrule.get());
 
    const std::vector<Real> & JxW = fe->get_JxW();
    const std::vector<Point> & q_point = fe->get_xyz();
    const std::vector<std::vector<Real>> * phi = &(fe->get_phi());
 
    std::vector<dof_id_type> dof_indices;
 
    unsigned int matsize = element_order + 1;
    if (dim > 1)
    {
      matsize *= (element_order + 2);
      matsize /= 2;
    }
    if (dim > 2)
    {
      matsize *= (element_order + 3);
      matsize /= 3;
    }
 
    DenseMatrix<Number> Kp(matsize, matsize);
    DenseVector<Number> F;
    DenseVector<Number> Pu_h;
 
    F.resize(matsize);
    Pu_h.resize(matsize);
 
 
    fe->reinit(elem);
 
    dof_map.dof_indices(elem, dof_indices, var);
    libmesh_assert_equal_to(dof_indices.size(), phi->size());
 
    const unsigned int n_dofs = cast_int<unsigned int>(dof_indices.size());
    const unsigned int n_qp   = qrule->n_points();
 
    for (unsigned int qp = 0; qp < n_qp; qp++)
    {
      std::vector<Real> psi = legepoly(dim, element_order, q_point[qp], matsize);
      const unsigned int psi_size = cast_int<unsigned int>(psi.size());
 
      for (unsigned int i = 0; i < matsize; i++)
        for (unsigned int j = 0; j < matsize; j++)
          Kp(i, j) += JxW[qp] * psi[i] * psi[j];
 
      Number u_h = libMesh::zero;
      for (unsigned int i = 0; i < n_dofs; i++)
        u_h += (*phi)[i][qp] * system.current_solution(dof_indices[i]);
 
      for (unsigned int i = 0; i < psi_size; i++)
        F(i) += JxW[qp] * u_h * psi[i];
    }
 
    Kp.lu_solve(F, Pu_h);
 
    // Generate index to total degree map. Total_degree_per_index[i] gives the degree ith basis function
    std::vector<unsigned int> total_degree_per_index;
 
    for (unsigned int poly_deg = 0; poly_deg <= static_cast<unsigned int>(element_order); ++poly_deg)
    {
      switch (dim)
      {
        case 3:
          for (int i = poly_deg; i >= 0; --i)
            for (int j = poly_deg - i; j >= 0; --j)
            {
              int k = poly_deg - i - j;
              total_degree_per_index.push_back(i + j + k);
            }
          break;
 
        case 2:
          for (int i = poly_deg; i >= 0; --i)
          {
            int j = poly_deg - i;
            total_degree_per_index.push_back(i + j);
          }
          break;
 
        case 1:
          total_degree_per_index.push_back(poly_deg);
          break;
 
        default:
          libmesh_error_msg("Invalid dimension dim " << dim);
      }
    }
 
    // Group coefficients |c_k| by degree k
    std::map<unsigned int, std::vector<Number>> coeff_by_degree;
 
    for (unsigned int i = 0; i < Pu_h.size(); ++i)
    {
      unsigned int degree = total_degree_per_index[i];
      // Excluding the constant mode i.e, zeroth order coefficient as we use ln|order| later
      // Constant order term often dominates the scale and doesn't inform about regularity.
      if (degree > 0)
        coeff_by_degree[degree].push_back(Pu_h(i));
    }
 
    // Compute L2 norm of each group |c_k|
    std::vector<unsigned int> degrees;
    std::vector<double> norms;
 
    for (const auto & pair : coeff_by_degree)
    {
      unsigned int deg = pair.first;
      const std::vector<Number> & coeffs = pair.second;
 
      double norm = 0.;
      for (const auto & c : coeffs)
        norm += std::norm(c);
 
      norm = std::sqrt(norm);
 
      degrees.push_back(deg);
      norms.push_back(norm);
    }
 
 
    // Collect log |c_k| and log |k|
    std::vector<double> log_norms, log_degrees;
 
    for (unsigned int i = 0; i < degrees.size(); ++i)
    {
      // filter tiny terms
      if (norms[i] > 1e-12)
      {
        log_degrees.push_back(std::log(static_cast<double>(degrees[i])));
        log_norms.push_back(std::log(norms[i]));
      }
    }
 
 
    // Fit log-log line - we use least-squares fit
    double Sx = 0., Sy = 0., Sxx = 0., Sxy = 0.;
    const size_t N = log_degrees.size();
 
    for (size_t i = 0; i < N; ++i)
    {
      Sx += log_degrees[i];
      Sy += log_norms[i];
      Sxx += log_degrees[i] * log_degrees[i];
      Sxy += log_degrees[i] * log_norms[i];
    }
 
    const double regularity = -compute_slope(N, Sx, Sy, Sxx, Sxy);
    // const double intercept = (Sy - slope * Sx) / N;
 
    Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
    smoothness_per_cell[e_id] = regularity;
  } // end variables loop
 
} // end element loop
 
} // End () operator definition

References libMesh::SmoothnessEstimator::_extra_order, libMesh::FEGenericBase< OutputType >::build(), libMesh::SmoothnessEstimator::compute_slope(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), dim, libMesh::DofMap::dof_indices(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::SmoothnessEstimator::legepoly(), libMesh::DenseMatrix< T >::lu_solve(), mesh, libMesh::MeshBase::mesh_dimension(), n_vars, libMesh::System::n_vars(), std::norm(), libMesh::FEType::order, libMesh::DenseVector< T >::resize(), libMesh::DenseVector< T >::size(), smoothness_estimator, smoothness_per_cell, libMesh::Threads::spin_mtx, system, libMesh::DofMap::variable_type(), and libMesh::zero.

Member Data Documentation

smoothness_estimator

const SmoothnessEstimator& libMesh::SmoothnessEstimator::EstimateSmoothness::smoothness_estimator

private

Definition at line 133 of file smoothness_estimator.h.

Referenced by operator()().

smoothness_per_cell

ErrorVector& libMesh::SmoothnessEstimator::EstimateSmoothness::smoothness_per_cell

private

Definition at line 134 of file smoothness_estimator.h.

Referenced by operator()().

system

const System& libMesh::SmoothnessEstimator::EstimateSmoothness::system

private

Definition at line 132 of file smoothness_estimator.h.

Referenced by operator()().

---

The documentation for this class was generated from the following files:

• include/error_estimation/smoothness_estimator.h
• src/error_estimation/smoothness_estimator.C