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 : // 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 vector in the positive x/y/z direction
27 8184809 : auto get_positive_vector = [](const Point & vec) -> Point {
28 8184809 : libmesh_error_msg_if(vec.norm() < TOLERANCE,
29 : "Can't define a positively-oriented vector with a zero vector.");
30 :
31 : // Choose sign such that direction vector points in the positive x/y/z
32 : // direction This helps to eliminate duplicate lines/planes
33 230558 : Point canonical{0, 0, 0};
34 : // Choose the canonical dimension to ensure the dot product below is nonzero
35 15283957 : for (const auto dim_id : make_range(3))
36 15714491 : if (!absolute_fuzzy_equals(vec(dim_id), 0.))
37 : {
38 8184809 : canonical(dim_id) = 1.;
39 7954251 : break;
40 : }
41 :
42 230558 : const auto dot_prod = vec * canonical;
43 230558 : libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));
44 :
45 12317843 : return (dot_prod > 0) ? vec.unit() : -vec.unit();
46 : };
47 :
48 13951004 : PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)
49 : {
50 13951004 : }
51 :
52 0 : bool PointConstraint::operator<(const PointConstraint & other) const
53 : {
54 0 : if (*this == other)
55 0 : return false;
56 :
57 0 : return _point < other.point();
58 : }
59 :
60 0 : bool PointConstraint::operator==(const PointConstraint & other) const
61 : {
62 0 : return _point.absolute_fuzzy_equals(other.point(), _tol);
63 : }
64 :
65 937727 : ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const
66 : {
67 : // using visit to resolve the variant to its actual type
68 : return std::visit(
69 2707499 : [&](auto && o) -> ConstraintVariant {
70 937727 : if (!o.contains_point(*this))
71 : // Point is not on the constraint
72 0 : return InvalidConstraint();
73 :
74 26423 : return *this;
75 : },
76 937727 : other);
77 : }
78 :
79 59285 : LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
80 59285 : : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
81 : {
82 59285 : libmesh_error_msg_if(_direction.norm() < _tol,
83 : "Can't define a line with zero magnitude direction vector.");
84 59285 : }
85 :
86 0 : bool LineConstraint::operator<(const LineConstraint & other) const
87 : {
88 0 : if (*this == other)
89 0 : return false;
90 :
91 0 : if (!(_direction.absolute_fuzzy_equals(other.direction(), _tol)))
92 0 : return _direction < other.direction();
93 0 : return (_direction * _point) < (other.direction() * other.point());
94 : }
95 :
96 0 : bool LineConstraint::operator==(const LineConstraint & other) const
97 : {
98 0 : if (!(_direction.absolute_fuzzy_equals(other.direction(), _tol)))
99 0 : return false;
100 0 : return this->contains_point(other.point());
101 : }
102 :
103 1420 : bool LineConstraint::contains_point(const PointConstraint & p) const
104 : {
105 : // If the point lies on the line, then the vector p - point is parallel to the
106 : // line In that case, the cross product of p - point with the line's direction
107 : // will be zero.
108 1420 : return _direction.cross(p.point() - _point).norm() < _tol;
109 : }
110 :
111 0 : bool LineConstraint::is_parallel(const LineConstraint & l) const
112 : {
113 0 : return _direction.absolute_fuzzy_equals(l.direction(), _tol);
114 : }
115 :
116 25773 : bool LineConstraint::is_parallel(const PlaneConstraint & p) const
117 : {
118 25773 : return _direction * p.normal() < _tol;
119 : }
120 :
121 49478 : ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const
122 : {
123 : // using visit to resolve the variant to its actual type
124 : return std::visit(
125 142860 : [&](auto && o) -> ConstraintVariant {
126 : // Use std::decay_t to strip references/const/volatile from o's type,
127 : // so we can match the actual stored variant type in std::is_same_v.
128 : using T = std::decay_t<decltype(o)>;
129 : if constexpr (std::is_same_v<T, LineConstraint>)
130 : {
131 49478 : if (*this == o)
132 0 : return *this;
133 :
134 0 : if (this->is_parallel(o))
135 : // Lines are parallel and do not intersect
136 0 : return InvalidConstraint();
137 :
138 : // Solve for t in the equation p1 + t·d1 = p2 + s·d2
139 : // The shortest vector between skew lines lies along the normal vector
140 : // (d1 × d2). Projecting the vector (p2 - p1) onto this normal gives a
141 : // scalar proportional to the distance. This is equivalent to solving:
142 : // ((p2 - p1) × d2) · (d1 × d2) = t · |d1 × d2|²
143 : // ⇒ t = ((delta × d2) · (d1 × d2)) / |d1 × d2|²
144 :
145 0 : const Point delta = o.point() - _point;
146 0 : const Point cross_d1_d2 = _direction.cross(o.direction());
147 0 : const Real cross_dot = (delta.cross(o.direction())) * cross_d1_d2;
148 0 : const Real denom = cross_d1_d2.norm_sq();
149 :
150 0 : const Real t = cross_dot / denom;
151 0 : const Point intersection = _point + t * _direction;
152 :
153 : // Verify that intersection lies on both lines
154 0 : if (o.direction().cross(intersection - o.point()).norm() > _tol)
155 : // Lines do not intersect at a single point
156 0 : return InvalidConstraint();
157 :
158 0 : return PointConstraint{intersection};
159 : }
160 :
161 : else if constexpr (std::is_same_v<T, PlaneConstraint>)
162 97563 : return o.intersect(*this);
163 :
164 : else if constexpr (std::is_same_v<T, PointConstraint>)
165 : {
166 0 : if (!this->contains_point(o))
167 : // Point is not on the line
168 0 : return InvalidConstraint();
169 :
170 0 : return o;
171 : }
172 :
173 : else
174 0 : libmesh_error_msg("Unsupported constraint type in Line::intersect.");
175 : },
176 49478 : other);
177 : }
178 :
179 8125524 : PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
180 8125524 : : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
181 : {
182 8125524 : libmesh_error_msg_if(_normal.norm() < _tol,
183 : "Can't define a plane with zero magnitude direction vector.");
184 8125524 : }
185 :
186 22103643 : bool PlaneConstraint::operator<(const PlaneConstraint & other) const
187 : {
188 22103643 : if (*this == other)
189 390759 : return false;
190 :
191 8464274 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
192 8231582 : return _normal < other.normal();
193 0 : return (_normal * _point) < (other.normal() * other.point());
194 : }
195 :
196 22155828 : bool PlaneConstraint::operator==(const PlaneConstraint & other) const
197 : {
198 22780745 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
199 232994 : return false;
200 13872061 : return this->contains_point(other.point());
201 : }
202 :
203 52185 : bool PlaneConstraint::is_parallel(const PlaneConstraint & p) const
204 : {
205 52185 : return _normal.absolute_fuzzy_equals(p.normal(), _tol);
206 : }
207 :
208 25773 : bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }
209 :
210 14857846 : bool PlaneConstraint::contains_point(const PointConstraint & p) const
211 : {
212 : // distance between the point and the plane
213 418535 : const Real dist = (p.point() - _point) * _normal;
214 14857846 : return std::abs(dist) < _tol;
215 : }
216 :
217 49478 : bool PlaneConstraint::contains_line(const LineConstraint & l) const
218 : {
219 49478 : const bool base_on_plane = this->contains_point(PointConstraint(l.point()));
220 49478 : const bool dir_orthogonal = std::abs(_normal * l.direction()) < _tol;
221 49478 : return base_on_plane && dir_orthogonal;
222 : }
223 :
224 101663 : ConstraintVariant PlaneConstraint::intersect(const ConstraintVariant & other) const
225 : {
226 : // using visit to resolve the variant to its actual type
227 : return std::visit(
228 293535 : [&](auto && o) -> ConstraintVariant {
229 : // Use std::decay_t to strip references/const/volatile from o's type,
230 : // so we can match the actual stored variant type in std::is_same_v.
231 : using T = std::decay_t<decltype(o)>;
232 : if constexpr (std::is_same_v<T, PlaneConstraint>)
233 : {
234 : // If planes are identical, return one of them
235 83648 : if (*this == o)
236 0 : return *this;
237 :
238 52185 : if (this->is_parallel(o))
239 : // Planes are parallel and do not intersect
240 0 : return InvalidConstraint();
241 :
242 : // Solve for a point on the intersection line of two planes.
243 : // Given planes:
244 : // Plane 1: n1 · (x - p1) = 0
245 : // Plane 2: n2 · (x - p2) = 0
246 : // The line of intersection has direction dir = n1 × n2.
247 : // To find a point on this line, we assume:
248 : // x = p1 + s·n1 = p2 + t·n2
249 : // ⇒ p1 - p2 = t·n2 - s·n1
250 : // Taking dot products with n1 and n2 leads to:
251 : // [-n1·n1 n1·n2] [s] = [n1 · (p1 - p2)]
252 : // [-n1·n2 n2·n2] [t] [n2 · (p1 - p2)]
253 :
254 50715 : const Point dir = this->_normal.cross(o.normal()); // direction of line of intersection
255 1470 : libmesh_assert(dir.norm() > _tol);
256 1470 : const Point w = _point - o.point();
257 :
258 : // Dot product terms used in 2x2 system
259 1470 : const Real n1_dot_n1 = _normal * _normal;
260 1470 : const Real n1_dot_n2 = _normal * o.normal();
261 1470 : const Real n2_dot_n2 = o.normal() * o.normal();
262 1470 : const Real n1_dot_w = _normal * w;
263 1470 : const Real n2_dot_w = o.normal() * w;
264 :
265 52185 : const Real denom = -(n1_dot_n1 * n2_dot_n2 - n1_dot_n2 * n1_dot_n2);
266 1470 : libmesh_assert(std::abs(denom) > _tol);
267 :
268 52185 : const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
269 1470 : const Point p0 = _point + s * _normal;
270 :
271 52185 : return LineConstraint{p0, dir};
272 : }
273 :
274 : else if constexpr (std::is_same_v<T, LineConstraint>)
275 : {
276 49478 : if (this->contains_line(o))
277 667 : return o;
278 :
279 25773 : if (this->is_parallel(o))
280 : // Line is parallel and does not intersect the plane
281 71 : return InvalidConstraint();
282 :
283 : // Solve for t in the parametric equation:
284 : // p(t) = point + t·d
285 : // such that this point also satisfies the plane equation:
286 : // n · (p(t) - p0) = 0
287 : // which leads to:
288 : // t = (n · (p0 - point)) / (n · d)
289 :
290 724 : const Real denom = _normal * o.direction();
291 724 : libmesh_assert(std::abs(denom) > _tol);
292 25702 : const Real t = (_normal * (_point - o.point())) / denom;
293 25702 : return PointConstraint{o.point() + t * o.direction()};
294 : }
295 :
296 : else if constexpr (std::is_same_v<T, PointConstraint>)
297 : {
298 0 : if (!this->contains_point(o))
299 : // Point is not on the plane
300 0 : return InvalidConstraint();
301 :
302 0 : return o;
303 : }
304 :
305 : else
306 0 : libmesh_error_msg("Unsupported constraint type in Plane::intersect.");
307 : },
308 101663 : other);
309 : }
310 :
311 852 : VariationalSmootherConstraint::VariationalSmootherConstraint(
312 852 : System & sys, const bool & preserve_subdomain_boundaries)
313 852 : : Constraint(), _sys(sys), _preserve_subdomain_boundaries(preserve_subdomain_boundaries)
314 : {
315 852 : }
316 :
317 804 : VariationalSmootherConstraint::~VariationalSmootherConstraint() = default;
318 :
319 852 : void VariationalSmootherConstraint::constrain()
320 : {
321 852 : const auto & mesh = _sys.get_mesh();
322 852 : const auto dim = mesh.mesh_dimension();
323 :
324 : // Only compute the node to elem map once
325 48 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
326 852 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
327 :
328 24 : const auto & boundary_info = mesh.get_boundary_info();
329 :
330 876 : const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
331 :
332 : // Identify/constrain subdomain boundary nodes, if requested
333 48 : std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
334 852 : if (_preserve_subdomain_boundaries)
335 : {
336 148255 : for (const auto * elem : mesh.active_element_ptr_range())
337 : {
338 2860 : const auto & sub_id1 = elem->subdomain_id();
339 312045 : for (const auto side : elem->side_index_range())
340 : {
341 261280 : const auto * neighbor = elem->neighbor_ptr(side);
342 261280 : if (neighbor == nullptr)
343 48714 : continue;
344 :
345 11896 : const auto & sub_id2 = neighbor->subdomain_id();
346 211154 : if (sub_id1 == sub_id2)
347 191406 : continue;
348 :
349 : // elem and neighbor are in different subdomains, and share nodes
350 : // that need to be constrained
351 117408 : for (const auto local_node_id : elem->nodes_on_side(side))
352 : {
353 105704 : const auto & node = mesh.node_ref(elem->node_id(local_node_id));
354 : // Make sure we haven't already processed this node
355 102808 : if (subdomain_boundary_map.count(node.id()))
356 77461 : continue;
357 :
358 : // Get the relevant nodal neighbors for the subdomain constraint
359 : const auto side_grouped_boundary_neighbors =
360 : get_neighbors_for_subdomain_constraint(
361 26061 : mesh, node, sub_id1, nodes_to_elem_map);
362 :
363 : // Determine which constraint should be imposed
364 : const auto subdomain_constraint =
365 25347 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
366 :
367 : // This subdomain boundary node does not lie on an external boundary,
368 : // go ahead and impose constraint
369 26061 : if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
370 16685 : this->impose_constraint(node, subdomain_constraint);
371 :
372 : // This subdomain boundary node could lie on an external boundary, save it
373 : // for later to combine with the external boundary constraint.
374 : // We also save constraints for non-boundary nodes so we don't try to
375 : // re-constrain the node when accessed from the neighboring elem.
376 : // See subdomain_boundary_map.count call above.
377 49980 : subdomain_boundary_map[node.id()] = subdomain_constraint;
378 :
379 : } // for local_node_id
380 :
381 : } // for side
382 804 : } // for elem
383 : }
384 :
385 : // Loop through boundary nodes and impose constraints
386 138166 : for (const auto & bid : boundary_node_ids)
387 : {
388 137314 : const auto & node = mesh.node_ref(bid);
389 :
390 : // Get the relevant nodal neighbors for the boundary constraint
391 : const auto side_grouped_boundary_neighbors = get_neighbors_for_boundary_constraint(
392 141182 : mesh, node, boundary_node_ids, boundary_info, nodes_to_elem_map);
393 :
394 : // Determine which constraint should be imposed
395 : const auto boundary_constraint =
396 141182 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
397 :
398 : // Check for the case where this boundary node is also part of a subdomain id boundary
399 137314 : if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
400 : {
401 8662 : const auto & subdomain_constraint = it->second;
402 : // Combine current boundary constraint with previously determined
403 : // subdomain_constraint
404 488 : auto combined_constraint = intersect_constraints(subdomain_constraint, boundary_constraint);
405 :
406 : // This will catch cases where constraints have no intersection
407 : // Fall back to fixed node constraint
408 8662 : if (std::holds_alternative<InvalidConstraint>(combined_constraint))
409 0 : combined_constraint = PointConstraint(node);
410 :
411 8662 : this->impose_constraint(node, combined_constraint);
412 : }
413 :
414 : else
415 128652 : this->impose_constraint(node, boundary_constraint);
416 :
417 : } // end bid
418 852 : }
419 :
420 29465 : void VariationalSmootherConstraint::fix_node(const Node & node)
421 : {
422 113671 : for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
423 : {
424 84206 : const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
425 4744 : DofConstraintRow constraint_row;
426 : // Leave the constraint row as all zeros so this dof is independent from other dofs
427 84206 : const auto constrained_value = node(d);
428 : // Simply constrain this dof to retain it's current value
429 84206 : _sys.get_dof_map().add_constraint_row(
430 : constrained_dof_index, constraint_row, constrained_value, true);
431 : }
432 29465 : }
433 :
434 92442 : void VariationalSmootherConstraint::constrain_node_to_plane(const Node & node,
435 : const Point & ref_normal_vec)
436 : {
437 92442 : const auto dim = _sys.get_mesh().mesh_dimension();
438 : // determine equation of plane: c_x * x + c_y * y + c_z * z + c = 0
439 5208 : std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
440 2604 : Real c = 0.;
441 :
442 : // We choose to constrain the dimension with the largest magnitude coefficient
443 : // This approach ensures the coefficients added to the constraint_row
444 : // (i.e., -c_xyz / c_max) have as small magnitude as possible
445 :
446 : // We initialize this to avoid maybe-uninitialized compiler error
447 2604 : unsigned int constrained_dim = 0;
448 :
449 : // Let's assert that we have a nonzero normal to ensure that constrained_dim
450 : // is always set
451 2604 : libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
452 :
453 2604 : Real max_abs_coef = 0.;
454 369768 : for (const auto d : make_range(dim))
455 : {
456 277326 : const auto coef = ref_normal_vec(d);
457 277326 : xyz_coefs.push_back(coef);
458 277326 : c -= coef * node(d);
459 :
460 7812 : const auto coef_abs = std::abs(coef);
461 277326 : if (coef_abs > max_abs_coef)
462 : {
463 2604 : max_abs_coef = coef_abs;
464 2604 : constrained_dim = d;
465 : }
466 : }
467 :
468 5208 : DofConstraintRow constraint_row;
469 369768 : for (const auto free_dim : make_range(dim))
470 : {
471 277326 : if (free_dim == constrained_dim)
472 92442 : continue;
473 :
474 184884 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
475 195300 : constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
476 : }
477 :
478 92442 : const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
479 92442 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
480 92442 : _sys.get_dof_map().add_constraint_row(
481 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
482 92442 : }
483 :
484 32092 : void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
485 : const Point & line_vec)
486 : {
487 32092 : const auto dim = _sys.get_mesh().mesh_dimension();
488 :
489 : // We will free the dimension most parallel to line_vec to keep the
490 : // constraint coefficients small
491 32996 : const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
492 32092 : auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
493 1808 : return std::abs(a) < std::abs(b);
494 1808 : });
495 32092 : const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
496 32092 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
497 :
498 : // A line is parameterized as r(t) = node + t * line_vec, so
499 : // x(t) = node(x) + t * line_vec(x)
500 : // y(t) = node(y) + t * line_vec(y)
501 : // z(t) = node(z) + t * line_vec(z)
502 : // Let's say we leave x free. Then t = (x(t) - node(x)) / line_vec(x)
503 : // Then y and z can be constrained as
504 : // y = node(y) + line_vec_y * (x(t) - node(x)) / line_vec(x)
505 : // = x(t) * line_vec(y) / line_vec(x) + (node(y) - node(x) * line_vec(y) / line_vec(x))
506 : // z = x(t) * line_vec(z) / line_vec(x) + (node(z) - node(x) * line_vec(z) / line_vec(x))
507 :
508 904 : libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
509 121836 : for (const auto constrained_dim : make_range(dim))
510 : {
511 89744 : if (constrained_dim == free_dim)
512 32092 : continue;
513 :
514 3248 : DofConstraintRow constraint_row;
515 57652 : constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
516 : const auto inhomogeneous_part =
517 57652 : node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
518 57652 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
519 57652 : _sys.get_dof_map().add_constraint_row(
520 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
521 : }
522 32092 : }
523 :
524 162661 : void VariationalSmootherConstraint::find_nodal_or_face_neighbors(
525 : const MeshBase & mesh,
526 : const Node & node,
527 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
528 : std::vector<const Node *> & neighbors)
529 : {
530 : // Find all the nodal neighbors... that is the nodes directly connected
531 : // to this node through one edge.
532 162661 : MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
533 :
534 : // If no neighbors are found, use all nodes on the containing side as
535 : // neighbors.
536 162661 : if (!neighbors.size())
537 : {
538 : // Grab the element containing node
539 26625 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
540 : // Find the element side containing node
541 97625 : for (const auto &side : elem->side_index_range())
542 : {
543 97625 : const auto &nodes_on_side = elem->nodes_on_side(side);
544 : const auto it =
545 116875 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
546 301125 : return elem->node_id(local_node_id) == node.id();
547 5500 : });
548 :
549 100375 : if (it != nodes_on_side.end())
550 : {
551 266250 : for (const auto &local_node_id : nodes_on_side)
552 : // No need to add node itself as a neighbor
553 246375 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
554 6750 : *node_ptr != node)
555 213000 : neighbors.push_back(node_ptr);
556 750 : break;
557 : }
558 : }
559 : }
560 4582 : libmesh_assert(neighbors.size());
561 162661 : }
562 :
563 : // Utility function to determine whether two nodes share a boundary ID.
564 : // The motivation for this is that a sliding boundary node on a triangular
565 : // element can have a neighbor boundary node in the same element that is not
566 : // part of the same boundary
567 : // Consider the below example with nodes A, C, D, F, G that comprise elements
568 : // E1, E2, E3, with boundaries B1, B2, B3, B4. To determine the constraint
569 : // equations for the sliding node C, the neighboring nodes A and D need to
570 : // be identified to construct the line that C is allowed to slide along.
571 : // Note that neighbors A and D both share the boundary B1 with C.
572 : // Without ensuring that neighbors share the same boundary as the current
573 : // node, a neighboring node that does not lie on the same boundary
574 : // (i.e. F and G) might be selected to define the constraining line,
575 : // resulting in an incorrect constraint.
576 : // Note that if, for example, boundaries B1 and B2 were to be combined
577 : // into boundary B12, the node F would share a boundary id with node C
578 : // and result in an incorrect constraint. It would be useful to design
579 : // additional checks to detect cases like this.
580 :
581 : // B3
582 : // G-----------F
583 : // | \ / |
584 : // B4 | \ E2 / | B2
585 : // | \ / |
586 : // | E1 \ / E3 |
587 : // A-----C-----D
588 : // B1
589 :
590 626504 : bool VariationalSmootherConstraint::nodes_share_boundary_id(const Node & boundary_node,
591 : const Node & neighbor_node,
592 : const Elem & containing_elem,
593 : const BoundaryInfo & boundary_info)
594 : {
595 17648 : bool nodes_share_bid = false;
596 :
597 : // Node ids local to containing_elem
598 608856 : const auto node_id = containing_elem.get_node_index(&boundary_node);
599 608856 : const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
600 :
601 2122261 : for (const auto side_id : containing_elem.side_index_range())
602 : {
603 : // We don't care about this side if it doesn't contain our boundary and neighbor nodes
604 3099505 : if (!(containing_elem.is_node_on_side(node_id, side_id) &&
605 978948 : containing_elem.is_node_on_side(neighbor_id, side_id)))
606 1298093 : continue;
607 :
608 : // If the current side, containing boundary_node and neighbor_node, lies on a boundary,
609 : // we can say that boundary_node and neighbor_node have a common boundary id.
610 23168 : std::vector<boundary_id_type> boundary_ids;
611 822464 : boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
612 822464 : if (boundary_ids.size())
613 : {
614 17600 : nodes_share_bid = true;
615 17600 : break;
616 : }
617 : }
618 626504 : return nodes_share_bid;
619 : }
620 :
621 : std::set<std::set<const Node *>>
622 25347 : VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
623 : const MeshBase & mesh,
624 : const Node & node,
625 : const subdomain_id_type sub_id,
626 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
627 : {
628 :
629 : // Find all the nodal neighbors... that is the nodes directly connected
630 : // to this node through one edge, or if none exists, use other nodes on the
631 : // containing face
632 1428 : std::vector<const Node *> neighbors;
633 25347 : VariationalSmootherConstraint::find_nodal_or_face_neighbors(
634 : mesh, node, nodes_to_elem_map, neighbors);
635 :
636 : // Each constituent set corresponds to neighbors sharing a face on the
637 : // subdomain boundary
638 714 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
639 :
640 125315 : for (const auto * neigh : neighbors)
641 : {
642 : // Determine whether the neighbor is on the subdomain boundary
643 : // First, find the common elements that both node and neigh belong to
644 99968 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
645 99968 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
646 2816 : const Elem * common_elem = nullptr;
647 634598 : for (const auto * neigh_elem : elems_containing_neigh)
648 : {
649 636456 : if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
650 30120 : elems_containing_node.end())
651 : // We should be able to find a common element on the sub_id boundary
652 534630 : && (neigh_elem->subdomain_id() == sub_id))
653 3388 : common_elem = neigh_elem;
654 : else
655 414356 : continue;
656 :
657 : // Now, determine whether node and neigh are on a side coincident
658 : // with the subdomain boundary
659 820192 : for (const auto common_side : common_elem->side_index_range())
660 : {
661 19716 : bool node_found_on_side = false;
662 19716 : bool neigh_found_on_side = false;
663 6833870 : for (const auto local_node_id : common_elem->nodes_on_side(common_side))
664 : {
665 6286468 : if (common_elem->node_id(local_node_id) == node.id())
666 6876 : node_found_on_side = true;
667 5870138 : else if (common_elem->node_id(local_node_id) == neigh->id())
668 9334 : neigh_found_on_side = true;
669 : }
670 :
671 705266 : if (!(node_found_on_side && neigh_found_on_side &&
672 10696 : common_elem->neighbor_ptr(common_side)))
673 515154 : continue;
674 :
675 4784 : const auto matched_side = common_side;
676 : // There could be multiple matched sides, so keep this next part
677 : // inside the common_side loop
678 : //
679 : // Does matched_side, containing both node and neigh, lie on the
680 : // sub_id subdomain boundary?
681 : const auto matched_neighbor_sub_id =
682 9568 : common_elem->neighbor_ptr(matched_side)->subdomain_id();
683 4784 : const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
684 :
685 169832 : if (is_matched_side_on_subdomain_boundary)
686 : {
687 : // Store all nodes that live on this side
688 108040 : const auto nodes_on_side = common_elem->nodes_on_side(common_side);
689 5920 : std::set<const Node *> node_ptrs_on_side;
690 1028932 : for (const auto local_node_id : nodes_on_side)
691 949876 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
692 108040 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
693 102120 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
694 :
695 2960 : continue;
696 : }
697 :
698 : } // for common_side
699 :
700 : } // for neigh_elem
701 : }
702 :
703 26061 : return side_grouped_boundary_neighbors;
704 : }
705 :
706 : std::set<std::set<const Node *>>
707 137314 : VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
708 : const MeshBase & mesh,
709 : const Node & node,
710 : const std::unordered_set<dof_id_type> & boundary_node_ids,
711 : const BoundaryInfo & boundary_info,
712 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
713 : {
714 :
715 : // Find all the nodal neighbors... that is the nodes directly connected
716 : // to this node through one edge, or if none exists, use other nodes on the
717 : // containing face
718 7736 : std::vector<const Node *> neighbors;
719 137314 : VariationalSmootherConstraint::find_nodal_or_face_neighbors(
720 : mesh, node, nodes_to_elem_map, neighbors);
721 :
722 : // Each constituent set corresponds to neighbors sharing a face on the
723 : // boundary
724 3868 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
725 :
726 603358 : for (const auto * neigh : neighbors)
727 : {
728 : const bool is_neighbor_boundary_node =
729 466044 : boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
730 466044 : if (!is_neighbor_boundary_node)
731 35328 : continue;
732 :
733 : // Determine whether nodes share a common boundary id
734 : // First, find the common element that both node and neigh belong to
735 429692 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
736 429692 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
737 12104 : const Elem * common_elem = nullptr;
738 1688593 : for (const auto * neigh_elem : elems_containing_neigh)
739 : {
740 : const bool is_neigh_common =
741 1294363 : std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
742 70924 : elems_containing_node.end();
743 1258901 : if (!is_neigh_common)
744 614583 : continue;
745 35296 : common_elem = neigh_elem;
746 : // Keep this in the neigh_elem loop because there can be multiple common
747 : // elements Now, determine whether node and neigh share a common boundary
748 : // id
749 626504 : const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
750 : node, *neigh, *common_elem, boundary_info);
751 626504 : if (nodes_have_common_bid)
752 : {
753 : // Find the side coinciding with the shared boundary
754 4332420 : for (const auto side : common_elem->side_index_range())
755 : {
756 : // We only care about external boundaries here, make sure side doesn't
757 : // have a neighbor
758 3812060 : if (common_elem->neighbor_ptr(side))
759 3035960 : continue;
760 :
761 30584 : bool node_found_on_side = false;
762 30584 : bool neigh_found_on_side = false;
763 1085732 : const auto nodes_on_side = common_elem->nodes_on_side(side);
764 9849972 : for (const auto local_node_id : nodes_on_side)
765 : {
766 9011120 : if (common_elem->node_id(local_node_id) == node.id())
767 20704 : node_found_on_side = true;
768 8029248 : else if (common_elem->node_id(local_node_id) == neigh->id())
769 22890 : neigh_found_on_side = true;
770 : }
771 1085732 : if (!(node_found_on_side && neigh_found_on_side))
772 11664 : continue;
773 :
774 37840 : std::set<const Node *> node_ptrs_on_side;
775 6052040 : for (const auto local_node_id : nodes_on_side)
776 5531940 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
777 690580 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
778 652740 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
779 : }
780 607200 : continue;
781 589600 : }
782 : }
783 : }
784 :
785 141182 : return side_grouped_boundary_neighbors;
786 : }
787 :
788 162661 : ConstraintVariant VariationalSmootherConstraint::determine_constraint(
789 : const Node & node,
790 : const unsigned int dim,
791 : const std::set<std::set<const Node *>> & side_grouped_boundary_neighbors)
792 : {
793 : // Determines the appropriate geometric constraint for a node based on its
794 : // neighbors.
795 :
796 : // Extract neighbors in flat vector
797 9164 : std::vector<const Node *> neighbors;
798 542653 : for (const auto & side : side_grouped_boundary_neighbors)
799 379992 : neighbors.insert(neighbors.end(), side.begin(), side.end());
800 :
801 : // Constrain the node to it's current location
802 162661 : if (dim == 1 || neighbors.size() == 1)
803 1125 : return PointConstraint{node};
804 :
805 161596 : if (dim == 2 || neighbors.size() == 2)
806 : {
807 : // Determine whether node + all neighbors are colinear
808 274 : bool all_colinear = true;
809 10001 : const Point ref_dir = (*neighbors[0] - node).unit();
810 19383 : for (const auto i : make_range(std::size_t(1), neighbors.size()))
811 : {
812 12629 : const Point delta = *(neighbors[i]) - node;
813 346 : libmesh_assert(delta.norm() >= TOLERANCE);
814 12283 : const Point dir = delta.unit();
815 12589 : if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
816 : {
817 74 : all_colinear = false;
818 2627 : break;
819 : }
820 : }
821 :
822 274 : if (all_colinear)
823 7500 : return LineConstraint{node, ref_dir};
824 :
825 2775 : return PointConstraint{node};
826 : }
827 :
828 : // dim == 3, neighbors.size() >= 3
829 8556 : std::set<PlaneConstraint> valid_planes;
830 513117 : for (const auto & side_nodes : side_grouped_boundary_neighbors)
831 : {
832 371424 : std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
833 2953032 : for (const auto i : index_range(side_nodes_vec))
834 : {
835 2664792 : const Point vec_i = (*side_nodes_vec[i] - node);
836 11045328 : for (const auto j : make_range(i))
837 : {
838 8691672 : const Point vec_j = (*side_nodes_vec[j] - node);
839 8215416 : Point candidate_normal = vec_i.cross(vec_j);
840 8453544 : if (candidate_normal.norm() <= TOLERANCE)
841 328020 : continue;
842 :
843 8583300 : const PlaneConstraint candidate_plane{node, candidate_normal};
844 7896636 : valid_planes.emplace(candidate_plane);
845 : }
846 : }
847 : }
848 :
849 : // Fall back to point constraint
850 151869 : if (valid_planes.empty())
851 0 : return PointConstraint(node);
852 :
853 : // Combine all the planes together to get a common constraint
854 4278 : auto it = valid_planes.begin();
855 8556 : ConstraintVariant current = *it++;
856 1182526 : for (; it != valid_planes.end(); ++it)
857 : {
858 1059766 : current = intersect_constraints(current, *it);
859 :
860 : // This will catch cases where constraints have no intersection
861 : // (i.e., the element surface is non-planar)
862 : // Fall back to fixed node constraint
863 1030728 : if (std::holds_alternative<InvalidConstraint>(current))
864 : {
865 73 : current = PointConstraint(node);
866 71 : break;
867 : }
868 : }
869 :
870 4278 : return current;
871 : }
872 :
873 : // Applies the computed constraint (PointConstraint, LineConstraint, or
874 : // PlaneConstraint) to a node.
875 153999 : void VariationalSmootherConstraint::impose_constraint(const Node & node,
876 : const ConstraintVariant & constraint)
877 : {
878 :
879 4338 : libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
880 : "Cannot impose constraint using InvalidConstraint.");
881 :
882 153999 : if (std::holds_alternative<PointConstraint>(constraint))
883 29465 : fix_node(node);
884 124534 : else if (std::holds_alternative<LineConstraint>(constraint))
885 32092 : constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
886 92442 : else if (std::holds_alternative<PlaneConstraint>(constraint))
887 92442 : constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
888 :
889 : else
890 0 : libmesh_assert_msg(false, "Unknown constraint type.");
891 153999 : }
892 :
893 : } // namespace libMesh
|
