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