libMesh::SmoothnessEstimator Class Reference

This class implements the Smoothness estimate. More...

#include <smoothness_estimator.h>

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

Public Member Functions
	SmoothnessEstimator ()
	Constructor. More...
	SmoothnessEstimator (const SmoothnessEstimator &)=default
	Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this class. More...
	SmoothnessEstimator (SmoothnessEstimator &&)=default
SmoothnessEstimator &	operator= (const SmoothnessEstimator &)=default
SmoothnessEstimator &	operator= (SmoothnessEstimator &&)=default
virtual	~SmoothnessEstimator ()=default
virtual void	estimate_smoothness (const System &system, ErrorVector &smoothness_per_cell, const NumericVector< Number > *solution_vector=nullptr)
	This function uses the Legendre expansion of solution to estimate coefficient decay to quantify the solution smoothness on each cell. More...
void	extra_quadrature_order (const int extraorder)
	Increases or decreases the order of the quadrature rule used for numerical integration. More...

Protected Member Functions
void	reduce_smoothness (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm)

Static Protected Member Functions
static std::vector< Real >	legepoly (const unsigned int dim, const Order order, const Point p, const unsigned int matsize)
static Real	compute_slope (int N, Real Sx, Real Sy, Real Sxx, Real Sxy)
	Computes slop in a linear regression. More...

Protected Attributes
int	_extra_order
	Extra order to use for quadrature rule. More...

Friends
class	EstimateSmoothness

Detailed Description

This class implements the Smoothness estimate.

Author

Arjun Kulathuvayal

Varis Carey

Benjamin S. Kirk

Date

2025

Definition at line 48 of file smoothness_estimator.h.

Constructor & Destructor Documentation

SmoothnessEstimator() [1/3]

libMesh::SmoothnessEstimator::SmoothnessEstimator

(

)

Constructor.

Definition at line 44 of file smoothness_estimator.C.

                                         :
    _extra_order(1)
{
}

SmoothnessEstimator() [2/3]

libMesh::SmoothnessEstimator::SmoothnessEstimator ( const SmoothnessEstimator &  )

default

Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this class.

SmoothnessEstimator() [3/3]

libMesh::SmoothnessEstimator::SmoothnessEstimator ( SmoothnessEstimator &&  )

default

~SmoothnessEstimator()

virtual libMesh::SmoothnessEstimator::~SmoothnessEstimator ( )

virtualdefault

Member Function Documentation

compute_slope()

Real libMesh::SmoothnessEstimator::compute_slope ( int  N, Real  Sx, Real  Sy, Real  Sxx, Real  Sxy  )

staticprotected

Computes slop in a linear regression.

Definition at line 124 of file smoothness_estimator.C.

References libMesh::Real.

Referenced by libMesh::SmoothnessEstimator::EstimateSmoothness::operator()().

{
    const Real denom = (N * Sxx - Sx * Sx);
    // Triggers for first order polynomials when there only 1 point
    // available at log |c_k| and log |k| space to fit.
    if (std::abs(denom) < std::numeric_limits<Real>::epsilon()) {
        return std::numeric_limits<Real>::max();
    }
    return (N * Sxy - Sx * Sy) / denom;
}

estimate_smoothness()

void libMesh::SmoothnessEstimator::estimate_smoothness ( const System &  system, ErrorVector &  smoothness_per_cell, const NumericVector< Number > *  solution_vector = nullptr  )

virtual

This function uses the Legendre expansion of solution to estimate coefficient decay to quantify the solution smoothness on each cell.

The estimated smoothness is output in the vector (for time being, we use ErrorVector) For element order 1, the least square fit of log|order| vs log |coefficients| fails. This leads to pure h refinement. smoothness_per_cell

Definition at line 143 of file smoothness_estimator.C.

References libMesh::ParallelObject::comm(), EstimateSmoothness, libMesh::System::get_mesh(), mesh, libMesh::Threads::parallel_for(), reduce_smoothness(), libMesh::System::solution, and libMesh::NumericVector< T >::swap().

Referenced by main(), and libMesh::HPCoarsenTest::select_refinement().

{
  LOG_SCOPE("estimate_smoothness()", "SmoothnessEstimator");

  // The current mesh
  const MeshBase & mesh = system.get_mesh();

  // Resize the smoothness_per_cell vector to be
  // the number of elements, initialize it to 0.
  smoothness_per_cell.resize (mesh.max_elem_id());
  std::fill (smoothness_per_cell.begin(), smoothness_per_cell.end(), 0.);

  // Prepare current_local_solution to localize a non-standard
  // solution vector if necessary
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number> * newsol =
        const_cast<NumericVector<Number> *>(solution_vector);
      System & sys = const_cast<System &>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }

  //------------------------------------------------------------
  // Iterate over all the active elements in the mesh
  // that live on this processor.
  Threads::parallel_for (ConstElemRange(mesh.active_local_elements_begin(),
                                        mesh.active_local_elements_end(),
                                        200),
                         EstimateSmoothness(system,
                                       *this,
                                       smoothness_per_cell)
                         );

  // Each processor has now computed the smoothness for its local
  // elements, and smoothness_per_cell contains 0 for all the
  // non-local elements.  Summing the vector will provide the true
  // value for each element, local or remote
  this->reduce_smoothness (smoothness_per_cell, system.comm());

  // If we used a non-standard solution before, now is the time to fix
  // the current_local_solution
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number> * newsol = const_cast<NumericVector<Number> *>(solution_vector);
      System & sys = const_cast<System &>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }
}

