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 4682521 : auto get_positive_vector = [](const Point & vec) -> Point {
28 4682521 : 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 131902 : Point canonical{0, 0, 0};
34 : // Choose the canonical dimension to ensure the dot product below is nonzero
35 9084308 : for (const auto dim_id : make_range(3))
36 9340204 : if (!absolute_fuzzy_equals(vec(dim_id), 0.))
37 : {
38 4682521 : canonical(dim_id) = 1.;
39 4550619 : break;
40 : }
41 :
42 131902 : const auto dot_prod = vec * canonical;
43 131902 : libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));
44 :
45 7036301 : return (dot_prod > 0) ? vec.unit() : -vec.unit();
46 : };
47 :
48 8660483 : PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)
49 : {
50 8660483 : }
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 1136 : bool PointConstraint::operator==(const PointConstraint & other) const
61 : {
62 1136 : return _point.absolute_fuzzy_equals(other.point(), _tol);
63 : }
64 :
65 77692 : ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const
66 : {
67 : // using visit to resolve the variant to its actual type
68 : return std::visit(
69 224324 : [&](auto && o) -> ConstraintVariant {
70 77692 : if (!o.contains_point(*this))
71 : // Point is not on the constraint
72 0 : return InvalidConstraint();
73 :
74 2190 : return *this;
75 : },
76 77692 : other);
77 : }
78 :
79 73414 : LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
80 73414 : : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
81 : {
82 73414 : libmesh_error_msg_if(_direction.norm() < _tol,
83 : "Can't define a line with zero magnitude direction vector.");
84 73414 : }
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 19170 : bool LineConstraint::is_parallel(const PlaneConstraint & p) const
117 : {
118 19170 : return _direction * p.normal() < _tol;
119 : }
120 :
121 24069 : ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const
122 : {
123 : // using visit to resolve the variant to its actual type
124 : return std::visit(
125 69495 : [&](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 24345 : 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 46900 : 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 24069 : other);
177 : }
178 :
179 4609107 : PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
180 4609107 : : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
181 : {
182 4609107 : libmesh_error_msg_if(_normal.norm() < _tol,
183 : "Can't define a plane with zero magnitude direction vector.");
184 4609107 : }
185 :
186 10699106 : bool PlaneConstraint::operator<(const PlaneConstraint & other) const
187 : {
188 10699106 : if (*this == other)
189 242596 : return false;
190 :
191 2145797 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
192 2086832 : return _normal < other.normal();
193 0 : return (_normal * _point) < (other.normal() * other.point());
194 : }
195 :
196 10762438 : bool PlaneConstraint::operator==(const PlaneConstraint & other) const
197 : {
198 11065791 : if (!(_normal.absolute_fuzzy_equals(other.normal(), _tol)))
199 60013 : return false;
200 8612274 : return this->contains_point(other.point());
201 : }
202 :
203 63332 : bool PlaneConstraint::is_parallel(const PlaneConstraint & p) const
204 : {
205 63332 : return _normal.absolute_fuzzy_equals(p.normal(), _tol);
206 : }
207 :
208 19170 : bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }
209 :
210 8712331 : bool PlaneConstraint::contains_point(const PointConstraint & p) const
211 : {
212 : // distance between the point and the plane
213 245416 : const Real dist = (p.point() - _point) * _normal;
214 8712331 : return std::abs(dist) < _tol;
215 : }
216 :
217 24637 : bool PlaneConstraint::contains_line(const LineConstraint & l) const
218 : {
219 24637 : const bool base_on_plane = this->contains_point(PointConstraint(l.point()));
220 24637 : const bool dir_orthogonal = std::abs(_normal * l.direction()) < _tol;
221 24637 : return base_on_plane && dir_orthogonal;
222 : }
223 :
224 87969 : ConstraintVariant PlaneConstraint::intersect(const ConstraintVariant & other) const
225 : {
226 : // using visit to resolve the variant to its actual type
227 : return std::visit(
228 253995 : [&](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 86986 : if (*this == o)
236 0 : return *this;
237 :
238 63332 : 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 61548 : const Point dir = this->_normal.cross(o.normal()); // direction of line of intersection
255 1784 : libmesh_assert(dir.norm() > _tol);
256 1784 : const Point w = _point - o.point();
257 :
258 : // Dot product terms used in 2x2 system
259 1784 : const Real n1_dot_n1 = _normal * _normal;
260 1784 : const Real n1_dot_n2 = _normal * o.normal();
261 1784 : const Real n2_dot_n2 = o.normal() * o.normal();
262 1784 : const Real n1_dot_w = _normal * w;
263 1784 : const Real n2_dot_w = o.normal() * w;
264 :
265 63332 : const Real denom = -(n1_dot_n1 * n2_dot_n2 - n1_dot_n2 * n1_dot_n2);
266 1784 : libmesh_assert(std::abs(denom) > _tol);
267 :
268 63332 : const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
269 1784 : const Point p0 = _point + s * _normal;
270 :
271 63332 : return LineConstraint{p0, dir};
272 : }
273 :
274 : else if constexpr (std::is_same_v<T, LineConstraint>)
275 : {
276 24637 : if (this->contains_line(o))
277 154 : return o;
278 :
279 19170 : if (this->is_parallel(o))
280 : // Line is parallel and does not intersect the plane
281 142 : 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 536 : const Real denom = _normal * o.direction();
291 536 : libmesh_assert(std::abs(denom) > _tol);
292 19028 : const Real t = (_normal * (_point - o.point())) / denom;
293 19028 : 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 87969 : other);
309 : }
310 :
311 0 : std::ostream &operator<<(std::ostream &os, const ConstraintVariant & c)
312 : {
313 0 : if (std::holds_alternative<PointConstraint>(c))
314 0 : os << "point constraint: point=" << std::get<PointConstraint>(c).point();
315 :
316 0 : else if (std::holds_alternative<LineConstraint>(c))
317 : {
318 0 : const auto & line = std::get<LineConstraint>(c);
319 0 : os << "line constraint: point=" << line.point() << ", direction=" << line.direction();
320 : }
321 :
322 0 : else if (std::holds_alternative<PlaneConstraint>(c))
323 : {
324 0 : const auto & plane = std::get<PlaneConstraint>(c);
325 0 : 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 0 : 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 :
346 : // Only compute the node to elem map once
347 136 : std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
348 2414 : MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
349 :
350 68 : const auto & boundary_info = mesh.get_boundary_info();
351 :
352 2482 : const auto boundary_node_ids = MeshTools::find_boundary_nodes(mesh);
353 :
354 : // Identify/constrain subdomain boundary nodes, if requested
355 136 : std::unordered_map<dof_id_type, ConstraintVariant> subdomain_boundary_map;
356 2414 : if (_preserve_subdomain_boundaries)
357 : {
358 551830 : for (const auto * elem : mesh.active_element_ptr_range())
359 : {
360 10680 : const auto & sub_id1 = elem->subdomain_id();
361 1001952 : for (const auto side : elem->side_index_range())
362 : {
363 812382 : const auto * neighbor = elem->neighbor_ptr(side);
364 812382 : if (neighbor == nullptr)
365 83490 : continue;
366 :
367 40928 : const auto & sub_id2 = neighbor->subdomain_id();
368 726472 : if (sub_id1 == sub_id2)
369 698832 : continue;
370 :
371 : // elem and neighbor are in different subdomains, and share nodes
372 : // that need to be constrained
373 48488 : for (const auto local_node_id : elem->nodes_on_side(side))
374 : {
375 42048 : const auto & node = mesh.node_ref(elem->node_id(local_node_id));
376 : // Make sure we haven't already processed this node
377 40896 : if (subdomain_boundary_map.count(node.id()))
378 28684 : continue;
379 :
380 : // Get the relevant nodal neighbors for the subdomain constraint
381 : const auto side_grouped_boundary_neighbors =
382 : get_neighbors_for_subdomain_constraint(
383 12556 : mesh, node, sub_id1, nodes_to_elem_map);
384 :
385 : // Determine which constraint should be imposed
386 : const auto subdomain_constraint =
387 12212 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
388 :
389 : // This subdomain boundary node does not lie on an external boundary,
390 : // go ahead and impose constraint
391 12556 : if (boundary_node_ids.find(node.id()) == boundary_node_ids.end())
392 : {
393 5254 : if (_verbosity > 20)
394 0 : libMesh::out << "Imposing subdomain constraint on "
395 0 : << std::endl << " " << node << std::endl
396 0 : << " " << subdomain_constraint << std::endl;
397 :
398 5254 : this->impose_constraint(node, subdomain_constraint);
399 : }
400 :
401 : // This subdomain boundary node could lie on an external boundary, save it
402 : // for later to combine with the external boundary constraint.
403 : // We also save constraints for non-boundary nodes so we don't try to
404 : // re-constrain the node when accessed from the neighboring elem.
405 : // See subdomain_boundary_map.count call above.
406 24080 : subdomain_boundary_map[node.id()] = subdomain_constraint;
407 :
408 : } // for local_node_id
409 :
410 : } // for side
411 2144 : } // for elem
412 : }
413 :
414 : // Loop through boundary nodes and impose constraints
415 146970 : for (const auto & bid : boundary_node_ids)
416 : {
417 144556 : const auto & node = mesh.node_ref(bid);
418 :
419 : // Get the relevant nodal neighbors for the boundary constraint
420 : const auto side_grouped_boundary_neighbors = get_neighbors_for_boundary_constraint(
421 148628 : mesh, node, boundary_node_ids, boundary_info, nodes_to_elem_map);
422 :
423 : // Determine which constraint should be imposed
424 : const auto boundary_constraint =
425 148628 : determine_constraint(node, dim, side_grouped_boundary_neighbors);
426 :
427 : // Check for the case where this boundary node is also part of a subdomain id boundary
428 144556 : if (const auto it = subdomain_boundary_map.find(bid); it != subdomain_boundary_map.end())
429 : {
430 6958 : const auto & subdomain_constraint = it->second;
431 : // Combine current boundary constraint with previously determined
432 : // subdomain_constraint
433 392 : auto combined_constraint = intersect_constraints(subdomain_constraint, boundary_constraint);
434 :
435 : // This will catch cases where constraints have no intersection
436 : // Fall back to fixed node constraint
437 6958 : if (std::holds_alternative<InvalidConstraint>(combined_constraint))
438 0 : combined_constraint = PointConstraint(node);
439 :
440 6958 : if (_verbosity > 20)
441 0 : libMesh::out << "Imposing boundary/subdomain constraint on "
442 0 : << std::endl << " " << node << std::endl
443 0 : << " " << combined_constraint << std::endl;
444 :
445 6958 : this->impose_constraint(node, combined_constraint);
446 : }
447 :
448 : else
449 : {
450 137598 : if (_verbosity > 20)
451 0 : libMesh::out << "Imposing boundary constraint on "
452 0 : << std::endl << node << std::endl
453 0 : << " " << boundary_constraint << std::endl;
454 :
455 137598 : this->impose_constraint(node, boundary_constraint);
456 : }
457 :
458 : } // end bid
459 2414 : }
460 :
461 22436 : void VariationalSmootherConstraint::fix_node(const Node & node)
462 : {
463 84845 : for (const auto d : make_range(_sys.get_mesh().mesh_dimension()))
464 : {
465 62409 : const auto constrained_dof_index = node.dof_number(_sys.number(), d, 0);
466 3516 : DofConstraintRow constraint_row;
467 : // Leave the constraint row as all zeros so this dof is independent from other dofs
468 62409 : const auto constrained_value = node(d);
469 : // Simply constrain this dof to retain it's current value
470 62409 : _sys.get_dof_map().add_constraint_row(
471 : constrained_dof_index, constraint_row, constrained_value, true);
472 : }
473 22436 : }
474 :
475 74834 : void VariationalSmootherConstraint::constrain_node_to_plane(const Node & node,
476 : const Point & ref_normal_vec)
477 : {
478 74834 : const auto dim = _sys.get_mesh().mesh_dimension();
479 : // determine equation of plane: c_x * x + c_y * y + c_z * z + c = 0
480 4216 : std::vector<Real> xyz_coefs; // vector to hold c_x, c_y, c_z
481 2108 : Real c = 0.;
482 :
483 : // We choose to constrain the dimension with the largest magnitude coefficient
484 : // This approach ensures the coefficients added to the constraint_row
485 : // (i.e., -c_xyz / c_max) have as small magnitude as possible
486 :
487 : // We initialize this to avoid maybe-uninitialized compiler error
488 2108 : unsigned int constrained_dim = 0;
489 :
490 : // Let's assert that we have a nonzero normal to ensure that constrained_dim
491 : // is always set
492 2108 : libmesh_assert(ref_normal_vec.norm() > TOLERANCE);
493 :
494 2108 : Real max_abs_coef = 0.;
495 299336 : for (const auto d : make_range(dim))
496 : {
497 224502 : const auto coef = ref_normal_vec(d);
498 224502 : xyz_coefs.push_back(coef);
499 224502 : c -= coef * node(d);
500 :
501 6324 : const auto coef_abs = std::abs(coef);
502 224502 : if (coef_abs > max_abs_coef)
503 : {
504 2108 : max_abs_coef = coef_abs;
505 2108 : constrained_dim = d;
506 : }
507 : }
508 :
509 4216 : DofConstraintRow constraint_row;
510 299336 : for (const auto free_dim : make_range(dim))
511 : {
512 224502 : if (free_dim == constrained_dim)
513 74834 : continue;
514 :
515 149668 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
516 158100 : constraint_row[free_dof_index] = -xyz_coefs[free_dim] / xyz_coefs[constrained_dim];
517 : }
518 :
519 74834 : const auto inhomogeneous_part = -c / xyz_coefs[constrained_dim];
520 74834 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
521 74834 : _sys.get_dof_map().add_constraint_row(
522 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
523 74834 : }
524 :
525 52540 : void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
526 : const Point & line_vec)
527 : {
528 52540 : const auto dim = _sys.get_mesh().mesh_dimension();
529 :
530 : // We will free the dimension most parallel to line_vec to keep the
531 : // constraint coefficients small
532 54020 : const std::vector<Real> line_vec_coefs{line_vec(0), line_vec(1), line_vec(2)};
533 52540 : auto it = std::max_element(line_vec_coefs.begin(), line_vec_coefs.end(), [](double a, double b) {
534 2960 : return std::abs(a) < std::abs(b);
535 2960 : });
536 52540 : const unsigned int free_dim = std::distance(line_vec_coefs.begin(), it);
537 52540 : const auto free_dof_index = node.dof_number(_sys.number(), free_dim, 0);
538 :
539 : // A line is parameterized as r(t) = node + t * line_vec, so
540 : // x(t) = node(x) + t * line_vec(x)
541 : // y(t) = node(y) + t * line_vec(y)
542 : // z(t) = node(z) + t * line_vec(z)
543 : // Let's say we leave x free. Then t = (x(t) - node(x)) / line_vec(x)
544 : // Then y and z can be constrained as
545 : // y = node(y) + line_vec_y * (x(t) - node(x)) / line_vec(x)
546 : // = x(t) * line_vec(y) / line_vec(x) + (node(y) - node(x) * line_vec(y) / line_vec(x))
547 : // z = x(t) * line_vec(z) / line_vec(x) + (node(z) - node(x) * line_vec(z) / line_vec(x))
548 :
549 1480 : libmesh_assert(!relative_fuzzy_equals(line_vec(free_dim), 0.));
550 201214 : for (const auto constrained_dim : make_range(dim))
551 : {
552 148674 : if (constrained_dim == free_dim)
553 52540 : continue;
554 :
555 5416 : DofConstraintRow constraint_row;
556 96134 : constraint_row[free_dof_index] = line_vec(constrained_dim) / line_vec(free_dim);
557 : const auto inhomogeneous_part =
558 96134 : node(constrained_dim) - node(free_dim) * line_vec(constrained_dim) / line_vec(free_dim);
559 96134 : const auto constrained_dof_index = node.dof_number(_sys.number(), constrained_dim, 0);
560 96134 : _sys.get_dof_map().add_constraint_row(
561 : constrained_dof_index, constraint_row, inhomogeneous_part, true);
562 : }
563 52540 : }
564 :
565 : // Utility function to determine whether two nodes share a boundary ID.
566 : // The motivation for this is that a sliding boundary node on a triangular
567 : // element can have a neighbor boundary node in the same element that is not
568 : // part of the same boundary
569 : // Consider the below example with nodes A, C, D, F, G that comprise elements
570 : // E1, E2, E3, with boundaries B1, B2, B3, B4. To determine the constraint
571 : // equations for the sliding node C, the neighboring nodes A and D need to
572 : // be identified to construct the line that C is allowed to slide along.
573 : // Note that neighbors A and D both share the boundary B1 with C.
574 : // Without ensuring that neighbors share the same boundary as the current
575 : // node, a neighboring node that does not lie on the same boundary
576 : // (i.e. F and G) might be selected to define the constraining line,
577 : // resulting in an incorrect constraint.
578 : // Note that if, for example, boundaries B1 and B2 were to be combined
579 : // into boundary B12, the node F would share a boundary id with node C
580 : // and result in an incorrect constraint. It would be useful to design
581 : // additional checks to detect cases like this.
582 :
583 : // B3
584 : // G-----------F
585 : // | \ / |
586 : // B4 | \ E2 / | B2
587 : // | \ / |
588 : // | E1 \ / E3 |
589 : // A-----C-----D
590 : // B1
591 :
592 1034612 : bool VariationalSmootherConstraint::nodes_share_boundary_id(const Node & boundary_node,
593 : const Node & neighbor_node,
594 : const Elem & containing_elem,
595 : const BoundaryInfo & boundary_info)
596 : {
597 29144 : bool nodes_share_bid = false;
598 :
599 : // Node ids local to containing_elem
600 1005468 : const auto node_id = containing_elem.get_node_index(&boundary_node);
601 1005468 : const auto neighbor_id = containing_elem.get_node_index(&neighbor_node);
602 :
603 3151690 : for (const auto side_id : containing_elem.side_index_range())
604 : {
605 : // We don't care about this side if it doesn't contain our boundary and neighbor nodes
606 4786181 : if (!(containing_elem.is_node_on_side(node_id, side_id) &&
607 1782029 : containing_elem.is_node_on_side(neighbor_id, side_id)))
608 1568532 : continue;
609 :
610 : // If the current side, containing boundary_node and neighbor_node, lies on a boundary,
611 : // we can say that boundary_node and neighbor_node have a common boundary id.
612 40440 : std::vector<boundary_id_type> boundary_ids;
613 1435620 : boundary_info.boundary_ids(&containing_elem, side_id, boundary_ids);
614 1435620 : if (boundary_ids.size())
615 : {
616 24988 : nodes_share_bid = true;
617 24988 : break;
618 : }
619 : }
620 1034612 : return nodes_share_bid;
621 : }
622 :
623 : std::set<std::set<const Node *>>
624 12212 : VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
625 : const MeshBase & mesh,
626 : const Node & node,
627 : const subdomain_id_type sub_id,
628 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
629 : {
630 :
631 : // Find all the nodal neighbors... that is the nodes directly connected
632 : // to this node through one edge, or if none exists, use other nodes on the
633 : // containing face
634 688 : std::vector<const Node *> neighbors;
635 12212 : MeshTools::find_nodal_or_face_neighbors(
636 : mesh, node, nodes_to_elem_map, neighbors);
637 :
638 : // Each constituent set corresponds to neighbors sharing a face on the
639 : // subdomain boundary
640 344 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
641 :
642 60066 : for (const auto * neigh : neighbors)
643 : {
644 : // Determine whether the neighbor is on the subdomain boundary
645 : // First, find the common elements that both node and neigh belong to
646 47854 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
647 47854 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
648 1348 : const Elem * common_elem = nullptr;
649 242536 : for (const auto * neigh_elem : elems_containing_neigh)
650 : {
651 238188 : if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
652 10968 : elems_containing_node.end())
653 : // We should be able to find a common element on the sub_id boundary
654 194682 : && (neigh_elem->subdomain_id() == sub_id))
655 1420 : common_elem = neigh_elem;
656 : else
657 144272 : continue;
658 :
659 : // Now, determine whether node and neigh are on a side coincident
660 : // with the subdomain boundary
661 289112 : for (const auto common_side : common_elem->side_index_range())
662 : {
663 6724 : bool node_found_on_side = false;
664 6724 : bool neigh_found_on_side = false;
665 2024544 : for (const auto local_node_id : common_elem->nodes_on_side(common_side))
666 : {
667 1829234 : if (common_elem->node_id(local_node_id) == node.id())
668 2872 : node_found_on_side = true;
669 1677162 : else if (common_elem->node_id(local_node_id) == neigh->id())
670 3436 : neigh_found_on_side = true;
671 : }
672 :
673 240946 : if (!(node_found_on_side && neigh_found_on_side &&
674 4488 : common_elem->neighbor_ptr(common_side)))
675 172776 : continue;
676 :
677 1716 : const auto matched_side = common_side;
678 : // There could be multiple matched sides, so keep this next part
679 : // inside the common_side loop
680 : //
681 : // Does matched_side, containing both node and neigh, lie on the
682 : // sub_id subdomain boundary?
683 : const auto matched_neighbor_sub_id =
684 3432 : common_elem->neighbor_ptr(matched_side)->subdomain_id();
685 1716 : const bool is_matched_side_on_subdomain_boundary = matched_neighbor_sub_id != sub_id;
686 :
687 60918 : if (is_matched_side_on_subdomain_boundary)
688 : {
689 : // Store all nodes that live on this side
690 43800 : const auto nodes_on_side = common_elem->nodes_on_side(common_side);
691 2400 : std::set<const Node *> node_ptrs_on_side;
692 361816 : for (const auto local_node_id : nodes_on_side)
693 328208 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
694 43800 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
695 41400 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
696 :
697 1200 : continue;
698 : }
699 :
700 : } // for common_side
701 :
702 : } // for neigh_elem
703 : }
704 :
705 12556 : return side_grouped_boundary_neighbors;
706 : }
707 :
708 : std::set<std::set<const Node *>>
709 144556 : VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
710 : const MeshBase & mesh,
711 : const Node & node,
712 : const std::unordered_set<dof_id_type> & boundary_node_ids,
713 : const BoundaryInfo & boundary_info,
714 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
715 : {
716 :
717 : // Find all the nodal neighbors... that is the nodes directly connected
718 : // to this node through one edge, or if none exists, use other nodes on the
719 : // containing face
720 8144 : std::vector<const Node *> neighbors;
721 144556 : MeshTools::find_nodal_or_face_neighbors(
722 : mesh, node, nodes_to_elem_map, neighbors);
723 :
724 : // Each constituent set corresponds to neighbors sharing a face on the
725 : // boundary
726 4072 : std::set<std::set<const Node *>> side_grouped_boundary_neighbors;
727 :
728 815364 : for (const auto * neigh : neighbors)
729 : {
730 : const bool is_neighbor_boundary_node =
731 670808 : boundary_node_ids.find(neigh->id()) != boundary_node_ids.end();
732 670808 : if (!is_neighbor_boundary_node)
733 124338 : continue;
734 :
735 : // Determine whether nodes share a common boundary id
736 : // First, find the common element that both node and neigh belong to
737 542866 : const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
738 542866 : const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
739 15292 : const Elem * common_elem = nullptr;
740 3368240 : for (const auto * neigh_elem : elems_containing_neigh)
741 : {
742 : const bool is_neigh_common =
743 2904962 : std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
744 159176 : elems_containing_node.end();
745 2825374 : if (!is_neigh_common)
746 1740318 : continue;
747 58288 : common_elem = neigh_elem;
748 : // Keep this in the neigh_elem loop because there can be multiple common
749 : // elements Now, determine whether node and neigh share a common boundary
750 : // id
751 1034612 : const bool nodes_have_common_bid = VariationalSmootherConstraint::nodes_share_boundary_id(
752 : node, *neigh, *common_elem, boundary_info);
753 1034612 : if (nodes_have_common_bid)
754 : {
755 : // Find the side coinciding with the shared boundary
756 5197910 : for (const auto side : common_elem->side_index_range())
757 : {
758 : // We only care about external boundaries here, make sure side doesn't
759 : // have a neighbor
760 4432268 : if (common_elem->neighbor_ptr(side))
761 3355744 : continue;
762 :
763 40440 : bool node_found_on_side = false;
764 40440 : bool neigh_found_on_side = false;
765 1435620 : const auto nodes_on_side = common_elem->nodes_on_side(side);
766 10338452 : for (const auto local_node_id : nodes_on_side)
767 : {
768 9153616 : if (common_elem->node_id(local_node_id) == node.id())
769 29120 : node_found_on_side = true;
770 7869072 : else if (common_elem->node_id(local_node_id) == neigh->id())
771 30488 : neigh_found_on_side = true;
772 : }
773 1435620 : if (!(node_found_on_side && neigh_found_on_side))
774 13536 : continue;
775 :
776 53808 : std::set<const Node *> node_ptrs_on_side;
777 6418684 : for (const auto local_node_id : nodes_on_side)
778 5617496 : node_ptrs_on_side.insert(common_elem->node_ptr(local_node_id));
779 981996 : node_ptrs_on_side.erase(node_ptrs_on_side.find(&node));
780 928188 : side_grouped_boundary_neighbors.insert(node_ptrs_on_side);
781 : }
782 862086 : continue;
783 837098 : }
784 : }
785 : }
786 :
787 148628 : return side_grouped_boundary_neighbors;
788 : }
789 :
790 156768 : ConstraintVariant VariationalSmootherConstraint::determine_constraint(
791 : const Node & node,
792 : const unsigned int dim,
793 : const std::set<std::set<const Node *>> & side_grouped_boundary_neighbors)
794 : {
795 : // Determines the appropriate geometric constraint for a node based on its
796 : // neighbors.
797 :
798 : // Extract neighbors in flat vector
799 8832 : std::vector<const Node *> neighbors;
800 623522 : for (const auto & side : side_grouped_boundary_neighbors)
801 466754 : neighbors.insert(neighbors.end(), side.begin(), side.end());
802 :
803 : // Constrain the node to it's current location
804 156768 : if (dim == 1 || neighbors.size() == 1)
805 1125 : return PointConstraint{node};
806 :
807 155703 : if (dim == 2 || neighbors.size() == 2)
808 : {
809 : // Determine whether node + all neighbors are colinear
810 370 : bool all_colinear = true;
811 13505 : const Point ref_dir = (*neighbors[0] - node).unit();
812 26696 : for (const auto i : make_range(std::size_t(1), neighbors.size()))
813 : {
814 17082 : const Point delta = *(neighbors[i]) - node;
815 468 : libmesh_assert(delta.norm() >= TOLERANCE);
816 16614 : const Point dir = delta.unit();
817 17008 : if (!dir.relative_fuzzy_equals(ref_dir) && !dir.relative_fuzzy_equals(-ref_dir))
818 : {
819 86 : all_colinear = false;
820 3053 : break;
821 : }
822 : }
823 :
824 370 : if (all_colinear)
825 10650 : return LineConstraint{node, ref_dir};
826 :
827 3225 : return PointConstraint{node};
828 : }
829 :
830 : // dim == 3, neighbors.size() >= 3
831 8032 : std::set<PlaneConstraint> valid_planes;
832 585324 : for (const auto & side_nodes : side_grouped_boundary_neighbors)
833 : {
834 455228 : std::vector<const Node *> side_nodes_vec(side_nodes.begin(), side_nodes.end());
835 2496644 : for (const auto i : index_range(side_nodes_vec))
836 : {
837 2111744 : const Point vec_i = (*side_nodes_vec[i] - node);
838 6983560 : for (const auto j : make_range(i))
839 : {
840 5068536 : const Point vec_j = (*side_nodes_vec[j] - node);
841 4790808 : Point candidate_normal = vec_i.cross(vec_j);
842 4929672 : if (candidate_normal.norm() <= TOLERANCE)
843 320565 : continue;
844 :
845 4868775 : const PlaneConstraint candidate_plane{node, candidate_normal};
846 4479273 : valid_planes.emplace(candidate_plane);
847 : }
848 : }
849 : }
850 :
851 : // Fall back to point constraint
852 142568 : if (valid_planes.empty())
853 0 : return PointConstraint(node);
854 :
855 : // Combine all the planes together to get a common constraint
856 4016 : auto it = valid_planes.begin();
857 8032 : ConstraintVariant current = *it++;
858 301413 : for (; it != valid_planes.end(); ++it)
859 : {
860 163463 : current = intersect_constraints(current, *it);
861 :
862 : // This will catch cases where constraints have no intersection
863 : // (i.e., the element surface is non-planar)
864 : // Fall back to fixed node constraint
865 158987 : if (std::holds_alternative<InvalidConstraint>(current))
866 : {
867 146 : current = PointConstraint(node);
868 142 : break;
869 : }
870 : }
871 :
872 4016 : return current;
873 : }
874 :
875 : // Applies the computed constraint (PointConstraint, LineConstraint, or
876 : // PlaneConstraint) to a node.
877 149810 : void VariationalSmootherConstraint::impose_constraint(const Node & node,
878 : const ConstraintVariant & constraint)
879 : {
880 :
881 4220 : libmesh_assert_msg(!std::holds_alternative<InvalidConstraint>(constraint),
882 : "Cannot impose constraint using InvalidConstraint.");
883 :
884 149810 : if (std::holds_alternative<PointConstraint>(constraint))
885 22436 : fix_node(node);
886 127374 : else if (std::holds_alternative<LineConstraint>(constraint))
887 52540 : constrain_node_to_line(node, std::get<LineConstraint>(constraint).direction());
888 74834 : else if (std::holds_alternative<PlaneConstraint>(constraint))
889 74834 : constrain_node_to_plane(node, std::get<PlaneConstraint>(constraint).normal());
890 :
891 : else
892 0 : libmesh_assert_msg(false, "Unknown constraint type.");
893 149810 : }
894 :
895 : } // namespace libMesh
|
