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