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

Functions

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...
 
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

◆ 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 70 of file MortarUtils.h.

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

84 {
85  const auto & primary_secondary_boundary_id_pair = amg.primarySecondaryBoundaryIDPair();
86 
87  const auto primary_boundary_id = primary_secondary_boundary_id_pair.first;
88  const auto secondary_boundary_id = primary_secondary_boundary_id_pair.second;
89 
90  // For 3D mortar get index for retrieving sub-element info
91  unsigned int secondary_sub_elem_index = 0, primary_sub_elem_index = 0;
92  if (amg.dim() == 3)
93  {
94  secondary_sub_elem_index = amg.mortarSegmentMesh().get_elem_integer_index("secondary_sub_elem");
95  primary_sub_elem_index = amg.mortarSegmentMesh().get_elem_integer_index("primary_sub_elem");
96  }
97 
98  // The mortar quadrature rule. Necessary for sizing the number of custom points for re-init'ing
99  // the secondary interior, primary interior, and secondary face elements
100  const auto & qrule_msm = assembly.qRuleMortar();
101 
102  // The element Jacobian times weights
103  const auto & JxW_msm = assembly.jxWMortar();
104 
105  // Set required material properties
106  std::unordered_set<unsigned int> needed_mat_props;
107  for (const auto & consumer : consumers)
108  {
109  const auto & mp_deps = consumer->getMatPropDependencies();
110  needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
111  }
112  fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/tid);
113 
114  // Loop through secondary elements, accumulating quadrature points for all corresponding mortar
115  // segments
116  for (const auto elem_to_msm : secondary_elems_to_mortar_segments)
117  {
118  const Elem * secondary_face_elem = subproblem.mesh().getMesh().elem_ptr(elem_to_msm->first);
119  // Set the secondary interior parent and side ids
120  const Elem * secondary_ip = secondary_face_elem->interior_parent();
121  unsigned int secondary_side_id = secondary_ip->which_side_am_i(secondary_face_elem);
122  const auto & secondary_ip_mats =
123  libmesh_map_find(secondary_ip_sub_to_mats, secondary_ip->subdomain_id());
124 
125  const auto & msm_elems = elem_to_msm->second;
126 
127  // Need to be able to check if there's edge dropping, in 3D we can't just compare
128 
129  // Map mortar segment integration points to primary and secondary sides
130  // Note points for segments will be held contiguously to allow reinit without moving
131  // cleaner way to do this would be with contiguously allocated 2D array but would either
132  // need to move into vector for calling
133  std::vector<Point> secondary_xi_pts, primary_xi_pts;
134 
135  std::vector<Real> JxW;
136 
137 #ifndef NDEBUG
138  unsigned int expected_length = 0;
139 #endif
140 
141  // Loop through contributing msm elements
142  for (const auto msm_elem : msm_elems)
143  {
144  // Initialize mortar segment quadrature and compute JxW
145  subproblem.reinitMortarElem(msm_elem, tid);
146 
147  // Get a reference to the MortarSegmentInfo for this Elem.
148  const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
149 
150  if (msm_elem->dim() == 1)
151  {
152  for (unsigned int qp = 0; qp < qrule_msm->n_points(); qp++)
153  {
154  const Real eta = qrule_msm->qp(qp)(0);
155 
156  // Map quadrature points to secondary side
157  const Real xi1_eta = 0.5 * (1 - eta) * msinfo.xi1_a + 0.5 * (1 + eta) * msinfo.xi1_b;
158  secondary_xi_pts.push_back(xi1_eta);
159 
160  // Map quadrature points to primary side
161  const Real xi2_eta = 0.5 * (1 - eta) * msinfo.xi2_a + 0.5 * (1 + eta) * msinfo.xi2_b;
162  primary_xi_pts.push_back(xi2_eta);
163  }
164  }
165  else
166  {
167  // Now that has_primary logic gone, could combine these to make more efficient
168  projectQPoints3d(msm_elem,
169  msinfo.secondary_elem,
170  secondary_sub_elem_index,
171  *qrule_msm,
172  secondary_xi_pts);
174  msm_elem, msinfo.primary_elem, primary_sub_elem_index, *qrule_msm, primary_xi_pts);
175  }
176 
177  // If edge dropping case we need JxW on the msm to compute dual shape functions
178  if (assembly.needDual())
179  std::copy(std::begin(JxW_msm), std::end(JxW_msm), std::back_inserter(JxW));
180 
181 #ifndef NDEBUG
182  // Verify that the expected number of quadrature points have been inserted
183  expected_length += qrule_msm->n_points();
184  mooseAssert(secondary_xi_pts.size() == expected_length,
185  "Fewer than expected secondary quadrature points");
186  mooseAssert(primary_xi_pts.size() == expected_length,
187  "Fewer than expected primary quadrature points");
188 
189  if (assembly.needDual())
190  mooseAssert(JxW.size() == expected_length, "Fewer than expected JxW values computed");
191 #endif
192  } // end loop over msm_elems
193 
194  // Reinit dual shape coeffs if dual shape functions needed
195  // lindsayad: is there any need to make sure we do this on both reference and displaced?
196  if (assembly.needDual())
197  assembly.reinitDual(secondary_face_elem, secondary_xi_pts, JxW);
198 
199  unsigned int n_segment = 0;
200 
201  // Loop through contributing msm elements, computing residual and Jacobian this time
202  for (const auto msm_elem : msm_elems)
203  {
204  n_segment++;
205 
206  // These will hold quadrature points for each segment
207  std::vector<Point> xi1_pts, xi2_pts;
208 
209  // Get a reference to the MortarSegmentInfo for this Elem.
210  const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
211 
212  // Set the primary interior parent and side ids
213  const Elem * primary_ip = msinfo.primary_elem->interior_parent();
214  unsigned int primary_side_id = primary_ip->which_side_am_i(msinfo.primary_elem);
215  const auto & primary_ip_mats =
216  libmesh_map_find(primary_ip_sub_to_mats, primary_ip->subdomain_id());
217 
218  // Compute a JxW for the actual mortar segment element (not the lower dimensional element on
219  // the secondary face!)
220  subproblem.reinitMortarElem(msm_elem, tid);
221 
222  // Extract previously computed mapped quadrature points for secondary and primary face
223  // elements
224  const unsigned int start = (n_segment - 1) * qrule_msm->n_points();
225  const unsigned int end = n_segment * qrule_msm->n_points();
226  xi1_pts.insert(
227  xi1_pts.begin(), secondary_xi_pts.begin() + start, secondary_xi_pts.begin() + end);
228  xi2_pts.insert(xi2_pts.begin(), primary_xi_pts.begin() + start, primary_xi_pts.begin() + end);
229 
230  const Elem * reinit_secondary_elem = secondary_ip;
231 
232  // If we're on the displaced mesh, we need to get the corresponding undisplaced elem before
233  // calling fe_problem.reinitElemFaceRef
234  if (displaced)
235  reinit_secondary_elem = fe_problem.mesh().elemPtr(reinit_secondary_elem->id());
236 
237  // NOTE to future developers: it can be tempting to try and change calls on fe_problem to
238  // calls on subproblem because it seems wasteful to reinit data on both the reference and
239  // displaced problems regardless of whether we are running on a "reference" or displaced
240  // mortar mesh. But making such changes opens a can of worms. For instance, a user may define
241  // constant material properties that are used by mortar constraints and not think about
242  // whether they should set `use_displaced_mesh` for the material (and indeed the user may want
243  // those constant properties usable by both reference and displaced consumer objects). If we
244  // reinit that material and we haven't reinit'd it's assembly member, then we will get things
245  // like segmentation faults. Moreover, one can easily imagine that a material may couple in
246  // variables so we need to make sure those are reinit'd too. So even though it's inefficient,
247  // it's safest to keep making calls on fe_problem instead of subproblem
248 
249  // reinit the variables/residuals/jacobians on the secondary interior
250  fe_problem.reinitElemFaceRef(
251  reinit_secondary_elem, secondary_side_id, TOLERANCE, &xi1_pts, nullptr, tid);
252 
253  const Elem * reinit_primary_elem = primary_ip;
254 
255  // If we're on the displaced mesh, we need to get the corresponding undisplaced elem before
256  // calling fe_problem.reinitElemFaceRef
257  if (displaced)
258  reinit_primary_elem = fe_problem.mesh().elemPtr(reinit_primary_elem->id());
259 
260  // reinit the variables/residuals/jacobians on the primary interior
261  fe_problem.reinitNeighborFaceRef(
262  reinit_primary_elem, primary_side_id, TOLERANCE, &xi2_pts, nullptr, tid);
263 
264  // reinit neighbor materials, but be careful not to execute stateful materials since
265  // conceptually they don't make sense with mortar (they're not interpolary)
266  fe_problem.reinitMaterialsNeighbor(primary_ip->subdomain_id(),
267  /*tid=*/tid,
268  /*swap_stateful=*/false,
269  &primary_ip_mats);
270 
271  // reinit the variables/residuals/jacobians on the lower dimensional element corresponding to
272  // the secondary face. This must be done last after the dof indices have been prepared for the
273  // secondary (element) and primary (neighbor)
274  subproblem.reinitLowerDElem(secondary_face_elem, /*tid=*/tid, &xi1_pts);
275 
276  // All this does currently is sets the neighbor/primary lower dimensional elem in Assembly and
277  // computes its volume for potential use in the MortarConstraints. Solution continuity
278  // stabilization for example relies on being able to access the volume
279  subproblem.reinitNeighborLowerDElem(msinfo.primary_elem, tid);
280 
281  // reinit higher-dimensional secondary face/boundary materials. Do this after we reinit
282  // lower-d variables in case we want to pull the lower-d variable values into the secondary
283  // face/boundary materials. Be careful not to execute stateful materials since conceptually
284  // they don't make sense with mortar (they're not interpolary)
285  fe_problem.reinitMaterialsFace(secondary_ip->subdomain_id(),
286  /*tid=*/tid,
287  /*swap_stateful=*/false,
288  &secondary_ip_mats);
289  fe_problem.reinitMaterialsBoundary(
290  secondary_boundary_id, /*tid=*/tid, /*swap_stateful=*/false, &secondary_boundary_mats);
291 
292  if (reinit_mortar_user_objects)
293  fe_problem.reinitMortarUserObjects(primary_boundary_id, secondary_boundary_id, displaced);
294 
295  act();
296 
297  } // End loop over msm segments on secondary face elem
298  } // End loop over (active) secondary elems
299 }
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:3134
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 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:3469
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
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
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:26
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:2270

