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