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