◆ 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 26 of file MortarUtils.C.

Referenced by loopOverMortarSegments().

31 {
32  const auto msm_elem_order = msm_elem->default_order();
33  const auto msm_elem_type = msm_elem->type();
34 
35  // Get normal to linearized element, could store and query but computation is easy
36  const Point e1 = msm_elem->point(0) - msm_elem->point(1);
37  const Point e2 = msm_elem->point(2) - msm_elem->point(1);
38  const Point normal = e2.cross(e1).unit();
39 
40  // Get sub-elem (for second order meshes, otherwise trivial)
41  const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
42  const ElemType primal_type = primal_elem->type();
43 
44  auto get_sub_elem_node_indices = [primal_type, sub_elem]() -> std::vector<unsigned int>
45  {
46  switch (primal_type)
47  {
48  case TRI3:
49  return {{0, 1, 2}};
50  case QUAD4:
51  return {{0, 1, 2, 3}};
52  case TRI6:
53  case TRI7:
54  switch (sub_elem)
55  {
56  case 0:
57  return {{0, 3, 5}};
58  case 1:
59  return {{3, 4, 5}};
60  case 2:
61  return {{3, 1, 4}};
62  case 3:
63  return {{5, 4, 2}};
64  default:
65  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
66  }
67  case QUAD8:
68  switch (sub_elem)
69  {
70  case 0:
71  return {{0, 4, 7}};
72  case 1:
73  return {{4, 1, 5}};
74  case 2:
75  return {{5, 2, 6}};
76  case 3:
77  return {{7, 6, 3}};
78  case 4:
79  return {{4, 5, 6, 7}};
80  default:
81  mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
82  }
83  case QUAD9:
84  switch (sub_elem)
85  {
86  case 0:
87  return {{0, 4, 8, 7}};
88  case 1:
89  return {{4, 1, 5, 8}};
90  case 2:
91  return {{8, 5, 2, 6}};
92  case 3:
93  return {{7, 8, 6, 3}};
94  default:
95  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
96  }
97  default:
98  mooseError("get_sub_elem_indices: Face element type: ",
99  libMesh::Utility::enum_to_string<ElemType>(primal_type),
100  " invalid for 3D mortar");
101  }
102  };
103 
104  // Transforms quadrature point from first order sub-elements (in case of second-order)
105  // to primal element
106  auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
107  {
108  switch (primal_type)
109  {
110  case TRI3:
111  return Point(nu, xi, 0);
112  case QUAD4:
113  return Point(nu, xi, 0);
114  case TRI6:
115  case TRI7:
116  switch (sub_elem)
117  {
118  case 0:
119  return Point(0.5 * nu, 0.5 * xi, 0);
120  case 1:
121  return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
122  case 2:
123  return Point(0.5 * (1 + nu), 0.5 * xi, 0);
124  case 3:
125  return Point(0.5 * nu, 0.5 * (1 + xi), 0);
126  default:
127  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
128  }
129  case QUAD8:
130  switch (sub_elem)
131  {
132  case 0:
133  return Point(nu - 1, xi - 1, 0);
134  case 1:
135  return Point(nu + xi, xi - 1, 0);
136  case 2:
137  return Point(1 - xi, nu + xi, 0);
138  case 3:
139  return Point(nu - 1, nu + xi, 0);
140  case 4:
141  return Point(0.5 * (nu - xi), 0.5 * (nu + xi), 0);
142  default:
143  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
144  }
145  case QUAD9:
146  switch (sub_elem)
147  {
148  case 0:
149  return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
150  case 1:
151  return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
152  case 2:
153  return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
154  case 3:
155  return Point(0.5 * (nu - 1), 0.5 * (xi + 1), 0);
156  default:
157  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
158  }
159  default:
160  mooseError("transform_qp: Face element type: ",
161  libMesh::Utility::enum_to_string<ElemType>(primal_type),
162  " invalid for 3D mortar");
163  }
164  };
165 
166  auto sub_element_type = [primal_type, sub_elem]()
167  {
168  switch (primal_type)
169  {
170  case TRI3:
171  case TRI6:
172  case TRI7:
173  return TRI3;
174  case QUAD4:
175  case QUAD9:
176  return QUAD4;
177  case QUAD8:
178  switch (sub_elem)
179  {
180  case 0:
181  case 1:
182  case 2:
183  case 3:
184  return TRI3;
185  case 4:
186  return QUAD4;
187  default:
188  mooseError("sub_element_type: Invalid sub_elem: ", sub_elem);
189  }
190  default:
191  mooseError("sub_element_type: Face element type: ",
192  libMesh::Utility::enum_to_string<ElemType>(primal_type),
193  " invalid for 3D mortar");
194  }
195  };
196 
197  // Get sub-elem node indices
198  auto sub_elem_node_indices = get_sub_elem_node_indices();
199 
200  // Loop through quadrature points on msm_elem
201  for (auto qp : make_range(qrule_msm.n_points()))
202  {
203  // Get physical point on msm_elem to project
204  Point x0;
205  for (auto n : make_range(msm_elem->n_nodes()))
206  x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
207  msm_elem_order,
208  n,
209  static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
210  msm_elem->point(n);
211 
212  // Use msm_elem quadrature point as initial guess
213  // (will be correct for aligned meshes)
214  Dual2 xi1{};
215  xi1.value() = qrule_msm.qp(qp)(0);
216  xi1.derivatives()[0] = 1.0;
217  Dual2 xi2{};
218  xi2.value() = qrule_msm.qp(qp)(1);
219  xi2.derivatives()[1] = 1.0;
220  VectorValue<Dual2> xi(xi1, xi2, 0);
221  unsigned int current_iterate = 0, max_iterates = 10;
222 
223  // Project qp from mortar segments to first order sub-elements (elements in case of first order
224  // geometry)
225  do
226  {
227  VectorValue<Dual2> x1;
228  for (auto n : make_range(sub_elem_node_indices.size()))
229  x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
230  primal_elem->point(sub_elem_node_indices[n]);
231  auto u = x1 - x0;
232 
233  VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
234  u(2) * normal(0) - u(0) * normal(2),
235  u(0) * normal(1) - u(1) * normal(0));
236 
237  Real projection_tolerance(1e-10);
238 
239  // Normalize tolerance with quantities involved in the projection.
240  // Absolute projection tolerance is loosened for displacements larger than those on the order
241  // of one. Tightening the tolerance for displacements of smaller orders causes this tolerance
242  // to not be reached in a number of tests.
243  if (!u.is_zero() && u.norm().value() > 1.0)
244  projection_tolerance *= u.norm().value();
245 
246  if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
247  break;
248 
249  RealEigenMatrix J(3, 2);
250  J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
251  F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
252  RealEigenVector f(3);
253  f << F(0).value(), F(1).value(), F(2).value();
254  const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
255 
256  xi(0) += dxi(0);
257  xi(1) += dxi(1);
258  } while (++current_iterate < max_iterates);
259 
260  if (current_iterate < max_iterates)
261  {
262  // Transfer quadrature point from sub-element to element and store
263  q_pts.push_back(transform_qp(xi(0).value(), xi(1).value()));
264 
265  // The following checks if quadrature point falls in correct domain.
266  // On small mortar segment elements with very distorted elements this can fail, instead of
267  // erroring simply truncate quadrature point, these points typically have very small
268  // contributions to integrals
269  auto & qp_back = q_pts.back();
270  if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
271  {
272  if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
273  qp_back(0) + qp_back(1) > (1 + TOLERANCE))
274  mooseException("Quadrature point: ", qp_back, " out of bounds, truncating.");
275  }
276  else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
277  primal_elem->type() == QUAD9)
278  {
279  if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
280  qp_back(1) < (-1 - TOLERANCE) || qp_back(1) > (1 + TOLERANCE))
281  mooseException("Quadrature point: ", qp_back, " out of bounds, truncating");
282  }
283  }
284  else
285  {
286  mooseError("Newton iteration for mortar quadrature mapping msm element: ",
287  msm_elem->id(),
288  " to elem: ",
289  primal_elem->id(),
290  " didn't converge. MSM element volume: ",
291  msm_elem->volume());
292  }
293  }
294 }
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:19
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:323
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
TRI3
QUAD4
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
TRI6
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > RealEigenMatrix
Definition: MooseTypes.h:149
auto norm(const T &a) -> decltype(std::abs(a))
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
TRI7
QUAD9
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
Definition: MooseTypes.h:146

