SBMUtils Namespace Reference

Classes
struct	ClassificationSubdomains
	The subdomain IDs a (possibly partial) element can be labeled with. More...
struct	ElementDomainOccupancy
	Results of testing an element against a domain. More...
struct	InterceptedSubdomainPolicy
	Policy for assigning partial elements to a dedicated intercepted subdomain. More...
struct	SubdomainOccupancy
	A candidate subdomain and its occupancy of an element. More...

Functions
libMesh::Point	distanceFrom (const SurfaceElement &surface_elem, const libMesh::Point &pt)
	Returns the vector from pt to the nearest point on the surface element: the normal projection if it falls inside the element, otherwise the nearest edge or vertex.
bool	checkWatertightnessFromRawElems (const std::vector< const Elem * > &bd_elements)
Real	activeElementFraction (const Elem &elem, Order qrule_order, const std::function< bool(const libMesh::Point &)> &is_active)
	Compute the fraction of an element's quadrature-weighted measure satisfying a predicate.
bool	isInactive (Real active_fraction, Real lambda)
	Return whether a partial element is inactive, i.e.
ElementDomainOccupancy	elementDomainOccupancy (const Elem &elem, Order qrule_order, const std::function< bool(const libMesh::Point &)> &is_in_domain)
	Test an element against the domain defined by is_in_domain.
ElementDomainOccupancy	elementDomainOccupancy (const Elem &elem, Order qrule_order, const std::function< bool(const libMesh::Point &)> &is_in_domain, const std::function< bool(const libMesh::Point &)> &is_outside_domain)
	As above, but use is_outside_domain to test nodes separately.
std::optional< SubdomainID >	selectSubdomainFromOccupancies (const std::vector< SubdomainOccupancy > &candidate_occupancies, const InterceptedSubdomainPolicy &intercepted_policy, Real lambda)
	Select a subdomain based on how multiple candidate domains occupy an element.
SubdomainID	classifySubdomainFromOccupancy (const ElementDomainOccupancy &occupancy, const ClassificationSubdomains &subdomain_id_settings, bool mark_intercepted, Real lambda)
	Classify an element from its retained-domain occupancy measurements.
std::vector< const Function * >	buildDistanceFunctions (const std::vector< FunctionName > &function_names, const FunctionInterface &function_provider)
	Build a list of distance functions based on names specified in input.
RealVectorValue	distanceVectorFromFunction (const Function *func, const libMesh::Point &pt, Real t)
	Compute the distance vector induced by a distance function.
RealVectorValue	trueNormalFromFunction (const Function *func, const libMesh::Point &pt, Real t)
	Compute the true boundary surface normal at the point on the boundary closest to pt.
RealVectorValue	closestDistanceVector (const std::vector< const Function * > &funcs, const libMesh::Point &pt, Real t)
	Scan all distance functions and return the closest distance vector.
RealVectorValue	closestTrueNormalVector (const std::vector< const Function * > &funcs, const libMesh::Point &pt, Real t)
	Scan all distance functions and return the corresponding normal vector.
Real	unionSignedDistance (const std::vector< const Function * > &funcs, Real t, const Point &p)
	Computes the union signed distance by taking the minimum of all signed distance functions.
Point	distanceFrom (const SurfaceElement &surface_elem, const Point &pt)

Function Documentation

activeElementFraction()

Real SBMUtils::activeElementFraction	(	const Elem &	elem,
		Order	qrule_order,
		const std::function< bool(const libMesh::Point &)> &	is_active
	)

Compute the fraction of an element's quadrature-weighted measure satisfying a predicate.

Definition at line 26 of file SBMUtils.C.

{
  const FEType fe_type(elem.default_order(), LAGRANGE);
  auto fe = FEBase::build(elem.dim(), fe_type);
  QGauss qrule(elem.dim(), qrule_order);
  const auto & q_points = fe->get_xyz();
  const auto & JxW = fe->get_JxW();
  fe->attach_quadrature_rule(&qrule);
  fe->reinit(&elem);
 
  Real active_measure = 0.0;
  Real total_measure = 0.0;
  for (const auto i : index_range(q_points))
  {
    if (is_active(q_points[i]))
      active_measure += JxW[i];
    total_measure += JxW[i];
  }
 
  return active_measure / total_measure;
}

Referenced by elementDomainOccupancy(), and TEST().

buildDistanceFunctions()

