https://mooseframework.inl.gov
AutomaticMortarGeneration.C
Go to the documentation of this file.
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 
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;
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 
66 {
67 public:
69  {
70  auto params = Output::validParams();
71  params.addPrivateParam<AutomaticMortarGeneration *>("_amg", nullptr);
72  params.addPrivateParam<MooseApp *>("_moose_app", nullptr);
73  params.set<std::string>("_type") = "MortarNodalGeometryOutput";
74  return params;
75  };
76 
78  : Output(params), _amg(*getCheckedPointerParam<AutomaticMortarGeneration *>("_amg"))
79  {
80  }
81 
82  void output() override
83  {
84  // Must call compute_nodal_geometry first!
85  if (_amg._secondary_node_to_nodal_normal.empty() ||
86  _amg._secondary_node_to_hh_nodal_tangents.empty())
87  mooseError("No entries found in the secondary node -> nodal geometry map.");
88 
89  auto & problem = _app.feProblem();
90  auto & subproblem = _amg._on_displaced
91  ? static_cast<SubProblem &>(*problem.getDisplacedProblem())
92  : static_cast<SubProblem &>(problem);
93  auto & nodal_normals_es = subproblem.es();
94 
95  const std::string nodal_normals_sys_name = "nodal_normals";
96 
97  if (!_nodal_normals_system)
98  {
99  for (const auto s : make_range(nodal_normals_es.n_systems()))
100  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  return;
104 
105  _nodal_normals_system =
106  &nodal_normals_es.template add_system<ExplicitSystem>(nodal_normals_sys_name);
107  _nnx_var_num = _nodal_normals_system->add_variable("nodal_normal_x", FEType(FIRST, LAGRANGE)),
108  _nny_var_num = _nodal_normals_system->add_variable("nodal_normal_y", FEType(FIRST, LAGRANGE));
109  _nnz_var_num = _nodal_normals_system->add_variable("nodal_normal_z", FEType(FIRST, LAGRANGE));
110 
111  _t1x_var_num =
112  _nodal_normals_system->add_variable("nodal_tangent_1_x", FEType(FIRST, LAGRANGE)),
113  _t1y_var_num =
114  _nodal_normals_system->add_variable("nodal_tangent_1_y", FEType(FIRST, LAGRANGE));
115  _t1z_var_num =
116  _nodal_normals_system->add_variable("nodal_tangent_1_z", FEType(FIRST, LAGRANGE));
117 
118  _t2x_var_num =
119  _nodal_normals_system->add_variable("nodal_tangent_2_x", FEType(FIRST, LAGRANGE)),
120  _t2y_var_num =
121  _nodal_normals_system->add_variable("nodal_tangent_2_y", FEType(FIRST, LAGRANGE));
122  _t2z_var_num =
123  _nodal_normals_system->add_variable("nodal_tangent_2_z", FEType(FIRST, LAGRANGE));
124  nodal_normals_es.reinit();
125  }
126 
127  const DofMap & dof_map = _nodal_normals_system->get_dof_map();
128  std::vector<dof_id_type> dof_indices_nnx, dof_indices_nny, dof_indices_nnz;
129  std::vector<dof_id_type> dof_indices_t1x, dof_indices_t1y, dof_indices_t1z;
130  std::vector<dof_id_type> dof_indices_t2x, dof_indices_t2y, dof_indices_t2z;
131 
132  for (MeshBase::const_element_iterator el = _amg._mesh.elements_begin(),
133  end_el = _amg._mesh.elements_end();
134  el != end_el;
135  ++el)
136  {
137  const Elem * elem = *el;
138 
139  // Get the nodal dofs for this Elem.
140  dof_map.dof_indices(elem, dof_indices_nnx, _nnx_var_num);
141  dof_map.dof_indices(elem, dof_indices_nny, _nny_var_num);
142  dof_map.dof_indices(elem, dof_indices_nnz, _nnz_var_num);
143 
144  dof_map.dof_indices(elem, dof_indices_t1x, _t1x_var_num);
145  dof_map.dof_indices(elem, dof_indices_t1y, _t1y_var_num);
146  dof_map.dof_indices(elem, dof_indices_t1z, _t1z_var_num);
147 
148  dof_map.dof_indices(elem, dof_indices_t2x, _t2x_var_num);
149  dof_map.dof_indices(elem, dof_indices_t2y, _t2y_var_num);
150  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  for (MooseIndex(elem->n_vertices()) n = 0; n < elem->n_vertices(); ++n)
157  {
158  auto it = _amg._secondary_node_to_nodal_normal.find(elem->node_ptr(n));
159  if (it != _amg._secondary_node_to_nodal_normal.end())
160  {
161  _nodal_normals_system->solution->set(dof_indices_nnx[n], it->second(0));
162  _nodal_normals_system->solution->set(dof_indices_nny[n], it->second(1));
163  _nodal_normals_system->solution->set(dof_indices_nnz[n], it->second(2));
164  }
165 
166  auto it_tangent = _amg._secondary_node_to_hh_nodal_tangents.find(elem->node_ptr(n));
167  if (it_tangent != _amg._secondary_node_to_hh_nodal_tangents.end())
168  {
169  _nodal_normals_system->solution->set(dof_indices_t1x[n], it_tangent->second[0](0));
170  _nodal_normals_system->solution->set(dof_indices_t1y[n], it_tangent->second[0](1));
171  _nodal_normals_system->solution->set(dof_indices_t1z[n], it_tangent->second[0](2));
172 
173  _nodal_normals_system->solution->set(dof_indices_t2x[n], it_tangent->second[1](0));
174  _nodal_normals_system->solution->set(dof_indices_t2y[n], it_tangent->second[1](1));
175  _nodal_normals_system->solution->set(dof_indices_t2z[n], it_tangent->second[1](2));
176  }
177 
178  } // end loop over nodes
179  } // end loop over elems
180 
181  // Finish assembly.
182  _nodal_normals_system->solution->close();
183 
184  std::set<std::string> sys_names = {nodal_normals_sys_name};
185 
186  // Write the nodal normals to file
187  ExodusII_IO nodal_normals_writer(_amg._mesh);
188 
189  // Default to non-HDF5 output for wider compatibility
190  nodal_normals_writer.set_hdf5_writing(false);
191 
192  nodal_normals_writer.write_equation_systems(
193  "nodal_geometry_only.e", nodal_normals_es, &sys_names);
194  }
195 
196 private:
199 
201 
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;
215 };
216 
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  const Real minimum_projection_angle)
227  : ConsoleStreamInterface(app),
228  _app(app),
229  _mesh(mesh_in),
230  _debug(debug),
231  _on_displaced(on_displaced),
232  _periodic(periodic),
233  _distributed(!_mesh.is_replicated()),
234  _correct_edge_dropping(correct_edge_dropping),
235  _minimum_projection_angle(minimum_projection_angle)
236 {
237  _primary_secondary_boundary_id_pairs.push_back(boundary_key);
238  _primary_requested_boundary_ids.insert(boundary_key.first);
239  _secondary_requested_boundary_ids.insert(boundary_key.second);
240  _primary_secondary_subdomain_id_pairs.push_back(subdomain_key);
241  _primary_boundary_subdomain_ids.insert(subdomain_key.first);
242  _secondary_boundary_subdomain_ids.insert(subdomain_key.second);
243 
244  if (_distributed)
245  _mortar_segment_mesh = std::make_unique<DistributedMesh>(_mesh.comm());
246  else
247  _mortar_segment_mesh = std::make_unique<ReplicatedMesh>(_mesh.comm());
248 }
249 
250 void
252 {
253  if (!_debug)
254  return;
255 
256  _output_params = std::make_unique<InputParameters>(MortarNodalGeometryOutput::validParams());
257  _output_params->set<AutomaticMortarGeneration *>("_amg") = this;
258  _output_params->set<FEProblemBase *>("_fe_problem_base") = &_app.feProblem();
259  _output_params->set<MooseApp *>("_moose_app") = &_app;
260  _output_params->set<std::string>("_object_name") =
261  "mortar_nodal_geometry_" +
262  std::to_string(_primary_secondary_boundary_id_pairs.front().first) +
263  std::to_string(_primary_secondary_boundary_id_pairs.front().second) + "_" +
264  (_on_displaced ? "displaced" : "undisplaced");
265  _output_params->finalize("MortarNodalGeometryOutput");
266  _app.getOutputWarehouse().addOutput(std::make_shared<MortarNodalGeometryOutput>(*_output_params));
267 }
268 
269 void
271 {
272  _mortar_segment_mesh->clear();
277  _msm_elem_to_info.clear();
278  _lower_elem_to_side_id.clear();
284  _secondary_ip_sub_ids.clear();
285  _primary_ip_sub_ids.clear();
286 }
287 
288 void
290 {
292  mooseError(
293  "Must specify secondary and primary boundary ids before building node-to-elem maps.");
294 
295  // Construct nodes_to_secondary_elem_map
296  for (const auto & secondary_elem :
297  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  if (!this->_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
301  continue;
302 
303  for (const auto & nd : secondary_elem->node_ref_range())
304  {
305  std::vector<const Elem *> & vec = _nodes_to_secondary_elem_map[nd.id()];
306  vec.push_back(secondary_elem);
307  }
308  }
309 
310  // Construct nodes_to_primary_elem_map
311  for (const auto & primary_elem :
312  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  if (!this->_primary_boundary_subdomain_ids.count(primary_elem->subdomain_id()))
316  continue;
317 
318  for (const auto & nd : primary_elem->node_ref_range())
319  {
320  std::vector<const Elem *> & vec = _nodes_to_primary_elem_map[nd.id()];
321  vec.push_back(primary_elem);
322  }
323  }
324 }
325 
326 std::vector<Point>
328 {
329  std::vector<Point> nodal_normals(secondary_elem.n_nodes());
330  for (const auto n : make_range(secondary_elem.n_nodes()))
331  nodal_normals[n] = _secondary_node_to_nodal_normal.at(secondary_elem.node_ptr(n));
332 
333  return nodal_normals;
334 }
335 
336 const Elem *
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  return _secondary_element_to_secondary_lowerd_element.at(secondary_elem_id);
344 }
345 
346 std::map<unsigned int, unsigned int>
348 {
349  std::map<unsigned int, unsigned int> secondary_ip_i_to_lower_secondary_i;
350  const Elem * const secondary_ip = lower_secondary_elem.interior_parent();
351  mooseAssert(secondary_ip, "This should be non-null");
352 
353  for (const auto i : make_range(lower_secondary_elem.n_nodes()))
354  {
355  const auto & nd = lower_secondary_elem.node_ref(i);
356  secondary_ip_i_to_lower_secondary_i[secondary_ip->get_node_index(&nd)] = i;
357  }
358 
359  return secondary_ip_i_to_lower_secondary_i;
360 }
361 
362 std::map<unsigned int, unsigned int>
364  const Elem & lower_primary_elem,
365  const Elem & primary_elem,
366  const Elem & /*lower_secondary_elem*/) const
367 {
368  std::map<unsigned int, unsigned int> primary_ip_i_to_lower_primary_i;
369 
370  for (const auto i : make_range(lower_primary_elem.n_nodes()))
371  {
372  const auto & nd = lower_primary_elem.node_ref(i);
373  primary_ip_i_to_lower_primary_i[primary_elem.get_node_index(&nd)] = i;
374  }
375 
376  return primary_ip_i_to_lower_primary_i;
377 }
378 
379 std::array<MooseUtils::SemidynamicVector<Point, 9>, 2>
381 {
382  // MetaPhysicL will check if we ran out of allocated space.
383  MooseUtils::SemidynamicVector<Point, 9> nodal_tangents_one(0);
384  MooseUtils::SemidynamicVector<Point, 9> nodal_tangents_two(0);
385 
386  for (const auto n : make_range(secondary_elem.n_nodes()))
387  {
388  const auto & tangent_vectors =
389  libmesh_map_find(_secondary_node_to_hh_nodal_tangents, secondary_elem.node_ptr(n));
390  nodal_tangents_one.push_back(tangent_vectors[0]);
391  nodal_tangents_two.push_back(tangent_vectors[1]);
392  }
393 
394  return {{nodal_tangents_one, nodal_tangents_two}};
395 }
396 
397 std::vector<Point>
399  const std::vector<Real> & oned_xi1_pts) const
400 {
401  std::vector<Point> xi1_pts(oned_xi1_pts.size());
402  for (const auto qp : index_range(oned_xi1_pts))
403  xi1_pts[qp] = oned_xi1_pts[qp];
404 
405  return getNormals(secondary_elem, xi1_pts);
406 }
407 
408 std::vector<Point>
410  const std::vector<Point> & xi1_pts) const
411 {
412  const auto mortar_dim = _mesh.mesh_dimension() - 1;
413  const auto num_qps = xi1_pts.size();
414  const auto nodal_normals = getNodalNormals(secondary_elem);
415  std::vector<Point> normals(num_qps);
416 
417  for (const auto n : make_range(secondary_elem.n_nodes()))
418  for (const auto qp : make_range(num_qps))
419  {
420  const auto phi =
421  (mortar_dim == 1)
422  ? Moose::fe_lagrange_1D_shape(secondary_elem.default_order(), n, xi1_pts[qp](0))
423  : Moose::fe_lagrange_2D_shape(secondary_elem.type(),
424  secondary_elem.default_order(),
425  n,
426  static_cast<const TypeVector<Real> &>(xi1_pts[qp]));
427  normals[qp] += phi * nodal_normals[n];
428  }
429 
430  if (_periodic)
431  for (auto & normal : normals)
432  normal *= -1;
433 
434  return normals;
435 }
436 
437 void
439 {
440  dof_id_type local_id_index = 0;
441  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  for (const auto & pr : _primary_secondary_boundary_id_pairs)
448  {
449  const auto primary_bnd_id = pr.first;
450  const auto secondary_bnd_id = pr.second;
451  const auto num_primary_nodes =
452  std::distance(_mesh.bid_nodes_begin(primary_bnd_id), _mesh.bid_nodes_end(primary_bnd_id));
453  const auto num_secondary_nodes = std::distance(_mesh.bid_nodes_begin(secondary_bnd_id),
454  _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  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  for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
467  end_el = _mesh.active_elements_end();
468  el != end_el;
469  ++el)
470  {
471  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  if (!this->_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
475  continue;
476 
477  std::vector<Node *> new_nodes;
478  for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
479  {
480  new_nodes.push_back(_mortar_segment_mesh->add_point(
481  secondary_elem->point(n), secondary_elem->node_id(n), secondary_elem->processor_id()));
482  Node * const new_node = new_nodes.back();
483  new_node->set_unique_id(new_node->id() + node_unique_id_offset);
484  }
485 
486  std::unique_ptr<Elem> new_elem;
487  if (secondary_elem->default_order() == SECOND)
488  new_elem = std::make_unique<Edge3>();
489  else
490  new_elem = std::make_unique<Edge2>();
491 
492  new_elem->processor_id() = secondary_elem->processor_id();
493  new_elem->subdomain_id() = secondary_elem->subdomain_id();
494  new_elem->set_id(local_id_index++);
495  new_elem->set_unique_id(new_elem->id());
496 
497  for (MooseIndex(new_elem->n_nodes()) n = 0; n < new_elem->n_nodes(); ++n)
498  new_elem->set_node(n, new_nodes[n]);
499 
500  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  MortarSegmentInfo msinfo;
504  msinfo.xi1_a = -1;
505  msinfo.xi1_b = +1;
506  msinfo.secondary_elem = secondary_elem;
507 
508  auto new_container_it0 = _secondary_node_and_elem_to_xi2_primary_elem.find(
509  std::make_pair(secondary_elem->node_ptr(0), secondary_elem)),
510  new_container_it1 = _secondary_node_and_elem_to_xi2_primary_elem.find(
511  std::make_pair(secondary_elem->node_ptr(1), secondary_elem));
512 
513  bool new_container_node0_found =
514  (new_container_it0 != _secondary_node_and_elem_to_xi2_primary_elem.end()),
515  new_container_node1_found =
516  (new_container_it1 != _secondary_node_and_elem_to_xi2_primary_elem.end());
517 
518  const Elem * node0_primary_candidate = nullptr;
519  const Elem * node1_primary_candidate = nullptr;
520 
521  if (new_container_node0_found)
522  {
523  const auto & xi2_primary_elem_pair = new_container_it0->second;
524  msinfo.xi2_a = xi2_primary_elem_pair.first;
525  node0_primary_candidate = xi2_primary_elem_pair.second;
526  }
527 
528  if (new_container_node1_found)
529  {
530  const auto & xi2_primary_elem_pair = new_container_it1->second;
531  msinfo.xi2_b = xi2_primary_elem_pair.first;
532  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  if (node0_primary_candidate == node1_primary_candidate)
540  msinfo.primary_elem = node0_primary_candidate;
541 
542  // Associate this MSM elem with the MortarSegmentInfo.
543  _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  _secondary_elems_to_mortar_segments[secondary_elem->id()].insert(new_elem_ptr);
548  }
549 
550  // 2.) Insert new nodes from primary side and split mortar segments as necessary.
551  for (const auto & pr : _primary_node_and_elem_to_xi1_secondary_elem)
552  {
553  auto key = pr.first;
554  auto val = pr.second;
555 
556  const Node * primary_node = std::get<1>(key);
557  Real xi1 = val.first;
558  const Elem * secondary_elem = val.second;
559 
560  // If this is an aligned node, we don't need to do anything.
561  if (std::abs(std::abs(xi1) - 1.) < _xi_tolerance)
562  continue;
563 
564  auto && order = secondary_elem->default_order();
565 
566  // Determine physical location of new point to be inserted.
567  Point new_pt(0);
568  for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
569  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  auto & mortar_segment_set = _secondary_elems_to_mortar_segments[secondary_elem->id()];
573  Elem * current_mortar_segment = nullptr;
574  MortarSegmentInfo * info = nullptr;
575 
576  for (const auto & mortar_segment_candidate : mortar_segment_set)
577  {
578  try
579  {
580  info = &_msm_elem_to_info.at(mortar_segment_candidate);
581  }
582  catch (std::out_of_range &)
583  {
584  mooseError("MortarSegmentInfo not found for the mortar segment candidate");
585  }
586  if (info->xi1_a <= xi1 && xi1 <= info->xi1_b)
587  {
588  current_mortar_segment = mortar_segment_candidate;
589  break;
590  }
591  }
592 
593  // Make sure we found one.
594  if (current_mortar_segment == nullptr)
595  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  if (info->xi1_a == xi1 || xi1 == info->xi1_b)
603  continue;
604 
605  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  _mortar_segment_mesh->add_point(new_pt, new_id, secondary_elem->processor_id());
610  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  const Point normal = getNormals(*secondary_elem, std::vector<Real>({xi1}))[0];
616 
617  // Get the set of primary_node neighbors.
618  if (this->_nodes_to_primary_elem_map.find(primary_node->id()) ==
619  this->_nodes_to_primary_elem_map.end())
620  mooseError("We should already have built this primary node to elem pair!");
621  const std::vector<const Elem *> & primary_node_neighbors =
622  this->_nodes_to_primary_elem_map[primary_node->id()];
623 
624  // Sanity check
625  if (primary_node_neighbors.size() == 0 || primary_node_neighbors.size() > 2)
626  mooseError("We must have either 1 or 2 primary side nodal neighbors, but we had ",
627  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  const Elem * left_primary_elem = primary_node_neighbors[0];
635  const Elem * right_primary_elem =
636  (primary_node_neighbors.size() == 2) ? primary_node_neighbors[1] : nullptr;
637 
639 
640  // Storage for z-component of cross products for determining
641  // orientation.
642  std::array<Real, 2> secondary_node_cps;
643  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  for (unsigned int nid = 0; nid < 2; ++nid)
647  secondary_node_cps[nid] = normal.cross(secondary_elem->point(nid) - new_pt)(2);
648 
649  for (MooseIndex(primary_node_neighbors) mnn = 0; mnn < primary_node_neighbors.size(); ++mnn)
650  {
651  const Elem * primary_neigh = primary_node_neighbors[mnn];
652  Point opposite = (primary_neigh->node_ptr(0) == primary_node) ? primary_neigh->point(1)
653  : primary_neigh->point(0);
654  Point cp = normal.cross(opposite - new_pt);
655  primary_node_cps[mnn] = cp(2);
656  }
657 
658  // We will verify that only 1 orientation is actually valid.
659  bool orientation1_valid = false, orientation2_valid = false;
660 
661  if (primary_node_neighbors.size() == 2)
662  {
663  // 2 primary neighbor case
664  orientation1_valid = (secondary_node_cps[0] * primary_node_cps[0] > 0.) &&
665  (secondary_node_cps[1] * primary_node_cps[1] > 0.);
666 
667  orientation2_valid = (secondary_node_cps[0] * primary_node_cps[1] > 0.) &&
668  (secondary_node_cps[1] * primary_node_cps[0] > 0.);
669  }
670  else if (primary_node_neighbors.size() == 1)
671  {
672  // 1 primary neighbor case
673  orientation1_valid = (secondary_node_cps[0] * primary_node_cps[0] > 0.);
674  orientation2_valid = (secondary_node_cps[1] * primary_node_cps[0] > 0.);
675  }
676  else
677  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  if (orientation1_valid && orientation2_valid)
684  throw MooseException(
685  "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  if (!orientation1_valid && !orientation2_valid)
692  {
693  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  continue;
701  }
702 
703  // Make an Elem on the left
704  std::unique_ptr<Elem> new_elem_left;
705  if (order == SECOND)
706  new_elem_left = std::make_unique<Edge3>();
707  else
708  new_elem_left = std::make_unique<Edge2>();
709 
710  new_elem_left->processor_id() = current_mortar_segment->processor_id();
711  new_elem_left->subdomain_id() = current_mortar_segment->subdomain_id();
712  new_elem_left->set_id(local_id_index++);
713  new_elem_left->set_unique_id(new_elem_left->id());
714  new_elem_left->set_node(0, current_mortar_segment->node_ptr(0));
715  new_elem_left->set_node(1, new_node);
716 
717  // Make an Elem on the right
718  std::unique_ptr<Elem> new_elem_right;
719  if (order == SECOND)
720  new_elem_right = std::make_unique<Edge3>();
721  else
722  new_elem_right = std::make_unique<Edge2>();
723 
724  new_elem_right->processor_id() = current_mortar_segment->processor_id();
725  new_elem_right->subdomain_id() = current_mortar_segment->subdomain_id();
726  new_elem_right->set_id(local_id_index++);
727  new_elem_right->set_unique_id(new_elem_right->id());
728  new_elem_right->set_node(0, new_node);
729  new_elem_right->set_node(1, current_mortar_segment->node_ptr(1));
730 
731  if (order == SECOND)
732  {
733  // left
734  Point left_interior_point(0);
735  Real left_interior_xi = (xi1 + info->xi1_a) / 2;
736 
737  // This is eta for the current mortar segment that we're splitting
738  Real current_left_interior_eta =
739  (2. * left_interior_xi - info->xi1_a - info->xi1_b) / (info->xi1_b - info->xi1_a);
740 
741  for (MooseIndex(current_mortar_segment->n_nodes()) n = 0;
742  n < current_mortar_segment->n_nodes();
743  ++n)
744  left_interior_point += Moose::fe_lagrange_1D_shape(order, n, current_left_interior_eta) *
745  current_mortar_segment->point(n);
746 
747  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  Node * const new_interior_node_left = _mortar_segment_mesh->add_point(
751  left_interior_point, new_interior_left_id, new_elem_left->processor_id());
752  new_elem_left->set_node(2, new_interior_node_left);
753  new_interior_node_left->set_unique_id(new_interior_left_id + node_unique_id_offset);
754 
755  // right
756  Point right_interior_point(0);
757  Real right_interior_xi = (xi1 + info->xi1_b) / 2;
758  // This is eta for the current mortar segment that we're splitting
759  Real current_right_interior_eta =
760  (2. * right_interior_xi - info->xi1_a - info->xi1_b) / (info->xi1_b - info->xi1_a);
761 
762  for (MooseIndex(current_mortar_segment->n_nodes()) n = 0;
763  n < current_mortar_segment->n_nodes();
764  ++n)
765  right_interior_point += Moose::fe_lagrange_1D_shape(order, n, current_right_interior_eta) *
766  current_mortar_segment->point(n);
767 
768  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  Node * const new_interior_node_right = _mortar_segment_mesh->add_point(
772  right_interior_point, new_interior_id_right, new_elem_right->processor_id());
773  new_elem_right->set_node(2, new_interior_node_right);
774  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  if (orientation2_valid)
779  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  if (left_primary_elem)
784  left_xi2 = (primary_node == left_primary_elem->node_ptr(0)) ? -1 : +1;
785  if (right_primary_elem)
786  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  auto msm_it = _msm_elem_to_info.find(current_mortar_segment);
794  if (msm_it == _msm_elem_to_info.end())
795  mooseError("MortarSegmentInfo not found for current_mortar_segment.");
796  MortarSegmentInfo current_msinfo = msm_it->second;
797 
798  // add_left
799  {
800  Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem_left.release());
801 
802  // Create new MortarSegmentInfo objects for new_elem_left
803  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  new_msinfo_left.xi1_a = current_msinfo.xi1_a;
809  new_msinfo_left.xi2_a = current_msinfo.xi2_a;
810  new_msinfo_left.secondary_elem = secondary_elem;
811  new_msinfo_left.xi1_b = xi1;
812  new_msinfo_left.xi2_b = left_xi2;
813  new_msinfo_left.primary_elem = left_primary_elem;
814 
815  // Add new msinfo objects to the map.
816  _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  mortar_segment_set.insert(msm_new_elem);
821  }
822 
823  // add_right
824  {
825  Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem_right.release());
826 
827  // Create new MortarSegmentInfo objects for new_elem_right
828  MortarSegmentInfo new_msinfo_right;
829 
830  new_msinfo_right.xi1_b = current_msinfo.xi1_b;
831  new_msinfo_right.xi2_b = current_msinfo.xi2_b;
832  new_msinfo_right.secondary_elem = secondary_elem;
833  new_msinfo_right.xi1_a = xi1;
834  new_msinfo_right.xi2_a = right_xi2;
835  new_msinfo_right.primary_elem = right_primary_elem;
836 
837  _msm_elem_to_info.emplace(msm_new_elem, new_msinfo_right);
838 
839  mortar_segment_set.insert(msm_new_elem);
840  }
841 
842  // Erase the MortarSegmentInfo object for current_mortar_segment from the map.
843  _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  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  _mortar_segment_mesh->delete_elem(current_mortar_segment);
852  }
853 
854  // Remove all MSM elements without a primary contribution
860  for (auto msm_elem : _mortar_segment_mesh->active_element_ptr_range())
861  {
862  MortarSegmentInfo & msinfo = libmesh_map_find(_msm_elem_to_info, msm_elem);
863  Elem * primary_elem = const_cast<Elem *>(msinfo.primary_elem);
864  if (primary_elem == nullptr || std::abs(msinfo.xi2_a) > 1.0 + TOLERANCE ||
865  std::abs(msinfo.xi2_b) > 1.0 + TOLERANCE)
866  {
867  // Erase from secondary to msms map
868  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  auto & msm_set = it->second;
872  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  if (msm_set.empty())
881 
882  // Erase msinfo
883  _msm_elem_to_info.erase(msm_elem);
884 
885  // Remove element from mortar segment mesh
886  _mortar_segment_mesh->delete_elem(msm_elem);
887  }
888  else
889  {
892  }
893  }
894 
895  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  for (const auto & element : _mortar_segment_mesh->element_ptr_range())
900  for (auto & n : element->node_ref_range())
901  msm_connected_nodes.insert(&n);
902 
903  for (const auto & node : _mortar_segment_mesh->node_ptr_range())
904  if (!msm_connected_nodes.count(node))
905  _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  _mortar_segment_mesh->cache_elem_data();
919 
920  // (Optionally) Write the mortar segment mesh to file for inspection
921  if (_debug)
922  {
923  ExodusII_IO mortar_segment_mesh_writer(*_mortar_segment_mesh);
924 
925  // Default to non-HDF5 output for wider compatibility
926  mortar_segment_mesh_writer.set_hdf5_writing(false);
927 
928  mortar_segment_mesh_writer.write("mortar_segment_mesh.e");
929  }
930 
932 }
933 
934 void
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  auto secondary_sub_elem = _mortar_segment_mesh->add_elem_integer("secondary_sub_elem");
940  auto primary_sub_elem = _mortar_segment_mesh->add_elem_integer("primary_sub_elem");
941 
942  dof_id_type local_id_index = 0;
943 
944  // Loop through mortar secondary and primary pairs to create mortar segment mesh between each
945  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
946  {
947  const auto primary_subd_id = pr.first;
948  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  NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, primary_subd_id);
952  subdomain_kd_tree_t kd_tree(
953  3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
954 
955  // Construct the KD tree.
956  kd_tree.buildIndex();
957 
958  // Define expression for getting sub-elements nodes (for sub-dividing secondary elements)
959  auto get_sub_elem_nodes = [](const ElemType type,
960  const unsigned int sub_elem) -> std::vector<unsigned int>
961  {
962  switch (type)
963  {
964  case TRI3:
965  return {{0, 1, 2}};
966  case QUAD4:
967  return {{0, 1, 2, 3}};
968  case TRI6:
969  case TRI7:
970  switch (sub_elem)
971  {
972  case 0:
973  return {{0, 3, 5}};
974  case 1:
975  return {{3, 4, 5}};
976  case 2:
977  return {{3, 1, 4}};
978  case 3:
979  return {{5, 4, 2}};
980  default:
981  mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
982  }
983  case QUAD8:
984  switch (sub_elem)
985  {
986  case 0:
987  return {{0, 4, 7}};
988  case 1:
989  return {{4, 1, 5}};
990  case 2:
991  return {{5, 2, 6}};
992  case 3:
993  return {{7, 6, 3}};
994  case 4:
995  return {{4, 5, 6, 7}};
996  default:
997  mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
998  }
999  case QUAD9:
1000  switch (sub_elem)
1001  {
1002  case 0:
1003  return {{0, 4, 8, 7}};
1004  case 1:
1005  return {{4, 1, 5, 8}};
1006  case 2:
1007  return {{8, 5, 2, 6}};
1008  case 3:
1009  return {{7, 8, 6, 3}};
1010  default:
1011  mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
1012  }
1013  default:
1014  mooseError("get_sub_elem_inds: Face element type: ",
1015  libMesh::Utility::enum_to_string<ElemType>(type),
1016  " invalid for 3D mortar");
1017  }
1018  };
1019 
1023  for (MeshBase::const_element_iterator el = _mesh.active_local_elements_begin(),
1024  end_el = _mesh.active_local_elements_end();
1025  el != end_el;
1026  ++el)
1027  {
1028  const Elem * secondary_side_elem = *el;
1029 
1030  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  if (secondary_side_elem->subdomain_id() != secondary_subd_id)
1034  continue;
1035 
1036  auto [secondary_elem_to_msm_map_it, insertion_happened] =
1037  _secondary_elems_to_mortar_segments.emplace(secondary_side_elem->id(),
1038  std::set<Elem *, CompareDofObjectsByID>{});
1039  libmesh_ignore(insertion_happened);
1040  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  secondary_side_elem->n_sub_elem());
1044  const auto nodal_normals = getNodalNormals(*secondary_side_elem);
1045 
1056  for (auto sel : make_range(secondary_side_elem->n_sub_elem()))
1057  {
1058  // Get indices of sub-element nodes in element
1059  auto sub_elem_nodes = get_sub_elem_nodes(secondary_side_elem->type(), sel);
1060 
1061  // Secondary sub-element center, normal, and nodes
1062  Point center;
1063  Point normal;
1064  std::vector<Point> nodes(sub_elem_nodes.size());
1065 
1066  // Loop through sub_element nodes, collect points and compute center and normal
1067  for (auto iv : make_range(sub_elem_nodes.size()))
1068  {
1069  const auto n = sub_elem_nodes[iv];
1070  nodes[iv] = secondary_side_elem->point(n);
1071  center += secondary_side_elem->point(n);
1072  normal += nodal_normals[n];
1073  }
1074  center /= sub_elem_nodes.size();
1075  normal = normal.unit();
1076 
1077  // Build and store linearized sub-elements for later use
1078  mortar_segment_helper[sel] = std::make_unique<MortarSegmentHelper>(nodes, center, normal);
1079  }
1080 
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  Point center_point;
1091  switch (secondary_side_elem->type())
1092  {
1093  case TRI3:
1094  case QUAD4:
1095  center_point = mortar_segment_helper[0]->center();
1096  query_pt = {{center_point(0), center_point(1), center_point(2)}};
1097  break;
1098  case TRI6:
1099  case TRI7:
1100  center_point = mortar_segment_helper[1]->center();
1101  query_pt = {{center_point(0), center_point(1), center_point(2)}};
1102  break;
1103  case QUAD8:
1104  center_point = mortar_segment_helper[4]->center();
1105  query_pt = {{center_point(0), center_point(1), center_point(2)}};
1106  break;
1107  case QUAD9:
1108  center_point = secondary_side_elem->point(8);
1109  query_pt = {{center_point(0), center_point(1), center_point(2)}};
1110  break;
1111  default:
1112  mooseError(
1113  "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  const std::size_t num_results = 3;
1120 
1121  // Initialize result_set and do the search.
1122  std::vector<size_t> ret_index(num_results);
1123  std::vector<Real> out_dist_sqr(num_results);
1124  nanoflann::KNNResultSet<Real> result_set(num_results);
1125  result_set.init(&ret_index[0], &out_dist_sqr[0]);
1126  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  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  bool primary_elem_found = false;
1134  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  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  this->_nodes_to_primary_elem_map.at(static_cast<dof_id_type>(ret_index[r]));
1147 
1148  // Uniquely add elems to candidate set
1149  for (auto elem : node_elems)
1150  primary_elem_candidates.insert(elem);
1151  }
1152 
1160  while (!primary_elem_candidates.empty())
1161  {
1162  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  if (processed_primary_elems.count(primary_elem_candidate))
1166  continue;
1167 
1168  // Initialize set of nodes used to construct mortar segment elements
1169  std::vector<Point> nodal_points;
1170 
1171  // Initialize map from mortar segment elements to nodes
1172  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  std::vector<std::pair<unsigned int, unsigned int>> sub_elem_map;
1176 
1181  for (auto p_el : make_range(primary_elem_candidate->n_sub_elem()))
1182  {
1183  // Get nodes of primary sub-elements
1184  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  std::vector<Point> primary_sub_elem(sub_elem_nodes.size());
1188  for (auto iv : make_range(sub_elem_nodes.size()))
1189  {
1190  const auto n = sub_elem_nodes[iv];
1191  primary_sub_elem[iv] = primary_elem_candidate->point(n);
1192  }
1193 
1194  // Loop through secondary sub-elements
1195  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  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  for (auto i = sub_elem_map.size(); i < elem_to_node_map.size(); ++i)
1209  sub_elem_map.push_back(std::make_pair(s_el, p_el));
1210  }
1211  }
1212 
1213  // Mark primary element as processed and remove from candidate list
1214  processed_primary_elems.insert(primary_elem_candidate);
1215  primary_elem_candidates.erase(primary_elem_candidate);
1216 
1217  // If overlap of polygons was non-trivial (created mortar segment elements)
1218  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  if (!primary_elem_found)
1224  {
1225  primary_elem_found = true;
1226  primary_elem_candidates.clear();
1227  }
1228 
1229  // Add neighbors to candidate list
1230  for (auto neighbor : primary_elem_candidate->neighbor_ptr_range())
1231  {
1232  // If not valid or not on lower dimensional secondary subdomain, skip
1233  if (neighbor == nullptr || neighbor->subdomain_id() != primary_subd_id)
1234  continue;
1235  // If already processed, skip
1236  if (processed_primary_elems.count(neighbor))
1237  continue;
1238  // Otherwise, add to candidates
1239  primary_elem_candidates.insert(neighbor);
1240  }
1241 
1245  std::vector<Node *> new_nodes;
1246  for (auto pt : nodal_points)
1247  new_nodes.push_back(_mortar_segment_mesh->add_point(
1248  pt, _mortar_segment_mesh->max_node_id(), secondary_side_elem->processor_id()));
1249 
1250  // Loop through triangular elements in map
1251  for (auto el : index_range(elem_to_node_map))
1252  {
1253  // Create new triangular element
1254  std::unique_ptr<Elem> new_elem;
1255  if (elem_to_node_map[el].size() == 3)
1256  new_elem = std::make_unique<Tri3>();
1257  else
1258  mooseError("Active mortar segments only supports TRI elements, 3 nodes expected "
1259  "but: ",
1260  elem_to_node_map[el].size(),
1261  " provided.");
1262 
1263  new_elem->processor_id() = secondary_side_elem->processor_id();
1264  new_elem->subdomain_id() = secondary_side_elem->subdomain_id();
1265  new_elem->set_id(local_id_index++);
1266 
1267  // Attach newly created nodes
1268  for (auto i : index_range(elem_to_node_map[el]))
1269  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  if (new_elem->volume() / secondary_volume < TOLERANCE)
1273  continue;
1274 
1275  // Add elements to mortar segment mesh
1276  Elem * msm_new_elem = _mortar_segment_mesh->add_elem(new_elem.release());
1277 
1278  msm_new_elem->set_extra_integer(secondary_sub_elem, sub_elem_map[el].first);
1279  msm_new_elem->set_extra_integer(primary_sub_elem, sub_elem_map[el].second);
1280 
1281  // Fill out mortar segment info
1282  MortarSegmentInfo msinfo;
1283  msinfo.secondary_elem = secondary_side_elem;
1284  msinfo.primary_elem = primary_elem_candidate;
1285 
1286  // Associate this MSM elem with the MortarSegmentInfo.
1287  _msm_elem_to_info.emplace(msm_new_elem, msinfo);
1288 
1289  // Add this mortar segment to the secondary elem to mortar segment map
1290  secondary_to_msm_element_set.insert(msm_new_elem);
1291 
1293  // Unlike for 2D, we always have a primary when building the mortar mesh so we don't
1294  // have to check for null
1296  }
1297  }
1298  // End loop through primary element candidates
1299  }
1300 
1301  for (auto sel : make_range(secondary_side_elem->n_sub_elem()))
1302  {
1303  // Check if any segments failed to project
1304  if (mortar_segment_helper[sel]->remainder() == 1.0)
1305  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  if (secondary_to_msm_element_set.empty())
1313  _secondary_elems_to_mortar_segments.erase(secondary_elem_to_msm_map_it);
1314  } // End loop through secondary elements
1315  } // End loop through mortar constraint pairs
1316 
1317  _mortar_segment_mesh->cache_elem_data();
1318 
1319  // Output mortar segment mesh
1320  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  for (const auto msm_el : _mortar_segment_mesh->active_local_element_ptr_range())
1325  if (msm_el->type() != TRI3)
1326  msm_el->subdomain_id()++;
1327 
1328  ExodusII_IO mortar_segment_mesh_writer(*_mortar_segment_mesh);
1329 
1330  // Default to non-HDF5 output for wider compatibility
1331  mortar_segment_mesh_writer.set_hdf5_writing(false);
1332 
1333  mortar_segment_mesh_writer.write("mortar_segment_mesh.e");
1334 
1335  // Undo increment
1336  for (const auto msm_el : _mortar_segment_mesh->active_local_element_ptr_range())
1337  if (msm_el->type() != TRI3)
1338  msm_el->subdomain_id()--;
1339  }
1340 
1342 
1343  // Print mortar segment mesh statistics
1344  if (_debug)
1345  {
1346  if (_mesh.n_processors() == 1)
1347  msmStatistics();
1348  else
1349  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 }
1353 
1354 void
1356 {
1357  std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type>>>
1358  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  for (const auto & pr : _msm_elem_to_info)
1367  {
1368  const Elem * secondary_elem = pr.second.secondary_elem;
1369  const Elem * primary_elem = pr.second.primary_elem;
1370 
1371  // LowerSecondary
1372  coupling_info[secondary_elem->processor_id()].emplace_back(
1373  secondary_elem->id(), secondary_elem->interior_parent()->id());
1374  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  _mortar_interface_coupling[secondary_elem->id()].insert(
1378  secondary_elem->interior_parent()->id());
1379 
1380  // LowerPrimary
1381  coupling_info[secondary_elem->processor_id()].emplace_back(
1382  secondary_elem->id(), primary_elem->interior_parent()->id());
1383  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  _mortar_interface_coupling[secondary_elem->id()].insert(
1387  primary_elem->interior_parent()->id());
1388 
1389  // Lower-LowerDimensionalPrimary
1390  coupling_info[secondary_elem->processor_id()].emplace_back(secondary_elem->id(),
1391  primary_elem->id());
1392  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  _mortar_interface_coupling[secondary_elem->id()].insert(primary_elem->id());
1396 
1397  // SecondaryLower
1398  coupling_info[secondary_elem->interior_parent()->processor_id()].emplace_back(
1399  secondary_elem->interior_parent()->id(), secondary_elem->id());
1400 
1401  // SecondaryPrimary
1402  coupling_info[secondary_elem->interior_parent()->processor_id()].emplace_back(
1403  secondary_elem->interior_parent()->id(), primary_elem->interior_parent()->id());
1404 
1405  // PrimaryLower
1406  coupling_info[primary_elem->interior_parent()->processor_id()].emplace_back(
1407  primary_elem->interior_parent()->id(), secondary_elem->id());
1408 
1409  // PrimarySecondary
1410  coupling_info[primary_elem->interior_parent()->processor_id()].emplace_back(
1411  primary_elem->interior_parent()->id(), secondary_elem->interior_parent()->id());
1412  }
1413 
1414  // Push the coupling information
1415  auto action_functor =
1416  [this](processor_id_type,
1417  const std::vector<std::pair<dof_id_type, dof_id_type>> & coupling_info)
1418  {
1419  for (auto [i, j] : coupling_info)
1420  _mortar_interface_coupling[i].insert(j);
1421  };
1422  TIMPI::push_parallel_vector_data(_mesh.comm(), coupling_info, action_functor);
1423 }
1424 
1425 void
1427 {
1428  // Print boundary pairs
1429  Moose::out << "Mortar Interface Statistics:" << std::endl;
1430 
1431  // Count number of elements on primary and secondary sides
1432  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1433  {
1434  const auto primary_subd_id = pr.first;
1435  const auto secondary_subd_id = pr.second;
1436 
1437  // Allocate statistics vectors for primary lower, secondary lower, and msm meshes
1438  StatisticsVector<Real> primary; // primary.reserve(mesh.n_elem());
1439  StatisticsVector<Real> secondary; // secondary.reserve(mesh.n_elem());
1440  StatisticsVector<Real> msm; // msm.reserve(mortar_segment_mesh->n_elem());
1441 
1442  for (auto * el : _mesh.active_element_ptr_range())
1443  {
1444  // Add secondary and primary elem volumes to statistics vector
1445  if (el->subdomain_id() == secondary_subd_id)
1446  secondary.push_back(el->volume());
1447  else if (el->subdomain_id() == primary_subd_id)
1448  primary.push_back(el->volume());
1449  }
1450 
1451  // Note: when we allow more than one primary secondary pair will need to make
1452  // separate mortar segment mesh for each
1453  for (auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
1454  {
1455  // Add msm elem volume to statistic vector
1456  msm.push_back(msm_elem->volume());
1457  }
1458 
1459  // Create table
1460  std::vector<std::string> col_names = {"mesh", "n_elems", "max", "min", "median"};
1461  std::vector<std::string> subds = {"secondary_lower", "primary_lower", "mortar_segment"};
1462  std::vector<size_t> n_elems = {secondary.size(), primary.size(), msm.size()};
1463  std::vector<Real> maxs = {secondary.maximum(), primary.maximum(), msm.maximum()};
1464  std::vector<Real> mins = {secondary.minimum(), primary.minimum(), msm.minimum()};
1465  std::vector<Real> medians = {secondary.median(), primary.median(), msm.median()};
1466 
1467  FormattedTable table;
1468  table.clear();
1469  for (auto i : index_range(subds))
1470  {
1471  table.addRow(i);
1472  table.addData<std::string>(col_names[0], subds[i]);
1473  table.addData<size_t>(col_names[1], n_elems[i]);
1474  table.addData<Real>(col_names[2], maxs[i]);
1475  table.addData<Real>(col_names[3], mins[i]);
1476  table.addData<Real>(col_names[4], medians[i]);
1477  }
1478 
1479  Moose::out << "secondary subdomain: " << secondary_subd_id
1480  << " \tprimary subdomain: " << primary_subd_id << std::endl;
1481  table.printTable(Moose::out, subds.size());
1482  }
1483 }
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
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  const Real tol = (dim() == 3) ? 0.1 : TOLERANCE;
1496 
1497  std::unordered_map<processor_id_type, std::set<dof_id_type>> proc_to_inactive_nodes_set;
1498  const auto my_pid = _mesh.processor_id();
1499 
1500  // List of inactive nodes on local secondary elements
1501  std::unordered_set<dof_id_type> inactive_node_ids;
1502 
1503  std::unordered_map<const Elem *, Real> active_volume{};
1504 
1505  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1506  for (const auto el : _mesh.active_subdomain_elements_ptr_range(pr.second))
1507  active_volume[el] = 0.;
1508 
1509  // Compute fraction of elements with corresponding primary elements
1510  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
1511  {
1512  const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
1513  const Elem * secondary_elem = msinfo.secondary_elem;
1514 
1515  active_volume[secondary_elem] += msm_elem->volume();
1516  }
1517 
1518  // Mark all inactive local nodes
1519  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1520  // Loop through all elements on my processor
1521  for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(pr.second))
1522  // If elem fully or partially dropped
1523  if (std::abs(active_volume[el] / el->volume() - 1.0) > tol)
1524  {
1525  // Add all nodes to list of inactive
1526  for (auto n : make_range(el->n_nodes()))
1527  inactive_node_ids.insert(el->node_id(n));
1528  }
1529 
1530  // Assemble list of procs that nodes contribute to
1531  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1532  {
1533  const auto secondary_subd_id = pr.second;
1534 
1535  // Loop through all elements not on my processor
1536  for (const auto el : _mesh.active_subdomain_elements_ptr_range(secondary_subd_id))
1537  {
1538  // Get processor_id
1539  const auto pid = el->processor_id();
1540 
1541  // If element is in my subdomain, skip
1542  if (pid == my_pid)
1543  continue;
1544 
1545  // If element on proc pid shares any of my inactive nodes, mark to send
1546  for (const auto n : make_range(el->n_nodes()))
1547  {
1548  const auto node_id = el->node_id(n);
1549  if (inactive_node_ids.find(node_id) != inactive_node_ids.end())
1550  proc_to_inactive_nodes_set[pid].insert(node_id);
1551  }
1552  }
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  std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_inactive_nodes_vector;
1559  for (const auto & proc_set : proc_to_inactive_nodes_set)
1560  proc_to_inactive_nodes_vector[proc_set.first].insert(
1561  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  auto action_functor = [this, &inactive_node_ids](const processor_id_type pid,
1567  const std::vector<dof_id_type> & sent_data)
1568  {
1569  if (pid == _mesh.processor_id())
1570  mooseError("Should not be communicating with self.");
1571  for (const auto pr : sent_data)
1572  inactive_node_ids.insert(pr);
1573  };
1574  TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_inactive_nodes_vector, action_functor);
1575  }
1576  _inactive_local_lm_nodes.clear();
1577  for (const auto node_id : inactive_node_ids)
1578  _inactive_local_lm_nodes.insert(_mesh.node_ptr(node_id));
1579 }
1580 
1581 void
1583 {
1585  {
1587  return;
1588  }
1589 
1590  std::unordered_map<processor_id_type, std::set<dof_id_type>> proc_to_active_nodes_set;
1591  const auto my_pid = _mesh.processor_id();
1592 
1593  // List of active nodes on local secondary elements
1594  std::unordered_set<dof_id_type> active_local_nodes;
1595 
1596  // Mark all active local nodes
1597  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
1598  {
1599  const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
1600  const Elem * secondary_elem = msinfo.secondary_elem;
1601 
1602  for (auto n : make_range(secondary_elem->n_nodes()))
1603  active_local_nodes.insert(secondary_elem->node_id(n));
1604  }
1605 
1606  // Assemble list of procs that nodes contribute to
1607  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1608  {
1609  const auto secondary_subd_id = pr.second;
1610 
1611  // Loop through all elements not on my processor
1612  for (const auto el : _mesh.active_subdomain_elements_ptr_range(secondary_subd_id))
1613  {
1614  // Get processor_id
1615  const auto pid = el->processor_id();
1616 
1617  // If element is in my subdomain, skip
1618  if (pid == my_pid)
1619  continue;
1620 
1621  // If element on proc pid shares any of my active nodes, mark to send
1622  for (const auto n : make_range(el->n_nodes()))
1623  {
1624  const auto node_id = el->node_id(n);
1625  if (active_local_nodes.find(node_id) != active_local_nodes.end())
1626  proc_to_active_nodes_set[pid].insert(node_id);
1627  }
1628  }
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  std::unordered_map<processor_id_type, std::vector<dof_id_type>> proc_to_active_nodes_vector;
1635  for (const auto & proc_set : proc_to_active_nodes_set)
1636  {
1637  proc_to_active_nodes_vector[proc_set.first].reserve(proc_to_active_nodes_set.size());
1638  for (const auto node_id : proc_set.second)
1639  proc_to_active_nodes_vector[proc_set.first].push_back(node_id);
1640  }
1641 
1642  // First push data
1643  auto action_functor = [this, &active_local_nodes](const processor_id_type pid,
1644  const std::vector<dof_id_type> & sent_data)
1645  {
1646  if (pid == _mesh.processor_id())
1647  mooseError("Should not be communicating with self.");
1648  active_local_nodes.insert(sent_data.begin(), sent_data.end());
1649  };
1650  TIMPI::push_parallel_vector_data(_mesh.comm(), proc_to_active_nodes_vector, action_functor);
1651  }
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  _inactive_local_lm_nodes.clear();
1656  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1657  for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(
1658  /*secondary_subd_id*/ pr.second))
1659  for (const auto n : make_range(el->n_nodes()))
1660  if (active_local_nodes.find(el->node_id(n)) == active_local_nodes.end())
1661  _inactive_local_lm_nodes.insert(el->node_ptr(n));
1662 }
1663 
1664 // Note: could be combined with previous routine, keeping separate for clarity (for now)
1665 void
1667 {
1668  // Mark all active secondary elements
1669  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  const Real tol = (dim() == 3) ? 0.1 : TOLERANCE;
1676 
1677  std::unordered_map<const Elem *, Real> active_volume;
1678 
1679  // Compute fraction of elements with corresponding primary elements
1681  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
1682  {
1683  const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
1684  const Elem * secondary_elem = msinfo.secondary_elem;
1685 
1686  active_volume[secondary_elem] += msm_elem->volume();
1687  }
1688  //****
1689 
1690  for (const auto msm_elem : _mortar_segment_mesh->active_local_element_ptr_range())
1691  {
1692  const MortarSegmentInfo & msinfo = _msm_elem_to_info.at(msm_elem);
1693  const Elem * secondary_elem = msinfo.secondary_elem;
1694 
1695  //****
1697  if (std::abs(active_volume[secondary_elem] / secondary_elem->volume() - 1.0) > tol)
1698  continue;
1699  //****
1700 
1701  active_local_elems.insert(secondary_elem);
1702  }
1703 
1704  // Take complement of active elements in active local subdomain to get inactive local elements
1705  _inactive_local_lm_elems.clear();
1706  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1707  for (const auto el : _mesh.active_local_subdomain_elements_ptr_range(
1708  /*secondary_subd_id*/ pr.second))
1709  if (active_local_elems.find(el) == active_local_elems.end())
1710  _inactive_local_lm_elems.insert(el);
1711 }
1712 
1713 void
1715 {
1716  // The dimension according to Mesh::mesh_dimension().
1717  const auto dim = _mesh.mesh_dimension();
1718 
1719  // A nodal lower-dimensional nodal quadrature rule to be used on faces.
1720  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  std::map<dof_id_type, std::vector<std::pair<Point, Real>>> node_to_normals_map;
1728 
1730  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  for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
1736  end_el = _mesh.active_elements_end();
1737  el != end_el;
1738  ++el)
1739  {
1740  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  if (!_secondary_boundary_subdomain_ids.count(secondary_elem->subdomain_id()))
1744  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  FEType nnx_fe_type(secondary_elem->default_order(), LAGRANGE);
1749  std::unique_ptr<FEBase> nnx_fe_face(FEBase::build(dim, nnx_fe_type));
1750  nnx_fe_face->attach_quadrature_rule(&qface);
1751  const std::vector<Point> & face_normals = nnx_fe_face->get_normals();
1752 
1753  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  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  _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  auto s = interior_parent->which_side_am_i(secondary_elem);
1770 
1771  // Reinit the face FE object on side s.
1772  nnx_fe_face->reinit(interior_parent, s);
1773 
1774  for (MooseIndex(secondary_elem->n_nodes()) n = 0; n < secondary_elem->n_nodes(); ++n)
1775  {
1776  auto & normals_and_weights_vec = node_to_normals_map[secondary_elem->node_id(n)];
1777  normals_and_weights_vec.push_back(std::make_pair(sign * face_normals[n], JxW[n]));
1778  }
1779  }
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  for (const auto & pr : node_to_normals_map)
1785  {
1786  // Compute normal vector
1787  const auto & node_id = pr.first;
1788  const auto & normals_and_weights_vec = pr.second;
1789 
1790  Point nodal_normal;
1791  for (const auto & norm_and_weight : normals_and_weights_vec)
1792  nodal_normal += norm_and_weight.first * norm_and_weight.second;
1793  nodal_normal = nodal_normal.unit();
1794 
1795  _secondary_node_to_nodal_normal[_mesh.node_ptr(node_id)] = nodal_normal;
1796 
1797  Point nodal_tangent_one;
1798  Point nodal_tangent_two;
1799  householderOrthogolization(nodal_normal, nodal_tangent_one, nodal_tangent_two);
1800 
1801  _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)][0] = nodal_tangent_one;
1802  _secondary_node_to_hh_nodal_tangents[_mesh.node_ptr(node_id)][1] = nodal_tangent_two;
1803  }
1804 }
1805 
1806 void
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  const Real nx = nodal_normal(0);
1815  const Real ny = nodal_normal(1);
1816  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  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  if (std::abs(h_vector(0)) < TOLERANCE)
1827  {
1828  nodal_tangent_one(0) = 0;
1829  nodal_tangent_one(1) = 1;
1830  nodal_tangent_one(2) = 0;
1831 
1832  nodal_tangent_two(0) = 0;
1833  nodal_tangent_two(1) = 0;
1834  nodal_tangent_two(2) = -1;
1835 
1836  return;
1837  }
1838 
1839  const Real h = h_vector.norm();
1840 
1841  nodal_tangent_one(0) = -2.0 * h_vector(0) * h_vector(1) / (h * h);
1842  nodal_tangent_one(1) = 1.0 - 2.0 * h_vector(1) * h_vector(1) / (h * h);
1843  nodal_tangent_one(2) = -2.0 * h_vector(1) * h_vector(2) / (h * h);
1844 
1845  nodal_tangent_two(0) = -2.0 * h_vector(0) * h_vector(2) / (h * h);
1846  nodal_tangent_two(1) = -2.0 * h_vector(1) * h_vector(2) / (h * h);
1847  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
1854 {
1855  // For each primary/secondary boundary id pair, call the
1856  // project_secondary_nodes_single_pair() helper function.
1857  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
1858  projectSecondaryNodesSinglePair(pr.first, pr.second);
1859 }
1860 
1861 void
1863  SubdomainID lower_dimensional_primary_subdomain_id,
1864  SubdomainID lower_dimensional_secondary_subdomain_id)
1865 {
1866  // Build the "subdomain" adaptor based KD Tree.
1867  NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, lower_dimensional_primary_subdomain_id);
1868  subdomain_kd_tree_t kd_tree(
1869  3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
1870 
1871  // Construct the KD tree.
1872  kd_tree.buildIndex();
1873 
1874  for (MeshBase::const_element_iterator el = _mesh.active_elements_begin(),
1875  end_el = _mesh.active_elements_end();
1876  el != end_el;
1877  ++el)
1878  {
1879  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  if (secondary_side_elem->subdomain_id() != lower_dimensional_secondary_subdomain_id)
1883  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  for (MooseIndex(secondary_side_elem->n_vertices()) n = 0; n < secondary_side_elem->n_vertices();
1890  ++n)
1891  {
1892  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  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  bool is_mapped = true;
1902  for (MooseIndex(secondary_node_neighbors) snn = 0; snn < secondary_node_neighbors.size();
1903  ++snn)
1904  {
1905  auto secondary_key = std::make_pair(secondary_node, secondary_node_neighbors[snn]);
1906  if (!_secondary_node_and_elem_to_xi2_primary_elem.count(secondary_key))
1907  {
1908  is_mapped = false;
1909  break;
1910  }
1911  }
1912 
1913  // Go to the next node if this one has already been mapped.
1914  if (is_mapped)
1915  continue;
1916 
1917  // Look up the new nodal normal value in the local storage, error if not found.
1918  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  {(*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  const std::size_t num_results = 3;
1929 
1930  // Initialize result_set and do the search.
1931  std::vector<size_t> ret_index(num_results);
1932  std::vector<Real> out_dist_sqr(num_results);
1933  nanoflann::KNNResultSet<Real> result_set(num_results);
1934  result_set.init(&ret_index[0], &out_dist_sqr[0]);
1935  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  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  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  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  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  for (MooseIndex(primary_elem_candidates) e = 0; e < primary_elem_candidates.size(); ++e)
1962  {
1963  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  if (rejected_primary_elem_candidates.count(primary_elem_candidate))
1967  continue;
1968 
1969  // Now generically solve for xi2
1970  auto && order = primary_elem_candidate->default_order();
1971  DualNumber<Real> xi2_dn{0, 1};
1972  unsigned int current_iterate = 0, max_iterates = 10;
1973 
1974  // Newton loop
1975  do
1976  {
1978  for (MooseIndex(primary_elem_candidate->n_nodes()) n = 0;
1979  n < primary_elem_candidate->n_nodes();
1980  ++n)
1981  x2 +=
1982  Moose::fe_lagrange_1D_shape(order, n, xi2_dn) * primary_elem_candidate->point(n);
1983  const auto u = x2 - (*secondary_node);
1984  const auto F = u(0) * nodal_normal(1) - u(1) * nodal_normal(0);
1985 
1986  if (std::abs(F) < _newton_tolerance)
1987  break;
1988 
1989  if (F.derivatives())
1990  {
1991  Real dxi2 = -F.value() / F.derivatives();
1992 
1993  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  current_iterate = max_iterates;
2000  } while (++current_iterate < max_iterates);
2001 
2002  Real xi2 = xi2_dn.value();
2003 
2004  // Check whether the projection worked. The last condition checks for obliqueness of the
2005  // projection
2006  if ((current_iterate < max_iterates) && (std::abs(xi2) <= 1. + _xi_tolerance) &&
2007  (std::abs(
2008  (primary_elem_candidate->point(0) - primary_elem_candidate->point(1)).unit() *
2009  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  if (std::abs(std::abs(xi2) - 1.) < _xi_tolerance * 5.0)
2026  {
2027  const Node * primary_node = (xi2 < 0) ? primary_elem_candidate->node_ptr(0)
2028  : primary_elem_candidate->node_ptr(1);
2029 
2030  const std::vector<const Elem *> & primary_node_neighbors =
2031  _nodes_to_primary_elem_map.at(primary_node->id());
2032 
2033  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  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  for (MooseIndex(secondary_node_neighbors) nn = 0;
2046  nn < secondary_node_neighbors.size();
2047  ++nn)
2048  {
2049  const Elem * secondary_neigh = secondary_node_neighbors[nn];
2050  Point opposite = (secondary_neigh->node_ptr(0) == secondary_node)
2051  ? secondary_neigh->point(1)
2052  : secondary_neigh->point(0);
2053  Point cp = nodal_normal.cross(opposite - *secondary_node);
2054  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  for (MooseIndex(primary_node_neighbors) nn = 0; nn < primary_node_neighbors.size();
2060  ++nn)
2061  {
2062  const Elem * primary_neigh = primary_node_neighbors[nn];
2063  Point opposite = (primary_neigh->node_ptr(0) == primary_node)
2064  ? primary_neigh->point(1)
2065  : primary_neigh->point(0);
2066  Point cp = nodal_normal.cross(opposite - *secondary_node);
2067  primary_node_neighbor_cps[nn] = cp(2);
2068  }
2069 
2070  // Associate secondary/primary elems on matching sides.
2071  bool found_match = false;
2072  for (MooseIndex(secondary_node_neighbors) snn = 0;
2073  snn < secondary_node_neighbors.size();
2074  ++snn)
2075  for (MooseIndex(primary_node_neighbors) mnn = 0;
2076  mnn < primary_node_neighbors.size();
2077  ++mnn)
2078  if (secondary_node_neighbor_cps[snn] * primary_node_neighbor_cps[mnn] > 0)
2079  {
2080  found_match = true;
2081  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  Real xi2 = (primary_node == primary_node_neighbors[mnn]->node_ptr(0)) ? -1 : +1;
2086  auto secondary_key =
2087  std::make_pair(secondary_node, secondary_node_neighbors[snn]);
2088  auto primary_val = std::make_pair(xi2, primary_node_neighbors[mnn]);
2089  _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  (secondary_node == secondary_node_neighbors[snn]->node_ptr(0)) ? -1 : +1;
2095 
2096  auto primary_key = std::make_tuple(
2097  primary_node->id(), primary_node, primary_node_neighbors[mnn]);
2098  auto secondary_val = std::make_pair(xi1, secondary_node_neighbors[snn]);
2100  secondary_val);
2101  }
2102 
2103  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  rejected_primary_elem_candidates.insert(primary_elem_candidate);
2108  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  if (secondary_node_neighbors.size() == 1 && primary_node_neighbors.size() == 2)
2115  for (auto it = primary_elems_mapped.begin(); it != primary_elems_mapped.end(); ++it)
2116  if (*it == false)
2117  {
2118  auto index = std::distance(primary_elems_mapped.begin(), it);
2120  std::make_tuple(
2121  primary_node->id(), primary_node, primary_node_neighbors[index]),
2122  std::make_pair(1, nullptr));
2123  }
2124  }
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  for (MooseIndex(secondary_node_neighbors) nn = 0;
2129  nn < secondary_node_neighbors.size();
2130  ++nn)
2131  {
2132  const Elem * neigh = secondary_node_neighbors[nn];
2133  for (MooseIndex(neigh->n_vertices()) nid = 0; nid < neigh->n_vertices(); ++nid)
2134  {
2135  const Node * neigh_node = neigh->node_ptr(nid);
2136  if (secondary_node == neigh_node)
2137  {
2138  auto key = std::make_pair(neigh_node, neigh);
2139  auto val = std::make_pair(xi2, primary_elem_candidate);
2141  }
2142  }
2143  }
2144  }
2145 
2146  projection_succeeded = true;
2147  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  rejected_primary_elem_candidates.insert(primary_elem_candidate);
2152  }
2153 
2154  if (projection_succeeded)
2155  break; // out of r-loop
2156  } // r-loop
2157 
2158  if (!projection_succeeded && _debug)
2159  _console << "Failed to find primary Elem into which secondary node "
2160  << static_cast<const Point &>(*secondary_node) << " was projected." << std::endl
2161  << std::endl;
2162  } // loop over side nodes
2163  } // end loop over lower-dimensional elements
2164 }
2165 
2166 // Inverse map primary nodes onto their corresponding secondary elements for each primary/secondary
2167 // pair.
2168 void
2170 {
2171  // For each primary/secondary boundary id pair, call the
2172  // project_primary_nodes_single_pair() helper function.
2173  for (const auto & pr : _primary_secondary_subdomain_id_pairs)
2174  projectPrimaryNodesSinglePair(pr.first, pr.second);
2175 }
2176 
2177 void
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  NanoflannMeshSubdomainAdaptor<3> mesh_adaptor(_mesh, lower_dimensional_secondary_subdomain_id);
2184  subdomain_kd_tree_t kd_tree(
2185  3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
2186 
2187  // Construct the KD tree for lower-dimensional elements in the volume mesh.
2188  kd_tree.buildIndex();
2189 
2190  std::unordered_set<dof_id_type> primary_nodes_visited;
2191 
2192  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  if (primary_side_elem->subdomain_id() != lower_dimensional_primary_subdomain_id)
2196  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  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  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  _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  std::make_tuple(primary_node->id(), primary_node, primary_node_neighbors[0]);
2216  if (!primary_nodes_visited.insert(primary_node->id()).second ||
2218  continue;
2219 
2220  // Data structure for performing Nanoflann searches.
2221  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  const size_t num_results = 3;
2228 
2229  // Initialize result_set and do the search.
2230  std::vector<size_t> ret_index(num_results);
2231  std::vector<Real> out_dist_sqr(num_results);
2232  nanoflann::KNNResultSet<Real> result_set(num_results);
2233  result_set.init(&ret_index[0], &out_dist_sqr[0]);
2234  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  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  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  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  _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  for (MooseIndex(secondary_elem_candidates) e = 0; e < secondary_elem_candidates.size(); ++e)
2260  {
2261  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  if (rejected_secondary_elem_candidates.count(secondary_elem_candidate))
2265  continue;
2266 
2267  std::vector<Point> nodal_normals(secondary_elem_candidate->n_nodes());
2268  for (const auto n : make_range(secondary_elem_candidate->n_nodes()))
2269  nodal_normals[n] =
2270  _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  DualNumber<Real> xi1_dn{0, 1}; // initial guess
2276  auto && order = secondary_elem_candidate->default_order();
2277  unsigned int current_iterate = 0, max_iterates = 10;
2278 
2279  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  {
2288  for (MooseIndex(secondary_elem_candidate->n_nodes()) n = 0;
2289  n < secondary_elem_candidate->n_nodes();
2290  ++n)
2291  {
2292  const auto phi = Moose::fe_lagrange_1D_shape(order, n, xi1_dn);
2293  x1 += phi * secondary_elem_candidate->point(n);
2294  normals += phi * nodal_normals[n];
2295  }
2296 
2297  const auto u = x1 - (*primary_node);
2298 
2299  const auto F = u(0) * normals(1) - u(1) * normals(0);
2300 
2301  if (std::abs(F) < _newton_tolerance)
2302  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  Real dxi1 = -F.value() / F.derivatives();
2308 
2309  xi1_dn += dxi1;
2310 
2311  normals = 0;
2312  } while (++current_iterate < max_iterates);
2313 
2314  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  if ((current_iterate < max_iterates) && (std::abs(xi1) <= 1. + _xi_tolerance) &&
2319  (std::abs((primary_side_elem->point(0) - primary_side_elem->point(1)).unit() *
2320  MetaPhysicL::raw_value(normals).unit()) <
2322  {
2323  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  throw MooseException("Nodes on primary and secondary surfaces are aligned. This is "
2332  "causing trouble when identifying projections from secondary "
2333  "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  const Elem * neigh = primary_node_neighbors[0];
2345  for (MooseIndex(neigh->n_vertices()) nid = 0; nid < neigh->n_vertices(); ++nid)
2346  {
2347  const Node * neigh_node = neigh->node_ptr(nid);
2348  if (primary_node == neigh_node)
2349  {
2350  auto key = std::make_tuple(neigh_node->id(), neigh_node, neigh);
2351  auto val = std::make_pair(xi1, secondary_elem_candidate);
2353  }
2354  }
2355  }
2356 
2357  projection_succeeded = true;
2358  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  rejected_secondary_elem_candidates.insert(secondary_elem_candidate);
2364  }
2365  } // end e-loop over candidate elems
2366 
2367  if (projection_succeeded)
2368  break; // out of r-loop
2369  } // r-loop
2370 
2371  if (!projection_succeeded && _debug)
2372  {
2373  _console << "Failed to find point from which primary node "
2374  << static_cast<const Point &>(*primary_node) << " was projected." << std::endl
2375  << std::endl;
2376  }
2377  } // loop over side nodes
2378  } // end loop over elements for finding where primary points would have projected from.
2379 }
2380 
2381 std::vector<AutomaticMortarGeneration::MortarFilterIter>
2383 {
2384  auto secondary_it = _nodes_to_secondary_elem_map.find(node.id());
2385  if (secondary_it == _nodes_to_secondary_elem_map.end())
2386  return {};
2387 
2388  const auto & secondary_elems = secondary_it->second;
2389  std::vector<MortarFilterIter> ret;
2390  ret.reserve(secondary_elems.size());
2391 
2392  for (const auto i : index_range(secondary_elems))
2393  {
2394  auto * const secondary_elem = secondary_elems[i];
2395  auto msm_it = _secondary_elems_to_mortar_segments.find(secondary_elem->id());
2396  if (msm_it == _secondary_elems_to_mortar_segments.end())
2397  // We may have removed this element key from this map
2398  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  ret.push_back(msm_it);
2407  }
2408 
2409  return ret;
2410 }
virtual T maximum() const
std::set< SubdomainID > _primary_boundary_subdomain_ids
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:42
const Elem * getSecondaryLowerdElemFromSecondaryElem(dof_id_type secondary_elem_id) const
Return lower dimensional secondary element given its interior parent.
T fe_lagrange_2D_shape(const libMesh::ElemType type, const Order order, const unsigned int i, const VectorType< T > &p)
ElemType
std::unique_ptr< FEGenericBase< Real > > build(const unsigned int dim, const FEType &fet)
unique_id_type & set_unique_id()
std::unordered_set< const Elem * > _inactive_local_lm_elems
List of inactive lagrange multiplier nodes (for elemental variables)
auto norm() const -> decltype(std::norm(Real()))
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
Definition: MooseUtils.h:380
QUAD8
unsigned int get_node_index(const Node *node_ptr) const
void computeInactiveLMNodes()
Get list of secondary nodes that don&#39;t contribute to interaction with any primary element...
std::vector< std::pair< BoundaryID, BoundaryID > > _primary_secondary_boundary_id_pairs
A list of primary/secondary boundary id pairs corresponding to each side of the mortar interface...
const Elem * interior_parent() const
MPI_Info info
void clear()
Clears the mortar segment mesh and accompanying data structures.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
std::set< SubdomainID > _secondary_boundary_subdomain_ids
The secondary/primary lower-dimensional boundary subdomain ids are the secondary/primary boundary ids...
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
std::unordered_map< dof_id_type, std::set< Elem *, CompareDofObjectsByID > > _secondary_elems_to_mortar_segments
We maintain a mapping from lower-dimensional secondary elements in the original mesh to (sets of) ele...
std::map< unsigned int, unsigned int > getPrimaryIpToLowerElementMap(const Elem &primary_elem, const Elem &primary_elem_ip, const Elem &lower_secondary_elem) const
Compute on-the-fly mapping from primary interior parent nodes to its corresponding lower dimensional ...
const bool _periodic
Whether this object will be generating a mortar segment mesh for periodic constraints.
void swap(std::vector< T > &data, const std::size_t idx0, const std::size_t idx1, const libMesh::Parallel::Communicator &comm)
Swap function for serial or distributed vector of data.
Definition: Shuffle.h:494
void addData(const std::string &name, const T &value)
Method for adding data to the output table.
unsigned int which_side_am_i(const Elem *e) const
Real _newton_tolerance
Newton solve tolerance for node projections.
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:336
void output() override
Overload this function with the desired output activities.
MortarNodalGeometryOutput(const InputParameters &params)
void push_back(const T &v)
Definition: MooseUtils.h:1120
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
std::unordered_map< const Elem *, MortarSegmentInfo > _msm_elem_to_info
Map between Elems in the mortar segment mesh and their info structs.
const bool _distributed
Whether the mortar segment mesh is distributed.
Base class for MOOSE-based applications.
Definition: MooseApp.h:96
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
void buildMortarSegmentMesh()
Builds the mortar segment mesh once the secondary and primary node projections have been completed...
const Parallel::Communicator & comm() const
void msmStatistics()
Outputs mesh statistics for mortar segment mesh.
void buildCouplingInformation()
build the _mortar_interface_coupling data
std::vector< Point > getNodalNormals(const Elem &secondary_elem) const
Special adaptor that works with subdomains of the Mesh.
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
static const Real invalid_xi
SECOND
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
FEProblemBase & feProblem() const
Definition: MooseApp.C:1961
void projectSecondaryNodes()
Project secondary nodes (find xi^(2) values) to the closest points on the primary surface...
TRI3
const std::unordered_map< dof_id_type, std::set< Elem *, CompareDofObjectsByID > > & secondariesToMortarSegments() const
Utility class template for a semidynamic vector with a maximum size N and a chosen dynamic size...
Definition: MooseUtils.h:1096
QUAD4
This class is a container/interface for the objects involved in automatic generation of mortar spaces...
const bool _debug
Whether to print debug output.
Based class for output objects.
Definition: Output.h:43
const Elem * primary_elem
void push_parallel_vector_data(const Communicator &comm, MapToVectors &&data, const ActionFunctor &act_on_data)
uint8_t processor_id_type
std::set< SubdomainID > _secondary_ip_sub_ids
All the secondary interior parent subdomain IDs associated with the mortar mesh.
virtual libMesh::EquationSystems & es()=0
processor_id_type n_processors() const
TypeVector< Real > unit() const
void libmesh_ignore(const Args &...)
void projectPrimaryNodes()
(Inverse) project primary nodes to the points on the secondary surface where they would have come fro...
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
const Node & node_ref(const unsigned int i) const
dof_id_type id() const
AutomaticMortarGeneration(MooseApp &app, MeshBase &mesh_in, const std::pair< BoundaryID, BoundaryID > &boundary_key, const std::pair< SubdomainID, SubdomainID > &subdomain_key, bool on_displaced, bool periodic, const bool debug, const bool correct_edge_dropping, const Real minimum_projection_angle)
Must be constructed with a reference to the Mesh we are generating mortar spaces for.
virtual unsigned int n_nodes() const=0
void buildNodeToElemMaps()
Once the secondary_requested_boundary_ids and primary_requested_boundary_ids containers have been fil...
std::unordered_map< dof_id_type, const Elem * > _secondary_element_to_secondary_lowerd_element
Map from full dimensional secondary element id to lower dimensional secondary element.
virtual void write_equation_systems(const std::string &fname, const EquationSystems &es, const std::set< std::string > *system_names=nullptr) override
TRI6
std::map< std::tuple< dof_id_type, const Node *, const Elem * >, std::pair< Real, const Elem * > > _primary_node_and_elem_to_xi1_secondary_elem
Same type of container, but for mapping (Primary Node ID, Primary Node, Primary Elem) -> (xi^(1)...
void projectPrimaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id, SubdomainID lower_dimensional_secondary_subdomain_id)
Helper function used internally by AutomaticMortarGeneration::project_primary_nodes().
std::unordered_map< dof_id_type, std::vector< const Elem * > > _nodes_to_primary_elem_map
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
const bool _correct_edge_dropping
Flag to enable regressed treatment of edge dropping where all LM DoFs on edge dropping element are st...
An inteface for the _console for outputting to the Console object.
T sign(T x)
Definition: MathUtils.h:84
std::array< MooseUtils::SemidynamicVector< Point, 9 >, 2 > getNodalTangents(const Elem &secondary_elem) const
Compute the two nodal tangents, which are built on-the-fly.
std::unique_ptr< InputParameters > _output_params
Storage for the input parameters used by the mortar nodal geometry output.
std::set< BoundaryID > _primary_requested_boundary_ids
The boundary ids corresponding to all the primary surfaces.
std::map< unsigned int, unsigned int > getSecondaryIpToLowerElementMap(const Elem &lower_secondary_elem) const
Compute on-the-fly mapping from secondary interior parent nodes to lower dimensional nodes...
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
nanoflann::KDTreeSingleIndexAdaptor< subdomain_adatper_t, NanoflannMeshSubdomainAdaptor< 3 >, 3 > subdomain_kd_tree_t
This class is used for building, formatting, and outputting tables of numbers.
void computeNodalGeometry()
Computes and stores the nodal normal/tangent vectors in a local data structure instead of using the E...
Real _xi_tolerance
Tolerance for checking projection xi values.
static InputParameters validParams()
void computeIncorrectEdgeDroppingInactiveLMNodes()
Computes inactive secondary nodes when incorrect edge dropping behavior is enabled (any node touching...
void buildMortarSegmentMesh3d()
Builds the mortar segment mesh once the secondary and primary node projections have been completed...
const Elem * secondary_elem
void projectSecondaryNodesSinglePair(SubdomainID lower_dimensional_primary_subdomain_id, SubdomainID lower_dimensional_secondary_subdomain_id)
Helper function responsible for projecting secondary nodes onto primary elements for a single primary...
Provides a way for users to bail out of the current solve.
void printTable(std::ostream &out, unsigned int last_n_entries=0)
Methods for dumping the table to the stream - either by filename or by stream handle.
virtual unsigned int n_vertices() const=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::unordered_map< std::pair< const Node *, const Elem * >, std::pair< Real, const Elem * > > _secondary_node_and_elem_to_xi2_primary_elem
Similar to the map above, but associates a (Secondary Node, Secondary Elem) pair to a (xi^(2)...
Holds xi^(1), xi^(2), and other data for a given mortar segment.
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:78
virtual SimpleRange< element_iterator > active_local_subdomain_elements_ptr_range(subdomain_id_type sid)=0
subdomain_id_type subdomain_id() const
const bool _on_displaced
Whether this object is on the displaced mesh.
const Node * node_ptr(const unsigned int i) const
auto norm_sq(const T &a) -> decltype(std::norm(a))
virtual void write(const std::string &fname) override
void initOutput()
initialize mortar-mesh based output
void addOutput(std::shared_ptr< Output > output)
Adds an existing output object to the warehouse.
virtual Real volume() const
std::unordered_map< const Node *, std::array< Point, 2 > > _secondary_node_to_hh_nodal_tangents
Container for storing the nodal tangent/binormal vectors associated with each secondary node (Househo...
IntRange< T > make_range(T beg, T end)
TRI7
unsigned int mesh_dimension() const
unsigned int level ElemType type std::set< subdomain_id_type > ss processor_id_type pid unsigned int level std::set< subdomain_id_type > virtual ss SimpleRange< element_iterator > active_subdomain_elements_ptr_range(subdomain_id_type sid)=0
virtual T minimum() const
T fe_lagrange_1D_shape(const Order order, const unsigned int i, const T &xi)
const Real _minimum_projection_angle
Parameter to control which angle (in degrees) is admissible for the creation of mortar segments...
QUAD9
MeshBase & _mesh
Reference to the mesh stored in equation_systems.
virtual const Point & point(const dof_id_type i) const=0
std::vector< std::pair< SubdomainID, SubdomainID > > _primary_secondary_subdomain_id_pairs
A list of primary/secondary subdomain id pairs corresponding to each side of the mortar interface...
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
SimpleRange< NeighborPtrIter > neighbor_ptr_range()
void addRow(Real time)
Force a new row in the table with the passed in time.
virtual unsigned int n_sub_elem() const=0
std::set< BoundaryID > _secondary_requested_boundary_ids
The boundary ids corresponding to all the secondary surfaces.
std::unordered_map< const Node *, Point > _secondary_node_to_nodal_normal
Container for storing the nodal normal vector associated with each secondary node.
virtual const Node * node_ptr(const dof_id_type i) const=0
std::unordered_set< const Node * > _inactive_local_lm_nodes
processor_id_type processor_id() const
virtual Order default_order() const=0
std::unordered_map< dof_id_type, std::unordered_set< dof_id_type > > _mortar_interface_coupling
Used by the AugmentSparsityOnInterface functor to determine whether a given Elem is coupled to any ot...
std::set< SubdomainID > _primary_ip_sub_ids
All the primary interior parent subdomain IDs associated with the mortar mesh.
SearchParams SearchParameters
AutomaticMortarGeneration & _amg
The mortar generation object that we will query for nodal normal and tangent information.
void set_hdf5_writing(bool write_hdf5)
processor_id_type processor_id() const
void householderOrthogolization(const Point &normal, Point &tangent_one, Point &tangent_two) const
Householder orthogonalization procedure to obtain proper basis for tangent and binormal vectors...
virtual ElemType type() const=0
dof_id_type node_id(const unsigned int i) const
std::vector< Point > getNormals(const Elem &secondary_elem, const std::vector< Point > &xi1_pts) const
Compute the normals at given reference points on a secondary element.
static InputParameters validParams()
Definition: Output.C:32
MooseApp & _app
The Moose app.
const Point & point(const unsigned int i) const
std::unique_ptr< MeshBase > _mortar_segment_mesh
1D Mesh of mortar segment elements which gets built by the call to build_mortar_segment_mesh().
auto index_range(const T &sizable)
void set_extra_integer(const unsigned int index, const dof_id_type value)
OutputWarehouse & getOutputWarehouse()
Get the OutputWarehouse objects.
Definition: MooseApp.C:2414
uint8_t dof_id_type
std::unordered_map< dof_id_type, std::vector< const Elem * > > _nodes_to_secondary_elem_map
Map from nodes to connected lower-dimensional elements on the secondary/primary subdomains.
std::unordered_map< const Elem *, unsigned int > _lower_elem_to_side_id
Keeps track of the mapping between lower-dimensional elements and the side_id of the interior_parent ...
const Real pi
void computeInactiveLMElems()
Get list of secondary elems without any corresponding primary elements.