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 2715324 : auto get_positive_vector = [](const Point & vec) -> Point {
28 2715324 : 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 76488 : Point canonical{0, 0, 0};
34 : // Choose the canonical dimension to ensure the dot product below is nonzero
35 5384782 : for (const auto dim_id : make_range(3))
36 5536466 : if (!absolute_fuzzy_equals(vec(dim_id), 0.))
37 : {
38 2715324 : canonical(dim_id) = 1.;
39 2638836 : break;
40 : }
41 :
42 76488 : const auto dot_prod = vec * canonical;
43 76488 : libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));
44 :
45 4078812 : return (dot_prod > 0) ? vec.unit() : -vec.unit();
46 : };
47 :
48 5133146 : PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)
49 : {
50 5133146 : }
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 426 : bool PointConstraint::operator==(const PointConstraint & other) const
61 : {
62 426 : return _point.absolute_fuzzy_equals(other.point(), _tol);
63 : }
64 :
65 21309 : ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const
66 : {
67 : // using visit to resolve the variant to its actual type
68 : return std::visit(
69 61525 : [&](auto && o) -> ConstraintVariant {
70 21309 : if (!o.contains_point(*this))
71 : // Point is not on the constraint
72 0 : return InvalidConstraint();
73 :
74 599 : return *this;
75 : },
76 21309 : other);
77 : }
78 :
79 37914 : LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
80 37914 : : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
81 : {
82 37914 : libmesh_error_msg_if(_direction.norm() < _tol,
83 : "Can't define a line with zero magnitude direction vector.");
84 37914 : }
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 994 : 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 994 : 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 8804 : bool LineConstraint::is_parallel(const PlaneConstraint & p) const
117 : {
118 8804 : return _direction * p.normal() < _tol;
119 : }
120 :
121 10153 : ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const
122 : {
123 : // using visit to resolve the variant to its actual type
124 : return std::visit(
125 29315 : [&](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 10153 : 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 20020 : 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 10153 : other);
177 : }
178 :
179 2677410 : PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
180 2677410 : : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
181 : {
182 2677410 : libmesh_error_msg_if(_normal.norm() < _tol,
183 : "Can't define a plane with zero magnitude direction vector.");
184 2677410 : }
185 :
186 6195380 : bool PlaneConstraint::operator<(const PlaneConstraint & other) const
187 : {
188 6195380 : if (*this == other)
189 143958 : return false;
190 :
191 1115650 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
192 1084954 : return _normal < other.normal();
193 0 : return (_normal * _point) < (other.normal() * other.point());
194 : }
195 :
196 6227614 : bool PlaneConstraint::operator==(const PlaneConstraint & other) const
197 : {
198 6403170 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
199 31582 : return false;
200 5110426 : return this->contains_point(other.point());
201 : }
202 :
203 32234 : bool PlaneConstraint::is_parallel(const PlaneConstraint & p) const
204 : {
205 32234 : return _normal.absolute_fuzzy_equals(p.normal(), _tol);
206 : }
207 :
208 8804 : bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }
209 :
210 5141320 : bool PlaneConstraint::contains_point(const PointConstraint & p) const
211 : {
212 : // distance between the point and the plane
213 144827 : const Real dist = (p.point() - _point) * _normal;
214 5141320 : return std::abs(dist) < _tol;
215 : }
216 :
217 11005 : bool PlaneConstraint::contains_line(const LineConstraint & l) const
218 : {
219 11005 : const bool base_on_plane = this->contains_point(PointConstraint(l.point()));
220 11005 : const bool dir_orthogonal = std::abs(_normal * l.direction()) < _tol;
221 11005 : return base_on_plane && dir_orthogonal;
222 : }
223 :
224 43239 : ConstraintVariant PlaneConstraint::intersect(const ConstraintVariant & other) const
225 : {
226 : // using visit to resolve the variant to its actual type
227 : return std::visit(
228 124845 : [&](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 43248 : if (*this == o)
236 0 : return *this;
237 :
238 32234 : 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 31326 : const Point dir = this->_normal.cross(o.normal()); // direction of line of intersection
255 908 : libmesh_assert(dir.norm() > _tol);
256 908 : const Point w = _point - o.point();
257 :
258 : // Dot product terms used in 2x2 system
259 908 : const Real n1_dot_n1 = _normal * _normal;
260 908 : const Real n1_dot_n2 = _normal * o.normal();
261 908 : const Real n2_dot_n2 = o.normal() * o.normal();
262 908 : const Real n1_dot_w = _normal * w;
263 908 : const Real n2_dot_w = o.normal() * w;
264 :
265 32234 : const Real denom = -(n1_dot_n1 * n2_dot_n2 - n1_dot_n2 * n1_dot_n2);
266 908 : libmesh_assert(std::abs(denom) > _tol);
267 :
268 32234 : const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
269 908 : const Point p0 = _point + s * _normal;
270 :
271 32234 : return LineConstraint{p0, dir};
272 : }
273 :
274 : else if constexpr (std::is_same_v<T, LineConstraint>)
275 : {
276 11005 : if (this->contains_line(o))
277 62 : return o;
278 :
279 8804 : if (this->is_parallel(o))
280 : // Line is parallel and does not intersect the plane
281 0 : 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 248 : const Real denom = _normal * o.direction();
291 248 : libmesh_assert(std::abs(denom) > _tol);
292 8804 : const Real t = (_normal * (_point - o.point())) / denom;
293 8804 : 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 43239 : other);
309 : }
310 :
311 1349 : VariationalSmootherConstraint::VariationalSmootherConstraint(
312 1349 : System & sys, const bool & preserve_subdomain_boundaries)
313 1349 : : Constraint(), _sys(sys), _preserve_subdomain_boundaries(preserve_subdomain_boundaries)
314 : {
315 1349 : }
316 :
317 2546 : VariationalSmootherConstraint::~VariationalSmootherConstraint() = default;
318 :
319 1349 : void VariationalSmootherConstraint::constrain()
320 : {
321 1349 : const auto & mesh = _sys.get_mesh();
322 1349 : const auto dim = mesh.mesh_dimension();
323 :
324 : // Only compute the node to elem map once
325 76 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
326 1349 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
327 :
328 38 : const auto & boundary_info = mesh.get_boundary_info();
329 :
330 1387 : const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
331 :
332 : // Identify/constrain subdomain boundary nodes, if requested
333 76 : std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
334 1349 : if (_preserve_subdomain_boundaries)
335 : {
336 149440 : for (const auto * elem : mesh.active_element_ptr_range())
337 : {
338 2864 : const auto & sub_id1 = elem->subdomain_id();
339 296496 : for (const auto side : elem->side_index_range())
340 : {
341 245660 : const auto * neighbor = elem->neighbor_ptr(side);
342 245660 : if (neighbor == nullptr)
343 30498 : continue;
344 :
345 12072 : const auto & sub_id2 = neighbor->subdomain_id();
346 214278 : if (sub_id1 == sub_id2)
347 203412 : continue;
348 :
349 : // elem and neighbor are in different subdomains, and share nodes
350 : // that need to be constrained
351 31096 : for (const auto local_node_id : elem->nodes_on_side(side))
352 : {
353 26718 : const auto & node = mesh.node_ref(elem->node_id(local_node_id));
354 : // Make sure we haven't already processed this node
355 25986 : if (subdomain_boundary_map.count(node.id()))
356 18531 : 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 7665 : mesh, node, sub_id1, nodes_to_elem_map);
362 :
363 : // Determine which constraint should be imposed
364 : const auto subdomain_constraint =
365 7455 : 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 7665 : if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
370 3053 : 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 14700 : subdomain_boundary_map[node.id()] = subdomain_constraint;
378 :
379 : } // for local_node_id
380 :
381 : } // for side
382 1273 : } // for elem
383 : }
384 :
385 : // Loop through boundary nodes and impose constraints
386 64539 : for (const auto & bid : boundary_node_ids)
387 : {
388 63190 : 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 64970 : 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 64970 : 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 63190 : if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
400 : {
401 4402 : const auto & subdomain_constraint = it->second;
402 : // Combine current boundary constraint with previously determined
403 : // subdomain_constraint
404 248 : 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 4402 : if (std::holds_alternative<InvalidConstraint>(combined_constraint))
409 0 : combined_constraint = PointConstraint(node);
410 :
411 4402 : this->impose_constraint(node, combined_constraint);
412 : }
413 :
414 : else
415 58788 : this->impose_constraint(node, boundary_constraint);
416 :
417 : } // end bid
418 1349 : }
419 :
420 11289 : void VariationalSmootherConstraint::fix_node(const Node & node)
421 : {
422 41748 : for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
423 : {
424 30459 : const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
425 1716 : DofConstraintRow constraint_row;
426 : // Leave the constraint row as all zeros so this dof is independent from other dofs
427 30459 : const auto constrained_value = node(d);
428 : // Simply constrain this dof to retain it's current value
429 30459 : _sys.get_dof_map().add_constraint_row(
430 : constrained_dof_index, constraint_row, constrained_value, true);
431 : }
432 11289 : }
433 :
434 26838 : void VariationalSmootherConstraint::constrain_node_to_plane(const Node & node,
435 : const Point & ref_normal_vec)
436 : {
437 26838 : 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 1512 : std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
440 756 : 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 756 : 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 756 : libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
452 :
453 756 : Real max_abs_coef = 0.;
454 107352 : for (const auto d : make_range(dim))
455 : {
456 80514 : const auto coef = ref_normal_vec(d);
457 80514 : xyz_coefs.push_back(coef);
458 80514 : c -= coef * node(d);
459 :
460 2268 : const auto coef_abs = std::abs(coef);
461 80514 : if (coef_abs > max_abs_coef)
462 : {
463 756 : max_abs_coef = coef_abs;
464 756 : constrained_dim = d;
465 : }
466 : }
467 :
468 1512 : DofConstraintRow constraint_row;
469 107352 : for (const auto free_dim : make_range(dim))
470 : {
471 80514 : if (free_dim == constrained_dim)
472 26838 : continue;
473 :
474 53676 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
475 56700 : constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
476 : }
477 :
478 26838 : const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
479 26838 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
480 26838 : _sys.get_dof_map().add_constraint_row(
481 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
482 26838 : }
483 :
484 28116 : void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
485 : const Point & line_vec)
486 : {
487 28116 : 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 28908 : const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
492 28116 : auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
493 1584 : return std::abs(a) < std::abs(b);
494 1584 : });
495 28116 : const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
496 28116 : 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 792 : libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
509 107352 : for (const auto constrained_dim : make_range(dim))
510 : {
511 79236 : if (constrained_dim == free_dim)
512 28116 : continue;
513 :
514 2880 : DofConstraintRow constraint_row;
515 51120 : constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
516 : const auto inhomogeneous_part =
517 51120 : node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
518 51120 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
519 51120 : _sys.get_dof_map().add_constraint_row(
520 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
521 : }
522 28116 : }
523 :
524 70645 : 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 70645 : 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 70645 : if (!neighbors.size())
537 : {
538 : // Grab the element containing node
539 8946 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
540 : // Find the element side containing node
541 35429 : for (const auto &side : elem->side_index_range())
542 : {
543 35429 : const auto &nodes_on_side = elem->nodes_on_side(side);
544 : const auto it =
545 41587 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
546 112674 : return elem->node_id(local_node_id) == node.id();
547 1996 : });
548 :
549 36427 : if (it != nodes_on_side.end())
550 : {
551 88608 : for (const auto &local_node_id : nodes_on_side)
552 : // No need to add node itself as a neighbor
553 81906 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
554 2244 : *node_ptr != node)
555 70716 : neighbors.push_back(node_ptr);
556 252 : break;
557 : }
558 : }
559 : }
560 1990 : libmesh_assert(neighbors.size());
561 70645 : }
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 373744 : 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 10528 : bool nodes_share_bid = false;
596 :
597 : // Node ids local to containing_elem
598 363216 : const auto node_id = containing_elem.get_node_index(&boundary_node);
599 363216 : const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
600 :
601 1640952 : 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 2365720 : if (!(containing_elem.is_node_on_side(node_id, side_id) &&
605 794064 : containing_elem.is_node_on_side(neighbor_id, side_id)))
606 969576 : 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 16960 : std::vector<boundary_id_type> boundary_ids;
611 602080 : boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
612 602080 : if (boundary_ids.size())
613 : {
614 8576 : nodes_share_bid = true;
615 8576 : break;
616 : }
617 : }
618 373744 : return nodes_share_bid;
619 : }
620 :
621 : std::set<std::set<const Node *>>
622 7455 : 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 420 : std::vector<const Node *> neighbors;
633 7455 : 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 210 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
639 :
640 34364 : 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 26909 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
645 26909 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
646 758 : const Elem * common_elem = nullptr;
647 122404 : for (const auto * neigh_elem : elems_containing_neigh)
648 : {
649 117162 : if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
650 5380 : elems_containing_node.end())
651 : // We should be able to find a common element on the sub_id boundary
652 95495 : && (neigh_elem->subdomain_id() == sub_id))
653 706 : common_elem = neigh_elem;
654 : else
655 70432 : continue;
656 :
657 : // Now, determine whether node and neigh are on a side coincident
658 : // with the subdomain boundary
659 156058 : for (const auto common_side : common_elem->side_index_range())
660 : {
661 3690 : bool node_found_on_side = false;
662 3690 : bool neigh_found_on_side = false;
663 1146080 : for (const auto local_node_id : common_elem->nodes_on_side(common_side))
664 : {
665 1039885 : if (common_elem->node_id(local_node_id) == node.id())
666 1368 : node_found_on_side = true;
667 962831 : else if (common_elem->node_id(local_node_id) == neigh->id())
668 1880 : neigh_found_on_side = true;
669 : }
670 :
671 132007 : if (!(node_found_on_side && neigh_found_on_side &&
672 2024 : common_elem->neighbor_ptr(common_side)))
673 101016 : continue;
674 :
675 762 : 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 1524 : common_elem->neighbor_ptr(matched_side)->subdomain_id();
683 762 : const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
684 :
685 27051 : if (is_matched_side_on_subdomain_boundary)
686 : {
687 : // Store all nodes that live on this side
688 22995 : const auto nodes_on_side = common_elem->nodes_on_side(common_side);
689 1260 : std::set<const Node *> node_ptrs_on_side;
690 196528 : for (const auto local_node_id : nodes_on_side)
691 179069 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
692 22995 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
693 21735 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
694 :
695 630 : continue;
696 : }
697 :
698 : } // for common_side
699 :
700 : } // for neigh_elem
701 : }
702 :
703 7665 : return side_grouped_boundary_neighbors;
704 : }
705 :
706 : std::set<std::set<const Node *>>
707 63190 : 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 3560 : std::vector<const Node *> neighbors;
719 63190 : 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 1780 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
725 :
726 313962 : for (const auto * neigh : neighbors)
727 : {
728 : const bool is_neighbor_boundary_node =
729 250772 : boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
730 250772 : if (!is_neighbor_boundary_node)
731 45540 : 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 203912 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
736 203912 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
737 5744 : const Elem * common_elem = nullptr;
738 1018921 : for (const auto * neigh_elem : elems_containing_neigh)
739 : {
740 : const bool is_neigh_common =
741 837967 : std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
742 45916 : elems_containing_node.end();
743 815009 : if (!is_neigh_common)
744 428835 : continue;
745 21056 : 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 373744 : const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
750 : node, *neigh, *common_elem, boundary_info);
751 373744 : if (nodes_have_common_bid)
752 : {
753 : // Find the side coinciding with the shared boundary
754 1959600 : 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 1701776 : if (common_elem->neighbor_ptr(side))
759 1323440 : continue;
760 :
761 15864 : bool node_found_on_side = false;
762 15864 : bool neigh_found_on_side = false;
763 563172 : const auto nodes_on_side = common_elem->nodes_on_side(side);
764 4698780 : for (const auto local_node_id : nodes_on_side)
765 : {
766 4252104 : if (common_elem->node_id(local_node_id) == node.id())
767 10328 : node_found_on_side = true;
768 3768964 : else if (common_elem->node_id(local_node_id) == neigh->id())
769 11514 : neigh_found_on_side = true;
770 : }
771 563172 : if (!(node_found_on_side && neigh_found_on_side))
772 6520 : continue;
773 :
774 18688 : std::set<const Node *> node_ptrs_on_side;
775 2635236 : for (const auto local_node_id : nodes_on_side)
776 2368412 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
777 341056 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
778 322368 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
779 : }
780 295872 : continue;
781 287296 : }
782 : }
783 : }
784 :
785 64970 : return side_grouped_boundary_neighbors;
786 : }
787 :
788 70645 : 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 3980 : std::vector<const Node *> neighbors;
798 244808 : for (const auto & side : side_grouped_boundary_neighbors)
799 174163 : neighbors.insert(neighbors.end(), side.begin(), side.end());
800 :
801 : // Constrain the node to it's current location
802 70645 : if (dim == 1 || neighbors.size() == 1)
803 1125 : return PointConstraint{node};
804 :
805 69580 : if (dim == 2 || neighbors.size() == 2)
806 : {
807 : // Determine whether node + all neighbors are colinear
808 212 : bool all_colinear = true;
809 7738 : const Point ref_dir = (*neighbors[0] - node).unit();
810 14058 : for (const auto i : make_range(std::size_t(1), neighbors.size()))
811 : {
812 8614 : const Point delta = *(neighbors[i]) - node;
813 236 : libmesh_assert(delta.norm() >= TOLERANCE);
814 8378 : const Point dir = delta.unit();
815 8598 : if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
816 : {
817 52 : all_colinear = false;
818 1846 : break;
819 : }
820 : }
821 :
822 212 : if (all_colinear)
823 6000 : return LineConstraint{node, ref_dir};
824 :
825 1950 : return PointConstraint{node};
826 : }
827 :
828 : // dim == 3, neighbors.size() >= 3
829 3496 : std::set<PlaneConstraint> valid_planes;
830 221165 : for (const auto & side_nodes : side_grouped_boundary_neighbors)
831 : {
832 163593 : std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
833 1126131 : for (const auto i : index_range(side_nodes_vec))
834 : {
835 994260 : const Point vec_i = (*side_nodes_vec[i] - node);
836 3829314 : for (const auto j : make_range(i))
837 : {
838 2942922 : const Point vec_j = (*side_nodes_vec[j] - node);
839 2781666 : Point candidate_normal = vec_i.cross(vec_j);
840 2862294 : if (candidate_normal.norm() <= TOLERANCE)
841 184884 : continue;
842 :
843 2828250 : const PlaneConstraint candidate_plane{node, candidate_normal};
844 2601990 : valid_planes.emplace(candidate_plane);
845 : }
846 : }
847 : }
848 :
849 : // Fall back to point constraint
850 62054 : if (valid_planes.empty())
851 0 : return PointConstraint(node);
852 :
853 : // Combine all the planes together to get a common constraint
854 1748 : auto it = valid_planes.begin();
855 3496 : ConstraintVariant current = *it++;
856 122200 : for (; it != valid_planes.end(); ++it)
857 : {
858 61842 : 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 60146 : if (std::holds_alternative<InvalidConstraint>(current))
864 : {
865 0 : current = PointConstraint(node);
866 0 : break;
867 : }
868 : }
869 :
870 1748 : return current;
871 : }
872 :
873 : // Applies the computed constraint (PointConstraint, LineConstraint, or
874 : // PlaneConstraint) to a node.
875 66243 : void VariationalSmootherConstraint::impose_constraint(const Node & node,
876 : const ConstraintVariant & constraint)
877 : {
878 :
879 1866 : libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
880 : "Cannot impose constraint using InvalidConstraint.");
881 :
882 66243 : if (std::holds_alternative<PointConstraint>(constraint))
883 11289 : fix_node(node);
884 54954 : else if (std::holds_alternative<LineConstraint>(constraint))
885 28116 : constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
886 26838 : else if (std::holds_alternative<PlaneConstraint>(constraint))
887 26838 : constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
888 :
889 : else
890 0 : libmesh_assert_msg(false, "Unknown constraint type.");
891 66243 : }
892 :
893 : } // namespace libMesh
|
