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 : #include "AutomaticMortarGeneration.h"
11 : #include "MortarSegmentInfo.h"
12 : #include "NanoflannMeshAdaptor.h"
13 : #include "MooseError.h"
14 : #include "MooseTypes.h"
15 : #include "MooseLagrangeHelpers.h"
16 : #include "MortarSegmentHelper.h"
17 : #include "MortarUtils.h"
18 : #include "FormattedTable.h"
19 : #include "FEProblemBase.h"
20 : #include "DisplacedProblem.h"
21 : #include "Output.h"
22 :
23 : #include "libmesh/mesh_tools.h"
24 : #include "libmesh/explicit_system.h"
25 : #include "libmesh/numeric_vector.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/node.h"
28 : #include "libmesh/dof_map.h"
29 : #include "libmesh/edge_edge2.h"
30 : #include "libmesh/edge_edge3.h"
31 : #include "libmesh/face_tri3.h"
32 : #include "libmesh/face_tri6.h"
33 : #include "libmesh/face_tri7.h"
34 : #include "libmesh/face_quad4.h"
35 : #include "libmesh/face_quad8.h"
36 : #include "libmesh/face_quad9.h"
37 : #include "libmesh/exodusII_io.h"
38 : #include "libmesh/quadrature_gauss.h"
39 : #include "libmesh/quadrature_nodal.h"
40 : #include "libmesh/distributed_mesh.h"
41 : #include "libmesh/replicated_mesh.h"
42 : #include "libmesh/enum_to_string.h"
43 : #include "libmesh/statistics.h"
44 : #include "libmesh/equation_systems.h"
45 :
46 : #include "metaphysicl/dualnumber.h"
47 :
48 : #include "timpi/communicator.h"
49 : #include "timpi/parallel_sync.h"
50 :
51 : #include <array>
52 : #include <algorithm>
53 : #include <cmath>
54 : #include <limits>
55 :
56 : using namespace libMesh;
57 : using MetaPhysicL::DualNumber;
58 :
59 : // Make newer nanoflann API spelling compatible with older nanoflann
60 : // versions
61 : #if NANOFLANN_VERSION < 0x150
62 : namespace nanoflann
63 : {
64 : typedef SearchParams SearchParameters;
65 : }
66 : #endif
67 :
68 : namespace
69 : {
70 : // QNodal on a parent side returns normals, weights, and physical points in the
71 : // parent-side quadrature ordering. That ordering is not guaranteed to match the
72 : // node ordering of the generated lower-dimensional secondary element, especially
73 : // for higher-order faces. Build the association geometrically so each
74 : // quadrature value is attached to the secondary node at the same physical point.
75 : std::vector<unsigned int>
76 28401 : nodalQuadraturePointToSecondaryNodeMap(const Elem & secondary_elem,
77 : const std::vector<Point> & q_points)
78 : {
79 28401 : const auto n_nodes = secondary_elem.n_nodes();
80 28401 : if (q_points.size() != n_nodes)
81 0 : mooseError("Nodal quadrature produced ",
82 0 : q_points.size(),
83 : " points for secondary mortar element ",
84 0 : secondary_elem.id(),
85 : " of type ",
86 0 : libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type()),
87 : ", but the element has ",
88 : n_nodes,
89 : " nodes.");
90 :
91 28401 : const auto invalid_node = std::numeric_limits<unsigned int>::max();
92 56802 : std::vector<unsigned int> qpoint_to_node(n_nodes, invalid_node);
93 28401 : std::vector<bool> node_used(n_nodes, false);
94 :
95 28401 : const Real element_size = secondary_elem.hmax();
96 : mooseAssert(element_size > 0,
97 : "Secondary mortar element "
98 : << secondary_elem.id() << " of type "
99 : << libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type())
100 : << " has a non-positive hmax and cannot be used for nodal quadrature point "
101 : "matching.");
102 :
103 : // The nodal quadrature locations and the generated secondary nodes are two floating-point
104 : // reconstructions of the same physical points. Scale the tolerance by element size so the
105 : // matching is insensitive to coordinate magnitude; the 100*TOLERANCE factor allows roundoff
106 : // from FE reinitialization and mesh generation while remaining far below a valid node spacing.
107 28401 : const Real matching_tol = 100 * TOLERANCE * element_size;
108 28401 : const Real matching_tol_sq = matching_tol * matching_tol;
109 :
110 : // Each nodal quadrature point should coincide with exactly one still-unused
111 : // secondary node. The unused-node search makes the mapping one-to-one and
112 : // avoids silently assigning two quadrature entries to the same node.
113 112905 : for (const auto qp : make_range(q_points.size()))
114 : {
115 84504 : unsigned int closest_node = invalid_node;
116 84504 : Real closest_dist_sq = std::numeric_limits<Real>::max();
117 84504 : Real second_closest_dist_sq = std::numeric_limits<Real>::max();
118 :
119 415844 : for (const auto n : make_range(n_nodes))
120 : {
121 331340 : if (node_used[n])
122 123418 : continue;
123 :
124 207922 : const Real dist_sq = (q_points[qp] - secondary_elem.point(n)).norm_sq();
125 207922 : if (dist_sq < closest_dist_sq)
126 : {
127 94833 : second_closest_dist_sq = closest_dist_sq;
128 94833 : closest_dist_sq = dist_sq;
129 94833 : closest_node = n;
130 : }
131 113089 : else if (dist_sq < second_closest_dist_sq)
132 69180 : second_closest_dist_sq = dist_sq;
133 : }
134 :
135 84504 : if (closest_node == invalid_node || closest_dist_sq > matching_tol_sq)
136 0 : mooseError("Could not match nodal quadrature point ",
137 : qp,
138 : " at ",
139 0 : q_points[qp],
140 : " to a node on secondary mortar element ",
141 0 : secondary_elem.id(),
142 : " of type ",
143 0 : libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type()),
144 : ". The nearest unmatched node distance is ",
145 0 : std::sqrt(closest_dist_sq),
146 : ", which exceeds the tolerance ",
147 : matching_tol,
148 : ".");
149 :
150 84504 : if (second_closest_dist_sq <= matching_tol_sq)
151 0 : mooseError("Nodal quadrature point ",
152 : qp,
153 : " at ",
154 0 : q_points[qp],
155 : " does not map uniquely to secondary mortar element ",
156 0 : secondary_elem.id(),
157 : " of type ",
158 0 : libMesh::Utility::enum_to_string<ElemType>(secondary_elem.type()),
159 : ". Two unmatched nodes are within the matching tolerance ",
160 : matching_tol,
161 : ".");
162 :
163 84504 : qpoint_to_node[qp] = closest_node;
164 84504 : node_used[closest_node] = true;
165 : }
166 :
167 : #ifdef DEBUG
168 : // In optimized builds the mapping above skips already matched nodes for speed. In debug builds,
169 : // audit the full candidate set to catch ambiguous geometry or accidental many-to-one matches.
170 : std::vector<unsigned int> node_to_qpoint(n_nodes, invalid_node);
171 : for (const auto qp : make_range(q_points.size()))
172 : {
173 : const auto mapped_node = qpoint_to_node[qp];
174 : mooseAssert(mapped_node != invalid_node && mapped_node < n_nodes,
175 : "Invalid secondary node mapping for nodal quadrature point " << qp << ".");
176 : mooseAssert(node_to_qpoint[mapped_node] == invalid_node,
177 : "Secondary node " << mapped_node << " on mortar element " << secondary_elem.id()
178 : << " was matched to both nodal quadrature point "
179 : << node_to_qpoint[mapped_node] << " and " << qp << ".");
180 : node_to_qpoint[mapped_node] = qp;
181 :
182 : // Check the qp -> node direction without excluding nodes already matched by previous qps.
183 : unsigned int candidate_count = 0;
184 : unsigned int candidate_node = invalid_node;
185 : for (const auto n : make_range(n_nodes))
186 : if ((q_points[qp] - secondary_elem.point(n)).norm_sq() <= matching_tol_sq)
187 : {
188 : ++candidate_count;
189 : candidate_node = n;
190 : }
191 :
192 : mooseAssert(candidate_count == 1,
193 : "Nodal quadrature point " << qp << " on mortar element " << secondary_elem.id()
194 : << " has " << candidate_count
195 : << " secondary node candidates within tolerance "
196 : << matching_tol << ".");
197 : mooseAssert(candidate_node == mapped_node,
198 : "Nodal quadrature point " << qp << " on mortar element " << secondary_elem.id()
199 : << " was matched to node " << mapped_node
200 : << ", but the full candidate search found node "
201 : << candidate_node << ".");
202 : }
203 :
204 : for (const auto n : make_range(n_nodes))
205 : {
206 : mooseAssert(node_to_qpoint[n] != invalid_node,
207 : "Secondary node " << n << " on mortar element " << secondary_elem.id()
208 : << " was not matched to a nodal quadrature point.");
209 :
210 : // Check the node -> qp direction so every secondary node is also uniquely represented.
211 : unsigned int candidate_count = 0;
212 : unsigned int candidate_qp = invalid_node;
213 : for (const auto qp : make_range(q_points.size()))
214 : if ((q_points[qp] - secondary_elem.point(n)).norm_sq() <= matching_tol_sq)
215 : {
216 : ++candidate_count;
217 : candidate_qp = qp;
218 : }
219 :
220 : mooseAssert(candidate_count == 1,
221 : "Secondary node " << n << " on mortar element " << secondary_elem.id() << " has "
222 : << candidate_count
223 : << " nodal quadrature point candidates within tolerance "
224 : << matching_tol << ".");
225 : mooseAssert(candidate_qp == node_to_qpoint[n],
226 : "Secondary node " << n << " on mortar element " << secondary_elem.id()
227 : << " was matched to nodal quadrature point " << node_to_qpoint[n]
228 : << ", but the full candidate search found point " << candidate_qp
229 : << ".");
230 : }
231 : #endif
232 :
233 56802 : return qpoint_to_node;
234 28401 : }
235 : }
236 :
237 : class MortarNodalGeometryOutput : public Output
238 : {
239 : public:
240 114 : static InputParameters validParams()
241 : {
242 114 : auto params = Output::validParams();
243 228 : params.addPrivateParam<AutomaticMortarGeneration *>("_amg", nullptr);
244 114 : params.addPrivateParam<MooseApp *>(MooseBase::app_param, nullptr);
245 114 : params.set<std::string>(MooseBase::type_param) = "MortarNodalGeometryOutput";
246 114 : return params;
247 0 : };
248 :
249 114 : MortarNodalGeometryOutput(const InputParameters & params)
250 456 : : Output(params), _amg(*getCheckedPointerParam<AutomaticMortarGeneration *>("_amg"))
251 : {
252 114 : }
253 :
254 204 : void output() override
255 : {
256 : // Must call compute_nodal_geometry first!
257 408 : if (_amg._secondary_node_to_nodal_normal.empty() ||
258 204 : _amg._secondary_node_to_hh_nodal_tangents.empty())
259 0 : mooseError("No entries found in the secondary node -> nodal geometry map.");
260 :
261 204 : auto & problem = _app.feProblem();
262 204 : auto & subproblem = _amg._on_displaced
263 0 : ? static_cast<SubProblem &>(*problem.getDisplacedProblem())
264 204 : : static_cast<SubProblem &>(problem);
265 204 : auto & nodal_normals_es = subproblem.es();
266 :
267 204 : const std::string nodal_normals_sys_name = "nodal_normals";
268 :
269 204 : if (!_nodal_normals_system)
270 : {
271 306 : for (const auto s : make_range(nodal_normals_es.n_systems()))
272 204 : if (!nodal_normals_es.get_system(s).is_initialized())
273 : // This is really early on in the simulation and the systems have not been initialized. We
274 : // thus need to avoid calling reinit on systems that haven't even had their first init yet
275 0 : return;
276 :
277 102 : _nodal_normals_system =
278 102 : &nodal_normals_es.template add_system<ExplicitSystem>(nodal_normals_sys_name);
279 102 : _nnx_var_num = _nodal_normals_system->add_variable("nodal_normal_x", FEType(FIRST, LAGRANGE)),
280 102 : _nny_var_num = _nodal_normals_system->add_variable("nodal_normal_y", FEType(FIRST, LAGRANGE));
281 102 : _nnz_var_num = _nodal_normals_system->add_variable("nodal_normal_z", FEType(FIRST, LAGRANGE));
282 :
283 102 : _t1x_var_num =
284 102 : _nodal_normals_system->add_variable("nodal_tangent_1_x", FEType(FIRST, LAGRANGE)),
285 102 : _t1y_var_num =
286 102 : _nodal_normals_system->add_variable("nodal_tangent_1_y", FEType(FIRST, LAGRANGE));
287 102 : _t1z_var_num =
288 102 : _nodal_normals_system->add_variable("nodal_tangent_1_z", FEType(FIRST, LAGRANGE));
289 :
290 102 : _t2x_var_num =
291 102 : _nodal_normals_system->add_variable("nodal_tangent_2_x", FEType(FIRST, LAGRANGE)),
292 102 : _t2y_var_num =
293 102 : _nodal_normals_system->add_variable("nodal_tangent_2_y", FEType(FIRST, LAGRANGE));
294 102 : _t2z_var_num =
295 102 : _nodal_normals_system->add_variable("nodal_tangent_2_z", FEType(FIRST, LAGRANGE));
296 102 : nodal_normals_es.reinit();
297 : }
298 :
299 204 : const DofMap & dof_map = _nodal_normals_system->get_dof_map();
300 204 : std::vector<dof_id_type> dof_indices_nnx, dof_indices_nny, dof_indices_nnz;
301 204 : std::vector<dof_id_type> dof_indices_t1x, dof_indices_t1y, dof_indices_t1z;
302 204 : std::vector<dof_id_type> dof_indices_t2x, dof_indices_t2y, dof_indices_t2z;
303 :
304 204 : for (MeshBase::const_element_iterator el = _amg._mesh.elements_begin(),
305 204 : end_el = _amg._mesh.elements_end();
306 82399 : el != end_el;
307 82195 : ++el)
308 : {
309 82195 : const Elem * elem = *el;
310 :
311 : // Get the nodal dofs for this Elem.
312 82195 : dof_map.dof_indices(elem, dof_indices_nnx, _nnx_var_num);
313 82195 : dof_map.dof_indices(elem, dof_indices_nny, _nny_var_num);
314 82195 : dof_map.dof_indices(elem, dof_indices_nnz, _nnz_var_num);
315 :
316 82195 : dof_map.dof_indices(elem, dof_indices_t1x, _t1x_var_num);
317 82195 : dof_map.dof_indices(elem, dof_indices_t1y, _t1y_var_num);
318 82195 : dof_map.dof_indices(elem, dof_indices_t1z, _t1z_var_num);
319 :
320 82195 : dof_map.dof_indices(elem, dof_indices_t2x, _t2x_var_num);
321 82195 : dof_map.dof_indices(elem, dof_indices_t2y, _t2y_var_num);
322 82195 : dof_map.dof_indices(elem, dof_indices_t2z, _t2z_var_num);
323 :
324 : //
325 :
326 : // For each node of the Elem, if it is in the secondary_node_to_nodal_normal
327 : // container, set the corresponding nodal normal dof values.
328 599171 : for (MooseIndex(elem->n_vertices()) n = 0; n < elem->n_vertices(); ++n)
329 : {
330 516976 : auto it = _amg._secondary_node_to_nodal_normal.find(elem->node_ptr(n));
331 516976 : if (it != _amg._secondary_node_to_nodal_normal.end())
332 : {
333 37928 : _nodal_normals_system->solution->set(dof_indices_nnx[n], it->second(0));
334 37928 : _nodal_normals_system->solution->set(dof_indices_nny[n], it->second(1));
335 37928 : _nodal_normals_system->solution->set(dof_indices_nnz[n], it->second(2));
336 : }
337 :
338 516976 : auto it_tangent = _amg._secondary_node_to_hh_nodal_tangents.find(elem->node_ptr(n));
339 516976 : if (it_tangent != _amg._secondary_node_to_hh_nodal_tangents.end())
340 : {
341 37928 : _nodal_normals_system->solution->set(dof_indices_t1x[n], it_tangent->second[0](0));
342 37928 : _nodal_normals_system->solution->set(dof_indices_t1y[n], it_tangent->second[0](1));
343 37928 : _nodal_normals_system->solution->set(dof_indices_t1z[n], it_tangent->second[0](2));
344 :
345 37928 : _nodal_normals_system->solution->set(dof_indices_t2x[n], it_tangent->second[1](0));
346 37928 : _nodal_normals_system->solution->set(dof_indices_t2y[n], it_tangent->second[1](1));
347 37928 : _nodal_normals_system->solution->set(dof_indices_t2z[n], it_tangent->second[1](2));
348 : }
349 :
350 : } // end loop over nodes
351 204 : } // end loop over elems
352 :
353 : // Finish assembly.
354 204 : _nodal_normals_system->solution->close();
355 :
356 612 : std::set<std::string> sys_names = {nodal_normals_sys_name};
357 :
358 : // Write the nodal normals to file
359 204 : ExodusII_IO nodal_normals_writer(_amg._mesh);
360 :
361 : // Default to non-HDF5 output for wider compatibility
362 204 : nodal_normals_writer.set_hdf5_writing(false);
363 :
364 204 : nodal_normals_writer.write_equation_systems(
365 : "nodal_geometry_only.e", nodal_normals_es, &sys_names);
366 408 : }
367 :
368 : private:
369 : /// The mortar generation object that we will query for nodal normal and tangent information
370 : AutomaticMortarGeneration & _amg;
371 :
372 : ///@{
373 : /** Member variables for geometry debug output */
374 : libMesh::System * _nodal_normals_system = nullptr;
375 : unsigned int _nnx_var_num;
376 : unsigned int _nny_var_num;
377 : unsigned int _nnz_var_num;
378 :
379 : unsigned int _t1x_var_num;
380 : unsigned int _t1y_var_num;
381 : unsigned int _t1z_var_num;
382 :
383 : unsigned int _t2x_var_num;
384 : unsigned int _t2y_var_num;
385 : unsigned int _t2z_var_num;
386 : ///@}
387 : };
388 :
389 1079 : AutomaticMortarGeneration::AutomaticMortarGeneration(
390 : MooseApp & app,
391 : MeshBase & mesh_in,
392 : const std::pair<BoundaryID, BoundaryID> & boundary_key,
393 : const std::pair<SubdomainID, SubdomainID> & subdomain_key,
394 : bool on_displaced,
395 : bool periodic,
396 : const bool debug,
397 : const bool correct_edge_dropping,
398 : const Real minimum_projection_angle,
399 : const Mortar3DSubpatchPlane mortar_3d_subpatch_plane,
400 : const MortarSegmentTriangulationMode triangulation_mode,
401 : const bool triangulate_triangles,
402 1079 : const Mortar3DQuadraturePointMapping mortar_3d_qp_mapping)
403 : : ConsoleStreamInterface(app),
404 1079 : _app(app),
405 1079 : _mesh(mesh_in),
406 1079 : _debug(debug),
407 1079 : _on_displaced(on_displaced),
408 1079 : _periodic(periodic),
409 : // 3D mortar always builds the mortar segment mesh distributedly (each rank adds only its local
410 : // secondary elements). For 2D, we ghost the entire mortar interface when displaced, so
411 : // displaced meshes are always replicated; otherwise follow the parent mesh.
412 1079 : _distributed(_mesh.mesh_dimension() == 3 ? true : (!_on_displaced && !_mesh.is_replicated())),
413 1079 : _correct_edge_dropping(correct_edge_dropping),
414 1079 : _minimum_projection_angle(minimum_projection_angle),
415 1079 : _mortar_3d_subpatch_plane(mortar_3d_subpatch_plane),
416 1079 : _triangulation_mode(triangulation_mode),
417 1079 : _triangulate_triangles(triangulate_triangles),
418 2158 : _mortar_3d_qp_mapping(mortar_3d_qp_mapping)
419 : {
420 1079 : _primary_secondary_boundary_id_pairs.push_back(boundary_key);
421 1079 : _primary_requested_boundary_ids.insert(boundary_key.first);
422 1079 : _secondary_requested_boundary_ids.insert(boundary_key.second);
423 1079 : _primary_secondary_subdomain_id_pairs.push_back(subdomain_key);
424 1079 : _primary_boundary_subdomain_ids.insert(subdomain_key.first);
425 1079 : _secondary_boundary_subdomain_ids.insert(subdomain_key.second);
426 :
427 1079 : if (_distributed)
428 : _mortar_segment_mesh =
429 448 : std::make_unique<DistributedMesh>(_mesh.comm(), _mesh.spatial_dimension());
430 : else
431 : _mortar_segment_mesh =
432 631 : std::make_unique<ReplicatedMesh>(_mesh.comm(), _mesh.spatial_dimension());
433 1079 : }
434 :
435 : std::string
436 219 : AutomaticMortarGeneration::mortarInterfaceName() const
437 : {
438 219 : std::vector<std::string> string_vec(_primary_secondary_boundary_id_pairs.size() * 2 + 1);
439 438 : for (const auto i : index_range(_primary_secondary_boundary_id_pairs))
440 : {
441 219 : const auto [primary_bnd_id, secondary_bnd_id] = _primary_secondary_boundary_id_pairs[i];
442 219 : string_vec[2 * i] = std::to_string(primary_bnd_id);
443 219 : string_vec[2 * i + 1] = std::to_string(secondary_bnd_id);
444 : }
445 219 : string_vec.back() = _on_displaced ? "displaced" : "undisplaced";
446 438 : return MooseUtils::join(string_vec, "_");
447 219 : }
448 :
449 : void
450 1079 : AutomaticMortarGeneration::initOutput()
451 : {
452 1079 : if (!_debug)
453 965 : return;
454 :
455 114 : _output_params = std::make_unique<InputParameters>(MortarNodalGeometryOutput::validParams());
456 228 : _output_params->set<AutomaticMortarGeneration *>("_amg") = this;
457 228 : _output_params->set<FEProblemBase *>("_fe_problem_base") = &_app.feProblem();
458 114 : _output_params->set<MooseApp *>(MooseBase::app_param) = &_app;
459 114 : _output_params->set<std::string>(MooseBase::name_param) =
460 228 : "mortar_nodal_geometry_" + mortarInterfaceName();
461 228 : _output_params->finalize("MortarNodalGeometryOutput");
462 114 : _app.getOutputWarehouse().addOutput(std::make_shared<MortarNodalGeometryOutput>(*_output_params));
463 : }
464 :
465 : void
466 4624 : AutomaticMortarGeneration::clear()
467 : {
468 4624 : _msm_elem_to_reference_points.clear();
469 4624 : _mortar_segment_mesh->clear();
470 4624 : _nodes_to_secondary_elem_map.clear();
471 4624 : _nodes_to_primary_elem_map.clear();
472 4624 : _secondary_node_and_elem_to_xi2_primary_elem.clear();
473 4624 : _primary_node_and_elem_to_xi1_secondary_elem.clear();
474 4624 : _msm_elem_to_info.clear();
475 4624 : _lower_elem_to_side_id.clear();
476 4624 : _mortar_interface_coupling.clear();
477 4624 : _secondary_node_to_nodal_normal.clear();
478 4624 : _secondary_node_to_hh_nodal_tangents.clear();
479 4624 : _secondary_element_to_secondary_lowerd_element.clear();
480 4624 : _secondary_elems_to_mortar_segments.clear();
481 4624 : _secondary_ip_sub_ids.clear();
482 4624 : _primary_ip_sub_ids.clear();
483 4624 : _projected_secondary_nodes.clear();
484 4624 : _failed_secondary_node_projections.clear();
485 4624 : }
486 :
487 : const MortarSegmentReferencePoints &
488 44034 : AutomaticMortarGeneration::mortarSegmentReferencePoints(const Elem & mortar_segment_elem) const
489 : {
490 44034 : if (_mortar_3d_qp_mapping != Mortar3DQuadraturePointMapping::REFERENCE_INTERPOLATION)
491 0 : mooseError("Mortar segment reference points were requested for mortar segment element ",
492 0 : mortar_segment_elem.id(),
493 : ", but the reference-interpolation mapping mode is not enabled.");
494 :
495 44034 : const auto reference_points_it = _msm_elem_to_reference_points.find(&mortar_segment_elem);
496 44034 : if (reference_points_it == _msm_elem_to_reference_points.end())
497 0 : mooseError("No reference-point record was found for mortar segment element ",
498 0 : mortar_segment_elem.id(),
499 : ". The mortar segment info and reference-point maps are not aligned.");
500 :
501 88068 : return reference_points_it->second;
502 : }
503 :
504 : void
505 4621 : AutomaticMortarGeneration::buildNodeToElemMaps()
506 : {
507 4621 : if (_secondary_requested_boundary_ids.empty() || _primary_requested_boundary_ids.empty())
508 0 : mooseError(
509 : "Must specify secondary and primary boundary ids before building node-to-elem maps.");
510 :
511 : // Construct nodes_to_secondary_elem_map
512 4621 : for (const auto & secondary_elem :
513 953522 : as_range(_mesh.active_elements_begin(), _mesh.active_elements_end()))
514 : {
515 : // If this is not one of the lower-dimensional secondary side elements, go on to the next one.
516 472140 : if (!this->_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
517 443739 : continue;
518 :
519 112905 : for (const auto & nd : secondary_elem->node_ref_range())
520 : {
521 84504 : std::vector<const Elem *> & vec = _nodes_to_secondary_elem_map[nd.id()];
522 84504 : vec.push_back(secondary_elem);
523 : }
524 4621 : }
525 :
526 : // Construct nodes_to_primary_elem_map
527 4621 : for (const auto & primary_elem :
528 953522 : as_range(_mesh.active_elements_begin(), _mesh.active_elements_end()))
529 : {
530 : // If this is not one of the lower-dimensional primary side elements, go on to the next one.
531 472140 : if (!this->_primary_boundary_subdomain_ids.count(primary_elem->subdomain_id()))
532 439090 : continue;
533 :
534 149756 : for (const auto & nd : primary_elem->node_ref_range())
535 : {
536 116706 : std::vector<const Elem *> & vec = _nodes_to_primary_elem_map[nd.id()];
537 116706 : vec.push_back(primary_elem);
538 : }
539 4621 : }
540 4621 : }
541 :
542 : std::vector<Point>
543 599704 : AutomaticMortarGeneration::getNodalNormals(const Elem & secondary_elem) const
544 : {
545 599704 : std::vector<Point> nodal_normals(secondary_elem.n_nodes());
546 4141840 : for (const auto n : make_range(secondary_elem.n_nodes()))
547 3542136 : nodal_normals[n] = _secondary_node_to_nodal_normal.at(secondary_elem.node_ptr(n));
548 :
549 599704 : return nodal_normals;
550 0 : }
551 :
552 : const Elem *
553 0 : AutomaticMortarGeneration::getSecondaryLowerdElemFromSecondaryElem(
554 : dof_id_type secondary_elem_id) const
555 : {
556 : mooseAssert(_secondary_element_to_secondary_lowerd_element.count(secondary_elem_id),
557 : "Map should locate secondary element");
558 :
559 0 : return _secondary_element_to_secondary_lowerd_element.at(secondary_elem_id);
560 : }
561 :
562 : std::map<unsigned int, unsigned int>
563 24126 : AutomaticMortarGeneration::getSecondaryIpToLowerElementMap(const Elem & lower_secondary_elem) const
564 : {
565 24126 : std::map<unsigned int, unsigned int> secondary_ip_i_to_lower_secondary_i;
566 24126 : const Elem * const secondary_ip = lower_secondary_elem.interior_parent();
567 : mooseAssert(secondary_ip, "This should be non-null");
568 :
569 72378 : for (const auto i : make_range(lower_secondary_elem.n_nodes()))
570 : {
571 48252 : const auto & nd = lower_secondary_elem.node_ref(i);
572 48252 : secondary_ip_i_to_lower_secondary_i[secondary_ip->get_node_index(&nd)] = i;
573 : }
574 :
575 24126 : return secondary_ip_i_to_lower_secondary_i;
576 0 : }
577 :
578 : std::map<unsigned int, unsigned int>
579 24126 : AutomaticMortarGeneration::getPrimaryIpToLowerElementMap(
580 : const Elem & lower_primary_elem,
581 : const Elem & primary_elem,
582 : const Elem & /*lower_secondary_elem*/) const
583 : {
584 24126 : std::map<unsigned int, unsigned int> primary_ip_i_to_lower_primary_i;
585 :
586 72378 : for (const auto i : make_range(lower_primary_elem.n_nodes()))
587 : {
588 48252 : const auto & nd = lower_primary_elem.node_ref(i);
589 48252 : primary_ip_i_to_lower_primary_i[primary_elem.get_node_index(&nd)] = i;
590 : }
591 :
592 24126 : return primary_ip_i_to_lower_primary_i;
593 0 : }
594 :
595 : std::array<MooseUtils::SemidynamicVector<Point, 9>, 2>
596 0 : AutomaticMortarGeneration::getNodalTangents(const Elem & secondary_elem) const
597 : {
598 : // MetaPhysicL will check if we ran out of allocated space.
599 0 : MooseUtils::SemidynamicVector<Point, 9> nodal_tangents_one(0);
600 0 : MooseUtils::SemidynamicVector<Point, 9> nodal_tangents_two(0);
601 :
602 0 : for (const auto n : make_range(secondary_elem.n_nodes()))
603 : {
604 : const auto & tangent_vectors =
605 0 : libmesh_map_find(_secondary_node_to_hh_nodal_tangents, secondary_elem.node_ptr(n));
606 0 : nodal_tangents_one.push_back(tangent_vectors[0]);
607 0 : nodal_tangents_two.push_back(tangent_vectors[1]);
608 : }
609 :
610 0 : return {{nodal_tangents_one, nodal_tangents_two}};
611 : }
612 :
613 : std::vector<Point>
614 12323 : AutomaticMortarGeneration::getNormals(const Elem & secondary_elem,
615 : const std::vector<Real> & oned_xi1_pts) const
616 : {
617 12323 : std::vector<Point> xi1_pts(oned_xi1_pts.size());
618 24646 : for (const auto qp : index_range(oned_xi1_pts))
619 12323 : xi1_pts[qp] = oned_xi1_pts[qp];
620 :
621 24646 : return getNormals(secondary_elem, xi1_pts);
622 12323 : }
623 :
624 : std::vector<Point>
625 592359 : AutomaticMortarGeneration::getNormals(const Elem & secondary_elem,
626 : const std::vector<Point> & xi1_pts) const
627 : {
628 592359 : const auto mortar_dim = _mesh.mesh_dimension() - 1;
629 592359 : const auto num_qps = xi1_pts.size();
630 592359 : const auto nodal_normals = getNodalNormals(secondary_elem);
631 592359 : std::vector<Point> normals(num_qps);
632 :
633 4101907 : for (const auto n : make_range(secondary_elem.n_nodes()))
634 25004724 : for (const auto qp : make_range(num_qps))
635 : {
636 : const auto phi =
637 : (mortar_dim == 1)
638 21495176 : ? Moose::fe_lagrange_1D_shape(secondary_elem.default_order(), n, xi1_pts[qp](0))
639 21152054 : : Moose::fe_lagrange_2D_shape(secondary_elem.type(),
640 21152054 : secondary_elem.default_order(),
641 : n,
642 21152054 : static_cast<const TypeVector<Real> &>(xi1_pts[qp]));
643 21495176 : normals[qp] += phi * nodal_normals[n];
644 : }
645 :
646 592359 : if (_periodic)
647 64266 : for (auto & normal : normals)
648 50605 : normal *= -1;
649 :
650 1184718 : return normals;
651 592359 : }
652 :
653 : void
654 4263 : AutomaticMortarGeneration::buildMortarSegmentMesh()
655 : {
656 : using std::abs;
657 :
658 4263 : dof_id_type local_id_index = 0;
659 4263 : std::size_t node_unique_id_offset = 0;
660 :
661 : // Create an offset by the maximum number of mortar segment elements that can be created *plus*
662 : // the number of lower-dimensional secondary subdomain elements. Recall that the number of mortar
663 : // segments created is a function of node projection, *and* that if we split elems we will delete
664 : // that elem which has already taken a unique id
665 8526 : for (const auto & pr : _primary_secondary_boundary_id_pairs)
666 : {
667 4263 : const auto primary_bnd_id = pr.first;
668 4263 : const auto secondary_bnd_id = pr.second;
669 : const auto num_primary_nodes =
670 8526 : std::distance(_mesh.bid_nodes_begin(primary_bnd_id), _mesh.bid_nodes_end(primary_bnd_id));
671 8526 : const auto num_secondary_nodes = std::distance(_mesh.bid_nodes_begin(secondary_bnd_id),
672 8526 : _mesh.bid_nodes_end(secondary_bnd_id));
673 : mooseAssert(num_primary_nodes,
674 : "There are no primary nodes on boundary ID "
675 : << primary_bnd_id << ". Does that bondary ID even exist on the mesh?");
676 : mooseAssert(num_secondary_nodes,
677 : "There are no secondary nodes on boundary ID "
678 : << secondary_bnd_id << ". Does that bondary ID even exist on the mesh?");
679 :
680 4263 : node_unique_id_offset += num_primary_nodes + 2 * num_secondary_nodes;
681 : }
682 :
683 : // 1.) Add all lower-dimensional secondary side elements as the "initial" mortar segments.
684 4263 : for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
685 4263 : end_el = _mesh.active_elements_end();
686 323189 : el != end_el;
687 318926 : ++el)
688 : {
689 318926 : const Elem * secondary_elem = *el;
690 :
691 : // If this is not one of the lower-dimensional secondary side elements, go on to the next one.
692 318926 : if (!this->_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
693 299218 : continue;
694 :
695 19708 : std::vector<Node *> new_nodes;
696 61134 : for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
697 : {
698 41426 : new_nodes.push_back(_mortar_segment_mesh->add_point(
699 : secondary_elem->point(n), secondary_elem->node_id(n), secondary_elem->processor_id()));
700 41426 : Node * const new_node = new_nodes.back();
701 41426 : new_node->set_unique_id(new_node->id() + node_unique_id_offset);
702 : }
703 :
704 19708 : std::unique_ptr<Elem> new_elem;
705 19708 : if (secondary_elem->default_order() == SECOND)
706 2010 : new_elem = std::make_unique<Edge3>();
707 : else
708 17698 : new_elem = std::make_unique<Edge2>();
709 :
710 19708 : new_elem->processor_id() = secondary_elem->processor_id();
711 19708 : new_elem->subdomain_id() = secondary_elem->subdomain_id();
712 19708 : new_elem->set_id(local_id_index++);
713 19708 : new_elem->set_unique_id(new_elem->id());
714 :
715 61134 : for (MooseIndex(new_elem->n_nodes()) n = 0; n < new_elem->n_nodes(); ++n)
716 41426 : new_elem->set_node(n, new_nodes[n]);
717 :
718 19708 : Elem * new_elem_ptr = _mortar_segment_mesh->add_elem(new_elem.release());
719 :
720 : // The xi^(1) values for this mortar segment are initially -1 and 1.
721 19708 : MortarSegmentInfo msinfo;
722 19708 : msinfo.xi1_a = -1;
723 19708 : msinfo.xi1_b = +1;
724 19708 : msinfo.secondary_elem = secondary_elem;
725 :
726 19708 : auto new_container_it0 = _secondary_node_and_elem_to_xi2_primary_elem.find(
727 19708 : std::make_pair(secondary_elem->node_ptr(0), secondary_elem)),
728 19708 : new_container_it1 = _secondary_node_and_elem_to_xi2_primary_elem.find(
729 19708 : std::make_pair(secondary_elem->node_ptr(1), secondary_elem));
730 :
731 : bool new_container_node0_found =
732 19708 : (new_container_it0 != _secondary_node_and_elem_to_xi2_primary_elem.end()),
733 : new_container_node1_found =
734 19708 : (new_container_it1 != _secondary_node_and_elem_to_xi2_primary_elem.end());
735 :
736 19708 : const Elem * node0_primary_candidate = nullptr;
737 19708 : const Elem * node1_primary_candidate = nullptr;
738 :
739 19708 : if (new_container_node0_found)
740 : {
741 16379 : const auto & xi2_primary_elem_pair = new_container_it0->second;
742 16379 : msinfo.xi2_a = xi2_primary_elem_pair.first;
743 16379 : node0_primary_candidate = xi2_primary_elem_pair.second;
744 : }
745 :
746 19708 : if (new_container_node1_found)
747 : {
748 19370 : const auto & xi2_primary_elem_pair = new_container_it1->second;
749 19370 : msinfo.xi2_b = xi2_primary_elem_pair.first;
750 19370 : node1_primary_candidate = xi2_primary_elem_pair.second;
751 : }
752 :
753 : // If both node0 and node1 agree on the primary element they are
754 : // projected into, then this mortar segment fits entirely within
755 : // a single primary element, and we can go ahead and set the
756 : // msinfo.primary_elem pointer now.
757 19708 : if (node0_primary_candidate == node1_primary_candidate)
758 7417 : msinfo.primary_elem = node0_primary_candidate;
759 :
760 : // Associate this MSM elem with the MortarSegmentInfo.
761 19708 : _msm_elem_to_info.emplace(new_elem_ptr, msinfo);
762 :
763 : // Maintain the mapping between secondary elems and mortar segment elems contained within them.
764 : // Initially, only the original secondary_elem is present.
765 19708 : _secondary_elems_to_mortar_segments[secondary_elem->id()].insert(new_elem_ptr);
766 23971 : }
767 :
768 : // 2.) Insert new nodes from primary side and split mortar segments as necessary.
769 24187 : for (const auto & pr : _primary_node_and_elem_to_xi1_secondary_elem)
770 : {
771 19924 : auto key = pr.first;
772 19924 : auto val = pr.second;
773 :
774 19924 : const Node * primary_node = std::get<1>(key);
775 19924 : Real xi1 = val.first;
776 19924 : const Elem * secondary_elem = val.second;
777 :
778 : // If this is an aligned node, we don't need to do anything.
779 19924 : if (abs(abs(xi1) - 1.) < _xi_tolerance)
780 7601 : continue;
781 :
782 12323 : auto && order = secondary_elem->default_order();
783 :
784 : // Determine physical location of new point to be inserted.
785 12323 : Point new_pt(0);
786 37501 : for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
787 25178 : new_pt += Moose::fe_lagrange_1D_shape(order, n, xi1) * secondary_elem->point(n);
788 :
789 : // Find the current mortar segment that will have to be split.
790 12323 : auto & mortar_segment_set = _secondary_elems_to_mortar_segments[secondary_elem->id()];
791 12323 : Elem * current_mortar_segment = nullptr;
792 12323 : MortarSegmentInfo * info = nullptr;
793 :
794 12323 : for (const auto & mortar_segment_candidate : mortar_segment_set)
795 : {
796 : try
797 : {
798 12323 : info = &_msm_elem_to_info.at(mortar_segment_candidate);
799 : }
800 0 : catch (std::out_of_range &)
801 : {
802 0 : mooseError("MortarSegmentInfo not found for the mortar segment candidate");
803 0 : }
804 12323 : if (info->xi1_a <= xi1 && xi1 <= info->xi1_b)
805 : {
806 12323 : current_mortar_segment = mortar_segment_candidate;
807 12323 : break;
808 : }
809 : }
810 :
811 : // Make sure we found one.
812 12323 : if (current_mortar_segment == nullptr)
813 0 : mooseError("Unable to find appropriate mortar segment during linear search!");
814 :
815 : // If node lands on endpoint of segment, don't split.
816 : // Jacob: This condition was getting missed by the < comparison a few lines above. To fix it I
817 : // just made it <= and put this condition in to handle equality different. It probably could be
818 : // done with a tolerance but the the toleranced equality is already handled later when we drop
819 : // segments with small volume.
820 12323 : if (info->xi1_a == xi1 || xi1 == info->xi1_b)
821 0 : continue;
822 :
823 12323 : const auto new_id = _mortar_segment_mesh->max_node_id();
824 : mooseAssert(_mortar_segment_mesh->comm().verify(new_id),
825 : "new_id must be the same on all processes");
826 : Node * const new_node =
827 12323 : _mortar_segment_mesh->add_point(new_pt, new_id, secondary_elem->processor_id());
828 12323 : new_node->set_unique_id(new_id + node_unique_id_offset);
829 :
830 : // Reconstruct the nodal normal at xi1. This will help us
831 : // determine the orientation of the primary elems relative to the
832 : // new mortar segments.
833 12323 : const Point normal = getNormals(*secondary_elem, std::vector<Real>({xi1}))[0];
834 :
835 : // Get the set of primary_node neighbors.
836 12323 : if (this->_nodes_to_primary_elem_map.find(primary_node->id()) ==
837 24646 : this->_nodes_to_primary_elem_map.end())
838 0 : mooseError("We should already have built this primary node to elem pair!");
839 : const std::vector<const Elem *> & primary_node_neighbors =
840 12323 : this->_nodes_to_primary_elem_map[primary_node->id()];
841 :
842 : // Sanity check
843 12323 : if (primary_node_neighbors.size() == 0 || primary_node_neighbors.size() > 2)
844 0 : mooseError("We must have either 1 or 2 primary side nodal neighbors, but we had ",
845 0 : primary_node_neighbors.size());
846 :
847 : // Primary Elem pointers which we will eventually assign to the
848 : // mortar segments being created. We start by assuming
849 : // primary_node_neighbor[0] is on the "left" and
850 : // primary_node_neighbor[1]/"nothing" is on the "right" and then
851 : // swap them if that's not the case.
852 12323 : const Elem * left_primary_elem = primary_node_neighbors[0];
853 : const Elem * right_primary_elem =
854 12323 : (primary_node_neighbors.size() == 2) ? primary_node_neighbors[1] : nullptr;
855 :
856 12323 : Real left_xi2 = MortarSegmentInfo::invalid_xi, right_xi2 = MortarSegmentInfo::invalid_xi;
857 :
858 : // Storage for z-component of cross products for determining
859 : // orientation.
860 : std::array<Real, 2> secondary_node_cps;
861 12323 : std::vector<Real> primary_node_cps(primary_node_neighbors.size());
862 :
863 : // Store z-component of left and right secondary node cross products with the nodal normal.
864 36969 : for (unsigned int nid = 0; nid < 2; ++nid)
865 24646 : secondary_node_cps[nid] = normal.cross(secondary_elem->point(nid) - new_pt)(2);
866 :
867 33978 : for (MooseIndex(primary_node_neighbors) mnn = 0; mnn < primary_node_neighbors.size(); ++mnn)
868 : {
869 21655 : const Elem * primary_neigh = primary_node_neighbors[mnn];
870 21655 : Point opposite = (primary_neigh->node_ptr(0) == primary_node) ? primary_neigh->point(1)
871 12323 : : primary_neigh->point(0);
872 21655 : Point cp = normal.cross(opposite - new_pt);
873 21655 : primary_node_cps[mnn] = cp(2);
874 : }
875 :
876 : // We will verify that only 1 orientation is actually valid.
877 12323 : bool orientation1_valid = false, orientation2_valid = false;
878 :
879 12323 : if (primary_node_neighbors.size() == 2)
880 : {
881 : // 2 primary neighbor case
882 9401 : orientation1_valid = (secondary_node_cps[0] * primary_node_cps[0] > 0.) &&
883 69 : (secondary_node_cps[1] * primary_node_cps[1] > 0.);
884 :
885 18595 : orientation2_valid = (secondary_node_cps[0] * primary_node_cps[1] > 0.) &&
886 9263 : (secondary_node_cps[1] * primary_node_cps[0] > 0.);
887 : }
888 2991 : else if (primary_node_neighbors.size() == 1)
889 : {
890 : // 1 primary neighbor case
891 2991 : orientation1_valid = (secondary_node_cps[0] * primary_node_cps[0] > 0.);
892 2991 : orientation2_valid = (secondary_node_cps[1] * primary_node_cps[0] > 0.);
893 : }
894 : else
895 0 : mooseError("Invalid primary node neighbors size ", primary_node_neighbors.size());
896 :
897 : // Verify that both orientations are not simultaneously valid/invalid. If they are not, then we
898 : // are going to throw an exception instead of erroring out since we can easily reach this point
899 : // if we have one bad linear solve. It's better in general to catch the error and then try a
900 : // smaller time-step
901 12323 : if (orientation1_valid && orientation2_valid)
902 : throw MooseException(
903 0 : "AutomaticMortarGeneration: Both orientations cannot simultaneously be valid.");
904 :
905 : // We are going to treat the case where both orientations are invalid as a case in which we
906 : // should not be splitting the mortar mesh to incorporate primary mesh elements.
907 : // In practice, this case has appeared for very oblique projections, so we assume these cases
908 : // will not be considered in mortar thermomechanical contact.
909 12323 : if (!orientation1_valid && !orientation2_valid)
910 : {
911 0 : mooseDoOnce(mooseWarning(
912 : "AutomaticMortarGeneration: Unable to determine valid secondary-primary orientation. "
913 : "Consequently we will consider projection of the primary node invalid and not split the "
914 : "mortar segment. "
915 : "This situation can indicate there are very oblique projections between primary (mortar) "
916 : "and secondary (non-mortar) surfaces for a good problem set up. It can also mean your "
917 : "time step is too large. This message is only printed once."));
918 0 : continue;
919 0 : }
920 :
921 : // Make an Elem on the left
922 12323 : std::unique_ptr<Elem> new_elem_left;
923 12323 : if (order == SECOND)
924 532 : new_elem_left = std::make_unique<Edge3>();
925 : else
926 11791 : new_elem_left = std::make_unique<Edge2>();
927 :
928 12323 : new_elem_left->processor_id() = current_mortar_segment->processor_id();
929 12323 : new_elem_left->subdomain_id() = current_mortar_segment->subdomain_id();
930 12323 : new_elem_left->set_id(local_id_index++);
931 12323 : new_elem_left->set_unique_id(new_elem_left->id());
932 12323 : new_elem_left->set_node(0, current_mortar_segment->node_ptr(0));
933 12323 : new_elem_left->set_node(1, new_node);
934 :
935 : // Make an Elem on the right
936 12323 : std::unique_ptr<Elem> new_elem_right;
937 12323 : if (order == SECOND)
938 532 : new_elem_right = std::make_unique<Edge3>();
939 : else
940 11791 : new_elem_right = std::make_unique<Edge2>();
941 :
942 12323 : new_elem_right->processor_id() = current_mortar_segment->processor_id();
943 12323 : new_elem_right->subdomain_id() = current_mortar_segment->subdomain_id();
944 12323 : new_elem_right->set_id(local_id_index++);
945 12323 : new_elem_right->set_unique_id(new_elem_right->id());
946 12323 : new_elem_right->set_node(0, new_node);
947 12323 : new_elem_right->set_node(1, current_mortar_segment->node_ptr(1));
948 :
949 12323 : if (order == SECOND)
950 : {
951 : // left
952 532 : Point left_interior_point(0);
953 532 : Real left_interior_xi = (xi1 + info->xi1_a) / 2;
954 :
955 : // This is eta for the current mortar segment that we're splitting
956 532 : Real current_left_interior_eta =
957 532 : (2. * left_interior_xi - info->xi1_a - info->xi1_b) / (info->xi1_b - info->xi1_a);
958 :
959 532 : for (MooseIndex(current_mortar_segment->n_nodes()) n = 0;
960 2128 : n < current_mortar_segment->n_nodes();
961 : ++n)
962 1596 : left_interior_point += Moose::fe_lagrange_1D_shape(order, n, current_left_interior_eta) *
963 1596 : current_mortar_segment->point(n);
964 :
965 532 : const auto new_interior_left_id = _mortar_segment_mesh->max_node_id();
966 : mooseAssert(_mortar_segment_mesh->comm().verify(new_interior_left_id),
967 : "new_id must be the same on all processes");
968 532 : Node * const new_interior_node_left = _mortar_segment_mesh->add_point(
969 532 : left_interior_point, new_interior_left_id, new_elem_left->processor_id());
970 532 : new_elem_left->set_node(2, new_interior_node_left);
971 532 : new_interior_node_left->set_unique_id(new_interior_left_id + node_unique_id_offset);
972 :
973 : // right
974 532 : Point right_interior_point(0);
975 532 : Real right_interior_xi = (xi1 + info->xi1_b) / 2;
976 : // This is eta for the current mortar segment that we're splitting
977 532 : Real current_right_interior_eta =
978 532 : (2. * right_interior_xi - info->xi1_a - info->xi1_b) / (info->xi1_b - info->xi1_a);
979 :
980 532 : for (MooseIndex(current_mortar_segment->n_nodes()) n = 0;
981 2128 : n < current_mortar_segment->n_nodes();
982 : ++n)
983 1596 : right_interior_point += Moose::fe_lagrange_1D_shape(order, n, current_right_interior_eta) *
984 1596 : current_mortar_segment->point(n);
985 :
986 532 : const auto new_interior_id_right = _mortar_segment_mesh->max_node_id();
987 : mooseAssert(_mortar_segment_mesh->comm().verify(new_interior_id_right),
988 : "new_id must be the same on all processes");
989 532 : Node * const new_interior_node_right = _mortar_segment_mesh->add_point(
990 532 : right_interior_point, new_interior_id_right, new_elem_right->processor_id());
991 532 : new_elem_right->set_node(2, new_interior_node_right);
992 532 : new_interior_node_right->set_unique_id(new_interior_id_right + node_unique_id_offset);
993 : }
994 :
995 : // If orientation 2 was valid, swap the left and right primaries.
996 12323 : if (orientation2_valid)
997 12254 : std::swap(left_primary_elem, right_primary_elem);
998 :
999 : // Now that we know left_primary_elem and right_primary_elem, we can determine left_xi2 and
1000 : // right_xi2.
1001 12323 : if (left_primary_elem)
1002 9332 : left_xi2 = (primary_node == left_primary_elem->node_ptr(0)) ? -1 : +1;
1003 12323 : if (right_primary_elem)
1004 12323 : right_xi2 = (primary_node == right_primary_elem->node_ptr(0)) ? -1 : +1;
1005 :
1006 : // Grab the MortarSegmentInfo object associated with this
1007 : // segment. We can use "at()" here since we want this to fail if
1008 : // current_mortar_segment is not found... Since we're going to
1009 : // erase this entry from the map momentarily, we make an actual
1010 : // copy rather than grabbing a reference.
1011 12323 : auto msm_it = _msm_elem_to_info.find(current_mortar_segment);
1012 12323 : if (msm_it == _msm_elem_to_info.end())
1013 0 : mooseError("MortarSegmentInfo not found for current_mortar_segment.");
1014 12323 : MortarSegmentInfo current_msinfo = msm_it->second;
1015 :
1016 : // add_left
1017 : {
1018 12323 : Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem_left.release());
1019 :
1020 : // Create new MortarSegmentInfo objects for new_elem_left
1021 12323 : MortarSegmentInfo new_msinfo_left;
1022 :
1023 : // The new MortarSegmentInfo info objects inherit their "outer"
1024 : // information from current_msinfo and the rest is determined by
1025 : // the Node being inserted.
1026 12323 : new_msinfo_left.xi1_a = current_msinfo.xi1_a;
1027 12323 : new_msinfo_left.xi2_a = current_msinfo.xi2_a;
1028 12323 : new_msinfo_left.secondary_elem = secondary_elem;
1029 12323 : new_msinfo_left.xi1_b = xi1;
1030 12323 : new_msinfo_left.xi2_b = left_xi2;
1031 12323 : new_msinfo_left.primary_elem = left_primary_elem;
1032 :
1033 : // Add new msinfo objects to the map.
1034 12323 : _msm_elem_to_info.emplace(msm_new_elem, new_msinfo_left);
1035 :
1036 : // We need to insert new_elem_left in
1037 : // the mortar_segment_set for this secondary_elem.
1038 12323 : mortar_segment_set.insert(msm_new_elem);
1039 : }
1040 :
1041 : // add_right
1042 : {
1043 12323 : Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem_right.release());
1044 :
1045 : // Create new MortarSegmentInfo objects for new_elem_right
1046 12323 : MortarSegmentInfo new_msinfo_right;
1047 :
1048 12323 : new_msinfo_right.xi1_b = current_msinfo.xi1_b;
1049 12323 : new_msinfo_right.xi2_b = current_msinfo.xi2_b;
1050 12323 : new_msinfo_right.secondary_elem = secondary_elem;
1051 12323 : new_msinfo_right.xi1_a = xi1;
1052 12323 : new_msinfo_right.xi2_a = right_xi2;
1053 12323 : new_msinfo_right.primary_elem = right_primary_elem;
1054 :
1055 12323 : _msm_elem_to_info.emplace(msm_new_elem, new_msinfo_right);
1056 :
1057 12323 : mortar_segment_set.insert(msm_new_elem);
1058 : }
1059 :
1060 : // Erase the MortarSegmentInfo object for current_mortar_segment from the map.
1061 12323 : _msm_elem_to_info.erase(msm_it);
1062 :
1063 : // current_mortar_segment must be erased from the
1064 : // mortar_segment_set since it has now been split.
1065 12323 : mortar_segment_set.erase(current_mortar_segment);
1066 :
1067 : // The original mortar segment has been split, so erase it from
1068 : // the mortar segment mesh.
1069 12323 : _mortar_segment_mesh->delete_elem(current_mortar_segment);
1070 12323 : }
1071 :
1072 : // Remove all MSM elements without a primary contribution
1073 : /**
1074 : * This was a change to how inactive LM DoFs are handled. Now mortar segment elements
1075 : * are not used in assembly if there is no corresponding primary element and inactive
1076 : * LM DoFs (those with no contribution to an active primary element) are zeroed.
1077 : */
1078 36294 : for (auto msm_elem : _mortar_segment_mesh->active_element_ptr_range())
1079 : {
1080 32031 : MortarSegmentInfo & msinfo = libmesh_map_find(_msm_elem_to_info, msm_elem);
1081 32031 : Elem * primary_elem = const_cast<Elem *>(msinfo.primary_elem);
1082 60733 : if (primary_elem == nullptr || abs(msinfo.xi2_a) > 1.0 + TOLERANCE ||
1083 28702 : abs(msinfo.xi2_b) > 1.0 + TOLERANCE)
1084 : {
1085 : // Erase from secondary to msms map
1086 3329 : auto it = _secondary_elems_to_mortar_segments.find(msinfo.secondary_elem->id());
1087 : mooseAssert(it != _secondary_elems_to_mortar_segments.end(),
1088 : "We should have found the element");
1089 3329 : auto & msm_set = it->second;
1090 3329 : msm_set.erase(msm_elem);
1091 : // We may be creating nodes with only one element neighbor where before this removal there
1092 : // were two. But the nodal normal used in computations will reflect the two-neighbor geometry.
1093 : // For a lower-d secondary mesh corner, that will imply the corner node will have a tilted
1094 : // normal vector (same for tangents) despite the mortar segment mesh not including its
1095 : // vertical neighboring element. It is the secondary element neighbors (not mortar segment
1096 : // mesh neighbors) that determine the nodal normal field.
1097 3329 : if (msm_set.empty())
1098 338 : _secondary_elems_to_mortar_segments.erase(it);
1099 :
1100 : // Erase msinfo
1101 3329 : _msm_elem_to_info.erase(msm_elem);
1102 :
1103 : // Remove element from mortar segment mesh
1104 3329 : _mortar_segment_mesh->delete_elem(msm_elem);
1105 : }
1106 : else
1107 : {
1108 28702 : _secondary_ip_sub_ids.insert(msinfo.secondary_elem->interior_parent()->subdomain_id());
1109 28702 : _primary_ip_sub_ids.insert(msinfo.primary_elem->interior_parent()->subdomain_id());
1110 : }
1111 4263 : }
1112 :
1113 4263 : std::unordered_set<Node *> msm_connected_nodes;
1114 :
1115 : // Deleting elements may produce isolated nodes.
1116 : // Loops for identifying and removing such nodes from mortar segment mesh.
1117 32965 : for (const auto & element : _mortar_segment_mesh->element_ptr_range())
1118 88648 : for (auto & n : element->node_ref_range())
1119 64209 : msm_connected_nodes.insert(&n);
1120 :
1121 43643 : for (const auto & node : _mortar_segment_mesh->node_ptr_range())
1122 39380 : if (!msm_connected_nodes.count(node))
1123 8124 : _mortar_segment_mesh->delete_node(node);
1124 :
1125 : #ifdef DEBUG
1126 : // Verify that all segments without primary contribution have been deleted
1127 : for (auto msm_elem : _mortar_segment_mesh->active_element_ptr_range())
1128 : {
1129 : const MortarSegmentInfo & msinfo = libmesh_map_find(_msm_elem_to_info, msm_elem);
1130 : mooseAssert(msinfo.primary_elem != nullptr,
1131 : "All mortar segment elements should have valid "
1132 : "primary element.");
1133 : }
1134 : #endif
1135 :
1136 4263 : _mortar_segment_mesh->cache_elem_data();
1137 :
1138 : // (Optionally) Write the mortar segment mesh to file for inspection
1139 4263 : if (_debug)
1140 12 : outputMortarMesh();
1141 :
1142 4263 : buildCouplingInformation();
1143 4263 : }
1144 :
1145 : void
1146 105 : AutomaticMortarGeneration::outputMortarMesh()
1147 : {
1148 105 : ExodusII_IO mortar_segment_mesh_writer(*_mortar_segment_mesh);
1149 :
1150 : // Default to non-HDF5 output for wider compatibility
1151 105 : mortar_segment_mesh_writer.set_hdf5_writing(false);
1152 :
1153 : std::array<std::string, 3> file_pieces = {
1154 105 : _app.getOutputFileBase(/*for_non_moose_build_output=*/true),
1155 : mortarInterfaceName(),
1156 210 : "mortar_segment_mesh.e"};
1157 105 : mortar_segment_mesh_writer.write(MooseUtils::join(file_pieces, "_"));
1158 105 : }
1159 :
1160 : void
1161 358 : AutomaticMortarGeneration::buildMortarSegmentMesh3d()
1162 : {
1163 358 : const bool use_reference_interpolation =
1164 358 : _mortar_3d_qp_mapping == Mortar3DQuadraturePointMapping::REFERENCE_INTERPOLATION;
1165 :
1166 : // Add an integer flag to mortar segment mesh to keep track of which subelem
1167 : // of second order primal elements mortar segments correspond to
1168 716 : auto secondary_sub_elem = _mortar_segment_mesh->add_elem_integer("secondary_sub_elem");
1169 716 : auto primary_sub_elem = _mortar_segment_mesh->add_elem_integer("primary_sub_elem");
1170 :
1171 : // Assign globally unique node/element IDs via an exclusive prefix scan: each rank's bound is
1172 : // local_secondary_sub_elems * visible_primary_sub_elems * 9, where 9 is the maximum nodes a
1173 : // single secondary/primary sub-element pair can produce (8-vertex clipped polygon + center).
1174 : // The result is cached and invalidated by meshChanged(), so the allgather only runs on topology
1175 : // changes, not on every displaced-mesh residual update.
1176 358 : if (!_msm_node_id_start.has_value())
1177 : {
1178 358 : dof_id_type local_secondary_sub_elems = 0, visible_primary_sub_elems = 0;
1179 716 : for (const auto & [primary_sub_id, secondary_sub_id] : _primary_secondary_subdomain_id_pairs)
1180 : {
1181 358 : for (const auto * const el :
1182 6868 : _mesh.active_local_subdomain_elements_ptr_range(secondary_sub_id))
1183 6510 : local_secondary_sub_elems += el->n_sub_elem();
1184 16890 : for (const auto * const el : _mesh.active_subdomain_elements_ptr_range(primary_sub_id))
1185 16890 : visible_primary_sub_elems += el->n_sub_elem();
1186 : }
1187 358 : const dof_id_type per_rank_bound = local_secondary_sub_elems * visible_primary_sub_elems * 9;
1188 358 : std::vector<dof_id_type> per_rank_bounds;
1189 358 : _mesh.comm().allgather(per_rank_bound, per_rank_bounds);
1190 358 : dof_id_type start = 0;
1191 468 : for (const auto r : make_range(_mesh.processor_id()))
1192 110 : start += per_rank_bounds[r];
1193 358 : _msm_node_id_start = start;
1194 358 : }
1195 358 : dof_id_type next_node_id = *_msm_node_id_start;
1196 : // Element IDs use the same starting offset: node and element IDs are separately numbered, and
1197 : // element count per clip (n triangles) is always <= node count (n+1), so per_rank_bound covers
1198 : // both.
1199 358 : dof_id_type next_elem_id = next_node_id;
1200 :
1201 : // Loop through mortar secondary and primary pairs to create mortar segment mesh between each
1202 713 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1203 : {
1204 358 : const auto primary_subd_id = pr.first;
1205 358 : const auto secondary_subd_id = pr.second;
1206 :
1207 : // Build k-d tree for use in Step 1.2 for primary interface coarse screening
1208 358 : NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, primary_subd_id);
1209 : subdomain_kd_tree_t kd_tree(
1210 358 : 3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
1211 :
1212 : // Construct the KD tree.
1213 358 : kd_tree.buildIndex();
1214 :
1215 : // Return the unoriented geometric normal of a linearized subpatch. These expressions are the
1216 : // TRI3 and QUAD4 mapping tangents evaluated at the reference center, equivalent to evaluating
1217 : // the first-order finite-element normal there without constructing a temporary element.
1218 131725 : auto get_sub_elem_geometric_normal = [](const std::vector<Point> & nodes)
1219 : {
1220 131725 : Point dxdxi;
1221 131725 : Point dxdeta;
1222 131725 : if (nodes.size() == 3)
1223 : {
1224 64056 : dxdxi = nodes[1] - nodes[0];
1225 64056 : dxdeta = nodes[2] - nodes[0];
1226 : }
1227 67669 : else if (nodes.size() == 4)
1228 : {
1229 : // Bilinear center tangents define one normal for the full quad instead of selecting one of
1230 : // the two diagonal triangle normals.
1231 67669 : dxdxi = 0.25 * (nodes[1] + nodes[2] - nodes[0] - nodes[3]);
1232 67669 : dxdeta = 0.25 * (nodes[2] + nodes[3] - nodes[0] - nodes[1]);
1233 : }
1234 : else
1235 0 : mooseError("GEOMETRIC_NORMAL 3D mortar subpatch plane construction only supports "
1236 : "triangular and quadrilateral subpatches, but received ",
1237 0 : nodes.size(),
1238 : " nodes.");
1239 :
1240 131725 : Point geometric_normal = dxdxi.cross(dxdeta);
1241 131725 : const auto normal_norm = geometric_normal.norm();
1242 : // The cross product has units of area, so compare it with the product of tangent lengths.
1243 : // Their ratio is the sine of the included angle and is independent of the mesh length scale.
1244 131725 : if (normal_norm <= TOLERANCE * dxdxi.norm() * dxdeta.norm())
1245 0 : mooseError("GEOMETRIC_NORMAL 3D mortar subpatch plane construction encountered a "
1246 : "degenerate subpatch.");
1247 :
1248 131725 : geometric_normal /= normal_norm;
1249 263450 : return geometric_normal;
1250 : };
1251 :
1252 : /**
1253 : * Step 1: Build mortar segments for all secondary elements
1254 : */
1255 358 : for (MeshBase::const_element_iterator el = _mesh.active_local_elements_begin(),
1256 358 : end_el = _mesh.active_local_elements_end();
1257 108596 : el != end_el;
1258 108238 : ++el)
1259 : {
1260 108241 : const Elem * secondary_side_elem = *el;
1261 :
1262 108241 : const Real secondary_volume = secondary_side_elem->volume();
1263 :
1264 : // If this Elem is not in the current secondary subdomain, go on to the next one.
1265 108241 : if (secondary_side_elem->subdomain_id() != secondary_subd_id)
1266 102092 : continue;
1267 :
1268 6149 : auto [secondary_elem_to_msm_map_it, insertion_happened] =
1269 6149 : _secondary_elems_to_mortar_segments.emplace(secondary_side_elem->id(),
1270 12298 : std::set<Elem *, CompareDofObjectsByID>{});
1271 6149 : libmesh_ignore(insertion_happened);
1272 6149 : auto & secondary_to_msm_element_set = secondary_elem_to_msm_map_it->second;
1273 :
1274 : std::vector<std::unique_ptr<MortarSegmentHelper>> mortar_segment_helper(
1275 6149 : secondary_side_elem->n_sub_elem());
1276 6149 : const auto nodal_normals = getNodalNormals(*secondary_side_elem);
1277 :
1278 : /**
1279 : * Step 1.1: Linearize secondary face elements
1280 : *
1281 : * For first order face elements (Tri3 and Quad4) elements are simply linearized around center
1282 : * For second order (Tri6 and Quad9) and third order (Tri7) face elements, elements are
1283 : * sub-divided into four first order elements then each of the sub-elements is linearized
1284 : * around their respective centers
1285 : * For Quad8 elements, they are sub-divided into one quad and four triangle elements and each
1286 : * sub-element is linearized around their respective centers
1287 : */
1288 17878 : for (auto sel : make_range(secondary_side_elem->n_sub_elem()))
1289 : {
1290 : // Get indices of sub-element nodes in element
1291 : const auto sub_elem_nodes =
1292 11729 : Moose::Mortar::getMortarSubElementNodeIndices(*secondary_side_elem, sel);
1293 :
1294 : // Secondary sub-element center, normal, and nodes
1295 11729 : Point center;
1296 11729 : Point normal;
1297 11729 : std::vector<Point> nodes(sub_elem_nodes.size());
1298 :
1299 : // Collect the sub-element points and evaluate its center and averaged nodal normal.
1300 52457 : for (auto iv : make_range(sub_elem_nodes.size()))
1301 : {
1302 40728 : const auto n = sub_elem_nodes[iv];
1303 40728 : nodes[iv] = secondary_side_elem->point(n);
1304 40728 : center += secondary_side_elem->point(n);
1305 40728 : normal += nodal_normals[n];
1306 : }
1307 11729 : center /= sub_elem_nodes.size();
1308 11729 : normal = normal.unit();
1309 :
1310 11729 : if (_mortar_3d_subpatch_plane == Mortar3DSubpatchPlane::GEOMETRIC_NORMAL)
1311 : {
1312 11249 : const Point averaged_normal = normal;
1313 11249 : normal = get_sub_elem_geometric_normal(nodes);
1314 11249 : if (normal * averaged_normal < 0)
1315 396 : normal *= -1;
1316 : }
1317 :
1318 11729 : if (use_reference_interpolation)
1319 : {
1320 704 : std::vector<Point> sub_elem_reference_points;
1321 704 : sub_elem_reference_points.reserve(sub_elem_nodes.size());
1322 3008 : for (const auto node_index : sub_elem_nodes)
1323 2304 : sub_elem_reference_points.push_back(secondary_side_elem->master_point(node_index));
1324 :
1325 704 : mortar_segment_helper[sel] =
1326 1408 : std::make_unique<MortarSegmentHelper>(std::move(nodes),
1327 704 : std::move(sub_elem_reference_points),
1328 : center,
1329 : normal,
1330 704 : _triangulation_mode,
1331 1408 : _triangulate_triangles);
1332 704 : }
1333 : else
1334 22050 : mortar_segment_helper[sel] = std::make_unique<MortarSegmentHelper>(
1335 22050 : std::move(nodes), center, normal, _triangulation_mode, _triangulate_triangles);
1336 11729 : }
1337 :
1338 : /**
1339 : * Step 1.2: Coarse screening using a k-d tree to find nodes on the primary interface that are
1340 : * 'close to' a center point of the secondary element.
1341 : */
1342 :
1343 : // Search point for performing Nanoflann (k-d tree) searches.
1344 : // In each case we use the center point of the original element (not sub-elements for second
1345 : // order elements). This is to do search for all sub-elements simultaneously
1346 : std::array<Real, 3> query_pt;
1347 6149 : Point center_point;
1348 6149 : switch (secondary_side_elem->type())
1349 : {
1350 4542 : case TRI3:
1351 : case QUAD4:
1352 4542 : center_point = mortar_segment_helper[0]->center();
1353 4542 : query_pt = {{center_point(0), center_point(1), center_point(2)}};
1354 4542 : break;
1355 608 : case TRI6:
1356 : case TRI7:
1357 608 : center_point = mortar_segment_helper[1]->center();
1358 608 : query_pt = {{center_point(0), center_point(1), center_point(2)}};
1359 608 : break;
1360 759 : case QUAD8:
1361 759 : center_point = mortar_segment_helper[4]->center();
1362 759 : query_pt = {{center_point(0), center_point(1), center_point(2)}};
1363 759 : break;
1364 240 : case QUAD9:
1365 240 : center_point = secondary_side_elem->point(8);
1366 240 : query_pt = {{center_point(0), center_point(1), center_point(2)}};
1367 240 : break;
1368 0 : default:
1369 0 : mooseError(
1370 0 : "Face element type: ", secondary_side_elem->type(), "not supported for 3D mortar");
1371 : }
1372 :
1373 : // The number of results we want to get. These results will only be used to find
1374 : // a single element with non-trivial overlap, after an element is identified a breadth
1375 : // first search is done on neighbors
1376 6149 : const std::size_t num_results = 3;
1377 :
1378 : // Initialize result_set and do the search.
1379 18444 : std::vector<size_t> ret_index(num_results);
1380 12295 : std::vector<Real> out_dist_sqr(num_results);
1381 6149 : nanoflann::KNNResultSet<Real> result_set(num_results);
1382 6149 : result_set.init(&ret_index[0], &out_dist_sqr[0]);
1383 6149 : kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParameters());
1384 :
1385 : // Initialize list of processed primary elements, we don't want to revisit processed elements
1386 12295 : std::set<const Elem *, CompareDofObjectsByID> processed_primary_elems;
1387 :
1388 : // Initialize candidate set and flag for switching between coarse screening and breadth-first
1389 : // search
1390 6149 : bool primary_elem_found = false;
1391 12295 : std::set<const Elem *, CompareDofObjectsByID> primary_elem_candidates;
1392 6149 : const bool use_geometric_subpatch_normals =
1393 6149 : _mortar_3d_subpatch_plane == Mortar3DSubpatchPlane::GEOMETRIC_NORMAL;
1394 : // In geometric mode the projection-angle cutoff also rejects near-orthogonal subpatch pairs.
1395 : // The absolute dot product below keeps opposing primary/secondary orientations admissible.
1396 6149 : const Real minimum_subpatch_normal_alignment =
1397 6149 : use_geometric_subpatch_normals ? std::sin(_minimum_projection_angle * libMesh::pi / 180.0)
1398 : : 0.0;
1399 :
1400 : // Loop candidate nodes (returned by Nanoflann) and add all adjoining elems to candidate set
1401 24596 : for (auto r : make_range(result_set.size()))
1402 : {
1403 : // Verify that the squared distance we compute is the same as nanoflann's
1404 : mooseAssert(abs((_mesh.point(ret_index[r]) - center_point).norm_sq() - out_dist_sqr[r]) <=
1405 : TOLERANCE,
1406 : "Lower-dimensional element squared distance verification failed.");
1407 :
1408 : // Get list of elems connected to node
1409 : std::vector<const Elem *> & node_elems =
1410 18447 : this->_nodes_to_primary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));
1411 :
1412 : // Uniquely add elems to candidate set
1413 87123 : for (auto elem : node_elems)
1414 68676 : primary_elem_candidates.insert(elem);
1415 : }
1416 :
1417 : /**
1418 : * Step 1.3: Loop through primary candidate nodes, create mortar segments
1419 : *
1420 : * Once an element with non-trivial projection onto secondary element identified, switch
1421 : * to breadth-first search (drop all current candidates and add only neighbors of elements
1422 : * with non-trivial overlap)
1423 : */
1424 75892 : while (!primary_elem_candidates.empty())
1425 : {
1426 69743 : const Elem * primary_elem_candidate = *primary_elem_candidates.begin();
1427 :
1428 : // If we've already processed this candidate, we don't need to check it again.
1429 69743 : if (processed_primary_elems.count(primary_elem_candidate))
1430 : {
1431 0 : primary_elem_candidates.erase(primary_elem_candidate);
1432 0 : continue;
1433 : }
1434 :
1435 : // Initialize set of nodes used to construct mortar segment elements
1436 69743 : std::vector<Point> nodal_points;
1437 :
1438 : // Initialize map from mortar segment elements to nodes
1439 69743 : std::vector<std::vector<unsigned int>> elem_to_node_map;
1440 :
1441 : // Initialize list of secondary and primary sub-elements that formed each mortar segment
1442 69743 : std::vector<std::pair<unsigned int, unsigned int>> sub_elem_map;
1443 69743 : std::vector<std::array<Point, 3>> elem_to_secondary_reference_points;
1444 69743 : std::vector<std::array<Point, 3>> elem_to_primary_reference_points;
1445 :
1446 : /**
1447 : * Step 1.3.2: Sub-divide primary element candidate, then project onto secondary
1448 : * sub-elements, perform polygon clipping, and triangulate to form mortar segments
1449 : */
1450 196302 : for (auto p_el : make_range(primary_elem_candidate->n_sub_elem()))
1451 : {
1452 : // Get nodes of primary sub-elements
1453 : const auto sub_elem_nodes =
1454 126559 : Moose::Mortar::getMortarSubElementNodeIndices(*primary_elem_candidate, p_el);
1455 :
1456 : // Get list of primary sub-element vertex nodes
1457 126559 : std::vector<Point> primary_sub_elem(sub_elem_nodes.size());
1458 573446 : for (auto iv : make_range(sub_elem_nodes.size()))
1459 : {
1460 446887 : const auto n = sub_elem_nodes[iv];
1461 446887 : primary_sub_elem[iv] = primary_elem_candidate->point(n);
1462 : }
1463 126559 : Point primary_sub_elem_normal;
1464 126559 : if (use_geometric_subpatch_normals)
1465 120476 : primary_sub_elem_normal = get_sub_elem_geometric_normal(primary_sub_elem);
1466 :
1467 126559 : std::vector<Point> sub_elem_reference_points;
1468 126559 : if (use_reference_interpolation)
1469 : {
1470 8632 : sub_elem_reference_points.reserve(sub_elem_nodes.size());
1471 36984 : for (const auto node_index : sub_elem_nodes)
1472 28352 : sub_elem_reference_points.push_back(primary_elem_candidate->master_point(node_index));
1473 : }
1474 :
1475 : // Loop through secondary sub-elements
1476 520314 : for (auto s_el : make_range(secondary_side_elem->n_sub_elem()))
1477 : {
1478 : // Nearby primary candidates can include adjacent corner faces. Those faces may clip to
1479 : // numerical slivers, which we do not consider valid face-to-face mortar pairs for this
1480 : // search.
1481 781427 : if (use_geometric_subpatch_normals &&
1482 387672 : std::abs(primary_sub_elem_normal * mortar_segment_helper[s_el]->normal()) <
1483 : minimum_subpatch_normal_alignment)
1484 375 : continue;
1485 :
1486 : // Mortar segment helpers were defined for each secondary sub-element, they will:
1487 : // 1. Project primary sub-element onto linearized secondary sub-element
1488 : // 2. Clip projected primary sub-element against secondary sub-element
1489 : // 3. Triangulate clipped polygon to form mortar segments
1490 : //
1491 : // Mortar segment helpers append a list of mortar segment nodes and connectivities that
1492 : // can be directly used to build mortar segments
1493 393380 : const auto segments_before_helper = elem_to_node_map.size();
1494 393380 : if (use_reference_interpolation)
1495 41288 : mortar_segment_helper[s_el]->getMortarSegments(primary_sub_elem,
1496 : sub_elem_reference_points,
1497 : nodal_points,
1498 : elem_to_node_map,
1499 : elem_to_secondary_reference_points,
1500 : elem_to_primary_reference_points,
1501 : TOLERANCE * secondary_volume);
1502 : else
1503 352092 : mortar_segment_helper[s_el]->getMortarSegments(
1504 : primary_sub_elem, nodal_points, elem_to_node_map);
1505 :
1506 : // Keep track of which secondary and primary sub-elements created segment
1507 548490 : for (auto i = segments_before_helper; i < elem_to_node_map.size(); ++i)
1508 155110 : sub_elem_map.push_back(std::make_pair(s_el, p_el));
1509 : }
1510 126559 : }
1511 :
1512 : // Mark primary element as processed and remove from candidate list
1513 69743 : processed_primary_elems.insert(primary_elem_candidate);
1514 69743 : primary_elem_candidates.erase(primary_elem_candidate);
1515 :
1516 : // If overlap of polygons was non-trivial (created mortar segment elements)
1517 69743 : if (!elem_to_node_map.empty())
1518 : {
1519 26622 : if (sub_elem_map.size() != elem_to_node_map.size())
1520 0 : mooseError("The mortar segment subpatch map is not aligned with the mortar segment "
1521 : "connectivity map.");
1522 28034 : if (use_reference_interpolation &&
1523 1412 : (elem_to_secondary_reference_points.size() != elem_to_node_map.size() ||
1524 706 : elem_to_primary_reference_points.size() != elem_to_node_map.size()))
1525 0 : mooseError("The mortar segment reference-point maps are not aligned with the mortar "
1526 : "segment connectivity map.");
1527 :
1528 : // Only overlap polygons large enough to become mortar segments may switch the candidate
1529 : // search to breadth first.
1530 26622 : bool seed_breadth_first_search = false;
1531 26622 : std::vector<bool> retained_mortar_segments(elem_to_node_map.size(), false);
1532 181732 : for (const auto el : index_range(elem_to_node_map))
1533 : {
1534 155110 : const auto & node_map = elem_to_node_map[el];
1535 155110 : if (node_map.size() != 3)
1536 0 : mooseError(
1537 : "Active mortar segments only supports TRI elements, 3 nodes expected but: ",
1538 0 : node_map.size(),
1539 : " provided.");
1540 :
1541 155110 : const Point e1 = nodal_points[node_map[1]] - nodal_points[node_map[0]];
1542 155110 : const Point e2 = nodal_points[node_map[2]] - nodal_points[node_map[0]];
1543 0 : retained_mortar_segments[el] =
1544 155110 : 0.5 * e1.cross(e2).norm() / secondary_volume >= TOLERANCE;
1545 155110 : seed_breadth_first_search = seed_breadth_first_search || retained_mortar_segments[el];
1546 : }
1547 :
1548 26622 : if (seed_breadth_first_search)
1549 : {
1550 : // If this is the first element with a qualifying overlap, set flag. Candidates will
1551 : // now be neighbors of elements that had qualifying overlap.
1552 26554 : if (!primary_elem_found)
1553 : {
1554 6146 : primary_elem_found = true;
1555 6146 : primary_elem_candidates.clear();
1556 : }
1557 :
1558 : // Add neighbors to candidate list
1559 130416 : for (auto neighbor : primary_elem_candidate->neighbor_ptr_range())
1560 : {
1561 : // If not valid or not on lower dimensional secondary subdomain, skip
1562 103862 : if (neighbor == nullptr || neighbor->subdomain_id() != primary_subd_id)
1563 6730 : continue;
1564 : // If already processed, skip
1565 97132 : if (processed_primary_elems.count(neighbor))
1566 33079 : continue;
1567 : // Otherwise, add to candidates
1568 64053 : primary_elem_candidates.insert(neighbor);
1569 : }
1570 : }
1571 :
1572 : /**
1573 : * Step 1.3.3: Create mortar segments and add to mortar segment mesh
1574 : */
1575 26622 : std::vector<Node *> new_nodes;
1576 : // Clipping can append points for triangles later rejected by the area tolerance. Add only
1577 : // points referenced by retained triangles so the mortar mesh has no orphan nodes.
1578 26622 : std::vector<bool> retained_nodes(nodal_points.size(), false);
1579 181732 : for (const auto el : index_range(elem_to_node_map))
1580 155110 : if (retained_mortar_segments[el])
1581 617952 : for (const auto node : elem_to_node_map[el])
1582 463464 : retained_nodes[node] = true;
1583 :
1584 26622 : new_nodes.resize(nodal_points.size(), nullptr);
1585 256752 : for (const auto node : index_range(nodal_points))
1586 230130 : if (retained_nodes[node])
1587 457782 : new_nodes[node] = _mortar_segment_mesh->add_point(
1588 228891 : nodal_points[node], next_node_id++, secondary_side_elem->processor_id());
1589 :
1590 : // Loop through triangular elements in map
1591 181732 : for (auto el : index_range(elem_to_node_map))
1592 : {
1593 155110 : if (!retained_mortar_segments[el])
1594 622 : continue;
1595 :
1596 154488 : std::unique_ptr<Elem> new_elem;
1597 154488 : if (elem_to_node_map[el].size() == 3)
1598 154488 : new_elem = std::make_unique<Tri3>();
1599 : else
1600 0 : mooseError("Active mortar segments only supports TRI elements, 3 nodes expected "
1601 : "but: ",
1602 0 : elem_to_node_map[el].size(),
1603 : " provided.");
1604 :
1605 154488 : new_elem->processor_id() = secondary_side_elem->processor_id();
1606 154488 : new_elem->subdomain_id() = secondary_side_elem->subdomain_id();
1607 154488 : new_elem->set_id(next_elem_id++);
1608 :
1609 : // Attach newly created nodes
1610 617952 : for (auto i : index_range(elem_to_node_map[el]))
1611 463464 : new_elem->set_node(i, new_nodes[elem_to_node_map[el][i]]);
1612 :
1613 : // If element is smaller than tolerance, don't add to msm
1614 154488 : if (new_elem->volume() / secondary_volume < TOLERANCE)
1615 0 : continue;
1616 :
1617 : // Add elements to mortar segment mesh
1618 154488 : Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem.release());
1619 :
1620 154488 : msm_new_elem->set_extra_integer(secondary_sub_elem, sub_elem_map[el].first);
1621 154488 : msm_new_elem->set_extra_integer(primary_sub_elem, sub_elem_map[el].second);
1622 :
1623 : // Fill out mortar segment info
1624 154488 : MortarSegmentInfo msinfo;
1625 154488 : msinfo.secondary_elem = secondary_side_elem;
1626 154488 : msinfo.primary_elem = primary_elem_candidate;
1627 :
1628 : // Associate this MSM elem with the MortarSegmentInfo.
1629 154488 : _msm_elem_to_info.emplace(msm_new_elem, msinfo);
1630 :
1631 : // Store reference data only for retained segments.
1632 154488 : if (use_reference_interpolation)
1633 : {
1634 14678 : MortarSegmentReferencePoints reference_points{elem_to_secondary_reference_points[el],
1635 14678 : elem_to_primary_reference_points[el]};
1636 14678 : _msm_elem_to_reference_points.emplace(msm_new_elem, reference_points);
1637 : }
1638 :
1639 : // Add this mortar segment to the secondary elem to mortar segment map
1640 154488 : secondary_to_msm_element_set.insert(msm_new_elem);
1641 :
1642 154488 : _secondary_ip_sub_ids.insert(msinfo.secondary_elem->interior_parent()->subdomain_id());
1643 : // Unlike for 2D, we always have a primary when building the mortar mesh so we don't
1644 : // have to check for null
1645 154488 : _primary_ip_sub_ids.insert(msinfo.primary_elem->interior_parent()->subdomain_id());
1646 154488 : }
1647 26622 : }
1648 : // End loop through primary element candidates
1649 69743 : }
1650 :
1651 6149 : if (use_geometric_subpatch_normals)
1652 : {
1653 : // A geometric corner filter may intentionally leave individual subpatches uncovered. Warn
1654 : // only when the complete secondary element failed to produce a retained segment.
1655 5669 : if (secondary_to_msm_element_set.empty())
1656 3 : mooseDoOnce(
1657 : mooseWarning("Some secondary elements on mortar interface were unable to identify"
1658 : " a corresponding primary element; this may be expected depending on"
1659 : " problem geometry but may indicate a failure of the element search"
1660 : " or projection"));
1661 : }
1662 : else
1663 960 : for (auto sel : make_range(secondary_side_elem->n_sub_elem()))
1664 480 : if (mortar_segment_helper[sel]->remainder() == 1.0)
1665 0 : mooseDoOnce(
1666 : mooseWarning("Some secondary elements on mortar interface were unable to identify"
1667 : " a corresponding primary element; this may be expected depending on"
1668 : " problem geometry but may indicate a failure of the element search"
1669 : " or projection"));
1670 :
1671 6146 : if (secondary_to_msm_element_set.empty())
1672 0 : _secondary_elems_to_mortar_segments.erase(secondary_elem_to_msm_map_it);
1673 6501 : } // End loop through secondary elements
1674 355 : } // End loop through mortar constraint pairs
1675 :
1676 : mooseAssert(!use_reference_interpolation ||
1677 : _msm_elem_to_reference_points.size() == _msm_elem_to_info.size(),
1678 : "Mortar segment info and reference-point maps must remain aligned.");
1679 :
1680 355 : _mortar_segment_mesh->cache_elem_data();
1681 :
1682 : // The mesh was built distributedly (each rank owns only its local elements), so mark it
1683 : // as such so MeshSerializer correctly gathers it to proc 0 for Exodus output.
1684 355 : _mortar_segment_mesh->set_distributed();
1685 :
1686 : // Output mortar segment mesh
1687 355 : if (_debug)
1688 : {
1689 : // If element is not triangular, increment subdomain id
1690 : // (ExodusII does not support mixed element types in a single subdomain)
1691 57915 : for (const auto msm_el : _mortar_segment_mesh->active_local_element_ptr_range())
1692 57822 : if (msm_el->type() != TRI3)
1693 93 : msm_el->subdomain_id()++;
1694 :
1695 93 : outputMortarMesh();
1696 :
1697 : // Undo increment
1698 57915 : for (const auto msm_el : _mortar_segment_mesh->active_local_element_ptr_range())
1699 57822 : if (msm_el->type() != TRI3)
1700 93 : msm_el->subdomain_id()--;
1701 : }
1702 :
1703 355 : buildCouplingInformation();
1704 :
1705 : // Print mortar segment mesh statistics
1706 355 : if (_debug)
1707 : {
1708 93 : msmStatistics();
1709 : }
1710 355 : }
1711 :
1712 : void
1713 4618 : AutomaticMortarGeneration::buildCouplingInformation()
1714 : {
1715 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type>>>
1716 4618 : coupling_info;
1717 :
1718 : // Loop over the msm_elem_to_info object and build a bi-directional
1719 : // multimap from secondary elements to the primary Elems which they are
1720 : // coupled to and vice-versa. This is used in the
1721 : // AugmentSparsityOnInterface functor to determine whether a given
1722 : // secondary Elem is coupled across the mortar interface to a primary
1723 : // element.
1724 187808 : for (const auto & pr : _msm_elem_to_info)
1725 : {
1726 183190 : const Elem * secondary_elem = pr.second.secondary_elem;
1727 183190 : const Elem * primary_elem = pr.second.primary_elem;
1728 :
1729 : // LowerSecondary
1730 183190 : coupling_info[secondary_elem->processor_id()].emplace_back(
1731 183190 : secondary_elem->id(), secondary_elem->interior_parent()->id());
1732 183190 : if (secondary_elem->processor_id() != _mesh.processor_id())
1733 : // We want to keep information for nonlocal lower-dimensional secondary element point
1734 : // neighbors for mortar nodal aux kernels
1735 7871 : _mortar_interface_coupling[secondary_elem->id()].insert(
1736 7871 : secondary_elem->interior_parent()->id());
1737 :
1738 : // LowerPrimary
1739 183190 : coupling_info[secondary_elem->processor_id()].emplace_back(
1740 183190 : secondary_elem->id(), primary_elem->interior_parent()->id());
1741 183190 : if (secondary_elem->processor_id() != _mesh.processor_id())
1742 : // We want to keep information for nonlocal lower-dimensional secondary element point
1743 : // neighbors for mortar nodal aux kernels
1744 7871 : _mortar_interface_coupling[secondary_elem->id()].insert(
1745 7871 : primary_elem->interior_parent()->id());
1746 :
1747 : // Lower-LowerDimensionalPrimary
1748 366380 : coupling_info[secondary_elem->processor_id()].emplace_back(secondary_elem->id(),
1749 183190 : primary_elem->id());
1750 183190 : if (secondary_elem->processor_id() != _mesh.processor_id())
1751 : // We want to keep information for nonlocal lower-dimensional secondary element point
1752 : // neighbors for mortar nodal aux kernels
1753 7871 : _mortar_interface_coupling[secondary_elem->id()].insert(primary_elem->id());
1754 :
1755 : // SecondaryLower
1756 183190 : coupling_info[secondary_elem->interior_parent()->processor_id()].emplace_back(
1757 183190 : secondary_elem->interior_parent()->id(), secondary_elem->id());
1758 :
1759 : // SecondaryPrimary
1760 183190 : coupling_info[secondary_elem->interior_parent()->processor_id()].emplace_back(
1761 183190 : secondary_elem->interior_parent()->id(), primary_elem->interior_parent()->id());
1762 :
1763 : // PrimaryLower
1764 183190 : coupling_info[primary_elem->interior_parent()->processor_id()].emplace_back(
1765 183190 : primary_elem->interior_parent()->id(), secondary_elem->id());
1766 :
1767 : // PrimarySecondary
1768 183190 : coupling_info[primary_elem->interior_parent()->processor_id()].emplace_back(
1769 183190 : primary_elem->interior_parent()->id(), secondary_elem->interior_parent()->id());
1770 : }
1771 :
1772 : // Push the coupling information
1773 : auto action_functor =
1774 7073 : [this](processor_id_type,
1775 : const std::vector<std::pair<dof_id_type, dof_id_type>> & coupling_info)
1776 : {
1777 1289403 : for (auto [i, j] : coupling_info)
1778 1282330 : _mortar_interface_coupling[i].insert(j);
1779 7073 : };
1780 4618 : TIMPI::push_parallel_vector_data(_mesh.comm(), coupling_info, action_functor);
1781 4618 : }
1782 :
1783 : std::vector<AutomaticMortarGeneration::MsmSubdomainStats>
1784 135 : AutomaticMortarGeneration::computeMsmStatistics()
1785 : {
1786 135 : std::vector<MsmSubdomainStats> result;
1787 135 : StatisticsVector<Real> primary;
1788 135 : StatisticsVector<Real> secondary;
1789 135 : StatisticsVector<Real> msm;
1790 135 : std::unordered_map<dof_id_type, Real> primary_elems_to_volume;
1791 :
1792 270 : for (const auto & [primary_subd_id, secondary_subd_id] : _primary_secondary_subdomain_id_pairs)
1793 : {
1794 135 : for (const auto * const secondary_el :
1795 4086 : _mesh.active_local_subdomain_element_ptr_range(secondary_subd_id))
1796 : {
1797 1908 : secondary.push_back(secondary_el->volume());
1798 : // We may not have projected onto a primary face in which case we may not have created mortar
1799 : // segments
1800 1908 : if (auto it = _secondary_elems_to_mortar_segments.find(secondary_el->id());
1801 1908 : it != _secondary_elems_to_mortar_segments.end())
1802 61778 : for (const auto * const msm_elem : it->second)
1803 : {
1804 59870 : msm.push_back(msm_elem->volume());
1805 59870 : const auto & msm_info = libmesh_map_find(_msm_elem_to_info, msm_elem);
1806 : // Now it's also possible that we didn't project onto a primary face and we *did* create
1807 : // mortar segments
1808 59870 : if (msm_info.primary_elem)
1809 : {
1810 59870 : if (msm_info.primary_elem->subdomain_id() != primary_subd_id)
1811 0 : mooseError("Unhandled primary-secondary pairing when computing mortar segment "
1812 : "statistics. This could happen if you have the same secondary "
1813 : "lower-dimensional subdomain ID paired with multiple lower-dimensional "
1814 : "primary subdomain IDs. Contact a MOOSE developer for help.");
1815 59870 : if (const auto [it, inserted] =
1816 59870 : primary_elems_to_volume.emplace(msm_info.primary_elem->id(), Real{});
1817 59870 : inserted)
1818 3945 : it->second = msm_info.primary_elem->volume();
1819 : else
1820 : mooseAssert(
1821 : MooseUtils::absoluteFuzzyEqual(it->second, msm_info.primary_elem->volume()),
1822 : "Volumes should be consistent");
1823 : }
1824 : }
1825 135 : }
1826 :
1827 135 : _mesh.comm().set_union(primary_elems_to_volume);
1828 135 : _mesh.comm().allgather(static_cast<std::vector<Real> &>(secondary));
1829 135 : _mesh.comm().allgather(static_cast<std::vector<Real> &>(msm));
1830 135 : primary.reserve(primary_elems_to_volume.size());
1831 5679 : for (const auto [_, volume] : primary_elems_to_volume)
1832 5544 : primary.push_back(volume);
1833 :
1834 : MsmSubdomainStats stats;
1835 135 : stats.primary_subd_id = primary_subd_id;
1836 135 : stats.secondary_subd_id = secondary_subd_id;
1837 135 : stats.secondary_lower_n_elems = secondary.size();
1838 135 : stats.secondary_lower_max_volume = secondary.maximum();
1839 135 : stats.secondary_lower_min_volume = secondary.minimum();
1840 135 : stats.secondary_lower_median_volume = secondary.median();
1841 135 : stats.primary_lower_n_elems = primary.size();
1842 135 : stats.primary_lower_max_volume = primary.maximum();
1843 135 : stats.primary_lower_min_volume = primary.minimum();
1844 135 : stats.primary_lower_median_volume = primary.median();
1845 135 : stats.msm_n_elems = msm.size();
1846 135 : stats.msm_max_volume = msm.maximum();
1847 135 : stats.msm_min_volume = msm.minimum();
1848 135 : stats.msm_median_volume = msm.median();
1849 135 : result.push_back(stats);
1850 :
1851 135 : primary.clear();
1852 135 : secondary.clear();
1853 135 : msm.clear();
1854 135 : primary_elems_to_volume.clear();
1855 : }
1856 :
1857 270 : return result;
1858 135 : }
1859 :
1860 : void
1861 93 : AutomaticMortarGeneration::msmStatistics()
1862 : {
1863 93 : const auto all_stats = computeMsmStatistics();
1864 :
1865 93 : if (_mesh.processor_id() != 0)
1866 27 : return;
1867 :
1868 66 : Moose::out << "Mortar Interface Statistics:" << std::endl;
1869 132 : for (const auto & stats : all_stats)
1870 : {
1871 132 : std::vector<std::string> col_names = {"mesh", "n_elems", "max", "min", "median"};
1872 132 : std::vector<std::string> subds = {"secondary_lower", "primary_lower", "mortar_segment"};
1873 : std::vector<size_t> n_elems = {
1874 132 : stats.secondary_lower_n_elems, stats.primary_lower_n_elems, stats.msm_n_elems};
1875 : std::vector<Real> maxs = {
1876 132 : stats.secondary_lower_max_volume, stats.primary_lower_max_volume, stats.msm_max_volume};
1877 : std::vector<Real> mins = {
1878 132 : stats.secondary_lower_min_volume, stats.primary_lower_min_volume, stats.msm_min_volume};
1879 66 : std::vector<Real> medians = {stats.secondary_lower_median_volume,
1880 66 : stats.primary_lower_median_volume,
1881 132 : stats.msm_median_volume};
1882 :
1883 66 : FormattedTable table;
1884 66 : table.clear();
1885 264 : for (auto i : index_range(subds))
1886 : {
1887 198 : table.addRow(i);
1888 198 : table.addData<std::string>(col_names[0], subds[i]);
1889 198 : table.addData<size_t>(col_names[1], n_elems[i]);
1890 198 : table.addData<Real>(col_names[2], maxs[i]);
1891 198 : table.addData<Real>(col_names[3], mins[i]);
1892 198 : table.addData<Real>(col_names[4], medians[i]);
1893 : }
1894 :
1895 66 : Moose::out << "secondary subdomain: " << stats.secondary_subd_id
1896 66 : << " \tprimary subdomain: " << stats.primary_subd_id << std::endl;
1897 66 : table.printTable(Moose::out, subds.size());
1898 66 : }
1899 93 : }
1900 :
1901 : // The blocks marked with **** are for regressing edge dropping treatment and should be removed
1902 : // eventually.
1903 : //****
1904 : // Compute inactve nodes when the old (incorrect) edge dropping treatemnt is enabled
1905 : void
1906 808 : AutomaticMortarGeneration::computeIncorrectEdgeDroppingInactiveLMNodes()
1907 : {
1908 : using std::abs;
1909 :
1910 : // Note that in 3D our trick to check whether an element has edge dropping needs loose tolerances
1911 : // since the mortar segments are on the linearized element and comparing the volume of the
1912 : // linearized element does not have the same volume as the warped element
1913 808 : const Real tol = (dim() == 3) ? 0.1 : TOLERANCE;
1914 :
1915 808 : std::unordered_map<processor_id_type, std::set<dof_id_type>> proc_to_inactive_nodes_set;
1916 808 : const auto my_pid = _mesh.processor_id();
1917 :
1918 : // List of inactive nodes on local secondary elements
1919 808 : std::unordered_set<dof_id_type> inactive_node_ids;
1920 :
1921 808 : std::unordered_map<const Elem *, Real> active_volume{};
1922 :
1923 1616 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1924 6633 : for (const auto el : _mesh.active_subdomain_elements_ptr_range(pr.second))
1925 6633 : active_volume[el] = 0.;
1926 :
1927 : // Compute fraction of elements with corresponding primary elements
1928 11599 : for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
1929 : {
1930 10791 : const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
1931 10791 : const Elem * secondary_elem = msinfo.secondary_elem;
1932 :
1933 10791 : active_volume[secondary_elem] += msm_elem->volume();
1934 808 : }
1935 :
1936 : // Mark all inactive local nodes
1937 1616 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1938 : // Loop through all elements on my processor
1939 9842 : for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(pr.second))
1940 : // If elem fully or partially dropped
1941 4517 : if (abs(active_volume[el] / el->volume() - 1.0) > tol)
1942 : {
1943 : // Add all nodes to list of inactive
1944 0 : for (auto n : make_range(el->n_nodes()))
1945 0 : inactive_node_ids.insert(el->node_id(n));
1946 808 : }
1947 :
1948 : // Assemble list of procs that nodes contribute to
1949 1616 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1950 : {
1951 808 : const auto secondary_subd_id = pr.second;
1952 :
1953 : // Loop through all elements not on my processor
1954 12458 : for (const auto el : _mesh.active_subdomain_elements_ptr_range(secondary_subd_id))
1955 : {
1956 : // Get processor_id
1957 5825 : const auto pid = el->processor_id();
1958 :
1959 : // If element is in my subdomain, skip
1960 5825 : if (pid == my_pid)
1961 4517 : continue;
1962 :
1963 : // If element on proc pid shares any of my inactive nodes, mark to send
1964 5935 : for (const auto n : make_range(el->n_nodes()))
1965 : {
1966 4627 : const auto node_id = el->node_id(n);
1967 4627 : if (inactive_node_ids.find(node_id) != inactive_node_ids.end())
1968 0 : proc_to_inactive_nodes_set[pid].insert(node_id);
1969 : }
1970 808 : }
1971 : }
1972 :
1973 : // Send list of inactive nodes
1974 : {
1975 : // Pack set into vector for sending (push_parallel_vector_data doesn't like sets)
1976 808 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_inactive_nodes_vector;
1977 808 : for (const auto & proc_set : proc_to_inactive_nodes_set)
1978 0 : proc_to_inactive_nodes_vector[proc_set.first].insert(
1979 0 : proc_to_inactive_nodes_vector[proc_set.first].end(),
1980 : proc_set.second.begin(),
1981 : proc_set.second.end());
1982 :
1983 : // First push data
1984 0 : auto action_functor = [this, &inactive_node_ids](const processor_id_type pid,
1985 : const std::vector<dof_id_type> & sent_data)
1986 : {
1987 0 : if (pid == _mesh.processor_id())
1988 0 : mooseError("Should not be communicating with self.");
1989 0 : for (const auto pr : sent_data)
1990 0 : inactive_node_ids.insert(pr);
1991 0 : };
1992 808 : TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_inactive_nodes_vector, action_functor);
1993 808 : }
1994 808 : _inactive_local_lm_nodes.clear();
1995 808 : for (const auto node_id : inactive_node_ids)
1996 0 : _inactive_local_lm_nodes.insert(_mesh.node_ptr(node_id));
1997 808 : }
1998 :
1999 : void
2000 4618 : AutomaticMortarGeneration::computeInactiveLMNodes()
2001 : {
2002 4618 : if (!_correct_edge_dropping)
2003 : {
2004 808 : computeIncorrectEdgeDroppingInactiveLMNodes();
2005 808 : return;
2006 : }
2007 :
2008 3810 : std::unordered_map<processor_id_type, std::set<dof_id_type>> proc_to_active_nodes_set;
2009 3810 : const auto my_pid = _mesh.processor_id();
2010 :
2011 : // List of active nodes on local secondary elements
2012 3810 : std::unordered_set<dof_id_type> active_local_nodes;
2013 :
2014 : // Mark all active local nodes
2015 332866 : for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
2016 : {
2017 164528 : const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
2018 164528 : const Elem * secondary_elem = msinfo.secondary_elem;
2019 :
2020 1151670 : for (auto n : make_range(secondary_elem->n_nodes()))
2021 987142 : active_local_nodes.insert(secondary_elem->node_id(n));
2022 3810 : }
2023 :
2024 : // Assemble list of procs that nodes contribute to
2025 7620 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
2026 : {
2027 3810 : const auto secondary_subd_id = pr.second;
2028 :
2029 : // Loop through all elements not on my processor
2030 48950 : for (const auto el : _mesh.active_subdomain_elements_ptr_range(secondary_subd_id))
2031 : {
2032 : // Get processor_id
2033 22570 : const auto pid = el->processor_id();
2034 :
2035 : // If element is in my subdomain, skip
2036 22570 : if (pid == my_pid)
2037 15963 : continue;
2038 :
2039 : // If element on proc pid shares any of my active nodes, mark to send
2040 26677 : for (const auto n : make_range(el->n_nodes()))
2041 : {
2042 20070 : const auto node_id = el->node_id(n);
2043 20070 : if (active_local_nodes.find(node_id) != active_local_nodes.end())
2044 354 : proc_to_active_nodes_set[pid].insert(node_id);
2045 : }
2046 3810 : }
2047 : }
2048 :
2049 : // Send list of active nodes
2050 : {
2051 : // Pack set into vector for sending (push_parallel_vector_data doesn't like sets)
2052 3810 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_active_nodes_vector;
2053 3984 : for (const auto & proc_set : proc_to_active_nodes_set)
2054 : {
2055 174 : proc_to_active_nodes_vector[proc_set.first].reserve(proc_to_active_nodes_set.size());
2056 470 : for (const auto node_id : proc_set.second)
2057 296 : proc_to_active_nodes_vector[proc_set.first].push_back(node_id);
2058 : }
2059 :
2060 : // First push data
2061 174 : auto action_functor = [this, &active_local_nodes](const processor_id_type pid,
2062 : const std::vector<dof_id_type> & sent_data)
2063 : {
2064 174 : if (pid == _mesh.processor_id())
2065 0 : mooseError("Should not be communicating with self.");
2066 174 : active_local_nodes.insert(sent_data.begin(), sent_data.end());
2067 3984 : };
2068 3810 : TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_active_nodes_vector, action_functor);
2069 3810 : }
2070 :
2071 : // Every proc has correct list of active local nodes, now take complement (list of inactive nodes)
2072 : // and store to use later to zero LM DoFs on inactive nodes
2073 3810 : _inactive_local_lm_nodes.clear();
2074 7620 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
2075 3810 : for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(
2076 39546 : /*secondary_subd_id*/ pr.second))
2077 61205 : for (const auto n : make_range(el->n_nodes()))
2078 45242 : if (active_local_nodes.find(el->node_id(n)) == active_local_nodes.end())
2079 4223 : _inactive_local_lm_nodes.insert(el->node_ptr(n));
2080 3810 : }
2081 :
2082 : // Note: could be combined with previous routine, keeping separate for clarity (for now)
2083 : void
2084 4618 : AutomaticMortarGeneration::computeInactiveLMElems()
2085 : {
2086 : // Mark all active secondary elements
2087 4618 : std::unordered_set<const Elem *> active_local_elems;
2088 :
2089 : //****
2090 : // Note that in 3D our trick to check whether an element has edge dropping needs loose tolerances
2091 : // since the mortar segments are on the linearized element and comparing the volume of the
2092 : // linearized element does not have the same volume as the warped element
2093 4618 : const Real tol = (dim() == 3) ? 0.1 : TOLERANCE;
2094 :
2095 4618 : std::unordered_map<const Elem *, Real> active_volume;
2096 :
2097 : // Compute fraction of elements with corresponding primary elements
2098 4618 : if (!_correct_edge_dropping)
2099 11599 : for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
2100 : {
2101 10791 : const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
2102 10791 : const Elem * secondary_elem = msinfo.secondary_elem;
2103 :
2104 10791 : active_volume[secondary_elem] += msm_elem->volume();
2105 808 : }
2106 : //****
2107 :
2108 355256 : for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
2109 : {
2110 175319 : const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
2111 175319 : const Elem * secondary_elem = msinfo.secondary_elem;
2112 :
2113 : //****
2114 175319 : if (!_correct_edge_dropping)
2115 10791 : if (abs(active_volume[secondary_elem] / secondary_elem->volume() - 1.0) > tol)
2116 0 : continue;
2117 : //****
2118 :
2119 175319 : active_local_elems.insert(secondary_elem);
2120 4618 : }
2121 :
2122 : // Take complement of active elements in active local subdomain to get inactive local elements
2123 4618 : _inactive_local_lm_elems.clear();
2124 9236 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
2125 4618 : for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(
2126 50196 : /*secondary_subd_id*/ pr.second))
2127 20480 : if (active_local_elems.find(el) == active_local_elems.end())
2128 4866 : _inactive_local_lm_elems.insert(el);
2129 4618 : }
2130 :
2131 : void
2132 4621 : AutomaticMortarGeneration::computeNodalGeometry()
2133 : {
2134 : // The dimension according to Mesh::mesh_dimension().
2135 4621 : const auto dim = _mesh.mesh_dimension();
2136 :
2137 : mooseAssert(dim == 2 || dim == 3,
2138 : "AutomaticMortarGeneration::computeNodalGeometry() is only valid for "
2139 : "mortar constraints on 2D or 3D meshes.");
2140 : // A nodal lower-dimensional nodal quadrature rule to be used on faces.
2141 4621 : QNodal qface(dim - 1);
2142 :
2143 : // A map from the node id to the attached elemental normals/weights evaluated at the node. Th
2144 : // length of the vector will correspond to the number of elements attached to the node. If it is a
2145 : // vertex node, for a 1D mortar mesh, the vector length will be two. If it is an interior node,
2146 : // the vector will be length 1. The first member of the pair is that element's normal at the node.
2147 : // The second member is that element's JxW at the node
2148 4621 : std::map<dof_id_type, std::vector<std::pair<Point, Real>>> node_to_normals_map;
2149 :
2150 : /// The _periodic flag tells us whether we want to inward vs outward facing normals
2151 4621 : Real sign = _periodic ? -1 : 1;
2152 :
2153 : // First loop over lower-dimensional secondary side elements and compute/save the outward normal
2154 : // for each one. We loop over all active elements currently, but this procedure could be
2155 : // parallelized as well.
2156 4621 : for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
2157 4621 : end_el = _mesh.active_elements_end();
2158 476761 : el != end_el;
2159 472140 : ++el)
2160 : {
2161 472140 : const Elem * secondary_elem = *el;
2162 :
2163 : // If this is not one of the lower-dimensional secondary side elements, go on to the next one.
2164 472140 : if (!_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
2165 443739 : continue;
2166 :
2167 : // We will create an FE object and attach the nodal quadrature rule such that we can get out the
2168 : // normals at the element nodes
2169 28401 : FEType nnx_fe_type(secondary_elem->default_order(), LAGRANGE);
2170 28401 : std::unique_ptr<FEBase> nnx_fe_face(FEBase::build(dim, nnx_fe_type));
2171 28401 : nnx_fe_face->attach_quadrature_rule(&qface);
2172 28401 : const auto & face_normals = nnx_fe_face->get_normals();
2173 28401 : const auto & face_points = nnx_fe_face->get_xyz();
2174 :
2175 28401 : const auto & JxW = nnx_fe_face->get_JxW();
2176 :
2177 : // Which side of the parent are we? We need to know this to know
2178 : // which side to reinit.
2179 28401 : const Elem * interior_parent = secondary_elem->interior_parent();
2180 : mooseAssert(interior_parent,
2181 : "No interior parent exists for element "
2182 : << secondary_elem->id()
2183 : << ". There may be a problem with your sideset set-up.");
2184 :
2185 : // Map to get lower dimensional element from interior parent on secondary surface
2186 : // This map can be used to provide a handle to methods in this class that need to
2187 : // operate on lower dimensional elements.
2188 28401 : _secondary_element_to_secondary_lowerd_element.emplace(interior_parent->id(), secondary_elem);
2189 :
2190 : // Look up which side of the interior parent secondary_elem is.
2191 28401 : auto s = interior_parent->which_side_am_i(secondary_elem);
2192 :
2193 : // Reinit the face FE object on side s.
2194 28401 : nnx_fe_face->reinit(interior_parent, s);
2195 :
2196 : // Match by physical location instead of assuming that parent-side nodal
2197 : // quadrature ordering and lower-dimensional side-element node ordering are
2198 : // identical.
2199 : const auto qpoint_to_secondary_node =
2200 28401 : nodalQuadraturePointToSecondaryNodeMap(*secondary_elem, face_points);
2201 :
2202 : mooseAssert(face_normals.size() == face_points.size() && JxW.size() == face_points.size(),
2203 : "Face nodal geometry vectors must have the same size.");
2204 :
2205 112905 : for (const auto qp : make_range(face_points.size()))
2206 : {
2207 84504 : const auto n = qpoint_to_secondary_node[qp];
2208 84504 : auto & normals_and_weights_vec = node_to_normals_map[secondary_elem->node_id(n)];
2209 84504 : normals_and_weights_vec.push_back(std::make_pair(sign * face_normals[qp], JxW[qp]));
2210 : }
2211 33022 : }
2212 :
2213 47673 : for (const auto & pr : node_to_normals_map)
2214 : {
2215 : // Compute normal vector
2216 43052 : const auto & node_id = pr.first;
2217 43052 : const auto & normals_and_weights_vec = pr.second;
2218 :
2219 43052 : Point nodal_normal;
2220 127556 : for (const auto & norm_and_weight : normals_and_weights_vec)
2221 84504 : nodal_normal += norm_and_weight.first * norm_and_weight.second;
2222 43052 : nodal_normal = nodal_normal.unit();
2223 :
2224 43052 : _secondary_node_to_nodal_normal[_mesh.node_ptr(node_id)] = nodal_normal;
2225 :
2226 43052 : Point nodal_tangent_one;
2227 43052 : Point nodal_tangent_two;
2228 43052 : householderOrthogolization(nodal_normal, nodal_tangent_one, nodal_tangent_two);
2229 :
2230 43052 : _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)][0] = nodal_tangent_one;
2231 43052 : _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)][1] = nodal_tangent_two;
2232 : }
2233 4621 : }
2234 :
2235 : void
2236 43052 : AutomaticMortarGeneration::householderOrthogolization(const Point & nodal_normal,
2237 : Point & nodal_tangent_one,
2238 : Point & nodal_tangent_two) const
2239 : {
2240 : using std::abs;
2241 :
2242 : mooseAssert(MooseUtils::absoluteFuzzyEqual(nodal_normal.norm(), 1),
2243 : "The input nodal normal should have unity norm");
2244 :
2245 43052 : const Real nx = nodal_normal(0);
2246 43052 : const Real ny = nodal_normal(1);
2247 43052 : const Real nz = nodal_normal(2);
2248 :
2249 : // See Lopes DS, Silva MT, Ambrosio JA. Tangent vectors to a 3-D surface normal: A geometric tool
2250 : // to find orthogonal vectors based on the Householder transformation. Computer-Aided Design. 2013
2251 : // Mar 1;45(3):683-94. We choose one definition of h_vector and deal with special case.
2252 43052 : const Point h_vector(nx + 1.0, ny, nz);
2253 :
2254 : // Avoid singularity of the equations at the end of routine by providing the solution to
2255 : // (nx,ny,nz)=(-1,0,0) Normal/tangent fields can be visualized by outputting nodal geometry mesh
2256 : // on a spherical problem.
2257 43052 : if (abs(h_vector(0)) < TOLERANCE)
2258 : {
2259 1878 : nodal_tangent_one(0) = 0;
2260 1878 : nodal_tangent_one(1) = 1;
2261 1878 : nodal_tangent_one(2) = 0;
2262 :
2263 1878 : nodal_tangent_two(0) = 0;
2264 1878 : nodal_tangent_two(1) = 0;
2265 1878 : nodal_tangent_two(2) = -1;
2266 :
2267 1878 : return;
2268 : }
2269 :
2270 41174 : const Real h = h_vector.norm();
2271 :
2272 41174 : nodal_tangent_one(0) = -2.0 * h_vector(0) * h_vector(1) / (h * h);
2273 41174 : nodal_tangent_one(1) = 1.0 - 2.0 * h_vector(1) * h_vector(1) / (h * h);
2274 41174 : nodal_tangent_one(2) = -2.0 * h_vector(1) * h_vector(2) / (h * h);
2275 :
2276 41174 : nodal_tangent_two(0) = -2.0 * h_vector(0) * h_vector(2) / (h * h);
2277 41174 : nodal_tangent_two(1) = -2.0 * h_vector(1) * h_vector(2) / (h * h);
2278 41174 : nodal_tangent_two(2) = 1.0 - 2.0 * h_vector(2) * h_vector(2) / (h * h);
2279 : }
2280 :
2281 : // Project secondary nodes onto their corresponding primary elements for each primary/secondary
2282 : // pair.
2283 : void
2284 4263 : AutomaticMortarGeneration::projectSecondaryNodes()
2285 : {
2286 : // For each primary/secondary boundary id pair, call the
2287 : // project_secondary_nodes_single_pair() helper function.
2288 8526 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
2289 4263 : projectSecondaryNodesSinglePair(pr.first, pr.second);
2290 4263 : }
2291 :
2292 : bool
2293 5153 : AutomaticMortarGeneration::processAlignedNodes(
2294 : const Node & secondary_node,
2295 : const Node & primary_node,
2296 : const std::vector<const Elem *> * secondary_node_neighbors,
2297 : const std::vector<const Elem *> * primary_node_neighbors,
2298 : const VectorValue<Real> & nodal_normal,
2299 : const Elem & candidate_element,
2300 : std::set<const Elem *> & rejected_elem_candidates)
2301 : {
2302 5153 : if (!secondary_node_neighbors)
2303 0 : secondary_node_neighbors = &libmesh_map_find(_nodes_to_secondary_elem_map, secondary_node.id());
2304 5153 : if (!primary_node_neighbors)
2305 5153 : primary_node_neighbors = &libmesh_map_find(_nodes_to_primary_elem_map, primary_node.id());
2306 :
2307 5153 : std::vector<bool> primary_elems_mapped(primary_node_neighbors->size(), false);
2308 :
2309 : // Add entries to secondary_node_and_elem_to_xi2_primary_elem container.
2310 : //
2311 : // First, determine "on left" vs. "on right" orientation of the nodal neighbors.
2312 : // There can be a max of 2 nodal neighbors, and we want to make sure that the
2313 : // secondary nodal neighbor on the "left" is associated with the primary nodal
2314 : // neighbor on the "left" and similarly for the "right". We use cross products to determine
2315 : // alignment. In the below diagram, 'x' denotes a node, and connected '|' are lower dimensional
2316 : // elements.
2317 : // x
2318 : // x |
2319 : // | |
2320 : // secondary x ----> x primary
2321 : // | |
2322 : // | x
2323 : // x
2324 : //
2325 : // Looking at the aligned nodes, the secondary node first, if we pick the top secondary lower
2326 : // dimensional element, then the cross product as written a few lines below points out of the
2327 : // screen towards you. (Point in the direction of the secondary nodal normal, and then curl your
2328 : // hand towards the secondary element's opposite node, then the thumb points in the direction of
2329 : // the cross product). Doing the same with the aligned primary node, if we pick the top primary
2330 : // element, then the cross product also points out of the screen. Because the cross products
2331 : // point in the same direction (positive dot product), then we know to associate the
2332 : // secondary-primary element pair. If we had picked the bottom primary element whose cross
2333 : // product points into the screen, then clearly the cross products point in the opposite
2334 : // direction and we don't have a match
2335 : std::array<Real, 2> secondary_node_neighbor_cps, primary_node_neighbor_cps;
2336 :
2337 13105 : for (const auto nn : index_range(*secondary_node_neighbors))
2338 : {
2339 7952 : const Elem * const secondary_neigh = (*secondary_node_neighbors)[nn];
2340 7952 : const Point opposite = (secondary_neigh->node_ptr(0) == &secondary_node)
2341 7952 : ? secondary_neigh->point(1)
2342 3980 : : secondary_neigh->point(0);
2343 7952 : const Point cp = nodal_normal.cross(opposite - secondary_node);
2344 7952 : secondary_node_neighbor_cps[nn] = cp(2);
2345 : }
2346 :
2347 12879 : for (const auto nn : index_range(*primary_node_neighbors))
2348 : {
2349 7726 : const Elem * const primary_neigh = (*primary_node_neighbors)[nn];
2350 7726 : const Point opposite = (primary_neigh->node_ptr(0) == &primary_node) ? primary_neigh->point(1)
2351 3980 : : primary_neigh->point(0);
2352 7726 : const Point cp = nodal_normal.cross(opposite - primary_node);
2353 7726 : primary_node_neighbor_cps[nn] = cp(2);
2354 : }
2355 :
2356 : // Associate secondary/primary elems on matching sides.
2357 5153 : bool found_match = false;
2358 13105 : for (const auto snn : index_range(*secondary_node_neighbors))
2359 21050 : for (const auto mnn : index_range(*primary_node_neighbors))
2360 13098 : if (secondary_node_neighbor_cps[snn] * primary_node_neighbor_cps[mnn] > 0)
2361 : {
2362 7714 : found_match = true;
2363 7714 : if (primary_elems_mapped[mnn])
2364 0 : continue;
2365 7714 : primary_elems_mapped[mnn] = true;
2366 :
2367 : // Figure out xi^(2) value by looking at which node primary_node is
2368 : // of the current primary node neighbor.
2369 7714 : const Real xi2 = (&primary_node == (*primary_node_neighbors)[mnn]->node_ptr(0)) ? -1 : +1;
2370 : const auto secondary_key =
2371 7714 : std::make_pair(&secondary_node, (*secondary_node_neighbors)[snn]);
2372 7714 : const auto primary_val = std::make_pair(xi2, (*primary_node_neighbors)[mnn]);
2373 7714 : _secondary_node_and_elem_to_xi2_primary_elem.emplace(secondary_key, primary_val);
2374 :
2375 : // Also map in the other direction.
2376 : const Real xi1 =
2377 7714 : (&secondary_node == (*secondary_node_neighbors)[snn]->node_ptr(0)) ? -1 : +1;
2378 :
2379 : const auto primary_key =
2380 7714 : std::make_tuple(primary_node.id(), &primary_node, (*primary_node_neighbors)[mnn]);
2381 7714 : const auto secondary_val = std::make_pair(xi1, (*secondary_node_neighbors)[snn]);
2382 7714 : _primary_node_and_elem_to_xi1_secondary_elem.emplace(primary_key, secondary_val);
2383 : }
2384 :
2385 5153 : if (!found_match)
2386 : {
2387 : // There could be coincident nodes and this might be a bad primary candidate (see
2388 : // issue #21680). Instead of giving up, let's try continuing
2389 12 : rejected_elem_candidates.insert(&candidate_element);
2390 12 : return false;
2391 : }
2392 :
2393 : // We need to handle the case where we've exactly projected a secondary node onto a
2394 : // primary node, but our secondary node is at one of the secondary boundary face endpoints and
2395 : // our primary node is not.
2396 5141 : if (secondary_node_neighbors->size() == 1 && primary_node_neighbors->size() == 2)
2397 0 : for (const auto i : index_range(primary_elems_mapped))
2398 0 : if (!primary_elems_mapped[i])
2399 : {
2400 0 : _primary_node_and_elem_to_xi1_secondary_elem.emplace(
2401 0 : std::make_tuple(primary_node.id(), &primary_node, (*primary_node_neighbors)[i]),
2402 0 : std::make_pair(1, nullptr));
2403 : }
2404 :
2405 5141 : return found_match;
2406 5153 : }
2407 :
2408 : void
2409 4263 : AutomaticMortarGeneration::projectSecondaryNodesSinglePair(
2410 : SubdomainID lower_dimensional_primary_subdomain_id,
2411 : SubdomainID lower_dimensional_secondary_subdomain_id)
2412 : {
2413 : using std::abs;
2414 :
2415 : // Build the "subdomain" adaptor based KD Tree.
2416 4263 : NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, lower_dimensional_primary_subdomain_id);
2417 : subdomain_kd_tree_t kd_tree(
2418 4263 : 3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
2419 :
2420 : // Construct the KD tree.
2421 4263 : kd_tree.buildIndex();
2422 :
2423 4263 : for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
2424 4263 : end_el = _mesh.active_elements_end();
2425 323189 : el != end_el;
2426 318926 : ++el)
2427 : {
2428 318926 : const Elem * secondary_side_elem = *el;
2429 :
2430 : // If this Elem is not in the current secondary subdomain, go on to the next one.
2431 318926 : if (secondary_side_elem->subdomain_id() != lower_dimensional_secondary_subdomain_id)
2432 299218 : continue;
2433 :
2434 : // For each node on the lower-dimensional element, find the nearest
2435 : // node on the primary side using the KDTree, then
2436 : // search in nearby elements for where it projects
2437 : // along the nodal normal direction.
2438 59124 : for (MooseIndex(secondary_side_elem->n_vertices()) n = 0; n < secondary_side_elem->n_vertices();
2439 : ++n)
2440 : {
2441 39416 : const Node * secondary_node = secondary_side_elem->node_ptr(n);
2442 :
2443 : // Get the nodal neighbors for secondary_node, so we can check whether we've
2444 : // already successfully projected it.
2445 : const std::vector<const Elem *> & secondary_node_neighbors =
2446 39416 : this->_nodes_to_secondary_elem_map.at(secondary_node->id());
2447 :
2448 : // Check whether we've already mapped this secondary node
2449 : // successfully for all of its nodal neighbors.
2450 39416 : bool is_mapped = true;
2451 69674 : for (MooseIndex(secondary_node_neighbors) snn = 0; snn < secondary_node_neighbors.size();
2452 : ++snn)
2453 : {
2454 54579 : auto secondary_key = std::make_pair(secondary_node, secondary_node_neighbors[snn]);
2455 54579 : if (!_secondary_node_and_elem_to_xi2_primary_elem.count(secondary_key))
2456 : {
2457 24321 : is_mapped = false;
2458 24321 : break;
2459 : }
2460 : }
2461 :
2462 : // Go to the next node if this one has already been mapped.
2463 39416 : if (is_mapped)
2464 15095 : continue;
2465 :
2466 : // Look up the new nodal normal value in the local storage, error if not found.
2467 24321 : Point nodal_normal = _secondary_node_to_nodal_normal.at(secondary_node);
2468 :
2469 : // Data structure for performing Nanoflann searches.
2470 : std::array<Real, 3> query_pt = {
2471 24321 : {(*secondary_node)(0), (*secondary_node)(1), (*secondary_node)(2)}};
2472 :
2473 : // The number of results we want to get. We'll look for a
2474 : // "few" nearest nodes, hopefully that is enough to let us
2475 : // figure out which lower-dimensional Elem on the primary
2476 : // side we are across from.
2477 24321 : const std::size_t num_results = 3;
2478 :
2479 : // Initialize result_set and do the search.
2480 48642 : std::vector<size_t> ret_index(num_results);
2481 24321 : std::vector<Real> out_dist_sqr(num_results);
2482 24321 : nanoflann::KNNResultSet<Real> result_set(num_results);
2483 24321 : result_set.init(&ret_index[0], &out_dist_sqr[0]);
2484 24321 : kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParameters());
2485 :
2486 : // If this flag gets set in the loop below, we can break out of the outer r-loop as well.
2487 24321 : bool projection_succeeded = false;
2488 :
2489 : // Once we've rejected a candidate for a given secondary_node,
2490 : // there's no reason to check it again.
2491 24321 : std::set<const Elem *> rejected_primary_elem_candidates;
2492 :
2493 : // Loop over the closest nodes, check whether
2494 : // the secondary node successfully projects into
2495 : // either of the closest neighbors, stop when
2496 : // the projection succeeds.
2497 34995 : for (MooseIndex(result_set) r = 0; r < result_set.size(); ++r)
2498 : {
2499 : // Verify that the squared distance we compute is the same as nanoflann'sFss
2500 : mooseAssert(abs((_mesh.point(ret_index[r]) - *secondary_node).norm_sq() -
2501 : out_dist_sqr[r]) <= TOLERANCE,
2502 : "Lower-dimensional element squared distance verification failed.");
2503 :
2504 : // Get a reference to the vector of lower dimensional elements from the
2505 : // nodes_to_primary_elem_map.
2506 : std::vector<const Elem *> & primary_elem_candidates =
2507 31441 : this->_nodes_to_primary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));
2508 :
2509 : // Search the Elems connected to this node on the primary mesh side.
2510 51139 : for (MooseIndex(primary_elem_candidates) e = 0; e < primary_elem_candidates.size(); ++e)
2511 : {
2512 40465 : const Elem * primary_elem_candidate = primary_elem_candidates[e];
2513 :
2514 : // If we've already rejected this candidate, we don't need to check it again.
2515 40465 : if (rejected_primary_elem_candidates.count(primary_elem_candidate))
2516 7120 : continue;
2517 :
2518 : // Now generically solve for xi2
2519 33357 : const auto order = primary_elem_candidate->default_order();
2520 33357 : DualNumber<Real> xi2_dn{0, 1};
2521 33357 : unsigned int current_iterate = 0, max_iterates = 10;
2522 :
2523 : // Newton loop
2524 : do
2525 : {
2526 65831 : VectorValue<DualNumber<Real>> x2(0);
2527 65831 : for (MooseIndex(primary_elem_candidate->n_nodes()) n = 0;
2528 203157 : n < primary_elem_candidate->n_nodes();
2529 : ++n)
2530 : x2 +=
2531 137326 : Moose::fe_lagrange_1D_shape(order, n, xi2_dn) * primary_elem_candidate->point(n);
2532 65831 : const auto u = x2 - (*secondary_node);
2533 65831 : const auto F = u(0) * nodal_normal(1) - u(1) * nodal_normal(0);
2534 :
2535 65831 : if (abs(F) < _newton_tolerance)
2536 33357 : break;
2537 :
2538 32474 : if (F.derivatives())
2539 : {
2540 32474 : Real dxi2 = -F.value() / F.derivatives();
2541 :
2542 32474 : xi2_dn += dxi2;
2543 : }
2544 : else
2545 : // It's possible that the secondary surface nodal normal is completely orthogonal to
2546 : // the primary surface normal, in which case the derivative is 0. We know in this case
2547 : // that the projection should be a failure
2548 0 : current_iterate = max_iterates;
2549 165019 : } while (++current_iterate < max_iterates);
2550 :
2551 33357 : Real xi2 = xi2_dn.value();
2552 :
2553 : // Check whether the projection worked. The last condition checks for obliqueness of the
2554 : // projection
2555 : //
2556 : // We are projecting on one side first and the other side second. If we make the
2557 : // tolerance bigger and remove the (5) factor we are going to continue to miss the
2558 : // second projection and fall into the exception message in
2559 : // projectPrimaryNodesSinglePair. What makes this modification to not fall in the
2560 : // exception is that we are projecting on one side more xi than in the other. There
2561 : // should be a better way of doing this by using actual distances and not parametric
2562 : // coordinates. But I believe making the tolerance uniformly larger or smaller won't do
2563 : // the trick here.
2564 54136 : if ((current_iterate < max_iterates) && (std::abs(xi2) <= 1. + 5 * _xi_tolerance) &&
2565 54136 : (abs((primary_elem_candidate->point(0) - primary_elem_candidate->point(1)).unit() *
2566 20779 : nodal_normal) < std::cos(_minimum_projection_angle * libMesh::pi / 180)))
2567 : {
2568 : // If xi2 == +1 or -1 then this secondary node mapped directly to a node on the primary
2569 : // surface. This isn't as unlikely as you might think, it will happen if the meshes
2570 : // on the interface start off being perfectly aligned. In this situation, we need to
2571 : // associate the secondary node with two different elements (and two corresponding
2572 : // xi^(2) values.
2573 20779 : if (abs(abs(xi2) - 1.) <= _xi_tolerance * 5.0)
2574 : {
2575 5153 : const Node * primary_node = (xi2 < 0) ? primary_elem_candidate->node_ptr(0)
2576 2926 : : primary_elem_candidate->node_ptr(1);
2577 : const bool created_mortar_segment =
2578 5153 : processAlignedNodes(*secondary_node,
2579 : *primary_node,
2580 : &secondary_node_neighbors,
2581 : nullptr,
2582 : nodal_normal,
2583 : *primary_elem_candidate,
2584 : rejected_primary_elem_candidates);
2585 :
2586 5153 : if (!created_mortar_segment)
2587 12 : continue;
2588 : }
2589 : else // Point falls somewhere in the middle of the Elem.
2590 : {
2591 : // Add two entries to secondary_node_and_elem_to_xi2_primary_elem.
2592 43774 : for (MooseIndex(secondary_node_neighbors) nn = 0;
2593 43774 : nn < secondary_node_neighbors.size();
2594 : ++nn)
2595 : {
2596 28148 : const Elem * neigh = secondary_node_neighbors[nn];
2597 84444 : for (MooseIndex(neigh->n_vertices()) nid = 0; nid < neigh->n_vertices(); ++nid)
2598 : {
2599 56296 : const Node * neigh_node = neigh->node_ptr(nid);
2600 56296 : if (secondary_node == neigh_node)
2601 : {
2602 28148 : auto key = std::make_pair(neigh_node, neigh);
2603 28148 : auto val = std::make_pair(xi2, primary_elem_candidate);
2604 28148 : _secondary_node_and_elem_to_xi2_primary_elem.emplace(key, val);
2605 : }
2606 : }
2607 : }
2608 : }
2609 :
2610 20767 : projection_succeeded = true;
2611 20767 : break; // out of e-loop
2612 : }
2613 : else
2614 : // The current secondary_node is not in this Elem, so keep track of the rejects.
2615 12578 : rejected_primary_elem_candidates.insert(primary_elem_candidate);
2616 33357 : }
2617 :
2618 31441 : if (projection_succeeded)
2619 20767 : break; // out of r-loop
2620 : } // r-loop
2621 :
2622 24321 : if (!projection_succeeded)
2623 : {
2624 3554 : _failed_secondary_node_projections.insert(secondary_node->id());
2625 3554 : if (_debug)
2626 0 : _console << "Failed to find primary Elem into which secondary node "
2627 0 : << static_cast<const Point &>(*secondary_node) << ", id '"
2628 0 : << secondary_node->id() << "', projects onto\n"
2629 0 : << std::endl;
2630 : }
2631 20767 : else if (_debug)
2632 48 : _projected_secondary_nodes.insert(secondary_node->id());
2633 24321 : } // loop over side nodes
2634 4263 : } // end loop over lower-dimensional elements
2635 :
2636 4263 : if (_distributed)
2637 : {
2638 96 : if (_debug)
2639 2 : _mesh.comm().set_union(_projected_secondary_nodes);
2640 96 : _mesh.comm().set_union(_failed_secondary_node_projections);
2641 : }
2642 :
2643 4263 : if (_debug)
2644 12 : _console << "\n"
2645 12 : << _projected_secondary_nodes.size() << " out of "
2646 12 : << _projected_secondary_nodes.size() + _failed_secondary_node_projections.size()
2647 12 : << " secondary nodes were successfully projected\n"
2648 12 : << std::endl;
2649 4263 : }
2650 :
2651 : // Inverse map primary nodes onto their corresponding secondary elements for each primary/secondary
2652 : // pair.
2653 : void
2654 4263 : AutomaticMortarGeneration::projectPrimaryNodes()
2655 : {
2656 : // For each primary/secondary boundary id pair, call the
2657 : // project_primary_nodes_single_pair() helper function.
2658 8526 : for (const auto & pr : _primary_secondary_subdomain_id_pairs)
2659 4263 : projectPrimaryNodesSinglePair(pr.first, pr.second);
2660 4263 : }
2661 :
2662 : void
2663 4263 : AutomaticMortarGeneration::projectPrimaryNodesSinglePair(
2664 : SubdomainID lower_dimensional_primary_subdomain_id,
2665 : SubdomainID lower_dimensional_secondary_subdomain_id)
2666 : {
2667 : using std::abs;
2668 :
2669 : // Build a Nanoflann object on the lower-dimensional secondary elements of the Mesh.
2670 4263 : NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, lower_dimensional_secondary_subdomain_id);
2671 : subdomain_kd_tree_t kd_tree(
2672 4263 : 3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
2673 :
2674 : // Construct the KD tree for lower-dimensional elements in the volume mesh.
2675 4263 : kd_tree.buildIndex();
2676 :
2677 4263 : std::unordered_set<dof_id_type> primary_nodes_visited;
2678 :
2679 323189 : for (const auto & primary_side_elem : _mesh.active_element_ptr_range())
2680 : {
2681 : // If this is not one of the lower-dimensional primary side elements, go on to the next one.
2682 318926 : if (primary_side_elem->subdomain_id() != lower_dimensional_primary_subdomain_id)
2683 302408 : continue;
2684 :
2685 : // For each node on this side, find the nearest node on the secondary side using the KDTree,
2686 : // then search in nearby elements for where it projects along the nodal normal direction.
2687 49554 : for (MooseIndex(primary_side_elem->n_vertices()) n = 0; n < primary_side_elem->n_vertices();
2688 : ++n)
2689 : {
2690 : // Get a pointer to this node.
2691 33036 : const Node * primary_node = primary_side_elem->node_ptr(n);
2692 :
2693 : // Get the nodal neighbors connected to this primary node.
2694 : const std::vector<const Elem *> & primary_node_neighbors =
2695 33036 : _nodes_to_primary_elem_map.at(primary_node->id());
2696 :
2697 : // Check whether we have already successfully inverse mapped this primary node (whether during
2698 : // secondary node projection or now during primary node projection) or we have already failed
2699 : // to inverse map this primary node (now during primary node projection), and then skip if
2700 : // either of those things is true
2701 : auto primary_key =
2702 33036 : std::make_tuple(primary_node->id(), primary_node, primary_node_neighbors[0]);
2703 53829 : if (!primary_nodes_visited.insert(primary_node->id()).second ||
2704 20793 : _primary_node_and_elem_to_xi1_secondary_elem.count(primary_key))
2705 17271 : continue;
2706 :
2707 : // Data structure for performing Nanoflann searches.
2708 15765 : Real query_pt[3] = {(*primary_node)(0), (*primary_node)(1), (*primary_node)(2)};
2709 :
2710 : // The number of results we want to get. We'll look for a
2711 : // "few" nearest nodes, hopefully that is enough to let us
2712 : // figure out which lower-dimensional Elem on the secondary side
2713 : // we are across from.
2714 15765 : const size_t num_results = 3;
2715 :
2716 : // Initialize result_set and do the search.
2717 31530 : std::vector<size_t> ret_index(num_results);
2718 15765 : std::vector<Real> out_dist_sqr(num_results);
2719 15765 : nanoflann::KNNResultSet<Real> result_set(num_results);
2720 15765 : result_set.init(&ret_index[0], &out_dist_sqr[0]);
2721 15765 : kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParameters());
2722 :
2723 : // If this flag gets set in the loop below, we can break out of the outer r-loop as well.
2724 15765 : bool projection_succeeded = false;
2725 :
2726 : // Once we've rejected a candidate for a given
2727 : // primary_node, there's no reason to check it
2728 : // again.
2729 15765 : std::set<const Elem *> rejected_secondary_elem_candidates;
2730 :
2731 : // Loop over the closest nodes, check whether the secondary node successfully projects into
2732 : // either of the closest neighbors, stop when the projection succeeds.
2733 26091 : for (MooseIndex(result_set) r = 0; r < result_set.size(); ++r)
2734 : {
2735 : // Verify that the squared distance we compute is the same as nanoflann's
2736 : mooseAssert(abs((_mesh.point(ret_index[r]) - *primary_node).norm_sq() - out_dist_sqr[r]) <=
2737 : TOLERANCE,
2738 : "Lower-dimensional element squared distance verification failed.");
2739 :
2740 : // Get a reference to the vector of lower dimensional elements from the
2741 : // nodes_to_secondary_elem_map.
2742 : const std::vector<const Elem *> & secondary_elem_candidates =
2743 22649 : _nodes_to_secondary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));
2744 :
2745 : // Print the Elems connected to this node on the secondary mesh side.
2746 44255 : for (MooseIndex(secondary_elem_candidates) e = 0; e < secondary_elem_candidates.size(); ++e)
2747 : {
2748 33929 : const Elem * secondary_elem_candidate = secondary_elem_candidates[e];
2749 :
2750 : // If we've already rejected this candidate, we don't need to check it again.
2751 33929 : if (rejected_secondary_elem_candidates.count(secondary_elem_candidate))
2752 6884 : continue;
2753 :
2754 27045 : std::vector<Point> nodal_normals(secondary_elem_candidate->n_nodes());
2755 82010 : for (const auto n : make_range(secondary_elem_candidate->n_nodes()))
2756 109930 : nodal_normals[n] =
2757 54965 : _secondary_node_to_nodal_normal.at(secondary_elem_candidate->node_ptr(n));
2758 :
2759 : // Use equation 2.4.6 from Bin Yang's dissertation to try and solve for
2760 : // the position on the secondary element where this primary came from. This
2761 : // requires a Newton iteration in general.
2762 27045 : DualNumber<Real> xi1_dn{0, 1}; // initial guess
2763 27045 : auto && order = secondary_elem_candidate->default_order();
2764 27045 : unsigned int current_iterate = 0, max_iterates = 10;
2765 :
2766 27045 : VectorValue<DualNumber<Real>> normals(0);
2767 :
2768 : // Newton iteration loop - this to converge in 1 iteration when it
2769 : // succeeds, and possibly two iterations when it converges to a
2770 : // xi outside the reference element. I don't know any reason why it should
2771 : // only take 1 iteration -- the Jacobian is not constant in general...
2772 : do
2773 : {
2774 53576 : VectorValue<DualNumber<Real>> x1(0);
2775 162303 : for (MooseIndex(secondary_elem_candidate->n_nodes()) n = 0;
2776 162303 : n < secondary_elem_candidate->n_nodes();
2777 : ++n)
2778 : {
2779 108727 : const auto phi = Moose::fe_lagrange_1D_shape(order, n, xi1_dn);
2780 108727 : x1 += phi * secondary_elem_candidate->point(n);
2781 108727 : normals += phi * nodal_normals[n];
2782 108727 : }
2783 :
2784 53576 : const auto u = x1 - (*primary_node);
2785 :
2786 53576 : const auto F = u(0) * normals(1) - u(1) * normals(0);
2787 :
2788 53576 : if (abs(F) < _newton_tolerance)
2789 27045 : break;
2790 :
2791 : // Unlike for projection of nodal normals onto primary surfaces, we should never have a
2792 : // case where the nodal normal is completely orthogonal to the secondary surface, so we
2793 : // do not have to guard against F.derivatives() == 0 here
2794 26531 : Real dxi1 = -F.value() / F.derivatives();
2795 :
2796 26531 : xi1_dn += dxi1;
2797 :
2798 26531 : normals = 0;
2799 134197 : } while (++current_iterate < max_iterates);
2800 :
2801 27045 : Real xi1 = xi1_dn.value();
2802 :
2803 : // Check for convergence to a valid solution... The last condition checks for obliqueness
2804 : // of the projection
2805 39368 : if ((current_iterate < max_iterates) && (abs(xi1) <= 1. + _xi_tolerance) &&
2806 12323 : (abs((primary_side_elem->point(0) - primary_side_elem->point(1)).unit() *
2807 39368 : MetaPhysicL::raw_value(normals).unit()) <
2808 12323 : std::cos(_minimum_projection_angle * libMesh::pi / 180.0)))
2809 : {
2810 12323 : if (abs(abs(xi1) - 1.) < _xi_tolerance)
2811 : {
2812 : // Special case: xi1=+/-1.
2813 : // It is unlikely that we get here, because this primary node should already
2814 : // have been mapped during the project_secondary_nodes() routine, but
2815 : // there is still a chance since the tolerances are applied to
2816 : // the xi coordinate and that value may be different on a primary element and a
2817 : // secondary element since they may have different sizes. It's also possible that we
2818 : // may reach this point if the solve has yielded a non-physical configuration such as
2819 : // one block being pushed way out into space
2820 0 : const Node & secondary_node = (xi1 < 0) ? secondary_elem_candidate->node_ref(0)
2821 0 : : secondary_elem_candidate->node_ref(1);
2822 0 : bool created_mortar_segment = false;
2823 :
2824 : // If we have failed to project this secondary node, let's try again now
2825 0 : if (_failed_secondary_node_projections.count(secondary_node.id()))
2826 0 : created_mortar_segment = processAlignedNodes(secondary_node,
2827 : *primary_node,
2828 : nullptr,
2829 : &primary_node_neighbors,
2830 0 : MetaPhysicL::raw_value(normals),
2831 : *secondary_elem_candidate,
2832 : rejected_secondary_elem_candidates);
2833 : else
2834 0 : rejected_secondary_elem_candidates.insert(secondary_elem_candidate);
2835 :
2836 0 : if (!created_mortar_segment)
2837 : // We used to throw an exception in this scope but now that we support processing
2838 : // aligned nodes within this primary node projection method, I don't see any harm in
2839 : // simply rejecting the secondary element candidate in the case of failure and
2840 : // continuing just as we do when projecting secondary nodes
2841 0 : continue;
2842 : }
2843 : else // somewhere in the middle of the Elem
2844 : {
2845 : // Add entry to primary_node_and_elem_to_xi1_secondary_elem
2846 : //
2847 : // Note: we originally duplicated the map values for the keys (node, left_neighbor)
2848 : // and (node, right_neighbor) but I don't think that should be necessary. Instead we
2849 : // just do it for neighbor 0, but really maybe we don't even need to do that since
2850 : // we can always look up the neighbors later given the Node... keeping it like this
2851 : // helps to maintain the "symmetry" of the two containers.
2852 12323 : const Elem * neigh = primary_node_neighbors[0];
2853 36969 : for (MooseIndex(neigh->n_vertices()) nid = 0; nid < neigh->n_vertices(); ++nid)
2854 : {
2855 24646 : const Node * neigh_node = neigh->node_ptr(nid);
2856 24646 : if (primary_node == neigh_node)
2857 : {
2858 12323 : auto key = std::make_tuple(neigh_node->id(), neigh_node, neigh);
2859 12323 : auto val = std::make_pair(xi1, secondary_elem_candidate);
2860 12323 : _primary_node_and_elem_to_xi1_secondary_elem.emplace(key, val);
2861 : }
2862 : }
2863 : }
2864 :
2865 12323 : projection_succeeded = true;
2866 12323 : break; // out of e-loop
2867 : }
2868 : else
2869 : {
2870 : // The current primary_point is not in this Elem, so keep track of the rejects.
2871 14722 : rejected_secondary_elem_candidates.insert(secondary_elem_candidate);
2872 : }
2873 51691 : } // end e-loop over candidate elems
2874 :
2875 22649 : if (projection_succeeded)
2876 12323 : break; // out of r-loop
2877 : } // r-loop
2878 :
2879 15765 : if (!projection_succeeded && _debug)
2880 : {
2881 0 : _console << "\nFailed to find point from which primary node "
2882 0 : << static_cast<const Point &>(*primary_node) << " was projected." << std::endl
2883 0 : << std::endl;
2884 : }
2885 15765 : } // loop over side nodes
2886 4263 : } // end loop over elements for finding where primary points would have projected from.
2887 4263 : }
2888 :
2889 : std::vector<AutomaticMortarGeneration::MortarFilterIter>
2890 595 : AutomaticMortarGeneration::secondariesToMortarSegments(const Node & node) const
2891 : {
2892 595 : auto secondary_it = _nodes_to_secondary_elem_map.find(node.id());
2893 595 : if (secondary_it == _nodes_to_secondary_elem_map.end())
2894 0 : return {};
2895 :
2896 595 : const auto & secondary_elems = secondary_it->second;
2897 595 : std::vector<MortarFilterIter> ret;
2898 595 : ret.reserve(secondary_elems.size());
2899 :
2900 1444 : for (const auto i : index_range(secondary_elems))
2901 : {
2902 849 : auto * const secondary_elem = secondary_elems[i];
2903 849 : auto msm_it = _secondary_elems_to_mortar_segments.find(secondary_elem->id());
2904 849 : if (msm_it == _secondary_elems_to_mortar_segments.end())
2905 : // We may have removed this element key from this map
2906 0 : continue;
2907 :
2908 : mooseAssert(secondary_elem->active(),
2909 : "We loop over active elements when building the mortar segment mesh, so we golly "
2910 : "well hope this is active.");
2911 : mooseAssert(!msm_it->second.empty(),
2912 : "We should have removed all secondaries from this map if they do not have any "
2913 : "mortar segments associated with them.");
2914 849 : ret.push_back(msm_it);
2915 : }
2916 :
2917 595 : return ret;
2918 595 : }
|
