libMesh
Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
libMesh::VariationalSmootherConstraint Class Reference

Constraint class for the VariationalMeshSmoother. More...

#include <variational_smoother_constraint.h>

Inheritance diagram for libMesh::VariationalSmootherConstraint:
[legend]

Public Member Functions

 VariationalSmootherConstraint (System &sys, const bool &preserve_subdomain_boundaries, const unsigned int verbosity)
 Constructor. More...
 
virtual ~VariationalSmootherConstraint () override
 
virtual void constrain () override
 Constraint function. More...
 

Private Member Functions

void fix_node (const Node &node)
 Constrain (i.e., fix) a node to not move during mesh smoothing. More...
 
void constrain_node_to_plane (const Node &node, const Point &ref_normal_vec)
 Constrain a node to remain in the given plane during mesh smoothing. More...
 
void constrain_node_to_line (const Node &node, const Point &line_vec)
 Constrain a node to remain on the given line during mesh smoothing. More...
 
void impose_constraint (const Node &node, const ConstraintVariant &constraint)
 Applies a given constraint to a node (e.g., fixing it, restricting it to a line or plane). More...
 

Static Private Member Functions

static bool nodes_share_boundary_id (const Node &boundary_node, const Node &neighbor_node, const Elem &containing_elem, const BoundaryInfo &boundary_info)
 Determines whether two neighboring nodes share a common boundary id. More...
 
static std::set< std::set< const Node * > > get_neighbors_for_subdomain_constraint (const MeshBase &mesh, const Node &node, const subdomain_id_type sub_id, const std::unordered_map< dof_id_type, std::vector< const Elem *>> &nodes_to_elem_map)
 Get the relevant nodal neighbors for a subdomain constraint. More...
 
static std::set< std::set< const Node * > > get_neighbors_for_boundary_constraint (const MeshBase &mesh, const Node &node, const std::unordered_set< dof_id_type > &boundary_node_ids, const BoundaryInfo &boundary_info, const std::unordered_map< dof_id_type, std::vector< const Elem *>> &nodes_to_elem_map)
 Get the relevant nodal neighbors for an external boundary constraint. More...
 
static ConstraintVariant determine_constraint (const Node &node, const unsigned int dim, const std::set< std::set< const Node *>> &side_grouped_boundary_neighbors)
 Determines the appropriate constraint (PointConstraint, LineConstraint, or PlaneConstraint) for a node based on its neighbors. More...
 

Private Attributes

const unsigned int _verbosity
 Verbosity setting. More...
 
System_sys
 
const bool _preserve_subdomain_boundaries
 Whether nodes on subdomain boundaries are subject to change via smoothing. More...
 

Detailed Description

Constraint class for the VariationalMeshSmoother.

Currently, all mesh boundary nodes are constrained to not move during smoothing. If requested (preserve_subdomain_boundaries = true), nodes on subdomain boundaries are also constrained to not move.

Definition at line 389 of file variational_smoother_constraint.h.

Constructor & Destructor Documentation

◆ VariationalSmootherConstraint()

libMesh::VariationalSmootherConstraint::VariationalSmootherConstraint ( System sys,
const bool &  preserve_subdomain_boundaries,
const unsigned int  verbosity 
)

Constructor.

Parameters
sysSystem to constrain.
preserve_subdomain_boundariesWhether to constrain nodes on subdomain boundaries to not move.

Definition at line 334 of file variational_smoother_constraint.C.

336  : Constraint(), _verbosity(verbosity), _sys(sys), _preserve_subdomain_boundaries(preserve_subdomain_boundaries) {
337 }
const unsigned int _verbosity
Verbosity setting.
const bool _preserve_subdomain_boundaries
Whether nodes on subdomain boundaries are subject to change via smoothing.

◆ ~VariationalSmootherConstraint()

libMesh::VariationalSmootherConstraint::~VariationalSmootherConstraint ( )
overridevirtualdefault

Member Function Documentation

◆ constrain()

void libMesh::VariationalSmootherConstraint::constrain ( )
overridevirtual

Constraint function.

This function will be called to constrain the system prior to a solve and must be provided by the user in a derived class.

