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