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