https://mooseframework.inl.gov
Functions
Moose::Mortar Namespace Reference

Functions

std::vector< unsigned intgetMortarSubElementNodeIndices (const Elem &parent_elem, unsigned int sub_elem)
 Return the node indices for a first-order sub-element of a parent face. More...
 
void projectQPoints3d (const Elem *msm_elem, const Elem *primal_elem, unsigned int sub_elem_index, const QBase &qrule_msm, std::vector< Point > &q_pts)
 3D projection operator for mapping qpoints on mortar segments to secondary or primary elements More...
 
void mapQPoints3dFromReference (const Elem &mortar_segment_elem, const MortarSegmentReferencePoints &reference_points, const QBase &qrule_msm, std::vector< Point > &secondary_q_pts, std::vector< Point > &primary_q_pts)
 3D mapping operator that interpolates stored parent reference points on each triangular mortar segment. More...
 
template<typename Iterators , typename Consumers , typename ActionFunctor >
void loopOverMortarSegments (const Iterators &secondary_elems_to_mortar_segments, Assembly &assembly, SubProblem &subproblem, FEProblemBase &fe_problem, const AutomaticMortarGeneration &amg, const bool displaced, const Consumers &consumers, const THREAD_ID tid, const std::map< SubdomainID, std::deque< MaterialBase *>> &secondary_ip_sub_to_mats, const std::map< SubdomainID, std::deque< MaterialBase *>> &primary_ip_sub_to_mats, const std::deque< MaterialBase *> &secondary_boundary_mats, const ActionFunctor act, const bool reinit_mortar_user_objects)
 This method will loop over pairs of secondary elements and their corresponding mortar segments, reinitialize all finite element shape functions, variables, and material properties, and then call a provided action function for each mortar segment. More...
 
template<typename Consumers >
void setupMortarMaterials (const Consumers &consumers, FEProblemBase &fe_problem, const AutomaticMortarGeneration &amg, const THREAD_ID tid, std::map< SubdomainID, std::deque< MaterialBase *>> &secondary_ip_sub_to_mats, std::map< SubdomainID, std::deque< MaterialBase *>> &primary_ip_sub_to_mats, std::deque< MaterialBase *> &secondary_boundary_mats)
 This function creates containers of materials necessary to execute the mortar method for a supplied set of consumers. More...
 

Function Documentation

◆ getMortarSubElementNodeIndices()

std::vector< unsigned int > Moose::Mortar::getMortarSubElementNodeIndices ( const Elem &  parent_elem,
unsigned int  sub_elem 
)

Return the node indices for a first-order sub-element of a parent face.

Definition at line 27 of file MortarUtils.C.

Referenced by AutomaticMortarGeneration::buildMortarSegmentMesh3d(), and projectQPoints3d().

28 {
29  if (sub_elem >= parent_elem.n_sub_elem())
30  mooseError("Invalid 3D mortar sub-element index ",
31  sub_elem,
32  " for parent element ",
33  parent_elem.id(),
34  " of type ",
35  libMesh::Utility::enum_to_string<ElemType>(parent_elem.type()),
36  ", which has ",
37  parent_elem.n_sub_elem(),
38  " sub-elements.");
39 
40  switch (parent_elem.type())
41  {
42  case TRI3:
43  return {0, 1, 2};
44  case QUAD4:
45  return {0, 1, 2, 3};
46  case TRI6:
47  case TRI7:
48  switch (sub_elem)
49  {
50  case 0:
51  return {0, 3, 5};
52  case 1:
53  return {3, 4, 5};
54  case 2:
55  return {3, 1, 4};
56  case 3:
57  return {5, 4, 2};
58  default:
59  mooseError("Invalid 3D mortar triangular sub-element index ", sub_elem, ".");
60  }
61  case QUAD8:
62  switch (sub_elem)
63  {
64  case 0:
65  return {0, 4, 7};
66  case 1:
67  return {4, 1, 5};
68  case 2:
69  return {5, 2, 6};
70  case 3:
71  return {7, 6, 3};
72  case 4:
73  return {4, 5, 6, 7};
74  default:
75  mooseError("Invalid 3D mortar QUAD8 sub-element index ", sub_elem, ".");
76  }
77  case QUAD9:
78  switch (sub_elem)
79  {
80  case 0:
81  return {0, 4, 8, 7};
82  case 1:
83  return {4, 1, 5, 8};
84  case 2:
85  return {8, 5, 2, 6};
86  case 3:
87  return {7, 8, 6, 3};
88  default:
89  mooseError("Invalid 3D mortar QUAD9 sub-element index ", sub_elem, ".");
90  }
91  default:
92  mooseError("Parent face element ",
93  parent_elem.id(),
94  " has unsupported type ",
95  libMesh::Utility::enum_to_string<ElemType>(parent_elem.type()),
96  " for 3D mortar sub-element topology.");
97  }
98 }
QUAD8
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:311
TRI3
QUAD4
TRI6
TRI7
QUAD9

