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

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

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

◆ 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:302
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 316 of file MortarUtils.h.

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

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