Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "Assembly.h"
13 : #include "FEProblemBase.h"
14 : #include "MaterialBase.h"
15 : #include "MaterialWarehouse.h"
16 : #include "AutomaticMortarGeneration.h"
17 :
18 : #include "libmesh/quadrature.h"
19 : #include "libmesh/elem.h"
20 : #include "libmesh/point.h"
21 :
22 : namespace Moose
23 : {
24 : namespace Mortar
25 : {
26 : /**
27 : * Return the node indices for a first-order sub-element of a parent face.
28 : */
29 : std::vector<unsigned int> getMortarSubElementNodeIndices(const Elem & parent_elem,
30 : unsigned int sub_elem);
31 :
32 : /**
33 : * 3D projection operator for mapping qpoints on mortar segments to secondary or primary elements
34 : * @param msm_elem The mortar segment element that we will be mapping quadrature points from
35 : * @param primal_elem The "persistent" mesh element (e.g. it exists on the simulation's MooseMesh)
36 : * that we will be mapping quadrature points for. This can be either an element on the secondary or
37 : * primary face
38 : * @param sub_elem_index We will call \p msm_elem->get_extra_integer(sub_elem_index) in the
39 : * implementation in order to determine which sub-element of the primal element the mortar segment
40 : * element corresponds to. This \p sub_elem_index should correspond to the secondary element index
41 : * if \p primal_elem is a secondary face element and the primary element index if \p primal_elem is
42 : * a primary face element
43 : * @param qrule_msm The rule that governs quadrature on the mortar segment element
44 : * @param q_pts The output of this function. This will correspond to the the (reference space)
45 : * quadrature points that we wish to evaluate shape functions, etc., at on the primal element
46 : */
47 : void projectQPoints3d(const Elem * msm_elem,
48 : const Elem * primal_elem,
49 : unsigned int sub_elem_index,
50 : const QBase & qrule_msm,
51 : std::vector<Point> & q_pts);
52 :
53 : /**
54 : * 3D mapping operator that interpolates stored parent reference points on each triangular mortar
55 : * segment.
56 : * @param mortar_segment_elem The triangular mortar segment carrying the quadrature rule
57 : * @param reference_points Parent-face reference points stored for the mortar segment
58 : * @param qrule_msm The rule that governs quadrature on the mortar segment element
59 : * @param secondary_q_pts Reference-space quadrature points on the secondary face
60 : * @param primary_q_pts Reference-space quadrature points on the primary face
61 : */
62 : void mapQPoints3dFromReference(const Elem & mortar_segment_elem,
63 : const MortarSegmentReferencePoints & reference_points,
64 : const QBase & qrule_msm,
65 : std::vector<Point> & secondary_q_pts,
66 : std::vector<Point> & primary_q_pts);
67 :
68 : /**
69 : * This method will loop over pairs of secondary elements and their corresponding mortar segments,
70 : * reinitialize all finite element shape functions, variables, and material properties, and then
71 : * call a provided action function for each mortar segment
72 : * @param secondary_elems_to_mortar_segments This is a container of iterators. Each iterator should
73 : * point to a pair. The first member of the pair should be a pointer to a secondary face element and
74 : * the second member of the pair should correspond to a container of mortar segment element pointers
75 : * that correspond to the secondary face element.
76 : * @param assembly The object we will to use to reinitalize finite element data
77 : * @param subproblem The object we will use to reinitialize variables
78 : * @param fe_problem The object we will use to reinitialize material properties
79 : * @param amg The mortar mesh generation object which holds all the mortar mesh data
80 : * @param displaced Whether the mortar mesh was built from a displaced parent mesh
81 : * @param consumers A container of objects that are going to be using all the data that we are
82 : * reinitializing within this function. This may be, for instance, a container of mortar constraints
83 : * or auxiliary kernels. This \p consumers parameter is important as it allows us to build up
84 : * variable and material property dependencies that we must make sure we reinit
85 : * @param act The action functor that we will call for each mortar segment after we have
86 : * reinitalized all of our prereq data. This functor may, for instance, call \p computeResidual or
87 : * \p computeJacobian on mortar constraints, or \p computeValue for an auxiliary kernel
88 : */
89 : template <typename Iterators, typename Consumers, typename ActionFunctor>
90 : void
91 13339 : loopOverMortarSegments(
92 : const Iterators & secondary_elems_to_mortar_segments,
93 : Assembly & assembly,
94 : SubProblem & subproblem,
95 : FEProblemBase & fe_problem,
96 : const AutomaticMortarGeneration & amg,
97 : const bool displaced,
98 : const Consumers & consumers,
99 : const THREAD_ID tid,
100 : const std::map<SubdomainID, std::deque<MaterialBase *>> & secondary_ip_sub_to_mats,
101 : const std::map<SubdomainID, std::deque<MaterialBase *>> & primary_ip_sub_to_mats,
102 : const std::deque<MaterialBase *> & secondary_boundary_mats,
103 : const ActionFunctor act,
104 : const bool reinit_mortar_user_objects)
105 : {
106 13339 : const auto & primary_secondary_boundary_id_pair = amg.primarySecondaryBoundaryIDPair();
107 :
108 13339 : const auto primary_boundary_id = primary_secondary_boundary_id_pair.first;
109 13339 : const auto secondary_boundary_id = primary_secondary_boundary_id_pair.second;
110 :
111 : // For 3D mortar get index for retrieving sub-element info
112 13339 : unsigned int secondary_sub_elem_index = 0, primary_sub_elem_index = 0;
113 13339 : if (amg.dim() == 3)
114 : {
115 1097 : secondary_sub_elem_index = amg.mortarSegmentMesh().get_elem_integer_index("secondary_sub_elem");
116 1097 : 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 13339 : const auto & qrule_msm = assembly.qRuleMortar();
122 :
123 : // The element Jacobian times weights
124 13339 : const auto & JxW_msm = assembly.jxWMortar();
125 :
126 : // Set required material properties
127 13339 : std::unordered_set<unsigned int> needed_mat_props;
128 29370 : for (const auto & consumer : consumers)
129 : {
130 16031 : const auto & mp_deps = consumer->getMatPropDependencies();
131 16031 : needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
132 : }
133 13339 : fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/tid);
134 :
135 : // Loop through secondary elements, accumulating quadrature points for all corresponding mortar
136 : // segments
137 80860 : for (const auto elem_to_msm : secondary_elems_to_mortar_segments)
138 : {
139 67521 : const Elem * secondary_face_elem = subproblem.mesh().getMesh().elem_ptr(elem_to_msm->first);
140 : // Set the secondary interior parent and side ids
141 67521 : const Elem * secondary_ip = secondary_face_elem->interior_parent();
142 67521 : unsigned int secondary_side_id = secondary_ip->which_side_am_i(secondary_face_elem);
143 : const auto & secondary_ip_mats =
144 67521 : libmesh_map_find(secondary_ip_sub_to_mats, secondary_ip->subdomain_id());
145 :
146 67521 : 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 67521 : std::vector<Point> secondary_xi_pts, primary_xi_pts;
155 :
156 67521 : std::vector<Real> JxW;
157 :
158 : #ifndef NDEBUG
159 : unsigned int expected_length = 0;
160 : #endif
161 :
162 : // Loop through contributing msm elements
163 631525 : for (const auto msm_elem : msm_elems)
164 : {
165 : // Initialize mortar segment quadrature and compute JxW
166 564004 : subproblem.reinitMortarElem(msm_elem, tid);
167 :
168 : // Get a reference to the MortarSegmentInfo for this Elem.
169 564004 : const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
170 :
171 564004 : if (msm_elem->dim() == 1)
172 : {
173 198004 : for (unsigned int qp = 0; qp < qrule_msm->n_points(); qp++)
174 : {
175 133816 : const Real eta = qrule_msm->qp(qp)(0);
176 :
177 : // Map quadrature points to secondary side
178 133816 : const Real xi1_eta = 0.5 * (1 - eta) * msinfo.xi1_a + 0.5 * (1 + eta) * msinfo.xi1_b;
179 133816 : secondary_xi_pts.push_back(xi1_eta);
180 :
181 : // Map quadrature points to primary side
182 133816 : const Real xi2_eta = 0.5 * (1 - eta) * msinfo.xi2_a + 0.5 * (1 + eta) * msinfo.xi2_b;
183 133816 : primary_xi_pts.push_back(xi2_eta);
184 : }
185 : }
186 : else
187 : {
188 499816 : if (amg.mortar3DQpMapping() == Mortar3DQuadraturePointMapping::REFERENCE_INTERPOLATION)
189 44034 : 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 455782 : projectQPoints3d(msm_elem,
198 455782 : msinfo.secondary_elem,
199 : secondary_sub_elem_index,
200 : *qrule_msm,
201 : secondary_xi_pts);
202 455782 : projectQPoints3d(
203 455782 : 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 564004 : if (assembly.needDual())
209 24284 : 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 67521 : if (assembly.needDual())
227 4057 : assembly.reinitDual(secondary_face_elem, secondary_xi_pts, JxW);
228 :
229 67521 : unsigned int n_segment = 0;
230 :
231 : // Loop through contributing msm elements, computing residual and Jacobian this time
232 631525 : for (const auto msm_elem : msm_elems)
233 : {
234 564004 : n_segment++;
235 :
236 : // These will hold quadrature points for each segment
237 564004 : std::vector<Point> xi1_pts, xi2_pts;
238 :
239 : // Get a reference to the MortarSegmentInfo for this Elem.
240 564004 : const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
241 :
242 : // Set the primary interior parent and side ids
243 564004 : const Elem * primary_ip = msinfo.primary_elem->interior_parent();
244 564004 : unsigned int primary_side_id = primary_ip->which_side_am_i(msinfo.primary_elem);
245 : const auto & primary_ip_mats =
246 564004 : 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 564004 : subproblem.reinitMortarElem(msm_elem, tid);
251 :
252 : // Extract previously computed mapped quadrature points for secondary and primary face
253 : // elements
254 564004 : const unsigned int start = (n_segment - 1) * qrule_msm->n_points();
255 564004 : const unsigned int end = n_segment * qrule_msm->n_points();
256 1692012 : xi1_pts.insert(
257 1692012 : xi1_pts.begin(), secondary_xi_pts.begin() + start, secondary_xi_pts.begin() + end);
258 564004 : xi2_pts.insert(xi2_pts.begin(), primary_xi_pts.begin() + start, primary_xi_pts.begin() + end);
259 :
260 564004 : 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 564004 : if (displaced)
265 20096 : 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 564004 : fe_problem.reinitElemFaceRef(
281 : reinit_secondary_elem, secondary_side_id, TOLERANCE, &xi1_pts, nullptr, tid);
282 :
283 564004 : 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 564004 : if (displaced)
288 20096 : reinit_primary_elem = fe_problem.mesh().elemPtr(reinit_primary_elem->id());
289 :
290 : // reinit the variables/residuals/jacobians on the primary interior
291 564004 : 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 564004 : 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 564004 : 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 564004 : 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 564004 : fe_problem.reinitMaterialsFace(secondary_ip->subdomain_id(),
316 : /*tid=*/tid,
317 : /*swap_stateful=*/false,
318 : &secondary_ip_mats);
319 564004 : fe_problem.reinitMaterialsBoundary(
320 : secondary_boundary_id, /*tid=*/tid, /*swap_stateful=*/false, &secondary_boundary_mats);
321 :
322 564004 : if (reinit_mortar_user_objects)
323 562808 : fe_problem.reinitMortarUserObjects(primary_boundary_id, secondary_boundary_id, displaced);
324 :
325 564004 : act();
326 :
327 : } // End loop over msm segments on secondary face elem
328 : } // End loop over (active) secondary elems
329 13339 : }
330 :
331 : /**
332 : * This function creates containers of materials necessary to execute the mortar method for a
333 : * supplied set of consumers
334 : * @param consumers The objects that we're building the material dependencies for. This could be a
335 : * container of mortar constraints or a "mortar" auxiliary kernel for example
336 : * @param fe_problem The finite element problem that we'll be querying for material warehouses
337 : * @param tid The thread ID that we will use for pulling the material warehouse
338 : * @param secondary_ip_sub_to_mats A map from the secondary interior parent subdomain IDs to the
339 : * required secondary \em block materials we will need to evaluate for the consumers
340 : * @param primary_ip_sub_to_mats A map from the primary interior parent subdomain IDs to the
341 : * required primary \em block materials we will need to evaluate for the consumers
342 : * @param secondary_boundary_mats The secondary \em boundary materials we will need to evaluate for
343 : * the consumers
344 : */
345 : template <typename Consumers>
346 : void
347 1300 : setupMortarMaterials(const Consumers & consumers,
348 : FEProblemBase & fe_problem,
349 : const AutomaticMortarGeneration & amg,
350 : const THREAD_ID tid,
351 : std::map<SubdomainID, std::deque<MaterialBase *>> & secondary_ip_sub_to_mats,
352 : std::map<SubdomainID, std::deque<MaterialBase *>> & primary_ip_sub_to_mats,
353 : std::deque<MaterialBase *> & secondary_boundary_mats)
354 : {
355 1300 : secondary_ip_sub_to_mats.clear();
356 1300 : primary_ip_sub_to_mats.clear();
357 1300 : secondary_boundary_mats.clear();
358 :
359 1300 : auto & mat_warehouse = fe_problem.getRegularMaterialsWarehouse();
360 1300 : auto get_required_sub_mats =
361 2648 : [&mat_warehouse, tid, &consumers](
362 : const SubdomainID sub_id,
363 : const Moose::MaterialDataType mat_data_type) -> std::deque<MaterialBase *>
364 : {
365 2648 : if (mat_warehouse[mat_data_type].hasActiveBlockObjects(sub_id, tid))
366 : {
367 736 : auto & sub_mats = mat_warehouse[mat_data_type].getActiveBlockObjects(sub_id, tid);
368 736 : return MaterialBase::buildRequiredMaterials(consumers, sub_mats, /*allow_stateful=*/false);
369 : }
370 : else
371 1912 : return {};
372 : };
373 :
374 : // Construct secondary *block* materials container
375 1300 : const auto & secondary_ip_sub_ids = amg.secondaryIPSubIDs();
376 2624 : for (const auto secondary_ip_sub : secondary_ip_sub_ids)
377 1324 : 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 1300 : const auto & primary_ip_sub_ids = amg.primaryIPSubIDs();
382 2624 : for (const auto primary_ip_sub : primary_ip_sub_ids)
383 1324 : 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 1300 : const auto & boundary_pr = amg.primarySecondaryBoundaryIDPair();
388 1300 : const auto secondary_boundary = boundary_pr.second;
389 1300 : if (mat_warehouse.hasActiveBoundaryObjects(secondary_boundary, tid))
390 : {
391 12 : auto & boundary_mats = mat_warehouse.getActiveBoundaryObjects(secondary_boundary, tid);
392 12 : secondary_boundary_mats =
393 : MaterialBase::buildRequiredMaterials(consumers, boundary_mats, /*allow_stateful=*/false);
394 : }
395 1300 : }
396 : }
397 : }
|