◆ loopOverMortarSegments()

template<typename Iterators , typename Consumers , typename ActionFunctor >
void Moose::Mortar::loopOverMortarSegments ( const Iterators &  secondary_elems_to_mortar_segments,
Assembly assembly,
SubProblem subproblem,
FEProblemBase fe_problem,
const AutomaticMortarGeneration amg,
const bool  displaced,
const Consumers &  consumers,
const THREAD_ID  tid,
const std::map< SubdomainID, std::deque< MaterialBase *>> &  secondary_ip_sub_to_mats,
const std::map< SubdomainID, std::deque< MaterialBase *>> &  primary_ip_sub_to_mats,
const std::deque< MaterialBase *> &  secondary_boundary_mats,
const ActionFunctor  act,
const bool  reinit_mortar_user_objects 
)

This method will loop over pairs of secondary elements and their corresponding mortar segments, reinitialize all finite element shape functions, variables, and material properties, and then call a provided action function for each mortar segment.

Parameters
secondary_elems_to_mortar_segmentsThis is a container of iterators. Each iterator should point to a pair. The first member of the pair should be a pointer to a secondary face element and the second member of the pair should correspond to a container of mortar segment element pointers that correspond to the secondary face element.
assemblyThe object we will to use to reinitalize finite element data
subproblemThe object we will use to reinitialize variables
fe_problemThe object we will use to reinitialize material properties
amgThe mortar mesh generation object which holds all the mortar mesh data
displacedWhether the mortar mesh was built from a displaced parent mesh
consumersA container of objects that are going to be using all the data that we are reinitializing within this function. This may be, for instance, a container of mortar constraints or auxiliary kernels. This consumers parameter is important as it allows us to build up variable and material property dependencies that we must make sure we reinit
actThe action functor that we will call for each mortar segment after we have reinitalized all of our prereq data. This functor may, for instance, call computeResidual or computeJacobian on mortar constraints, or computeValue for an auxiliary kernel

Definition at line 91 of file MortarUtils.h.

Referenced by MortarNodalAuxKernelTempl< ComputeValueType >::compute(), MortarUserObjectThread::operator()(), and ComputeMortarFunctor::operator()().