Implements libMesh::System::Constraint.

Definition at line 341 of file variational_smoother_constraint.C.

References _preserve_subdomain_boundaries, _sys, _verbosity, libMesh::MeshTools::build_nodes_to_elem_map(), determine_constraint(), dim, libMesh::MeshTools::find_boundary_nodes(), libMesh::System::get_mesh(), get_neighbors_for_boundary_constraint(), get_neighbors_for_subdomain_constraint(), impose_constraint(), libMesh::intersect_constraints(), mesh, and libMesh::out.

342 {
343  const auto & mesh = _sys.get_mesh();
344  const auto dim = mesh.mesh_dimension();
345  const auto & proc_id = mesh.processor_id();
346 
347  // Only compute the node to elem map once
348  std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
349  MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
350 
351  const auto & boundary_info = mesh.get_boundary_info();
352 
353  const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
354 
355  // Identify/constrain subdomain boundary nodes, if requested
356  std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
358  {
359  for (const auto * elem : mesh.active_element_ptr_range())
360  {
361  const auto & sub_id1 = elem->subdomain_id();
362  for (const auto side : elem->side_index_range())
363  {
364  const auto * neighbor = elem->neighbor_ptr(side);
365  if (neighbor == nullptr)
366  continue;
367 
368  const auto & sub_id2 = neighbor->subdomain_id();
369  if (sub_id1 == sub_id2)
370  continue;
371 
372  // elem and neighbor are in different subdomains, and share nodes
373  // that need to be constrained
374  for (const auto local_node_id : elem->nodes_on_side(side))
375  {
376  const auto & node = mesh.node_ref(elem->node_id(local_node_id));
377  // Make sure we haven't already processed this node.
378  // We leave the responsibility of computing the constraint to
379  // the proc that owns the node.
380  if (subdomain_boundary_map.count(node.id()) ||
381  node.processor_id() != proc_id)
382  continue;
383 
384  // Get the relevant nodal neighbors for the subdomain constraint
385  const auto side_grouped_boundary_neighbors =
387  mesh, node, sub_id1, nodes_to_elem_map);
388 
389  // Determine which constraint should be imposed
390  const auto subdomain_constraint =
391  determine_constraint(node, dim, side_grouped_boundary_neighbors);
392 
393  // This subdomain boundary node does not lie on an external boundary,
394  // go ahead and impose constraint
395  if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
396  {
397  if (_verbosity > 20)
398  libMesh::out << "Imposing subdomain constraint on "
399  << std::endl << " " << node << std::endl
400  << " " << subdomain_constraint << std::endl;
401 
402  this->impose_constraint(node, subdomain_constraint);
403  }
404 
405  // This subdomain boundary node could lie on an external boundary, save it
406  // for later to combine with the external boundary constraint.
407  // We also save constraints for non-boundary nodes so we don't try to
408  // re-constrain the node when accessed from the neighboring elem.
409  // See subdomain_boundary_map.count call above.
410  subdomain_boundary_map[node.id()] = subdomain_constraint;
411 
412  } // for local_node_id
413 
414  } // for side
415  } // for elem
416  }
417 
418  // Loop through boundary nodes and impose constraints
419  for (const auto & bid : boundary_node_ids)
420  {
421  const auto & node = mesh.node_ref(bid);
422  // We leave the responsibility of computing the constraint to the proc
423  // that owns the node.
424  if (node.processor_id() != proc_id)
425  continue;
426 
427  // Get the relevant nodal neighbors for the boundary constraint
428  const auto side_grouped_boundary_neighbors = get_neighbors_for_boundary_constraint(
429  mesh, node, boundary_node_ids, boundary_info, nodes_to_elem_map);
430 
431  // Determine which constraint should be imposed
432  const auto boundary_constraint =
433  determine_constraint(node, dim, side_grouped_boundary_neighbors);
434 
435  // Check for the case where this boundary node is also part of a subdomain id boundary
436  if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
437  {
438  const auto & subdomain_constraint = it->second;
439  // Combine current boundary constraint with previously determined
440  // subdomain_constraint
441  auto combined_constraint = intersect_constraints(subdomain_constraint, boundary_constraint);
442 
443  // This will catch cases where constraints have no intersection
444  // Fall back to fixed node constraint
445  if (std::holds_alternative<InvalidConstraint>(combined_constraint))
446  combined_constraint = PointConstraint(node);
447 
448  if (_verbosity > 20)
449  libMesh::out << "Imposing boundary/subdomain constraint on "
450  << std::endl << " " << node << std::endl
451  << " " << combined_constraint << std::endl;
452 
453  this->impose_constraint(node, combined_constraint);
454  }
455 
456  else
457  {
458  if (_verbosity > 20)
459  libMesh::out << "Imposing boundary constraint on "
460  << std::endl << node << std::endl
461  << " " << boundary_constraint << std::endl;
462 
463  this->impose_constraint(node, boundary_constraint);
464  }
465 
466  } // end bid
467 }
unsigned int dim
MeshBase & mesh
void build_nodes_to_elem_map(const MeshBase &mesh, std::vector< std::vector< dof_id_type >> &nodes_to_elem_map)
After calling this function the input vector nodes_to_elem_map will contain the node to element conne...
Definition: mesh_tools.C:456
const MeshBase & get_mesh() const
Definition: system.h:2401
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:524
ConstraintVariant intersect_constraints(const ConstraintVariant &a, const ConstraintVariant &b)
Dispatch intersection between two constraint variants.
static std::set< std::set< const Node * > > get_neighbors_for_boundary_constraint(const MeshBase &mesh, const Node &node, const std::unordered_set< dof_id_type > &boundary_node_ids, const BoundaryInfo &boundary_info, const std::unordered_map< dof_id_type, std::vector< const Elem *>> &nodes_to_elem_map)
Get the relevant nodal neighbors for an external boundary constraint.
static std::set< std::set< const Node * > > get_neighbors_for_subdomain_constraint(const MeshBase &mesh, const Node &node, const subdomain_id_type sub_id, const std::unordered_map< dof_id_type, std::vector< const Elem *>> &nodes_to_elem_map)
Get the relevant nodal neighbors for a subdomain constraint.
const unsigned int _verbosity
Verbosity setting.
const bool _preserve_subdomain_boundaries
Whether nodes on subdomain boundaries are subject to change via smoothing.
OStreamProxy out
static ConstraintVariant determine_constraint(const Node &node, const unsigned int dim, const std::set< std::set< const Node *>> &side_grouped_boundary_neighbors)
Determines the appropriate constraint (PointConstraint, LineConstraint, or PlaneConstraint) for a nod...
void impose_constraint(const Node &node, const ConstraintVariant &constraint)
Applies a given constraint to a node (e.g., fixing it, restricting it to a line or plane)...