std::vector< const Function * > SBMUtils::buildDistanceFunctions	(	const std::vector< FunctionName > &	function_names,
		const FunctionInterface &	function_provider
	)

Build a list of distance functions based on names specified in input.

Definition at line 172 of file SBMUtils.C.

{
  std::vector<const Function *> funcs;
  funcs.reserve(function_names.size());
 
  for (const auto & name : function_names)
  {
    const Function * func = &function_provider.getFunctionByName(name);
    if (!dynamic_cast<const MooseParsedFunction *>(func) &&
        !dynamic_cast<const UnsignedDistanceToSurfaceMesh *>(func) &&
        !dynamic_cast<const SignedDistanceToSurfaceMesh *>(func))
    {
      mooseError("SBM distance helpers only support ParsedFunction, "
                 "UnsignedDistanceToSurfaceMesh, or SignedDistanceToSurfaceMesh types. Offending "
                 "function: ",
                 name);
    }
    funcs.emplace_back(func);
  }
 
  return funcs;
}

Referenced by BoundaryShortestDistanceToSurface::BoundaryShortestDistanceToSurface(), and ShortestDistanceToSurface::ShortestDistanceToSurface().

checkWatertightnessFromRawElems()

bool SBMUtils::checkWatertightnessFromRawElems

(

const std::vector< const Elem * > &

bd_elements

)

Definition at line 161 of file SBMUtils.C.

{
  for (const auto * el : bd_elements)
    for (unsigned int s = 0; s < el->n_sides(); ++s)
      if (!el->neighbor_ptr(s))
        return false;
 
  return true;
}

classifySubdomainFromOccupancy()

SubdomainID SBMUtils::classifySubdomainFromOccupancy	(	const ElementDomainOccupancy &	occupancy,
		const ClassificationSubdomains &	subdomain_id_settings,
		bool	mark_intercepted,
		Real	lambda
	)

Classify an element from its retained-domain occupancy measurements.

An element whose nodes and measure are entirely in the domain is inside; one whose nodes and measure are entirely outside the domain is outside. A remaining (partial) element is assigned the intercepted subdomain when mark_intercepted is true, and otherwise is resolved by the lambda threshold (see isInactive).

Definition at line 138 of file SBMUtils.C.

{
  // Same-side nodes do not rule out a surface crossing the element or enclosing a region
  // within it. Quadrature sampling detects most such cases, but may miss very small regions.
  // Exact endpoint comparisons are intentional: activeElementFraction returns exactly zero
  // or one when no or all quadrature points are in the retained domain, respectively.
  if (occupancy.all_nodes_in_domain && occupancy.fraction_of_element_in_domain == 1.0)
    return subdomain_id_settings.inside;
 
  if (occupancy.all_nodes_outside_domain && occupancy.fraction_of_element_in_domain == 0.0)
    return subdomain_id_settings.outside;
 
  if (mark_intercepted)
    return subdomain_id_settings.intercepted;
 
  return isInactive(occupancy.fraction_of_element_in_domain, lambda) ? subdomain_id_settings.outside
                                                                     : subdomain_id_settings.inside;
}

Referenced by InterceptedElementModifier::computeSubdomainID(), and SubdomainInterceptedGenerator::generate().

closestDistanceVector()

RealVectorValue SBMUtils::closestDistanceVector	(	const std::vector< const Function * > &	funcs,
		const libMesh::Point &	pt,
		Real	t
	)

Scan all distance functions and return the closest distance vector.

Definition at line 240 of file SBMUtils.C.

{
  Real min_dist = std::numeric_limits<Real>::max();
  RealVectorValue closest_dist_vec;
 
  for (const auto & func : funcs)
  {
    const auto dist_vec = distanceVectorFromFunction(func, pt, t);
    const auto dist = dist_vec.norm();
    if (dist < min_dist)
    {
      min_dist = dist;
      closest_dist_vec = dist_vec;
    }
  }
 
  return closest_dist_vec;
}

Referenced by ShortestDistanceToSurface::distanceVector(), and TEST_F().

closestTrueNormalVector()

RealVectorValue SBMUtils::closestTrueNormalVector	(	const std::vector< const Function * > &	funcs,
		const libMesh::Point &	pt,
		Real	t
	)

Scan all distance functions and return the corresponding normal vector.

Definition at line 262 of file SBMUtils.C.