105 {
106  const auto & primary_secondary_boundary_id_pair = amg.primarySecondaryBoundaryIDPair();
107 
108  const auto primary_boundary_id = primary_secondary_boundary_id_pair.first;
109  const auto secondary_boundary_id = primary_secondary_boundary_id_pair.second;
110 
111  // For 3D mortar get index for retrieving sub-element info
112  unsigned int secondary_sub_elem_index = 0, primary_sub_elem_index = 0;
113  if (amg.dim() == 3)
114  {
115  secondary_sub_elem_index = amg.mortarSegmentMesh().get_elem_integer_index("secondary_sub_elem");
116  primary_sub_elem_index = amg.mortarSegmentMesh().get_elem_integer_index("primary_sub_elem");
117  }
118 
119  // The mortar quadrature rule. Necessary for sizing the number of custom points for re-init'ing
120  // the secondary interior, primary interior, and secondary face elements
121  const auto & qrule_msm = assembly.qRuleMortar();
122 
123  // The element Jacobian times weights
124  const auto & JxW_msm = assembly.jxWMortar();
125 
126  // Set required material properties
127  std::unordered_set<unsigned int> needed_mat_props;
128  for (const auto & consumer : consumers)
129  {
130  const auto & mp_deps = consumer->getMatPropDependencies();
131  needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
132  }
133  fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/tid);
134 
135  // Loop through secondary elements, accumulating quadrature points for all corresponding mortar
136  // segments
137  for (const auto elem_to_msm : secondary_elems_to_mortar_segments)
138  {
139  const Elem * secondary_face_elem = subproblem.mesh().getMesh().elem_ptr(elem_to_msm->first);
140  // Set the secondary interior parent and side ids
141  const Elem * secondary_ip = secondary_face_elem->interior_parent();
142  unsigned int secondary_side_id = secondary_ip->which_side_am_i(secondary_face_elem);
143  const auto & secondary_ip_mats =
144  libmesh_map_find(secondary_ip_sub_to_mats, secondary_ip->subdomain_id());
145 
146  const auto & msm_elems = elem_to_msm->second;
147 
148  // Need to be able to check if there's edge dropping, in 3D we can't just compare
149 
150  // Map mortar segment integration points to primary and secondary sides
151  // Note points for segments will be held contiguously to allow reinit without moving
152  // cleaner way to do this would be with contiguously allocated 2D array but would either
153  // need to move into vector for calling
154  std::vector<Point> secondary_xi_pts, primary_xi_pts;
155 
156  std::vector<Real> JxW;
157 
158 #ifndef NDEBUG
159  unsigned int expected_length = 0;
160 #endif
161 
162  // Loop through contributing msm elements
163  for (const auto msm_elem : msm_elems)
164  {
165  // Initialize mortar segment quadrature and compute JxW
166  subproblem.reinitMortarElem(msm_elem, tid);
167 
168  // Get a reference to the MortarSegmentInfo for this Elem.
169  const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
170 
171  if (msm_elem->dim() == 1)
172  {
173  for (unsigned int qp = 0; qp < qrule_msm->n_points(); qp++)
174  {
175  const Real eta = qrule_msm->qp(qp)(0);
176 
177  // Map quadrature points to secondary side
178  const Real xi1_eta = 0.5 * (1 - eta) * msinfo.xi1_a + 0.5 * (1 + eta) * msinfo.xi1_b;
179  secondary_xi_pts.push_back(xi1_eta);
180 
181  // Map quadrature points to primary side
182  const Real xi2_eta = 0.5 * (1 - eta) * msinfo.xi2_a + 0.5 * (1 + eta) * msinfo.xi2_b;
183  primary_xi_pts.push_back(xi2_eta);
184  }
185  }
186  else
187  {
188  if (amg.mortar3DQpMapping() == Mortar3DQuadraturePointMapping::REFERENCE_INTERPOLATION)
189  mapQPoints3dFromReference(*msm_elem,
190  amg.mortarSegmentReferencePoints(*msm_elem),
191  *qrule_msm,
192  secondary_xi_pts,
193  primary_xi_pts);
194  else
195  {
196  // Map independently because the parent-face linearizations differ.
197  projectQPoints3d(msm_elem,
198  msinfo.secondary_elem,
199  secondary_sub_elem_index,
200  *qrule_msm,
201  secondary_xi_pts);
203  msm_elem, msinfo.primary_elem, primary_sub_elem_index, *qrule_msm, primary_xi_pts);
204  }
205  }
206 
207  // If edge dropping case we need JxW on the msm to compute dual shape functions
208  if (assembly.needDual())
209  std::copy(std::begin(JxW_msm), std::end(JxW_msm), std::back_inserter(JxW));
210 
211 #ifndef NDEBUG
212  // Verify that the expected number of quadrature points have been inserted
213  expected_length += qrule_msm->n_points();
214  mooseAssert(secondary_xi_pts.size() == expected_length,
215  "Fewer than expected secondary quadrature points");
216  mooseAssert(primary_xi_pts.size() == expected_length,
217  "Fewer than expected primary quadrature points");
218 
219  if (assembly.needDual())
220  mooseAssert(JxW.size() == expected_length, "Fewer than expected JxW values computed");
221 #endif
222  } // end loop over msm_elems
223 
224  // Reinit dual shape coeffs if dual shape functions needed
225  // lindsayad: is there any need to make sure we do this on both reference and displaced?
226  if (assembly.needDual())
227  assembly.reinitDual(secondary_face_elem, secondary_xi_pts, JxW);
228 
229  unsigned int n_segment = 0;
230 
231  // Loop through contributing msm elements, computing residual and Jacobian this time
232  for (const auto msm_elem : msm_elems)
233  {
234  n_segment++;
235 
236  // These will hold quadrature points for each segment
237  std::vector<Point> xi1_pts, xi2_pts;
238 
239  // Get a reference to the MortarSegmentInfo for this Elem.
240  const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
241 
242  // Set the primary interior parent and side ids
243  const Elem * primary_ip = msinfo.primary_elem->interior_parent();
244  unsigned int primary_side_id = primary_ip->which_side_am_i(msinfo.primary_elem);
245  const auto & primary_ip_mats =
246  libmesh_map_find(primary_ip_sub_to_mats, primary_ip->subdomain_id());
247 
248  // Compute a JxW for the actual mortar segment element (not the lower dimensional element on
249  // the secondary face!)
250  subproblem.reinitMortarElem(msm_elem, tid);
251 
252  // Extract previously computed mapped quadrature points for secondary and primary face
253  // elements
254  const unsigned int start = (n_segment - 1) * qrule_msm->n_points();
255  const unsigned int end = n_segment * qrule_msm->n_points();
256  xi1_pts.insert(
257  xi1_pts.begin(), secondary_xi_pts.begin() + start, secondary_xi_pts.begin() + end);
258  xi2_pts.insert(xi2_pts.begin(), primary_xi_pts.begin() + start, primary_xi_pts.begin() + end);
259 
260  const Elem * reinit_secondary_elem = secondary_ip;
261 
262  // If we're on the displaced mesh, we need to get the corresponding undisplaced elem before
263  // calling fe_problem.reinitElemFaceRef
264  if (displaced)
265  reinit_secondary_elem = fe_problem.mesh().elemPtr(reinit_secondary_elem->id());
266 
267  // NOTE to future developers: it can be tempting to try and change calls on fe_problem to
268  // calls on subproblem because it seems wasteful to reinit data on both the reference and
269  // displaced problems regardless of whether we are running on a "reference" or displaced
270  // mortar mesh. But making such changes opens a can of worms. For instance, a user may define
271  // constant material properties that are used by mortar constraints and not think about
272  // whether they should set `use_displaced_mesh` for the material (and indeed the user may want
273  // those constant properties usable by both reference and displaced consumer objects). If we
274  // reinit that material and we haven't reinit'd it's assembly member, then we will get things
275  // like segmentation faults. Moreover, one can easily imagine that a material may couple in
276  // variables so we need to make sure those are reinit'd too. So even though it's inefficient,
277  // it's safest to keep making calls on fe_problem instead of subproblem
278 
279  // reinit the variables/residuals/jacobians on the secondary interior
280  fe_problem.reinitElemFaceRef(
281  reinit_secondary_elem, secondary_side_id, TOLERANCE, &xi1_pts, nullptr, tid);
282 
283  const Elem * reinit_primary_elem = primary_ip;
284 
285  // If we're on the displaced mesh, we need to get the corresponding undisplaced elem before
286  // calling fe_problem.reinitElemFaceRef
287  if (displaced)
288  reinit_primary_elem = fe_problem.mesh().elemPtr(reinit_primary_elem->id());
289 
290  // reinit the variables/residuals/jacobians on the primary interior
291  fe_problem.reinitNeighborFaceRef(
292  reinit_primary_elem, primary_side_id, TOLERANCE, &xi2_pts, nullptr, tid);
293 
294  // reinit neighbor materials, but be careful not to execute stateful materials since
295  // conceptually they don't make sense with mortar (they're not interpolary)
296  fe_problem.reinitMaterialsNeighbor(primary_ip->subdomain_id(),
297  /*tid=*/tid,
298  /*swap_stateful=*/false,
299  &primary_ip_mats);
300 
301  // reinit the variables/residuals/jacobians on the lower dimensional element corresponding to
302  // the secondary face. This must be done last after the dof indices have been prepared for the
303  // secondary (element) and primary (neighbor)
304  subproblem.reinitLowerDElem(secondary_face_elem, /*tid=*/tid, &xi1_pts);
305 
306  // All this does currently is sets the neighbor/primary lower dimensional elem in Assembly and
307  // computes its volume for potential use in the MortarConstraints. Solution continuity
308  // stabilization for example relies on being able to access the volume
309  subproblem.reinitNeighborLowerDElem(msinfo.primary_elem, tid);
310 
311  // reinit higher-dimensional secondary face/boundary materials. Do this after we reinit
312  // lower-d variables in case we want to pull the lower-d variable values into the secondary
313  // face/boundary materials. Be careful not to execute stateful materials since conceptually
314  // they don't make sense with mortar (they're not interpolary)
315  fe_problem.reinitMaterialsFace(secondary_ip->subdomain_id(),
316  /*tid=*/tid,
317  /*swap_stateful=*/false,
318  &secondary_ip_mats);
319  fe_problem.reinitMaterialsBoundary(
320  secondary_boundary_id, /*tid=*/tid, /*swap_stateful=*/false, &secondary_boundary_mats);
321 
322  if (reinit_mortar_user_objects)
323  fe_problem.reinitMortarUserObjects(primary_boundary_id, secondary_boundary_id, displaced);
324 
325  act();
326 
327  } // End loop over msm segments on secondary face elem
328  } // End loop over (active) secondary elems
329 }
virtual MooseMesh & mesh()=0
void setActiveMaterialProperties(const std::unordered_set< unsigned int > &mat_prop_ids, const THREAD_ID tid)
Record and set the material properties required by the current computing thread.
void reinitNeighborLowerDElem(const Elem *elem, const THREAD_ID tid=0)
reinitialize a neighboring lower dimensional element
Definition: SubProblem.C:988
const std::pair< BoundaryID, BoundaryID > & primarySecondaryBoundaryIDPair() const
virtual void reinitLowerDElem(const Elem *lower_d_elem, const THREAD_ID tid, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)
Definition: SubProblem.C:958
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:3213
void reinitMortarUserObjects(BoundaryID primary_boundary_id, BoundaryID secondary_boundary_id, bool displaced)
Call reinit on mortar user objects with matching primary boundary ID, secondary boundary ID...
void reinitMaterialsBoundary(BoundaryID boundary_id, const THREAD_ID tid, bool swap_stateful=true, const std::deque< MaterialBase *> *reinit_mats=nullptr)
reinit materials on a boundary
const libMesh::QBase *const & qRuleMortar() const
Returns a reference to the quadrature rule for the mortar segments.
Definition: Assembly.h:712
void mapQPoints3dFromReference(const Elem &mortar_segment_elem, const MortarSegmentReferencePoints &reference_points, const QBase &qrule_msm, std::vector< Point > &secondary_q_pts, std::vector< Point > &primary_q_pts)
3D mapping operator that interpolates stored parent reference points on each triangular mortar segmen...
Definition: MortarUtils.C:101
void reinitMaterialsFace(SubdomainID blk_id, const THREAD_ID tid, bool swap_stateful=true, const std::deque< MaterialBase *> *reinit_mats=nullptr)
reinit materials on element faces
virtual void reinitElemFaceRef(const Elem *elem, unsigned int side, Real tolerance, const std::vector< Point > *const pts, const std::vector< Real > *const weights=nullptr, const THREAD_ID tid=0) override
reinitialize FE objects on a given element on a given side at a given set of reference points and the...
const Elem * primary_elem
bool needDual() const
Indicates whether dual shape functions are used (computation is now repeated on each element so expen...
Definition: Assembly.h:623
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3548
void reinitMaterialsNeighbor(SubdomainID blk_id, const THREAD_ID tid, bool swap_stateful=true, const std::deque< MaterialBase *> *reinit_mats=nullptr)
reinit materials on the neighboring element face
const MortarSegmentReferencePoints & mortarSegmentReferencePoints(const Elem &mortar_segment_elem) const
Return the parent-face reference coordinates for a mortar segment.
virtual void reinitNeighborFaceRef(const Elem *neighbor_elem, unsigned int neighbor_side, Real tolerance, const std::vector< Point > *const pts, const std::vector< Real > *const weights=nullptr, const THREAD_ID tid=0) override
reinitialize FE objects on a given neighbor element on a given side at a given set of reference point...
const Elem * secondary_elem
Mortar3DQuadraturePointMapping mortar3DQpMapping() const
Return the 3D mortar quadrature-point mapping method.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Holds xi^(1), xi^(2), and other data for a given mortar segment.
const std::unordered_map< const Elem *, MortarSegmentInfo > & mortarSegmentMeshElemToInfo() const
virtual MooseMesh & mesh() override
const MeshBase & mortarSegmentMesh() const
const std::vector< Real > & jxWMortar() const
Returns a reference to JxW for mortar segment elements.
Definition: Assembly.h:707
void reinitMortarElem(const Elem *elem, const THREAD_ID tid=0)
Reinit a mortar element to obtain a valid JxW.
Definition: SubProblem.C:995
void projectQPoints3d(const Elem *msm_elem, const Elem *primal_elem, unsigned int sub_elem_index, const QBase &qrule_msm, std::vector< Point > &q_pts)
3D projection operator for mapping qpoints on mortar segments to secondary or primary elements ...
Definition: MortarUtils.C:130
void reinitDual(const Elem *elem, const std::vector< Point > &pts, const std::vector< Real > &JxW)
Reintialize dual basis coefficients based on a customized quadrature rule.
Definition: Assembly.C:2274

◆ mapQPoints3dFromReference()

void Moose::Mortar::mapQPoints3dFromReference ( const Elem &  mortar_segment_elem,
const MortarSegmentReferencePoints reference_points,
const QBase &  qrule_msm,
std::vector< Point > &  secondary_q_pts,
std::vector< Point > &  primary_q_pts 
)

3D mapping operator that interpolates stored parent reference points on each triangular mortar segment.

Parameters
mortar_segment_elemThe triangular mortar segment carrying the quadrature rule
reference_pointsParent-face reference points stored for the mortar segment
qrule_msmThe rule that governs quadrature on the mortar segment element
secondary_q_ptsReference-space quadrature points on the secondary face
primary_q_ptsReference-space quadrature points on the primary face

Definition at line 101 of file MortarUtils.C.

Referenced by loopOverMortarSegments().

106 {
107  mooseAssert(mortar_segment_elem.type() == TRI3,
108  "Reference interpolation expects triangular mortar segments.");
109  const FEType fe_type(FIRST, LAGRANGE);
110 
111  for (const auto qp : make_range(qrule_msm.n_points()))
112  {
113  Point secondary_qp;
114  Point primary_qp;
115 
116  for (const auto n : index_range(reference_points.secondary_reference_points))
117  {
118  const auto phi =
119  FEInterface::shape(fe_type, &mortar_segment_elem, n, qrule_msm.qp(qp), false);
120  secondary_qp += phi * reference_points.secondary_reference_points[n];
121  primary_qp += phi * reference_points.primary_reference_points[n];
122  }
123 
124  secondary_q_pts.push_back(secondary_qp);
125  primary_q_pts.push_back(primary_qp);
126  }
127 }
std::array< Point, 3 > secondary_reference_points
TRI3
std::array< Point, 3 > primary_reference_points
IntRange< T > make_range(T beg, T end)
auto index_range(const T &sizable)

◆ projectQPoints3d()

void Moose::Mortar::projectQPoints3d ( const Elem *  msm_elem,
const Elem *  primal_elem,
unsigned int  sub_elem_index,
const QBase &  qrule_msm,
std::vector< Point > &  q_pts 
)

3D projection operator for mapping qpoints on mortar segments to secondary or primary elements

Parameters
msm_elemThe mortar segment element that we will be mapping quadrature points from
primal_elemThe "persistent" mesh element (e.g. it exists on the simulation's MooseMesh) that we will be mapping quadrature points for. This can be either an element on the secondary or primary face
sub_elem_indexWe will call msm_elem->get_extra_integer(sub_elem_index) in the implementation in order to determine which sub-element of the primal element the mortar segment element corresponds to. This sub_elem_index should correspond to the secondary element index if primal_elem is a secondary face element and the primary element index if primal_elem is a primary face element
qrule_msmThe rule that governs quadrature on the mortar segment element
q_ptsThe output of this function. This will correspond to the the (reference space) quadrature points that we wish to evaluate shape functions, etc., at on the primal element

Definition at line 130 of file MortarUtils.C.

Referenced by loopOverMortarSegments().

135 {
136  const auto msm_elem_order = msm_elem->default_order();
137  const auto msm_elem_type = msm_elem->type();
138 
139  // Get normal to linearized element, could store and query but computation is easy
140  const Point e1 = msm_elem->point(0) - msm_elem->point(1);
141  const Point e2 = msm_elem->point(2) - msm_elem->point(1);
142  const Point normal = e2.cross(e1).unit();
143 
144  // Get sub-elem (for second order meshes, otherwise trivial)
145  const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
146  const ElemType primal_type = primal_elem->type();
147 
148  // Transforms quadrature point from first order sub-elements (in case of second-order)
149  // to primal element
150  auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
151  {
152  switch (primal_type)
153  {
154  case TRI3:
155  return Point(nu, xi, 0);
156  case QUAD4:
157  return Point(nu, xi, 0);
158  case TRI6:
159  case TRI7:
160  switch (sub_elem)
161  {
162  case 0:
163  return Point(0.5 * nu, 0.5 * xi, 0);
164  case 1:
165  return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
166  case 2:
167  return Point(0.5 * (1 + nu), 0.5 * xi, 0);
168  case 3:
169  return Point(0.5 * nu, 0.5 * (1 + xi), 0);
170  default:
171  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
172  }
173  case QUAD8:
174  switch (sub_elem)
175  {
176  case 0:
177  return Point(nu - 1, xi - 1, 0);
178  case 1:
179  return Point(nu + xi, xi - 1, 0);
180  case 2:
181  return Point(1 - xi, nu + xi, 0);
182  case 3:
183  return Point(nu - 1, nu + xi, 0);
184  case 4:
185  return Point(0.5 * (nu - xi), 0.5 * (nu + xi), 0);
186  default:
187  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
188  }
189  case QUAD9:
190  switch (sub_elem)
191  {
192  case 0:
193  return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
194  case 1:
195  return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
196  case 2:
197  return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
198  case 3:
199  return Point(0.5 * (nu - 1), 0.5 * (xi + 1), 0);
200  default:
201  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
202  }
203  default:
204  mooseError("transform_qp: Face element type: ",
205  libMesh::Utility::enum_to_string<ElemType>(primal_type),
206  " invalid for 3D mortar");
207  }
208  };
209 
210  auto sub_element_type = [primal_type, sub_elem]()
211  {
212  switch (primal_type)
213  {
214  case TRI3:
215  case TRI6:
216  case TRI7:
217  return TRI3;
218  case QUAD4:
219  case QUAD9:
220  return QUAD4;
221  case QUAD8:
222  switch (sub_elem)
223  {
224  case 0:
225  case 1:
226  case 2:
227  case 3:
228  return TRI3;
229  case 4:
230  return QUAD4;
231  default:
232  mooseError("sub_element_type: Invalid sub_elem: ", sub_elem);
233  }
234  default:
235  mooseError("sub_element_type: Face element type: ",
236  libMesh::Utility::enum_to_string<ElemType>(primal_type),
237  " invalid for 3D mortar");
238  }
239  };
240 
241  // Get sub-elem node indices
242  const auto sub_elem_node_indices = getMortarSubElementNodeIndices(*primal_elem, sub_elem);
243 
244  // Loop through quadrature points on msm_elem
245  for (auto qp : make_range(qrule_msm.n_points()))
246  {
247  // Get physical point on msm_elem to project
248  Point x0;
249  for (auto n : make_range(msm_elem->n_nodes()))
250  x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
251  msm_elem_order,
252  n,
253  static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
254  msm_elem->point(n);
255 
256  // Use msm_elem quadrature point as initial guess
257  // (will be correct for aligned meshes)
258  Dual2 xi1{};
259  xi1.value() = qrule_msm.qp(qp)(0);
260  xi1.derivatives()[0] = 1.0;
261  Dual2 xi2{};
262  xi2.value() = qrule_msm.qp(qp)(1);
263  xi2.derivatives()[1] = 1.0;
264  VectorValue<Dual2> xi(xi1, xi2, 0);
265  unsigned int current_iterate = 0, max_iterates = 10;
266 
267  // Project qp from mortar segments to first order sub-elements (elements in case of first order
268  // geometry)
269  do
270  {
271  VectorValue<Dual2> x1;
272  for (auto n : make_range(sub_elem_node_indices.size()))
273  x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
274  primal_elem->point(sub_elem_node_indices[n]);
275  auto u = x1 - x0;
276 
277  VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
278  u(2) * normal(0) - u(0) * normal(2),
279  u(0) * normal(1) - u(1) * normal(0));
280 
281  Real projection_tolerance(1e-10);
282 
283  // Normalize tolerance with quantities involved in the projection.
284  // Absolute projection tolerance is loosened for displacements larger than those on the order
285  // of one. Tightening the tolerance for displacements of smaller orders causes this tolerance
286  // to not be reached in a number of tests.
287  if (!u.is_zero() && u.norm().value() > 1.0)
288  projection_tolerance *= u.norm().value();
289 
290  if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
291  break;
292 
293  RealEigenMatrix J(3, 2);
294  J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
295  F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
296  RealEigenVector f(3);
297  f << F(0).value(), F(1).value(), F(2).value();
298  const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
299 
300  xi(0) += dxi(0);
301  xi(1) += dxi(1);
302  } while (++current_iterate < max_iterates);
303 
304  if (current_iterate < max_iterates)
305  {
306  // Transfer quadrature point from sub-element to element and store
307  q_pts.push_back(transform_qp(xi(0).value(), xi(1).value()));
308 
309  // The following checks if quadrature point falls in correct domain.
310  // On small mortar segment elements with very distorted elements this can fail, instead of
311  // erroring simply truncate quadrature point, these points typically have very small
312  // contributions to integrals
313  auto & qp_back = q_pts.back();
314  if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
315  {
316  if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
317  qp_back(0) + qp_back(1) > (1 + TOLERANCE))
318  mooseException("Quadrature point: ", qp_back, " out of bounds, truncating.");
319  }
320  else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
321  primal_elem->type() == QUAD9)
322  {
323  if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
324  qp_back(1) < (-1 - TOLERANCE) || qp_back(1) > (1 + TOLERANCE))
325  mooseException("Quadrature point: ", qp_back, " out of bounds, truncating");
326  }
327  }
328  else
329  {
330  mooseError("Newton iteration for mortar quadrature mapping msm element: ",
331  msm_elem->id(),
332  " to elem: ",
333  primal_elem->id(),
334  " didn't converge. MSM element volume: ",
335  msm_elem->volume());
336  }
337  }
338 }
T fe_lagrange_2D_shape(const libMesh::ElemType type, const Order order, const unsigned int i, const VectorType< T > &p)
ElemType
QUAD8
DualNumber< Real, NumberArray< 2, Real > > Dual2
Definition: MortarUtils.C:20
std::vector< unsigned int > getMortarSubElementNodeIndices(const Elem &parent_elem, unsigned int sub_elem)
Return the node indices for a first-order sub-element of a parent face.
Definition: MortarUtils.C:27
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:311
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:100
TRI3
QUAD4
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
TRI6
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > RealEigenMatrix
Definition: MooseTypes.h:151
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
auto norm(const T &a)
IntRange< T > make_range(T beg, T end)
TRI7
QUAD9
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
Definition: MooseTypes.h:147