◆ constrain_node_to_line()

void libMesh::VariationalSmootherConstraint::constrain_node_to_line ( const Node node,
const Point line_vec 
)
private

Constrain a node to remain on the given line during mesh smoothing.

Parameters
nodeNode to constrain
line_vecvector parallel to the constraining line. This, along with the coordinates of node, are used to define the constraining line.

Definition at line 533 of file variational_smoother_constraint.C.

References _sys, libMesh::DofMap::add_constraint_row(), b, dim, distance(), libMesh::DofObject::dof_number(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::mesh_dimension(), libMesh::System::number(), and libMesh::relative_fuzzy_equals().

Referenced by impose_constraint().

535 {
536  const auto dim = _sys.get_mesh().mesh_dimension();
537 
538  // We will free the dimension most parallel to line_vec to keep the
539  // constraint coefficients small
540  const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
541  auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
542  return std::abs(a) < std::abs(b);
543  });
544  const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
545  const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
546 
547  // A line is parameterized as r(t) = node + t * line_vec, so
548  // x(t) = node(x) + t * line_vec(x)
549  // y(t) = node(y) + t * line_vec(y)
550  // z(t) = node(z) + t * line_vec(z)
551  // Let's say we leave x free. Then t = (x(t) - node(x)) / line_vec(x)
552  // Then y and z can be constrained as
553  // y = node(y) + line_vec_y * (x(t) - node(x)) / line_vec(x)
554  // = x(t) * line_vec(y) / line_vec(x) + (node(y) - node(x) * line_vec(y) / line_vec(x))
555  // z = x(t) * line_vec(z) / line_vec(x) + (node(z) - node(x) * line_vec(z) / line_vec(x))
556 
557  libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
558  for (const auto constrained_dim : make_range(dim))
559  {
560  if (constrained_dim == free_dim)
561  continue;
562 
563  DofConstraintRow constraint_row;
564  constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
565  const auto inhomogeneous_part =
566  node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
567  const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
569  constrained_dof_index, constraint_row, inhomogeneous_part, true);
570  }
571 }
unsigned int dim
const MeshBase & get_mesh() const
Definition: system.h:2401
Real distance(const Point &p)
unsigned int number() const
Definition: system.h:2393
void add_constraint_row(const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
libmesh_assert(ctx)
static const Real b
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition: dof_map.h:100
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:176
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
Definition: fuzzy_equals.h:78
const DofMap & get_dof_map() const
Definition: system.h:2417

◆ constrain_node_to_plane()

void libMesh::VariationalSmootherConstraint::constrain_node_to_plane ( const Node node,
const Point ref_normal_vec 
)
private

Constrain a node to remain in the given plane during mesh smoothing.

Parameters
nodeNode to constrain
ref_normal_vecReference normal vector to the constraining plane. This, along with the coordinates of node, are used to define the constraining plane.

Definition at line 483 of file variational_smoother_constraint.C.

References _sys, libMesh::DofMap::add_constraint_row(), dim, libMesh::DofObject::dof_number(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::mesh_dimension(), libMesh::TypeVector< T >::norm(), libMesh::System::number(), libMesh::Real, and libMesh::TOLERANCE.

Referenced by impose_constraint().

485 {
486  const auto dim = _sys.get_mesh().mesh_dimension();
487  // determine equation of plane: c_x * x + c_y * y + c_z * z + c = 0
488  std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
489  Real c = 0.;
490 
491  // We choose to constrain the dimension with the largest magnitude coefficient
492  // This approach ensures the coefficients added to the constraint_row
493  // (i.e., -c_xyz / c_max) have as small magnitude as possible
494 
495  // We initialize this to avoid maybe-uninitialized compiler error
496  unsigned int constrained_dim = 0;
497 
498  // Let's assert that we have a nonzero normal to ensure that constrained_dim
499  // is always set
500  libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
501 
502  Real max_abs_coef = 0.;
503  for (const auto d : make_range(dim))
504  {
505  const auto coef = ref_normal_vec(d);
506  xyz_coefs.push_back(coef);
507  c -= coef * node(d);
508 
509  const auto coef_abs = std::abs(coef);
510  if (coef_abs > max_abs_coef)
511  {
512  max_abs_coef = coef_abs;
513  constrained_dim = d;
514  }
515  }
516 
517  DofConstraintRow constraint_row;
518  for (const auto free_dim : make_range(dim))
519  {
520  if (free_dim == constrained_dim)
521  continue;
522 
523  const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
524  constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
525  }
526 
527  const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
528  const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
530  constrained_dof_index, constraint_row, inhomogeneous_part, true);
531 }
static constexpr Real TOLERANCE
unsigned int dim
const MeshBase & get_mesh() const
Definition: system.h:2401
unsigned int number() const
Definition: system.h:2393
void add_constraint_row(const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
libmesh_assert(ctx)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition: dof_map.h:100
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:176
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
const DofMap & get_dof_map() const
Definition: system.h:2417

◆ determine_constraint()

ConstraintVariant libMesh::VariationalSmootherConstraint::determine_constraint ( const Node node,
const unsigned int  dim,
const std::set< std::set< const Node *>> &  side_grouped_boundary_neighbors 
)
staticprivate

Determines the appropriate constraint (PointConstraint, LineConstraint, or PlaneConstraint) for a node based on its neighbors.

Parameters
nodeThe node to constrain.
dimThe mesh dimension.
Returns
The best-fit constraint for the given geometry.

Definition at line 798 of file variational_smoother_constraint.C.

References libMesh::TypeVector< T >::cross(), dim, libMesh::index_range(), libMesh::intersect_constraints(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::TypeVector< T >::norm(), libMesh::TOLERANCE, and libMesh::TypeVector< T >::unit().

Referenced by constrain().

802 {
803  // Determines the appropriate geometric constraint for a node based on its
804  // neighbors.
805 
806  // Extract neighbors in flat vector
807  std::vector<const Node *> neighbors;
808  for (const auto & side : side_grouped_boundary_neighbors)
809  neighbors.insert(neighbors.end(), side.begin(), side.end());
810 
811  // Constrain the node to it's current location
812  // Note that neighbors.size() <= 1 is used here instead of == 1 so we don't
813  // crash when trying to access *neighbors[0] a few lines below in the
814  // dim == 2 case
815  if (dim == 1 || neighbors.size() <= 1)
816  return PointConstraint{node};
817 
818  if (dim == 2 || neighbors.size() == 2)
819  {
820  // Determine whether node + all neighbors are colinear
821  bool all_colinear = true;
822  const Point ref_dir = (*neighbors[0] - node).unit();
823  for (const auto i : make_range(std::size_t(1), neighbors.size()))
824  {
825  const Point delta = *(neighbors[i]) - node;
826  libmesh_assert(delta.norm() >= TOLERANCE);
827  const Point dir = delta.unit();
828  if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
829  {
830  all_colinear = false;
831  break;
832  }
833  }
834 
835  if (all_colinear)
836  return LineConstraint{node, ref_dir};
837 
838  return PointConstraint{node};
839  }
840 
841  // dim == 3, neighbors.size() >= 3
842  std::set<PlaneConstraint> valid_planes;
843  for (const auto & side_nodes : side_grouped_boundary_neighbors)
844  {
845  std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
846  for (const auto i : index_range(side_nodes_vec))
847  {
848  const Point vec_i = (*side_nodes_vec[i] - node);
849  for (const auto j : make_range(i))
850  {
851  const Point vec_j = (*side_nodes_vec[j] - node);
852  Point candidate_normal = vec_i.cross(vec_j);
853  if (candidate_normal.norm() <= TOLERANCE)
854  continue;
855 
856  const PlaneConstraint candidate_plane{node, candidate_normal};
857  valid_planes.emplace(candidate_plane);
858  }
859  }
860  }
861 
862  // Fall back to point constraint
863  if (valid_planes.empty())
864  return PointConstraint(node);
865 
866  // Combine all the planes together to get a common constraint
867  auto it = valid_planes.begin();
868  ConstraintVariant current = *it++;
869  for (; it != valid_planes.end(); ++it)
870  {
871  current = intersect_constraints(current, *it);
872 
873  // This will catch cases where constraints have no intersection
874  // (i.e., the element surface is non-planar)
875  // Fall back to fixed node constraint
876  if (std::holds_alternative<InvalidConstraint>(current))
877  {
878  current = PointConstraint(node);
879  break;
880  }
881  }
882 
883  return current;
884 }
static constexpr Real TOLERANCE
unsigned int dim
ConstraintVariant intersect_constraints(const ConstraintVariant &a, const ConstraintVariant &b)
Dispatch intersection between two constraint variants.
libmesh_assert(ctx)
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:176
std::variant< PointConstraint, LineConstraint, PlaneConstraint, InvalidConstraint > ConstraintVariant
Type used to store a constraint that may be a PlaneConstraint, LineConstraint, or PointConstraint...
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:153

◆ fix_node()

void libMesh::VariationalSmootherConstraint::fix_node ( const Node node)
private

Constrain (i.e., fix) a node to not move during mesh smoothing.

Parameters
nodeNode to fix.

Definition at line 469 of file variational_smoother_constraint.C.

References _sys, libMesh::DofMap::add_constraint_row(), libMesh::DofObject::dof_number(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::make_range(), libMesh::MeshBase::mesh_dimension(), and libMesh::System::number().

Referenced by impose_constraint().

470 {
471  for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
472  {
473  const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
474  DofConstraintRow constraint_row;
475  // Leave the constraint row as all zeros so this dof is independent from other dofs
476  const auto constrained_value = node(d);
477  // Simply constrain this dof to retain it's current value
479  constrained_dof_index, constraint_row, constrained_value, true);
480  }
481 }
const MeshBase & get_mesh() const
Definition: system.h:2401
unsigned int number() const
Definition: system.h:2393
void add_constraint_row(const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition: dof_map.h:100
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:176
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
const DofMap & get_dof_map() const
Definition: system.h:2417

◆ get_neighbors_for_boundary_constraint()

std::set< std::set< const Node * > > libMesh::VariationalSmootherConstraint::get_neighbors_for_boundary_constraint ( const MeshBase mesh,
const Node node,
const std::unordered_set< dof_id_type > &  boundary_node_ids,
const BoundaryInfo boundary_info,
const std::unordered_map< dof_id_type, std::vector< const Elem *>> &  nodes_to_elem_map 
)
staticprivate

Get the relevant nodal neighbors for an external boundary constraint.

Parameters
meshThe mesh being smoothed.
nodeThe node (on the external boundary) being constrained.
boundary_node_idsThe set of mesh's external boundary node ids.
boundary_infoThe mesh's BoundaryInfo.
nodes_to_elem_mapA mapping from node id to containing element ids.
Returns
A set of node pointer sets containing nodal neighbors to 'node' on the external boundary. The subsets are grouped by element faces that form the external boundary. Note that 'node' itself does not appear in this set.

Definition at line 717 of file variational_smoother_constraint.C.

References libMesh::MeshTools::find_nodal_or_face_neighbors(), libMesh::DofObject::id(), mesh, and nodes_share_boundary_id().

Referenced by constrain().

723 {
724 
725  // Find all the nodal neighbors... that is the nodes directly connected
726  // to this node through one edge, or if none exists, use other nodes on the
727  // containing face
728  std::vector<const Node *> neighbors;
730  mesh, node, nodes_to_elem_map, neighbors);
731 
732  // Each constituent set corresponds to neighbors sharing a face on the
733  // boundary
734  std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
735 
736  for (const auto * neigh : neighbors)
737  {
738  const bool is_neighbor_boundary_node =
739  boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
740  if (!is_neighbor_boundary_node)
741  continue;
742 
743  // Determine whether nodes share a common boundary id
744  // First, find the common element that both node and neigh belong to
745  const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
746  const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
747  const Elem * common_elem = nullptr;
748  for (const auto * neigh_elem : elems_containing_neigh)
749  {
750  const bool is_neigh_common =
751  std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
752  elems_containing_node.end();
753  if (!is_neigh_common)
754  continue;
755  common_elem = neigh_elem;
756  // Keep this in the neigh_elem loop because there can be multiple common
757  // elements Now, determine whether node and neigh share a common boundary
758  // id
759  const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
760  node, *neigh, *common_elem, boundary_info);
761  if (nodes_have_common_bid)
762  {
763  // Find the side coinciding with the shared boundary
764  for (const auto side : common_elem->side_index_range())
765  {
766  // We only care about external boundaries here, make sure side doesn't
767  // have a neighbor
768  if (common_elem->neighbor_ptr(side))
769  continue;
770 
771  bool node_found_on_side = false;
772  bool neigh_found_on_side = false;
773  const auto nodes_on_side = common_elem->nodes_on_side(side);
774  for (const auto local_node_id : nodes_on_side)
775  {
776  if (common_elem->node_id(local_node_id) == node.id())
777  node_found_on_side = true;
778  else if (common_elem->node_id(local_node_id) == neigh->id())
779  neigh_found_on_side = true;
780  }
781  if (!(node_found_on_side && neigh_found_on_side))
782  continue;
783 
784  std::set<const Node *> node_ptrs_on_side;
785  for (const auto local_node_id : nodes_on_side)
786  node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
787  node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
788  side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
789  }
790  continue;
791  }
792  }
793  }
794 
795  return side_grouped_boundary_neighbors;
796 }
MeshBase & mesh
static bool nodes_share_boundary_id(const Node &boundary_node, const Node &neighbor_node, const Elem &containing_elem, const BoundaryInfo &boundary_info)
Determines whether two neighboring nodes share a common boundary id.
void find_nodal_or_face_neighbors(const MeshBase &mesh, const Node &node, const std::unordered_map< dof_id_type, std::vector< const Elem *>> &nodes_to_elem_map, std::vector< const Node *> &neighbors)
Given a mesh and a node in the mesh, the vector will be filled with every node directly attached to t...
Definition: mesh_tools.C:1092

◆ get_neighbors_for_subdomain_constraint()

std::set< std::set< const Node * > > libMesh::VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint ( const MeshBase mesh,
const Node node,
const subdomain_id_type  sub_id,
const std::unordered_map< dof_id_type, std::vector< const Elem *>> &  nodes_to_elem_map 
)
staticprivate

Get the relevant nodal neighbors for a subdomain constraint.

Parameters
meshThe mesh being smoothed.
nodeThe node (on the subdomain boundary) being constrained.
sub_idThe subdomain id of the block on one side of the subdomain boundary.
nodes_to_elem_mapA mapping from node id to containing element ids.
Returns
A set of node pointer sets containing nodal neighbors to 'node' on the sub_id1-sub_id2 boundary. The subsets are grouped by element faces that form the subdomain boundary. Note that 'node' itself does not appear in this set.

Definition at line 632 of file variational_smoother_constraint.C.

References libMesh::MeshTools::find_nodal_or_face_neighbors(), libMesh::DofObject::id(), mesh, libMesh::Elem::neighbor_ptr(), and libMesh::Elem::subdomain_id().

Referenced by constrain().

637 {
638 
639  // Find all the nodal neighbors... that is the nodes directly connected
640  // to this node through one edge, or if none exists, use other nodes on the
641  // containing face
642  std::vector<const Node *> neighbors;
644  mesh, node, nodes_to_elem_map, neighbors);
645 
646  // Each constituent set corresponds to neighbors sharing a face on the
647  // subdomain boundary
648  std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
649 
650  for (const auto * neigh : neighbors)
651  {
652  // Determine whether the neighbor is on the subdomain boundary
653  // First, find the common elements that both node and neigh belong to
654  const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
655  const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
656  const Elem * common_elem = nullptr;
657  for (const auto * neigh_elem : elems_containing_neigh)
658  {
659  if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
660  elems_containing_node.end())
661  // We should be able to find a common element on the sub_id boundary
662  && (neigh_elem->subdomain_id() == sub_id))
663  common_elem = neigh_elem;
664  else
665  continue;
666 
667  // Now, determine whether node and neigh are on a side coincident
668  // with the subdomain boundary
669  for (const auto common_side : common_elem->side_index_range())
670  {
671  bool node_found_on_side = false;
672  bool neigh_found_on_side = false;
673  for (const auto local_node_id : common_elem->nodes_on_side(common_side))
674  {
675  if (common_elem->node_id(local_node_id) == node.id())
676  node_found_on_side = true;
677  else if (common_elem->node_id(local_node_id) == neigh->id())
678  neigh_found_on_side = true;
679  }
680 
681  if (!(node_found_on_side && neigh_found_on_side &&
682  common_elem->neighbor_ptr(common_side)))
683  continue;
684 
685  const auto matched_side = common_side;
686  // There could be multiple matched sides, so keep this next part
687  // inside the common_side loop
688  //
689  // Does matched_side, containing both node and neigh, lie on the
690  // sub_id subdomain boundary?
691  const auto matched_neighbor_sub_id =
692  common_elem->neighbor_ptr(matched_side)->subdomain_id();
693  const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
694 
695  if (is_matched_side_on_subdomain_boundary)
696  {
697  // Store all nodes that live on this side
698  const auto nodes_on_side = common_elem->nodes_on_side(common_side);
699  std::set<const Node *> node_ptrs_on_side;
700  for (const auto local_node_id : nodes_on_side)
701  node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
702  node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
703  side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
704 
705  continue;
706  }
707 
708  } // for common_side
709 
710  } // for neigh_elem
711  }
712 
713  return side_grouped_boundary_neighbors;
714 }
MeshBase & mesh
void find_nodal_or_face_neighbors(const MeshBase &mesh, const Node &node, const std::unordered_map< dof_id_type, std::vector< const Elem *>> &nodes_to_elem_map, std::vector< const Node *> &neighbors)
Given a mesh and a node in the mesh, the vector will be filled with every node directly attached to t...
Definition: mesh_tools.C:1092