◆ 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 317 of file MortarUtils.h.

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

324 {
325  secondary_ip_sub_to_mats.clear();
326  primary_ip_sub_to_mats.clear();
327  secondary_boundary_mats.clear();
328 
329  auto & mat_warehouse = fe_problem.getRegularMaterialsWarehouse();
330  auto get_required_sub_mats =
331  [&mat_warehouse, tid, &consumers](
332  const SubdomainID sub_id,
333  const Moose::MaterialDataType mat_data_type) -> std::deque<MaterialBase *>
334  {
335  if (mat_warehouse[mat_data_type].hasActiveBlockObjects(sub_id, tid))
336  {
337  auto & sub_mats = mat_warehouse[mat_data_type].getActiveBlockObjects(sub_id, tid);
338  return MaterialBase::buildRequiredMaterials(consumers, sub_mats, /*allow_stateful=*/false);
339  }
340  else
341  return {};
342  };
343 
344  // Construct secondary *block* materials container
345  const auto & secondary_ip_sub_ids = amg.secondaryIPSubIDs();
346  for (const auto secondary_ip_sub : secondary_ip_sub_ids)
347  secondary_ip_sub_to_mats.emplace(
348  secondary_ip_sub, get_required_sub_mats(secondary_ip_sub, Moose::FACE_MATERIAL_DATA));
349 
350  // Construct primary *block* materials container
351  const auto & primary_ip_sub_ids = amg.primaryIPSubIDs();
352  for (const auto primary_ip_sub : primary_ip_sub_ids)
353  primary_ip_sub_to_mats.emplace(
354  primary_ip_sub, get_required_sub_mats(primary_ip_sub, Moose::NEIGHBOR_MATERIAL_DATA));
355 
356  // Construct secondary *boundary* materials container
357  const auto & boundary_pr = amg.primarySecondaryBoundaryIDPair();
358  const auto secondary_boundary = boundary_pr.second;
359  if (mat_warehouse.hasActiveBoundaryObjects(secondary_boundary, tid))
360  {
361  auto & boundary_mats = mat_warehouse.getActiveBoundaryObjects(secondary_boundary, tid);
362  secondary_boundary_mats =
363  MaterialBase::buildRequiredMaterials(consumers, boundary_mats, /*allow_stateful=*/false);
364  }
365 }
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:692
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:534