{
  Real min_dist = std::numeric_limits<Real>::max();
  RealVectorValue closest_normal_vec;
 
  for (const auto & func : funcs)
  {
    const auto dist_vec = distanceVectorFromFunction(func, pt, t);
    const auto dist = dist_vec.norm();
    if (dist < min_dist)
    {
      min_dist = dist;
      closest_normal_vec = trueNormalFromFunction(func, pt, t);
    }
  }
 
  return closest_normal_vec;
}

Referenced by TEST_F(), and ShortestDistanceToSurface::trueNormal().

distanceFrom() [1/2]

libMesh::Point SBMUtils::distanceFrom	(	const SurfaceElement &	surface_elem,
		const libMesh::Point &	pt
	)

Returns the vector from pt to the nearest point on the surface element: the normal projection if it falls inside the element, otherwise the nearest edge or vertex.

Precondition (mooseAssert in debug builds): the element's sides are EDGE2 or NODEELEM (the only side types this routine handles).

Referenced by UnsignedDistanceToSurfaceMesh::distanceVectorToSurface(), SBMInterfaceManager::queryInterface(), TEST(), TEST(), TEST(), TEST(), and TEST().

distanceFrom() [2/2]

Point SBMUtils::distanceFrom	(	const SurfaceElement &	surface_elem,
		const Point &	pt
	)

Definition at line 23 of file SBMSurfaceDistance.C.

{
  const Elem & element = surface_elem.elem();
  const Point & normal = surface_elem.normal();
 
  // Precondition: this routine only handles EDGE2 and NODEELEM element sides.
  // libMesh standard elements have uniform side types, so probing side(0) is
  // sufficient. Validated once here as a debug assert; callers build these from
  // supported surface meshes.
  mooseAssert(
      [&element]()
      {
        if (element.n_sides() == 0)
          return true;
        const auto t = element.build_side_ptr(0)->type();
        return t == EDGE2 || t == NODEELEM;
      }(),
      "SBMUtils::distanceFrom only handles EDGE2 and NODEELEM element sides.");
 
  // (a) Project pt onto the normal direction
  const auto vec_to_first = element.point(0) - pt;
  const auto scale = vec_to_first * normal;
  const auto projection = normal * scale;
 
  // Check if projection point lands inside the geometry
  if (element.contains_point(pt + projection))
    return projection;
 
  // (b) Point to closest edge or node
  Real min_dist = std::numeric_limits<Real>::max();
  Point closest_vec;
 
  const unsigned int n_edges = element.n_sides();
  for (unsigned int j = 0; j < n_edges; ++j)
  {
    std::unique_ptr<const Elem> curr_edge = element.build_side_ptr(j);
 
    switch (curr_edge->type())
    {
      case EDGE2:
      {
        const Point & p1 = *curr_edge->node_ptr(0);
        const Point & p2 = *curr_edge->node_ptr(1);
 
        const Point edge = p2 - p1;
        Real t = ((pt - p1) * edge) / (edge * edge);
        t = std::clamp(t, 0.0, 1.0);
        const Point proj = p1 + t * edge;
        const Real dist = (pt - proj).norm();
 
        if (dist < min_dist)
        {
          min_dist = dist;
          closest_vec = proj - pt;
        }
        break;
      }
 
      case NODEELEM:
      {
        const Point & p = *curr_edge->node_ptr(0);
        const Real dist = (pt - p).norm();
        if (dist < min_dist)
        {
          min_dist = dist;
          closest_vec = p - pt;
        }
        break;
      }
 
      default:
        mooseAssert(false, "unreachable: side type validated by the precondition above");
    }
  }
 
  return closest_vec;
}

distanceVectorFromFunction()

RealVectorValue SBMUtils::distanceVectorFromFunction	(	const Function *	func,
		const libMesh::Point &	pt,
		Real	t
	)

Compute the distance vector induced by a distance function.

The provided function must represent a distance to the true boundary: either a signed distance function (ParsedFunction) or an unsigned distance function provided by UnsignedDistanceToSurfaceMesh.

Definition at line 197 of file SBMUtils.C.

{
  mooseAssert(dynamic_cast<const MooseParsedFunction *>(func) ||
                  dynamic_cast<const UnsignedDistanceToSurfaceMesh *>(func) ||
                  dynamic_cast<const SignedDistanceToSurfaceMesh *>(func),
              "Function was not a valid distance strategy, the only "
              "supported types are ParsedFunction, UnsignedDistanceToSurfaceMesh, or "
              "SignedDistanceToSurfaceMesh.");
 
  const Real phi = func->value(t, pt);
  const RealVectorValue grad_phi = func->gradient(t, pt);
  const Real grad_norm = grad_phi.norm();
 
  if (grad_norm <= libMesh::TOLERANCE)
    return RealVectorValue(0.0, 0.0, 0.0);
 
  return -(phi / grad_norm) * grad_phi;
}