◆ impose_constraint()

void libMesh::VariationalSmootherConstraint::impose_constraint ( const Node node,
const ConstraintVariant constraint 
)
private

Applies a given constraint to a node (e.g., fixing it, restricting it to a line or plane).

Parameters
nodeThe node to constrain.
constraintThe geometric constraint variant to apply.
Exceptions
libMesh::logicErrorIf the constraint cannot be imposed.

Definition at line 888 of file variational_smoother_constraint.C.

References constrain_node_to_line(), constrain_node_to_plane(), and fix_node().

Referenced by constrain().

890 {
891 
892  libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
893  "Cannot impose constraint using InvalidConstraint.");
894 
895  if (std::holds_alternative<PointConstraint>(constraint))
896  fix_node(node);
897  else if (std::holds_alternative<LineConstraint>(constraint))
898  constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
899  else if (std::holds_alternative<PlaneConstraint>(constraint))
900  constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
901 
902  else
903  libmesh_assert_msg(false, "Unknown constraint type.");
904 }
void constrain_node_to_plane(const Node &node, const Point &ref_normal_vec)
Constrain a node to remain in the given plane during mesh smoothing.
void constrain_node_to_line(const Node &node, const Point &line_vec)
Constrain a node to remain on the given line during mesh smoothing.
void fix_node(const Node &node)
Constrain (i.e., fix) a node to not move during mesh smoothing.

