Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local Includes
19 : #include "libmesh/variational_smoother_constraint.h"
20 : #include "libmesh/mesh_tools.h"
21 : #include "libmesh/boundary_info.h"
22 :
23 : namespace libMesh
24 : {
25 :
26 : // Helper function to orient a unit vector in the positive x/y/z direction
27 795061 : auto get_positive_vector = [](const Point & vec) -> Point {
28 : // unit() rounds to within a few ulps of length 1, so TOLERANCE only
29 : // rejects vectors that were never normalized
30 66021 : libmesh_assert(absolute_fuzzy_equals(vec.norm(), 1., TOLERANCE));
31 :
32 : // Choose sign such that direction vector points in the positive x/y/z
33 : // direction This helps to eliminate duplicate lines/planes
34 66021 : Point canonical{0, 0, 0};
35 : // Choose the canonical dimension to ensure the dot product below is nonzero
36 1542106 : for (const auto dim_id : make_range(3))
37 1670170 : if (!absolute_fuzzy_equals(vec(dim_id), 0.))
38 : {
39 795061 : canonical(dim_id) = 1.;
40 729040 : break;
41 : }
42 :
43 66021 : const auto dot_prod = vec * canonical;
44 66021 : libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));
45 :
46 795061 : return (dot_prod > 0) ? vec : Point(-vec);
47 : };
48 :
49 1469398 : PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)
50 : {
51 1469398 : }
52 :
53 0 : bool PointConstraint::operator<(const PointConstraint & other) const
54 : {
55 0 : if (*this == other)
56 0 : return false;
57 :
58 0 : return _point < other.point();
59 : }
60 :
61 192 : bool PointConstraint::operator==(const PointConstraint & other) const
62 : {
63 192 : return _point.absolute_fuzzy_equals(other.point(), _tol);
64 : }
65 :
66 13868 : ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const
67 : {
68 : // using visit to resolve the variant to its actual type
69 : return std::visit(
70 37040 : [&](auto && o) -> ConstraintVariant {
71 13868 : if (!o.contains_point(*this))
72 : // Point is not on the constraint
73 0 : return InvalidConstraint();
74 :
75 1141 : return *this;
76 : },
77 13868 : other);
78 : }
79 :
80 13270 : LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
81 13270 : : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
82 : {
83 1098 : libmesh_assert_msg(_direction.norm() >= _tol,
84 : "Can't define a line with zero magnitude direction vector.");
85 13270 : }
86 :
87 0 : bool LineConstraint::operator<(const LineConstraint & other) const
88 : {
89 0 : if (*this == other)
90 0 : return false;
91 :
92 0 : if (!(_direction.absolute_fuzzy_equals(other.direction(), _tol)))
93 0 : return _direction < other.direction();
94 0 : return (_direction * _point) < (other.direction() * other.point());
95 : }
96 :
97 48 : bool LineConstraint::operator==(const LineConstraint & other) const
98 : {
99 52 : if (!(_direction.absolute_fuzzy_equals(other.direction(), _tol)))
100 4 : return false;
101 0 : return this->contains_point(other.point());
102 : }
103 :
104 599 : bool LineConstraint::contains_point(const PointConstraint & p) const
105 : {
106 : // If the point lies on the line, then the vector p - point is parallel to the
107 : // line In that case, the cross product of p - point with the line's direction
108 : // will be zero.
109 599 : return _direction.cross(p.point() - _point).norm() < _tol;
110 : }
111 :
112 48 : bool LineConstraint::is_parallel(const LineConstraint & l) const
113 : {
114 48 : return _direction.absolute_fuzzy_equals(l.direction(), _tol);
115 : }
116 :
117 3453 : bool LineConstraint::is_parallel(const PlaneConstraint & p) const
118 : {
119 3453 : return std::abs(_direction * p.normal()) < _tol;
120 : }
121 :
122 4139 : ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const
123 : {
124 : // using visit to resolve the variant to its actual type
125 : return std::visit(
126 11053 : [&](auto && o) -> ConstraintVariant {
127 : // Use std::decay_t to strip references/const/volatile from o's type,
128 : // so we can match the actual stored variant type in std::is_same_v.
129 : using T = std::decay_t<decltype(o)>;
130 : if constexpr (std::is_same_v<T, LineConstraint>)
131 : {
132 4183 : if (*this == o)
133 0 : return *this;
134 :
135 48 : if (this->is_parallel(o))
136 : // Lines are parallel and do not intersect
137 0 : return InvalidConstraint();
138 :
139 : // Solve for t in the equation p1 + t·d1 = p2 + s·d2
140 : // The shortest vector between skew lines lies along the normal vector
141 : // (d1 × d2). Projecting the vector (p2 - p1) onto this normal gives a
142 : // scalar proportional to the distance. This is equivalent to solving:
143 : // ((p2 - p1) × d2) · (d1 × d2) = t · |d1 × d2|²
144 : // ⇒ t = ((delta × d2) · (d1 × d2)) / |d1 × d2|²
145 :
146 4 : const Point delta = o.point() - _point;
147 44 : const Point cross_d1_d2 = _direction.cross(o.direction());
148 44 : const Real cross_dot = (delta.cross(o.direction())) * cross_d1_d2;
149 4 : const Real denom = cross_d1_d2.norm_sq();
150 :
151 48 : const Real t = cross_dot / denom;
152 4 : const Point intersection = _point + t * _direction;
153 :
154 : // Verify that intersection lies on both lines
155 48 : if (o.direction().cross(intersection - o.point()).norm() > _tol)
156 : // Lines do not intersect at a single point
157 0 : return InvalidConstraint();
158 :
159 48 : return PointConstraint{intersection};
160 : }
161 :
162 : else if constexpr (std::is_same_v<T, PlaneConstraint>)
163 7845 : return o.intersect(*this);
164 :
165 : else if constexpr (std::is_same_v<T, PointConstraint>)
166 : {
167 0 : if (!this->contains_point(o))
168 : // Point is not on the line
169 0 : return InvalidConstraint();
170 :
171 0 : return o;
172 : }
173 :
174 : else
175 0 : libmesh_error_msg("Unsupported constraint type in Line::intersect.");
176 : },
177 4139 : other);
178 : }
179 :
180 781791 : PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
181 781791 : : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
182 : {
183 64923 : libmesh_assert_msg(_normal.norm() >= _tol,
184 : "Can't define a plane with zero magnitude direction vector.");
185 781791 : }
186 :
187 1817863 : bool PlaneConstraint::operator<(const PlaneConstraint & other) const
188 : {
189 1817863 : if (*this == other)
190 121302 : return false;
191 :
192 386796 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
193 357424 : return _normal < other.normal();
194 0 : return (_normal * _point) < (other.normal() * other.point());
195 : }
196 :
197 1828709 : bool PlaneConstraint::operator==(const PlaneConstraint & other) const
198 : {
199 1980279 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
200 30379 : return false;
201 1460439 : return this->contains_point(other.point());
202 : }
203 :
204 10846 : bool PlaneConstraint::is_parallel(const PlaneConstraint & p) const
205 : {
206 : // |n1 x n2| is the sine of the angle between the unit normals, and its
207 : // square is the determinant that intersect() divides by, so any pair of
208 : // planes that fails this test can be intersected.
209 10846 : return _normal.cross(p.normal()).norm() <= _tol;
210 : }
211 :
212 3382 : bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }
213 :
214 1477893 : bool PlaneConstraint::contains_point(const PointConstraint & p) const
215 : {
216 : // distance between the point and the plane
217 122734 : const Real dist = (p.point() - _point) * _normal;
218 1477893 : return std::abs(dist) < _tol;
219 : }
220 :
221 4306 : bool PlaneConstraint::contains_line(const LineConstraint & l) const
222 : {
223 4306 : const bool base_on_plane = this->contains_point(PointConstraint(l.point()));
224 4306 : const bool dir_orthogonal = std::abs(_normal * l.direction()) < _tol;
225 4306 : return base_on_plane && dir_orthogonal;
226 : }
227 :
228 15081 : ConstraintVariant PlaneConstraint::intersect(const ConstraintVariant & other) const
229 : {
230 : // using visit to resolve the variant to its actual type
231 : return std::visit(
232 40263 : [&](auto && o) -> ConstraintVariant {
233 : // Use std::decay_t to strip references/const/volatile from o's type,
234 : // so we can match the actual stored variant type in std::is_same_v.
235 : using T = std::decay_t<decltype(o)>;
236 : if constexpr (std::is_same_v<T, PlaneConstraint>)
237 : {
238 : // If planes are identical, return one of them
239 16498 : if (*this == o)
240 0 : return *this;
241 :
242 10775 : if (this->is_parallel(o))
243 : // Planes are parallel and do not intersect
244 0 : return InvalidConstraint();
245 :
246 : // Solve for a point on the intersection line of two planes.
247 : // Given planes:
248 : // Plane 1: n1 · (x - p1) = 0
249 : // Plane 2: n2 · (x - p2) = 0
250 : // The line of intersection has direction dir = n1 × n2.
251 : // To find a point on this line, we assume:
252 : // x = p1 + s·n1 = p2 + t·n2
253 : // ⇒ p1 - p2 = t·n2 - s·n1
254 : // Taking dot products with n1 and n2 leads to:
255 : // [-n1·n1 n1·n2] [s] = [n1 · (p1 - p2)]
256 : // [-n1·n2 n2·n2] [t] [n2 · (p1 - p2)]
257 :
258 9881 : const Point dir = this->_normal.cross(o.normal()); // direction of line of intersection
259 894 : libmesh_assert(dir.norm() > _tol);
260 894 : const Point w = _point - o.point();
261 :
262 : // Dot product terms used in 2x2 system
263 894 : const Real n1_dot_n2 = _normal * o.normal();
264 894 : const Real n2_dot_n2 = o.normal() * o.normal();
265 894 : const Real n1_dot_w = _normal * w;
266 894 : const Real n2_dot_w = o.normal() * w;
267 :
268 : // By Lagrange's identity the determinant
269 : // -(n1.n1 * n2.n2 - (n1.n2)^2) equals -|n1 x n2|^2, which is
270 : // nonzero because the planes are not parallel. Computing it from
271 : // dir avoids the cancellation in 1 - cos^2 for nearly parallel
272 : // planes.
273 894 : const Real denom = -dir.norm_sq();
274 :
275 10775 : const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
276 894 : const Point p0 = _point + s * _normal;
277 :
278 10775 : return LineConstraint{p0, dir.unit()};
279 : }
280 :
281 : else if constexpr (std::is_same_v<T, LineConstraint>)
282 : {
283 4306 : if (this->contains_line(o))
284 77 : return o;
285 :
286 3382 : if (this->is_parallel(o))
287 : // Line is parallel and does not intersect the plane
288 0 : return InvalidConstraint();
289 :
290 : // Solve for t in the parametric equation:
291 : // p(t) = point + t·d
292 : // such that this point also satisfies the plane equation:
293 : // n · (p(t) - p0) = 0
294 : // which leads to:
295 : // t = (n · (p0 - point)) / (n · d)
296 :
297 274 : const Real denom = _normal * o.direction();
298 274 : libmesh_assert(std::abs(denom) > _tol);
299 3382 : const Real t = (_normal * (_point - o.point())) / denom;
300 3382 : return PointConstraint{o.point() + t * o.direction()};
301 : }
302 :
303 : else if constexpr (std::is_same_v<T, PointConstraint>)
304 : {
305 0 : if (!this->contains_point(o))
306 : // Point is not on the plane
307 0 : return InvalidConstraint();
308 :
309 0 : return o;
310 : }
311 :
312 : else
313 0 : libmesh_error_msg("Unsupported constraint type in Plane::intersect.");
314 : },
315 15081 : other);
316 : }
317 :
318 2520 : std::ostream &operator<<(std::ostream &os, const ConstraintVariant & c)
319 : {
320 2520 : if (std::holds_alternative<PointConstraint>(c))
321 264 : os << "point constraint: point=" << std::get<PointConstraint>(c).point();
322 :
323 2376 : else if (std::holds_alternative<LineConstraint>(c))
324 : {
325 48 : const auto & line = std::get<LineConstraint>(c);
326 624 : os << "line constraint: point=" << line.point() << ", direction=" << line.direction();
327 : }
328 :
329 1800 : else if (std::holds_alternative<PlaneConstraint>(c))
330 : {
331 150 : const auto & plane = std::get<PlaneConstraint>(c);
332 1950 : os << "plane constraint: point=" << plane.point() << ", normal=" << plane.normal();
333 : }
334 :
335 0 : else if (std::holds_alternative<InvalidConstraint>(c))
336 0 : os << "invalid constraint";
337 :
338 2520 : return os;
339 : }
340 :
341 2556 : VariationalSmootherConstraint::VariationalSmootherConstraint(
342 2556 : System & sys, const bool & preserve_subdomain_boundaries, const unsigned int verbosity)
343 2556 : : Constraint(), _verbosity(verbosity), _sys(sys), _preserve_subdomain_boundaries(preserve_subdomain_boundaries) {
344 2556 : }
345 :
346 4824 : VariationalSmootherConstraint::~VariationalSmootherConstraint() = default;
347 :
348 2556 : void VariationalSmootherConstraint::constrain()
349 : {
350 2556 : const auto & mesh = _sys.get_mesh();
351 2556 : const auto dim = mesh.mesh_dimension();
352 144 : const auto & proc_id = mesh.processor_id();
353 :
354 : // Only compute the node to elem map once
355 144 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
356 2556 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
357 :
358 72 : const auto & boundary_info = mesh.get_boundary_info();
359 :
360 2628 : const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
361 :
362 : // Identify/constrain subdomain boundary nodes, if requested
363 144 : std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
364 2556 : if (_preserve_subdomain_boundaries)
365 : {
366 260477 : for (const auto * elem : mesh.active_element_ptr_range())
367 : {
368 10896 : const auto & sub_id1 = elem->subdomain_id();
369 490089 : for (const auto side : elem->side_index_range())
370 : {
371 397586 : const auto * neighbor = elem->neighbor_ptr(side);
372 397586 : if (neighbor == nullptr)
373 42510 : continue;
374 :
375 41360 : const auto & sub_id2 = neighbor->subdomain_id();
376 352520 : if (sub_id1 == sub_id2)
377 323386 : continue;
378 :
379 : // elem and neighbor are in different subdomains, and share nodes
380 : // that need to be constrained
381 47282 : for (const auto local_node_id : elem->nodes_on_side(side))
382 : {
383 39564 : const auto & node = mesh.node_ref(elem->node_id(local_node_id));
384 : // Make sure we haven't already processed this node.
385 : // We leave the responsibility of computing the constraint to
386 : // the proc that owns the node.
387 39094 : if (subdomain_boundary_map.count(node.id()) ||
388 1748 : node.processor_id() != proc_id)
389 35796 : continue;
390 :
391 : // Get the relevant nodal neighbors for the subdomain constraint
392 : const auto side_grouped_boundary_neighbors =
393 : get_neighbors_for_subdomain_constraint(
394 2626 : mesh, node, sub_id1, nodes_to_elem_map);
395 :
396 : // Determine which constraint should be imposed
397 : const auto subdomain_constraint =
398 2424 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
399 :
400 : // This subdomain boundary node does not lie on an external boundary,
401 : // go ahead and impose constraint
402 2626 : if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
403 : {
404 912 : if (_verbosity > 20)
405 0 : libMesh::out << "Imposing subdomain constraint on "
406 0 : << std::endl << " " << node << std::endl
407 0 : << " " << subdomain_constraint << std::endl;
408 :
409 912 : this->impose_constraint(node, subdomain_constraint);
410 : }
411 :
412 : // This subdomain boundary node could lie on an external boundary, save it
413 : // for later to combine with the external boundary constraint.
414 : // We also save constraints for non-boundary nodes so we don't try to
415 : // re-constrain the node when accessed from the neighboring elem.
416 : // See subdomain_boundary_map.count call above.
417 4646 : subdomain_boundary_map[node.id()] = subdomain_constraint;
418 :
419 : } // for local_node_id
420 :
421 : } // for side
422 2278 : } // for elem
423 : }
424 :
425 : // Loop through boundary nodes and impose constraints
426 104020 : for (const auto & bid : boundary_node_ids)
427 : {
428 101464 : const auto & node = mesh.node_ref(bid);
429 : // We leave the responsibility of computing the constraint to the proc
430 : // that owns the node.
431 101464 : if (node.processor_id() != proc_id)
432 76216 : continue;
433 :
434 : // Get the relevant nodal neighbors for the boundary constraint
435 : const auto side_grouped_boundary_neighbors = get_neighbors_for_boundary_constraint(
436 27352 : mesh, node, boundary_node_ids, boundary_info, nodes_to_elem_map);
437 :
438 : // Determine which constraint should be imposed
439 : const auto boundary_constraint =
440 27352 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
441 :
442 : // Check for the case where this boundary node is also part of a subdomain id boundary
443 25248 : if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
444 : {
445 1512 : const auto & subdomain_constraint = it->second;
446 : // Combine current boundary constraint with previously determined
447 : // subdomain_constraint
448 252 : auto combined_constraint = intersect_constraints(subdomain_constraint, boundary_constraint);
449 :
450 : // This will catch cases where constraints have no intersection
451 : // Fall back to fixed node constraint
452 1512 : if (std::holds_alternative<InvalidConstraint>(combined_constraint))
453 0 : combined_constraint = PointConstraint(node);
454 :
455 1512 : if (_verbosity > 20)
456 0 : libMesh::out << "Imposing boundary/subdomain constraint on "
457 0 : << std::endl << " " << node << std::endl
458 0 : << " " << combined_constraint << std::endl;
459 :
460 1512 : this->impose_constraint(node, combined_constraint);
461 : }
462 :
463 : else
464 : {
465 23736 : if (_verbosity > 20)
466 210 : libMesh::out << "Imposing boundary constraint on "
467 210 : << std::endl << node << std::endl
468 210 : << " " << boundary_constraint << std::endl;
469 :
470 23736 : this->impose_constraint(node, boundary_constraint);
471 : }
472 :
473 : } // end bid
474 2556 : }
475 :
476 4248 : void VariationalSmootherConstraint::fix_node(const Node & node)
477 : {
478 15708 : for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
479 : {
480 11460 : const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
481 1910 : DofConstraintRow constraint_row;
482 : // Leave the constraint row as all zeros so this dof is independent from other dofs
483 11460 : const auto constrained_value = node(d);
484 : // Simply constrain this dof to retain it's current value
485 11460 : _sys.get_dof_map().add_constraint_row(
486 : constrained_dof_index, constraint_row, constrained_value, true);
487 : }
488 4248 : }
489 :
490 12648 : void VariationalSmootherConstraint::constrain_node_to_plane(const Node & node,
491 : const Point & ref_normal_vec)
492 : {
493 12648 : const auto dim = _sys.get_mesh().mesh_dimension();
494 : // determine equation of plane: c_x * x + c_y * y + c_z * z + c = 0
495 2108 : std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
496 1054 : Real c = 0.;
497 :
498 : // We choose to constrain the dimension with the largest magnitude coefficient
499 : // This approach ensures the coefficients added to the constraint_row
500 : // (i.e., -c_xyz / c_max) have as small magnitude as possible
501 :
502 : // We initialize this to avoid maybe-uninitialized compiler error
503 1054 : unsigned int constrained_dim = 0;
504 :
505 : // Let's assert that we have a nonzero normal to ensure that constrained_dim
506 : // is always set
507 1054 : libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
508 :
509 1054 : Real max_abs_coef = 0.;
510 50592 : for (const auto d : make_range(dim))
511 : {
512 37944 : const auto coef = ref_normal_vec(d);
513 37944 : xyz_coefs.push_back(coef);
514 37944 : c -= coef * node(d);
515 :
516 3162 : const auto coef_abs = std::abs(coef);
517 37944 : if (coef_abs > max_abs_coef)
518 : {
519 1054 : max_abs_coef = coef_abs;
520 1054 : constrained_dim = d;
521 : }
522 : }
523 :
524 2108 : DofConstraintRow constraint_row;
525 50592 : for (const auto free_dim : make_range(dim))
526 : {
527 37944 : if (free_dim == constrained_dim)
528 12648 : continue;
529 :
530 25296 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
531 29512 : constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
532 : }
533 :
534 12648 : const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
535 12648 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
536 12648 : _sys.get_dof_map().add_constraint_row(
537 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
538 12648 : }
539 :
540 9264 : void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
541 : const Point & line_vec)
542 : {
543 9264 : const auto dim = _sys.get_mesh().mesh_dimension();
544 :
545 : // We will free the dimension most parallel to line_vec to keep the
546 : // constraint coefficients small
547 10036 : const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
548 9264 : auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
549 1544 : return std::abs(a) < std::abs(b);
550 1544 : });
551 9264 : const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
552 9264 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
553 :
554 : // A line is parameterized as r(t) = node + t * line_vec, so
555 : // x(t) = node(x) + t * line_vec(x)
556 : // y(t) = node(y) + t * line_vec(y)
557 : // z(t) = node(z) + t * line_vec(z)
558 : // Let's say we leave x free. Then t = (x(t) - node(x)) / line_vec(x)
559 : // Then y and z can be constrained as
560 : // y = node(y) + line_vec_y * (x(t) - node(x)) / line_vec(x)
561 : // = x(t) * line_vec(y) / line_vec(x) + (node(y) - node(x) * line_vec(y) / line_vec(x))
562 : // z = x(t) * line_vec(z) / line_vec(x) + (node(z) - node(x) * line_vec(z) / line_vec(x))
563 :
564 772 : libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
565 35160 : for (const auto constrained_dim : make_range(dim))
566 : {
567 25896 : if (constrained_dim == free_dim)
568 9264 : continue;
569 :
570 2772 : DofConstraintRow constraint_row;
571 16632 : constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
572 : const auto inhomogeneous_part =
573 16632 : node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
574 16632 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
575 16632 : _sys.get_dof_map().add_constraint_row(
576 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
577 : }
578 9264 : }
579 :
580 : // Utility function to determine whether two nodes share a boundary ID.
581 : // The motivation for this is that a sliding boundary node on a triangular
582 : // element can have a neighbor boundary node in the same element that is not
583 : // part of the same boundary
584 : // Consider the below example with nodes A, C, D, F, G that comprise elements
585 : // E1, E2, E3, with boundaries B1, B2, B3, B4. To determine the constraint
586 : // equations for the sliding node C, the neighboring nodes A and D need to
587 : // be identified to construct the line that C is allowed to slide along.
588 : // Note that neighbors A and D both share the boundary B1 with C.
589 : // Without ensuring that neighbors share the same boundary as the current
590 : // node, a neighboring node that does not lie on the same boundary
591 : // (i.e. F and G) might be selected to define the constraining line,
592 : // resulting in an incorrect constraint.
593 : // Note that if, for example, boundaries B1 and B2 were to be combined
594 : // into boundary B12, the node F would share a boundary id with node C
595 : // and result in an incorrect constraint. It would be useful to design
596 : // additional checks to detect cases like this.
597 :
598 : // B3
599 : // G-----------F
600 : // | \ / |
601 : // B4 | \ E2 / | B2
602 : // | \ / |
603 : // | E1 \ / E3 |
604 : // A-----C-----D
605 : // B1
606 :
607 177648 : bool VariationalSmootherConstraint::nodes_share_boundary_id(const Node & boundary_node,
608 : const Node & neighbor_node,
609 : const Elem & containing_elem,
610 : const BoundaryInfo & boundary_info)
611 : {
612 14804 : bool nodes_share_bid = false;
613 :
614 : // Node ids local to containing_elem
615 162844 : const auto node_id = containing_elem.get_node_index(&boundary_node);
616 162844 : const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
617 :
618 540600 : for (const auto side_id : containing_elem.side_index_range())
619 : {
620 : // We don't care about this side if it doesn't contain our boundary and neighbor nodes
621 820092 : if (!(containing_elem.is_node_on_side(node_id, side_id) &&
622 305580 : containing_elem.is_node_on_side(neighbor_id, side_id)))
623 269088 : continue;
624 :
625 : // If the current side, containing boundary_node and neighbor_node, lies on a boundary,
626 : // we can say that boundary_node and neighbor_node have a common boundary id.
627 20452 : std::vector<boundary_id_type> boundary_ids;
628 245424 : boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
629 245424 : if (boundary_ids.size())
630 : {
631 12630 : nodes_share_bid = true;
632 12630 : break;
633 : }
634 : }
635 177648 : return nodes_share_bid;
636 : }
637 :
638 : std::set<std::set<const Node *>>
639 2424 : VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
640 : const MeshBase & mesh,
641 : const Node & node,
642 : const subdomain_id_type sub_id,
643 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
644 : {
645 :
646 : // Find all the nodal neighbors... that is the nodes directly connected
647 : // to this node through one edge, or if none exists, use other nodes on the
648 : // containing face
649 404 : std::vector<const Node *> neighbors;
650 2424 : MeshTools::find_nodal_or_face_neighbors(
651 : mesh, node, nodes_to_elem_map, neighbors);
652 :
653 : // Each constituent set corresponds to neighbors sharing a face on the
654 : // subdomain boundary
655 202 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
656 :
657 12096 : for (const auto * neigh : neighbors)
658 : {
659 : // Determine whether the neighbor is on the subdomain boundary
660 : // First, find the common elements that both node and neigh belong to
661 9672 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
662 9672 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
663 806 : const Elem * common_elem = nullptr;
664 47415 : for (const auto * neigh_elem : elems_containing_neigh)
665 : {
666 43957 : if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
667 6364 : elems_containing_node.end())
668 : // We should be able to find a common element on the sub_id boundary
669 37743 : && (neigh_elem->subdomain_id() == sub_id))
670 838 : common_elem = neigh_elem;
671 : else
672 27509 : continue;
673 :
674 : // Now, determine whether node and neigh are on a side coincident
675 : // with the subdomain boundary
676 56124 : for (const auto common_side : common_elem->side_index_range())
677 : {
678 3746 : bool node_found_on_side = false;
679 3746 : bool neigh_found_on_side = false;
680 367136 : for (const auto local_node_id : common_elem->nodes_on_side(common_side))
681 : {
682 343326 : if (common_elem->node_id(local_node_id) == node.id())
683 1692 : node_found_on_side = true;
684 296696 : else if (common_elem->node_id(local_node_id) == neigh->id())
685 1974 : neigh_found_on_side = true;
686 : }
687 :
688 47140 : if (!(node_found_on_side && neigh_found_on_side &&
689 2500 : common_elem->neighbor_ptr(common_side)))
690 31328 : continue;
691 :
692 954 : const auto matched_side = common_side;
693 : // There could be multiple matched sides, so keep this next part
694 : // inside the common_side loop
695 : //
696 : // Does matched_side, containing both node and neigh, lie on the
697 : // sub_id subdomain boundary?
698 : const auto matched_neighbor_sub_id =
699 1908 : common_elem->neighbor_ptr(matched_side)->subdomain_id();
700 954 : const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
701 :
702 11770 : if (is_matched_side_on_subdomain_boundary)
703 : {
704 : // Store all nodes that live on this side
705 8670 : const auto nodes_on_side = common_elem->nodes_on_side(common_side);
706 1296 : std::set<const Node *> node_ptrs_on_side;
707 65080 : for (const auto local_node_id : nodes_on_side)
708 61650 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
709 8670 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
710 7374 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
711 :
712 648 : continue;
713 : }
714 :
715 : } // for common_side
716 :
717 : } // for neigh_elem
718 : }
719 :
720 2626 : return side_grouped_boundary_neighbors;
721 : }
722 :
723 : std::set<std::set<const Node *>>
724 25248 : VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
725 : const MeshBase & mesh,
726 : const Node & node,
727 : const std::unordered_set<dof_id_type> & boundary_node_ids,
728 : const BoundaryInfo & boundary_info,
729 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
730 : {
731 :
732 : // Find all the nodal neighbors... that is the nodes directly connected
733 : // to this node through one edge, or if none exists, use other nodes on the
734 : // containing face
735 4208 : std::vector<const Node *> neighbors;
736 25248 : MeshTools::find_nodal_or_face_neighbors(
737 : mesh, node, nodes_to_elem_map, neighbors);
738 :
739 : // Each constituent set corresponds to neighbors sharing a face on the
740 : // boundary
741 2104 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
742 :
743 141192 : for (const auto * neigh : neighbors)
744 : {
745 : const bool is_neighbor_boundary_node =
746 115944 : boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
747 115944 : if (!is_neighbor_boundary_node)
748 20152 : continue;
749 :
750 : // Determine whether nodes share a common boundary id
751 : // First, find the common element that both node and neigh belong to
752 93960 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
753 93960 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
754 7830 : const Elem * common_elem = nullptr;
755 565736 : for (const auto * neigh_elem : elems_containing_neigh)
756 : {
757 : const bool is_neigh_common =
758 511994 : std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
759 80436 : elems_containing_node.end();
760 471776 : if (!is_neigh_common)
761 268714 : continue;
762 29608 : common_elem = neigh_elem;
763 : // Keep this in the neigh_elem loop because there can be multiple common
764 : // elements Now, determine whether node and neigh share a common boundary
765 : // id
766 177648 : const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
767 : node, *neigh, *common_elem, boundary_info);
768 177648 : if (nodes_have_common_bid)
769 : {
770 : // Find the side coinciding with the shared boundary
771 885528 : for (const auto side : common_elem->side_index_range())
772 : {
773 : // We only care about external boundaries here, make sure side doesn't
774 : // have a neighbor
775 795132 : if (common_elem->neighbor_ptr(side))
776 570912 : continue;
777 :
778 20412 : bool node_found_on_side = false;
779 20412 : bool neigh_found_on_side = false;
780 244944 : const auto nodes_on_side = common_elem->nodes_on_side(side);
781 1754256 : for (const auto local_node_id : nodes_on_side)
782 : {
783 1635088 : if (common_elem->node_id(local_node_id) == node.id())
784 14712 : node_found_on_side = true;
785 1332768 : else if (common_elem->node_id(local_node_id) == neigh->id())
786 15396 : neigh_found_on_side = true;
787 : }
788 244944 : if (!(node_found_on_side && neigh_found_on_side))
789 6824 : continue;
790 :
791 27176 : std::set<const Node *> node_ptrs_on_side;
792 1089744 : for (const auto local_node_id : nodes_on_side)
793 1003912 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
794 176644 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
795 149468 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
796 : }
797 138930 : continue;
798 126300 : }
799 : }
800 : }
801 :
802 27352 : return side_grouped_boundary_neighbors;
803 : }
804 :
805 27672 : ConstraintVariant VariationalSmootherConstraint::determine_constraint(
806 : const Node & node,
807 : const unsigned int dim,
808 : const std::set<std::set<const Node *>> & side_grouped_boundary_neighbors)
809 : {
810 : // Determines the appropriate geometric constraint for a node based on its
811 : // neighbors.
812 :
813 : // Extract neighbors in flat vector
814 4612 : std::vector<const Node *> neighbors;
815 108902 : for (const auto & side : side_grouped_boundary_neighbors)
816 81230 : neighbors.insert(neighbors.end(), side.begin(), side.end());
817 :
818 : // Constrain the node to it's current location
819 : // Note that neighbors.size() <= 1 is used here instead of == 1 so we don't
820 : // crash when trying to access *neighbors[0] a few lines below in the
821 : // dim == 2 case
822 27672 : if (dim == 1 || neighbors.size() <= 1)
823 490 : return PointConstraint{node};
824 :
825 27252 : if (dim == 2 || neighbors.size() == 2)
826 : {
827 : // Determine whether node + all neighbors are colinear
828 263 : bool all_colinear = true;
829 3419 : const Point ref_dir = (*neighbors[0] - node).unit();
830 6171 : for (const auto i : make_range(std::size_t(1), neighbors.size()))
831 : {
832 4059 : const Point delta = *(neighbors[i]) - node;
833 312 : libmesh_assert(delta.norm() >= TOLERANCE);
834 3747 : const Point dir = delta.unit();
835 4022 : if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
836 : {
837 61 : all_colinear = false;
838 732 : break;
839 : }
840 : }
841 :
842 263 : if (all_colinear)
843 2828 : return LineConstraint{node, ref_dir};
844 :
845 854 : return PointConstraint{node};
846 : }
847 :
848 : // dim == 3, neighbors.size() >= 3
849 4016 : std::set<PlaneConstraint> valid_planes;
850 99040 : for (const auto & side_nodes : side_grouped_boundary_neighbors)
851 : {
852 81180 : std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
853 422916 : for (const auto i : index_range(side_nodes_vec))
854 : {
855 376900 : const Point vec_i = (*side_nodes_vec[i] - node);
856 1183902 : for (const auto j : make_range(i))
857 : {
858 905362 : const Point vec_j = (*side_nodes_vec[j] - node);
859 766498 : Point candidate_normal = vec_i.cross(vec_j);
860 835930 : if (candidate_normal.norm() <= TOLERANCE)
861 54352 : continue;
862 :
863 781578 : const PlaneConstraint candidate_plane{node, candidate_normal.unit()};
864 716661 : valid_planes.emplace(candidate_plane);
865 : }
866 : }
867 : }
868 :
869 : // Fall back to point constraint
870 24096 : if (valid_planes.empty())
871 0 : return PointConstraint(node);
872 :
873 : // Combine all the planes together to get a common constraint
874 2008 : auto it = valid_planes.begin();
875 4016 : ConstraintVariant current = *it++;
876 51368 : for (; it != valid_planes.end(); ++it)
877 : {
878 29530 : current = intersect_constraints(current, *it);
879 :
880 : // This will catch cases where constraints have no intersection
881 : // (i.e., the element surface is non-planar)
882 : // Fall back to fixed node constraint
883 27272 : if (std::holds_alternative<InvalidConstraint>(current))
884 : {
885 0 : current = PointConstraint(node);
886 0 : break;
887 : }
888 : }
889 :
890 2008 : return current;
891 : }
892 :
893 : // Applies the computed constraint (PointConstraint, LineConstraint, or
894 : // PlaneConstraint) to a node.
895 26160 : void VariationalSmootherConstraint::impose_constraint(const Node & node,
896 : const ConstraintVariant & constraint)
897 : {
898 :
899 2180 : libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
900 : "Cannot impose constraint using InvalidConstraint.");
901 :
902 26160 : if (std::holds_alternative<PointConstraint>(constraint))
903 4248 : fix_node(node);
904 21912 : else if (std::holds_alternative<LineConstraint>(constraint))
905 9264 : constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
906 12648 : else if (std::holds_alternative<PlaneConstraint>(constraint))
907 12648 : constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
908 :
909 : else
910 0 : libmesh_assert_msg(false, "Unknown constraint type.");
911 26160 : }
912 :
913 : } // namespace libMesh
|