Referenced by closestDistanceVector(), closestTrueNormalVector(), ShortestDistanceToSurface::distanceVectorByFunc(), ShortestDistanceToSurface::distanceVectorByIndex(), BoundaryShortestDistanceToSurface::execute(), TEST_F(), TEST_F(), and trueNormalFromFunction().

elementDomainOccupancy() [1/2]

ElementDomainOccupancy SBMUtils::elementDomainOccupancy	(	const Elem &	elem,
		Order	qrule_order,
		const std::function< bool(const libMesh::Point &)> &	is_in_domain
	)

Test an element against the domain defined by is_in_domain.

The result says whether all nodes are in or out of the domain and estimates the fraction of the element inside it using quadrature points. is_in_domain decides if a point belongs to the domain.

Definition at line 51 of file SBMUtils.C.

{
  const auto is_outside_domain = [&is_in_domain](const Point & point)
  { return !is_in_domain(point); };
  return elementDomainOccupancy(elem, qrule_order, is_in_domain, is_outside_domain);
}

Referenced by InterceptedElementModifier::computeSubdomainID(), SubdomainElementModifier::computeSubdomainID(), elementDomainOccupancy(), SubdomainGrainIDGenerator::generate(), SubdomainInterceptedGenerator::generate(), TEST(), and TEST().

elementDomainOccupancy() [2/2]

ElementDomainOccupancy SBMUtils::elementDomainOccupancy	(	const Elem &	elem,
		Order	qrule_order,
		const std::function< bool(const libMesh::Point &)> &	is_in_domain,
		const std::function< bool(const libMesh::Point &)> &	is_outside_domain
	)

As above, but use is_outside_domain to test nodes separately.

This allows boundary nodes to be neither inside nor outside the domain.

Definition at line 61 of file SBMUtils.C.

{
  bool all_nodes_in_domain = true;
  bool all_nodes_outside_domain = true;
  for (const auto node : make_range(elem.n_nodes()))
  {
    const auto & point = elem.point(node);
    all_nodes_in_domain = all_nodes_in_domain && is_in_domain(point);
    all_nodes_outside_domain = all_nodes_outside_domain && is_outside_domain(point);
  }
 
  return {all_nodes_in_domain,
          all_nodes_outside_domain,
          activeElementFraction(elem, qrule_order, is_in_domain)};
}

isInactive()

bool SBMUtils::isInactive	(	Real	active_fraction,
		Real	lambda
	)

Return whether a partial element is inactive, i.e.

its inactive fraction exceeds lambda.

The endpoint handling preserves the convention that lambda zero rejects and lambda one accepts a partially active element.

Definition at line 81 of file SBMUtils.C.

{
  if (MooseUtils::absoluteFuzzyEqual(lambda, 0))
    return true;
  if (MooseUtils::absoluteFuzzyEqual(lambda, 1))
    return false;
 
  return MooseUtils::absoluteFuzzyGreaterThan(1.0 - active_fraction, lambda);
}

Referenced by classifySubdomainFromOccupancy(), selectSubdomainFromOccupancies(), and TEST().

selectSubdomainFromOccupancies()

std::optional< SubdomainID > SBMUtils::selectSubdomainFromOccupancies	(	const std::vector< SubdomainOccupancy > &	candidate_occupancies,
		const InterceptedSubdomainPolicy &	intercepted_policy,
		Real	lambda
	)

Select a subdomain based on how multiple candidate domains occupy an element.

A fully occupied candidate takes precedence. Otherwise, the candidate with the largest domain fraction is selected, with fuzzy ties resolved by the lowest subdomain ID. A partial candidate is accepted according to the intercepted-subdomain policy and lambda.

Parameters

candidate_occupancies	Candidate subdomain IDs and their measured element occupancies
intercepted_policy	Whether partial elements receive a dedicated subdomain ID, and that ID
lambda	Maximum fraction outside the selected domain for accepting a partial element

Returns

The selected subdomain ID, or std::nullopt if no candidate subdomain is selected

Definition at line 92 of file SBMUtils.C.