◆ nodes_share_boundary_id()

bool libMesh::VariationalSmootherConstraint::nodes_share_boundary_id ( const Node boundary_node,
const Node neighbor_node,
const Elem containing_elem,
const BoundaryInfo boundary_info 
)
staticprivate

Determines whether two neighboring nodes share a common boundary id.

Parameters
boundary_nodeThe first of the two prospective nodes.
neighbor_nodeThe second of the two prospective nodes.
containing_elemThe element containing node1 and node2.
boundary_infoThe mesh's BoundaryInfo.
Returns
nodes_share_bid Whether node1 and node2 share a common boundary id.

Definition at line 600 of file variational_smoother_constraint.C.

References libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::get_node_index(), libMesh::Elem::is_node_on_side(), side_id, and libMesh::Elem::side_index_range().

Referenced by get_neighbors_for_boundary_constraint().

604 {
605  bool nodes_share_bid = false;
606 
607  // Node ids local to containing_elem
608  const auto node_id = containing_elem.get_node_index(&boundary_node);
609  const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
610 
611  for (const auto side_id : containing_elem.side_index_range())
612  {
613  // We don't care about this side if it doesn't contain our boundary and neighbor nodes
614  if (!(containing_elem.is_node_on_side(node_id, side_id) &&
615  containing_elem.is_node_on_side(neighbor_id, side_id)))
616  continue;
617 
618  // If the current side, containing boundary_node and neighbor_node, lies on a boundary,
619  // we can say that boundary_node and neighbor_node have a common boundary id.
620  std::vector<boundary_id_type> boundary_ids;
621  boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
622  if (boundary_ids.size())
623  {
624  nodes_share_bid = true;
625  break;
626  }
627  }
628  return nodes_share_bid;
629 }
const boundary_id_type side_id

Member Data Documentation

◆ _preserve_subdomain_boundaries

const bool libMesh::VariationalSmootherConstraint::_preserve_subdomain_boundaries
private

Whether nodes on subdomain boundaries are subject to change via smoothing.

Definition at line 410 of file variational_smoother_constraint.h.

Referenced by constrain().

◆ _sys

System& libMesh::VariationalSmootherConstraint::_sys
private

◆ _verbosity

const unsigned int libMesh::VariationalSmootherConstraint::_verbosity
private

Verbosity setting.

The verbosity levels and the corresponding information output are as follows:

verbosity = 0: No information.

20 < verbosity: Prints:

  • Constraint information on all boundary and subdomain nodes

Definition at line 403 of file variational_smoother_constraint.h.

Referenced by constrain().


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