Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #ifndef LIBMESH_VARIATIONAL_SMOOTHER_CONSTRAINT_H
19 : #define LIBMESH_VARIATIONAL_SMOOTHER_CONSTRAINT_H
20 :
21 : // Local Includes
22 : #include "libmesh/system.h"
23 : #include "libmesh/dof_map.h"
24 :
25 : // C++ includes
26 : #include <variant>
27 :
28 : namespace libMesh
29 : {
30 :
31 : // Forward declarations
32 : class PointConstraint;
33 : class LineConstraint;
34 : class PlaneConstraint;
35 : class InvalidConstraint;
36 :
37 : /**
38 : * Type used to store a constraint that may be a PlaneConstraint,
39 : * LineConstraint, or PointConstraint. std::variant is an alternative to using
40 : * the classic polymorphic approach where these constraints inherit from an
41 : * common base class.
42 : */
43 : using ConstraintVariant = std::variant<PointConstraint, LineConstraint,
44 : PlaneConstraint, InvalidConstraint>;
45 :
46 : /**
47 : * Represents a fixed point constraint.
48 : */
49 : class PointConstraint
50 : {
51 :
52 : public:
53 714 : PointConstraint() = default;
54 :
55 : /**
56 : * Constructor
57 : * @param point The point defining the constraint.
58 : * @param tol The tolerance to use for numerical comparisons.
59 : */
60 14295970 : PointConstraint(const Point &point, const Real &tol = TOLERANCE * TOLERANCE);
61 :
62 : /**
63 : * Comparison operator for ordering PointConstraint objects.
64 : * A PointConstraint is considered less than another if its location
65 : * is lexicographically less than the other's location.
66 : * @param other The PointConstraint to compare with.
67 : * @param tol The tolerance to use for numerical comparisons.
68 : * @return True if this PointConstraint is less than the other.
69 : */
70 : bool operator<(const PointConstraint &other) const;
71 :
72 : /**
73 : * Equality operator.
74 : * @param other The PointConstraint to compare with.
75 : * @return True if both PointConstraints have the same location.
76 : */
77 : bool operator==(const PointConstraint &other) const;
78 :
79 : /**
80 : * Query whether a point lies on another point.
81 : * @param p The point in question
82 : * @return bool indicating whether p lies on this point.
83 : */
84 0 : bool contains_point(const PointConstraint &p) const { return *this == p; }
85 :
86 : /**
87 : * Computes the intersection of this point with another constraint.
88 : * Handles intersection with PointConstraint, LineConstraint, or
89 : * PlaneConstraint.
90 : * @param other The constraint to intersect with.
91 : * @return The most specific ConstraintVariant that satisfies both
92 : * constraints. constraints. If no intersection exists, return an
93 : * InvalidConstraint.
94 : */
95 : ConstraintVariant intersect(const ConstraintVariant &other) const;
96 :
97 : /**
98 : * Const getter for the _point attribute
99 : */
100 418575 : const Point &point() const { return _point; }
101 :
102 : /**
103 : * Const getter for the _tol attribute
104 : */
105 : const Real &tol() const { return _tol; }
106 :
107 : private:
108 : // Life is easier if we don't make this const
109 : /**
110 : * Location of constraint
111 : */
112 : Point _point;
113 :
114 : /**
115 : * Tolerance to use for numerical comparisons
116 : */
117 : Real _tol;
118 : };
119 :
120 : /**
121 : * Represents a line constraint defined by a base point and direction vector.
122 : */
123 : class LineConstraint
124 : {
125 : public:
126 : LineConstraint() = default;
127 :
128 : /**
129 : * Constructor
130 : * @param point A point on the constraining line.
131 : * @param direction the direction of the constraining line.
132 : * @param tol The tolerance to use for numerical comparisons.
133 : */
134 : LineConstraint(const Point &point, const Point &direction,
135 60955 : const Real &tol = TOLERANCE * TOLERANCE);
136 :
137 : /**
138 : * Comparison operator for ordering LineConstraint objects.
139 : * The comparison is primarily based on the direction vector. If the direction
140 : * vectors are equal (within tolerance), the tie is broken using the dot
141 : * product of the direction with the base point.
142 : * @param other The LineConstraint to compare with.
143 : * @return True if this LineConstraint is less than the other.
144 : */
145 : bool operator<(const LineConstraint &other) const;
146 :
147 : /**
148 : * Equality operator.
149 : * @param other The LineConstraint to compare with.
150 : * @return True if both LineConstraints represent the same line.
151 : */
152 : bool operator==(const LineConstraint &other) const;
153 :
154 : /**
155 : * Query whether a point lies on the line.
156 : * @param p The point in question
157 : * @return bool indicating whether p lies on the line.
158 : */
159 : bool contains_point(const PointConstraint &p) const;
160 :
161 : /**
162 : * Query whether a line is parallel to this line
163 : * @param l The line in question
164 : * @return bool indicating whether l is parallel to this line.
165 : */
166 : bool is_parallel(const LineConstraint &l) const;
167 :
168 : /**
169 : * Query whether a plane is parallel to this line
170 : * @param p The plane in question
171 : * @return bool indicating whether p is parallel to this line.
172 : */
173 : bool is_parallel(const PlaneConstraint &p) const;
174 :
175 : /**
176 : * Computes the intersection of this line with another constraint.
177 : * Handles intersection with LineConstraint, PlaneConstraint, or
178 : * PointConstraint.
179 : * @param other The constraint to intersect with.
180 : * @return The most specific ConstraintVariant that satisfies both
181 : * constraints. constraints. If no intersection exists, return an
182 : * InvalidConstraint.
183 : */
184 : ConstraintVariant intersect(const ConstraintVariant &other) const;
185 :
186 : /**
187 : * Const getter for the _point attribute
188 : */
189 50926 : const Point &point() const { return _point; }
190 :
191 : /**
192 : * Const getter for the _direction attribute
193 : */
194 34933 : const Point &direction() const { return _direction; }
195 :
196 : /**
197 : * Const getter for the _tol attribute
198 : */
199 : const Real &tol() const { return _tol; }
200 :
201 : private:
202 : // Life is easier if we don't make these const
203 : /**
204 : * A point on the constraining line
205 : */
206 : Point _point;
207 :
208 : /**
209 : * Direction of the constraining line
210 : */
211 : Point _direction;
212 :
213 : /**
214 : * Tolerance to use for numerical comparisons
215 : */
216 : Real _tol;
217 : };
218 :
219 : /**
220 : * Represents a plane constraint defined by a point and normal vector.
221 : */
222 : class PlaneConstraint
223 : {
224 :
225 : public:
226 : PlaneConstraint() = default;
227 :
228 : /**
229 : * Constructor
230 : * @param point A point on the constraining plane.
231 : * @param normal the direction normal to the constraining plane.
232 : * @param tol The tolerance to use for numerical comparisons.
233 : */
234 : PlaneConstraint(const Point &point, const Point &normal,
235 8125524 : const Real &tol = TOLERANCE * TOLERANCE);
236 :
237 : /**
238 : * Comparison operator for ordering PlaneConstraint objects.
239 : * The comparison is primarily based on the normal vector. If the normal
240 : * vectors are equal (within tolerance), the tie is broken using the dot
241 : * product of the normal with the point on the plane.
242 : * @param other The PlaneConstraint to compare with.
243 : * @return True if this PlaneConstraint is less than the other.
244 : */
245 : bool operator<(const PlaneConstraint &other) const;
246 :
247 : /**
248 : * Equality operator.
249 : * @param other The PlaneConstraint to compare with.
250 : * @return True if both PlaneConstraints represent the same plane.
251 : */
252 : bool operator==(const PlaneConstraint &other) const;
253 :
254 : /**
255 : * Query whether a point lies on the plane.
256 : * @param p The point in question
257 : * @return bool indicating whether p lies on the plane.
258 : */
259 : bool contains_point(const PointConstraint &p) const;
260 :
261 : /**
262 : * Query whether a line lies on the plane.
263 : * @param l The line in question
264 : * @return bool indicating whether l lies on the plane.
265 : */
266 : bool contains_line(const LineConstraint &l) const;
267 :
268 : /**
269 : * Query whether a plane is parallel to this plane
270 : * @param p The plane in question
271 : * @return bool indicating whether p is parallel to this plane.
272 : */
273 : bool is_parallel(const PlaneConstraint &p) const;
274 :
275 : /**
276 : * Query whether a line is parallel to this plane
277 : * @param l The line in question
278 : * @return bool indicating whether l is parallel to this plane.
279 : */
280 : bool is_parallel(const LineConstraint &l) const;
281 :
282 : /**
283 : * Computes the intersection of this plane with another constraint.
284 : * Handles intersection with PlaneConstraint, LineConstraint, or
285 : * PointConstraint.
286 : * @param other The constraint to intersect with.
287 : * @return The most specific ConstraintVariant that satisfies both
288 : * constraints. constraints. If no intersection exists, return an
289 : * InvalidConstraint.
290 : */
291 : ConstraintVariant intersect(const ConstraintVariant &other) const;
292 :
293 : /**
294 : * Const getter for the _point attribute
295 : */
296 13873531 : const Point &point() const { return _point; }
297 :
298 : /**
299 : * Const getter for the _normal attribute
300 : */
301 1188789 : const Point &normal() const { return _normal; }
302 :
303 : /**
304 : * Const getter for the _tol attribute
305 : */
306 : const Real &tol() const { return _tol; }
307 :
308 : private:
309 : // Life is easier if we don't make these const
310 : /**
311 : * A point on the constraining plane
312 : */
313 : Point _point;
314 :
315 : /**
316 : * The direction normal to the constraining plane
317 : */
318 : Point _normal;
319 :
320 : /**
321 : * Tolerance to use for numerical comparisons
322 : */
323 : Real _tol;
324 : };
325 :
326 : /**
327 : * Represents an invalid constraint (i.e., when the two constraints don't
328 : * intersect)
329 : */
330 268 : class InvalidConstraint
331 : {
332 :
333 : public:
334 2 : InvalidConstraint()
335 71 : : _err_msg("We should never get here! The InvalidConstraint object should be "
336 : "detected and replaced with a valid ConstraintVariant prior to calling "
337 4 : "any class methods.")
338 : {
339 2 : }
340 :
341 : /**
342 : * Dummy intersect method that should never be called.
343 : */
344 0 : ConstraintVariant intersect(const ConstraintVariant &) const {
345 0 : libmesh_assert_msg(false, _err_msg);
346 0 : return *this;
347 : }
348 :
349 : /**
350 : * Dummy contains_point method that should never be called.
351 : */
352 0 : bool contains_point(const PointConstraint &) const {
353 0 : libmesh_assert_msg(false, _err_msg);
354 0 : return false;
355 : }
356 :
357 : private:
358 : std::string _err_msg;
359 : };
360 :
361 : /**
362 : * Dispatch intersection between two constraint variants.
363 : * Resolves to the appropriate method based on the type of the first operand.
364 : * @param a First constraint.
365 : * @param b Constraint to combine with a.
366 : * @return Combination (intersection) of constraint a and b.
367 : */
368 29286 : inline ConstraintVariant intersect_constraints(const ConstraintVariant &a,
369 : const ConstraintVariant &b) {
370 : // std::visit applies the visitor v (a Callable that can be called with any
371 : // combination of types from Variants) to the active value inside a
372 : // std::Variant. This circumvents the issue that the literal ConstraintVariant
373 : // type does not have a method called 'intersect' (but the types defining
374 : // ConstraintVariant do)
375 : return std::visit(
376 2020212 : [](const auto &lhs, const auto &rhs) -> ConstraintVariant {
377 2049494 : return lhs.intersect(rhs);
378 : },
379 1039390 : a, b);
380 : }
381 :
382 : /**
383 : * Constraint class for the VariationalMeshSmoother.
384 : *
385 : * Currently, all mesh boundary nodes are constrained to not move during
386 : * smoothing. If requested (preserve_subdomain_boundaries = true), nodes on
387 : * subdomain boundaries are also constrained to not move.
388 : */
389 48 : class VariationalSmootherConstraint : public System::Constraint
390 : {
391 : private:
392 :
393 : System & _sys;
394 :
395 : /**
396 : * Whether nodes on subdomain boundaries are subject to change via smoothing
397 : */
398 : const bool _preserve_subdomain_boundaries;
399 :
400 : /**
401 : * Constrain (i.e., fix) a node to not move during mesh smoothing.
402 : * @param node Node to fix.
403 : */
404 : void fix_node(const Node & node);
405 :
406 : /**
407 : * Constrain a node to remain in the given plane during mesh smoothing.
408 : * @param node Node to constrain
409 : * @param ref_normal_vec Reference normal vector to the constraining plane.
410 : * This, along with the coordinates of node, are used to define the
411 : * constraining plane.
412 : */
413 : void constrain_node_to_plane(const Node & node, const Point & ref_normal_vec);
414 :
415 : /**
416 : * Constrain a node to remain on the given line during mesh smoothing.
417 : * @param node Node to constrain
418 : * @param line_vec vector parallel to the constraining line.
419 : * This, along with the coordinates of node, are used to define the
420 : * constraining line.
421 : */
422 : void constrain_node_to_line(const Node & node, const Point & line_vec);
423 :
424 : /**
425 : * Given a mesh and a node in the mesh, the vector will be filled with
426 : * every node directly attached to the given one. IF NO neighbors are found,
427 : * all nodes on the containing side are treated as neighbors. This is useful
428 : * when the node does not lie on an edge, such as the central face node in
429 : * HEX27 elements.
430 : */
431 : static void find_nodal_or_face_neighbors(
432 : const MeshBase & mesh,
433 : const Node & node,
434 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
435 : std::vector<const Node *> & neighbors);
436 :
437 : /**
438 : * Determines whether two neighboring nodes share a common boundary id.
439 : * @param boundary_node The first of the two prospective nodes.
440 : * @param neighbor_node The second of the two prospective nodes.
441 : * @param containing_elem The element containing node1 and node2.
442 : * @param boundary_info The mesh's BoundaryInfo.
443 : * @return nodes_share_bid Whether node1 and node2 share a common boundary id.
444 : */
445 : static bool nodes_share_boundary_id(
446 : const Node & boundary_node,
447 : const Node & neighbor_node,
448 : const Elem & containing_elem,
449 : const BoundaryInfo & boundary_info);
450 :
451 : /**
452 : * Get the relevant nodal neighbors for a subdomain constraint.
453 : * @param mesh The mesh being smoothed.
454 : * @param node The node (on the subdomain boundary) being constrained.
455 : * @param sub_id The subdomain id of the block on one side of the subdomain
456 : * boundary.
457 : * @param nodes_to_elem_map A mapping from node id to containing element ids.
458 : * @return A set of node pointer sets containing nodal neighbors to 'node' on
459 : * the sub_id1-sub_id2 boundary. The subsets are grouped by element faces
460 : * that form the subdomain boundary. Note that 'node' itself does not appear
461 : * in this set.
462 : */
463 : static std::set<std::set<const Node *>>
464 : get_neighbors_for_subdomain_constraint(
465 : const MeshBase &mesh, const Node &node, const subdomain_id_type sub_id,
466 : const std::unordered_map<dof_id_type, std::vector<const Elem *>>
467 : &nodes_to_elem_map);
468 :
469 : /**
470 : * Get the relevant nodal neighbors for an external boundary constraint.
471 : * @param mesh The mesh being smoothed.
472 : * @param node The node (on the external boundary) being constrained.
473 : * @param boundary_node_ids The set of mesh's external boundary node ids.
474 : * @param boundary_info The mesh's BoundaryInfo.
475 : * @param nodes_to_elem_map A mapping from node id to containing element ids.
476 : * @return A set of node pointer sets containing nodal neighbors to 'node' on
477 : * the external boundary. The subsets are grouped by element faces that form
478 : * the external boundary. Note that 'node' itself does not appear in this
479 : * set.
480 : */
481 : static std::set<std::set<const Node *>> get_neighbors_for_boundary_constraint(
482 : const MeshBase &mesh, const Node &node,
483 : const std::unordered_set<dof_id_type> &boundary_node_ids,
484 : const BoundaryInfo &boundary_info,
485 : const std::unordered_map<dof_id_type, std::vector<const Elem *>>
486 : &nodes_to_elem_map);
487 :
488 : /**
489 : * Determines the appropriate constraint (PointConstraint, LineConstraint, or
490 : * PlaneConstraint) for a node based on its neighbors.
491 : * @param node The node to constrain.
492 : * @param dim The mesh dimension.
493 : * @return The best-fit constraint for the given geometry.
494 : */
495 : static ConstraintVariant determine_constraint(
496 : const Node &node, const unsigned int dim,
497 : const std::set<std::set<const Node *>> &side_grouped_boundary_neighbors);
498 :
499 : /**
500 : * Applies a given constraint to a node (e.g., fixing it, restricting it to a
501 : * line or plane).
502 : * @param node The node to constrain.
503 : * @param constraint The geometric constraint variant to apply.
504 : * @throw libMesh::logicError If the constraint cannot be imposed.
505 : */
506 : void impose_constraint(const Node &node, const ConstraintVariant &constraint);
507 :
508 : public:
509 : /**
510 : * Constructor
511 : * @param sys System to constrain.
512 : * @param preserve_subdomain_boundaries Whether to constrain nodes on
513 : * subdomain boundaries to not move.
514 : */
515 : VariationalSmootherConstraint(System & sys, const bool & preserve_subdomain_boundaries);
516 :
517 : virtual ~VariationalSmootherConstraint() override;
518 :
519 : virtual void constrain() override;
520 : };
521 :
522 :
523 : } // namespace libMesh
524 :
525 : #endif // LIBMESH_VARIATIONAL_SMOOTHER_CONSTRAINT_H
|