extra_quadrature_order()

void libMesh::SmoothnessEstimator::extra_quadrature_order ( const int  extraorder )

inline

Increases or decreases the order of the quadrature rule used for numerical integration.

The default extraorder is 1, because properly integrating L2 error requires integrating the squares of terms with order p+1, and 2p+2 is 1 higher than what we default to using for reasonable mass matrix integration.

Definition at line 86 of file smoothness_estimator.h.

References _extra_order.

  { _extra_order = extraorder; }

legepoly()

std::vector< Real > libMesh::SmoothnessEstimator::legepoly ( const unsigned int  dim, const Order  order, const Point  p, const unsigned int  matsize  )

staticprotected

Returns

The Legendre polynomial basis function values at a point (x,y,z).

Definition at line 51 of file smoothness_estimator.C.

References dim, and libMesh::Real.

Referenced by libMesh::SmoothnessEstimator::EstimateSmoothness::operator()().

{
  std::vector<Real> psi;
  psi.reserve(matsize);

  // Evaluate 1D Legendre polynomials at x, y, z up to order
  const int npols = static_cast<int>(order) + 1;

  std::vector<Real> Lx(npols, 0.), Ly(npols, 1.), Lz(npols, 1.);
  const Real x = p(0);
  Lx[0] = 1.;
  if (npols > 1)
    Lx[1] = x;
  for (int n = 2; n < npols; ++n)
    Lx[n] = ((2. * n - 1) * x * Lx[n - 1] - (n - 1) * Lx[n - 2]) / n;

  if (dim > 1)
  {
    const Real y = p(1);
    Ly[0] = 1.;
    if (npols > 1)
      Ly[1] = y;
    for (int n = 2; n < npols; ++n)
      Ly[n] = ((2. * n - 1) * y * Ly[n - 1] - (n - 1) * Ly[n - 2]) / n;
  }

  if (dim > 2)
  {
    const Real z = p(2);
    Lz[0] = 1.;
    if (npols > 1)
      Lz[1] = z;
    for (int n = 2; n < npols; ++n)
      Lz[n] = ((2. * n - 1) * z * Lz[n - 1] - (n - 1) * Lz[n - 2]) / n;
  }

  // Loop over total degree and build tensor-product Legendre basis
  for (unsigned int poly_deg = 0; poly_deg <= static_cast<unsigned int>(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;
            psi.push_back(Lx[i] * Ly[j] * Lz[k]);
          }
        break;

      case 2:
        for (int i = poly_deg; i >= 0; --i)
        {
          int j = poly_deg - i;
          psi.push_back(Lx[i] * Ly[j]);
        }
        break;

      case 1:
        psi.push_back(Lx[poly_deg]);
        break;

      default:
        libmesh_error_msg("Invalid dimension dim " << dim);
    }
  }

  return psi;
}

operator=() [1/2]

SmoothnessEstimator& libMesh::SmoothnessEstimator::operator= ( const SmoothnessEstimator &  )

default

operator=() [2/2]

SmoothnessEstimator& libMesh::SmoothnessEstimator::operator= ( SmoothnessEstimator &&  )

default

reduce_smoothness()

void libMesh::SmoothnessEstimator::reduce_smoothness ( std::vector< ErrorVectorReal > &  error_per_cell, const Parallel::Communicator &  comm  )

protected

Definition at line 135 of file smoothness_estimator.C.

References TIMPI::Communicator::sum().

Referenced by estimate_smoothness().

{
  // Aggregates element-wise contributions computed
  // in parallel across all processors
  comm.sum(smoothness_per_cell);
}

Friends And Related Function Documentation

EstimateSmoothness

friend class EstimateSmoothness

friend

Definition at line 137 of file smoothness_estimator.h.

Referenced by estimate_smoothness().

Member Data Documentation

_extra_order

int libMesh::SmoothnessEstimator::_extra_order

protected

Extra order to use for quadrature rule.

Definition at line 102 of file smoothness_estimator.h.

Referenced by extra_quadrature_order(), and libMesh::SmoothnessEstimator::EstimateSmoothness::operator()().

---

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

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