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