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 : #include "Assembly.h"
11 :
12 : // MOOSE includes
13 : #include "SubProblem.h"
14 : #include "ArbitraryQuadrature.h"
15 : #include "SystemBase.h"
16 : #include "MooseTypes.h"
17 : #include "MooseMesh.h"
18 : #include "MooseVariableFE.h"
19 : #include "MooseVariableScalar.h"
20 : #include "XFEMInterface.h"
21 : #include "DisplacedSystem.h"
22 : #include "MooseMeshUtils.h"
23 :
24 : // libMesh
25 : #include "libmesh/coupling_matrix.h"
26 : #include "libmesh/dof_map.h"
27 : #include "libmesh/elem.h"
28 : #include "libmesh/equation_systems.h"
29 : #include "libmesh/fe_interface.h"
30 : #include "libmesh/node.h"
31 : #include "libmesh/quadrature_gauss.h"
32 : #include "libmesh/sparse_matrix.h"
33 : #include "libmesh/tensor_value.h"
34 : #include "libmesh/vector_value.h"
35 : #include "libmesh/fe.h"
36 : #include "libmesh/static_condensation.h"
37 :
38 : using namespace libMesh;
39 :
40 : template <typename P, typename C>
41 : void
42 1823172886 : coordTransformFactor(const SubProblem & s,
43 : const SubdomainID sub_id,
44 : const P & point,
45 : C & factor,
46 : const SubdomainID neighbor_sub_id)
47 : {
48 1823172886 : coordTransformFactor(s.mesh(), sub_id, point, factor, neighbor_sub_id);
49 1823172886 : }
50 :
51 : template <typename P, typename C>
52 : void
53 1826490689 : coordTransformFactor(const MooseMesh & mesh,
54 : const SubdomainID sub_id,
55 : const P & point,
56 : C & factor,
57 : const SubdomainID libmesh_dbg_var(neighbor_sub_id))
58 : {
59 : mooseAssert(neighbor_sub_id != libMesh::Elem::invalid_subdomain_id
60 : ? mesh.getCoordSystem(sub_id) == mesh.getCoordSystem(neighbor_sub_id)
61 : : true,
62 : "Coordinate systems must be the same between element and neighbor");
63 1826490689 : const auto coord_type = mesh.getCoordSystem(sub_id);
64 :
65 1826490689 : if (coord_type == Moose::COORD_RZ)
66 : {
67 3858246 : if (mesh.usingGeneralAxisymmetricCoordAxes())
68 : {
69 1065929 : const auto & axis = mesh.getGeneralAxisymmetricCoordAxis(sub_id);
70 1065929 : MooseMeshUtils::coordTransformFactorRZGeneral(point, axis, factor);
71 : }
72 : else
73 2792317 : MooseMeshUtils::coordTransformFactor(
74 : point, factor, coord_type, mesh.getAxisymmetricRadialCoord());
75 : }
76 : else
77 1822632443 : MooseMeshUtils::coordTransformFactor(point, factor, coord_type, libMesh::invalid_uint);
78 1826490689 : }
79 :
80 66571 : Assembly::Assembly(SystemBase & sys, THREAD_ID tid)
81 66571 : : _sys(sys),
82 133142 : _subproblem(_sys.subproblem()),
83 66571 : _displaced(dynamic_cast<DisplacedSystem *>(&sys) ? true : false),
84 66571 : _nonlocal_cm(_subproblem.nonlocalCouplingMatrix(_sys.number())),
85 66571 : _computing_residual(_subproblem.currentlyComputingResidual()),
86 66571 : _computing_jacobian(_subproblem.currentlyComputingJacobian()),
87 66571 : _computing_residual_and_jacobian(_subproblem.currentlyComputingResidualAndJacobian()),
88 66571 : _dof_map(_sys.dofMap()),
89 66571 : _tid(tid),
90 66571 : _mesh(sys.mesh()),
91 66571 : _mesh_dimension(_mesh.dimension()),
92 66571 : _helper_type(_mesh.hasSecondOrderElements() ? SECOND : FIRST, LAGRANGE),
93 66571 : _user_added_fe_of_helper_type(false),
94 66571 : _user_added_fe_face_of_helper_type(false),
95 66571 : _user_added_fe_face_neighbor_of_helper_type(false),
96 66571 : _user_added_fe_neighbor_of_helper_type(false),
97 66571 : _user_added_fe_lower_of_helper_type(false),
98 66571 : _building_helpers(false),
99 66571 : _current_qrule(nullptr),
100 66571 : _current_qrule_volume(nullptr),
101 66571 : _current_qrule_arbitrary(nullptr),
102 66571 : _coord_type(Moose::COORD_XYZ),
103 66571 : _current_qrule_face(nullptr),
104 66571 : _current_qface_arbitrary(nullptr),
105 66571 : _current_qrule_neighbor(nullptr),
106 66571 : _need_JxW_neighbor(false),
107 66571 : _qrule_msm(nullptr),
108 66571 : _custom_mortar_qrule(false),
109 66571 : _current_qrule_lower(nullptr),
110 :
111 66571 : _current_elem(nullptr),
112 66571 : _current_elem_volume(0),
113 66571 : _current_side(0),
114 66571 : _current_side_elem(nullptr),
115 66571 : _current_side_volume(0),
116 66571 : _current_neighbor_elem(nullptr),
117 66571 : _current_neighbor_side(0),
118 66571 : _current_neighbor_side_elem(nullptr),
119 66571 : _need_neighbor_elem_volume(false),
120 66571 : _current_neighbor_volume(0),
121 66571 : _current_node(nullptr),
122 66571 : _current_neighbor_node(nullptr),
123 66571 : _current_elem_volume_computed(false),
124 66571 : _current_side_volume_computed(false),
125 :
126 66571 : _current_lower_d_elem(nullptr),
127 66571 : _current_neighbor_lower_d_elem(nullptr),
128 66571 : _need_lower_d_elem_volume(false),
129 66571 : _need_neighbor_lower_d_elem_volume(false),
130 66571 : _need_dual(false),
131 :
132 66571 : _residual_vector_tags(_subproblem.getVectorTags(Moose::VECTOR_TAG_RESIDUAL)),
133 66571 : _cached_residual_values(2), // The 2 is for TIME and NONTIME
134 66571 : _cached_residual_rows(2), // The 2 is for TIME and NONTIME
135 66571 : _max_cached_residuals(0),
136 66571 : _max_cached_jacobians(0),
137 :
138 66571 : _block_diagonal_matrix(false),
139 66571 : _calculate_xyz(false),
140 66571 : _calculate_face_xyz(false),
141 66571 : _calculate_curvatures(false),
142 66571 : _calculate_ad_coord(false),
143 66571 : _have_p_refinement(false),
144 266284 : _sc(nullptr)
145 : {
146 66571 : const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
147 66571 : _building_helpers = true;
148 : // Build fe's for the helpers
149 66571 : buildFE(FEType(helper_order, LAGRANGE));
150 66571 : buildFaceFE(FEType(helper_order, LAGRANGE));
151 66571 : buildNeighborFE(FEType(helper_order, LAGRANGE));
152 66571 : buildFaceNeighborFE(FEType(helper_order, LAGRANGE));
153 66571 : buildLowerDFE(FEType(helper_order, LAGRANGE));
154 66571 : _building_helpers = false;
155 :
156 : // Build an FE helper object for this type for each dimension up to the dimension of the current
157 : // mesh
158 258653 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
159 : {
160 192082 : _holder_fe_helper[dim] = _fe[dim][FEType(helper_order, LAGRANGE)];
161 192082 : _holder_fe_face_helper[dim] = _fe_face[dim][FEType(helper_order, LAGRANGE)];
162 192082 : _holder_fe_face_neighbor_helper[dim] = _fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)];
163 192082 : _holder_fe_neighbor_helper[dim] = _fe_neighbor[dim][FEType(helper_order, LAGRANGE)];
164 : }
165 :
166 192082 : for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
167 125511 : _holder_fe_lower_helper[dim] = _fe_lower[dim][FEType(helper_order, LAGRANGE)];
168 :
169 : // request phi, dphi, xyz, JxW, etc. data
170 66571 : helpersRequestData();
171 :
172 : // For 3D mortar, mortar segments are always TRI3 elements so we want FIRST LAGRANGE regardless
173 : // of discretization
174 66571 : _fe_msm = (_mesh_dimension == 2)
175 133142 : ? FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(helper_order, LAGRANGE))
176 66571 : : FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(FIRST, LAGRANGE));
177 : // This FE object should not take part in p-refinement
178 66571 : _fe_msm->add_p_level_in_reinit(false);
179 66571 : _JxW_msm = &_fe_msm->get_JxW();
180 : // Prerequest xyz so that it is computed for _fe_msm so that it can be used for calculating
181 : // _coord_msm
182 66571 : _fe_msm->get_xyz();
183 :
184 66571 : _extra_elem_ids.resize(_mesh.getMesh().n_elem_integers() + 1);
185 66571 : _neighbor_extra_elem_ids.resize(_mesh.getMesh().n_elem_integers() + 1);
186 66571 : }
187 :
188 122386 : Assembly::~Assembly()
189 : {
190 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
191 435440 : for (auto & it : _fe[dim])
192 258259 : delete it.second;
193 :
194 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
195 435440 : for (auto & it : _fe_face[dim])
196 258259 : delete it.second;
197 :
198 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
199 435440 : for (auto & it : _fe_neighbor[dim])
200 258259 : delete it.second;
201 :
202 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
203 435440 : for (auto & it : _fe_face_neighbor[dim])
204 258259 : delete it.second;
205 :
206 177181 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
207 278596 : for (auto & it : _fe_lower[dim])
208 162608 : delete it.second;
209 :
210 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
211 180190 : for (auto & it : _vector_fe[dim])
212 3009 : delete it.second;
213 :
214 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
215 180190 : for (auto & it : _vector_fe_face[dim])
216 3009 : delete it.second;
217 :
218 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
219 180190 : for (auto & it : _vector_fe_neighbor[dim])
220 3009 : delete it.second;
221 :
222 238374 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
223 180190 : for (auto & it : _vector_fe_face_neighbor[dim])
224 3009 : delete it.second;
225 :
226 177181 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
227 117536 : for (auto & it : _vector_fe_lower[dim])
228 1548 : delete it.second;
229 :
230 137803 : for (auto & it : _ad_grad_phi_data)
231 76610 : it.second.release();
232 :
233 62402 : for (auto & it : _ad_vector_grad_phi_data)
234 1209 : it.second.release();
235 :
236 131809 : for (auto & it : _ad_grad_phi_data_face)
237 70616 : it.second.release();
238 :
239 62402 : for (auto & it : _ad_vector_grad_phi_data_face)
240 1209 : it.second.release();
241 :
242 61193 : _current_physical_points.release();
243 :
244 61193 : _coord.release();
245 61193 : _coord_neighbor.release();
246 61193 : _coord_msm.release();
247 :
248 61193 : _ad_JxW.release();
249 61193 : _ad_q_points.release();
250 61193 : _ad_JxW_face.release();
251 61193 : _ad_normals.release();
252 61193 : _ad_q_points_face.release();
253 61193 : _curvatures.release();
254 61193 : _ad_curvatures.release();
255 61193 : _ad_coord.release();
256 :
257 61193 : delete _qrule_msm;
258 122386 : }
259 :
260 : const MooseArray<Real> &
261 93 : Assembly::JxWNeighbor() const
262 : {
263 93 : _need_JxW_neighbor = true;
264 93 : return _current_JxW_neighbor;
265 : }
266 :
267 : void
268 440628 : Assembly::buildFE(FEType type) const
269 : {
270 440628 : if (!_building_helpers && type == _helper_type)
271 195257 : _user_added_fe_of_helper_type = true;
272 :
273 440628 : if (!_fe_shape_data[type])
274 95647 : _fe_shape_data[type] = std::make_unique<FEShapeData>();
275 :
276 : // Build an FE object for this type for each dimension up to the dimension of the current mesh
277 1759972 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
278 : {
279 1319344 : if (!_fe[dim][type])
280 278159 : _fe[dim][type] = FEGenericBase<Real>::build(dim, type).release();
281 :
282 1319344 : _fe[dim][type]->get_phi();
283 1319344 : _fe[dim][type]->get_dphi();
284 : // Pre-request xyz. We have always computed xyz, but due to
285 : // recent optimizations in libmesh, we now need to explicity
286 : // request it, since apps (Yak) may rely on it being computed.
287 1319344 : _fe[dim][type]->get_xyz();
288 1319344 : if (_need_second_derivative.count(type))
289 80632 : _fe[dim][type]->get_d2phi();
290 : }
291 440628 : }
292 :
293 : void
294 414552 : Assembly::buildFaceFE(FEType type) const
295 : {
296 414552 : if (!_building_helpers && type == _helper_type)
297 190110 : _user_added_fe_face_of_helper_type = true;
298 :
299 414552 : if (!_fe_shape_data_face[type])
300 95647 : _fe_shape_data_face[type] = std::make_unique<FEShapeData>();
301 :
302 : // Build an FE object for this type for each dimension up to the dimension of the current mesh
303 1655608 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
304 : {
305 1241056 : if (!_fe_face[dim][type])
306 278159 : _fe_face[dim][type] = FEGenericBase<Real>::build(dim, type).release();
307 :
308 1241056 : _fe_face[dim][type]->get_phi();
309 1241056 : _fe_face[dim][type]->get_dphi();
310 1241056 : if (_need_second_derivative.count(type))
311 16198 : _fe_face[dim][type]->get_d2phi();
312 : }
313 414552 : }
314 :
315 : void
316 409196 : Assembly::buildNeighborFE(FEType type) const
317 : {
318 409196 : if (!_building_helpers && type == _helper_type)
319 189998 : _user_added_fe_neighbor_of_helper_type = true;
320 :
321 409196 : if (!_fe_shape_data_neighbor[type])
322 95647 : _fe_shape_data_neighbor[type] = std::make_unique<FEShapeData>();
323 :
324 : // Build an FE object for this type for each dimension up to the dimension of the current mesh
325 1634171 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
326 : {
327 1224975 : if (!_fe_neighbor[dim][type])
328 278159 : _fe_neighbor[dim][type] = FEGenericBase<Real>::build(dim, type).release();
329 :
330 1224975 : _fe_neighbor[dim][type]->get_phi();
331 1224975 : _fe_neighbor[dim][type]->get_dphi();
332 1224975 : if (_need_second_derivative_neighbor.count(type))
333 117 : _fe_neighbor[dim][type]->get_d2phi();
334 : }
335 409196 : }
336 :
337 : void
338 409196 : Assembly::buildFaceNeighborFE(FEType type) const
339 : {
340 409196 : if (!_building_helpers && type == _helper_type)
341 189998 : _user_added_fe_face_neighbor_of_helper_type = true;
342 :
343 409196 : if (!_fe_shape_data_face_neighbor[type])
344 95647 : _fe_shape_data_face_neighbor[type] = std::make_unique<FEShapeData>();
345 :
346 : // Build an FE object for this type for each dimension up to the dimension of the current mesh
347 1634171 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
348 : {
349 1224975 : if (!_fe_face_neighbor[dim][type])
350 278159 : _fe_face_neighbor[dim][type] = FEGenericBase<Real>::build(dim, type).release();
351 :
352 1224975 : _fe_face_neighbor[dim][type]->get_phi();
353 1224975 : _fe_face_neighbor[dim][type]->get_dphi();
354 1224975 : if (_need_second_derivative_neighbor.count(type))
355 117 : _fe_face_neighbor[dim][type]->get_d2phi();
356 : }
357 409196 : }
358 :
359 : void
360 720463 : Assembly::buildLowerDFE(FEType type) const
361 : {
362 720463 : if (!_building_helpers && type == _helper_type)
363 379456 : _user_added_fe_lower_of_helper_type = true;
364 :
365 720463 : if (!_fe_shape_data_lower[type])
366 90821 : _fe_shape_data_lower[type] = std::make_unique<FEShapeData>();
367 :
368 : // Build an FE object for this type for each dimension up to the dimension of
369 : // the current mesh minus one (because this is for lower-dimensional
370 : // elements!)
371 2178546 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
372 : {
373 1458083 : if (!_fe_lower[dim][type])
374 175206 : _fe_lower[dim][type] = FEGenericBase<Real>::build(dim, type).release();
375 :
376 1458083 : _fe_lower[dim][type]->get_phi();
377 1458083 : _fe_lower[dim][type]->get_dphi();
378 1458083 : if (_need_second_derivative.count(type))
379 0 : _fe_lower[dim][type]->get_d2phi();
380 : }
381 720463 : }
382 :
383 : void
384 540 : Assembly::buildLowerDDualFE(FEType type) const
385 : {
386 540 : if (!_fe_shape_data_dual_lower[type])
387 135 : _fe_shape_data_dual_lower[type] = std::make_unique<FEShapeData>();
388 :
389 : // Build an FE object for this type for each dimension up to the dimension of
390 : // the current mesh minus one (because this is for lower-dimensional
391 : // elements!)
392 1672 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
393 : {
394 1132 : if (!_fe_lower[dim][type])
395 0 : _fe_lower[dim][type] = FEGenericBase<Real>::build(dim, type).release();
396 :
397 1132 : _fe_lower[dim][type]->get_dual_phi();
398 1132 : _fe_lower[dim][type]->get_dual_dphi();
399 1132 : if (_need_second_derivative.count(type))
400 0 : _fe_lower[dim][type]->get_dual_d2phi();
401 : }
402 540 : }
403 :
404 : void
405 6284 : Assembly::buildVectorLowerDFE(FEType type) const
406 : {
407 6284 : if (!_vector_fe_shape_data_lower[type])
408 1259 : _vector_fe_shape_data_lower[type] = std::make_unique<VectorFEShapeData>();
409 :
410 : // Build an FE object for this type for each dimension up to the dimension of
411 : // the current mesh minus one (because this is for lower-dimensional
412 : // elements!)
413 6284 : unsigned int dim = ((type.family == LAGRANGE_VEC) || (type.family == MONOMIAL_VEC)) ? 0 : 2;
414 6284 : const auto ending_dim = cast_int<unsigned int>(_mesh_dimension - 1);
415 6284 : if (ending_dim < dim)
416 1716 : return;
417 13536 : for (; dim <= ending_dim; dim++)
418 : {
419 8968 : if (!_vector_fe_lower[dim][type])
420 1663 : _vector_fe_lower[dim][type] = FEVectorBase::build(dim, type).release();
421 :
422 8968 : _vector_fe_lower[dim][type]->get_phi();
423 8968 : _vector_fe_lower[dim][type]->get_dphi();
424 8968 : if (_need_second_derivative.count(type))
425 0 : _vector_fe_lower[dim][type]->get_d2phi();
426 : }
427 : }
428 :
429 : void
430 0 : Assembly::buildVectorDualLowerDFE(FEType type) const
431 : {
432 0 : if (!_vector_fe_shape_data_dual_lower[type])
433 0 : _vector_fe_shape_data_dual_lower[type] = std::make_unique<VectorFEShapeData>();
434 :
435 : // Build an FE object for this type for each dimension up to the dimension of
436 : // the current mesh minus one (because this is for lower-dimensional
437 : // elements!)
438 0 : unsigned int dim = ((type.family == LAGRANGE_VEC) || (type.family == MONOMIAL_VEC)) ? 0 : 2;
439 0 : const auto ending_dim = cast_int<unsigned int>(_mesh_dimension - 1);
440 0 : if (ending_dim < dim)
441 0 : return;
442 0 : for (; dim <= ending_dim; dim++)
443 : {
444 0 : if (!_vector_fe_lower[dim][type])
445 0 : _vector_fe_lower[dim][type] = FEVectorBase::build(dim, type).release();
446 :
447 0 : _vector_fe_lower[dim][type]->get_dual_phi();
448 0 : _vector_fe_lower[dim][type]->get_dual_dphi();
449 0 : if (_need_second_derivative.count(type))
450 0 : _vector_fe_lower[dim][type]->get_dual_d2phi();
451 : }
452 : }
453 :
454 : void
455 529455 : Assembly::buildVectorFE(const FEType type) const
456 : {
457 529455 : if (!_vector_fe_shape_data[type])
458 1259 : _vector_fe_shape_data[type] = std::make_unique<VectorFEShapeData>();
459 :
460 : // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
461 : unsigned int min_dim;
462 529455 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
463 2968 : type.family == L2_RAVIART_THOMAS)
464 526795 : min_dim = 2;
465 : else
466 2660 : min_dim = 0;
467 :
468 : // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
469 : // current mesh
470 1520683 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
471 : {
472 991228 : if (!_vector_fe[dim][type])
473 3174 : _vector_fe[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
474 :
475 991228 : _vector_fe[dim][type]->get_phi();
476 991228 : _vector_fe[dim][type]->get_dphi();
477 991228 : if (_need_curl.count(type))
478 63894 : _vector_fe[dim][type]->get_curl_phi();
479 991228 : if (_need_div.count(type))
480 919204 : _vector_fe[dim][type]->get_div_phi();
481 991228 : _vector_fe[dim][type]->get_xyz();
482 : }
483 529455 : }
484 :
485 : void
486 3663 : Assembly::buildVectorFaceFE(const FEType type) const
487 : {
488 3663 : if (!_vector_fe_shape_data_face[type])
489 1259 : _vector_fe_shape_data_face[type] = std::make_unique<VectorFEShapeData>();
490 :
491 : // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
492 : unsigned int min_dim;
493 3663 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
494 2464 : type.family == L2_RAVIART_THOMAS)
495 1353 : min_dim = 2;
496 : else
497 2310 : min_dim = 0;
498 :
499 : // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
500 : // current mesh
501 12492 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
502 : {
503 8829 : if (!_vector_fe_face[dim][type])
504 3174 : _vector_fe_face[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
505 :
506 8829 : _vector_fe_face[dim][type]->get_phi();
507 8829 : _vector_fe_face[dim][type]->get_dphi();
508 8829 : if (_need_curl.count(type))
509 283 : _vector_fe_face[dim][type]->get_curl_phi();
510 8829 : if (_need_face_div.count(type))
511 416 : _vector_fe_face[dim][type]->get_div_phi();
512 : }
513 3663 : }
514 :
515 : void
516 3142 : Assembly::buildVectorNeighborFE(const FEType type) const
517 : {
518 3142 : if (!_vector_fe_shape_data_neighbor[type])
519 1259 : _vector_fe_shape_data_neighbor[type] = std::make_unique<VectorFEShapeData>();
520 :
521 : // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
522 : unsigned int min_dim;
523 3142 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
524 2464 : type.family == L2_RAVIART_THOMAS)
525 832 : min_dim = 2;
526 : else
527 2310 : min_dim = 0;
528 :
529 : // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
530 : // current mesh
531 11272 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
532 : {
533 8130 : if (!_vector_fe_neighbor[dim][type])
534 3174 : _vector_fe_neighbor[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
535 :
536 8130 : _vector_fe_neighbor[dim][type]->get_phi();
537 8130 : _vector_fe_neighbor[dim][type]->get_dphi();
538 8130 : if (_need_curl.count(type))
539 0 : _vector_fe_neighbor[dim][type]->get_curl_phi();
540 8130 : if (_need_neighbor_div.count(type))
541 0 : _vector_fe_neighbor[dim][type]->get_div_phi();
542 : }
543 3142 : }
544 :
545 : void
546 3663 : Assembly::buildVectorFaceNeighborFE(const FEType type) const
547 : {
548 3663 : if (!_vector_fe_shape_data_face_neighbor[type])
549 1259 : _vector_fe_shape_data_face_neighbor[type] = std::make_unique<VectorFEShapeData>();
550 :
551 : // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
552 : unsigned int min_dim;
553 3663 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
554 2464 : type.family == L2_RAVIART_THOMAS)
555 1353 : min_dim = 2;
556 : else
557 2310 : min_dim = 0;
558 :
559 : // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
560 : // current mesh
561 12492 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
562 : {
563 8829 : if (!_vector_fe_face_neighbor[dim][type])
564 3174 : _vector_fe_face_neighbor[dim][type] =
565 6348 : FEGenericBase<VectorValue<Real>>::build(dim, type).release();
566 :
567 8829 : _vector_fe_face_neighbor[dim][type]->get_phi();
568 8829 : _vector_fe_face_neighbor[dim][type]->get_dphi();
569 8829 : if (_need_curl.count(type))
570 283 : _vector_fe_face_neighbor[dim][type]->get_curl_phi();
571 8829 : if (_need_face_neighbor_div.count(type))
572 0 : _vector_fe_face_neighbor[dim][type]->get_div_phi();
573 : }
574 3663 : }
575 :
576 : void
577 90 : Assembly::bumpVolumeQRuleOrder(Order volume_order, SubdomainID block)
578 : {
579 90 : auto & qdefault = _qrules[Moose::ANY_BLOCK_ID];
580 : mooseAssert(qdefault.size() > 0, "default quadrature must be initialized before order bumps");
581 :
582 90 : unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
583 90 : auto & qvec = _qrules[block];
584 90 : if (qvec.size() != ndims || !qvec[0].vol)
585 52 : createQRules(qdefault[0].vol->type(),
586 26 : qdefault[0].arbitrary_vol->get_order(),
587 : volume_order,
588 26 : qdefault[0].face->get_order(),
589 : block);
590 64 : else if (qvec[0].vol->get_order() < volume_order)
591 0 : createQRules(qvec[0].vol->type(),
592 0 : qvec[0].arbitrary_vol->get_order(),
593 : volume_order,
594 0 : qvec[0].face->get_order(),
595 : block);
596 : // otherwise do nothing - quadrature order is already as high as requested
597 90 : }
598 :
599 : void
600 15 : Assembly::bumpAllQRuleOrder(Order order, SubdomainID block)
601 : {
602 15 : auto & qdefault = _qrules[Moose::ANY_BLOCK_ID];
603 : mooseAssert(qdefault.size() > 0, "default quadrature must be initialized before order bumps");
604 :
605 15 : unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
606 15 : auto & qvec = _qrules[block];
607 15 : if (qvec.size() != ndims || !qvec[0].vol)
608 13 : createQRules(qdefault[0].vol->type(), order, order, order, block);
609 2 : else if (qvec[0].vol->get_order() < order || qvec[0].face->get_order() < order)
610 0 : createQRules(qvec[0].vol->type(),
611 0 : std::max(order, qvec[0].arbitrary_vol->get_order()),
612 0 : std::max(order, qvec[0].vol->get_order()),
613 0 : std::max(order, qvec[0].face->get_order()),
614 : block);
615 : // otherwise do nothing - quadrature order is already as high as requested
616 15 : }
617 :
618 : void
619 66176 : Assembly::createQRules(QuadratureType type,
620 : Order order,
621 : Order volume_order,
622 : Order face_order,
623 : SubdomainID block,
624 : bool allow_negative_qweights)
625 : {
626 66176 : auto & qvec = _qrules[block];
627 66176 : unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
628 66176 : if (qvec.size() != ndims)
629 66176 : qvec.resize(ndims);
630 :
631 257077 : for (unsigned int i = 0; i < qvec.size(); i++)
632 : {
633 190901 : int dim = i;
634 190901 : auto & q = qvec[dim];
635 190901 : q.vol = QBase::build(type, dim, volume_order);
636 190901 : q.vol->allow_rules_with_negative_weights = allow_negative_qweights;
637 190901 : q.face = QBase::build(type, dim - 1, face_order);
638 190901 : q.face->allow_rules_with_negative_weights = allow_negative_qweights;
639 190901 : q.fv_face = QBase::build(QMONOMIAL, dim - 1, CONSTANT);
640 190901 : q.fv_face->allow_rules_with_negative_weights = allow_negative_qweights;
641 190901 : q.neighbor = std::make_unique<ArbitraryQuadrature>(dim - 1, face_order);
642 190901 : q.neighbor->allow_rules_with_negative_weights = allow_negative_qweights;
643 190901 : q.arbitrary_vol = std::make_unique<ArbitraryQuadrature>(dim, order);
644 190901 : q.arbitrary_vol->allow_rules_with_negative_weights = allow_negative_qweights;
645 190901 : q.arbitrary_face = std::make_unique<ArbitraryQuadrature>(dim - 1, face_order);
646 190901 : q.arbitrary_face->allow_rules_with_negative_weights = allow_negative_qweights;
647 : }
648 :
649 66176 : delete _qrule_msm;
650 66176 : _custom_mortar_qrule = false;
651 66176 : _qrule_msm = QBase::build(type, _mesh_dimension - 1, face_order).release();
652 66176 : _qrule_msm->allow_rules_with_negative_weights = allow_negative_qweights;
653 66176 : _fe_msm->attach_quadrature_rule(_qrule_msm);
654 66176 : }
655 :
656 : void
657 225641 : Assembly::setVolumeQRule(QBase * qrule, unsigned int dim)
658 : {
659 225641 : _current_qrule = qrule;
660 :
661 225641 : if (qrule) // Don't set a NULL qrule
662 : {
663 545192 : for (auto & it : _fe[dim])
664 319551 : it.second->attach_quadrature_rule(qrule);
665 229601 : for (auto & it : _vector_fe[dim])
666 3960 : it.second->attach_quadrature_rule(qrule);
667 225641 : if (!_unique_fe_helper.empty())
668 : {
669 : mooseAssert(dim < _unique_fe_helper.size(), "We should not be indexing out of bounds");
670 214 : _unique_fe_helper[dim]->attach_quadrature_rule(qrule);
671 : }
672 : }
673 225641 : }
674 :
675 : void
676 563870 : Assembly::setFaceQRule(QBase * qrule, unsigned int dim)
677 : {
678 563870 : _current_qrule_face = qrule;
679 :
680 1334156 : for (auto & it : _fe_face[dim])
681 770286 : it.second->attach_quadrature_rule(qrule);
682 564494 : for (auto & it : _vector_fe_face[dim])
683 624 : it.second->attach_quadrature_rule(qrule);
684 563870 : if (!_unique_fe_face_helper.empty())
685 : {
686 : mooseAssert(dim < _unique_fe_face_helper.size(), "We should not be indexing out of bounds");
687 200 : _unique_fe_face_helper[dim]->attach_quadrature_rule(qrule);
688 : }
689 563870 : }
690 :
691 : void
692 263515 : Assembly::setLowerQRule(QBase * qrule, unsigned int dim)
693 : {
694 : // The lower-dimensional quadrature rule matches the face quadrature rule
695 263515 : setFaceQRule(qrule, dim);
696 :
697 263515 : _current_qrule_lower = qrule;
698 :
699 616116 : for (auto & it : _fe_lower[dim])
700 352601 : it.second->attach_quadrature_rule(qrule);
701 263515 : for (auto & it : _vector_fe_lower[dim])
702 0 : it.second->attach_quadrature_rule(qrule);
703 263515 : if (!_unique_fe_lower_helper.empty())
704 : {
705 : mooseAssert(dim < _unique_fe_lower_helper.size(), "We should not be indexing out of bounds");
706 0 : _unique_fe_lower_helper[dim]->attach_quadrature_rule(qrule);
707 : }
708 263515 : }
709 :
710 : void
711 27021639 : Assembly::setNeighborQRule(QBase * qrule, unsigned int dim)
712 : {
713 27021639 : _current_qrule_neighbor = qrule;
714 :
715 80577481 : for (auto & it : _fe_face_neighbor[dim])
716 53555842 : it.second->attach_quadrature_rule(qrule);
717 27047504 : for (auto & it : _vector_fe_face_neighbor[dim])
718 25865 : it.second->attach_quadrature_rule(qrule);
719 27021639 : if (!_unique_fe_face_neighbor_helper.empty())
720 : {
721 : mooseAssert(dim < _unique_fe_face_neighbor_helper.size(),
722 : "We should not be indexing out of bounds");
723 62324 : _unique_fe_face_neighbor_helper[dim]->attach_quadrature_rule(qrule);
724 : }
725 27021639 : }
726 :
727 : void
728 62685 : Assembly::clearCachedQRules()
729 : {
730 62685 : _current_qrule = nullptr;
731 62685 : _current_qrule_face = nullptr;
732 62685 : _current_qrule_lower = nullptr;
733 62685 : _current_qrule_neighbor = nullptr;
734 62685 : }
735 :
736 : void
737 0 : Assembly::setMortarQRule(Order order)
738 : {
739 0 : if (order != _qrule_msm->get_order())
740 : {
741 : // If custom mortar qrule has not yet been specified
742 0 : if (!_custom_mortar_qrule)
743 : {
744 0 : _custom_mortar_qrule = true;
745 0 : const unsigned int dim = _qrule_msm->get_dim();
746 0 : const QuadratureType type = _qrule_msm->type();
747 0 : delete _qrule_msm;
748 :
749 0 : _qrule_msm = QBase::build(type, dim, order).release();
750 0 : _fe_msm->attach_quadrature_rule(_qrule_msm);
751 : }
752 : else
753 0 : mooseError("Mortar quadrature_order: ",
754 : order,
755 : " does not match previously specified quadrature_order: ",
756 0 : _qrule_msm->get_order(),
757 : ". Quadrature_order (when specified) must match for all mortar constraints.");
758 : }
759 0 : }
760 :
761 : void
762 394339425 : Assembly::reinitFE(const Elem * elem)
763 : {
764 394339425 : unsigned int dim = elem->dim();
765 :
766 888329751 : for (const auto & it : _fe[dim])
767 : {
768 493990354 : FEBase & fe = *it.second;
769 493990354 : const FEType & fe_type = it.first;
770 :
771 493990354 : _current_fe[fe_type] = &fe;
772 :
773 493990354 : FEShapeData & fesd = *_fe_shape_data[fe_type];
774 :
775 493990354 : fe.reinit(elem);
776 :
777 493990326 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_phi()));
778 493990326 : fesd._grad_phi.shallowCopy(
779 493990326 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_dphi()));
780 493990326 : if (_need_second_derivative.count(fe_type))
781 59941 : fesd._second_phi.shallowCopy(
782 59941 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_d2phi()));
783 : }
784 399959602 : for (const auto & it : _vector_fe[dim])
785 : {
786 5620205 : FEVectorBase & fe = *it.second;
787 5620205 : const FEType & fe_type = it.first;
788 :
789 5620205 : _current_vector_fe[fe_type] = &fe;
790 :
791 5620205 : VectorFEShapeData & fesd = *_vector_fe_shape_data[fe_type];
792 :
793 5620205 : fe.reinit(elem);
794 :
795 5620205 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_phi()));
796 5620205 : fesd._grad_phi.shallowCopy(
797 5620205 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_dphi()));
798 5620205 : if (_need_second_derivative.count(fe_type))
799 0 : fesd._second_phi.shallowCopy(
800 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe.get_d2phi()));
801 5620205 : if (_need_curl.count(fe_type))
802 2433834 : fesd._curl_phi.shallowCopy(
803 2433834 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_curl_phi()));
804 5620205 : if (_need_div.count(fe_type))
805 1228120 : fesd._div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_div_phi()));
806 : }
807 394339397 : if (!_unique_fe_helper.empty())
808 : {
809 : mooseAssert(dim < _unique_fe_helper.size(), "We should be in bounds here");
810 667582 : _unique_fe_helper[dim]->reinit(elem);
811 : }
812 :
813 : // During that last loop the helper objects will have been reinitialized as well
814 : // We need to dig out the q_points and JxW from it.
815 394339397 : _current_q_points.shallowCopy(
816 394339397 : const_cast<std::vector<Point> &>(_holder_fe_helper[dim]->get_xyz()));
817 394339397 : _current_JxW.shallowCopy(const_cast<std::vector<Real> &>(_holder_fe_helper[dim]->get_JxW()));
818 :
819 394339397 : if (_subproblem.haveADObjects())
820 : {
821 33195038 : auto n_qp = _current_qrule->n_points();
822 33195038 : resizeADMappingObjects(n_qp, dim);
823 33195038 : if (_displaced)
824 : {
825 575152 : const auto & qw = _current_qrule->get_weights();
826 3165489 : for (unsigned int qp = 0; qp != n_qp; qp++)
827 2590337 : computeSinglePointMapAD(elem, qw, qp, _holder_fe_helper[dim]);
828 : }
829 : else
830 : {
831 106563719 : for (unsigned qp = 0; qp < n_qp; ++qp)
832 73943833 : _ad_JxW[qp] = _current_JxW[qp];
833 32619886 : if (_calculate_xyz)
834 59813164 : for (unsigned qp = 0; qp < n_qp; ++qp)
835 47438565 : _ad_q_points[qp] = _current_q_points[qp];
836 : }
837 :
838 89399961 : for (const auto & it : _fe[dim])
839 : {
840 56204923 : FEBase & fe = *it.second;
841 56204923 : auto fe_type = it.first;
842 56204923 : auto num_shapes = FEInterface::n_shape_functions(fe_type, elem);
843 56204923 : auto & grad_phi = _ad_grad_phi_data[fe_type];
844 :
845 56204923 : grad_phi.resize(num_shapes);
846 219905495 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
847 163700572 : grad_phi[i].resize(n_qp);
848 :
849 56204923 : if (_displaced)
850 865825 : computeGradPhiAD(elem, n_qp, grad_phi, &fe);
851 : else
852 : {
853 55339098 : const auto & regular_grad_phi = _fe_shape_data[fe_type]->_grad_phi;
854 215959497 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
855 675698696 : for (unsigned qp = 0; qp < n_qp; ++qp)
856 515078297 : grad_phi[i][qp] = regular_grad_phi[i][qp];
857 : }
858 : }
859 36315248 : for (const auto & it : _vector_fe[dim])
860 : {
861 3120210 : FEVectorBase & fe = *it.second;
862 3120210 : auto fe_type = it.first;
863 3120210 : auto num_shapes = FEInterface::n_shape_functions(fe_type, elem);
864 3120210 : auto & grad_phi = _ad_vector_grad_phi_data[fe_type];
865 :
866 3120210 : grad_phi.resize(num_shapes);
867 18652890 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
868 15532680 : grad_phi[i].resize(n_qp);
869 :
870 3120210 : if (_displaced)
871 0 : computeGradPhiAD(elem, n_qp, grad_phi, &fe);
872 : else
873 : {
874 3120210 : const auto & regular_grad_phi = _vector_fe_shape_data[fe_type]->_grad_phi;
875 18652890 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
876 79399800 : for (unsigned qp = 0; qp < n_qp; ++qp)
877 63867120 : grad_phi[i][qp] = regular_grad_phi[i][qp];
878 : }
879 : }
880 : }
881 :
882 394339397 : auto n = numExtraElemIntegers();
883 397010005 : for (auto i : make_range(n))
884 2670608 : _extra_elem_ids[i] = _current_elem->get_extra_integer(i);
885 394339397 : _extra_elem_ids[n] = _current_elem->subdomain_id();
886 :
887 394339397 : if (_xfem != nullptr)
888 0 : modifyWeightsDueToXFEM(elem);
889 394339397 : }
890 :
891 : template <typename OutputType>
892 : void
893 997281 : Assembly::computeGradPhiAD(const Elem * elem,
894 : unsigned int n_qp,
895 : ADTemplateVariablePhiGradient<OutputType> & grad_phi,
896 : FEGenericBase<OutputType> * fe)
897 : {
898 : // This function relies on the fact that FE::reinit has already been called. FE::reinit will
899 : // importantly have already called FEMap::init_shape_functions which will have computed
900 : // these quantities at the integration/quadrature points: dphidxi,
901 : // dphideta, and dphidzeta (e.g. \nabla phi w.r.t. reference coordinates). These *phi* quantities
902 : // are independent of mesh displacements when using a quadrature rule.
903 : //
904 : // Note that a user could have specified custom integration points (e.g. independent of a
905 : // quadrature rule) which could very well depend on displacements. In that case even the *phi*
906 : // quantities from the above paragraph would be a function of the displacements and we would be
907 : // missing that derivative information in the calculations below
908 :
909 997281 : auto dim = elem->dim();
910 997281 : const auto & dphidxi = fe->get_dphidxi();
911 997281 : const auto & dphideta = fe->get_dphideta();
912 997281 : const auto & dphidzeta = fe->get_dphidzeta();
913 997281 : auto num_shapes = grad_phi.size();
914 :
915 997281 : switch (dim)
916 : {
917 0 : case 0:
918 : {
919 0 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
920 0 : for (unsigned qp = 0; qp < n_qp; ++qp)
921 0 : grad_phi[i][qp] = 0;
922 0 : break;
923 : }
924 :
925 51957 : case 1:
926 : {
927 139340 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
928 274853 : for (unsigned qp = 0; qp < n_qp; ++qp)
929 : {
930 187470 : grad_phi[i][qp].slice(0) = dphidxi[i][qp] * _ad_dxidx_map[qp];
931 187470 : grad_phi[i][qp].slice(1) = dphidxi[i][qp] * _ad_dxidy_map[qp];
932 187470 : grad_phi[i][qp].slice(2) = dphidxi[i][qp] * _ad_dxidz_map[qp];
933 : }
934 51957 : break;
935 : }
936 :
937 945324 : case 2:
938 : {
939 4624758 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
940 21534379 : for (unsigned qp = 0; qp < n_qp; ++qp)
941 : {
942 35709890 : grad_phi[i][qp].slice(0) =
943 35709890 : dphidxi[i][qp] * _ad_dxidx_map[qp] + dphideta[i][qp] * _ad_detadx_map[qp];
944 35709890 : grad_phi[i][qp].slice(1) =
945 35709890 : dphidxi[i][qp] * _ad_dxidy_map[qp] + dphideta[i][qp] * _ad_detady_map[qp];
946 35709890 : grad_phi[i][qp].slice(2) =
947 35709890 : dphidxi[i][qp] * _ad_dxidz_map[qp] + dphideta[i][qp] * _ad_detadz_map[qp];
948 : }
949 945324 : break;
950 : }
951 :
952 0 : case 3:
953 : {
954 0 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
955 0 : for (unsigned qp = 0; qp < n_qp; ++qp)
956 : {
957 0 : grad_phi[i][qp].slice(0) = dphidxi[i][qp] * _ad_dxidx_map[qp] +
958 0 : dphideta[i][qp] * _ad_detadx_map[qp] +
959 0 : dphidzeta[i][qp] * _ad_dzetadx_map[qp];
960 0 : grad_phi[i][qp].slice(1) = dphidxi[i][qp] * _ad_dxidy_map[qp] +
961 0 : dphideta[i][qp] * _ad_detady_map[qp] +
962 0 : dphidzeta[i][qp] * _ad_dzetady_map[qp];
963 0 : grad_phi[i][qp].slice(2) = dphidxi[i][qp] * _ad_dxidz_map[qp] +
964 0 : dphideta[i][qp] * _ad_detadz_map[qp] +
965 0 : dphidzeta[i][qp] * _ad_dzetadz_map[qp];
966 : }
967 0 : break;
968 : }
969 : }
970 997281 : }
971 :
972 : void
973 37284428 : Assembly::resizeADMappingObjects(unsigned int n_qp, unsigned int dim)
974 : {
975 37284428 : _ad_dxyzdxi_map.resize(n_qp);
976 37284428 : _ad_dxidx_map.resize(n_qp);
977 37284428 : _ad_dxidy_map.resize(n_qp); // 1D element may live in 2D ...
978 37284428 : _ad_dxidz_map.resize(n_qp); // ... or 3D
979 :
980 37284428 : if (dim > 1)
981 : {
982 22752064 : _ad_dxyzdeta_map.resize(n_qp);
983 22752064 : _ad_detadx_map.resize(n_qp);
984 22752064 : _ad_detady_map.resize(n_qp);
985 22752064 : _ad_detadz_map.resize(n_qp);
986 :
987 22752064 : if (dim > 2)
988 : {
989 468926 : _ad_dxyzdzeta_map.resize(n_qp);
990 468926 : _ad_dzetadx_map.resize(n_qp);
991 468926 : _ad_dzetady_map.resize(n_qp);
992 468926 : _ad_dzetadz_map.resize(n_qp);
993 : }
994 : }
995 :
996 37284428 : _ad_jac.resize(n_qp);
997 37284428 : _ad_JxW.resize(n_qp);
998 37284428 : if (_calculate_xyz)
999 16407608 : _ad_q_points.resize(n_qp);
1000 37284428 : }
1001 :
1002 : void
1003 2697665 : Assembly::computeSinglePointMapAD(const Elem * elem,
1004 : const std::vector<Real> & qw,
1005 : unsigned p,
1006 : FEBase * fe)
1007 : {
1008 : // This function relies on the fact that FE::reinit has already been called. FE::reinit will
1009 : // importantly have already called FEMap::init_reference_to_physical_map which will have computed
1010 : // these quantities at the integration/quadrature points: phi_map, dphidxi_map,
1011 : // dphideta_map, and dphidzeta_map (e.g. phi and \nabla phi w.r.t reference coordinates). *_map is
1012 : // used to denote that quantities are in reference to a mapping Lagrange FE object. The FE<Dim,
1013 : // LAGRANGE> objects used for mapping will in general have an order matching the order of the
1014 : // mesh. These *phi*_map quantities are independent of mesh displacements when using a quadrature
1015 : // rule.
1016 : //
1017 : // Note that a user could have specified custom integration points (e.g. independent of a
1018 : // quadrature rule) which could very well depend on displacements. In that case even the *phi*_map
1019 : // quantities from the above paragraph would be a function of the displacements and we would be
1020 : // missing that derivative information in the calculations below
1021 : //
1022 : // Important quantities calculated by this method:
1023 : // - _ad_JxW;
1024 : // - _ad_q_points;
1025 : // And the following quantities are important because they are used in the computeGradPhiAD method
1026 : // to calculate the shape function gradients with respect to the physical coordinates
1027 : // dphi/dphys = dphi/dref * dref/dphys:
1028 : // - _ad_dxidx_map;
1029 : // - _ad_dxidy_map;
1030 : // - _ad_dxidz_map;
1031 : // - _ad_detadx_map;
1032 : // - _ad_detady_map;
1033 : // - _ad_detadz_map;
1034 : // - _ad_dzetadx_map;
1035 : // - _ad_dzetady_map;
1036 : // - _ad_dzetadz_map;
1037 : //
1038 : // Some final notes. This method will be called both when we are reinit'ing in the volume and on
1039 : // faces. When reinit'ing on faces, computation of _ad_JxW will be garbage because we will be
1040 : // using dummy quadrature weights. _ad_q_points computation is also currently extraneous during
1041 : // face reinit because we compute _ad_q_points_face in the computeFaceMap method. However,
1042 : // computation of dref/dphys is absolutely necessary (and the reason we call this method for the
1043 : // face case) for both volume and face reinit
1044 :
1045 2697665 : auto dim = elem->dim();
1046 2697665 : const auto & elem_nodes = elem->get_nodes();
1047 2697665 : auto num_shapes = FEInterface::n_shape_functions(fe->get_fe_type(), elem);
1048 2697665 : const auto & phi_map = fe->get_fe_map().get_phi_map();
1049 2697665 : const auto & dphidxi_map = fe->get_fe_map().get_dphidxi_map();
1050 2697665 : const auto & dphideta_map = fe->get_fe_map().get_dphideta_map();
1051 2697665 : const auto & dphidzeta_map = fe->get_fe_map().get_dphidzeta_map();
1052 2697665 : const auto sys_num = _sys.number();
1053 : const bool do_derivatives =
1054 2697665 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
1055 :
1056 2697665 : switch (dim)
1057 : {
1058 0 : case 0:
1059 : {
1060 0 : _ad_jac[p] = 1.0;
1061 0 : _ad_JxW[p] = qw[p];
1062 0 : if (_calculate_xyz)
1063 0 : _ad_q_points[p] = *elem_nodes[0];
1064 0 : break;
1065 : }
1066 :
1067 62684 : case 1:
1068 : {
1069 62684 : if (_calculate_xyz)
1070 12502 : _ad_q_points[p].zero();
1071 :
1072 62684 : _ad_dxyzdxi_map[p].zero();
1073 :
1074 195972 : for (std::size_t i = 0; i < num_shapes; i++)
1075 : {
1076 : libmesh_assert(elem_nodes[i]);
1077 133288 : const Node & node = *elem_nodes[i];
1078 133288 : libMesh::VectorValue<ADReal> elem_point = node;
1079 133288 : if (do_derivatives)
1080 50772 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1081 2714 : if (node.n_dofs(sys_num, disp_num))
1082 5428 : Moose::derivInsert(
1083 2714 : elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1084 :
1085 133288 : _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1086 :
1087 133288 : if (_calculate_xyz)
1088 32924 : _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1089 133288 : }
1090 :
1091 62684 : _ad_jac[p] = _ad_dxyzdxi_map[p].norm();
1092 :
1093 62684 : if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
1094 : {
1095 : static bool failing = false;
1096 0 : if (!failing)
1097 : {
1098 0 : failing = true;
1099 0 : elem->print_info(libMesh::err);
1100 0 : libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
1101 : << p << " in element " << elem->id());
1102 : }
1103 : else
1104 0 : return;
1105 : }
1106 :
1107 62684 : const auto jacm2 = 1. / _ad_jac[p] / _ad_jac[p];
1108 62684 : _ad_dxidx_map[p] = jacm2 * _ad_dxyzdxi_map[p](0);
1109 62684 : _ad_dxidy_map[p] = jacm2 * _ad_dxyzdxi_map[p](1);
1110 62684 : _ad_dxidz_map[p] = jacm2 * _ad_dxyzdxi_map[p](2);
1111 :
1112 62684 : _ad_JxW[p] = _ad_jac[p] * qw[p];
1113 :
1114 62684 : break;
1115 62684 : }
1116 :
1117 2634981 : case 2:
1118 : {
1119 2634981 : if (_calculate_xyz)
1120 1903561 : _ad_q_points[p].zero();
1121 2634981 : _ad_dxyzdxi_map[p].zero();
1122 2634981 : _ad_dxyzdeta_map[p].zero();
1123 :
1124 16203878 : for (std::size_t i = 0; i < num_shapes; i++)
1125 : {
1126 : libmesh_assert(elem_nodes[i]);
1127 13568897 : const Node & node = *elem_nodes[i];
1128 13568897 : libMesh::VectorValue<ADReal> elem_point = node;
1129 13568897 : if (do_derivatives)
1130 1505152 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1131 146496 : if (node.n_dofs(sys_num, disp_num))
1132 292992 : Moose::derivInsert(
1133 146496 : elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1134 :
1135 13568897 : _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1136 13568897 : _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
1137 :
1138 13568897 : if (_calculate_xyz)
1139 11200209 : _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1140 13568897 : }
1141 :
1142 2634981 : const auto &dx_dxi = _ad_dxyzdxi_map[p](0), &dx_deta = _ad_dxyzdeta_map[p](0),
1143 2634981 : &dy_dxi = _ad_dxyzdxi_map[p](1), &dy_deta = _ad_dxyzdeta_map[p](1),
1144 2634981 : &dz_dxi = _ad_dxyzdxi_map[p](2), &dz_deta = _ad_dxyzdeta_map[p](2);
1145 :
1146 2634981 : const auto g11 = (dx_dxi * dx_dxi + dy_dxi * dy_dxi + dz_dxi * dz_dxi);
1147 :
1148 2634981 : const auto g12 = (dx_dxi * dx_deta + dy_dxi * dy_deta + dz_dxi * dz_deta);
1149 :
1150 2634981 : const auto & g21 = g12;
1151 :
1152 2634981 : const auto g22 = (dx_deta * dx_deta + dy_deta * dy_deta + dz_deta * dz_deta);
1153 :
1154 2634981 : auto det = (g11 * g22 - g12 * g21);
1155 :
1156 2634981 : if (det.value() <= -TOLERANCE * TOLERANCE)
1157 : {
1158 : static bool failing = false;
1159 0 : if (!failing)
1160 : {
1161 0 : failing = true;
1162 0 : elem->print_info(libMesh::err);
1163 0 : libmesh_error_msg("ERROR: negative Jacobian " << det << " at point index " << p
1164 : << " in element " << elem->id());
1165 : }
1166 : else
1167 0 : return;
1168 : }
1169 2634981 : else if (det.value() <= 0.)
1170 0 : det.value() = TOLERANCE * TOLERANCE;
1171 :
1172 2634981 : const auto inv_det = 1. / det;
1173 2634981 : _ad_jac[p] = std::sqrt(det);
1174 :
1175 2634981 : _ad_JxW[p] = _ad_jac[p] * qw[p];
1176 :
1177 2634981 : const auto g11inv = g22 * inv_det;
1178 2634981 : const auto g12inv = -g12 * inv_det;
1179 2634981 : const auto g21inv = -g21 * inv_det;
1180 2634981 : const auto g22inv = g11 * inv_det;
1181 :
1182 2634981 : _ad_dxidx_map[p] = g11inv * dx_dxi + g12inv * dx_deta;
1183 2634981 : _ad_dxidy_map[p] = g11inv * dy_dxi + g12inv * dy_deta;
1184 2634981 : _ad_dxidz_map[p] = g11inv * dz_dxi + g12inv * dz_deta;
1185 :
1186 2634981 : _ad_detadx_map[p] = g21inv * dx_dxi + g22inv * dx_deta;
1187 2634981 : _ad_detady_map[p] = g21inv * dy_dxi + g22inv * dy_deta;
1188 2634981 : _ad_detadz_map[p] = g21inv * dz_dxi + g22inv * dz_deta;
1189 :
1190 2634981 : break;
1191 13174905 : }
1192 :
1193 0 : case 3:
1194 : {
1195 0 : if (_calculate_xyz)
1196 0 : _ad_q_points[p].zero();
1197 0 : _ad_dxyzdxi_map[p].zero();
1198 0 : _ad_dxyzdeta_map[p].zero();
1199 0 : _ad_dxyzdzeta_map[p].zero();
1200 :
1201 0 : for (std::size_t i = 0; i < num_shapes; i++)
1202 : {
1203 : libmesh_assert(elem_nodes[i]);
1204 0 : const Node & node = *elem_nodes[i];
1205 0 : libMesh::VectorValue<ADReal> elem_point = node;
1206 0 : if (do_derivatives)
1207 0 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1208 0 : if (node.n_dofs(sys_num, disp_num))
1209 0 : Moose::derivInsert(
1210 0 : elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1211 :
1212 0 : _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1213 0 : _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
1214 0 : _ad_dxyzdzeta_map[p].add_scaled(elem_point, dphidzeta_map[i][p]);
1215 :
1216 0 : if (_calculate_xyz)
1217 0 : _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1218 0 : }
1219 :
1220 0 : const auto &dx_dxi = _ad_dxyzdxi_map[p](0), &dy_dxi = _ad_dxyzdxi_map[p](1),
1221 0 : &dz_dxi = _ad_dxyzdxi_map[p](2), &dx_deta = _ad_dxyzdeta_map[p](0),
1222 0 : &dy_deta = _ad_dxyzdeta_map[p](1), &dz_deta = _ad_dxyzdeta_map[p](2),
1223 0 : &dx_dzeta = _ad_dxyzdzeta_map[p](0), &dy_dzeta = _ad_dxyzdzeta_map[p](1),
1224 0 : &dz_dzeta = _ad_dxyzdzeta_map[p](2);
1225 :
1226 0 : _ad_jac[p] = (dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) +
1227 0 : dy_dxi * (dz_deta * dx_dzeta - dx_deta * dz_dzeta) +
1228 0 : dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta));
1229 :
1230 0 : if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
1231 : {
1232 : static bool failing = false;
1233 0 : if (!failing)
1234 : {
1235 0 : failing = true;
1236 0 : elem->print_info(libMesh::err);
1237 0 : libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
1238 : << p << " in element " << elem->id());
1239 : }
1240 : else
1241 0 : return;
1242 : }
1243 :
1244 0 : _ad_JxW[p] = _ad_jac[p] * qw[p];
1245 :
1246 0 : const auto inv_jac = 1. / _ad_jac[p];
1247 :
1248 0 : _ad_dxidx_map[p] = (dy_deta * dz_dzeta - dz_deta * dy_dzeta) * inv_jac;
1249 0 : _ad_dxidy_map[p] = (dz_deta * dx_dzeta - dx_deta * dz_dzeta) * inv_jac;
1250 0 : _ad_dxidz_map[p] = (dx_deta * dy_dzeta - dy_deta * dx_dzeta) * inv_jac;
1251 :
1252 0 : _ad_detadx_map[p] = (dz_dxi * dy_dzeta - dy_dxi * dz_dzeta) * inv_jac;
1253 0 : _ad_detady_map[p] = (dx_dxi * dz_dzeta - dz_dxi * dx_dzeta) * inv_jac;
1254 0 : _ad_detadz_map[p] = (dy_dxi * dx_dzeta - dx_dxi * dy_dzeta) * inv_jac;
1255 :
1256 0 : _ad_dzetadx_map[p] = (dy_dxi * dz_deta - dz_dxi * dy_deta) * inv_jac;
1257 0 : _ad_dzetady_map[p] = (dz_dxi * dx_deta - dx_dxi * dz_deta) * inv_jac;
1258 0 : _ad_dzetadz_map[p] = (dx_dxi * dy_deta - dy_dxi * dx_deta) * inv_jac;
1259 :
1260 0 : break;
1261 0 : }
1262 :
1263 0 : default:
1264 0 : libmesh_error_msg("Invalid dim = " << dim);
1265 : }
1266 : }
1267 :
1268 : void
1269 11953620 : Assembly::reinitFEFace(const Elem * elem, unsigned int side)
1270 : {
1271 11953620 : unsigned int dim = elem->dim();
1272 :
1273 37969179 : for (const auto & it : _fe_face[dim])
1274 : {
1275 26015559 : FEBase & fe_face = *it.second;
1276 26015559 : const FEType & fe_type = it.first;
1277 26015559 : FEShapeData & fesd = *_fe_shape_data_face[fe_type];
1278 26015559 : fe_face.reinit(elem, side);
1279 26015559 : _current_fe_face[fe_type] = &fe_face;
1280 :
1281 26015559 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face.get_phi()));
1282 26015559 : fesd._grad_phi.shallowCopy(
1283 26015559 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_dphi()));
1284 26015559 : if (_need_second_derivative.count(fe_type))
1285 16216 : fesd._second_phi.shallowCopy(
1286 16216 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
1287 : }
1288 14723830 : for (const auto & it : _vector_fe_face[dim])
1289 : {
1290 2770210 : FEVectorBase & fe_face = *it.second;
1291 2770210 : const FEType & fe_type = it.first;
1292 :
1293 2770210 : _current_vector_fe_face[fe_type] = &fe_face;
1294 :
1295 2770210 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face[fe_type];
1296 :
1297 2770210 : fe_face.reinit(elem, side);
1298 :
1299 2770210 : fesd._phi.shallowCopy(
1300 2770210 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_phi()));
1301 2770210 : fesd._grad_phi.shallowCopy(
1302 2770210 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_dphi()));
1303 2770210 : if (_need_second_derivative.count(fe_type))
1304 0 : fesd._second_phi.shallowCopy(
1305 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
1306 2770210 : if (_need_curl.count(fe_type))
1307 486265 : fesd._curl_phi.shallowCopy(
1308 486265 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
1309 2770210 : if (_need_face_div.count(fe_type))
1310 37632 : fesd._div_phi.shallowCopy(
1311 37632 : const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
1312 : }
1313 11953620 : if (!_unique_fe_face_helper.empty())
1314 : {
1315 : mooseAssert(dim < _unique_fe_face_helper.size(), "We should be in bounds here");
1316 107270 : _unique_fe_face_helper[dim]->reinit(elem, side);
1317 : }
1318 :
1319 : // During that last loop the helper objects will have been reinitialized as well
1320 : // We need to dig out the q_points and JxW from it.
1321 11953620 : _current_q_points_face.shallowCopy(
1322 11953620 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_xyz()));
1323 11953620 : _current_JxW_face.shallowCopy(
1324 11953620 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_JxW()));
1325 11953620 : _current_normals.shallowCopy(
1326 11953620 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_normals()));
1327 :
1328 11953620 : _mapped_normals.resize(_current_normals.size(), Eigen::Map<RealDIMValue>(nullptr));
1329 38984236 : for (unsigned int i = 0; i < _current_normals.size(); i++)
1330 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
1331 27030616 : new (&_mapped_normals[i]) Eigen::Map<RealDIMValue>(const_cast<Real *>(&_current_normals[i](0)));
1332 :
1333 11953620 : if (_calculate_curvatures)
1334 416 : _curvatures.shallowCopy(
1335 416 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_curvatures()));
1336 :
1337 11953620 : computeADFace(*elem, side);
1338 :
1339 11953620 : if (_xfem != nullptr)
1340 0 : modifyFaceWeightsDueToXFEM(elem, side);
1341 :
1342 11953620 : auto n = numExtraElemIntegers();
1343 12096740 : for (auto i : make_range(n))
1344 143120 : _extra_elem_ids[i] = _current_elem->get_extra_integer(i);
1345 11953620 : _extra_elem_ids[n] = _current_elem->subdomain_id();
1346 11953620 : }
1347 :
1348 : void
1349 53788 : Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std::vector<Real> & qw)
1350 : {
1351 : // Important quantities calculated by this method:
1352 : // - _ad_JxW_face
1353 : // - _ad_q_points_face
1354 : // - _ad_normals
1355 : // - _ad_curvatures
1356 :
1357 53788 : const Elem & side_elem = _compute_face_map_side_elem_builder(elem, side);
1358 53788 : const auto dim = elem.dim();
1359 53788 : const auto n_qp = qw.size();
1360 53788 : const auto & dpsidxi_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsidxi();
1361 53788 : const auto & dpsideta_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsideta();
1362 53788 : const auto & psi_map = _holder_fe_face_helper[dim]->get_fe_map().get_psi();
1363 53788 : std::vector<std::vector<Real>> const * d2psidxi2_map = nullptr;
1364 53788 : std::vector<std::vector<Real>> const * d2psidxideta_map = nullptr;
1365 53788 : std::vector<std::vector<Real>> const * d2psideta2_map = nullptr;
1366 53788 : const auto sys_num = _sys.number();
1367 53788 : const bool do_derivatives = ADReal::do_derivatives && sys_num == _subproblem.currentNlSysNum();
1368 :
1369 53788 : if (_calculate_curvatures)
1370 : {
1371 0 : d2psidxi2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxi2();
1372 0 : d2psidxideta_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxideta();
1373 0 : d2psideta2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psideta2();
1374 : }
1375 :
1376 53788 : switch (dim)
1377 : {
1378 248 : case 1:
1379 : {
1380 248 : if (!n_qp)
1381 0 : break;
1382 :
1383 248 : if (side_elem.node_id(0) == elem.node_id(0))
1384 248 : _ad_normals[0] = Point(-1.);
1385 : else
1386 0 : _ad_normals[0] = Point(1.);
1387 :
1388 248 : VectorValue<ADReal> side_point;
1389 248 : if (_calculate_face_xyz)
1390 : {
1391 248 : const Node & node = side_elem.node_ref(0);
1392 248 : side_point = node;
1393 :
1394 248 : if (do_derivatives)
1395 112 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1396 112 : Moose::derivInsert(
1397 56 : side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1398 : }
1399 :
1400 496 : for (const auto p : make_range(n_qp))
1401 : {
1402 248 : if (_calculate_face_xyz)
1403 : {
1404 248 : _ad_q_points_face[p].zero();
1405 248 : _ad_q_points_face[p].add_scaled(side_point, psi_map[0][p]);
1406 : }
1407 :
1408 248 : _ad_normals[p] = _ad_normals[0];
1409 248 : _ad_JxW_face[p] = 1.0 * qw[p];
1410 : }
1411 :
1412 248 : break;
1413 248 : }
1414 :
1415 53540 : case 2:
1416 : {
1417 53540 : _ad_dxyzdxi_map.resize(n_qp);
1418 53540 : if (_calculate_curvatures)
1419 0 : _ad_d2xyzdxi2_map.resize(n_qp);
1420 :
1421 160620 : for (const auto p : make_range(n_qp))
1422 107080 : _ad_dxyzdxi_map[p].zero();
1423 53540 : if (_calculate_face_xyz)
1424 99504 : for (const auto p : make_range(n_qp))
1425 66336 : _ad_q_points_face[p].zero();
1426 53540 : if (_calculate_curvatures)
1427 0 : for (const auto p : make_range(n_qp))
1428 0 : _ad_d2xyzdxi2_map[p].zero();
1429 :
1430 : const auto n_mapping_shape_functions =
1431 53540 : FE<2, LAGRANGE>::n_dofs(&side_elem, side_elem.default_order());
1432 :
1433 189420 : for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1434 : {
1435 135880 : const Node & node = side_elem.node_ref(i);
1436 135880 : VectorValue<ADReal> side_point = node;
1437 :
1438 135880 : if (do_derivatives)
1439 38272 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1440 0 : Moose::derivInsert(
1441 0 : side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1442 :
1443 407640 : for (const auto p : make_range(n_qp))
1444 271760 : _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1445 135880 : if (_calculate_face_xyz)
1446 285408 : for (const auto p : make_range(n_qp))
1447 190272 : _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1448 135880 : if (_calculate_curvatures)
1449 0 : for (const auto p : make_range(n_qp))
1450 0 : _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1451 135880 : }
1452 :
1453 160620 : for (const auto p : make_range(n_qp))
1454 : {
1455 107080 : _ad_normals[p] =
1456 214160 : (VectorValue<ADReal>(_ad_dxyzdxi_map[p](1), -_ad_dxyzdxi_map[p](0), 0.)).unit();
1457 107080 : const auto the_jac = _ad_dxyzdxi_map[p].norm();
1458 107080 : _ad_JxW_face[p] = the_jac * qw[p];
1459 107080 : if (_calculate_curvatures)
1460 : {
1461 0 : const auto numerator = _ad_d2xyzdxi2_map[p] * _ad_normals[p];
1462 0 : const auto denominator = _ad_dxyzdxi_map[p].norm_sq();
1463 : libmesh_assert_not_equal_to(denominator, 0);
1464 0 : _ad_curvatures[p] = numerator / denominator;
1465 0 : }
1466 107080 : }
1467 :
1468 53540 : break;
1469 : }
1470 :
1471 0 : case 3:
1472 : {
1473 0 : _ad_dxyzdxi_map.resize(n_qp);
1474 0 : _ad_dxyzdeta_map.resize(n_qp);
1475 0 : if (_calculate_curvatures)
1476 : {
1477 0 : _ad_d2xyzdxi2_map.resize(n_qp);
1478 0 : _ad_d2xyzdxideta_map.resize(n_qp);
1479 0 : _ad_d2xyzdeta2_map.resize(n_qp);
1480 : }
1481 :
1482 0 : for (const auto p : make_range(n_qp))
1483 : {
1484 0 : _ad_dxyzdxi_map[p].zero();
1485 0 : _ad_dxyzdeta_map[p].zero();
1486 : }
1487 0 : if (_calculate_face_xyz)
1488 0 : for (const auto p : make_range(n_qp))
1489 0 : _ad_q_points_face[p].zero();
1490 0 : if (_calculate_curvatures)
1491 0 : for (const auto p : make_range(n_qp))
1492 : {
1493 0 : _ad_d2xyzdxi2_map[p].zero();
1494 0 : _ad_d2xyzdxideta_map[p].zero();
1495 0 : _ad_d2xyzdeta2_map[p].zero();
1496 : }
1497 :
1498 : const unsigned int n_mapping_shape_functions =
1499 0 : FE<3, LAGRANGE>::n_dofs(&side_elem, side_elem.default_order());
1500 :
1501 0 : for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1502 : {
1503 0 : const Node & node = side_elem.node_ref(i);
1504 0 : VectorValue<ADReal> side_point = node;
1505 :
1506 0 : if (do_derivatives)
1507 0 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1508 0 : Moose::derivInsert(
1509 0 : side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1510 :
1511 0 : for (const auto p : make_range(n_qp))
1512 : {
1513 0 : _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1514 0 : _ad_dxyzdeta_map[p].add_scaled(side_point, dpsideta_map[i][p]);
1515 : }
1516 0 : if (_calculate_face_xyz)
1517 0 : for (const auto p : make_range(n_qp))
1518 0 : _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1519 0 : if (_calculate_curvatures)
1520 0 : for (const auto p : make_range(n_qp))
1521 : {
1522 0 : _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1523 0 : _ad_d2xyzdxideta_map[p].add_scaled(side_point, (*d2psidxideta_map)[i][p]);
1524 0 : _ad_d2xyzdeta2_map[p].add_scaled(side_point, (*d2psideta2_map)[i][p]);
1525 : }
1526 0 : }
1527 :
1528 0 : for (const auto p : make_range(n_qp))
1529 : {
1530 0 : _ad_normals[p] = _ad_dxyzdxi_map[p].cross(_ad_dxyzdeta_map[p]).unit();
1531 :
1532 0 : const auto &dxdxi = _ad_dxyzdxi_map[p](0), &dxdeta = _ad_dxyzdeta_map[p](0),
1533 0 : &dydxi = _ad_dxyzdxi_map[p](1), &dydeta = _ad_dxyzdeta_map[p](1),
1534 0 : &dzdxi = _ad_dxyzdxi_map[p](2), &dzdeta = _ad_dxyzdeta_map[p](2);
1535 :
1536 0 : const auto g11 = (dxdxi * dxdxi + dydxi * dydxi + dzdxi * dzdxi);
1537 :
1538 0 : const auto g12 = (dxdxi * dxdeta + dydxi * dydeta + dzdxi * dzdeta);
1539 :
1540 0 : const auto & g21 = g12;
1541 :
1542 0 : const auto g22 = (dxdeta * dxdeta + dydeta * dydeta + dzdeta * dzdeta);
1543 :
1544 0 : const auto the_jac = std::sqrt(g11 * g22 - g12 * g21);
1545 :
1546 0 : _ad_JxW_face[p] = the_jac * qw[p];
1547 :
1548 0 : if (_calculate_curvatures)
1549 : {
1550 0 : const auto L = -_ad_d2xyzdxi2_map[p] * _ad_normals[p];
1551 0 : const auto M = -_ad_d2xyzdxideta_map[p] * _ad_normals[p];
1552 0 : const auto N = -_ad_d2xyzdeta2_map[p] * _ad_normals[p];
1553 0 : const auto E = _ad_dxyzdxi_map[p].norm_sq();
1554 0 : const auto F = _ad_dxyzdxi_map[p] * _ad_dxyzdeta_map[p];
1555 0 : const auto G = _ad_dxyzdeta_map[p].norm_sq();
1556 :
1557 0 : const auto numerator = E * N - 2. * F * M + G * L;
1558 0 : const auto denominator = E * G - F * F;
1559 : libmesh_assert_not_equal_to(denominator, 0.);
1560 0 : _ad_curvatures[p] = 0.5 * numerator / denominator;
1561 0 : }
1562 0 : }
1563 :
1564 0 : break;
1565 : }
1566 :
1567 0 : default:
1568 0 : mooseError("Invalid dimension dim = ", dim);
1569 : }
1570 53788 : }
1571 :
1572 : void
1573 4024111 : Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1574 : {
1575 4024111 : unsigned int neighbor_dim = neighbor->dim();
1576 :
1577 : // reinit neighbor face
1578 11784063 : for (const auto & it : _fe_face_neighbor[neighbor_dim])
1579 : {
1580 7759952 : FEBase & fe_face_neighbor = *it.second;
1581 7759952 : FEType fe_type = it.first;
1582 7759952 : FEShapeData & fesd = *_fe_shape_data_face_neighbor[fe_type];
1583 :
1584 7759952 : fe_face_neighbor.reinit(neighbor, &reference_points);
1585 :
1586 7759952 : _current_fe_face_neighbor[fe_type] = &fe_face_neighbor;
1587 :
1588 7759952 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_phi()));
1589 7759952 : fesd._grad_phi.shallowCopy(
1590 7759952 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor.get_dphi()));
1591 7759952 : if (_need_second_derivative_neighbor.count(fe_type))
1592 8640 : fesd._second_phi.shallowCopy(
1593 8640 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
1594 : }
1595 4049976 : for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
1596 : {
1597 25865 : FEVectorBase & fe_face_neighbor = *it.second;
1598 25865 : const FEType & fe_type = it.first;
1599 :
1600 25865 : _current_vector_fe_face_neighbor[fe_type] = &fe_face_neighbor;
1601 :
1602 25865 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face_neighbor[fe_type];
1603 :
1604 25865 : fe_face_neighbor.reinit(neighbor, &reference_points);
1605 :
1606 25865 : fesd._phi.shallowCopy(
1607 25865 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor.get_phi()));
1608 25865 : fesd._grad_phi.shallowCopy(
1609 25865 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_dphi()));
1610 25865 : if (_need_second_derivative.count(fe_type))
1611 0 : fesd._second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
1612 0 : fe_face_neighbor.get_d2phi()));
1613 25865 : if (_need_curl.count(fe_type))
1614 105 : fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
1615 105 : fe_face_neighbor.get_curl_phi()));
1616 25865 : if (_need_face_neighbor_div.count(fe_type))
1617 0 : fesd._div_phi.shallowCopy(
1618 0 : const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
1619 : }
1620 4024111 : if (!_unique_fe_face_neighbor_helper.empty())
1621 : {
1622 : mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
1623 : "We should be in bounds here");
1624 62324 : _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
1625 : }
1626 :
1627 4024111 : _current_q_points_face_neighbor.shallowCopy(
1628 4024111 : const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
1629 4024111 : }
1630 :
1631 : void
1632 19880 : Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1633 : {
1634 19880 : unsigned int neighbor_dim = neighbor->dim();
1635 :
1636 : // reinit neighbor face
1637 39760 : for (const auto & it : _fe_neighbor[neighbor_dim])
1638 : {
1639 19880 : FEBase & fe_neighbor = *it.second;
1640 19880 : FEType fe_type = it.first;
1641 19880 : FEShapeData & fesd = *_fe_shape_data_neighbor[fe_type];
1642 :
1643 19880 : fe_neighbor.reinit(neighbor, &reference_points);
1644 :
1645 19880 : _current_fe_neighbor[fe_type] = &fe_neighbor;
1646 :
1647 19880 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_phi()));
1648 19880 : fesd._grad_phi.shallowCopy(
1649 19880 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_neighbor.get_dphi()));
1650 19880 : if (_need_second_derivative_neighbor.count(fe_type))
1651 0 : fesd._second_phi.shallowCopy(
1652 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_d2phi()));
1653 : }
1654 19880 : for (const auto & it : _vector_fe_neighbor[neighbor_dim])
1655 : {
1656 0 : FEVectorBase & fe_neighbor = *it.second;
1657 0 : const FEType & fe_type = it.first;
1658 :
1659 0 : _current_vector_fe_neighbor[fe_type] = &fe_neighbor;
1660 :
1661 0 : VectorFEShapeData & fesd = *_vector_fe_shape_data_neighbor[fe_type];
1662 :
1663 0 : fe_neighbor.reinit(neighbor, &reference_points);
1664 :
1665 0 : fesd._phi.shallowCopy(
1666 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_phi()));
1667 0 : fesd._grad_phi.shallowCopy(
1668 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_dphi()));
1669 0 : if (_need_second_derivative.count(fe_type))
1670 0 : fesd._second_phi.shallowCopy(
1671 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_neighbor.get_d2phi()));
1672 0 : if (_need_curl.count(fe_type))
1673 0 : fesd._curl_phi.shallowCopy(
1674 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_curl_phi()));
1675 0 : if (_need_neighbor_div.count(fe_type))
1676 0 : fesd._div_phi.shallowCopy(
1677 0 : const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_div_phi()));
1678 : }
1679 19880 : if (!_unique_fe_neighbor_helper.empty())
1680 : {
1681 : mooseAssert(neighbor_dim < _unique_fe_neighbor_helper.size(), "We should be in bounds here");
1682 0 : _unique_fe_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
1683 : }
1684 19880 : }
1685 :
1686 : void
1687 4043991 : Assembly::reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1688 : {
1689 4043991 : unsigned int neighbor_dim = neighbor->dim();
1690 : mooseAssert(_current_neighbor_subdomain_id == neighbor->subdomain_id(),
1691 : "Neighbor subdomain ID has not been correctly set");
1692 :
1693 : ArbitraryQuadrature * neighbor_rule =
1694 4043991 : qrules(neighbor_dim, _current_neighbor_subdomain_id).neighbor.get();
1695 4043991 : neighbor_rule->setPoints(reference_points);
1696 4043991 : setNeighborQRule(neighbor_rule, neighbor_dim);
1697 :
1698 4043991 : _current_neighbor_elem = neighbor;
1699 : mooseAssert(_current_neighbor_subdomain_id == _current_neighbor_elem->subdomain_id(),
1700 : "current neighbor subdomain has been set incorrectly");
1701 :
1702 : // Calculate the volume of the neighbor
1703 4043991 : if (_need_neighbor_elem_volume)
1704 : {
1705 2245686 : unsigned int dim = neighbor->dim();
1706 2245686 : FEBase & fe = *_holder_fe_neighbor_helper[dim];
1707 2245686 : QBase * qrule = qrules(dim).vol.get();
1708 :
1709 2245686 : fe.attach_quadrature_rule(qrule);
1710 2245686 : fe.reinit(neighbor);
1711 :
1712 2245686 : const std::vector<Real> & JxW = fe.get_JxW();
1713 2245686 : MooseArray<Point> q_points;
1714 2245686 : q_points.shallowCopy(const_cast<std::vector<Point> &>(fe.get_xyz()));
1715 :
1716 2245686 : setCoordinateTransformation(qrule, q_points, _coord_neighbor, _current_neighbor_subdomain_id);
1717 :
1718 2245686 : _current_neighbor_volume = 0.;
1719 19551269 : for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1720 17305583 : _current_neighbor_volume += JxW[qp] * _coord_neighbor[qp];
1721 2245686 : }
1722 :
1723 4043991 : auto n = numExtraElemIntegers();
1724 4124271 : for (auto i : make_range(n))
1725 80280 : _neighbor_extra_elem_ids[i] = _current_neighbor_elem->get_extra_integer(i);
1726 4043991 : _neighbor_extra_elem_ids[n] = _current_neighbor_elem->subdomain_id();
1727 4043991 : }
1728 :
1729 : template <typename Points, typename Coords>
1730 : void
1731 425457703 : Assembly::setCoordinateTransformation(const QBase * qrule,
1732 : const Points & q_points,
1733 : Coords & coord,
1734 : SubdomainID sub_id)
1735 : {
1736 :
1737 : mooseAssert(qrule, "The quadrature rule is null in Assembly::setCoordinateTransformation");
1738 425457703 : auto n_points = qrule->n_points();
1739 : mooseAssert(n_points == q_points.size(),
1740 : "The number of points in the quadrature rule doesn't match the number of passed-in "
1741 : "points in Assembly::setCoordinateTransformation");
1742 :
1743 : // Make sure to honor the name of this method and set the _coord_type member because users may
1744 : // make use of the const Moose::CoordinateSystem & coordTransformation() { return _coord_type; }
1745 : // API. MaterialBase for example uses it
1746 425457703 : _coord_type = _subproblem.getCoordSystem(sub_id);
1747 :
1748 425457703 : coord.resize(n_points);
1749 2248630589 : for (unsigned int qp = 0; qp < n_points; qp++)
1750 1823172886 : coordTransformFactor(_subproblem, sub_id, q_points[qp], coord[qp]);
1751 425457703 : }
1752 :
1753 : void
1754 394339397 : Assembly::computeCurrentElemVolume()
1755 : {
1756 394339397 : if (_current_elem_volume_computed)
1757 0 : return;
1758 :
1759 394339397 : setCoordinateTransformation(
1760 394339397 : _current_qrule, _current_q_points, _coord, _current_elem->subdomain_id());
1761 394339397 : if (_calculate_ad_coord)
1762 12751991 : setCoordinateTransformation(
1763 12751991 : _current_qrule, _ad_q_points, _ad_coord, _current_elem->subdomain_id());
1764 :
1765 394339397 : _current_elem_volume = 0.;
1766 2113796998 : for (unsigned int qp = 0; qp < _current_qrule->n_points(); qp++)
1767 1719457601 : _current_elem_volume += _current_JxW[qp] * _coord[qp];
1768 :
1769 394339397 : _current_elem_volume_computed = true;
1770 : }
1771 :
1772 : void
1773 11953620 : Assembly::computeCurrentFaceVolume()
1774 : {
1775 11953620 : if (_current_side_volume_computed)
1776 0 : return;
1777 :
1778 11953620 : setCoordinateTransformation(
1779 11953620 : _current_qrule_face, _current_q_points_face, _coord, _current_elem->subdomain_id());
1780 11953620 : if (_calculate_ad_coord)
1781 3637669 : setCoordinateTransformation(
1782 3637669 : _current_qrule_face, _ad_q_points_face, _ad_coord, _current_elem->subdomain_id());
1783 :
1784 11953620 : _current_side_volume = 0.;
1785 38984236 : for (unsigned int qp = 0; qp < _current_qrule_face->n_points(); qp++)
1786 27030616 : _current_side_volume += _current_JxW_face[qp] * _coord[qp];
1787 :
1788 11953620 : _current_side_volume_computed = true;
1789 : }
1790 :
1791 : void
1792 351016 : Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points)
1793 : {
1794 351016 : _current_elem = elem;
1795 351016 : _current_neighbor_elem = nullptr;
1796 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1797 : "current subdomain has been set incorrectly");
1798 351016 : _current_elem_volume_computed = false;
1799 :
1800 351016 : FEMap::inverse_map(elem->dim(), elem, physical_points, _temp_reference_points);
1801 :
1802 351016 : reinit(elem, _temp_reference_points);
1803 :
1804 : // Save off the physical points
1805 351016 : _current_physical_points = physical_points;
1806 351016 : }
1807 :
1808 : void
1809 394062908 : Assembly::setVolumeQRule(const Elem * const elem)
1810 : {
1811 394062908 : unsigned int elem_dimension = elem->dim();
1812 394062908 : _current_qrule_volume = qrules(elem_dimension).vol.get();
1813 : // Make sure the qrule is the right one
1814 394062908 : if (_current_qrule != _current_qrule_volume)
1815 194730 : setVolumeQRule(_current_qrule_volume, elem_dimension);
1816 394062908 : }
1817 :
1818 : void
1819 393988409 : Assembly::reinit(const Elem * elem)
1820 : {
1821 393988409 : _current_elem = elem;
1822 393988409 : _current_neighbor_elem = nullptr;
1823 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1824 : "current subdomain has been set incorrectly");
1825 393988409 : _current_elem_volume_computed = false;
1826 393988409 : if (_sc)
1827 883838 : _sc->set_current_elem(*elem);
1828 :
1829 393988409 : setVolumeQRule(elem);
1830 393988409 : reinitFE(elem);
1831 :
1832 393988381 : computeCurrentElemVolume();
1833 393988381 : }
1834 :
1835 : void
1836 351016 : Assembly::reinit(const Elem * elem, const std::vector<Point> & reference_points)
1837 : {
1838 351016 : _current_elem = elem;
1839 351016 : _current_neighbor_elem = nullptr;
1840 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1841 : "current subdomain has been set incorrectly");
1842 351016 : _current_elem_volume_computed = false;
1843 :
1844 351016 : unsigned int elem_dimension = _current_elem->dim();
1845 :
1846 351016 : _current_qrule_arbitrary = qrules(elem_dimension).arbitrary_vol.get();
1847 :
1848 : // Make sure the qrule is the right one
1849 351016 : if (_current_qrule != _current_qrule_arbitrary)
1850 30911 : setVolumeQRule(_current_qrule_arbitrary, elem_dimension);
1851 :
1852 351016 : _current_qrule_arbitrary->setPoints(reference_points);
1853 :
1854 351016 : reinitFE(elem);
1855 :
1856 351016 : computeCurrentElemVolume();
1857 351016 : }
1858 :
1859 : void
1860 24516206 : Assembly::reinitFVFace(const FaceInfo & fi)
1861 : {
1862 24516206 : _current_elem = &fi.elem();
1863 24516206 : _current_neighbor_elem = fi.neighborPtr();
1864 24516206 : _current_side = fi.elemSideID();
1865 24516206 : _current_neighbor_side = fi.neighborSideID();
1866 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1867 : "current subdomain has been set incorrectly");
1868 :
1869 24516206 : _current_elem_volume_computed = false;
1870 24516206 : _current_side_volume_computed = false;
1871 :
1872 24516206 : prepareResidual();
1873 24516206 : prepareNeighbor();
1874 24516206 : prepareJacobianBlock();
1875 :
1876 24516206 : unsigned int dim = _current_elem->dim();
1877 24516206 : if (_current_qrule_face != qrules(dim).fv_face.get())
1878 : {
1879 5104 : setFaceQRule(qrules(dim).fv_face.get(), dim);
1880 : // The order of the element that is used for initing here doesn't matter since this will just
1881 : // be used for constant monomials (which only need a single integration point)
1882 5104 : if (dim == 3)
1883 12 : _current_qrule_face->init(QUAD4, /* p_level = */ 0, /* simple_type_only = */ true);
1884 : else
1885 5092 : _current_qrule_face->init(EDGE2, /* p_level = */ 0, /* simple_type_only = */ true);
1886 : }
1887 :
1888 24516206 : _current_side_elem = &_current_side_elem_builder(*_current_elem, _current_side);
1889 :
1890 : mooseAssert(_current_qrule_face->n_points() == 1,
1891 : "Our finite volume quadrature rule should always yield a single point");
1892 :
1893 : // We've initialized the reference points. Now we need to compute the physical location of the
1894 : // quadrature points. We do not do any FE initialization so we cannot simply copy over FE
1895 : // results like we do in reinitFEFace. Instead we handle the computation of the physical
1896 : // locations manually
1897 24516206 : _current_q_points_face.resize(1);
1898 24516206 : const auto & ref_points = _current_qrule_face->get_points();
1899 24516206 : const auto & ref_point = ref_points[0];
1900 24516206 : auto physical_point = FEMap::map(_current_side_elem->dim(), _current_side_elem, ref_point);
1901 24516206 : _current_q_points_face[0] = physical_point;
1902 :
1903 24516206 : if (_current_neighbor_elem)
1904 : {
1905 : mooseAssert(_current_neighbor_subdomain_id == _current_neighbor_elem->subdomain_id(),
1906 : "current neighbor subdomain has been set incorrectly");
1907 : // Now handle the neighbor qrule/qpoints
1908 : ArbitraryQuadrature * const neighbor_rule =
1909 22694094 : qrules(_current_neighbor_elem->dim(), _current_neighbor_subdomain_id).neighbor.get();
1910 : // Here we are setting a reference point that is correct for the neighbor *side* element. It
1911 : // would be wrong if this reference point is used for a volumetric FE reinit with the neighbor
1912 22694094 : neighbor_rule->setPoints(ref_points);
1913 22694094 : setNeighborQRule(neighbor_rule, _current_neighbor_elem->dim());
1914 22694094 : _current_q_points_face_neighbor.resize(1);
1915 22694094 : _current_q_points_face_neighbor[0] = std::move(physical_point);
1916 : }
1917 24516206 : }
1918 :
1919 : QBase *
1920 11953620 : Assembly::qruleFace(const Elem * elem, unsigned int side)
1921 : {
1922 24117378 : return qruleFaceHelper<QBase>(elem, side, [](QRules & q) { return q.face.get(); });
1923 : }
1924 :
1925 : ArbitraryQuadrature *
1926 283554 : Assembly::qruleArbitraryFace(const Elem * elem, unsigned int side)
1927 : {
1928 567108 : return qruleFaceHelper<ArbitraryQuadrature>(
1929 850662 : elem, side, [](QRules & q) { return q.arbitrary_face.get(); });
1930 : }
1931 :
1932 : void
1933 11953620 : Assembly::setFaceQRule(const Elem * const elem, const unsigned int side)
1934 : {
1935 11953620 : const auto elem_dimension = elem->dim();
1936 : //// Make sure the qrule is the right one
1937 11953620 : auto rule = qruleFace(elem, side);
1938 11953620 : if (_current_qrule_face != rule)
1939 11697 : setFaceQRule(rule, elem_dimension);
1940 11953620 : }
1941 :
1942 : void
1943 11953620 : Assembly::reinit(const Elem * const elem, const unsigned int side)
1944 : {
1945 11953620 : _current_elem = elem;
1946 11953620 : _current_neighbor_elem = nullptr;
1947 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1948 : "current subdomain has been set incorrectly");
1949 11953620 : _current_side = side;
1950 11953620 : _current_elem_volume_computed = false;
1951 11953620 : _current_side_volume_computed = false;
1952 :
1953 11953620 : _current_side_elem = &_current_side_elem_builder(*elem, side);
1954 :
1955 11953620 : setFaceQRule(elem, side);
1956 11953620 : reinitFEFace(elem, side);
1957 :
1958 11953620 : computeCurrentFaceVolume();
1959 11953620 : }
1960 :
1961 : void
1962 0 : Assembly::reinit(const Elem * elem, unsigned int side, const std::vector<Point> & reference_points)
1963 : {
1964 0 : _current_elem = elem;
1965 0 : _current_neighbor_elem = nullptr;
1966 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1967 : "current subdomain has been set incorrectly");
1968 0 : _current_side = side;
1969 0 : _current_elem_volume_computed = false;
1970 0 : _current_side_volume_computed = false;
1971 :
1972 0 : unsigned int elem_dimension = _current_elem->dim();
1973 :
1974 0 : _current_qrule_arbitrary_face = qruleArbitraryFace(elem, side);
1975 :
1976 : // Make sure the qrule is the right one
1977 0 : if (_current_qrule_face != _current_qrule_arbitrary_face)
1978 0 : setFaceQRule(_current_qrule_arbitrary_face, elem_dimension);
1979 :
1980 0 : _current_qrule_arbitrary->setPoints(reference_points);
1981 :
1982 0 : _current_side_elem = &_current_side_elem_builder(*elem, side);
1983 :
1984 0 : reinitFEFace(elem, side);
1985 :
1986 0 : computeCurrentFaceVolume();
1987 0 : }
1988 :
1989 : void
1990 106537100 : Assembly::reinit(const Node * node)
1991 : {
1992 106537100 : _current_node = node;
1993 106537100 : _current_neighbor_node = NULL;
1994 106537100 : }
1995 :
1996 : void
1997 3876795 : Assembly::reinitElemAndNeighbor(const Elem * elem,
1998 : unsigned int side,
1999 : const Elem * neighbor,
2000 : unsigned int neighbor_side,
2001 : const std::vector<Point> * neighbor_reference_points)
2002 : {
2003 3876795 : _current_neighbor_side = neighbor_side;
2004 :
2005 3876795 : reinit(elem, side);
2006 :
2007 3876795 : unsigned int neighbor_dim = neighbor->dim();
2008 :
2009 3876795 : if (neighbor_reference_points)
2010 64908 : _current_neighbor_ref_points = *neighbor_reference_points;
2011 : else
2012 3811887 : FEMap::inverse_map(
2013 7623774 : neighbor_dim, neighbor, _current_q_points_face.stdVector(), _current_neighbor_ref_points);
2014 :
2015 3876795 : _current_neighbor_side_elem = &_current_neighbor_side_elem_builder(*neighbor, neighbor_side);
2016 :
2017 3876795 : reinitFEFaceNeighbor(neighbor, _current_neighbor_ref_points);
2018 3876795 : reinitNeighbor(neighbor, _current_neighbor_ref_points);
2019 3876795 : }
2020 :
2021 : void
2022 283554 : Assembly::reinitElemFaceRef(const Elem * elem,
2023 : unsigned int elem_side,
2024 : Real tolerance,
2025 : const std::vector<Point> * const pts,
2026 : const std::vector<Real> * const weights)
2027 : {
2028 283554 : _current_elem = elem;
2029 :
2030 283554 : unsigned int elem_dim = elem->dim();
2031 :
2032 : // Attach the quadrature rules
2033 283554 : if (pts)
2034 : {
2035 283554 : auto face_rule = qruleArbitraryFace(elem, elem_side);
2036 283554 : face_rule->setPoints(*pts);
2037 283554 : setFaceQRule(face_rule, elem_dim);
2038 : }
2039 : else
2040 : {
2041 0 : auto rule = qruleFace(elem, elem_side);
2042 0 : if (_current_qrule_face != rule)
2043 0 : setFaceQRule(rule, elem_dim);
2044 : }
2045 :
2046 : // reinit face
2047 671376 : for (const auto & it : _fe_face[elem_dim])
2048 : {
2049 387822 : FEBase & fe_face = *it.second;
2050 387822 : FEType fe_type = it.first;
2051 387822 : FEShapeData & fesd = *_fe_shape_data_face[fe_type];
2052 :
2053 387822 : fe_face.reinit(elem, elem_side, tolerance, pts, weights);
2054 :
2055 387822 : _current_fe_face[fe_type] = &fe_face;
2056 :
2057 387822 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face.get_phi()));
2058 387822 : fesd._grad_phi.shallowCopy(
2059 387822 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_face.get_dphi()));
2060 387822 : if (_need_second_derivative_neighbor.count(fe_type))
2061 0 : fesd._second_phi.shallowCopy(
2062 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
2063 : }
2064 283554 : for (const auto & it : _vector_fe_face[elem_dim])
2065 : {
2066 0 : FEVectorBase & fe_face = *it.second;
2067 0 : const FEType & fe_type = it.first;
2068 :
2069 0 : _current_vector_fe_face[fe_type] = &fe_face;
2070 :
2071 0 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face[fe_type];
2072 :
2073 0 : fe_face.reinit(elem, elem_side, tolerance, pts, weights);
2074 :
2075 0 : fesd._phi.shallowCopy(
2076 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_phi()));
2077 0 : fesd._grad_phi.shallowCopy(
2078 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_dphi()));
2079 0 : if (_need_second_derivative.count(fe_type))
2080 0 : fesd._second_phi.shallowCopy(
2081 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
2082 0 : if (_need_curl.count(fe_type))
2083 0 : fesd._curl_phi.shallowCopy(
2084 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
2085 0 : if (_need_face_div.count(fe_type))
2086 0 : fesd._div_phi.shallowCopy(
2087 0 : const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
2088 : }
2089 283554 : if (!_unique_fe_face_helper.empty())
2090 : {
2091 : mooseAssert(elem_dim < _unique_fe_face_helper.size(), "We should be in bounds here");
2092 0 : _unique_fe_face_helper[elem_dim]->reinit(elem, elem_side, tolerance, pts, weights);
2093 : }
2094 :
2095 : // During that last loop the helper objects will have been reinitialized
2096 283554 : _current_q_points_face.shallowCopy(
2097 283554 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_xyz()));
2098 283554 : _current_normals.shallowCopy(
2099 283554 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_normals()));
2100 283554 : _current_tangents.shallowCopy(const_cast<std::vector<std::vector<Point>> &>(
2101 283554 : _holder_fe_face_helper[elem_dim]->get_tangents()));
2102 : // Note that if the user did pass in points and not weights to this method, JxW will be garbage
2103 : // and should not be used
2104 283554 : _current_JxW_face.shallowCopy(
2105 283554 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_JxW()));
2106 283554 : if (_calculate_curvatures)
2107 0 : _curvatures.shallowCopy(
2108 0 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_curvatures()));
2109 :
2110 283554 : computeADFace(*elem, elem_side);
2111 283554 : }
2112 :
2113 : void
2114 12237174 : Assembly::computeADFace(const Elem & elem, const unsigned int side)
2115 : {
2116 12237174 : const auto dim = elem.dim();
2117 :
2118 12237174 : if (_subproblem.haveADObjects())
2119 : {
2120 4089390 : auto n_qp = _current_qrule_face->n_points();
2121 4089390 : resizeADMappingObjects(n_qp, dim);
2122 4089390 : _ad_normals.resize(n_qp);
2123 4089390 : _ad_JxW_face.resize(n_qp);
2124 4089390 : if (_calculate_face_xyz)
2125 3655617 : _ad_q_points_face.resize(n_qp);
2126 4089390 : if (_calculate_curvatures)
2127 416 : _ad_curvatures.resize(n_qp);
2128 :
2129 4089390 : if (_displaced)
2130 : {
2131 53788 : const auto & qw = _current_qrule_face->get_weights();
2132 53788 : computeFaceMap(elem, side, qw);
2133 53788 : const std::vector<Real> dummy_qw(n_qp, 1.);
2134 :
2135 161116 : for (unsigned int qp = 0; qp != n_qp; qp++)
2136 107328 : computeSinglePointMapAD(&elem, dummy_qw, qp, _holder_fe_face_helper[dim]);
2137 53788 : }
2138 : else
2139 : {
2140 12836956 : for (unsigned qp = 0; qp < n_qp; ++qp)
2141 : {
2142 8801354 : _ad_JxW_face[qp] = _current_JxW_face[qp];
2143 8801354 : _ad_normals[qp] = _current_normals[qp];
2144 : }
2145 4035602 : if (_calculate_face_xyz)
2146 11156515 : for (unsigned qp = 0; qp < n_qp; ++qp)
2147 7534314 : _ad_q_points_face[qp] = _current_q_points_face[qp];
2148 4035602 : if (_calculate_curvatures)
2149 832 : for (unsigned qp = 0; qp < n_qp; ++qp)
2150 416 : _ad_curvatures[qp] = _curvatures[qp];
2151 : }
2152 :
2153 13010932 : for (const auto & it : _fe_face[dim])
2154 : {
2155 8921542 : FEBase & fe = *it.second;
2156 8921542 : auto fe_type = it.first;
2157 8921542 : auto num_shapes = FEInterface::n_shape_functions(fe_type, &elem);
2158 8921542 : auto & grad_phi = _ad_grad_phi_data_face[fe_type];
2159 :
2160 8921542 : grad_phi.resize(num_shapes);
2161 66260802 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2162 57339260 : grad_phi[i].resize(n_qp);
2163 :
2164 8921542 : const auto & regular_grad_phi = _fe_shape_data_face[fe_type]->_grad_phi;
2165 :
2166 8921542 : if (_displaced)
2167 131456 : computeGradPhiAD(&elem, n_qp, grad_phi, &fe);
2168 : else
2169 65442702 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2170 187861180 : for (unsigned qp = 0; qp < n_qp; ++qp)
2171 131208564 : grad_phi[i][qp] = regular_grad_phi[i][qp];
2172 : }
2173 4588203 : for (const auto & it : _vector_fe_face[dim])
2174 : {
2175 498813 : FEVectorBase & fe = *it.second;
2176 498813 : auto fe_type = it.first;
2177 498813 : auto num_shapes = FEInterface::n_shape_functions(fe_type, &elem);
2178 498813 : auto & grad_phi = _ad_vector_grad_phi_data_face[fe_type];
2179 :
2180 498813 : grad_phi.resize(num_shapes);
2181 2621137 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2182 2122324 : grad_phi[i].resize(n_qp);
2183 :
2184 498813 : const auto & regular_grad_phi = _vector_fe_shape_data_face[fe_type]->_grad_phi;
2185 :
2186 498813 : if (_displaced)
2187 0 : computeGradPhiAD(&elem, n_qp, grad_phi, &fe);
2188 : else
2189 2621137 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2190 6480120 : for (unsigned qp = 0; qp < n_qp; ++qp)
2191 4357796 : grad_phi[i][qp] = regular_grad_phi[i][qp];
2192 : }
2193 : }
2194 12237174 : }
2195 :
2196 : void
2197 283554 : Assembly::reinitNeighborFaceRef(const Elem * neighbor,
2198 : unsigned int neighbor_side,
2199 : Real tolerance,
2200 : const std::vector<Point> * const pts,
2201 : const std::vector<Real> * const weights)
2202 : {
2203 283554 : _current_neighbor_elem = neighbor;
2204 :
2205 283554 : unsigned int neighbor_dim = neighbor->dim();
2206 :
2207 : ArbitraryQuadrature * neighbor_rule =
2208 283554 : qrules(neighbor_dim, neighbor->subdomain_id()).neighbor.get();
2209 283554 : neighbor_rule->setPoints(*pts);
2210 :
2211 : // Attach this quadrature rule to all the _fe_face_neighbor FE objects. This
2212 : // has to have garbage quadrature weights but that's ok because we never
2213 : // actually use the JxW coming from these FE reinit'd objects, e.g. we use the
2214 : // JxW coming from the element face reinit for DGKernels or we use the JxW
2215 : // coming from reinit of the mortar segment element in the case of mortar
2216 283554 : setNeighborQRule(neighbor_rule, neighbor_dim);
2217 :
2218 : // reinit neighbor face
2219 671376 : for (const auto & it : _fe_face_neighbor[neighbor_dim])
2220 : {
2221 387822 : FEBase & fe_face_neighbor = *it.second;
2222 387822 : FEType fe_type = it.first;
2223 387822 : FEShapeData & fesd = *_fe_shape_data_face_neighbor[fe_type];
2224 :
2225 387822 : fe_face_neighbor.reinit(neighbor, neighbor_side, tolerance, pts, weights);
2226 :
2227 387822 : _current_fe_face_neighbor[fe_type] = &fe_face_neighbor;
2228 :
2229 387822 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_phi()));
2230 387822 : fesd._grad_phi.shallowCopy(
2231 387822 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor.get_dphi()));
2232 387822 : if (_need_second_derivative_neighbor.count(fe_type))
2233 0 : fesd._second_phi.shallowCopy(
2234 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
2235 : }
2236 283554 : for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
2237 : {
2238 0 : FEVectorBase & fe_face_neighbor = *it.second;
2239 0 : const FEType & fe_type = it.first;
2240 :
2241 0 : _current_vector_fe_face_neighbor[fe_type] = &fe_face_neighbor;
2242 :
2243 0 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face_neighbor[fe_type];
2244 :
2245 0 : fe_face_neighbor.reinit(neighbor, neighbor_side, tolerance, pts, weights);
2246 :
2247 0 : fesd._phi.shallowCopy(
2248 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor.get_phi()));
2249 0 : fesd._grad_phi.shallowCopy(
2250 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_dphi()));
2251 0 : if (_need_second_derivative.count(fe_type))
2252 0 : fesd._second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
2253 0 : fe_face_neighbor.get_d2phi()));
2254 0 : if (_need_curl.count(fe_type))
2255 0 : fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
2256 0 : fe_face_neighbor.get_curl_phi()));
2257 0 : if (_need_face_neighbor_div.count(fe_type))
2258 0 : fesd._div_phi.shallowCopy(
2259 0 : const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
2260 : }
2261 283554 : if (!_unique_fe_face_neighbor_helper.empty())
2262 : {
2263 : mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
2264 : "We should be in bounds here");
2265 0 : _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(
2266 : neighbor, neighbor_side, tolerance, pts, weights);
2267 : }
2268 : // During that last loop the helper objects will have been reinitialized as well
2269 : // We need to dig out the q_points from it
2270 283554 : _current_q_points_face_neighbor.shallowCopy(
2271 283554 : const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
2272 283554 : }
2273 :
2274 : void
2275 4057 : Assembly::reinitDual(const Elem * elem,
2276 : const std::vector<Point> & pts,
2277 : const std::vector<Real> & JxW)
2278 : {
2279 4057 : const unsigned int elem_dim = elem->dim();
2280 : mooseAssert(elem_dim == _mesh_dimension - 1,
2281 : "Dual shape functions should only be computed on lower dimensional face elements");
2282 :
2283 8390 : for (const auto & it : _fe_lower[elem_dim])
2284 : {
2285 4333 : FEBase & fe_lower = *it.second;
2286 : // We use customized quadrature rule for integration along the mortar segment elements
2287 4333 : fe_lower.set_calculate_default_dual_coeff(false);
2288 4333 : fe_lower.reinit_dual_shape_coeffs(elem, pts, JxW);
2289 : }
2290 4057 : }
2291 :
2292 : void
2293 299795 : Assembly::reinitLowerDElem(const Elem * elem,
2294 : const std::vector<Point> * const pts,
2295 : const std::vector<Real> * const weights)
2296 : {
2297 299795 : _current_lower_d_elem = elem;
2298 :
2299 299795 : const unsigned int elem_dim = elem->dim();
2300 : mooseAssert(elem_dim < _mesh_dimension,
2301 : "The lower dimensional element should truly be a lower dimensional element");
2302 :
2303 299795 : if (pts)
2304 : {
2305 : // Lower rule matches the face rule for the higher dimensional element
2306 263086 : ArbitraryQuadrature * lower_rule = qrules(elem_dim + 1).arbitrary_face.get();
2307 :
2308 : // This also sets the quadrature weights to unity
2309 263086 : lower_rule->setPoints(*pts);
2310 :
2311 263086 : if (weights)
2312 0 : lower_rule->setWeights(*weights);
2313 :
2314 263086 : setLowerQRule(lower_rule, elem_dim);
2315 : }
2316 36709 : else if (_current_qrule_lower != qrules(elem_dim + 1).face.get())
2317 429 : setLowerQRule(qrules(elem_dim + 1).face.get(), elem_dim);
2318 :
2319 740084 : for (const auto & it : _fe_lower[elem_dim])
2320 : {
2321 440289 : FEBase & fe_lower = *it.second;
2322 440289 : FEType fe_type = it.first;
2323 :
2324 440289 : fe_lower.reinit(elem);
2325 :
2326 440289 : if (FEShapeData * fesd = _fe_shape_data_lower[fe_type].get())
2327 : {
2328 440289 : fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower.get_phi()));
2329 440289 : fesd->_grad_phi.shallowCopy(
2330 440289 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower.get_dphi()));
2331 440289 : if (_need_second_derivative_neighbor.count(fe_type))
2332 0 : fesd->_second_phi.shallowCopy(
2333 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_d2phi()));
2334 : }
2335 :
2336 : // Dual shape functions need to be computed after primal basis being initialized
2337 440289 : if (FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type].get())
2338 : {
2339 12142 : fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower.get_dual_phi()));
2340 12142 : fesd->_grad_phi.shallowCopy(
2341 12142 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower.get_dual_dphi()));
2342 12142 : if (_need_second_derivative_neighbor.count(fe_type))
2343 0 : fesd->_second_phi.shallowCopy(
2344 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_dual_d2phi()));
2345 : }
2346 : }
2347 299795 : if (!_unique_fe_lower_helper.empty())
2348 : {
2349 : mooseAssert(elem_dim < _unique_fe_lower_helper.size(), "We should be in bounds here");
2350 0 : _unique_fe_lower_helper[elem_dim]->reinit(elem);
2351 : }
2352 :
2353 299795 : if (!_need_lower_d_elem_volume)
2354 199345 : return;
2355 :
2356 100450 : if (pts && !weights)
2357 : {
2358 : // We only have dummy weights so the JxWs computed during our FE reinits are meaningless and
2359 : // we cannot use them
2360 :
2361 97282 : if (_subproblem.getCoordSystem(elem->subdomain_id()) == Moose::CoordinateSystemType::COORD_XYZ)
2362 : // We are in a Cartesian coordinate system and we can just use the element volume method
2363 : // which has fast computation for certain element types
2364 97282 : _current_lower_d_elem_volume = elem->volume();
2365 : else
2366 : // We manually compute the volume taking the curvilinear coordinate transformations into
2367 : // account
2368 0 : _current_lower_d_elem_volume = elementVolume(elem);
2369 : }
2370 : else
2371 : {
2372 : // During that last loop the helper objects will have been reinitialized as well
2373 3168 : FEBase & helper_fe = *_holder_fe_lower_helper[elem_dim];
2374 3168 : const auto & physical_q_points = helper_fe.get_xyz();
2375 3168 : const auto & JxW = helper_fe.get_JxW();
2376 3168 : MooseArray<Real> coord;
2377 3168 : setCoordinateTransformation(
2378 3168 : _current_qrule_lower, physical_q_points, coord, elem->subdomain_id());
2379 3168 : _current_lower_d_elem_volume = 0;
2380 46944 : for (const auto qp : make_range(_current_qrule_lower->n_points()))
2381 43776 : _current_lower_d_elem_volume += JxW[qp] * coord[qp];
2382 3168 : }
2383 : }
2384 :
2385 : void
2386 263086 : Assembly::reinitNeighborLowerDElem(const Elem * elem)
2387 : {
2388 : mooseAssert(elem->dim() < _mesh_dimension,
2389 : "You should be calling reinitNeighborLowerDElem on a lower dimensional element");
2390 :
2391 263086 : _current_neighbor_lower_d_elem = elem;
2392 :
2393 263086 : if (!_need_neighbor_lower_d_elem_volume)
2394 165804 : return;
2395 :
2396 97282 : if (_subproblem.getCoordSystem(elem->subdomain_id()) == Moose::CoordinateSystemType::COORD_XYZ)
2397 : // We are in a Cartesian coordinate system and we can just use the element volume method which
2398 : // has fast computation for certain element types
2399 97282 : _current_neighbor_lower_d_elem_volume = elem->volume();
2400 : else
2401 : // We manually compute the volume taking the curvilinear coordinate transformations into
2402 : // account
2403 0 : _current_neighbor_lower_d_elem_volume = elementVolume(elem);
2404 : }
2405 :
2406 : void
2407 526172 : Assembly::reinitMortarElem(const Elem * elem)
2408 : {
2409 : mooseAssert(elem->dim() == _mesh_dimension - 1,
2410 : "You should be calling reinitMortarElem on a lower dimensional element");
2411 :
2412 526172 : _fe_msm->reinit(elem);
2413 526172 : _msm_elem = elem;
2414 :
2415 526172 : MooseArray<Point> array_q_points;
2416 526172 : array_q_points.shallowCopy(const_cast<std::vector<Point> &>(_fe_msm->get_xyz()));
2417 526172 : setCoordinateTransformation(_qrule_msm, array_q_points, _coord_msm, elem->subdomain_id());
2418 526172 : }
2419 :
2420 : void
2421 147316 : Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
2422 : unsigned int neighbor_side,
2423 : const std::vector<Point> & physical_points)
2424 : {
2425 147316 : unsigned int neighbor_dim = neighbor->dim();
2426 147316 : FEMap::inverse_map(neighbor_dim, neighbor, physical_points, _current_neighbor_ref_points);
2427 :
2428 147316 : if (_need_JxW_neighbor)
2429 : {
2430 : mooseAssert(
2431 : physical_points.size() == 1,
2432 : "If reinitializing with more than one point, then I am dubious of your use case. Perhaps "
2433 : "you are performing a DG type method and you are reinitializing using points from the "
2434 : "element face. In such a case your neighbor JxW must have its index order 'match' the "
2435 : "element JxW index order, e.g. imagining a vertical 1D face with two quadrature points, "
2436 : "if "
2437 : "index 0 for elem JxW corresponds to the 'top' quadrature point, then index 0 for "
2438 : "neighbor "
2439 : "JxW must also correspond to the 'top' quadrature point. And libMesh/MOOSE has no way to "
2440 : "guarantee that with multiple quadrature points.");
2441 :
2442 96052 : _current_neighbor_side_elem = &_current_neighbor_side_elem_builder(*neighbor, neighbor_side);
2443 :
2444 : // With a single point our size-1 JxW should just be the element volume
2445 96052 : _current_JxW_neighbor.resize(1);
2446 96052 : _current_JxW_neighbor[0] = _current_neighbor_side_elem->volume();
2447 : }
2448 :
2449 147316 : reinitFEFaceNeighbor(neighbor, _current_neighbor_ref_points);
2450 147316 : reinitNeighbor(neighbor, _current_neighbor_ref_points);
2451 :
2452 : // Save off the physical points
2453 147316 : _current_physical_points = physical_points;
2454 147316 : }
2455 :
2456 : void
2457 19880 : Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
2458 : const std::vector<Point> & physical_points)
2459 : {
2460 19880 : unsigned int neighbor_dim = neighbor->dim();
2461 19880 : FEMap::inverse_map(neighbor_dim, neighbor, physical_points, _current_neighbor_ref_points);
2462 :
2463 19880 : reinitFENeighbor(neighbor, _current_neighbor_ref_points);
2464 19880 : reinitNeighbor(neighbor, _current_neighbor_ref_points);
2465 : // Save off the physical points
2466 19880 : _current_physical_points = physical_points;
2467 19880 : }
2468 :
2469 : void
2470 63135 : Assembly::init(const CouplingMatrix * cm)
2471 : {
2472 63135 : _cm = cm;
2473 :
2474 63135 : unsigned int n_vars = _sys.nVariables();
2475 :
2476 63135 : _cm_ss_entry.clear();
2477 63135 : _cm_sf_entry.clear();
2478 63135 : _cm_fs_entry.clear();
2479 63135 : _cm_ff_entry.clear();
2480 :
2481 63135 : auto & vars = _sys.getVariables(_tid);
2482 :
2483 63135 : _block_diagonal_matrix = true;
2484 127183 : for (auto & ivar : vars)
2485 : {
2486 64048 : auto i = ivar->number();
2487 64048 : if (i >= _component_block_diagonal.size())
2488 64048 : _component_block_diagonal.resize(i + 1, true);
2489 :
2490 64048 : auto ivar_start = _cm_ff_entry.size();
2491 129709 : for (unsigned int k = 0; k < ivar->count(); ++k)
2492 : {
2493 65661 : unsigned int iv = i + k;
2494 153922 : for (const auto & j : ConstCouplingRow(iv, *_cm))
2495 : {
2496 88261 : if (_sys.isScalarVariable(j))
2497 : {
2498 1090 : auto & jvar = _sys.getScalarVariable(_tid, j);
2499 1090 : _cm_fs_entry.push_back(std::make_pair(ivar, &jvar));
2500 1090 : _block_diagonal_matrix = false;
2501 : }
2502 : else
2503 : {
2504 87171 : auto & jvar = _sys.getVariable(_tid, j);
2505 87171 : auto pair = std::make_pair(ivar, &jvar);
2506 87171 : auto c = ivar_start;
2507 : // check if the pair has been pushed or not
2508 87171 : bool has_pair = false;
2509 114105 : for (; c < _cm_ff_entry.size(); ++c)
2510 32627 : if (_cm_ff_entry[c] == pair)
2511 : {
2512 5693 : has_pair = true;
2513 5693 : break;
2514 : }
2515 87171 : if (!has_pair)
2516 81478 : _cm_ff_entry.push_back(pair);
2517 : // only set having diagonal matrix to false when ivar and jvar numbers are different
2518 : // Note: for array variables, since we save the entire local Jacobian of all components,
2519 : // even there are couplings among components of the same array variable, we still
2520 : // do not set the flag to false.
2521 87171 : if (i != jvar.number())
2522 19896 : _block_diagonal_matrix = false;
2523 67275 : else if (iv != j)
2524 1614 : _component_block_diagonal[i] = false;
2525 : }
2526 : }
2527 : }
2528 : }
2529 :
2530 63135 : auto & scalar_vars = _sys.getScalarVariables(_tid);
2531 :
2532 64561 : for (auto & ivar : scalar_vars)
2533 : {
2534 1426 : auto i = ivar->number();
2535 1426 : if (i >= _component_block_diagonal.size())
2536 1352 : _component_block_diagonal.resize(i + 1, true);
2537 :
2538 4078 : for (const auto & j : ConstCouplingRow(i, *_cm))
2539 2652 : if (_sys.isScalarVariable(j))
2540 : {
2541 1562 : auto & jvar = _sys.getScalarVariable(_tid, j);
2542 1562 : _cm_ss_entry.push_back(std::make_pair(ivar, &jvar));
2543 : }
2544 : else
2545 : {
2546 1090 : auto & jvar = _sys.getVariable(_tid, j);
2547 1090 : _cm_sf_entry.push_back(std::make_pair(ivar, &jvar));
2548 : }
2549 : }
2550 :
2551 63135 : if (_block_diagonal_matrix && scalar_vars.size() != 0)
2552 443 : _block_diagonal_matrix = false;
2553 :
2554 63135 : auto num_vector_tags = _residual_vector_tags.size();
2555 :
2556 63135 : _sub_Re.resize(num_vector_tags);
2557 63135 : _sub_Rn.resize(num_vector_tags);
2558 63135 : _sub_Rl.resize(num_vector_tags);
2559 224372 : for (MooseIndex(_sub_Re) i = 0; i < _sub_Re.size(); i++)
2560 : {
2561 161237 : _sub_Re[i].resize(n_vars);
2562 161237 : _sub_Rn[i].resize(n_vars);
2563 161237 : _sub_Rl[i].resize(n_vars);
2564 : }
2565 :
2566 63135 : _cached_residual_values.resize(num_vector_tags);
2567 63135 : _cached_residual_rows.resize(num_vector_tags);
2568 :
2569 63135 : auto num_matrix_tags = _subproblem.numMatrixTags();
2570 :
2571 63135 : _cached_jacobian_values.resize(num_matrix_tags);
2572 63135 : _cached_jacobian_rows.resize(num_matrix_tags);
2573 63135 : _cached_jacobian_cols.resize(num_matrix_tags);
2574 :
2575 : // Element matrices
2576 63135 : _sub_Kee.resize(num_matrix_tags);
2577 63135 : _sub_Keg.resize(num_matrix_tags);
2578 63135 : _sub_Ken.resize(num_matrix_tags);
2579 63135 : _sub_Kne.resize(num_matrix_tags);
2580 63135 : _sub_Knn.resize(num_matrix_tags);
2581 63135 : _sub_Kll.resize(num_matrix_tags);
2582 63135 : _sub_Kle.resize(num_matrix_tags);
2583 63135 : _sub_Kln.resize(num_matrix_tags);
2584 63135 : _sub_Kel.resize(num_matrix_tags);
2585 63135 : _sub_Knl.resize(num_matrix_tags);
2586 :
2587 63135 : _jacobian_block_used.resize(num_matrix_tags);
2588 63135 : _jacobian_block_neighbor_used.resize(num_matrix_tags);
2589 63135 : _jacobian_block_lower_used.resize(num_matrix_tags);
2590 63135 : _jacobian_block_nonlocal_used.resize(num_matrix_tags);
2591 :
2592 191692 : for (MooseIndex(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
2593 : {
2594 128557 : _sub_Keg[tag].resize(n_vars);
2595 128557 : _sub_Ken[tag].resize(n_vars);
2596 128557 : _sub_Kne[tag].resize(n_vars);
2597 128557 : _sub_Knn[tag].resize(n_vars);
2598 128557 : _sub_Kee[tag].resize(n_vars);
2599 128557 : _sub_Kll[tag].resize(n_vars);
2600 128557 : _sub_Kle[tag].resize(n_vars);
2601 128557 : _sub_Kln[tag].resize(n_vars);
2602 128557 : _sub_Kel[tag].resize(n_vars);
2603 128557 : _sub_Knl[tag].resize(n_vars);
2604 :
2605 128557 : _jacobian_block_used[tag].resize(n_vars);
2606 128557 : _jacobian_block_neighbor_used[tag].resize(n_vars);
2607 128557 : _jacobian_block_lower_used[tag].resize(n_vars);
2608 128557 : _jacobian_block_nonlocal_used[tag].resize(n_vars);
2609 266295 : for (MooseIndex(n_vars) i = 0; i < n_vars; ++i)
2610 : {
2611 137738 : if (!_block_diagonal_matrix)
2612 : {
2613 32056 : _sub_Kee[tag][i].resize(n_vars);
2614 32056 : _sub_Keg[tag][i].resize(n_vars);
2615 32056 : _sub_Ken[tag][i].resize(n_vars);
2616 32056 : _sub_Kne[tag][i].resize(n_vars);
2617 32056 : _sub_Knn[tag][i].resize(n_vars);
2618 32056 : _sub_Kll[tag][i].resize(n_vars);
2619 32056 : _sub_Kle[tag][i].resize(n_vars);
2620 32056 : _sub_Kln[tag][i].resize(n_vars);
2621 32056 : _sub_Kel[tag][i].resize(n_vars);
2622 32056 : _sub_Knl[tag][i].resize(n_vars);
2623 :
2624 32056 : _jacobian_block_used[tag][i].resize(n_vars);
2625 32056 : _jacobian_block_neighbor_used[tag][i].resize(n_vars);
2626 32056 : _jacobian_block_lower_used[tag][i].resize(n_vars);
2627 32056 : _jacobian_block_nonlocal_used[tag][i].resize(n_vars);
2628 : }
2629 : else
2630 : {
2631 105682 : _sub_Kee[tag][i].resize(1);
2632 105682 : _sub_Keg[tag][i].resize(1);
2633 105682 : _sub_Ken[tag][i].resize(1);
2634 105682 : _sub_Kne[tag][i].resize(1);
2635 105682 : _sub_Knn[tag][i].resize(1);
2636 105682 : _sub_Kll[tag][i].resize(1);
2637 105682 : _sub_Kle[tag][i].resize(1);
2638 105682 : _sub_Kln[tag][i].resize(1);
2639 105682 : _sub_Kel[tag][i].resize(1);
2640 105682 : _sub_Knl[tag][i].resize(1);
2641 :
2642 105682 : _jacobian_block_used[tag][i].resize(1);
2643 105682 : _jacobian_block_neighbor_used[tag][i].resize(1);
2644 105682 : _jacobian_block_lower_used[tag][i].resize(1);
2645 105682 : _jacobian_block_nonlocal_used[tag][i].resize(1);
2646 : }
2647 : }
2648 : }
2649 63135 : }
2650 :
2651 : void
2652 78 : Assembly::initNonlocalCoupling()
2653 : {
2654 78 : _cm_nonlocal_entry.clear();
2655 :
2656 78 : auto & vars = _sys.getVariables(_tid);
2657 :
2658 234 : for (auto & ivar : vars)
2659 : {
2660 156 : auto i = ivar->number();
2661 156 : auto ivar_start = _cm_nonlocal_entry.size();
2662 312 : for (unsigned int k = 0; k < ivar->count(); ++k)
2663 : {
2664 156 : unsigned int iv = i + k;
2665 268 : for (const auto & j : ConstCouplingRow(iv, _nonlocal_cm))
2666 112 : if (!_sys.isScalarVariable(j))
2667 : {
2668 112 : auto & jvar = _sys.getVariable(_tid, j);
2669 112 : auto pair = std::make_pair(ivar, &jvar);
2670 112 : auto c = ivar_start;
2671 : // check if the pair has been pushed or not
2672 112 : bool has_pair = false;
2673 146 : for (; c < _cm_nonlocal_entry.size(); ++c)
2674 34 : if (_cm_nonlocal_entry[c] == pair)
2675 : {
2676 0 : has_pair = true;
2677 0 : break;
2678 : }
2679 112 : if (!has_pair)
2680 112 : _cm_nonlocal_entry.push_back(pair);
2681 : }
2682 : }
2683 : }
2684 78 : }
2685 :
2686 : void
2687 82784409 : Assembly::prepareJacobianBlock()
2688 : {
2689 218911243 : for (const auto & it : _cm_ff_entry)
2690 : {
2691 136126834 : MooseVariableFEBase & ivar = *(it.first);
2692 136126834 : MooseVariableFEBase & jvar = *(it.second);
2693 :
2694 136126834 : unsigned int vi = ivar.number();
2695 136126834 : unsigned int vj = jvar.number();
2696 :
2697 136126834 : unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2698 :
2699 411034138 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2700 : {
2701 274907304 : jacobianBlock(vi, vj, LocalDataKey{}, tag)
2702 274907304 : .resize(ivar.dofIndices().size() * ivar.count(), jvar.dofIndices().size() * jcount);
2703 274907304 : jacobianBlockUsed(tag, vi, vj, false);
2704 : }
2705 : }
2706 82784409 : }
2707 :
2708 : void
2709 418601017 : Assembly::prepareResidual()
2710 : {
2711 418601017 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2712 886015799 : for (const auto & var : vars)
2713 1748462738 : for (auto & tag_Re : _sub_Re)
2714 1281047956 : tag_Re[var->number()].resize(var->dofIndices().size() * var->count());
2715 418601017 : }
2716 :
2717 : void
2718 542568 : Assembly::prepare()
2719 : {
2720 542568 : prepareJacobianBlock();
2721 542568 : prepareResidual();
2722 542568 : }
2723 :
2724 : void
2725 10364 : Assembly::prepareNonlocal()
2726 : {
2727 20992 : for (const auto & it : _cm_nonlocal_entry)
2728 : {
2729 10628 : MooseVariableFEBase & ivar = *(it.first);
2730 10628 : MooseVariableFEBase & jvar = *(it.second);
2731 :
2732 10628 : unsigned int vi = ivar.number();
2733 10628 : unsigned int vj = jvar.number();
2734 :
2735 10628 : unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2736 :
2737 31884 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2738 31884 : tag < _jacobian_block_nonlocal_used.size();
2739 : tag++)
2740 : {
2741 21256 : jacobianBlockNonlocal(vi, vj, LocalDataKey{}, tag)
2742 21256 : .resize(ivar.dofIndices().size() * ivar.count(), jvar.allDofIndices().size() * jcount);
2743 21256 : jacobianBlockNonlocalUsed(tag, vi, vj, false);
2744 : }
2745 : }
2746 10364 : }
2747 :
2748 : void
2749 1488 : Assembly::prepareVariable(MooseVariableFEBase * var)
2750 : {
2751 5624 : for (const auto & it : _cm_ff_entry)
2752 : {
2753 4136 : MooseVariableFEBase & ivar = *(it.first);
2754 4136 : MooseVariableFEBase & jvar = *(it.second);
2755 :
2756 4136 : unsigned int vi = ivar.number();
2757 4136 : unsigned int vj = jvar.number();
2758 :
2759 4136 : unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2760 :
2761 4136 : if (vi == var->number() || vj == var->number())
2762 : {
2763 8640 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2764 : {
2765 5760 : jacobianBlock(vi, vj, LocalDataKey{}, tag)
2766 5760 : .resize(ivar.dofIndices().size() * ivar.count(), jvar.dofIndices().size() * jcount);
2767 5760 : jacobianBlockUsed(tag, vi, vj, false);
2768 : }
2769 : }
2770 : }
2771 :
2772 5160 : for (auto & tag_Re : _sub_Re)
2773 3672 : tag_Re[var->number()].resize(var->dofIndices().size() * var->count());
2774 1488 : }
2775 :
2776 : void
2777 0 : Assembly::prepareVariableNonlocal(MooseVariableFEBase * var)
2778 : {
2779 0 : for (const auto & it : _cm_nonlocal_entry)
2780 : {
2781 0 : MooseVariableFEBase & ivar = *(it.first);
2782 0 : MooseVariableFEBase & jvar = *(it.second);
2783 :
2784 0 : unsigned int vi = ivar.number();
2785 0 : unsigned int vj = jvar.number();
2786 :
2787 0 : unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2788 :
2789 0 : if (vi == var->number() || vj == var->number())
2790 : {
2791 0 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2792 0 : tag < _jacobian_block_nonlocal_used.size();
2793 : tag++)
2794 : {
2795 0 : jacobianBlockNonlocal(vi, vj, LocalDataKey{}, tag)
2796 0 : .resize(ivar.dofIndices().size() * ivar.count(), jvar.allDofIndices().size() * jcount);
2797 0 : jacobianBlockNonlocalUsed(tag, vi, vj);
2798 : }
2799 : }
2800 : }
2801 0 : }
2802 :
2803 : void
2804 28843751 : Assembly::prepareNeighbor()
2805 : {
2806 88480614 : for (const auto & it : _cm_ff_entry)
2807 : {
2808 59636863 : MooseVariableFEBase & ivar = *(it.first);
2809 59636863 : MooseVariableFEBase & jvar = *(it.second);
2810 :
2811 59636863 : unsigned int vi = ivar.number();
2812 59636863 : unsigned int vj = jvar.number();
2813 :
2814 59636863 : unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2815 :
2816 178918077 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
2817 178918077 : tag < _jacobian_block_neighbor_used.size();
2818 : tag++)
2819 : {
2820 119281214 : jacobianBlockNeighbor(Moose::ElementNeighbor, vi, vj, LocalDataKey{}, tag)
2821 119281214 : .resize(ivar.dofIndices().size() * ivar.count(),
2822 119281214 : jvar.dofIndicesNeighbor().size() * jcount);
2823 :
2824 119281214 : jacobianBlockNeighbor(Moose::NeighborElement, vi, vj, LocalDataKey{}, tag)
2825 119281214 : .resize(ivar.dofIndicesNeighbor().size() * ivar.count(),
2826 119281214 : jvar.dofIndices().size() * jcount);
2827 :
2828 119281214 : jacobianBlockNeighbor(Moose::NeighborNeighbor, vi, vj, LocalDataKey{}, tag)
2829 119281214 : .resize(ivar.dofIndicesNeighbor().size() * ivar.count(),
2830 119281214 : jvar.dofIndicesNeighbor().size() * jcount);
2831 :
2832 119281214 : jacobianBlockNeighborUsed(tag, vi, vj, false);
2833 : }
2834 : }
2835 :
2836 28843751 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2837 67261704 : for (const auto & var : vars)
2838 136723390 : for (auto & tag_Rn : _sub_Rn)
2839 98305437 : tag_Rn[var->number()].resize(var->dofIndicesNeighbor().size() * var->count());
2840 28843751 : }
2841 :
2842 : void
2843 299795 : Assembly::prepareLowerD()
2844 : {
2845 1285175 : for (const auto & it : _cm_ff_entry)
2846 : {
2847 985380 : MooseVariableFEBase & ivar = *(it.first);
2848 985380 : MooseVariableFEBase & jvar = *(it.second);
2849 :
2850 985380 : unsigned int vi = ivar.number();
2851 985380 : unsigned int vj = jvar.number();
2852 :
2853 985380 : unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2854 :
2855 2956140 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
2856 : tag++)
2857 : {
2858 : // To cover all possible cases we should have 9 combinations below for every 2-permutation
2859 : // of Lower,Secondary,Primary. However, 4 cases will in general be covered by calls to
2860 : // prepare() and prepareNeighbor(). These calls will cover SecondarySecondary
2861 : // (ElementElement), SecondaryPrimary (ElementNeighbor), PrimarySecondary (NeighborElement),
2862 : // and PrimaryPrimary (NeighborNeighbor). With these covered we only need to prepare the 5
2863 : // remaining below
2864 :
2865 : // derivatives w.r.t. lower dimensional residuals
2866 1970760 : jacobianBlockMortar(Moose::LowerLower, vi, vj, LocalDataKey{}, tag)
2867 1970760 : .resize(ivar.dofIndicesLower().size() * ivar.count(),
2868 1970760 : jvar.dofIndicesLower().size() * jcount);
2869 :
2870 1970760 : jacobianBlockMortar(Moose::LowerSecondary, vi, vj, LocalDataKey{}, tag)
2871 1970760 : .resize(ivar.dofIndicesLower().size() * ivar.count(),
2872 1970760 : jvar.dofIndices().size() * jvar.count());
2873 :
2874 1970760 : jacobianBlockMortar(Moose::LowerPrimary, vi, vj, LocalDataKey{}, tag)
2875 1970760 : .resize(ivar.dofIndicesLower().size() * ivar.count(),
2876 1970760 : jvar.dofIndicesNeighbor().size() * jvar.count());
2877 :
2878 : // derivatives w.r.t. interior secondary residuals
2879 1970760 : jacobianBlockMortar(Moose::SecondaryLower, vi, vj, LocalDataKey{}, tag)
2880 1970760 : .resize(ivar.dofIndices().size() * ivar.count(),
2881 1970760 : jvar.dofIndicesLower().size() * jvar.count());
2882 :
2883 : // derivatives w.r.t. interior primary residuals
2884 1970760 : jacobianBlockMortar(Moose::PrimaryLower, vi, vj, LocalDataKey{}, tag)
2885 1970760 : .resize(ivar.dofIndicesNeighbor().size() * ivar.count(),
2886 1970760 : jvar.dofIndicesLower().size() * jvar.count());
2887 :
2888 1970760 : jacobianBlockLowerUsed(tag, vi, vj, false);
2889 : }
2890 : }
2891 :
2892 299795 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2893 805435 : for (const auto & var : vars)
2894 1590640 : for (auto & tag_Rl : _sub_Rl)
2895 1085000 : tag_Rl[var->number()].resize(var->dofIndicesLower().size() * var->count());
2896 299795 : }
2897 :
2898 : void
2899 0 : Assembly::prepareBlock(unsigned int ivar,
2900 : unsigned int jvar,
2901 : const std::vector<dof_id_type> & dof_indices)
2902 : {
2903 0 : const auto & iv = _sys.getVariable(_tid, ivar);
2904 0 : const auto & jv = _sys.getVariable(_tid, jvar);
2905 0 : const unsigned int ivn = iv.number();
2906 0 : const unsigned int jvn = jv.number();
2907 0 : const unsigned int icount = iv.count();
2908 0 : unsigned int jcount = jv.count();
2909 0 : if (ivn == jvn && _component_block_diagonal[ivn])
2910 0 : jcount = 1;
2911 :
2912 0 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2913 : {
2914 0 : jacobianBlock(ivn, jvn, LocalDataKey{}, tag)
2915 0 : .resize(dof_indices.size() * icount, dof_indices.size() * jcount);
2916 0 : jacobianBlockUsed(tag, ivn, jvn, false);
2917 : }
2918 :
2919 0 : for (auto & tag_Re : _sub_Re)
2920 0 : tag_Re[ivn].resize(dof_indices.size() * icount);
2921 0 : }
2922 :
2923 : void
2924 0 : Assembly::prepareBlockNonlocal(unsigned int ivar,
2925 : unsigned int jvar,
2926 : const std::vector<dof_id_type> & idof_indices,
2927 : const std::vector<dof_id_type> & jdof_indices)
2928 : {
2929 0 : const auto & iv = _sys.getVariable(_tid, ivar);
2930 0 : const auto & jv = _sys.getVariable(_tid, jvar);
2931 0 : const unsigned int ivn = iv.number();
2932 0 : const unsigned int jvn = jv.number();
2933 0 : const unsigned int icount = iv.count();
2934 0 : unsigned int jcount = jv.count();
2935 0 : if (ivn == jvn && _component_block_diagonal[ivn])
2936 0 : jcount = 1;
2937 :
2938 0 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2939 0 : tag < _jacobian_block_nonlocal_used.size();
2940 : tag++)
2941 : {
2942 0 : jacobianBlockNonlocal(ivn, jvn, LocalDataKey{}, tag)
2943 0 : .resize(idof_indices.size() * icount, jdof_indices.size() * jcount);
2944 :
2945 0 : jacobianBlockNonlocalUsed(tag, ivn, jvn, false);
2946 : }
2947 0 : }
2948 :
2949 : void
2950 8600239 : Assembly::prepareScalar()
2951 : {
2952 8600239 : const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
2953 8890286 : for (const auto & ivar : vars)
2954 : {
2955 290047 : auto idofs = ivar->dofIndices().size();
2956 :
2957 1149856 : for (auto & tag_Re : _sub_Re)
2958 859809 : tag_Re[ivar->number()].resize(idofs);
2959 :
2960 727700 : for (const auto & jvar : vars)
2961 : {
2962 437653 : auto jdofs = jvar->dofIndices().size();
2963 :
2964 1316869 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2965 : {
2966 879216 : jacobianBlock(ivar->number(), jvar->number(), LocalDataKey{}, tag).resize(idofs, jdofs);
2967 879216 : jacobianBlockUsed(tag, ivar->number(), jvar->number(), false);
2968 : }
2969 : }
2970 : }
2971 8600239 : }
2972 :
2973 : void
2974 186723 : Assembly::prepareOffDiagScalar()
2975 : {
2976 186723 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2977 186723 : const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
2978 :
2979 378761 : for (const auto & ivar : scalar_vars)
2980 : {
2981 192038 : auto idofs = ivar->dofIndices().size();
2982 :
2983 473087 : for (const auto & jvar : vars)
2984 : {
2985 281049 : auto jdofs = jvar->dofIndices().size() * jvar->count();
2986 843147 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2987 : {
2988 562098 : jacobianBlock(ivar->number(), jvar->number(), LocalDataKey{}, tag).resize(idofs, jdofs);
2989 562098 : jacobianBlockUsed(tag, ivar->number(), jvar->number(), false);
2990 :
2991 562098 : jacobianBlock(jvar->number(), ivar->number(), LocalDataKey{}, tag).resize(jdofs, idofs);
2992 562098 : jacobianBlockUsed(tag, jvar->number(), ivar->number(), false);
2993 : }
2994 : }
2995 : }
2996 186723 : }
2997 :
2998 : template <typename T>
2999 : void
3000 130886798 : Assembly::copyShapes(MooseVariableField<T> & v)
3001 : {
3002 130886798 : phi(v).shallowCopy(v.phi());
3003 130886798 : gradPhi(v).shallowCopy(v.gradPhi());
3004 130886798 : if (v.computingSecond())
3005 26616 : secondPhi(v).shallowCopy(v.secondPhi());
3006 130886798 : }
3007 :
3008 : void
3009 130886798 : Assembly::copyShapes(unsigned int var)
3010 : {
3011 130886798 : auto & v = _sys.getVariable(_tid, var);
3012 130886798 : if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3013 : {
3014 127890134 : auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3015 127890134 : copyShapes(v);
3016 : }
3017 2996664 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3018 : {
3019 1107654 : auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3020 1107654 : copyShapes(v);
3021 : }
3022 1889010 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3023 : {
3024 1889010 : auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3025 1889010 : copyShapes(v);
3026 1889010 : if (v.computingCurl())
3027 50512 : curlPhi(v).shallowCopy(v.curlPhi());
3028 1889010 : if (v.computingDiv())
3029 449280 : divPhi(v).shallowCopy(v.divPhi());
3030 : }
3031 : else
3032 0 : mooseError("Unsupported variable field type!");
3033 130886798 : }
3034 :
3035 : template <typename T>
3036 : void
3037 659288 : Assembly::copyFaceShapes(MooseVariableField<T> & v)
3038 : {
3039 659288 : phiFace(v).shallowCopy(v.phiFace());
3040 659288 : gradPhiFace(v).shallowCopy(v.gradPhiFace());
3041 659288 : if (v.computingSecond())
3042 5208 : secondPhiFace(v).shallowCopy(v.secondPhiFace());
3043 659288 : }
3044 :
3045 : void
3046 659288 : Assembly::copyFaceShapes(unsigned int var)
3047 : {
3048 659288 : auto & v = _sys.getVariable(_tid, var);
3049 659288 : if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3050 : {
3051 585047 : auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3052 585047 : copyFaceShapes(v);
3053 : }
3054 74241 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3055 : {
3056 12832 : auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3057 12832 : copyFaceShapes(v);
3058 : }
3059 61409 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3060 : {
3061 61409 : auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3062 61409 : copyFaceShapes(v);
3063 61409 : if (v.computingCurl())
3064 6380 : _vector_curl_phi_face.shallowCopy(v.curlPhi());
3065 61409 : if (v.computingDiv())
3066 18816 : _vector_div_phi_face.shallowCopy(v.divPhi());
3067 : }
3068 : else
3069 0 : mooseError("Unsupported variable field type!");
3070 659288 : }
3071 :
3072 : template <typename T>
3073 : void
3074 194810 : Assembly::copyNeighborShapes(MooseVariableField<T> & v)
3075 : {
3076 194810 : if (v.usesPhiNeighbor())
3077 : {
3078 194810 : phiFaceNeighbor(v).shallowCopy(v.phiFaceNeighbor());
3079 194810 : phiNeighbor(v).shallowCopy(v.phiNeighbor());
3080 : }
3081 194810 : if (v.usesGradPhiNeighbor())
3082 : {
3083 194810 : gradPhiFaceNeighbor(v).shallowCopy(v.gradPhiFaceNeighbor());
3084 194810 : gradPhiNeighbor(v).shallowCopy(v.gradPhiNeighbor());
3085 : }
3086 194810 : if (v.usesSecondPhiNeighbor())
3087 : {
3088 0 : secondPhiFaceNeighbor(v).shallowCopy(v.secondPhiFaceNeighbor());
3089 0 : secondPhiNeighbor(v).shallowCopy(v.secondPhiNeighbor());
3090 : }
3091 194810 : }
3092 :
3093 : void
3094 194810 : Assembly::copyNeighborShapes(unsigned int var)
3095 : {
3096 194810 : auto & v = _sys.getVariable(_tid, var);
3097 194810 : if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3098 : {
3099 190935 : auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3100 190935 : copyNeighborShapes(v);
3101 : }
3102 3875 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3103 : {
3104 3600 : auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3105 3600 : copyNeighborShapes(v);
3106 : }
3107 275 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3108 : {
3109 275 : auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3110 275 : copyNeighborShapes(v);
3111 : }
3112 : else
3113 0 : mooseError("Unsupported variable field type!");
3114 194810 : }
3115 :
3116 : DenseMatrix<Number> &
3117 358632242 : Assembly::jacobianBlockNeighbor(
3118 : Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
3119 : {
3120 358632242 : if (type == Moose::ElementElement)
3121 0 : jacobianBlockUsed(tag, ivar, jvar, true);
3122 : else
3123 358632242 : jacobianBlockNeighborUsed(tag, ivar, jvar, true);
3124 :
3125 358632242 : if (_block_diagonal_matrix)
3126 : {
3127 132101427 : switch (type)
3128 : {
3129 0 : default:
3130 : case Moose::ElementElement:
3131 0 : return _sub_Kee[tag][ivar][0];
3132 44034903 : case Moose::ElementNeighbor:
3133 44034903 : return _sub_Ken[tag][ivar][0];
3134 44028681 : case Moose::NeighborElement:
3135 44028681 : return _sub_Kne[tag][ivar][0];
3136 44037843 : case Moose::NeighborNeighbor:
3137 44037843 : return _sub_Knn[tag][ivar][0];
3138 : }
3139 : }
3140 : else
3141 : {
3142 226530815 : switch (type)
3143 : {
3144 0 : default:
3145 : case Moose::ElementElement:
3146 0 : return _sub_Kee[tag][ivar][jvar];
3147 75510501 : case Moose::ElementNeighbor:
3148 75510501 : return _sub_Ken[tag][ivar][jvar];
3149 75509813 : case Moose::NeighborElement:
3150 75509813 : return _sub_Kne[tag][ivar][jvar];
3151 75510501 : case Moose::NeighborNeighbor:
3152 75510501 : return _sub_Knn[tag][ivar][jvar];
3153 : }
3154 : }
3155 : }
3156 :
3157 : DenseMatrix<Number> &
3158 10591604 : Assembly::jacobianBlockMortar(Moose::ConstraintJacobianType type,
3159 : unsigned int ivar,
3160 : unsigned int jvar,
3161 : LocalDataKey,
3162 : TagID tag)
3163 : {
3164 10591604 : jacobianBlockLowerUsed(tag, ivar, jvar, true);
3165 10591604 : if (_block_diagonal_matrix)
3166 : {
3167 1926080 : switch (type)
3168 : {
3169 329552 : default:
3170 : case Moose::LowerLower:
3171 329552 : return _sub_Kll[tag][ivar][0];
3172 329168 : case Moose::LowerSecondary:
3173 329168 : return _sub_Kle[tag][ivar][0];
3174 328976 : case Moose::LowerPrimary:
3175 328976 : return _sub_Kln[tag][ivar][0];
3176 357192 : case Moose::SecondaryLower:
3177 357192 : return _sub_Kel[tag][ivar][0];
3178 56048 : case Moose::SecondarySecondary:
3179 56048 : return _sub_Kee[tag][ivar][0];
3180 56048 : case Moose::SecondaryPrimary:
3181 56048 : return _sub_Ken[tag][ivar][0];
3182 357000 : case Moose::PrimaryLower:
3183 357000 : return _sub_Knl[tag][ivar][0];
3184 56048 : case Moose::PrimarySecondary:
3185 56048 : return _sub_Kne[tag][ivar][0];
3186 56048 : case Moose::PrimaryPrimary:
3187 56048 : return _sub_Knn[tag][ivar][0];
3188 : }
3189 : }
3190 : else
3191 : {
3192 8665524 : switch (type)
3193 : {
3194 1714644 : default:
3195 : case Moose::LowerLower:
3196 1714644 : return _sub_Kll[tag][ivar][jvar];
3197 1713764 : case Moose::LowerSecondary:
3198 1713764 : return _sub_Kle[tag][ivar][jvar];
3199 1706196 : case Moose::LowerPrimary:
3200 1706196 : return _sub_Kln[tag][ivar][jvar];
3201 1715924 : case Moose::SecondaryLower:
3202 1715924 : return _sub_Kel[tag][ivar][jvar];
3203 26660 : case Moose::SecondarySecondary:
3204 26660 : return _sub_Kee[tag][ivar][jvar];
3205 26660 : case Moose::SecondaryPrimary:
3206 26660 : return _sub_Ken[tag][ivar][jvar];
3207 1708356 : case Moose::PrimaryLower:
3208 1708356 : return _sub_Knl[tag][ivar][jvar];
3209 26660 : case Moose::PrimarySecondary:
3210 26660 : return _sub_Kne[tag][ivar][jvar];
3211 26660 : case Moose::PrimaryPrimary:
3212 26660 : return _sub_Knn[tag][ivar][jvar];
3213 : }
3214 : }
3215 : }
3216 :
3217 : void
3218 1098238072 : Assembly::processLocalResidual(DenseVector<Number> & res_block,
3219 : std::vector<dof_id_type> & dof_indices,
3220 : const std::vector<Real> & scaling_factor,
3221 : bool is_nodal)
3222 : {
3223 : // For an array variable, ndof is the number of dofs of the zero-th component and
3224 : // ntdof is the number of dofs of all components.
3225 : // For standard or vector variables, ndof will be the same as ntdof.
3226 1098238072 : auto ndof = dof_indices.size();
3227 1098238072 : auto ntdof = res_block.size();
3228 1098238072 : auto count = ntdof / ndof;
3229 : mooseAssert(count == scaling_factor.size(), "Inconsistent of number of components");
3230 : mooseAssert(count * ndof == ntdof, "Inconsistent of number of components");
3231 1098238072 : if (count > 1)
3232 : {
3233 : // expanding dof indices
3234 4940709 : dof_indices.resize(ntdof);
3235 4940709 : unsigned int p = 0;
3236 17148087 : for (MooseIndex(count) j = 0; j < count; ++j)
3237 67883010 : for (MooseIndex(ndof) i = 0; i < ndof; ++i)
3238 : {
3239 55675632 : dof_indices[p] = dof_indices[i] + (is_nodal ? j : j * ndof);
3240 55675632 : res_block(p) *= scaling_factor[j];
3241 55675632 : ++p;
3242 : }
3243 : }
3244 : else
3245 : {
3246 1093297363 : if (scaling_factor[0] != 1.0)
3247 4901982 : res_block *= scaling_factor[0];
3248 : }
3249 :
3250 1098238072 : _dof_map.constrain_element_vector(res_block, dof_indices, false);
3251 1098238072 : }
3252 :
3253 : void
3254 14563758 : Assembly::addResidualBlock(NumericVector<Number> & residual,
3255 : DenseVector<Number> & res_block,
3256 : const std::vector<dof_id_type> & dof_indices,
3257 : const std::vector<Real> & scaling_factor,
3258 : bool is_nodal)
3259 : {
3260 14563758 : if (dof_indices.size() > 0 && res_block.size())
3261 : {
3262 7468928 : _temp_dof_indices = dof_indices;
3263 7468928 : _tmp_Re = res_block;
3264 7468928 : processLocalResidual(_tmp_Re, _temp_dof_indices, scaling_factor, is_nodal);
3265 7468928 : residual.add_vector(_tmp_Re, _temp_dof_indices);
3266 : }
3267 14563758 : }
3268 :
3269 : void
3270 1165837821 : Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
3271 : std::vector<dof_id_type> & cached_residual_rows,
3272 : DenseVector<Number> & res_block,
3273 : const std::vector<dof_id_type> & dof_indices,
3274 : const std::vector<Real> & scaling_factor,
3275 : bool is_nodal)
3276 : {
3277 1165837821 : if (dof_indices.size() > 0 && res_block.size())
3278 : {
3279 1090769144 : _temp_dof_indices = dof_indices;
3280 1090769144 : _tmp_Re = res_block;
3281 1090769144 : processLocalResidual(_tmp_Re, _temp_dof_indices, scaling_factor, is_nodal);
3282 :
3283 5181152016 : for (MooseIndex(_tmp_Re) i = 0; i < _tmp_Re.size(); i++)
3284 : {
3285 4090382872 : cached_residual_values.push_back(_tmp_Re(i));
3286 4090382872 : cached_residual_rows.push_back(_temp_dof_indices[i]);
3287 : }
3288 : }
3289 :
3290 1165837821 : res_block.zero();
3291 1165837821 : }
3292 :
3293 : void
3294 0 : Assembly::setResidualBlock(NumericVector<Number> & residual,
3295 : DenseVector<Number> & res_block,
3296 : const std::vector<dof_id_type> & dof_indices,
3297 : const std::vector<Real> & scaling_factor,
3298 : bool is_nodal)
3299 : {
3300 0 : if (dof_indices.size() > 0)
3301 : {
3302 0 : std::vector<dof_id_type> di(dof_indices);
3303 0 : _tmp_Re = res_block;
3304 0 : processLocalResidual(_tmp_Re, di, scaling_factor, is_nodal);
3305 0 : residual.insert(_tmp_Re, di);
3306 0 : }
3307 0 : }
3308 :
3309 : void
3310 782880 : Assembly::addResidual(const VectorTag & vector_tag)
3311 : {
3312 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3313 : "Non-residual tag in Assembly::addResidual");
3314 :
3315 782880 : auto & tag_Re = _sub_Re[vector_tag._type_id];
3316 782880 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3317 782880 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3318 2255232 : for (const auto & var : vars)
3319 1472352 : addResidualBlock(residual,
3320 1472352 : tag_Re[var->number()],
3321 1472352 : var->dofIndices(),
3322 1472352 : var->arrayScalingFactor(),
3323 1472352 : var->isNodal());
3324 782880 : }
3325 :
3326 : void
3327 270736 : Assembly::addResidual(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3328 : {
3329 1053616 : for (const auto & vector_tag : vector_tags)
3330 782880 : if (_sys.hasVector(vector_tag._id))
3331 782880 : addResidual(vector_tag);
3332 270736 : }
3333 :
3334 : void
3335 5172304 : Assembly::addResidualNeighbor(const VectorTag & vector_tag)
3336 : {
3337 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3338 : "Non-residual tag in Assembly::addResidualNeighbor");
3339 :
3340 5172304 : auto & tag_Rn = _sub_Rn[vector_tag._type_id];
3341 5172304 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3342 5172304 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3343 11638556 : for (const auto & var : vars)
3344 6466252 : addResidualBlock(residual,
3345 6466252 : tag_Rn[var->number()],
3346 6466252 : var->dofIndicesNeighbor(),
3347 6466252 : var->arrayScalingFactor(),
3348 6466252 : var->isNodal());
3349 5172304 : }
3350 :
3351 : void
3352 2080743 : Assembly::addResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3353 : {
3354 7253047 : for (const auto & vector_tag : vector_tags)
3355 5172304 : if (_sys.hasVector(vector_tag._id))
3356 5172304 : addResidualNeighbor(vector_tag);
3357 2080743 : }
3358 :
3359 : void
3360 5130576 : Assembly::addResidualLower(const VectorTag & vector_tag)
3361 : {
3362 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3363 : "Non-residual tag in Assembly::addResidualLower");
3364 :
3365 5130576 : auto & tag_Rl = _sub_Rl[vector_tag._type_id];
3366 5130576 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3367 5130576 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3368 11531280 : for (const auto & var : vars)
3369 6400704 : addResidualBlock(residual,
3370 6400704 : tag_Rl[var->number()],
3371 6400704 : var->dofIndicesLower(),
3372 6400704 : var->arrayScalingFactor(),
3373 6400704 : var->isNodal());
3374 5130576 : }
3375 :
3376 : void
3377 2055577 : Assembly::addResidualLower(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3378 : {
3379 7186153 : for (const auto & vector_tag : vector_tags)
3380 5130576 : if (_sys.hasVector(vector_tag._id))
3381 5130576 : addResidualLower(vector_tag);
3382 2055577 : }
3383 :
3384 : // private method, so no key required
3385 : void
3386 171809 : Assembly::addResidualScalar(const VectorTag & vector_tag)
3387 : {
3388 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3389 : "Non-residual tag in Assembly::addResidualScalar");
3390 :
3391 : // add the scalar variables residuals
3392 171809 : auto & tag_Re = _sub_Re[vector_tag._type_id];
3393 171809 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3394 171809 : const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
3395 396259 : for (const auto & var : vars)
3396 224450 : addResidualBlock(
3397 224450 : residual, tag_Re[var->number()], var->dofIndices(), var->arrayScalingFactor(), false);
3398 171809 : }
3399 :
3400 : void
3401 57299 : Assembly::addResidualScalar(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3402 : {
3403 229108 : for (const auto & vector_tag : vector_tags)
3404 171809 : if (_sys.hasVector(vector_tag._id))
3405 171809 : addResidualScalar(vector_tag);
3406 57299 : }
3407 :
3408 : void
3409 338499112 : Assembly::cacheResidual(GlobalDataKey, const std::vector<VectorTag> & tags)
3410 : {
3411 338499112 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3412 721895973 : for (const auto & var : vars)
3413 1425962564 : for (const auto & vector_tag : tags)
3414 1042565703 : if (_sys.hasVector(vector_tag._id))
3415 1042565703 : cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3416 1042565703 : _cached_residual_rows[vector_tag._type_id],
3417 1042565703 : _sub_Re[vector_tag._type_id][var->number()],
3418 1042565703 : var->dofIndices(),
3419 1042565703 : var->arrayScalingFactor(),
3420 1042565703 : var->isNodal());
3421 338499112 : }
3422 :
3423 : // private method, so no key required
3424 : void
3425 117727606 : Assembly::cacheResidual(dof_id_type dof, Real value, TagID tag_id)
3426 : {
3427 117727606 : const VectorTag & tag = _subproblem.getVectorTag(tag_id);
3428 :
3429 117727606 : _cached_residual_values[tag._type_id].push_back(value);
3430 117727606 : _cached_residual_rows[tag._type_id].push_back(dof);
3431 117727606 : }
3432 :
3433 : // private method, so no key required
3434 : void
3435 117296362 : Assembly::cacheResidual(dof_id_type dof, Real value, const std::set<TagID> & tags)
3436 : {
3437 235023968 : for (auto & tag : tags)
3438 117727606 : cacheResidual(dof, value, tag);
3439 117296362 : }
3440 :
3441 : void
3442 0 : Assembly::cacheResidualNodes(const DenseVector<Number> & res,
3443 : const std::vector<dof_id_type> & dof_index,
3444 : LocalDataKey,
3445 : TagID tag)
3446 : {
3447 : // Add the residual value and dof_index to cached_residual_values and cached_residual_rows
3448 : // respectively.
3449 : // This is used by NodalConstraint.C to cache the residual calculated for primary and secondary
3450 : // node.
3451 0 : const VectorTag & vector_tag = _subproblem.getVectorTag(tag);
3452 0 : for (MooseIndex(dof_index) i = 0; i < dof_index.size(); ++i)
3453 : {
3454 0 : _cached_residual_values[vector_tag._type_id].push_back(res(i));
3455 0 : _cached_residual_rows[vector_tag._type_id].push_back(dof_index[i]);
3456 : }
3457 0 : }
3458 :
3459 : void
3460 35796406 : Assembly::cacheResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & tags)
3461 : {
3462 35796406 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3463 83907714 : for (const auto & var : vars)
3464 170794490 : for (const auto & vector_tag : tags)
3465 122683182 : if (_sys.hasVector(vector_tag._id))
3466 122683182 : cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3467 122683182 : _cached_residual_rows[vector_tag._type_id],
3468 122683182 : _sub_Rn[vector_tag._type_id][var->number()],
3469 122683182 : var->dofIndicesNeighbor(),
3470 122683182 : var->arrayScalingFactor(),
3471 122683182 : var->isNodal());
3472 35796406 : }
3473 :
3474 : void
3475 181094 : Assembly::cacheResidualLower(GlobalDataKey, const std::vector<VectorTag> & tags)
3476 : {
3477 181094 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3478 459598 : for (const auto & var : vars)
3479 867440 : for (const auto & vector_tag : tags)
3480 588936 : if (_sys.hasVector(vector_tag._id))
3481 588936 : cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3482 588936 : _cached_residual_rows[vector_tag._type_id],
3483 588936 : _sub_Rl[vector_tag._type_id][var->number()],
3484 588936 : var->dofIndicesLower(),
3485 588936 : var->arrayScalingFactor(),
3486 588936 : var->isNodal());
3487 181094 : }
3488 :
3489 : void
3490 53017553 : Assembly::addCachedResiduals(GlobalDataKey, const std::vector<VectorTag> & tags)
3491 : {
3492 191638945 : for (const auto & vector_tag : tags)
3493 : {
3494 138621392 : if (!_sys.hasVector(vector_tag._id))
3495 : {
3496 0 : _cached_residual_values[vector_tag._type_id].clear();
3497 0 : _cached_residual_rows[vector_tag._type_id].clear();
3498 0 : continue;
3499 : }
3500 138621392 : addCachedResidualDirectly(_sys.getVector(vector_tag._id), GlobalDataKey{}, vector_tag);
3501 : }
3502 53017553 : }
3503 :
3504 : void
3505 10870 : Assembly::clearCachedResiduals(GlobalDataKey)
3506 : {
3507 42552 : for (const auto & vector_tag : _residual_vector_tags)
3508 31682 : clearCachedResiduals(vector_tag);
3509 10870 : }
3510 :
3511 : // private method, so no key required
3512 : void
3513 131065695 : Assembly::clearCachedResiduals(const VectorTag & vector_tag)
3514 : {
3515 131065695 : auto & values = _cached_residual_values[vector_tag._type_id];
3516 131065695 : auto & rows = _cached_residual_rows[vector_tag._type_id];
3517 :
3518 : mooseAssert(values.size() == rows.size(),
3519 : "Number of cached residuals and number of rows must match!");
3520 :
3521 : // Keep track of the largest size so we can use it to reserve and avoid
3522 : // as much dynamic allocation as possible
3523 131065695 : if (_max_cached_residuals < values.size())
3524 45082 : _max_cached_residuals = values.size();
3525 :
3526 : // Clear both vectors (keeps the capacity the same)
3527 131065695 : values.clear();
3528 131065695 : rows.clear();
3529 : // And then reserve: use 2 as a fudge factor to *really* avoid dynamic allocation!
3530 131065695 : values.reserve(_max_cached_residuals * 2);
3531 131065695 : rows.reserve(_max_cached_residuals * 2);
3532 131065695 : }
3533 :
3534 : void
3535 138642204 : Assembly::addCachedResidualDirectly(NumericVector<Number> & residual,
3536 : GlobalDataKey,
3537 : const VectorTag & vector_tag)
3538 : {
3539 138642204 : const auto & values = _cached_residual_values[vector_tag._type_id];
3540 138642204 : const auto & rows = _cached_residual_rows[vector_tag._type_id];
3541 :
3542 : mooseAssert(values.size() == rows.size(),
3543 : "Number of cached residuals and number of rows must match!");
3544 :
3545 138642204 : if (!values.empty())
3546 : {
3547 131034013 : residual.add_vector(values, rows);
3548 131034013 : clearCachedResiduals(vector_tag);
3549 : }
3550 138642204 : }
3551 :
3552 : void
3553 0 : Assembly::setResidual(NumericVector<Number> & residual, GlobalDataKey, const VectorTag & vector_tag)
3554 : {
3555 0 : auto & tag_Re = _sub_Re[vector_tag._type_id];
3556 0 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3557 0 : for (const auto & var : vars)
3558 0 : setResidualBlock(residual,
3559 0 : tag_Re[var->number()],
3560 0 : var->dofIndices(),
3561 0 : var->arrayScalingFactor(),
3562 0 : var->isNodal());
3563 0 : }
3564 :
3565 : void
3566 0 : Assembly::setResidualNeighbor(NumericVector<Number> & residual,
3567 : GlobalDataKey,
3568 : const VectorTag & vector_tag)
3569 : {
3570 0 : auto & tag_Rn = _sub_Rn[vector_tag._type_id];
3571 0 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3572 0 : for (const auto & var : vars)
3573 0 : setResidualBlock(residual,
3574 0 : tag_Rn[var->number()],
3575 0 : var->dofIndicesNeighbor(),
3576 0 : var->arrayScalingFactor(),
3577 0 : var->isNodal());
3578 0 : }
3579 :
3580 : // private method, so no key required
3581 : void
3582 459818 : Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
3583 : DenseMatrix<Number> & jac_block,
3584 : const MooseVariableBase & ivar,
3585 : const MooseVariableBase & jvar,
3586 : const std::vector<dof_id_type> & idof_indices,
3587 : const std::vector<dof_id_type> & jdof_indices)
3588 : {
3589 459818 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3590 86636 : return;
3591 373182 : if (jac_block.n() == 0 || jac_block.m() == 0)
3592 0 : return;
3593 :
3594 373182 : auto & scaling_factor = ivar.arrayScalingFactor();
3595 :
3596 760092 : for (unsigned int i = 0; i < ivar.count(); ++i)
3597 : {
3598 386910 : unsigned int iv = ivar.number();
3599 931704 : for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3600 : {
3601 544794 : unsigned int jv = jvar.number();
3602 544794 : if (jt < jv || jt >= jv + jvar.count())
3603 132348 : continue;
3604 412446 : unsigned int j = jt - jv;
3605 :
3606 412446 : auto di = ivar.componentDofIndices(idof_indices, i);
3607 412446 : auto dj = jvar.componentDofIndices(jdof_indices, j);
3608 412446 : auto indof = di.size();
3609 412446 : auto jndof = dj.size();
3610 :
3611 412446 : unsigned int jj = j;
3612 412446 : if (iv == jv && _component_block_diagonal[iv])
3613 : // here i must be equal to j
3614 335208 : jj = 0;
3615 :
3616 412446 : auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3617 412446 : if (scaling_factor[i] != 1.0)
3618 6108 : sub *= scaling_factor[i];
3619 :
3620 : // If we're computing the jacobian for automatically scaling variables we do not want
3621 : // to constrain the element matrix because it introduces 1s on the diagonal for the
3622 : // constrained dofs
3623 412446 : if (!_sys.computingScalingJacobian())
3624 412376 : _dof_map.constrain_element_matrix(sub, di, dj, false);
3625 :
3626 412446 : jacobian.add_matrix(sub, di, dj);
3627 412446 : }
3628 : }
3629 : }
3630 :
3631 : // private method, so no key required
3632 : void
3633 60452113 : Assembly::cacheJacobianBlock(DenseMatrix<Number> & jac_block,
3634 : const MooseVariableBase & ivar,
3635 : const MooseVariableBase & jvar,
3636 : const std::vector<dof_id_type> & idof_indices,
3637 : const std::vector<dof_id_type> & jdof_indices,
3638 : TagID tag)
3639 : {
3640 60452113 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3641 394052 : return;
3642 60058061 : if (jac_block.n() == 0 || jac_block.m() == 0)
3643 0 : return;
3644 60058061 : if (!_sys.hasMatrix(tag))
3645 0 : return;
3646 :
3647 60058061 : auto & scaling_factor = ivar.arrayScalingFactor();
3648 :
3649 120943150 : for (unsigned int i = 0; i < ivar.count(); ++i)
3650 : {
3651 60885089 : unsigned int iv = ivar.number();
3652 153175077 : for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3653 : {
3654 92289988 : unsigned int jv = jvar.number();
3655 92289988 : if (jt < jv || jt >= jv + jvar.count())
3656 28939755 : continue;
3657 63350233 : unsigned int j = jt - jv;
3658 :
3659 63350233 : auto di = ivar.componentDofIndices(idof_indices, i);
3660 63350233 : auto dj = jvar.componentDofIndices(jdof_indices, j);
3661 63350233 : auto indof = di.size();
3662 63350233 : auto jndof = dj.size();
3663 :
3664 63350233 : unsigned int jj = j;
3665 63350233 : if (iv == jv && _component_block_diagonal[iv])
3666 : // here i must be equal to j
3667 50135997 : jj = 0;
3668 :
3669 63350233 : auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3670 63350233 : if (scaling_factor[i] != 1.0)
3671 1794194 : sub *= scaling_factor[i];
3672 :
3673 : // If we're computing the jacobian for automatically scaling variables we do not want
3674 : // to constrain the element matrix because it introduces 1s on the diagonal for the
3675 : // constrained dofs
3676 63350233 : if (!_sys.computingScalingJacobian())
3677 63127198 : _dof_map.constrain_element_matrix(sub, di, dj, false);
3678 :
3679 342621978 : for (MooseIndex(di) i = 0; i < di.size(); i++)
3680 1684008149 : for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3681 : {
3682 1404736404 : _cached_jacobian_values[tag].push_back(sub(i, j));
3683 1404736404 : _cached_jacobian_rows[tag].push_back(di[i]);
3684 1404736404 : _cached_jacobian_cols[tag].push_back(dj[j]);
3685 : }
3686 63350233 : }
3687 : }
3688 :
3689 60058061 : jac_block.zero();
3690 : }
3691 :
3692 : // private method, so no key required
3693 : void
3694 770 : Assembly::cacheJacobianBlockNonzero(DenseMatrix<Number> & jac_block,
3695 : const MooseVariableBase & ivar,
3696 : const MooseVariableBase & jvar,
3697 : const std::vector<dof_id_type> & idof_indices,
3698 : const std::vector<dof_id_type> & jdof_indices,
3699 : TagID tag)
3700 : {
3701 770 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3702 0 : return;
3703 770 : if (jac_block.n() == 0 || jac_block.m() == 0)
3704 0 : return;
3705 770 : if (!_sys.hasMatrix(tag))
3706 0 : return;
3707 :
3708 770 : auto & scaling_factor = ivar.arrayScalingFactor();
3709 :
3710 1540 : for (unsigned int i = 0; i < ivar.count(); ++i)
3711 : {
3712 770 : unsigned int iv = ivar.number();
3713 2254 : for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3714 : {
3715 1484 : unsigned int jv = jvar.number();
3716 1484 : if (jt < jv || jt >= jv + jvar.count())
3717 714 : continue;
3718 770 : unsigned int j = jt - jv;
3719 :
3720 770 : auto di = ivar.componentDofIndices(idof_indices, i);
3721 770 : auto dj = jvar.componentDofIndices(jdof_indices, j);
3722 770 : auto indof = di.size();
3723 770 : auto jndof = dj.size();
3724 :
3725 770 : unsigned int jj = j;
3726 770 : if (iv == jv && _component_block_diagonal[iv])
3727 : // here i must be equal to j
3728 196 : jj = 0;
3729 :
3730 770 : auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3731 770 : if (scaling_factor[i] != 1.0)
3732 322 : sub *= scaling_factor[i];
3733 :
3734 770 : _dof_map.constrain_element_matrix(sub, di, dj, false);
3735 :
3736 3850 : for (MooseIndex(di) i = 0; i < di.size(); i++)
3737 232400 : for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3738 229320 : if (sub(i, j) != 0.0) // no storage allocated for unimplemented jacobian terms,
3739 : // maintaining maximum sparsity possible
3740 : {
3741 32648 : _cached_jacobian_values[tag].push_back(sub(i, j));
3742 32648 : _cached_jacobian_rows[tag].push_back(di[i]);
3743 32648 : _cached_jacobian_cols[tag].push_back(dj[j]);
3744 : }
3745 770 : }
3746 : }
3747 :
3748 770 : jac_block.zero();
3749 : }
3750 :
3751 : void
3752 6706460 : Assembly::cacheJacobianBlock(DenseMatrix<Number> & jac_block,
3753 : const std::vector<dof_id_type> & idof_indices,
3754 : const std::vector<dof_id_type> & jdof_indices,
3755 : Real scaling_factor,
3756 : LocalDataKey,
3757 : TagID tag)
3758 : {
3759 : // Only cache data when the matrix exists
3760 13392220 : if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m() &&
3761 6685760 : _sys.hasMatrix(tag))
3762 : {
3763 6685760 : std::vector<dof_id_type> di(idof_indices);
3764 6685760 : std::vector<dof_id_type> dj(jdof_indices);
3765 :
3766 : // If we're computing the jacobian for automatically scaling variables we do not want to
3767 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
3768 : // dofs
3769 6685760 : if (!_sys.computingScalingJacobian())
3770 6685760 : _dof_map.constrain_element_matrix(jac_block, di, dj, false);
3771 :
3772 6685760 : if (scaling_factor != 1.0)
3773 7521 : jac_block *= scaling_factor;
3774 :
3775 39563256 : for (MooseIndex(di) i = 0; i < di.size(); i++)
3776 217463582 : for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3777 : {
3778 184586086 : _cached_jacobian_values[tag].push_back(jac_block(i, j));
3779 184586086 : _cached_jacobian_rows[tag].push_back(di[i]);
3780 184586086 : _cached_jacobian_cols[tag].push_back(dj[j]);
3781 : }
3782 6685760 : }
3783 6706460 : jac_block.zero();
3784 6706460 : }
3785 :
3786 : Real
3787 0 : Assembly::elementVolume(const Elem * elem) const
3788 : {
3789 0 : FEType fe_type(elem->default_order(), LAGRANGE);
3790 0 : std::unique_ptr<FEBase> fe(FEBase::build(elem->dim(), fe_type));
3791 :
3792 : // references to the quadrature points and weights
3793 0 : const std::vector<Real> & JxW = fe->get_JxW();
3794 0 : const std::vector<Point> & q_points = fe->get_xyz();
3795 :
3796 : // The default quadrature rule should integrate the mass matrix,
3797 : // thus it should be plenty to compute the volume
3798 0 : QGauss qrule(elem->dim(), fe_type.default_quadrature_order());
3799 0 : fe->attach_quadrature_rule(&qrule);
3800 0 : fe->reinit(elem);
3801 :
3802 : // perform a sanity check to ensure that size of quad rule and size of q_points is
3803 : // identical
3804 : mooseAssert(qrule.n_points() == q_points.size(),
3805 : "The number of points in the quadrature rule doesn't match the number of passed-in "
3806 : "points in Assembly::setCoordinateTransformation");
3807 :
3808 : // compute the coordinate transformation
3809 0 : Real vol = 0;
3810 0 : for (unsigned int qp = 0; qp < qrule.n_points(); ++qp)
3811 : {
3812 : Real coord;
3813 0 : coordTransformFactor(_subproblem, elem->subdomain_id(), q_points[qp], coord);
3814 0 : vol += JxW[qp] * coord;
3815 : }
3816 0 : return vol;
3817 0 : }
3818 :
3819 : void
3820 3831499 : Assembly::addCachedJacobian(GlobalDataKey)
3821 : {
3822 : #ifndef NDEBUG
3823 : if (!_subproblem.checkNonlocalCouplingRequirement())
3824 : {
3825 : mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
3826 : "Error: Cached data sizes MUST be the same!");
3827 : for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3828 : mooseAssert(_cached_jacobian_rows[i].size() == _cached_jacobian_cols[i].size(),
3829 : "Error: Cached data sizes MUST be the same for a given tag!");
3830 : }
3831 : #endif
3832 :
3833 11563925 : for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3834 7732430 : if (_sys.hasMatrix(i))
3835 1823566469 : for (MooseIndex(_cached_jacobian_rows[i]) j = 0; j < _cached_jacobian_rows[i].size(); j++)
3836 3639439762 : _sys.getMatrix(i).add(_cached_jacobian_rows[i][j],
3837 1819719881 : _cached_jacobian_cols[i][j],
3838 1819719881 : _cached_jacobian_values[i][j]);
3839 :
3840 11563921 : for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3841 : {
3842 7732426 : if (!_sys.hasMatrix(i))
3843 3885838 : continue;
3844 :
3845 3846588 : if (_max_cached_jacobians < _cached_jacobian_values[i].size())
3846 46876 : _max_cached_jacobians = _cached_jacobian_values[i].size();
3847 :
3848 : // Try to be more efficient from now on
3849 : // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
3850 3846588 : _cached_jacobian_values[i].clear();
3851 3846588 : _cached_jacobian_values[i].reserve(_max_cached_jacobians * 2);
3852 :
3853 3846588 : _cached_jacobian_rows[i].clear();
3854 3846588 : _cached_jacobian_rows[i].reserve(_max_cached_jacobians * 2);
3855 :
3856 3846588 : _cached_jacobian_cols[i].clear();
3857 3846588 : _cached_jacobian_cols[i].reserve(_max_cached_jacobians * 2);
3858 : }
3859 3831495 : }
3860 :
3861 : inline void
3862 103486 : Assembly::addJacobianCoupledVarPair(const MooseVariableBase & ivar, const MooseVariableBase & jvar)
3863 : {
3864 103486 : auto i = ivar.number();
3865 103486 : auto j = jvar.number();
3866 311494 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
3867 208008 : if (jacobianBlockUsed(tag, i, j) && _sys.hasMatrix(tag))
3868 55316 : addJacobianBlock(_sys.getMatrix(tag),
3869 110632 : jacobianBlock(i, j, LocalDataKey{}, tag),
3870 : ivar,
3871 : jvar,
3872 55316 : ivar.dofIndices(),
3873 55316 : jvar.dofIndices());
3874 103486 : }
3875 :
3876 : void
3877 37388 : Assembly::addJacobian(GlobalDataKey)
3878 : {
3879 110174 : for (const auto & it : _cm_ff_entry)
3880 72786 : addJacobianCoupledVarPair(*it.first, *it.second);
3881 :
3882 37388 : for (const auto & it : _cm_sf_entry)
3883 0 : addJacobianCoupledVarPair(*it.first, *it.second);
3884 :
3885 37388 : for (const auto & it : _cm_fs_entry)
3886 0 : addJacobianCoupledVarPair(*it.first, *it.second);
3887 37388 : }
3888 :
3889 : void
3890 0 : Assembly::addJacobianNonlocal(GlobalDataKey)
3891 : {
3892 0 : for (const auto & it : _cm_nonlocal_entry)
3893 : {
3894 0 : auto ivar = it.first;
3895 0 : auto jvar = it.second;
3896 0 : auto i = ivar->number();
3897 0 : auto j = jvar->number();
3898 0 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
3899 0 : tag < _jacobian_block_nonlocal_used.size();
3900 : tag++)
3901 0 : if (jacobianBlockNonlocalUsed(tag, i, j) && _sys.hasMatrix(tag))
3902 0 : addJacobianBlock(_sys.getMatrix(tag),
3903 0 : jacobianBlockNonlocal(i, j, LocalDataKey{}, tag),
3904 : *ivar,
3905 : *jvar,
3906 0 : ivar->dofIndices(),
3907 : jvar->allDofIndices());
3908 : }
3909 0 : }
3910 :
3911 : void
3912 8244 : Assembly::addJacobianNeighbor(GlobalDataKey)
3913 : {
3914 38364 : for (const auto & it : _cm_ff_entry)
3915 : {
3916 30120 : auto ivar = it.first;
3917 30120 : auto jvar = it.second;
3918 30120 : auto i = ivar->number();
3919 30120 : auto j = jvar->number();
3920 90744 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3921 90744 : tag < _jacobian_block_neighbor_used.size();
3922 : tag++)
3923 60624 : if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
3924 : {
3925 22711 : addJacobianBlock(_sys.getMatrix(tag),
3926 22711 : jacobianBlockNeighbor(Moose::ElementNeighbor, i, j, LocalDataKey{}, tag),
3927 : *ivar,
3928 : *jvar,
3929 22711 : ivar->dofIndices(),
3930 22711 : jvar->dofIndicesNeighbor());
3931 :
3932 22711 : addJacobianBlock(_sys.getMatrix(tag),
3933 22711 : jacobianBlockNeighbor(Moose::NeighborElement, i, j, LocalDataKey{}, tag),
3934 : *ivar,
3935 : *jvar,
3936 22711 : ivar->dofIndicesNeighbor(),
3937 22711 : jvar->dofIndices());
3938 :
3939 22711 : addJacobianBlock(_sys.getMatrix(tag),
3940 45422 : jacobianBlockNeighbor(Moose::NeighborNeighbor, i, j, LocalDataKey{}, tag),
3941 : *ivar,
3942 : *jvar,
3943 22711 : ivar->dofIndicesNeighbor(),
3944 22711 : jvar->dofIndicesNeighbor());
3945 : }
3946 : }
3947 8244 : }
3948 :
3949 : void
3950 115562 : Assembly::addJacobianNeighborLowerD(GlobalDataKey)
3951 : {
3952 300559 : for (const auto & it : _cm_ff_entry)
3953 : {
3954 184997 : auto ivar = it.first;
3955 184997 : auto jvar = it.second;
3956 184997 : auto i = ivar->number();
3957 184997 : auto j = jvar->number();
3958 560079 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
3959 : tag++)
3960 375082 : if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
3961 : {
3962 7824 : addJacobianBlock(_sys.getMatrix(tag),
3963 7824 : jacobianBlockMortar(Moose::LowerLower, i, j, LocalDataKey{}, tag),
3964 : *ivar,
3965 : *jvar,
3966 7824 : ivar->dofIndicesLower(),
3967 7824 : jvar->dofIndicesLower());
3968 :
3969 7824 : addJacobianBlock(_sys.getMatrix(tag),
3970 7824 : jacobianBlockMortar(Moose::LowerSecondary, i, j, LocalDataKey{}, tag),
3971 : *ivar,
3972 : *jvar,
3973 7824 : ivar->dofIndicesLower(),
3974 7824 : jvar->dofIndicesNeighbor());
3975 :
3976 7824 : addJacobianBlock(_sys.getMatrix(tag),
3977 7824 : jacobianBlockMortar(Moose::LowerPrimary, i, j, LocalDataKey{}, tag),
3978 : *ivar,
3979 : *jvar,
3980 7824 : ivar->dofIndicesLower(),
3981 7824 : jvar->dofIndices());
3982 :
3983 7824 : addJacobianBlock(_sys.getMatrix(tag),
3984 7824 : jacobianBlockMortar(Moose::SecondaryLower, i, j, LocalDataKey{}, tag),
3985 : *ivar,
3986 : *jvar,
3987 7824 : ivar->dofIndicesNeighbor(),
3988 7824 : jvar->dofIndicesLower());
3989 :
3990 7824 : addJacobianBlock(_sys.getMatrix(tag),
3991 15648 : jacobianBlockMortar(Moose::PrimaryLower, i, j, LocalDataKey{}, tag),
3992 : *ivar,
3993 : *jvar,
3994 7824 : ivar->dofIndices(),
3995 7824 : jvar->dofIndicesLower());
3996 : }
3997 :
3998 560079 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3999 560079 : tag < _jacobian_block_neighbor_used.size();
4000 : tag++)
4001 375082 : if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
4002 : {
4003 92379 : addJacobianBlock(_sys.getMatrix(tag),
4004 92379 : jacobianBlockNeighbor(Moose::ElementNeighbor, i, j, LocalDataKey{}, tag),
4005 : *ivar,
4006 : *jvar,
4007 92379 : ivar->dofIndices(),
4008 92379 : jvar->dofIndicesNeighbor());
4009 :
4010 92379 : addJacobianBlock(_sys.getMatrix(tag),
4011 92379 : jacobianBlockNeighbor(Moose::NeighborElement, i, j, LocalDataKey{}, tag),
4012 : *ivar,
4013 : *jvar,
4014 92379 : ivar->dofIndicesNeighbor(),
4015 92379 : jvar->dofIndices());
4016 :
4017 92379 : addJacobianBlock(_sys.getMatrix(tag),
4018 184758 : jacobianBlockNeighbor(Moose::NeighborNeighbor, i, j, LocalDataKey{}, tag),
4019 : *ivar,
4020 : *jvar,
4021 92379 : ivar->dofIndicesNeighbor(),
4022 92379 : jvar->dofIndicesNeighbor());
4023 : }
4024 : }
4025 115562 : }
4026 :
4027 : void
4028 4816 : Assembly::addJacobianLowerD(GlobalDataKey)
4029 : {
4030 41840 : for (const auto & it : _cm_ff_entry)
4031 : {
4032 37024 : auto ivar = it.first;
4033 37024 : auto jvar = it.second;
4034 37024 : auto i = ivar->number();
4035 37024 : auto j = jvar->number();
4036 111072 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
4037 : tag++)
4038 74048 : if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
4039 : {
4040 6704 : addJacobianBlock(_sys.getMatrix(tag),
4041 6704 : jacobianBlockMortar(Moose::LowerLower, i, j, LocalDataKey{}, tag),
4042 : *ivar,
4043 : *jvar,
4044 6704 : ivar->dofIndicesLower(),
4045 6704 : jvar->dofIndicesLower());
4046 :
4047 6704 : addJacobianBlock(_sys.getMatrix(tag),
4048 6704 : jacobianBlockMortar(Moose::LowerSecondary, i, j, LocalDataKey{}, tag),
4049 : *ivar,
4050 : *jvar,
4051 6704 : ivar->dofIndicesLower(),
4052 6704 : jvar->dofIndices());
4053 :
4054 6704 : addJacobianBlock(_sys.getMatrix(tag),
4055 13408 : jacobianBlockMortar(Moose::SecondaryLower, i, j, LocalDataKey{}, tag),
4056 : *ivar,
4057 : *jvar,
4058 6704 : ivar->dofIndices(),
4059 6704 : jvar->dofIndicesLower());
4060 : }
4061 : }
4062 4816 : }
4063 :
4064 : void
4065 51041992 : Assembly::cacheJacobian(GlobalDataKey)
4066 : {
4067 129043144 : for (const auto & it : _cm_ff_entry)
4068 78001152 : cacheJacobianCoupledVarPair(*it.first, *it.second);
4069 :
4070 51293406 : for (const auto & it : _cm_fs_entry)
4071 251414 : cacheJacobianCoupledVarPair(*it.first, *it.second);
4072 :
4073 51293406 : for (const auto & it : _cm_sf_entry)
4074 251414 : cacheJacobianCoupledVarPair(*it.first, *it.second);
4075 51041992 : }
4076 :
4077 : // private method, so no key required
4078 : void
4079 78503980 : Assembly::cacheJacobianCoupledVarPair(const MooseVariableBase & ivar,
4080 : const MooseVariableBase & jvar)
4081 : {
4082 78503980 : auto i = ivar.number();
4083 78503980 : auto j = jvar.number();
4084 237498600 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
4085 158994620 : if (jacobianBlockUsed(tag, i, j) && _sys.hasMatrix(tag))
4086 60026527 : cacheJacobianBlock(jacobianBlock(i, j, LocalDataKey{}, tag),
4087 : ivar,
4088 : jvar,
4089 60026527 : ivar.dofIndices(),
4090 60026527 : jvar.dofIndices(),
4091 : tag);
4092 78503980 : }
4093 :
4094 : void
4095 5166 : Assembly::cacheJacobianNonlocal(GlobalDataKey)
4096 : {
4097 10458 : for (const auto & it : _cm_nonlocal_entry)
4098 : {
4099 5292 : auto ivar = it.first;
4100 5292 : auto jvar = it.second;
4101 5292 : auto i = ivar->number();
4102 5292 : auto j = jvar->number();
4103 15876 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
4104 15876 : tag < _jacobian_block_nonlocal_used.size();
4105 : tag++)
4106 10584 : if (jacobianBlockNonlocalUsed(tag, i, j) && _sys.hasMatrix(tag))
4107 1540 : cacheJacobianBlockNonzero(jacobianBlockNonlocal(i, j, LocalDataKey{}, tag),
4108 : *ivar,
4109 : *jvar,
4110 770 : ivar->dofIndices(),
4111 : jvar->allDofIndices(),
4112 : tag);
4113 : }
4114 5166 : }
4115 :
4116 : void
4117 9290 : Assembly::cacheJacobianNeighbor(GlobalDataKey)
4118 : {
4119 32404 : for (const auto & it : _cm_ff_entry)
4120 : {
4121 23114 : auto ivar = it.first;
4122 23114 : auto jvar = it.second;
4123 23114 : auto i = ivar->number();
4124 23114 : auto j = jvar->number();
4125 :
4126 69342 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
4127 69342 : tag < _jacobian_block_neighbor_used.size();
4128 : tag++)
4129 46228 : if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
4130 : {
4131 7726 : cacheJacobianBlock(jacobianBlockNeighbor(Moose::ElementNeighbor, i, j, LocalDataKey{}, tag),
4132 : *ivar,
4133 : *jvar,
4134 7726 : ivar->dofIndices(),
4135 7726 : jvar->dofIndicesNeighbor(),
4136 : tag);
4137 7726 : cacheJacobianBlock(jacobianBlockNeighbor(Moose::NeighborElement, i, j, LocalDataKey{}, tag),
4138 : *ivar,
4139 : *jvar,
4140 7726 : ivar->dofIndicesNeighbor(),
4141 7726 : jvar->dofIndices(),
4142 : tag);
4143 7726 : cacheJacobianBlock(
4144 15452 : jacobianBlockNeighbor(Moose::NeighborNeighbor, i, j, LocalDataKey{}, tag),
4145 : *ivar,
4146 : *jvar,
4147 7726 : ivar->dofIndicesNeighbor(),
4148 7726 : jvar->dofIndicesNeighbor(),
4149 : tag);
4150 : }
4151 : }
4152 9290 : }
4153 :
4154 : void
4155 84056 : Assembly::cacheJacobianMortar(GlobalDataKey)
4156 : {
4157 344172 : for (const auto & it : _cm_ff_entry)
4158 : {
4159 260116 : auto ivar = it.first;
4160 260116 : auto jvar = it.second;
4161 260116 : auto i = ivar->number();
4162 260116 : auto j = jvar->number();
4163 780348 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
4164 : tag++)
4165 520232 : if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
4166 : {
4167 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::LowerLower, i, j, LocalDataKey{}, tag),
4168 : *ivar,
4169 : *jvar,
4170 44712 : ivar->dofIndicesLower(),
4171 44712 : jvar->dofIndicesLower(),
4172 : tag);
4173 :
4174 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::LowerSecondary, i, j, LocalDataKey{}, tag),
4175 : *ivar,
4176 : *jvar,
4177 44712 : ivar->dofIndicesLower(),
4178 44712 : jvar->dofIndices(),
4179 : tag);
4180 :
4181 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::LowerPrimary, i, j, LocalDataKey{}, tag),
4182 : *ivar,
4183 : *jvar,
4184 44712 : ivar->dofIndicesLower(),
4185 44712 : jvar->dofIndicesNeighbor(),
4186 : tag);
4187 :
4188 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::SecondaryLower, i, j, LocalDataKey{}, tag),
4189 : *ivar,
4190 : *jvar,
4191 44712 : ivar->dofIndices(),
4192 44712 : jvar->dofIndicesLower(),
4193 : tag);
4194 :
4195 44712 : cacheJacobianBlock(
4196 44712 : jacobianBlockMortar(Moose::SecondarySecondary, i, j, LocalDataKey{}, tag),
4197 : *ivar,
4198 : *jvar,
4199 44712 : ivar->dofIndices(),
4200 44712 : jvar->dofIndices(),
4201 : tag);
4202 :
4203 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::SecondaryPrimary, i, j, LocalDataKey{}, tag),
4204 : *ivar,
4205 : *jvar,
4206 44712 : ivar->dofIndices(),
4207 44712 : jvar->dofIndicesNeighbor(),
4208 : tag);
4209 :
4210 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::PrimaryLower, i, j, LocalDataKey{}, tag),
4211 : *ivar,
4212 : *jvar,
4213 44712 : ivar->dofIndicesNeighbor(),
4214 44712 : jvar->dofIndicesLower(),
4215 : tag);
4216 :
4217 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::PrimarySecondary, i, j, LocalDataKey{}, tag),
4218 : *ivar,
4219 : *jvar,
4220 44712 : ivar->dofIndicesNeighbor(),
4221 44712 : jvar->dofIndices(),
4222 : tag);
4223 :
4224 44712 : cacheJacobianBlock(jacobianBlockMortar(Moose::PrimaryPrimary, i, j, LocalDataKey{}, tag),
4225 : *ivar,
4226 : *jvar,
4227 44712 : ivar->dofIndicesNeighbor(),
4228 44712 : jvar->dofIndicesNeighbor(),
4229 : tag);
4230 : }
4231 : }
4232 84056 : }
4233 :
4234 : void
4235 65544 : Assembly::addJacobianBlockTags(SparseMatrix<Number> & jacobian,
4236 : unsigned int ivar,
4237 : unsigned int jvar,
4238 : const DofMap & dof_map,
4239 : std::vector<dof_id_type> & dof_indices,
4240 : GlobalDataKey,
4241 : const std::set<TagID> & tags)
4242 : {
4243 164376 : for (auto tag : tags)
4244 98832 : addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices, GlobalDataKey{}, tag);
4245 65544 : }
4246 :
4247 : void
4248 98832 : Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
4249 : unsigned int ivar,
4250 : unsigned int jvar,
4251 : const DofMap & dof_map,
4252 : std::vector<dof_id_type> & dof_indices,
4253 : GlobalDataKey,
4254 : TagID tag)
4255 : {
4256 98832 : if (dof_indices.size() == 0)
4257 0 : return;
4258 98832 : if (!(*_cm)(ivar, jvar))
4259 0 : return;
4260 :
4261 98832 : auto & iv = _sys.getVariable(_tid, ivar);
4262 98832 : auto & jv = _sys.getVariable(_tid, jvar);
4263 98832 : auto & scaling_factor = iv.arrayScalingFactor();
4264 :
4265 98832 : const unsigned int ivn = iv.number();
4266 98832 : const unsigned int jvn = jv.number();
4267 98832 : auto & ke = jacobianBlock(ivn, jvn, LocalDataKey{}, tag);
4268 :
4269 : // It is guaranteed by design iv.number <= ivar since iv is obtained
4270 : // through SystemBase::getVariable with ivar.
4271 : // Most of times ivar will just be equal to iv.number except for array variables,
4272 : // where ivar could be a number for a component of an array variable but calling
4273 : // getVariable will return the array variable that has the number of the 0th component.
4274 : // It is the same for jvar.
4275 98832 : const unsigned int i = ivar - ivn;
4276 98832 : const unsigned int j = jvar - jvn;
4277 :
4278 : // DoF indices are independently given
4279 98832 : auto di = dof_indices;
4280 98832 : auto dj = dof_indices;
4281 :
4282 98832 : auto indof = di.size();
4283 98832 : auto jndof = dj.size();
4284 :
4285 98832 : unsigned int jj = j;
4286 98832 : if (ivar == jvar && _component_block_diagonal[ivn])
4287 94304 : jj = 0;
4288 :
4289 98832 : auto sub = ke.sub_matrix(i * indof, indof, jj * jndof, jndof);
4290 : // If we're computing the jacobian for automatically scaling variables we do not want to
4291 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4292 : // dofs
4293 98832 : if (!_sys.computingScalingJacobian())
4294 98832 : dof_map.constrain_element_matrix(sub, di, dj, false);
4295 :
4296 98832 : if (scaling_factor[i] != 1.0)
4297 0 : sub *= scaling_factor[i];
4298 :
4299 98832 : jacobian.add_matrix(sub, di, dj);
4300 98832 : }
4301 :
4302 : void
4303 0 : Assembly::addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
4304 : const unsigned int ivar,
4305 : const unsigned int jvar,
4306 : const DofMap & dof_map,
4307 : const std::vector<dof_id_type> & idof_indices,
4308 : const std::vector<dof_id_type> & jdof_indices,
4309 : GlobalDataKey,
4310 : const TagID tag)
4311 : {
4312 0 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
4313 0 : return;
4314 0 : if (jacobian.n() == 0 || jacobian.m() == 0)
4315 0 : return;
4316 0 : if (!(*_cm)(ivar, jvar))
4317 0 : return;
4318 :
4319 0 : auto & iv = _sys.getVariable(_tid, ivar);
4320 0 : auto & jv = _sys.getVariable(_tid, jvar);
4321 0 : auto & scaling_factor = iv.arrayScalingFactor();
4322 :
4323 0 : const unsigned int ivn = iv.number();
4324 0 : const unsigned int jvn = jv.number();
4325 0 : auto & keg = jacobianBlockNonlocal(ivn, jvn, LocalDataKey{}, tag);
4326 :
4327 : // It is guaranteed by design iv.number <= ivar since iv is obtained
4328 : // through SystemBase::getVariable with ivar.
4329 : // Most of times ivar will just be equal to iv.number except for array variables,
4330 : // where ivar could be a number for a component of an array variable but calling
4331 : // getVariable will return the array variable that has the number of the 0th component.
4332 : // It is the same for jvar.
4333 0 : const unsigned int i = ivar - ivn;
4334 0 : const unsigned int j = jvar - jvn;
4335 :
4336 : // DoF indices are independently given
4337 0 : auto di = idof_indices;
4338 0 : auto dj = jdof_indices;
4339 :
4340 0 : auto indof = di.size();
4341 0 : auto jndof = dj.size();
4342 :
4343 0 : unsigned int jj = j;
4344 0 : if (ivar == jvar && _component_block_diagonal[ivn])
4345 0 : jj = 0;
4346 :
4347 0 : auto sub = keg.sub_matrix(i * indof, indof, jj * jndof, jndof);
4348 : // If we're computing the jacobian for automatically scaling variables we do not want to
4349 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4350 : // dofs
4351 0 : if (!_sys.computingScalingJacobian())
4352 0 : dof_map.constrain_element_matrix(sub, di, dj, false);
4353 :
4354 0 : if (scaling_factor[i] != 1.0)
4355 0 : sub *= scaling_factor[i];
4356 :
4357 0 : jacobian.add_matrix(sub, di, dj);
4358 0 : }
4359 :
4360 : void
4361 0 : Assembly::addJacobianBlockNonlocalTags(SparseMatrix<Number> & jacobian,
4362 : const unsigned int ivar,
4363 : const unsigned int jvar,
4364 : const DofMap & dof_map,
4365 : const std::vector<dof_id_type> & idof_indices,
4366 : const std::vector<dof_id_type> & jdof_indices,
4367 : GlobalDataKey,
4368 : const std::set<TagID> & tags)
4369 : {
4370 0 : for (auto tag : tags)
4371 0 : addJacobianBlockNonlocal(
4372 0 : jacobian, ivar, jvar, dof_map, idof_indices, jdof_indices, GlobalDataKey{}, tag);
4373 0 : }
4374 :
4375 : void
4376 2688 : Assembly::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
4377 : const unsigned int ivar,
4378 : const unsigned int jvar,
4379 : const DofMap & dof_map,
4380 : std::vector<dof_id_type> & dof_indices,
4381 : std::vector<dof_id_type> & neighbor_dof_indices,
4382 : GlobalDataKey,
4383 : const TagID tag)
4384 : {
4385 2688 : if (dof_indices.size() == 0 && neighbor_dof_indices.size() == 0)
4386 0 : return;
4387 2688 : if (!(*_cm)(ivar, jvar))
4388 0 : return;
4389 :
4390 2688 : auto & iv = _sys.getVariable(_tid, ivar);
4391 2688 : auto & jv = _sys.getVariable(_tid, jvar);
4392 2688 : auto & scaling_factor = iv.arrayScalingFactor();
4393 :
4394 2688 : const unsigned int ivn = iv.number();
4395 2688 : const unsigned int jvn = jv.number();
4396 2688 : auto & ken = jacobianBlockNeighbor(Moose::ElementNeighbor, ivn, jvn, LocalDataKey{}, tag);
4397 2688 : auto & kne = jacobianBlockNeighbor(Moose::NeighborElement, ivn, jvn, LocalDataKey{}, tag);
4398 2688 : auto & knn = jacobianBlockNeighbor(Moose::NeighborNeighbor, ivn, jvn, LocalDataKey{}, tag);
4399 :
4400 : // It is guaranteed by design iv.number <= ivar since iv is obtained
4401 : // through SystemBase::getVariable with ivar.
4402 : // Most of times ivar will just be equal to iv.number except for array variables,
4403 : // where ivar could be a number for a component of an array variable but calling
4404 : // getVariable will return the array variable that has the number of the 0th component.
4405 : // It is the same for jvar.
4406 2688 : const unsigned int i = ivar - ivn;
4407 2688 : const unsigned int j = jvar - jvn;
4408 : // DoF indices are independently given
4409 2688 : auto dc = dof_indices;
4410 2688 : auto dn = neighbor_dof_indices;
4411 2688 : auto cndof = dc.size();
4412 2688 : auto nndof = dn.size();
4413 :
4414 2688 : unsigned int jj = j;
4415 2688 : if (ivar == jvar && _component_block_diagonal[ivn])
4416 2688 : jj = 0;
4417 :
4418 2688 : auto suben = ken.sub_matrix(i * cndof, cndof, jj * nndof, nndof);
4419 2688 : auto subne = kne.sub_matrix(i * nndof, nndof, jj * cndof, cndof);
4420 2688 : auto subnn = knn.sub_matrix(i * nndof, nndof, jj * nndof, nndof);
4421 :
4422 : // If we're computing the jacobian for automatically scaling variables we do not want to
4423 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4424 : // dofs
4425 2688 : if (!_sys.computingScalingJacobian())
4426 : {
4427 2688 : dof_map.constrain_element_matrix(suben, dc, dn, false);
4428 2688 : dof_map.constrain_element_matrix(subne, dn, dc, false);
4429 2688 : dof_map.constrain_element_matrix(subnn, dn, dn, false);
4430 : }
4431 :
4432 2688 : if (scaling_factor[i] != 1.0)
4433 : {
4434 0 : suben *= scaling_factor[i];
4435 0 : subne *= scaling_factor[i];
4436 0 : subnn *= scaling_factor[i];
4437 : }
4438 :
4439 2688 : jacobian.add_matrix(suben, dc, dn);
4440 2688 : jacobian.add_matrix(subne, dn, dc);
4441 2688 : jacobian.add_matrix(subnn, dn, dn);
4442 2688 : }
4443 :
4444 : void
4445 1344 : Assembly::addJacobianNeighborTags(SparseMatrix<Number> & jacobian,
4446 : const unsigned int ivar,
4447 : const unsigned int jvar,
4448 : const DofMap & dof_map,
4449 : std::vector<dof_id_type> & dof_indices,
4450 : std::vector<dof_id_type> & neighbor_dof_indices,
4451 : GlobalDataKey,
4452 : const std::set<TagID> & tags)
4453 : {
4454 4032 : for (const auto tag : tags)
4455 2688 : addJacobianNeighbor(
4456 5376 : jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices, GlobalDataKey{}, tag);
4457 1344 : }
4458 :
4459 : void
4460 13840 : Assembly::addJacobianScalar(GlobalDataKey)
4461 : {
4462 31574 : for (const auto & it : _cm_ss_entry)
4463 17734 : addJacobianCoupledVarPair(*it.first, *it.second);
4464 13840 : }
4465 :
4466 : void
4467 34764 : Assembly::addJacobianOffDiagScalar(unsigned int ivar, GlobalDataKey)
4468 : {
4469 34764 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
4470 34764 : MooseVariableScalar & var_i = _sys.getScalarVariable(_tid, ivar);
4471 47730 : for (const auto & var_j : vars)
4472 12966 : addJacobianCoupledVarPair(var_i, *var_j);
4473 34764 : }
4474 :
4475 : void
4476 238614908 : Assembly::cacheJacobian(
4477 : numeric_index_type i, numeric_index_type j, Real value, LocalDataKey, TagID tag)
4478 : {
4479 238614908 : _cached_jacobian_rows[tag].push_back(i);
4480 238614908 : _cached_jacobian_cols[tag].push_back(j);
4481 238614908 : _cached_jacobian_values[tag].push_back(value);
4482 238614908 : }
4483 :
4484 : void
4485 238562004 : Assembly::cacheJacobian(numeric_index_type i,
4486 : numeric_index_type j,
4487 : Real value,
4488 : LocalDataKey,
4489 : const std::set<TagID> & tags)
4490 : {
4491 492102516 : for (auto tag : tags)
4492 253540512 : if (_sys.hasMatrix(tag))
4493 238614908 : cacheJacobian(i, j, value, LocalDataKey{}, tag);
4494 238562004 : }
4495 :
4496 : void
4497 394261 : Assembly::setCachedJacobian(GlobalDataKey)
4498 : {
4499 1194593 : for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4500 800332 : if (_sys.hasMatrix(tag))
4501 : {
4502 : // First zero the rows (including the diagonals) to prepare for
4503 : // setting the cached values.
4504 395648 : _sys.getMatrix(tag).zero_rows(_cached_jacobian_rows[tag], 0.0);
4505 :
4506 : // TODO: Use SparseMatrix::set_values() for efficiency
4507 8635269 : for (MooseIndex(_cached_jacobian_values) i = 0; i < _cached_jacobian_values[tag].size(); ++i)
4508 16479242 : _sys.getMatrix(tag).set(_cached_jacobian_rows[tag][i],
4509 8239621 : _cached_jacobian_cols[tag][i],
4510 8239621 : _cached_jacobian_values[tag][i]);
4511 : }
4512 :
4513 394261 : clearCachedJacobian();
4514 394261 : }
4515 :
4516 : void
4517 0 : Assembly::zeroCachedJacobian(GlobalDataKey)
4518 : {
4519 0 : for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4520 0 : if (_sys.hasMatrix(tag))
4521 0 : _sys.getMatrix(tag).zero_rows(_cached_jacobian_rows[tag], 0.0);
4522 :
4523 0 : clearCachedJacobian();
4524 0 : }
4525 :
4526 : void
4527 394261 : Assembly::clearCachedJacobian()
4528 : {
4529 1194593 : for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4530 : {
4531 800332 : _cached_jacobian_rows[tag].clear();
4532 800332 : _cached_jacobian_cols[tag].clear();
4533 800332 : _cached_jacobian_values[tag].clear();
4534 : }
4535 394261 : }
4536 :
4537 : void
4538 0 : Assembly::modifyWeightsDueToXFEM(const Elem * elem)
4539 : {
4540 : mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
4541 :
4542 0 : if (_current_qrule == _current_qrule_arbitrary)
4543 0 : return;
4544 :
4545 0 : MooseArray<Real> xfem_weight_multipliers;
4546 0 : if (_xfem->getXFEMWeights(xfem_weight_multipliers, elem, _current_qrule, _current_q_points))
4547 : {
4548 : mooseAssert(xfem_weight_multipliers.size() == _current_JxW.size(),
4549 : "Size of weight multipliers in xfem doesn't match number of quadrature points");
4550 0 : for (unsigned i = 0; i < xfem_weight_multipliers.size(); i++)
4551 0 : _current_JxW[i] = _current_JxW[i] * xfem_weight_multipliers[i];
4552 :
4553 0 : xfem_weight_multipliers.release();
4554 : }
4555 0 : }
4556 :
4557 : void
4558 0 : Assembly::modifyFaceWeightsDueToXFEM(const Elem * elem, unsigned int side)
4559 : {
4560 : mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
4561 :
4562 0 : if (_current_qrule_face == _current_qrule_arbitrary)
4563 0 : return;
4564 :
4565 0 : MooseArray<Real> xfem_face_weight_multipliers;
4566 0 : if (_xfem->getXFEMFaceWeights(
4567 0 : xfem_face_weight_multipliers, elem, _current_qrule_face, _current_q_points_face, side))
4568 : {
4569 : mooseAssert(xfem_face_weight_multipliers.size() == _current_JxW_face.size(),
4570 : "Size of weight multipliers in xfem doesn't match number of quadrature points");
4571 0 : for (unsigned i = 0; i < xfem_face_weight_multipliers.size(); i++)
4572 0 : _current_JxW_face[i] = _current_JxW_face[i] * xfem_face_weight_multipliers[i];
4573 :
4574 0 : xfem_face_weight_multipliers.release();
4575 : }
4576 0 : }
4577 :
4578 : void
4579 412 : Assembly::hasScalingVector()
4580 : {
4581 412 : _scaling_vector = &_sys.getVector("scaling_factors");
4582 412 : }
4583 :
4584 : void
4585 0 : Assembly::modifyArbitraryWeights(const std::vector<Real> & weights)
4586 : {
4587 : mooseAssert(_current_qrule == _current_qrule_arbitrary, "Rule should be arbitrary");
4588 : mooseAssert(weights.size() == _current_physical_points.size(), "Size mismatch");
4589 :
4590 0 : for (MooseIndex(weights.size()) i = 0; i < weights.size(); ++i)
4591 0 : _current_JxW[i] = weights[i];
4592 0 : }
4593 :
4594 : template <>
4595 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4596 1571 : Assembly::fePhi<VectorValue<Real>>(FEType type) const
4597 : {
4598 1571 : buildVectorFE(type);
4599 1571 : return _vector_fe_shape_data[type]->_phi;
4600 : }
4601 :
4602 : template <>
4603 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4604 1571 : Assembly::feGradPhi<VectorValue<Real>>(FEType type) const
4605 : {
4606 1571 : buildVectorFE(type);
4607 1571 : return _vector_fe_shape_data[type]->_grad_phi;
4608 : }
4609 :
4610 : template <>
4611 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4612 0 : Assembly::feSecondPhi<VectorValue<Real>>(FEType type) const
4613 : {
4614 0 : _need_second_derivative.insert(type);
4615 0 : buildVectorFE(type);
4616 0 : return _vector_fe_shape_data[type]->_second_phi;
4617 : }
4618 :
4619 : template <>
4620 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4621 3142 : Assembly::fePhiLower<VectorValue<Real>>(FEType type) const
4622 : {
4623 3142 : buildVectorLowerDFE(type);
4624 3142 : return _vector_fe_shape_data_lower[type]->_phi;
4625 : }
4626 :
4627 : template <>
4628 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4629 0 : Assembly::feDualPhiLower<VectorValue<Real>>(FEType type) const
4630 : {
4631 0 : buildVectorDualLowerDFE(type);
4632 0 : return _vector_fe_shape_data_dual_lower[type]->_phi;
4633 : }
4634 :
4635 : template <>
4636 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4637 3142 : Assembly::feGradPhiLower<VectorValue<Real>>(FEType type) const
4638 : {
4639 3142 : buildVectorLowerDFE(type);
4640 3142 : return _vector_fe_shape_data_lower[type]->_grad_phi;
4641 : }
4642 :
4643 : template <>
4644 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4645 0 : Assembly::feGradDualPhiLower<VectorValue<Real>>(FEType type) const
4646 : {
4647 0 : buildVectorDualLowerDFE(type);
4648 0 : return _vector_fe_shape_data_dual_lower[type]->_grad_phi;
4649 : }
4650 :
4651 : template <>
4652 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4653 1571 : Assembly::fePhiFace<VectorValue<Real>>(FEType type) const
4654 : {
4655 1571 : buildVectorFaceFE(type);
4656 1571 : return _vector_fe_shape_data_face[type]->_phi;
4657 : }
4658 :
4659 : template <>
4660 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4661 1571 : Assembly::feGradPhiFace<VectorValue<Real>>(FEType type) const
4662 : {
4663 1571 : buildVectorFaceFE(type);
4664 1571 : return _vector_fe_shape_data_face[type]->_grad_phi;
4665 : }
4666 :
4667 : template <>
4668 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4669 0 : Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type) const
4670 : {
4671 0 : _need_second_derivative.insert(type);
4672 0 : buildVectorFaceFE(type);
4673 :
4674 : // If we're building for a face we probably need to build for a
4675 : // neighbor while _need_second_derivative is set;
4676 : // onInterface/reinitNeighbor/etc don't distinguish
4677 0 : buildVectorFaceNeighborFE(type);
4678 :
4679 0 : return _vector_fe_shape_data_face[type]->_second_phi;
4680 : }
4681 :
4682 : template <>
4683 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4684 1571 : Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type) const
4685 : {
4686 1571 : buildVectorNeighborFE(type);
4687 1571 : return _vector_fe_shape_data_neighbor[type]->_phi;
4688 : }
4689 :
4690 : template <>
4691 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4692 1571 : Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type) const
4693 : {
4694 1571 : buildVectorNeighborFE(type);
4695 1571 : return _vector_fe_shape_data_neighbor[type]->_grad_phi;
4696 : }
4697 :
4698 : template <>
4699 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4700 0 : Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type) const
4701 : {
4702 0 : _need_second_derivative_neighbor.insert(type);
4703 0 : buildVectorNeighborFE(type);
4704 0 : return _vector_fe_shape_data_neighbor[type]->_second_phi;
4705 : }
4706 :
4707 : template <>
4708 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4709 1571 : Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4710 : {
4711 1571 : buildVectorFaceNeighborFE(type);
4712 1571 : return _vector_fe_shape_data_face_neighbor[type]->_phi;
4713 : }
4714 :
4715 : template <>
4716 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4717 1571 : Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4718 : {
4719 1571 : buildVectorFaceNeighborFE(type);
4720 1571 : return _vector_fe_shape_data_face_neighbor[type]->_grad_phi;
4721 : }
4722 :
4723 : template <>
4724 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4725 0 : Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4726 : {
4727 0 : _need_second_derivative_neighbor.insert(type);
4728 0 : buildVectorFaceNeighborFE(type);
4729 0 : return _vector_fe_shape_data_face_neighbor[type]->_second_phi;
4730 : }
4731 :
4732 : template <>
4733 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4734 57245 : Assembly::feCurlPhi<VectorValue<Real>>(FEType type) const
4735 : {
4736 57245 : _need_curl.insert(type);
4737 57245 : buildVectorFE(type);
4738 57245 : return _vector_fe_shape_data[type]->_curl_phi;
4739 : }
4740 :
4741 : template <>
4742 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4743 209 : Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type) const
4744 : {
4745 209 : _need_curl.insert(type);
4746 209 : buildVectorFaceFE(type);
4747 :
4748 : // If we're building for a face we probably need to build for a
4749 : // neighbor while _need_curl is set;
4750 : // onInterface/reinitNeighbor/etc don't distinguish
4751 209 : buildVectorFaceNeighborFE(type);
4752 :
4753 209 : return _vector_fe_shape_data_face[type]->_curl_phi;
4754 : }
4755 :
4756 : template <>
4757 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4758 0 : Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type) const
4759 : {
4760 0 : _need_curl.insert(type);
4761 0 : buildVectorNeighborFE(type);
4762 0 : return _vector_fe_shape_data_neighbor[type]->_curl_phi;
4763 : }
4764 :
4765 : template <>
4766 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4767 0 : Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4768 : {
4769 0 : _need_curl.insert(type);
4770 0 : buildVectorFaceNeighborFE(type);
4771 :
4772 0 : return _vector_fe_shape_data_face_neighbor[type]->_curl_phi;
4773 : }
4774 :
4775 : template <>
4776 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4777 469068 : Assembly::feDivPhi<VectorValue<Real>>(FEType type) const
4778 : {
4779 469068 : _need_div.insert(type);
4780 469068 : buildVectorFE(type);
4781 469068 : return _vector_fe_shape_data[type]->_div_phi;
4782 : }
4783 :
4784 : template <>
4785 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4786 312 : Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const
4787 : {
4788 312 : _need_face_div.insert(type);
4789 312 : buildVectorFaceFE(type);
4790 :
4791 : // If we're building for a face we probably need to build for a
4792 : // neighbor while _need_face_div is set;
4793 : // onInterface/reinitNeighbor/etc don't distinguish
4794 312 : buildVectorFaceNeighborFE(type);
4795 :
4796 312 : return _vector_fe_shape_data_face[type]->_div_phi;
4797 : }
4798 :
4799 : template <>
4800 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4801 0 : Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const
4802 : {
4803 0 : _need_neighbor_div.insert(type);
4804 0 : buildVectorNeighborFE(type);
4805 0 : return _vector_fe_shape_data_neighbor[type]->_div_phi;
4806 : }
4807 :
4808 : template <>
4809 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4810 0 : Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4811 : {
4812 0 : _need_face_neighbor_div.insert(type);
4813 0 : buildVectorFaceNeighborFE(type);
4814 0 : return _vector_fe_shape_data_face_neighbor[type]->_div_phi;
4815 : }
4816 :
4817 : const MooseArray<ADReal> &
4818 26 : Assembly::adCurvatures() const
4819 : {
4820 26 : _calculate_curvatures = true;
4821 26 : const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
4822 26 : const FEType helper_type(helper_order, LAGRANGE);
4823 : // Must prerequest the second derivatives. Sadly because there is only one
4824 : // _need_second_derivative map for both volumetric and face FE objects we must request both here
4825 26 : feSecondPhi<Real>(helper_type);
4826 26 : feSecondPhiFace<Real>(helper_type);
4827 26 : return _ad_curvatures;
4828 : }
4829 :
4830 : void
4831 66805 : Assembly::helpersRequestData()
4832 : {
4833 259602 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
4834 : {
4835 192797 : _holder_fe_helper[dim]->get_phi();
4836 192797 : _holder_fe_helper[dim]->get_dphi();
4837 192797 : _holder_fe_helper[dim]->get_xyz();
4838 192797 : _holder_fe_helper[dim]->get_JxW();
4839 :
4840 192797 : _holder_fe_face_helper[dim]->get_phi();
4841 192797 : _holder_fe_face_helper[dim]->get_dphi();
4842 192797 : _holder_fe_face_helper[dim]->get_xyz();
4843 192797 : _holder_fe_face_helper[dim]->get_JxW();
4844 192797 : _holder_fe_face_helper[dim]->get_normals();
4845 :
4846 192797 : _holder_fe_face_neighbor_helper[dim]->get_xyz();
4847 192797 : _holder_fe_face_neighbor_helper[dim]->get_JxW();
4848 192797 : _holder_fe_face_neighbor_helper[dim]->get_normals();
4849 :
4850 192797 : _holder_fe_neighbor_helper[dim]->get_xyz();
4851 192797 : _holder_fe_neighbor_helper[dim]->get_JxW();
4852 : }
4853 :
4854 192797 : for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
4855 : {
4856 : // We need these computations in order to compute correct lower-d element volumes in
4857 : // curvilinear coordinates
4858 125992 : _holder_fe_lower_helper[dim]->get_xyz();
4859 125992 : _holder_fe_lower_helper[dim]->get_JxW();
4860 : }
4861 66805 : }
4862 :
4863 : void
4864 234 : Assembly::havePRefinement(const std::unordered_set<FEFamily> & disable_families)
4865 : {
4866 234 : if (_have_p_refinement)
4867 : // Already performed tasks for p-refinement
4868 0 : return;
4869 :
4870 234 : const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
4871 234 : const FEType helper_type(helper_order, LAGRANGE);
4872 : auto process_fe =
4873 9984 : [&disable_families](const unsigned int num_dimensionalities, auto & fe_container)
4874 : {
4875 2340 : if (!disable_families.empty())
4876 7540 : for (const auto dim : make_range(num_dimensionalities))
4877 : {
4878 5590 : auto fe_container_it = fe_container.find(dim);
4879 5590 : if (fe_container_it != fe_container.end())
4880 11349 : for (auto & [fe_type, fe_ptr] : fe_container_it->second)
4881 7644 : if (disable_families.count(fe_type.family))
4882 2665 : fe_ptr->add_p_level_in_reinit(false);
4883 : }
4884 2574 : };
4885 1170 : auto process_fe_and_helpers = [process_fe, &helper_type](auto & unique_helper_container,
4886 : auto & helper_container,
4887 : const unsigned int num_dimensionalities,
4888 : const bool user_added_helper_type,
4889 5707 : auto & fe_container)
4890 : {
4891 1170 : unique_helper_container.resize(num_dimensionalities);
4892 4511 : for (const auto dim : make_range(num_dimensionalities))
4893 : {
4894 3341 : auto & unique_helper = unique_helper_container[dim];
4895 3341 : unique_helper = FEGenericBase<Real>::build(dim, helper_type);
4896 : // don't participate in p-refinement
4897 3341 : unique_helper->add_p_level_in_reinit(false);
4898 3341 : helper_container[dim] = unique_helper.get();
4899 :
4900 : // If the user did not request the helper type then we should erase it from our FE container
4901 : // so that they're not penalized (in the "we should be able to do p-refinement sense") for
4902 : // our perhaps silly helpers
4903 3341 : if (!user_added_helper_type)
4904 : {
4905 2366 : auto & fe_container_dim = libmesh_map_find(fe_container, dim);
4906 2366 : auto fe_it = fe_container_dim.find(helper_type);
4907 : mooseAssert(fe_it != fe_container_dim.end(), "We should have the helper type");
4908 2366 : delete fe_it->second;
4909 2366 : fe_container_dim.erase(fe_it);
4910 : }
4911 : }
4912 :
4913 1170 : process_fe(num_dimensionalities, fe_container);
4914 1170 : };
4915 :
4916 : // Handle scalar field families
4917 234 : process_fe_and_helpers(_unique_fe_helper,
4918 234 : _holder_fe_helper,
4919 234 : _mesh_dimension + 1,
4920 234 : _user_added_fe_of_helper_type,
4921 234 : _fe);
4922 234 : process_fe_and_helpers(_unique_fe_face_helper,
4923 234 : _holder_fe_face_helper,
4924 234 : _mesh_dimension + 1,
4925 234 : _user_added_fe_face_of_helper_type,
4926 234 : _fe_face);
4927 234 : process_fe_and_helpers(_unique_fe_face_neighbor_helper,
4928 234 : _holder_fe_face_neighbor_helper,
4929 234 : _mesh_dimension + 1,
4930 234 : _user_added_fe_face_neighbor_of_helper_type,
4931 234 : _fe_face_neighbor);
4932 234 : process_fe_and_helpers(_unique_fe_neighbor_helper,
4933 234 : _holder_fe_neighbor_helper,
4934 234 : _mesh_dimension + 1,
4935 234 : _user_added_fe_neighbor_of_helper_type,
4936 234 : _fe_neighbor);
4937 234 : process_fe_and_helpers(_unique_fe_lower_helper,
4938 234 : _holder_fe_lower_helper,
4939 : _mesh_dimension,
4940 234 : _user_added_fe_lower_of_helper_type,
4941 234 : _fe_lower);
4942 : // Handle vector field families
4943 234 : process_fe(_mesh_dimension + 1, _vector_fe);
4944 234 : process_fe(_mesh_dimension + 1, _vector_fe_face);
4945 234 : process_fe(_mesh_dimension + 1, _vector_fe_neighbor);
4946 234 : process_fe(_mesh_dimension + 1, _vector_fe_face_neighbor);
4947 234 : process_fe(_mesh_dimension, _vector_fe_lower);
4948 :
4949 234 : helpersRequestData();
4950 :
4951 234 : _have_p_refinement = true;
4952 : }
4953 :
4954 : template void coordTransformFactor<Point, Real>(const SubProblem & s,
4955 : SubdomainID sub_id,
4956 : const Point & point,
4957 : Real & factor,
4958 : SubdomainID neighbor_sub_id);
4959 : template void coordTransformFactor<ADPoint, ADReal>(const SubProblem & s,
4960 : SubdomainID sub_id,
4961 : const ADPoint & point,
4962 : ADReal & factor,
4963 : SubdomainID neighbor_sub_id);
4964 : template void coordTransformFactor<Point, Real>(const MooseMesh & mesh,
4965 : SubdomainID sub_id,
4966 : const Point & point,
4967 : Real & factor,
4968 : SubdomainID neighbor_sub_id);
4969 : template void coordTransformFactor<ADPoint, ADReal>(const MooseMesh & mesh,
4970 : SubdomainID sub_id,
4971 : const ADPoint & point,
4972 : ADReal & factor,
4973 : SubdomainID neighbor_sub_id);
4974 :
4975 : template <>
4976 : const MooseArray<Moose::GenericType<Point, false>> &
4977 23 : Assembly::genericQPoints<false>() const
4978 : {
4979 23 : return qPoints();
4980 : }
4981 :
4982 : template <>
4983 : const MooseArray<Moose::GenericType<Point, true>> &
4984 11 : Assembly::genericQPoints<true>() const
4985 : {
4986 11 : return adQPoints();
4987 : }
|