◆ setupMortarMaterials()

template<typename Consumers >
void Moose::Mortar::setupMortarMaterials ( const Consumers &  consumers,
FEProblemBase fe_problem,
const AutomaticMortarGeneration amg,
const THREAD_ID  tid,
std::map< SubdomainID, std::deque< MaterialBase *>> &  secondary_ip_sub_to_mats,
std::map< SubdomainID, std::deque< MaterialBase *>> &  primary_ip_sub_to_mats,
std::deque< MaterialBase *> &  secondary_boundary_mats 
)

This function creates containers of materials necessary to execute the mortar method for a supplied set of consumers.

Parameters
consumersThe objects that we're building the material dependencies for. This could be a container of mortar constraints or a "mortar" auxiliary kernel for example
fe_problemThe finite element problem that we'll be querying for material warehouses
tidThe thread ID that we will use for pulling the material warehouse
secondary_ip_sub_to_matsA map from the secondary interior parent subdomain IDs to the required secondary block materials we will need to evaluate for the consumers
primary_ip_sub_to_matsA map from the primary interior parent subdomain IDs to the required primary block materials we will need to evaluate for the consumers
secondary_boundary_matsThe secondary boundary materials we will need to evaluate for the consumers

Definition at line 347 of file MortarUtils.h.

Referenced by ComputeMortarFunctor::ComputeMortarFunctor(), MortarNodalAuxKernelTempl< ComputeValueType >::initialSetup(), MortarUserObjectThread::MortarUserObjectThread(), and ComputeMortarFunctor::setupMortarMaterials().

