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 793986 : auto get_positive_vector = [](const Point & vec) -> Point {
28 793986 : 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 65951 : Point canonical{0, 0, 0};
34 : // Choose the canonical dimension to ensure the dot product below is nonzero
35 1540244 : for (const auto dim_id : make_range(3))
36 1668192 : if (!absolute_fuzzy_equals(vec(dim_id), 0.))
37 : {
38 793986 : canonical(dim_id) = 1.;
39 728035 : break;
40 : }
41 :
42 65951 : const auto dot_prod = vec * canonical;
43 65951 : libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));
44 :
45 1193945 : return (dot_prod > 0) ? vec.unit() : -vec.unit();
46 : };
47 :
48 1468559 : PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)
49 : {
50 1468559 : }
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 192 : bool PointConstraint::operator==(const PointConstraint & other) const
61 : {
62 192 : return _point.absolute_fuzzy_equals(other.point(), _tol);
63 : }
64 :
65 13306 : ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const
66 : {
67 : // using visit to resolve the variant to its actual type
68 : return std::visit(
69 35542 : [&](auto && o) -> ConstraintVariant {
70 13306 : if (!o.contains_point(*this))
71 : // Point is not on the constraint
72 0 : return InvalidConstraint();
73 :
74 1093 : return *this;
75 : },
76 13306 : other);
77 : }
78 :
79 12408 : LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
80 12408 : : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
81 : {
82 12408 : libmesh_error_msg_if(_direction.norm() < _tol,
83 : "Can't define a line with zero magnitude direction vector.");
84 12408 : }
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 48 : bool LineConstraint::operator==(const LineConstraint & other) const
97 : {
98 52 : if (!(_direction.absolute_fuzzy_equals(other.direction(), _tol)))
99 4 : return false;
100 0 : return this->contains_point(other.point());
101 : }
102 :
103 192 : 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 192 : return _direction.cross(p.point() - _point).norm() < _tol;
109 : }
110 :
111 48 : bool LineConstraint::is_parallel(const LineConstraint & l) const
112 : {
113 48 : return _direction.absolute_fuzzy_equals(l.direction(), _tol);
114 : }
115 :
116 3240 : bool LineConstraint::is_parallel(const PlaneConstraint & p) const
117 : {
118 3240 : return _direction * p.normal() < _tol;
119 : }
120 :
121 4068 : ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const
122 : {
123 : // using visit to resolve the variant to its actual type
124 : return std::visit(
125 10848 : [&](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 4112 : if (*this == o)
132 0 : return *this;
133 :
134 48 : 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 4 : const Point delta = o.point() - _point;
146 44 : const Point cross_d1_d2 = _direction.cross(o.direction());
147 44 : const Real cross_dot = (delta.cross(o.direction())) * cross_d1_d2;
148 4 : const Real denom = cross_d1_d2.norm_sq();
149 :
150 48 : const Real t = cross_dot / denom;
151 4 : const Point intersection = _point + t * _direction;
152 :
153 : // Verify that intersection lies on both lines
154 48 : if (o.direction().cross(intersection - o.point()).norm() > _tol)
155 : // Lines do not intersect at a single point
156 0 : return InvalidConstraint();
157 :
158 48 : return PointConstraint{intersection};
159 : }
160 :
161 : else if constexpr (std::is_same_v<T, PlaneConstraint>)
162 7705 : 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 4068 : other);
177 : }
178 :
179 781578 : PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
180 781578 : : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
181 : {
182 781578 : libmesh_error_msg_if(_normal.norm() < _tol,
183 : "Can't define a plane with zero magnitude direction vector.");
184 781578 : }
185 :
186 1817672 : bool PlaneConstraint::operator<(const PlaneConstraint & other) const
187 : {
188 1817672 : if (*this == other)
189 121302 : return false;
190 :
191 386802 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
192 357261 : return _normal < other.normal();
193 0 : return (_normal * _point) < (other.normal() * other.point());
194 : }
195 :
196 1828376 : bool PlaneConstraint::operator==(const PlaneConstraint & other) const
197 : {
198 1980107 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
199 30261 : return false;
200 1460411 : return this->contains_point(other.point());
201 : }
202 :
203 10704 : bool PlaneConstraint::is_parallel(const PlaneConstraint & p) const
204 : {
205 10704 : return _normal.absolute_fuzzy_equals(p.normal(), _tol);
206 : }
207 :
208 3240 : bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }
209 :
210 1477497 : bool PlaneConstraint::contains_point(const PointConstraint & p) const
211 : {
212 : // distance between the point and the plane
213 122710 : const Real dist = (p.point() - _point) * _normal;
214 1477497 : return std::abs(dist) < _tol;
215 : }
216 :
217 4164 : bool PlaneConstraint::contains_line(const LineConstraint & l) const
218 : {
219 4164 : const bool base_on_plane = this->contains_point(PointConstraint(l.point()));
220 4164 : const bool dir_orthogonal = std::abs(_normal * l.direction()) < _tol;
221 4164 : return base_on_plane && dir_orthogonal;
222 : }
223 :
224 14868 : ConstraintVariant PlaneConstraint::intersect(const ConstraintVariant & other) const
225 : {
226 : // using visit to resolve the variant to its actual type
227 : return std::visit(
228 39648 : [&](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 16233 : if (*this == o)
236 0 : return *this;
237 :
238 10704 : 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 9812 : const Point dir = this->_normal.cross(o.normal()); // direction of line of intersection
255 892 : libmesh_assert(dir.norm() > _tol);
256 892 : const Point w = _point - o.point();
257 :
258 : // Dot product terms used in 2x2 system
259 892 : const Real n1_dot_n1 = _normal * _normal;
260 892 : const Real n1_dot_n2 = _normal * o.normal();
261 892 : const Real n2_dot_n2 = o.normal() * o.normal();
262 892 : const Real n1_dot_w = _normal * w;
263 892 : const Real n2_dot_w = o.normal() * w;
264 :
265 10704 : const Real denom = -(n1_dot_n1 * n2_dot_n2 - n1_dot_n2 * n1_dot_n2);
266 892 : libmesh_assert(std::abs(denom) > _tol);
267 :
268 10704 : const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
269 892 : const Point p0 = _point + s * _normal;
270 :
271 10704 : return LineConstraint{p0, dir};
272 : }
273 :
274 : else if constexpr (std::is_same_v<T, LineConstraint>)
275 : {
276 4164 : if (this->contains_line(o))
277 77 : return o;
278 :
279 3240 : if (this->is_parallel(o))
280 : // Line is parallel and does not intersect the plane
281 24 : 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 268 : const Real denom = _normal * o.direction();
291 268 : libmesh_assert(std::abs(denom) > _tol);
292 3216 : const Real t = (_normal * (_point - o.point())) / denom;
293 3216 : 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 14868 : other);
309 : }
310 :
311 2520 : std::ostream &operator<<(std::ostream &os, const ConstraintVariant & c)
312 : {
313 2520 : if (std::holds_alternative<PointConstraint>(c))
314 264 : os << "point constraint: point=" << std::get<PointConstraint>(c).point();
315 :
316 2376 : else if (std::holds_alternative<LineConstraint>(c))
317 : {
318 48 : const auto & line = std::get<LineConstraint>(c);
319 624 : os << "line constraint: point=" << line.point() << ", direction=" << line.direction();
320 : }
321 :
322 1800 : else if (std::holds_alternative<PlaneConstraint>(c))
323 : {
324 150 : const auto & plane = std::get<PlaneConstraint>(c);
325 1950 : os << "plane constraint: point=" << plane.point() << ", normal=" << plane.normal();
326 : }
327 :
328 0 : else if (std::holds_alternative<InvalidConstraint>(c))
329 0 : os << "invalid constraint";
330 :
331 2520 : return os;
332 : }
333 :
334 2414 : VariationalSmootherConstraint::VariationalSmootherConstraint(
335 2414 : System & sys, const bool & preserve_subdomain_boundaries, const unsigned int verbosity)
336 2414 : : Constraint(), _verbosity(verbosity), _sys(sys), _preserve_subdomain_boundaries(preserve_subdomain_boundaries) {
337 2414 : }
338 :
339 4556 : VariationalSmootherConstraint::~VariationalSmootherConstraint() = default;
340 :
341 2414 : void VariationalSmootherConstraint::constrain()
342 : {
343 2414 : const auto & mesh = _sys.get_mesh();
344 2414 : const auto dim = mesh.mesh_dimension();
345 136 : const auto & proc_id = mesh.processor_id();
346 :
347 : // Only compute the node to elem map once
348 136 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
349 2414 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
350 :
351 68 : const auto & boundary_info = mesh.get_boundary_info();
352 :
353 2482 : const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
354 :
355 : // Identify/constrain subdomain boundary nodes, if requested
356 136 : std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
357 2414 : if (_preserve_subdomain_boundaries)
358 : {
359 249127 : for (const auto * elem : mesh.active_element_ptr_range())
360 : {
361 10680 : const auto & sub_id1 = elem->subdomain_id();
362 473759 : for (const auto side : elem->side_index_range())
363 : {
364 385090 : const auto * neighbor = elem->neighbor_ptr(side);
365 385090 : if (neighbor == nullptr)
366 37818 : continue;
367 :
368 40928 : const auto & sub_id2 = neighbor->subdomain_id();
369 344852 : if (sub_id1 == sub_id2)
370 319246 : continue;
371 :
372 : // elem and neighbor are in different subdomains, and share nodes
373 : // that need to be constrained
374 36962 : for (const auto local_node_id : elem->nodes_on_side(side))
375 : {
376 32556 : const auto & node = mesh.node_ref(elem->node_id(local_node_id));
377 : // Make sure we haven't already processed this node.
378 : // We leave the responsibility of computing the constraint to
379 : // the proc that owns the node.
380 32152 : if (subdomain_boundary_map.count(node.id()) ||
381 1496 : node.processor_id() != proc_id)
382 29340 : continue;
383 :
384 : // Get the relevant nodal neighbors for the subdomain constraint
385 : const auto side_grouped_boundary_neighbors =
386 : get_neighbors_for_subdomain_constraint(
387 2236 : mesh, node, sub_id1, nodes_to_elem_map);
388 :
389 : // Determine which constraint should be imposed
390 : const auto subdomain_constraint =
391 2064 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
392 :
393 : // This subdomain boundary node does not lie on an external boundary,
394 : // go ahead and impose constraint
395 2236 : if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
396 : {
397 888 : if (_verbosity > 20)
398 0 : libMesh::out << "Imposing subdomain constraint on "
399 0 : << std::endl << " " << node << std::endl
400 0 : << " " << subdomain_constraint << std::endl;
401 :
402 888 : this->impose_constraint(node, subdomain_constraint);
403 : }
404 :
405 : // This subdomain boundary node could lie on an external boundary, save it
406 : // for later to combine with the external boundary constraint.
407 : // We also save constraints for non-boundary nodes so we don't try to
408 : // re-constrain the node when accessed from the neighboring elem.
409 : // See subdomain_boundary_map.count call above.
410 3956 : subdomain_boundary_map[node.id()] = subdomain_constraint;
411 :
412 : } // for local_node_id
413 :
414 : } // for side
415 2144 : } // for elem
416 : }
417 :
418 : // Loop through boundary nodes and impose constraints
419 99050 : for (const auto & bid : boundary_node_ids)
420 : {
421 96636 : const auto & node = mesh.node_ref(bid);
422 : // We leave the responsibility of computing the constraint to the proc
423 : // that owns the node.
424 96636 : if (node.processor_id() != proc_id)
425 72204 : continue;
426 :
427 : // Get the relevant nodal neighbors for the boundary constraint
428 : const auto side_grouped_boundary_neighbors = get_neighbors_for_boundary_constraint(
429 26468 : mesh, node, boundary_node_ids, boundary_info, nodes_to_elem_map);
430 :
431 : // Determine which constraint should be imposed
432 : const auto boundary_constraint =
433 26468 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
434 :
435 : // Check for the case where this boundary node is also part of a subdomain id boundary
436 24432 : if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
437 : {
438 1176 : const auto & subdomain_constraint = it->second;
439 : // Combine current boundary constraint with previously determined
440 : // subdomain_constraint
441 196 : auto combined_constraint = intersect_constraints(subdomain_constraint, boundary_constraint);
442 :
443 : // This will catch cases where constraints have no intersection
444 : // Fall back to fixed node constraint
445 1176 : if (std::holds_alternative<InvalidConstraint>(combined_constraint))
446 0 : combined_constraint = PointConstraint(node);
447 :
448 1176 : if (_verbosity > 20)
449 0 : libMesh::out << "Imposing boundary/subdomain constraint on "
450 0 : << std::endl << " " << node << std::endl
451 0 : << " " << combined_constraint << std::endl;
452 :
453 1176 : this->impose_constraint(node, combined_constraint);
454 : }
455 :
456 : else
457 : {
458 23256 : if (_verbosity > 20)
459 210 : libMesh::out << "Imposing boundary constraint on "
460 210 : << std::endl << node << std::endl
461 210 : << " " << boundary_constraint << std::endl;
462 :
463 23256 : this->impose_constraint(node, boundary_constraint);
464 : }
465 :
466 : } // end bid
467 2414 : }
468 :
469 3792 : void VariationalSmootherConstraint::fix_node(const Node & node)
470 : {
471 14340 : for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
472 : {
473 10548 : const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
474 1758 : DofConstraintRow constraint_row;
475 : // Leave the constraint row as all zeros so this dof is independent from other dofs
476 10548 : const auto constrained_value = node(d);
477 : // Simply constrain this dof to retain it's current value
478 10548 : _sys.get_dof_map().add_constraint_row(
479 : constrained_dof_index, constraint_row, constrained_value, true);
480 : }
481 3792 : }
482 :
483 12648 : void VariationalSmootherConstraint::constrain_node_to_plane(const Node & node,
484 : const Point & ref_normal_vec)
485 : {
486 12648 : const auto dim = _sys.get_mesh().mesh_dimension();
487 : // determine equation of plane: c_x * x + c_y * y + c_z * z + c = 0
488 2108 : std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
489 1054 : Real c = 0.;
490 :
491 : // We choose to constrain the dimension with the largest magnitude coefficient
492 : // This approach ensures the coefficients added to the constraint_row
493 : // (i.e., -c_xyz / c_max) have as small magnitude as possible
494 :
495 : // We initialize this to avoid maybe-uninitialized compiler error
496 1054 : unsigned int constrained_dim = 0;
497 :
498 : // Let's assert that we have a nonzero normal to ensure that constrained_dim
499 : // is always set
500 1054 : libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
501 :
502 1054 : Real max_abs_coef = 0.;
503 50592 : for (const auto d : make_range(dim))
504 : {
505 37944 : const auto coef = ref_normal_vec(d);
506 37944 : xyz_coefs.push_back(coef);
507 37944 : c -= coef * node(d);
508 :
509 3162 : const auto coef_abs = std::abs(coef);
510 37944 : if (coef_abs > max_abs_coef)
511 : {
512 1054 : max_abs_coef = coef_abs;
513 1054 : constrained_dim = d;
514 : }
515 : }
516 :
517 2108 : DofConstraintRow constraint_row;
518 50592 : for (const auto free_dim : make_range(dim))
519 : {
520 37944 : if (free_dim == constrained_dim)
521 12648 : continue;
522 :
523 25296 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
524 29512 : constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
525 : }
526 :
527 12648 : const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
528 12648 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
529 12648 : _sys.get_dof_map().add_constraint_row(
530 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
531 12648 : }
532 :
533 8880 : void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
534 : const Point & line_vec)
535 : {
536 8880 : const auto dim = _sys.get_mesh().mesh_dimension();
537 :
538 : // We will free the dimension most parallel to line_vec to keep the
539 : // constraint coefficients small
540 9620 : const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
541 8880 : auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
542 1480 : return std::abs(a) < std::abs(b);
543 1480 : });
544 8880 : const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
545 8880 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
546 :
547 : // A line is parameterized as r(t) = node + t * line_vec, so
548 : // x(t) = node(x) + t * line_vec(x)
549 : // y(t) = node(y) + t * line_vec(y)
550 : // z(t) = node(z) + t * line_vec(z)
551 : // Let's say we leave x free. Then t = (x(t) - node(x)) / line_vec(x)
552 : // Then y and z can be constrained as
553 : // y = node(y) + line_vec_y * (x(t) - node(x)) / line_vec(x)
554 : // = x(t) * line_vec(y) / line_vec(x) + (node(y) - node(x) * line_vec(y) / line_vec(x))
555 : // z = x(t) * line_vec(z) / line_vec(x) + (node(z) - node(x) * line_vec(z) / line_vec(x))
556 :
557 740 : libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
558 34008 : for (const auto constrained_dim : make_range(dim))
559 : {
560 25128 : if (constrained_dim == free_dim)
561 8880 : continue;
562 :
563 2708 : DofConstraintRow constraint_row;
564 16248 : constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
565 : const auto inhomogeneous_part =
566 16248 : node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
567 16248 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
568 16248 : _sys.get_dof_map().add_constraint_row(
569 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
570 : }
571 8880 : }
572 :
573 : // Utility function to determine whether two nodes share a boundary ID.
574 : // The motivation for this is that a sliding boundary node on a triangular
575 : // element can have a neighbor boundary node in the same element that is not
576 : // part of the same boundary
577 : // Consider the below example with nodes A, C, D, F, G that comprise elements
578 : // E1, E2, E3, with boundaries B1, B2, B3, B4. To determine the constraint
579 : // equations for the sliding node C, the neighboring nodes A and D need to
580 : // be identified to construct the line that C is allowed to slide along.
581 : // Note that neighbors A and D both share the boundary B1 with C.
582 : // Without ensuring that neighbors share the same boundary as the current
583 : // node, a neighboring node that does not lie on the same boundary
584 : // (i.e. F and G) might be selected to define the constraining line,
585 : // resulting in an incorrect constraint.
586 : // Note that if, for example, boundaries B1 and B2 were to be combined
587 : // into boundary B12, the node F would share a boundary id with node C
588 : // and result in an incorrect constraint. It would be useful to design
589 : // additional checks to detect cases like this.
590 :
591 : // B3
592 : // G-----------F
593 : // | \ / |
594 : // B4 | \ E2 / | B2
595 : // | \ / |
596 : // | E1 \ / E3 |
597 : // A-----C-----D
598 : // B1
599 :
600 174864 : bool VariationalSmootherConstraint::nodes_share_boundary_id(const Node & boundary_node,
601 : const Node & neighbor_node,
602 : const Elem & containing_elem,
603 : const BoundaryInfo & boundary_info)
604 : {
605 14572 : bool nodes_share_bid = false;
606 :
607 : // Node ids local to containing_elem
608 160292 : const auto node_id = containing_elem.get_node_index(&boundary_node);
609 160292 : const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
610 :
611 532680 : for (const auto side_id : containing_elem.side_index_range())
612 : {
613 : // We don't care about this side if it doesn't contain our boundary and neighbor nodes
614 808932 : if (!(containing_elem.is_node_on_side(node_id, side_id) &&
615 301188 : containing_elem.is_node_on_side(neighbor_id, side_id)))
616 265104 : continue;
617 :
618 : // If the current side, containing boundary_node and neighbor_node, lies on a boundary,
619 : // we can say that boundary_node and neighbor_node have a common boundary id.
620 20220 : std::vector<boundary_id_type> boundary_ids;
621 242640 : boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
622 242640 : if (boundary_ids.size())
623 : {
624 12494 : nodes_share_bid = true;
625 12494 : break;
626 : }
627 : }
628 174864 : return nodes_share_bid;
629 : }
630 :
631 : std::set<std::set<const Node *>>
632 2064 : VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
633 : const MeshBase & mesh,
634 : const Node & node,
635 : const subdomain_id_type sub_id,
636 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
637 : {
638 :
639 : // Find all the nodal neighbors... that is the nodes directly connected
640 : // to this node through one edge, or if none exists, use other nodes on the
641 : // containing face
642 344 : std::vector<const Node *> neighbors;
643 2064 : MeshTools::find_nodal_or_face_neighbors(
644 : mesh, node, nodes_to_elem_map, neighbors);
645 :
646 : // Each constituent set corresponds to neighbors sharing a face on the
647 : // subdomain boundary
648 172 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
649 :
650 10152 : for (const auto * neigh : neighbors)
651 : {
652 : // Determine whether the neighbor is on the subdomain boundary
653 : // First, find the common elements that both node and neigh belong to
654 8088 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
655 8088 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
656 674 : const Elem * common_elem = nullptr;
657 40551 : for (const auto * neigh_elem : elems_containing_neigh)
658 : {
659 37709 : if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
660 5484 : elems_containing_node.end())
661 : // We should be able to find a common element on the sub_id boundary
662 32463 : && (neigh_elem->subdomain_id() == sub_id))
663 710 : common_elem = neigh_elem;
664 : else
665 23765 : continue;
666 :
667 : // Now, determine whether node and neigh are on a side coincident
668 : // with the subdomain boundary
669 49980 : for (const auto common_side : common_elem->side_index_range())
670 : {
671 3362 : bool node_found_on_side = false;
672 3362 : bool neigh_found_on_side = false;
673 352928 : for (const auto local_node_id : common_elem->nodes_on_side(common_side))
674 : {
675 333342 : if (common_elem->node_id(local_node_id) == node.id())
676 1436 : node_found_on_side = true;
677 290552 : else if (common_elem->node_id(local_node_id) == neigh->id())
678 1718 : neigh_found_on_side = true;
679 : }
680 :
681 42404 : if (!(node_found_on_side && neigh_found_on_side &&
682 2244 : common_elem->neighbor_ptr(common_side)))
683 28160 : continue;
684 :
685 858 : const auto matched_side = common_side;
686 : // There could be multiple matched sides, so keep this next part
687 : // inside the common_side loop
688 : //
689 : // Does matched_side, containing both node and neigh, lie on the
690 : // sub_id subdomain boundary?
691 : const auto matched_neighbor_sub_id =
692 1716 : common_elem->neighbor_ptr(matched_side)->subdomain_id();
693 858 : const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
694 :
695 10618 : if (is_matched_side_on_subdomain_boundary)
696 : {
697 : // Store all nodes that live on this side
698 8046 : const auto nodes_on_side = common_elem->nodes_on_side(common_side);
699 1200 : std::set<const Node *> node_ptrs_on_side;
700 63352 : for (const auto local_node_id : nodes_on_side)
701 60402 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
702 8046 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
703 6846 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
704 :
705 600 : continue;
706 : }
707 :
708 : } // for common_side
709 :
710 : } // for neigh_elem
711 : }
712 :
713 2236 : return side_grouped_boundary_neighbors;
714 : }
715 :
716 : std::set<std::set<const Node *>>
717 24432 : VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
718 : const MeshBase & mesh,
719 : const Node & node,
720 : const std::unordered_set<dof_id_type> & boundary_node_ids,
721 : const BoundaryInfo & boundary_info,
722 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
723 : {
724 :
725 : // Find all the nodal neighbors... that is the nodes directly connected
726 : // to this node through one edge, or if none exists, use other nodes on the
727 : // containing face
728 4072 : std::vector<const Node *> neighbors;
729 24432 : MeshTools::find_nodal_or_face_neighbors(
730 : mesh, node, nodes_to_elem_map, neighbors);
731 :
732 : // Each constituent set corresponds to neighbors sharing a face on the
733 : // boundary
734 2036 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
735 :
736 137808 : for (const auto * neigh : neighbors)
737 : {
738 : const bool is_neighbor_boundary_node =
739 113376 : boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
740 113376 : if (!is_neighbor_boundary_node)
741 19822 : continue;
742 :
743 : // Determine whether nodes share a common boundary id
744 : // First, find the common element that both node and neigh belong to
745 91752 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
746 91752 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
747 7646 : const Elem * common_elem = nullptr;
748 558440 : for (const auto * neigh_elem : elems_containing_neigh)
749 : {
750 : const bool is_neigh_common =
751 506482 : std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
752 79588 : elems_containing_node.end();
753 466688 : if (!is_neigh_common)
754 266602 : continue;
755 29144 : common_elem = neigh_elem;
756 : // Keep this in the neigh_elem loop because there can be multiple common
757 : // elements Now, determine whether node and neigh share a common boundary
758 : // id
759 174864 : const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
760 : node, *neigh, *common_elem, boundary_info);
761 174864 : if (nodes_have_common_bid)
762 : {
763 : // Find the side coinciding with the shared boundary
764 878520 : for (const auto side : common_elem->side_index_range())
765 : {
766 : // We only care about external boundaries here, make sure side doesn't
767 : // have a neighbor
768 789308 : if (common_elem->neighbor_ptr(side))
769 567168 : continue;
770 :
771 20220 : bool node_found_on_side = false;
772 20220 : bool neigh_found_on_side = false;
773 242640 : const auto nodes_on_side = common_elem->nodes_on_side(side);
774 1747344 : for (const auto local_node_id : nodes_on_side)
775 : {
776 1630096 : if (common_elem->node_id(local_node_id) == node.id())
777 14560 : node_found_on_side = true;
778 1329984 : else if (common_elem->node_id(local_node_id) == neigh->id())
779 15244 : neigh_found_on_side = true;
780 : }
781 242640 : if (!(node_found_on_side && neigh_found_on_side))
782 6768 : continue;
783 :
784 26904 : std::set<const Node *> node_ptrs_on_side;
785 1084848 : for (const auto local_node_id : nodes_on_side)
786 1000376 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
787 174876 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
788 147972 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
789 : }
790 137434 : continue;
791 124940 : }
792 : }
793 : }
794 :
795 26468 : return side_grouped_boundary_neighbors;
796 : }
797 :
798 26496 : ConstraintVariant VariationalSmootherConstraint::determine_constraint(
799 : const Node & node,
800 : const unsigned int dim,
801 : const std::set<std::set<const Node *>> & side_grouped_boundary_neighbors)
802 : {
803 : // Determines the appropriate geometric constraint for a node based on its
804 : // neighbors.
805 :
806 : // Extract neighbors in flat vector
807 4416 : std::vector<const Node *> neighbors;
808 105518 : for (const auto & side : side_grouped_boundary_neighbors)
809 79022 : neighbors.insert(neighbors.end(), side.begin(), side.end());
810 :
811 : // Constrain the node to it's current location
812 26496 : if (dim == 1 || neighbors.size() == 1)
813 210 : return PointConstraint{node};
814 :
815 26316 : if (dim == 2 || neighbors.size() == 2)
816 : {
817 : // Determine whether node + all neighbors are colinear
818 185 : bool all_colinear = true;
819 2405 : const Point ref_dir = (*neighbors[0] - node).unit();
820 4516 : for (const auto i : make_range(std::size_t(1), neighbors.size()))
821 : {
822 3046 : const Point delta = *(neighbors[i]) - node;
823 234 : libmesh_assert(delta.norm() >= TOLERANCE);
824 2812 : const Point dir = delta.unit();
825 3009 : if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
826 : {
827 43 : all_colinear = false;
828 516 : break;
829 : }
830 : }
831 :
832 185 : if (all_colinear)
833 1988 : return LineConstraint{node, ref_dir};
834 :
835 602 : return PointConstraint{node};
836 : }
837 :
838 : // dim == 3, neighbors.size() >= 3
839 4016 : std::set<PlaneConstraint> valid_planes;
840 99040 : for (const auto & side_nodes : side_grouped_boundary_neighbors)
841 : {
842 81180 : std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
843 422916 : for (const auto i : index_range(side_nodes_vec))
844 : {
845 376900 : const Point vec_i = (*side_nodes_vec[i] - node);
846 1183902 : for (const auto j : make_range(i))
847 : {
848 905362 : const Point vec_j = (*side_nodes_vec[j] - node);
849 766498 : Point candidate_normal = vec_i.cross(vec_j);
850 835930 : if (candidate_normal.norm() <= TOLERANCE)
851 54352 : continue;
852 :
853 911412 : const PlaneConstraint candidate_plane{node, candidate_normal};
854 716661 : valid_planes.emplace(candidate_plane);
855 : }
856 : }
857 : }
858 :
859 : // Fall back to point constraint
860 24096 : if (valid_planes.empty())
861 0 : return PointConstraint(node);
862 :
863 : // Combine all the planes together to get a common constraint
864 2008 : auto it = valid_planes.begin();
865 4016 : ConstraintVariant current = *it++;
866 51118 : for (; it != valid_planes.end(); ++it)
867 : {
868 29286 : current = intersect_constraints(current, *it);
869 :
870 : // This will catch cases where constraints have no intersection
871 : // (i.e., the element surface is non-planar)
872 : // Fall back to fixed node constraint
873 27046 : if (std::holds_alternative<InvalidConstraint>(current))
874 : {
875 26 : current = PointConstraint(node);
876 24 : break;
877 : }
878 : }
879 :
880 2008 : return current;
881 : }
882 :
883 : // Applies the computed constraint (PointConstraint, LineConstraint, or
884 : // PlaneConstraint) to a node.
885 25320 : void VariationalSmootherConstraint::impose_constraint(const Node & node,
886 : const ConstraintVariant & constraint)
887 : {
888 :
889 2110 : libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
890 : "Cannot impose constraint using InvalidConstraint.");
891 :
892 25320 : if (std::holds_alternative<PointConstraint>(constraint))
893 3792 : fix_node(node);
894 21528 : else if (std::holds_alternative<LineConstraint>(constraint))
895 8880 : constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
896 12648 : else if (std::holds_alternative<PlaneConstraint>(constraint))
897 12648 : constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
898 :
899 : else
900 0 : libmesh_assert_msg(false, "Unknown constraint type.");
901 25320 : }
902 :
903 : } // namespace libMesh
|