{
  std::optional<SubdomainID> fully_inside_subdomain;
  std::optional<SubdomainID> best_partial_subdomain;
  Real max_fraction = 0.0;
 
  for (const auto & candidate : candidate_occupancies)
  {
    const auto & occupancy = candidate.occupancy;
    if (occupancy.all_nodes_in_domain && occupancy.fraction_of_element_in_domain == 1.0)
    {
      if (!fully_inside_subdomain || candidate.subdomain_id < *fully_inside_subdomain)
        fully_inside_subdomain = candidate.subdomain_id;
    }
    else if ((!occupancy.all_nodes_outside_domain ||
              occupancy.fraction_of_element_in_domain > 0.0) &&
             occupancy.fraction_of_element_in_domain >= max_fraction)
      max_fraction = occupancy.fraction_of_element_in_domain;
  }
 
  if (fully_inside_subdomain)
    return fully_inside_subdomain;
 
  for (const auto & candidate : candidate_occupancies)
    if ((!candidate.occupancy.all_nodes_outside_domain ||
         candidate.occupancy.fraction_of_element_in_domain > 0.0) &&
        MooseUtils::absoluteFuzzyEqual(candidate.occupancy.fraction_of_element_in_domain,
                                       max_fraction) &&
        (!best_partial_subdomain || candidate.subdomain_id < *best_partial_subdomain))
      best_partial_subdomain = candidate.subdomain_id;
 
  if (!best_partial_subdomain)
    return std::nullopt;
 
  if (intercepted_policy.mark_intercepted)
    return intercepted_policy.subdomain_id;
 
  if (isInactive(max_fraction, lambda))
    return std::nullopt;
 
  return best_partial_subdomain;
}

Referenced by SubdomainElementModifier::computeSubdomainID(), SubdomainGrainIDGenerator::generate(), TEST(), TEST(), and TEST().

trueNormalFromFunction()

RealVectorValue SBMUtils::trueNormalFromFunction	(	const Function *	func,
		const libMesh::Point &	pt,
		Real	t
	)

Compute the true boundary surface normal at the point on the boundary closest to pt.

Definition at line 217 of file SBMUtils.C.

{
  if (const auto * parsed = dynamic_cast<const MooseParsedFunction *>(func))
  {
    const auto proj_pt = pt + distanceVectorFromFunction(func, pt, t);
    const RealVectorValue grad_phi = parsed->gradient(t, proj_pt);
    const Real grad_norm = grad_phi.norm();
    if (grad_norm <= libMesh::TOLERANCE)
      return RealVectorValue(0.0, 0.0, 0.0);
    return grad_phi / grad_norm;
  }
  else
  {
    const auto * mesh_func = dynamic_cast<const UnsignedDistanceToSurfaceMesh *>(func);
    if (!mesh_func)
      mesh_func = dynamic_cast<const SignedDistanceToSurfaceMesh *>(func);
 
    mooseAssert(mesh_func, "Function was not a valid distance strategy");
    return mesh_func->surfaceNormal(pt);
  }
}

Referenced by closestTrueNormalVector(), BoundaryShortestDistanceToSurface::execute(), TEST_F(), ShortestDistanceToSurface::trueNormalByFunc(), and ShortestDistanceToSurface::trueNormalByIndex().

unionSignedDistance()

Real SBMUtils::unionSignedDistance	(	const std::vector< const Function * > &	funcs,
		Real	t,
		const Point &	p
	)

Computes the union signed distance by taking the minimum of all signed distance functions.

This definition is valid for overlapping geometries and preserves the correct union boundary.

Definition at line 284 of file SBMUtils.C.

{
  // Ensure all distance functions are valid signed distance strategies
  for (const auto * func : funcs)
  {
    if (!dynamic_cast<const MooseParsedFunction *>(func) &&
        !dynamic_cast<const SignedDistanceToSurfaceMesh *>(func))
      mooseError("Signed distance requested but function was not a valid signed distance strategy. "
                 "Valid types are MooseParsedFunction or SignedDistanceToSurfaceMesh. "
                 "Offending function: ",
                 func->name());
  }
 
  mooseAssert(!funcs.empty(), "unionSignedDistance requires at least one function.");
 
  // Union signed distance: min of all signed distances
  Real min_value = funcs[0]->value(t, p);
  for (const auto i : make_range(std::size_t(1), funcs.size()))
  {
    const Real val = funcs[i]->value(t, p);
    if (val < min_value)
      min_value = val;
  }
 
  return min_value;
}

Referenced by TEST_F(), TEST_F(), and ShortestDistanceToSurface::value().