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 3801766 : auto get_positive_vector = [](const Point & vec) -> Point {
28 3801766 : 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 107092 : Point canonical{0, 0, 0};
34 : // Choose the canonical dimension to ensure the dot product below is nonzero
35 7397135 : for (const auto dim_id : make_range(3))
36 7605505 : if (!absolute_fuzzy_equals(vec(dim_id), 0.))
37 : {
38 3801766 : canonical(dim_id) = 1.;
39 3694674 : break;
40 : }
41 :
42 107092 : const auto dot_prod = vec * canonical;
43 107092 : libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));
44 :
45 5715138 : return (dot_prod > 0) ? vec.unit() : -vec.unit();
46 : };
47 :
48 7136960 : PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)
49 : {
50 7136960 : }
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 710 : bool PointConstraint::operator==(const PointConstraint & other) const
61 : {
62 710 : return _point.absolute_fuzzy_equals(other.point(), _tol);
63 : }
64 :
65 26464 : ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const
66 : {
67 : // using visit to resolve the variant to its actual type
68 : return std::visit(
69 76402 : [&](auto && o) -> ConstraintVariant {
70 26464 : if (!o.contains_point(*this))
71 : // Point is not on the constraint
72 0 : return InvalidConstraint();
73 :
74 749 : return *this;
75 : },
76 26464 : other);
77 : }
78 :
79 60847 : LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
80 60847 : : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
81 : {
82 60847 : libmesh_error_msg_if(_direction.norm() < _tol,
83 : "Can't define a line with zero magnitude direction vector.");
84 60847 : }
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 284 : bool LineConstraint::operator==(const LineConstraint & other) const
97 : {
98 292 : if (!(_direction.absolute_fuzzy_equals(other.direction(), _tol)))
99 8 : return false;
100 0 : return this->contains_point(other.point());
101 : }
102 :
103 1136 : 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 1136 : return _direction.cross(p.point() - _point).norm() < _tol;
109 : }
110 :
111 284 : bool LineConstraint::is_parallel(const LineConstraint & l) const
112 : {
113 284 : return _direction.absolute_fuzzy_equals(l.direction(), _tol);
114 : }
115 :
116 13419 : bool LineConstraint::is_parallel(const PlaneConstraint & p) const
117 : {
118 13419 : return _direction * p.normal() < _tol;
119 : }
120 :
121 15549 : ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const
122 : {
123 : // using visit to resolve the variant to its actual type
124 : return std::visit(
125 44895 : [&](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 15825 : if (*this == o)
132 0 : return *this;
133 :
134 284 : 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 8 : const Point delta = o.point() - _point;
146 276 : const Point cross_d1_d2 = _direction.cross(o.direction());
147 276 : const Real cross_dot = (delta.cross(o.direction())) * cross_d1_d2;
148 8 : const Real denom = cross_d1_d2.norm_sq();
149 :
150 284 : const Real t = cross_dot / denom;
151 8 : const Point intersection = _point + t * _direction;
152 :
153 : // Verify that intersection lies on both lines
154 284 : if (o.direction().cross(intersection - o.point()).norm() > _tol)
155 : // Lines do not intersect at a single point
156 0 : return InvalidConstraint();
157 :
158 284 : return PointConstraint{intersection};
159 : }
160 :
161 : else if constexpr (std::is_same_v<T, PlaneConstraint>)
162 30100 : 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 15549 : other);
177 : }
178 :
179 3740919 : PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
180 3740919 : : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
181 : {
182 3740919 : libmesh_error_msg_if(_normal.norm() < _tol,
183 : "Can't define a plane with zero magnitude direction vector.");
184 3740919 : }
185 :
186 8624326 : bool PlaneConstraint::operator<(const PlaneConstraint & other) const
187 : {
188 8624326 : if (*this == other)
189 200078 : return false;
190 :
191 1564129 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
192 1521304 : return _normal < other.normal();
193 0 : return (_normal * _point) < (other.normal() * other.point());
194 : }
195 :
196 8675091 : bool PlaneConstraint::operator==(const PlaneConstraint & other) const
197 : {
198 8919430 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
199 44370 : return false;
200 7103022 : return this->contains_point(other.point());
201 : }
202 :
203 50765 : bool PlaneConstraint::is_parallel(const PlaneConstraint & p) const
204 : {
205 50765 : return _normal.absolute_fuzzy_equals(p.normal(), _tol);
206 : }
207 :
208 13419 : bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }
209 :
210 7143757 : bool PlaneConstraint::contains_point(const PointConstraint & p) const
211 : {
212 : // distance between the point and the plane
213 201229 : const Real dist = (p.point() - _point) * _normal;
214 7143757 : return std::abs(dist) < _tol;
215 : }
216 :
217 16117 : bool PlaneConstraint::contains_line(const LineConstraint & l) const
218 : {
219 16117 : const bool base_on_plane = this->contains_point(PointConstraint(l.point()));
220 16117 : const bool dir_orthogonal = std::abs(_normal * l.direction()) < _tol;
221 16117 : return base_on_plane && dir_orthogonal;
222 : }
223 :
224 66882 : ConstraintVariant PlaneConstraint::intersect(const ConstraintVariant & other) const
225 : {
226 : // using visit to resolve the variant to its actual type
227 : return std::visit(
228 193110 : [&](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 67580 : if (*this == o)
236 0 : return *this;
237 :
238 50765 : 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 49335 : const Point dir = this->_normal.cross(o.normal()); // direction of line of intersection
255 1430 : libmesh_assert(dir.norm() > _tol);
256 1430 : const Point w = _point - o.point();
257 :
258 : // Dot product terms used in 2x2 system
259 1430 : const Real n1_dot_n1 = _normal * _normal;
260 1430 : const Real n1_dot_n2 = _normal * o.normal();
261 1430 : const Real n2_dot_n2 = o.normal() * o.normal();
262 1430 : const Real n1_dot_w = _normal * w;
263 1430 : const Real n2_dot_w = o.normal() * w;
264 :
265 50765 : const Real denom = -(n1_dot_n1 * n2_dot_n2 - n1_dot_n2 * n1_dot_n2);
266 1430 : libmesh_assert(std::abs(denom) > _tol);
267 :
268 50765 : const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
269 1430 : const Point p0 = _point + s * _normal;
270 :
271 50765 : return LineConstraint{p0, dir};
272 : }
273 :
274 : else if constexpr (std::is_same_v<T, LineConstraint>)
275 : {
276 16117 : if (this->contains_line(o))
277 76 : return o;
278 :
279 13419 : 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 378 : const Real denom = _normal * o.direction();
291 378 : libmesh_assert(std::abs(denom) > _tol);
292 13419 : const Real t = (_normal * (_point - o.point())) / denom;
293 13419 : 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 66882 : other);
309 : }
310 :
311 2130 : VariationalSmootherConstraint::VariationalSmootherConstraint(
312 2130 : System & sys, const bool & preserve_subdomain_boundaries)
313 2130 : : Constraint(), _sys(sys), _preserve_subdomain_boundaries(preserve_subdomain_boundaries)
314 : {
315 2130 : }
316 :
317 4020 : VariationalSmootherConstraint::~VariationalSmootherConstraint() = default;
318 :
319 2130 : void VariationalSmootherConstraint::constrain()
320 : {
321 2130 : const auto & mesh = _sys.get_mesh();
322 2130 : const auto dim = mesh.mesh_dimension();
323 :
324 : // Only compute the node to elem map once
325 120 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
326 2130 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
327 :
328 60 : const auto & boundary_info = mesh.get_boundary_info();
329 :
330 2190 : const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
331 :
332 : // Identify/constrain subdomain boundary nodes, if requested
333 120 : std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
334 2130 : if (_preserve_subdomain_boundaries)
335 : {
336 187190 : for (const auto * elem : mesh.active_element_ptr_range())
337 : {
338 3576 : const auto & sub_id1 = elem->subdomain_id();
339 371472 : for (const auto side : elem->side_index_range())
340 : {
341 307998 : const auto * neighbor = elem->neighbor_ptr(side);
342 307998 : if (neighbor == nullptr)
343 47058 : continue;
344 :
345 14624 : const auto & sub_id2 = neighbor->subdomain_id();
346 259576 : if (sub_id1 == sub_id2)
347 246744 : continue;
348 :
349 : // elem and neighbor are in different subdomains, and share nodes
350 : // that need to be constrained
351 34808 : for (const auto local_node_id : elem->nodes_on_side(side))
352 : {
353 29784 : const auto & node = mesh.node_ref(elem->node_id(local_node_id));
354 : // Make sure we haven't already processed this node
355 28968 : if (subdomain_boundary_map.count(node.id()))
356 20235 : 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 8979 : mesh, node, sub_id1, nodes_to_elem_map);
362 :
363 : // Determine which constraint should be imposed
364 : const auto subdomain_constraint =
365 8733 : 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 8979 : if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
370 3479 : 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 17220 : subdomain_boundary_map[node.id()] = subdomain_constraint;
378 :
379 : } // for local_node_id
380 :
381 : } // for side
382 1876 : } // for elem
383 : }
384 :
385 : // Loop through boundary nodes and impose constraints
386 108630 : for (const auto & bid : boundary_node_ids)
387 : {
388 106500 : 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 109500 : 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 109500 : 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 106500 : if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
400 : {
401 5254 : const auto & subdomain_constraint = it->second;
402 : // Combine current boundary constraint with previously determined
403 : // subdomain_constraint
404 296 : 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 5254 : if (std::holds_alternative<InvalidConstraint>(combined_constraint))
409 0 : combined_constraint = PointConstraint(node);
410 :
411 5254 : this->impose_constraint(node, combined_constraint);
412 : }
413 :
414 : else
415 101246 : this->impose_constraint(node, boundary_constraint);
416 :
417 : } // end bid
418 2130 : }
419 :
420 17111 : void VariationalSmootherConstraint::fix_node(const Node & node)
421 : {
422 63545 : for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
423 : {
424 46434 : const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
425 2616 : DofConstraintRow constraint_row;
426 : // Leave the constraint row as all zeros so this dof is independent from other dofs
427 46434 : const auto constrained_value = node(d);
428 : // Simply constrain this dof to retain it's current value
429 46434 : _sys.get_dof_map().add_constraint_row(
430 : constrained_dof_index, constraint_row, constrained_value, true);
431 : }
432 17111 : }
433 :
434 47144 : void VariationalSmootherConstraint::constrain_node_to_plane(const Node & node,
435 : const Point & ref_normal_vec)
436 : {
437 47144 : 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 2656 : std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
440 1328 : 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 1328 : 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 1328 : libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
452 :
453 1328 : Real max_abs_coef = 0.;
454 188576 : for (const auto d : make_range(dim))
455 : {
456 141432 : const auto coef = ref_normal_vec(d);
457 141432 : xyz_coefs.push_back(coef);
458 141432 : c -= coef * node(d);
459 :
460 3984 : const auto coef_abs = std::abs(coef);
461 141432 : if (coef_abs > max_abs_coef)
462 : {
463 1328 : max_abs_coef = coef_abs;
464 1328 : constrained_dim = d;
465 : }
466 : }
467 :
468 2656 : DofConstraintRow constraint_row;
469 188576 : for (const auto free_dim : make_range(dim))
470 : {
471 141432 : if (free_dim == constrained_dim)
472 47144 : continue;
473 :
474 94288 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
475 99600 : constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
476 : }
477 :
478 47144 : const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
479 47144 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
480 47144 : _sys.get_dof_map().add_constraint_row(
481 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
482 47144 : }
483 :
484 45724 : void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
485 : const Point & line_vec)
486 : {
487 45724 : 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 47012 : const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
492 45724 : auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
493 2576 : return std::abs(a) < std::abs(b);
494 2576 : });
495 45724 : const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
496 45724 : 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 1288 : libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
509 173950 : for (const auto constrained_dim : make_range(dim))
510 : {
511 128226 : if (constrained_dim == free_dim)
512 45724 : continue;
513 :
514 4648 : DofConstraintRow constraint_row;
515 82502 : constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
516 : const auto inhomogeneous_part =
517 82502 : node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
518 82502 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
519 82502 : _sys.get_dof_map().add_constraint_row(
520 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
521 : }
522 45724 : }
523 :
524 115233 : 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 115233 : 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 115233 : if (!neighbors.size())
537 : {
538 : // Grab the element containing node
539 11857 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
540 : // Find the element side containing node
541 43665 : for (const auto &side : elem->side_index_range())
542 : {
543 43665 : const auto &nodes_on_side = elem->nodes_on_side(side);
544 : const auto it =
545 51227 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
546 139935 : return elem->node_id(local_node_id) == node.id();
547 2460 : });
548 :
549 44895 : if (it != nodes_on_side.end())
550 : {
551 116014 : for (const auto &local_node_id : nodes_on_side)
552 : // No need to add node itself as a neighbor
553 107091 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
554 2934 : *node_ptr != node)
555 92300 : neighbors.push_back(node_ptr);
556 334 : break;
557 : }
558 : }
559 : }
560 3246 : libmesh_assert(neighbors.size());
561 115233 : }
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 673364 : 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 18968 : bool nodes_share_bid = false;
596 :
597 : // Node ids local to containing_elem
598 654396 : const auto node_id = containing_elem.get_node_index(&boundary_node);
599 654396 : const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
600 :
601 2572330 : 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 3797861 : if (!(containing_elem.is_node_on_side(node_id, side_id) &&
605 1318541 : containing_elem.is_node_on_side(neighbor_id, side_id)))
606 1459476 : 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 28728 : std::vector<boundary_id_type> boundary_ids;
611 1019844 : boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
612 1019844 : if (boundary_ids.size())
613 : {
614 16348 : nodes_share_bid = true;
615 16348 : break;
616 : }
617 : }
618 673364 : return nodes_share_bid;
619 : }
620 :
621 : std::set<std::set<const Node *>>
622 8733 : 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 492 : std::vector<const Node *> neighbors;
633 8733 : 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 246 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
639 :
640 41961 : 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 33228 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
645 33228 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
646 936 : const Elem * common_elem = nullptr;
647 137598 : for (const auto * neigh_elem : elems_containing_neigh)
648 : {
649 131514 : if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
650 5880 : elems_containing_node.end())
651 : // We should be able to find a common element on the sub_id boundary
652 104370 : && (neigh_elem->subdomain_id() == sub_id))
653 872 : common_elem = neigh_elem;
654 : else
655 73414 : continue;
656 :
657 : // Now, determine whether node and neigh are on a side coincident
658 : // with the subdomain boundary
659 191842 : for (const auto common_side : common_elem->side_index_range())
660 : {
661 4532 : bool node_found_on_side = false;
662 4532 : bool neigh_found_on_side = false;
663 1399824 : for (const auto local_node_id : common_elem->nodes_on_side(common_side))
664 : {
665 1269178 : if (common_elem->node_id(local_node_id) == node.id())
666 1692 : node_found_on_side = true;
667 1174340 : else if (common_elem->node_id(local_node_id) == neigh->id())
668 2092 : neigh_found_on_side = true;
669 : }
670 :
671 162178 : if (!(node_found_on_side && neigh_found_on_side &&
672 2584 : common_elem->neighbor_ptr(common_side)))
673 123924 : continue;
674 :
675 940 : 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 1880 : common_elem->neighbor_ptr(matched_side)->subdomain_id();
683 940 : const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
684 :
685 33370 : if (is_matched_side_on_subdomain_boundary)
686 : {
687 : // Store all nodes that live on this side
688 28908 : const auto nodes_on_side = common_elem->nodes_on_side(common_side);
689 1584 : std::set<const Node *> node_ptrs_on_side;
690 245944 : for (const auto local_node_id : nodes_on_side)
691 223964 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
692 28908 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
693 27324 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
694 :
695 792 : continue;
696 : }
697 :
698 : } // for common_side
699 :
700 : } // for neigh_elem
701 : }
702 :
703 8979 : return side_grouped_boundary_neighbors;
704 : }
705 :
706 : std::set<std::set<const Node *>>
707 106500 : 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 6000 : std::vector<const Node *> neighbors;
719 106500 : 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 3000 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
725 :
726 553516 : for (const auto * neigh : neighbors)
727 : {
728 : const bool is_neighbor_boundary_node =
729 447016 : boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
730 447016 : if (!is_neighbor_boundary_node)
731 65826 : 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 379282 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
736 379282 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
737 10684 : const Elem * common_elem = nullptr;
738 1684688 : for (const auto * neigh_elem : elems_containing_neigh)
739 : {
740 : const bool is_neigh_common =
741 1342178 : std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
742 73544 : elems_containing_node.end();
743 1305406 : if (!is_neigh_common)
744 614238 : continue;
745 37936 : 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 673364 : const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
750 : node, *neigh, *common_elem, boundary_info);
751 673364 : if (nodes_have_common_bid)
752 : {
753 : // Find the side coinciding with the shared boundary
754 3664310 : 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 3170828 : if (common_elem->neighbor_ptr(side))
759 2435584 : continue;
760 :
761 31800 : bool node_found_on_side = false;
762 31800 : bool neigh_found_on_side = false;
763 1128900 : const auto nodes_on_side = common_elem->nodes_on_side(side);
764 8620820 : for (const auto local_node_id : nodes_on_side)
765 : {
766 7702960 : if (common_elem->node_id(local_node_id) == node.id())
767 20480 : node_found_on_side = true;
768 6764880 : else if (common_elem->node_id(local_node_id) == neigh->id())
769 21848 : neigh_found_on_side = true;
770 : }
771 1128900 : if (!(node_found_on_side && neigh_found_on_side))
772 13536 : continue;
773 :
774 36528 : std::set<const Node *> node_ptrs_on_side;
775 4701052 : for (const auto local_node_id : nodes_on_side)
776 4166840 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
777 666636 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
778 630108 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
779 : }
780 564006 : continue;
781 547658 : }
782 : }
783 : }
784 :
785 109500 : return side_grouped_boundary_neighbors;
786 : }
787 :
788 115233 : 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 6492 : std::vector<const Node *> neighbors;
798 429905 : for (const auto & side : side_grouped_boundary_neighbors)
799 314672 : neighbors.insert(neighbors.end(), side.begin(), side.end());
800 :
801 : // Constrain the node to it's current location
802 115233 : if (dim == 1 || neighbors.size() == 1)
803 1125 : return PointConstraint{node};
804 :
805 114168 : if (dim == 2 || neighbors.size() == 2)
806 : {
807 : // Determine whether node + all neighbors are colinear
808 370 : bool all_colinear = true;
809 13505 : const Point ref_dir = (*neighbors[0] - node).unit();
810 26696 : for (const auto i : make_range(std::size_t(1), neighbors.size()))
811 : {
812 17082 : const Point delta = *(neighbors[i]) - node;
813 468 : libmesh_assert(delta.norm() >= TOLERANCE);
814 16614 : const Point dir = delta.unit();
815 17008 : if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
816 : {
817 86 : all_colinear = false;
818 3053 : break;
819 : }
820 : }
821 :
822 370 : if (all_colinear)
823 10650 : return LineConstraint{node, ref_dir};
824 :
825 3225 : return PointConstraint{node};
826 : }
827 :
828 : // dim == 3, neighbors.size() >= 3
829 5692 : std::set<PlaneConstraint> valid_planes;
830 391707 : for (const auto & side_nodes : side_grouped_boundary_neighbors)
831 : {
832 298862 : std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
833 1800134 : for (const auto i : index_range(side_nodes_vec))
834 : {
835 1551980 : const Point vec_i = (*side_nodes_vec[i] - node);
836 5507470 : for (const auto j : make_range(i))
837 : {
838 4110630 : const Point vec_j = (*side_nodes_vec[j] - node);
839 3885390 : Point candidate_normal = vec_i.cross(vec_j);
840 3998010 : if (candidate_normal.norm() <= TOLERANCE)
841 257091 : continue;
842 :
843 3951675 : const PlaneConstraint candidate_plane{node, candidate_normal};
844 3635541 : valid_planes.emplace(candidate_plane);
845 : }
846 : }
847 : }
848 :
849 : // Fall back to point constraint
850 101033 : if (valid_planes.empty())
851 0 : return PointConstraint(node);
852 :
853 : // Combine all the planes together to get a common constraint
854 2846 : auto it = valid_planes.begin();
855 5692 : ConstraintVariant current = *it++;
856 189409 : for (; it != valid_planes.end(); ++it)
857 : {
858 90866 : 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 88376 : if (std::holds_alternative<InvalidConstraint>(current))
864 : {
865 0 : current = PointConstraint(node);
866 0 : break;
867 : }
868 : }
869 :
870 2846 : return current;
871 : }
872 :
873 : // Applies the computed constraint (PointConstraint, LineConstraint, or
874 : // PlaneConstraint) to a node.
875 109979 : void VariationalSmootherConstraint::impose_constraint(const Node & node,
876 : const ConstraintVariant & constraint)
877 : {
878 :
879 3098 : libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
880 : "Cannot impose constraint using InvalidConstraint.");
881 :
882 109979 : if (std::holds_alternative<PointConstraint>(constraint))
883 17111 : fix_node(node);
884 92868 : else if (std::holds_alternative<LineConstraint>(constraint))
885 45724 : constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
886 47144 : else if (std::holds_alternative<PlaneConstraint>(constraint))
887 47144 : constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
888 :
889 : else
890 0 : libmesh_assert_msg(false, "Unknown constraint type.");
891 109979 : }
892 :
893 : } // namespace libMesh
|
