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 : * 3D projection operator for mapping qpoints on mortar segments to secondary or primary elements
28 : * @param msm_elem The mortar segment element that we will be mapping quadrature points from
29 : * @param primal_elem The "persistent" mesh element (e.g. it exists on the simulation's MooseMesh)
30 : * that we will be mapping quadrature points for. This can be either an element on the secondary or
31 : * primary face
32 : * @param sub_elem_index We will call \p msm_elem->get_extra_integer(sub_elem_index) in the
33 : * implementation in order to determine which sub-element of the primal element the mortar segment
34 : * element corresponds to. This \p sub_elem_index should correspond to the secondary element index
35 : * if \p primal_elem is a secondary face element and the primary element index if \p primal_elem is
36 : * a primary face element
37 : * @param qrule_msm The rule that governs quadrature on the mortar segment element
38 : * @param q_pts The output of this function. This will correspond to the the (reference space)
39 : * quadrature points that we wish to evaluate shape functions, etc., at on the primal element
40 : */
41 : void projectQPoints3d(const Elem * msm_elem,
42 : const Elem * primal_elem,
43 : unsigned int sub_elem_index,
44 : const QBase & qrule_msm,
45 : std::vector<Point> & q_pts);
46 :
47 : /**
48 : * This method will loop over pairs of secondary elements and their corresponding mortar segments,
49 : * reinitialize all finite element shape functions, variables, and material properties, and then
50 : * call a provided action function for each mortar segment
51 : * @param secondary_elems_to_mortar_segments This is a container of iterators. Each iterator should
52 : * point to a pair. The first member of the pair should be a pointer to a secondary face element and
53 : * the second member of the pair should correspond to a container of mortar segment element pointers
54 : * that correspond to the secondary face element.
55 : * @param assembly The object we will to use to reinitalize finite element data
56 : * @param subproblem The object we will use to reinitialize variables
57 : * @param fe_problem The object we will use to reinitialize material properties
58 : * @param amg The mortar mesh generation object which holds all the mortar mesh data
59 : * @param displaced Whether the mortar mesh was built from a displaced parent mesh
60 : * @param consumers A container of objects that are going to be using all the data that we are
61 : * reinitializing within this function. This may be, for instance, a container of mortar constraints
62 : * or auxiliary kernels. This \p consumers parameter is important as it allows us to build up
63 : * variable and material property dependencies that we must make sure we reinit
64 : * @param act The action functor that we will call for each mortar segment after we have
65 : * reinitalized all of our prereq data. This functor may, for instance, call \p computeResidual or
66 : * \p computeJacobian on mortar constraints, or \p computeValue for an auxiliary kernel
67 : */
68 : template <typename Iterators, typename Consumers, typename ActionFunctor>
69 : void
70 14518 : loopOverMortarSegments(
71 : const Iterators & secondary_elems_to_mortar_segments,
72 : Assembly & assembly,
73 : SubProblem & subproblem,
74 : FEProblemBase & fe_problem,
75 : const AutomaticMortarGeneration & amg,
76 : const bool displaced,
77 : const Consumers & consumers,
78 : const THREAD_ID tid,
79 : const std::map<SubdomainID, std::deque<MaterialBase *>> & secondary_ip_sub_to_mats,
80 : const std::map<SubdomainID, std::deque<MaterialBase *>> & primary_ip_sub_to_mats,
81 : const std::deque<MaterialBase *> & secondary_boundary_mats,
82 : const ActionFunctor act,
83 : const bool reinit_mortar_user_objects)
84 : {
85 14518 : const auto & primary_secondary_boundary_id_pair = amg.primarySecondaryBoundaryIDPair();
86 :
87 14518 : const auto primary_boundary_id = primary_secondary_boundary_id_pair.first;
88 14518 : const auto secondary_boundary_id = primary_secondary_boundary_id_pair.second;
89 :
90 : // For 3D mortar get index for retrieving sub-element info
91 14518 : unsigned int secondary_sub_elem_index = 0, primary_sub_elem_index = 0;
92 14518 : if (amg.dim() == 3)
93 : {
94 795 : secondary_sub_elem_index = amg.mortarSegmentMesh().get_elem_integer_index("secondary_sub_elem");
95 795 : 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 14518 : const auto & qrule_msm = assembly.qRuleMortar();
101 :
102 : // The element Jacobian times weights
103 14518 : const auto & JxW_msm = assembly.jxWMortar();
104 :
105 : // Set required material properties
106 14518 : std::unordered_set<unsigned int> needed_mat_props;
107 31916 : for (const auto & consumer : consumers)
108 : {
109 17398 : const auto & mp_deps = consumer->getMatPropDependencies();
110 17398 : needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
111 : }
112 14518 : fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/tid);
113 :
114 : // Loop through secondary elements, accumulating quadrature points for all corresponding mortar
115 : // segments
116 85331 : for (const auto elem_to_msm : secondary_elems_to_mortar_segments)
117 : {
118 70813 : const Elem * secondary_face_elem = subproblem.mesh().getMesh().elem_ptr(elem_to_msm->first);
119 : // Set the secondary interior parent and side ids
120 70813 : const Elem * secondary_ip = secondary_face_elem->interior_parent();
121 70813 : unsigned int secondary_side_id = secondary_ip->which_side_am_i(secondary_face_elem);
122 : const auto & secondary_ip_mats =
123 70813 : libmesh_map_find(secondary_ip_sub_to_mats, secondary_ip->subdomain_id());
124 :
125 70813 : 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 70813 : std::vector<Point> secondary_xi_pts, primary_xi_pts;
134 :
135 70813 : std::vector<Real> JxW;
136 :
137 : #ifndef NDEBUG
138 : unsigned int expected_length = 0;
139 : #endif
140 :
141 : // Loop through contributing msm elements
142 355242 : for (const auto msm_elem : msm_elems)
143 : {
144 : // Initialize mortar segment quadrature and compute JxW
145 284429 : subproblem.reinitMortarElem(msm_elem, tid);
146 :
147 : // Get a reference to the MortarSegmentInfo for this Elem.
148 284429 : const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
149 :
150 284429 : if (msm_elem->dim() == 1)
151 : {
152 236591 : for (unsigned int qp = 0; qp < qrule_msm->n_points(); qp++)
153 : {
154 158388 : const Real eta = qrule_msm->qp(qp)(0);
155 :
156 : // Map quadrature points to secondary side
157 158388 : const Real xi1_eta = 0.5 * (1 - eta) * msinfo.xi1_a + 0.5 * (1 + eta) * msinfo.xi1_b;
158 158388 : secondary_xi_pts.push_back(xi1_eta);
159 :
160 : // Map quadrature points to primary side
161 158388 : const Real xi2_eta = 0.5 * (1 - eta) * msinfo.xi2_a + 0.5 * (1 + eta) * msinfo.xi2_b;
162 158388 : 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 206226 : projectQPoints3d(msm_elem,
169 206226 : msinfo.secondary_elem,
170 : secondary_sub_elem_index,
171 : *qrule_msm,
172 : secondary_xi_pts);
173 206226 : projectQPoints3d(
174 206226 : 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 284429 : if (assembly.needDual())
179 27376 : 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 70813 : if (assembly.needDual())
197 4587 : assembly.reinitDual(secondary_face_elem, secondary_xi_pts, JxW);
198 :
199 70813 : unsigned int n_segment = 0;
200 :
201 : // Loop through contributing msm elements, computing residual and Jacobian this time
202 355242 : for (const auto msm_elem : msm_elems)
203 : {
204 284429 : n_segment++;
205 :
206 : // These will hold quadrature points for each segment
207 284429 : std::vector<Point> xi1_pts, xi2_pts;
208 :
209 : // Get a reference to the MortarSegmentInfo for this Elem.
210 284429 : const MortarSegmentInfo & msinfo = amg.mortarSegmentMeshElemToInfo().at(msm_elem);
211 :
212 : // Set the primary interior parent and side ids
213 284429 : const Elem * primary_ip = msinfo.primary_elem->interior_parent();
214 284429 : unsigned int primary_side_id = primary_ip->which_side_am_i(msinfo.primary_elem);
215 : const auto & primary_ip_mats =
216 284429 : 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 284429 : subproblem.reinitMortarElem(msm_elem, tid);
221 :
222 : // Extract previously computed mapped quadrature points for secondary and primary face
223 : // elements
224 284429 : const unsigned int start = (n_segment - 1) * qrule_msm->n_points();
225 284429 : const unsigned int end = n_segment * qrule_msm->n_points();
226 853287 : xi1_pts.insert(
227 853287 : xi1_pts.begin(), secondary_xi_pts.begin() + start, secondary_xi_pts.begin() + end);
228 284429 : xi2_pts.insert(xi2_pts.begin(), primary_xi_pts.begin() + start, primary_xi_pts.begin() + end);
229 :
230 284429 : 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 284429 : if (displaced)
235 22778 : 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 284429 : fe_problem.reinitElemFaceRef(
251 : reinit_secondary_elem, secondary_side_id, TOLERANCE, &xi1_pts, nullptr, tid);
252 :
253 284429 : 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 284429 : if (displaced)
258 22778 : reinit_primary_elem = fe_problem.mesh().elemPtr(reinit_primary_elem->id());
259 :
260 : // reinit the variables/residuals/jacobians on the primary interior
261 284429 : 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 284429 : 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 284429 : 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 284429 : 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 284429 : fe_problem.reinitMaterialsFace(secondary_ip->subdomain_id(),
286 : /*tid=*/tid,
287 : /*swap_stateful=*/false,
288 : &secondary_ip_mats);
289 284429 : fe_problem.reinitMaterialsBoundary(
290 : secondary_boundary_id, /*tid=*/tid, /*swap_stateful=*/false, &secondary_boundary_mats);
291 :
292 284429 : if (reinit_mortar_user_objects)
293 283073 : fe_problem.reinitMortarUserObjects(primary_boundary_id, secondary_boundary_id, displaced);
294 :
295 284429 : act();
296 :
297 : } // End loop over msm segments on secondary face elem
298 : } // End loop over (active) secondary elems
299 14518 : }
300 :
301 : /**
302 : * This function creates containers of materials necessary to execute the mortar method for a
303 : * supplied set of consumers
304 : * @param consumers The objects that we're building the material dependencies for. This could be a
305 : * container of mortar constraints or a "mortar" auxiliary kernel for example
306 : * @param fe_problem The finite element problem that we'll be querying for material warehouses
307 : * @param tid The thread ID that we will use for pulling the material warehouse
308 : * @param secondary_ip_sub_to_mats A map from the secondary interior parent subdomain IDs to the
309 : * required secondary \em block materials we will need to evaluate for the consumers
310 : * @param primary_ip_sub_to_mats A map from the primary interior parent subdomain IDs to the
311 : * required primary \em block materials we will need to evaluate for the consumers
312 : * @param secondary_boundary_mats The secondary \em boundary materials we will need to evaluate for
313 : * the consumers
314 : */
315 : template <typename Consumers>
316 : void
317 1281 : setupMortarMaterials(const Consumers & consumers,
318 : FEProblemBase & fe_problem,
319 : const AutomaticMortarGeneration & amg,
320 : const THREAD_ID tid,
321 : std::map<SubdomainID, std::deque<MaterialBase *>> & secondary_ip_sub_to_mats,
322 : std::map<SubdomainID, std::deque<MaterialBase *>> & primary_ip_sub_to_mats,
323 : std::deque<MaterialBase *> & secondary_boundary_mats)
324 : {
325 1281 : secondary_ip_sub_to_mats.clear();
326 1281 : primary_ip_sub_to_mats.clear();
327 1281 : secondary_boundary_mats.clear();
328 :
329 1281 : auto & mat_warehouse = fe_problem.getRegularMaterialsWarehouse();
330 1281 : auto get_required_sub_mats =
331 2714 : [&mat_warehouse, tid, &consumers](
332 : const SubdomainID sub_id,
333 : const Moose::MaterialDataType mat_data_type) -> std::deque<MaterialBase *>
334 : {
335 2714 : if (mat_warehouse[mat_data_type].hasActiveBlockObjects(sub_id, tid))
336 : {
337 786 : auto & sub_mats = mat_warehouse[mat_data_type].getActiveBlockObjects(sub_id, tid);
338 786 : return MaterialBase::buildRequiredMaterials(consumers, sub_mats, /*allow_stateful=*/false);
339 : }
340 : else
341 1928 : return {};
342 : };
343 :
344 : // Construct secondary *block* materials container
345 1281 : const auto & secondary_ip_sub_ids = amg.secondaryIPSubIDs();
346 2638 : for (const auto secondary_ip_sub : secondary_ip_sub_ids)
347 1357 : 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 1281 : const auto & primary_ip_sub_ids = amg.primaryIPSubIDs();
352 2638 : for (const auto primary_ip_sub : primary_ip_sub_ids)
353 1357 : 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 1281 : const auto & boundary_pr = amg.primarySecondaryBoundaryIDPair();
358 1281 : const auto secondary_boundary = boundary_pr.second;
359 1281 : if (mat_warehouse.hasActiveBoundaryObjects(secondary_boundary, tid))
360 : {
361 13 : auto & boundary_mats = mat_warehouse.getActiveBoundaryObjects(secondary_boundary, tid);
362 13 : secondary_boundary_mats =
363 : MaterialBase::buildRequiredMaterials(consumers, boundary_mats, /*allow_stateful=*/false);
364 : }
365 1281 : }
366 : }
367 : }
|
