Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "MortarSegmentInfo.h"
13 : #include "MooseHashing.h"
14 : #include "ConsoleStreamInterface.h"
15 : #include "MooseError.h"
16 : #include "MooseUtils.h"
17 :
18 : // libMesh includes
19 : #include "libmesh/id_types.h"
20 : #include "libmesh/equation_systems.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/int_range.h"
23 :
24 : // C++ includes
25 : #include <optional>
26 : #include <set>
27 : #include <memory>
28 : #include <vector>
29 : #include <unordered_map>
30 :
31 : // Forward declarations
32 : namespace libMesh
33 : {
34 : class MeshBase;
35 : class System;
36 : }
37 : class GetPot;
38 :
39 : // Using statements
40 : using libMesh::boundary_id_type;
41 : using libMesh::CompareDofObjectsByID;
42 : using libMesh::dof_id_type;
43 : using libMesh::Elem;
44 : using libMesh::MeshBase;
45 : using libMesh::Node;
46 : using libMesh::Point;
47 : using libMesh::Real;
48 : using libMesh::subdomain_id_type;
49 :
50 : typedef boundary_id_type BoundaryID;
51 : typedef subdomain_id_type SubdomainID;
52 :
53 : /**
54 : * This class is a container/interface for the objects involved in
55 : * automatic generation of mortar spaces.
56 : */
57 : class AutomaticMortarGeneration : public ConsoleStreamInterface
58 : {
59 : public:
60 : /**
61 : * The name of the nodal normals system. We store this in one place
62 : * so it's easy to change later.
63 : */
64 : const static std::string system_name;
65 :
66 : /**
67 : * Must be constructed with a reference to the Mesh we are
68 : * generating mortar spaces for.
69 : */
70 : AutomaticMortarGeneration(MooseApp & app,
71 : MeshBase & mesh_in,
72 : const std::pair<BoundaryID, BoundaryID> & boundary_key,
73 : const std::pair<SubdomainID, SubdomainID> & subdomain_key,
74 : bool on_displaced,
75 : bool periodic,
76 : const bool debug,
77 : const bool correct_edge_dropping,
78 : const Real minimum_projection_angle,
79 : const MortarSegmentTriangulationMode triangulation_mode,
80 : const bool triangulate_triangles);
81 :
82 : /**
83 : * Once the secondary_requested_boundary_ids and
84 : * primary_requested_boundary_ids containers have been filled in,
85 : * call this function to build node-to-Elem maps for the
86 : * lower-dimensional elements.
87 : */
88 : void buildNodeToElemMaps();
89 :
90 : /**
91 : * Computes and stores the nodal normal/tangent vectors in a local data
92 : * structure instead of using the ExplicitSystem/NumericVector
93 : * approach. This design was triggered by the way that the
94 : * GhostingFunctor operates, but I think it is a better/more
95 : * efficient way to do it anyway.
96 : */
97 : void computeNodalGeometry();
98 :
99 : /**
100 : * Project secondary nodes (find xi^(2) values) to the closest points on
101 : * the primary surface.
102 : * Inputs:
103 : * - The nodal normals values
104 : * - mesh
105 : * - nodes_to_primary_elem_map
106 : *
107 : * Outputs:
108 : * - secondary_node_and_elem_to_xi2_primary_elem
109 : *
110 : * Defined in the file project_secondary_nodes.C.
111 : */
112 : void projectSecondaryNodes();
113 :
114 : /**
115 : * (Inverse) project primary nodes to the points on the secondary surface
116 : * where they would have come from (find (xi^(1) values)).
117 : *
118 : * Inputs:
119 : * - The nodal normals values
120 : * - mesh
121 : * - nodes_to_secondary_elem_map
122 : *
123 : * Outputs:
124 : * - primary_node_and_elem_to_xi1_secondary_elem
125 : *
126 : * Defined in the file project_primary_nodes.C.
127 : */
128 : void projectPrimaryNodes();
129 :
130 : /**
131 : * Builds the mortar segment mesh once the secondary and primary node
132 : * projections have been completed.
133 : *
134 : * Inputs:
135 : * - mesh
136 : * - primary_node_and_elem_to_xi1_secondary_elem
137 : * - secondary_node_and_elem_to_xi2_primary_elem
138 : * - nodes_to_primary_elem_map
139 : *
140 : * Outputs:
141 : * - mortar_segment_mesh
142 : * - msm_elem_to_info
143 : *
144 : * Defined in the file build_mortar_segment_mesh.C.
145 : */
146 : void buildMortarSegmentMesh();
147 :
148 : /**
149 : * Builds the mortar segment mesh once the secondary and primary node
150 : * projections have been completed.
151 : *
152 : * Inputs:
153 : * - mesh
154 : *
155 : * Outputs:
156 : * - mortar_segment_mesh
157 : * - msm_elem_to_info
158 : */
159 : void buildMortarSegmentMesh3d();
160 :
161 : /**
162 : * Statistics for one primary-secondary subdomain pair.
163 : * Secondary/primary lower-d stats reflect local data (all data for replicated meshes).
164 : */
165 : struct MsmSubdomainStats
166 : {
167 : SubdomainID primary_subd_id;
168 : SubdomainID secondary_subd_id;
169 : std::size_t secondary_lower_n_elems;
170 : Real secondary_lower_max_volume;
171 : Real secondary_lower_min_volume;
172 : Real secondary_lower_median_volume;
173 : std::size_t primary_lower_n_elems;
174 : Real primary_lower_max_volume;
175 : Real primary_lower_min_volume;
176 : Real primary_lower_median_volume;
177 : std::size_t msm_n_elems;
178 : Real msm_max_volume;
179 : Real msm_min_volume;
180 : Real msm_median_volume;
181 : };
182 :
183 : /**
184 : * Computes mortar segment mesh statistics and returns one entry per subdomain pair.
185 : * Must be called collectively on all ranks.
186 : */
187 : std::vector<MsmSubdomainStats> computeMsmStatistics();
188 :
189 : /**
190 : * Prints mortar segment mesh statistics to console (calls computeMsmStatistics internally)
191 : */
192 : void msmStatistics();
193 :
194 : /**
195 : * Clears the mortar segment mesh and accompanying data structures
196 : */
197 : void clear();
198 :
199 : /**
200 : * Invalidates the cached MSM node/element ID starting offset so that the next call to
201 : * buildMortarSegmentMesh3d() recomputes it via allgather. Call this when mesh topology changes.
202 : */
203 172 : void meshChanged() { _msm_node_id_start = std::nullopt; }
204 :
205 : /**
206 : * returns whether this object is on the displaced mesh
207 : */
208 : bool onDisplaced() const { return _on_displaced; }
209 :
210 : /**
211 : * @return The nodal normals associated with the provided \p secondary_elem
212 : */
213 : std::vector<Point> getNodalNormals(const Elem & secondary_elem) const;
214 :
215 : /**
216 : * Compute the two nodal tangents, which are built on-the-fly.
217 : * @return The nodal tangents associated with the provided \p secondary_elem
218 : */
219 : std::array<MooseUtils::SemidynamicVector<Point, 9>, 2>
220 : getNodalTangents(const Elem & secondary_elem) const;
221 :
222 : /**
223 : * Compute on-the-fly mapping from secondary interior parent nodes to lower dimensional nodes
224 : * @return The map from secondary interior parent nodes to lower dimensional nodes
225 : */
226 : std::map<unsigned int, unsigned int>
227 : getSecondaryIpToLowerElementMap(const Elem & lower_secondary_elem) const;
228 :
229 : /**
230 : * Compute on-the-fly mapping from primary interior parent nodes to its corresponding lower
231 : * dimensional nodes
232 : * @return The map from primary interior parent nodes to its corresponding lower dimensional
233 : * nodes
234 : */
235 : std::map<unsigned int, unsigned int>
236 : getPrimaryIpToLowerElementMap(const Elem & primary_elem,
237 : const Elem & primary_elem_ip,
238 : const Elem & lower_secondary_elem) const;
239 :
240 : /**
241 : * Compute the normals at given reference points on a secondary element
242 : * @param secondary_elem The secondary element used to query for associated nodal normals
243 : * @param xi1_pts The reference points on the secondary element to evaluate the normals at. The
244 : * points should only be non-zero in the zeroth entry because right now our mortar mesh elements
245 : * are always 1D
246 : * @return The normals
247 : */
248 : std::vector<Point> getNormals(const Elem & secondary_elem,
249 : const std::vector<Point> & xi1_pts) const;
250 :
251 : /**
252 : * Compute the normals at given reference points on a secondary element
253 : * @param secondary_elem The secondary element used to query for associated nodal normals
254 : * @param 1d_xi1_pts The reference points on the secondary element to evaluate the normals at. The
255 : * "points" are single reals corresponding to xi because right now our mortar mesh elements are
256 : * always 1D
257 : * @return The normals
258 : */
259 : std::vector<Point> getNormals(const Elem & secondary_elem,
260 : const std::vector<Real> & oned_xi1_pts) const;
261 :
262 : /**
263 : * Return lower dimensional secondary element given its interior parent. Helpful outside the
264 : * mortar generation to locate mortar-related quantities.
265 : * @param secondary_elem_id The secondary interior parent element id used to query for associated
266 : * lower dimensional element
267 : * @return The corresponding lower dimensional secondary element
268 : */
269 : const Elem * getSecondaryLowerdElemFromSecondaryElem(dof_id_type secondary_elem_id) const;
270 :
271 : /**
272 : * Get list of secondary nodes that don't contribute to interaction with any primary element.
273 : * Used to enforce zero values on inactive DoFs of nodal variables.
274 : */
275 : void computeInactiveLMNodes();
276 :
277 : /**
278 : * Computes inactive secondary nodes when incorrect edge dropping behavior is enabled
279 : * (any node touching a partially or fully dropped element is dropped)
280 : */
281 : void computeIncorrectEdgeDroppingInactiveLMNodes();
282 :
283 : /**
284 : * Get list of secondary elems without any corresponding primary elements.
285 : * Used to enforce zero values on inactive DoFs of elemental variables.
286 : */
287 : void computeInactiveLMElems();
288 :
289 : /**
290 : * @return The mortar interface coupling
291 : */
292 : const std::unordered_map<dof_id_type, std::unordered_set<dof_id_type>> &
293 405074 : mortarInterfaceCoupling() const
294 : {
295 405074 : return _mortar_interface_coupling;
296 : }
297 :
298 : /**
299 : * @return The primary-secondary boundary ID pair
300 : */
301 : const std::pair<BoundaryID, BoundaryID> & primarySecondaryBoundaryIDPair() const;
302 :
303 : /**
304 : * @return The mortar segment mesh
305 : */
306 2008 : const MeshBase & mortarSegmentMesh() const { return *_mortar_segment_mesh; }
307 :
308 : /**
309 : * @return The mortar segment element to corresponding information
310 : */
311 949900 : const std::unordered_map<const Elem *, MortarSegmentInfo> & mortarSegmentMeshElemToInfo() const
312 : {
313 949900 : return _msm_elem_to_info;
314 : }
315 :
316 23165 : int dim() const { return _mesh.mesh_dimension(); }
317 :
318 : /**
319 : * @return The set of nodes on which mortar constraints are not active
320 : */
321 21771 : const std::unordered_set<const Node *> & getInactiveLMNodes() const
322 : {
323 21771 : return _inactive_local_lm_nodes;
324 : }
325 :
326 : /**
327 : * @return The list of secondary elems on which mortar constraint is not active
328 : */
329 15321 : const std::unordered_set<const Elem *> & getInactiveLMElems() const
330 : {
331 15321 : return _inactive_local_lm_elems;
332 : }
333 :
334 15321 : bool incorrectEdgeDropping() const { return !_correct_edge_dropping; }
335 :
336 : using MortarFilterIter =
337 : std::unordered_map<dof_id_type, std::set<Elem *, CompareDofObjectsByID>>::const_iterator;
338 :
339 : /**
340 : * @return A vector of iterators that point to the lower dimensional secondary elements and their
341 : * associated mortar segment elements that would have nonzero values for a Lagrange shape function
342 : * associated with the provided node. This method may return an empty container if the node is
343 : * away from the mortar mesh
344 : */
345 : std::vector<MortarFilterIter> secondariesToMortarSegments(const Node & node) const;
346 :
347 : /**
348 : * @return the lower dimensional secondary element ids and their associated mortar segment
349 : * elements
350 : */
351 : const std::unordered_map<dof_id_type, std::set<Elem *, CompareDofObjectsByID>> &
352 12651 : secondariesToMortarSegments() const
353 : {
354 12651 : return _secondary_elems_to_mortar_segments;
355 : }
356 :
357 : /**
358 : * @return All the secondary interior parent subdomain IDs associated with the mortar mesh
359 : */
360 1248 : const std::set<SubdomainID> & secondaryIPSubIDs() const { return _secondary_ip_sub_ids; }
361 :
362 : /**
363 : * @return All the primary interior parent subdomain IDs associated with the mortar mesh
364 : */
365 1248 : const std::set<SubdomainID> & primaryIPSubIDs() const { return _primary_ip_sub_ids; }
366 :
367 : /**
368 : * @return Map from node id to secondary lower-d element pointer
369 : */
370 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodesToSecondaryElem() const
371 : {
372 : return _nodes_to_secondary_elem_map;
373 : }
374 :
375 : /**
376 : * initialize mortar-mesh based output
377 : */
378 : void initOutput();
379 :
380 : private:
381 : /**
382 : * Write the mortar segment mesh to exodus
383 : */
384 : void outputMortarMesh();
385 :
386 : /**
387 : * @returns A string uniquely identifying this mortar interface
388 : */
389 : std::string mortarInterfaceName() const;
390 :
391 : /**
392 : * build the \p _mortar_interface_coupling data
393 : */
394 : void buildCouplingInformation();
395 :
396 : /// The Moose app
397 : MooseApp & _app;
398 :
399 : /// Reference to the mesh stored in equation_systems.
400 : MeshBase & _mesh;
401 :
402 : /// The boundary ids corresponding to all the secondary surfaces.
403 : std::set<BoundaryID> _secondary_requested_boundary_ids;
404 :
405 : /// The boundary ids corresponding to all the primary surfaces.
406 : std::set<BoundaryID> _primary_requested_boundary_ids;
407 :
408 : /// A list of primary/secondary boundary id pairs corresponding to each
409 : /// side of the mortar interface.
410 : std::vector<std::pair<BoundaryID, BoundaryID>> _primary_secondary_boundary_id_pairs;
411 :
412 : /// Map from nodes to connected lower-dimensional elements on the secondary/primary subdomains.
413 : std::unordered_map<dof_id_type, std::vector<const Elem *>> _nodes_to_secondary_elem_map;
414 : std::unordered_map<dof_id_type, std::vector<const Elem *>> _nodes_to_primary_elem_map;
415 :
416 : /// Similar to the map above, but associates a (Secondary Node, Secondary Elem)
417 : /// pair to a (xi^(2), primary Elem) pair. This allows a single secondary node, which is
418 : /// potentially connected to two elements on the secondary side, to be associated with
419 : /// multiple primary Elem/xi^(2) values to handle the case where the primary and secondary
420 : /// nodes are "matching".
421 : /// In this configuration:
422 : ///
423 : /// A B
424 : /// o-----o-----o (secondary orientation ->)
425 : /// |
426 : /// v
427 : /// ------x------ (primary orientation <-)
428 : /// C D
429 : ///
430 : /// The entries in the map should be:
431 : /// (Elem A, Node 1) -> (Elem C, xi^(2)=-1)
432 : /// (Elem B, Node 0) -> (Elem D, xi^(2)=+1)
433 : std::unordered_map<std::pair<const Node *, const Elem *>, std::pair<Real, const Elem *>>
434 : _secondary_node_and_elem_to_xi2_primary_elem;
435 :
436 : /// Same type of container, but for mapping (Primary Node ID, Primary Node,
437 : /// Primary Elem) -> (xi^(1), Secondary Elem) where they are inverse-projected along
438 : /// the nodal normal direction. Note that the first item of the key, the primary
439 : /// node ID, is important for storing the key-value pairs in a consistent order
440 : /// across processes, e.g. this container has to be ordered!
441 : std::map<std::tuple<dof_id_type, const Node *, const Elem *>, std::pair<Real, const Elem *>>
442 : _primary_node_and_elem_to_xi1_secondary_elem;
443 :
444 : /// 1D Mesh of mortar segment elements which gets built by the call
445 : /// to build_mortar_segment_mesh().
446 : std::unique_ptr<MeshBase> _mortar_segment_mesh;
447 :
448 : /// Map between Elems in the mortar segment mesh and their info
449 : /// structs. This gets filled in by the call to
450 : /// build_mortar_segment_mesh().
451 : std::unordered_map<const Elem *, MortarSegmentInfo> _msm_elem_to_info;
452 :
453 : /// Keeps track of the mapping between lower-dimensional elements and
454 : /// the side_id of the interior_parent which they are.
455 : std::unordered_map<const Elem *, unsigned int> _lower_elem_to_side_id;
456 :
457 : /// A list of primary/secondary subdomain id pairs corresponding to each
458 : /// side of the mortar interface.
459 : std::vector<std::pair<SubdomainID, SubdomainID>> _primary_secondary_subdomain_id_pairs;
460 :
461 : /// The secondary/primary lower-dimensional boundary subdomain ids are the
462 : /// secondary/primary *boundary* ids
463 : std::set<SubdomainID> _secondary_boundary_subdomain_ids;
464 : std::set<SubdomainID> _primary_boundary_subdomain_ids;
465 :
466 : /// Used by the AugmentSparsityOnInterface functor to determine
467 : /// whether a given Elem is coupled to any others across the gap, and
468 : /// to explicitly set up the dependence between interior_parent()
469 : /// elements on the secondary side and their lower-dimensional sides
470 : /// which are on the interface. This latter type of coupling must be
471 : /// explicitly declared when there is no primary_elem for a given
472 : /// mortar segment and you are using e.g. a P^1-P^0 discretization
473 : /// which does not induce the coupling automatically.
474 : std::unordered_map<dof_id_type, std::unordered_set<dof_id_type>> _mortar_interface_coupling;
475 :
476 : /// Container for storing the nodal normal vector associated with each secondary node.
477 : std::unordered_map<const Node *, Point> _secondary_node_to_nodal_normal;
478 :
479 : /// Container for storing the nodal tangent/binormal vectors associated with each secondary node
480 : /// (Householder approach).
481 : std::unordered_map<const Node *, std::array<Point, 2>> _secondary_node_to_hh_nodal_tangents;
482 :
483 : /// Map from full dimensional secondary element id to lower dimensional secondary element
484 : std::unordered_map<dof_id_type, const Elem *> _secondary_element_to_secondary_lowerd_element;
485 :
486 : // List of inactive lagrange multiplier nodes (for nodal variables)
487 : std::unordered_set<const Node *> _inactive_local_lm_nodes;
488 :
489 : /// List of inactive lagrange multiplier nodes (for elemental variables)
490 : std::unordered_set<const Elem *> _inactive_local_lm_elems;
491 :
492 : /// We maintain a mapping from lower-dimensional secondary elements in the original mesh to (sets
493 : /// of) elements in mortar_segment_mesh. This allows us to quickly determine which elements need
494 : /// to be split.
495 : std::unordered_map<dof_id_type, std::set<Elem *, CompareDofObjectsByID>>
496 : _secondary_elems_to_mortar_segments;
497 :
498 : /// All the secondary interior parent subdomain IDs associated with the mortar mesh
499 : std::set<SubdomainID> _secondary_ip_sub_ids;
500 :
501 : /// All the primary interior parent subdomain IDs associated with the mortar mesh
502 : std::set<SubdomainID> _primary_ip_sub_ids;
503 :
504 : /**
505 : * Helper function responsible for projecting secondary nodes
506 : * onto primary elements for a single primary/secondary pair. Called by the class member
507 : * AutomaticMortarGeneration::project_secondary_nodes().
508 : */
509 : void projectSecondaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id,
510 : SubdomainID lower_dimensional_secondary_subdomain_id);
511 :
512 : /**
513 : * Helper function used internally by AutomaticMortarGeneration::project_primary_nodes().
514 : */
515 : void projectPrimaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id,
516 : SubdomainID lower_dimensional_secondary_subdomain_id);
517 :
518 : /**
519 : * Householder orthogonalization procedure to obtain proper basis for tangent and binormal vectors
520 : */
521 : void
522 : householderOrthogolization(const Point & normal, Point & tangent_one, Point & tangent_two) const;
523 :
524 : /**
525 : * Process aligned nodes
526 : * @returns whether mortar segment(s) were created
527 : */
528 : bool processAlignedNodes(const Node & secondary_node,
529 : const Node & primary_node,
530 : const std::vector<const Elem *> * secondary_node_neighbors,
531 : const std::vector<const Elem *> * primary_node_neighbors,
532 : const VectorValue<Real> & nodal_normal,
533 : const Elem & candidate_element,
534 : std::set<const Elem *> & rejected_element_candidates);
535 :
536 : /// Whether to print debug output
537 : const bool _debug;
538 :
539 : /// Whether this object is on the displaced mesh
540 : const bool _on_displaced;
541 :
542 : /// Whether this object will be generating a mortar segment mesh for periodic constraints
543 : const bool _periodic;
544 :
545 : /// Whether the mortar segment mesh is distributed
546 : const bool _distributed;
547 :
548 : /// Newton solve tolerance for node projections
549 : Real _newton_tolerance = 1e-12;
550 :
551 : /// Tolerance for checking projection xi values. Usually we are checking whether we projected onto
552 : /// a certain element (in which case -1 <= xi <= 1) or whether we should have *already* projected
553 : /// a primary node (in which case we error if abs(xi) is sufficiently close to 1)
554 : Real _xi_tolerance = 1e-6;
555 :
556 : /// Flag to enable regressed treatment of edge dropping where all LM DoFs on edge dropping element
557 : /// are strongly set to 0.
558 : const bool _correct_edge_dropping;
559 :
560 : /// Parameter to control which angle (in degrees) is admissible for the creation of mortar segments.
561 : /// If set to a value close to zero, very oblique projections are allowed, which can result in mortar
562 : /// segments solving physics not meaningfully and overprojection of primary nodes onto the mortar
563 : /// segment mesh in extreme cases. This parameter is mostly intended for mortar mesh debugging purposes in 2D.
564 : const Real _minimum_projection_angle;
565 :
566 : /// Triangulation mode used for clipped 3D mortar polygons.
567 : const MortarSegmentTriangulationMode _triangulation_mode;
568 :
569 : /// Whether already-triangular clipped polygons should still be centroid-subdivided.
570 : const bool _triangulate_triangles;
571 :
572 : /// Storage for the input parameters used by the mortar nodal geometry output
573 : std::unique_ptr<InputParameters> _output_params;
574 :
575 : /// Cached per-rank starting ID for 3D MSM nodes/elements. nullopt forces recomputation on next
576 : /// buildMortarSegmentMesh3d() call. Reset by meshChanged() on topology change so the allgather
577 : /// is skipped for displaced-mesh residual updates that only move nodes.
578 : std::optional<dof_id_type> _msm_node_id_start;
579 :
580 : /// Debugging container for printing information about fraction of successful projections for
581 : /// secondary nodes. If !_debug then this should always be empty
582 : std::unordered_set<dof_id_type> _projected_secondary_nodes;
583 :
584 : /// Secondary nodes that failed to project
585 : std::unordered_set<dof_id_type> _failed_secondary_node_projections;
586 :
587 : friend class MortarNodalGeometryOutput;
588 : friend class AugmentSparsityOnInterface;
589 : };
590 :
591 : inline const std::pair<BoundaryID, BoundaryID> &
592 14494 : AutomaticMortarGeneration::primarySecondaryBoundaryIDPair() const
593 : {
594 : mooseAssert(_primary_secondary_boundary_id_pairs.size() == 1,
595 : "We currently only support a single boundary pair per mortar generation object");
596 :
597 14494 : return _primary_secondary_boundary_id_pairs.front();
598 : }
|