354 {
355  secondary_ip_sub_to_mats.clear();
356  primary_ip_sub_to_mats.clear();
357  secondary_boundary_mats.clear();
358 
359  auto & mat_warehouse = fe_problem.getRegularMaterialsWarehouse();
360  auto get_required_sub_mats =
361  [&mat_warehouse, tid, &consumers](
362  const SubdomainID sub_id,
363  const Moose::MaterialDataType mat_data_type) -> std::deque<MaterialBase *>
364  {
365  if (mat_warehouse[mat_data_type].hasActiveBlockObjects(sub_id, tid))
366  {
367  auto & sub_mats = mat_warehouse[mat_data_type].getActiveBlockObjects(sub_id, tid);
368  return MaterialBase::buildRequiredMaterials(consumers, sub_mats, /*allow_stateful=*/false);
369  }
370  else
371  return {};
372  };
373 
374  // Construct secondary *block* materials container
375  const auto & secondary_ip_sub_ids = amg.secondaryIPSubIDs();
376  for (const auto secondary_ip_sub : secondary_ip_sub_ids)
377  secondary_ip_sub_to_mats.emplace(
378  secondary_ip_sub, get_required_sub_mats(secondary_ip_sub, Moose::FACE_MATERIAL_DATA));
379 
380  // Construct primary *block* materials container
381  const auto & primary_ip_sub_ids = amg.primaryIPSubIDs();
382  for (const auto primary_ip_sub : primary_ip_sub_ids)
383  primary_ip_sub_to_mats.emplace(
384  primary_ip_sub, get_required_sub_mats(primary_ip_sub, Moose::NEIGHBOR_MATERIAL_DATA));
385 
386  // Construct secondary *boundary* materials container
387  const auto & boundary_pr = amg.primarySecondaryBoundaryIDPair();
388  const auto secondary_boundary = boundary_pr.second;
389  if (mat_warehouse.hasActiveBoundaryObjects(secondary_boundary, tid))
390  {
391  auto & boundary_mats = mat_warehouse.getActiveBoundaryObjects(secondary_boundary, tid);
392  secondary_boundary_mats =
393  MaterialBase::buildRequiredMaterials(consumers, boundary_mats, /*allow_stateful=*/false);
394  }
395 }
const std::pair< BoundaryID, BoundaryID > & primarySecondaryBoundaryIDPair() const
const std::map< SubdomainID, std::vector< std::shared_ptr< T > > > & getActiveBlockObjects(THREAD_ID tid=0) const
MaterialDataType
MaterialData types.
Definition: MooseTypes.h:740
const std::set< SubdomainID > & primaryIPSubIDs() const
const MaterialWarehouse & getRegularMaterialsWarehouse() const
const std::set< SubdomainID > & secondaryIPSubIDs() const
static std::deque< MaterialBase * > buildRequiredMaterials(const Consumers &mat_consumers, const std::vector< std::shared_ptr< MaterialBase >> &mats, const bool allow_stateful)
Build the materials required by a set of consumer objects.
Definition: MaterialBase.h:538