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 :
37 : #include <algorithm>
38 :
39 : using namespace libMesh;
40 :
41 : template <typename P, typename C>
42 : void
43 1807666843 : coordTransformFactor(const SubProblem & s,
44 : const SubdomainID sub_id,
45 : const P & point,
46 : C & factor,
47 : const SubdomainID neighbor_sub_id)
48 : {
49 1807666843 : coordTransformFactor(s.mesh(), sub_id, point, factor, neighbor_sub_id);
50 1807666843 : }
51 :
52 : template <typename P, typename C>
53 : void
54 1810930003 : coordTransformFactor(const MooseMesh & mesh,
55 : const SubdomainID sub_id,
56 : const P & point,
57 : C & factor,
58 : const SubdomainID libmesh_dbg_var(neighbor_sub_id))
59 : {
60 : mooseAssert(neighbor_sub_id != libMesh::Elem::invalid_subdomain_id
61 : ? mesh.getCoordSystem(sub_id) == mesh.getCoordSystem(neighbor_sub_id)
62 : : true,
63 : "Coordinate systems must be the same between element and neighbor");
64 1810930003 : const auto coord_type = mesh.getCoordSystem(sub_id);
65 :
66 1810930003 : if (coord_type == Moose::COORD_RZ)
67 : {
68 3487480 : if (mesh.usingGeneralAxisymmetricCoordAxes())
69 : {
70 1041955 : const auto & axis = mesh.getGeneralAxisymmetricCoordAxis(sub_id);
71 1041955 : MooseMeshUtils::coordTransformFactorRZGeneral(point, axis, factor);
72 : }
73 : else
74 2445525 : MooseMeshUtils::coordTransformFactor(
75 : point, factor, coord_type, mesh.getAxisymmetricRadialCoord());
76 : }
77 : else
78 1807442523 : MooseMeshUtils::coordTransformFactor(point, factor, coord_type, libMesh::invalid_uint);
79 1810930003 : }
80 :
81 72852 : Assembly::Assembly(SystemBase & sys, THREAD_ID tid)
82 72852 : : _sys(sys),
83 145704 : _subproblem(_sys.subproblem()),
84 72852 : _displaced(dynamic_cast<DisplacedSystem *>(&sys) ? true : false),
85 72852 : _nonlocal_cm(_subproblem.nonlocalCouplingMatrix(_sys.number())),
86 72852 : _computing_residual(_subproblem.currentlyComputingResidual()),
87 72852 : _computing_jacobian(_subproblem.currentlyComputingJacobian()),
88 72852 : _computing_residual_and_jacobian(_subproblem.currentlyComputingResidualAndJacobian()),
89 72852 : _dof_map(_sys.dofMap()),
90 72852 : _tid(tid),
91 72852 : _mesh(sys.mesh()),
92 72852 : _mesh_dimension(_mesh.dimension()),
93 72852 : _helper_type(_mesh.hasSecondOrderElements() ? SECOND : FIRST, LAGRANGE),
94 72852 : _user_added_fe_of_helper_type(false),
95 72852 : _user_added_fe_face_of_helper_type(false),
96 72852 : _user_added_fe_face_neighbor_of_helper_type(false),
97 72852 : _user_added_fe_neighbor_of_helper_type(false),
98 72852 : _user_added_fe_lower_of_helper_type(false),
99 72852 : _building_helpers(false),
100 72852 : _current_qrule(nullptr),
101 72852 : _current_qrule_volume(nullptr),
102 72852 : _current_qrule_arbitrary(nullptr),
103 72852 : _coord_type(Moose::COORD_XYZ),
104 72852 : _current_qrule_face(nullptr),
105 72852 : _current_qface_arbitrary(nullptr),
106 72852 : _current_qrule_neighbor(nullptr),
107 72852 : _need_JxW_neighbor(false),
108 72852 : _qrule_msm(nullptr),
109 72852 : _custom_mortar_qrule(false),
110 72852 : _current_qrule_lower(nullptr),
111 :
112 72852 : _current_elem(nullptr),
113 72852 : _current_elem_volume(0),
114 72852 : _current_side(0),
115 72852 : _current_side_elem(nullptr),
116 72852 : _current_side_volume(0),
117 72852 : _current_neighbor_elem(nullptr),
118 72852 : _current_neighbor_side(0),
119 72852 : _current_neighbor_side_elem(nullptr),
120 72852 : _need_neighbor_elem_volume(false),
121 72852 : _current_neighbor_volume(0),
122 72852 : _current_node(nullptr),
123 72852 : _current_neighbor_node(nullptr),
124 72852 : _current_elem_volume_computed(false),
125 72852 : _current_side_volume_computed(false),
126 :
127 72852 : _current_lower_d_elem(nullptr),
128 72852 : _current_neighbor_lower_d_elem(nullptr),
129 72852 : _need_lower_d_elem_volume(false),
130 72852 : _need_neighbor_lower_d_elem_volume(false),
131 72852 : _need_dual(false),
132 :
133 72852 : _residual_vector_tags(_subproblem.getVectorTags(Moose::VECTOR_TAG_RESIDUAL)),
134 145704 : _cached_residual_values(2), // The 2 is for TIME and NONTIME
135 145704 : _cached_residual_rows(2), // The 2 is for TIME and NONTIME
136 72852 : _max_cached_residuals(0),
137 72852 : _max_cached_jacobians(0),
138 :
139 72852 : _block_diagonal_matrix(false),
140 72852 : _calculate_xyz(false),
141 72852 : _calculate_face_xyz(false),
142 72852 : _calculate_curvatures(false),
143 72852 : _calculate_ad_coord(false),
144 291408 : _have_p_refinement(false)
145 : {
146 72852 : const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
147 72852 : _building_helpers = true;
148 : // Build fe's for the helpers
149 72852 : buildFE(FEType(helper_order, LAGRANGE));
150 72852 : buildFaceFE(FEType(helper_order, LAGRANGE));
151 72852 : buildNeighborFE(FEType(helper_order, LAGRANGE));
152 72852 : buildFaceNeighborFE(FEType(helper_order, LAGRANGE));
153 72852 : buildLowerDFE(FEType(helper_order, LAGRANGE));
154 72852 : _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 285263 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
159 : {
160 212411 : _holder_fe_helper[dim] = _fe[dim][FEType(helper_order, LAGRANGE)];
161 212411 : _holder_fe_face_helper[dim] = _fe_face[dim][FEType(helper_order, LAGRANGE)];
162 212411 : _holder_fe_face_neighbor_helper[dim] = _fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)];
163 212411 : _holder_fe_neighbor_helper[dim] = _fe_neighbor[dim][FEType(helper_order, LAGRANGE)];
164 : }
165 :
166 212411 : for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
167 139559 : _holder_fe_lower_helper[dim] = _fe_lower[dim][FEType(helper_order, LAGRANGE)];
168 :
169 : // request phi, dphi, xyz, JxW, etc. data
170 72852 : helpersRequestData();
171 :
172 : // For 3D mortar, mortar segments are always TRI3 elements so we want FIRST LAGRANGE regardless
173 : // of discretization
174 72852 : _fe_msm = (_mesh_dimension == 2)
175 145704 : ? FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(helper_order, LAGRANGE))
176 72852 : : FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(FIRST, LAGRANGE));
177 : // This FE object should not take part in p-refinement
178 72852 : _fe_msm->add_p_level_in_reinit(false);
179 72852 : _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 72852 : _fe_msm->get_xyz();
183 :
184 72852 : _extra_elem_ids.resize(_mesh.getMesh().n_elem_integers() + 1);
185 72852 : _neighbor_extra_elem_ids.resize(_mesh.getMesh().n_elem_integers() + 1);
186 72852 : }
187 :
188 137198 : Assembly::~Assembly()
189 : {
190 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
191 487943 : for (auto & it : _fe[dim])
192 287173 : delete it.second;
193 :
194 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
195 487943 : for (auto & it : _fe_face[dim])
196 287173 : delete it.second;
197 :
198 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
199 487943 : for (auto & it : _fe_neighbor[dim])
200 287173 : delete it.second;
201 :
202 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
203 487943 : for (auto & it : _fe_face_neighbor[dim])
204 287173 : delete it.second;
205 :
206 200770 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
207 315435 : for (auto & it : _fe_lower[dim])
208 183264 : delete it.second;
209 :
210 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
211 204329 : for (auto & it : _vector_fe[dim])
212 3559 : delete it.second;
213 :
214 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
215 204329 : for (auto & it : _vector_fe_face[dim])
216 3559 : delete it.second;
217 :
218 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
219 204329 : for (auto & it : _vector_fe_neighbor[dim])
220 3559 : delete it.second;
221 :
222 269369 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
223 204329 : for (auto & it : _vector_fe_face_neighbor[dim])
224 3559 : delete it.second;
225 :
226 200770 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
227 134126 : for (auto & it : _vector_fe_lower[dim])
228 1955 : delete it.second;
229 :
230 149997 : for (auto & it : _ad_grad_phi_data)
231 81398 : it.second.release();
232 :
233 70017 : for (auto & it : _ad_vector_grad_phi_data)
234 1418 : it.second.release();
235 :
236 145483 : for (auto & it : _ad_grad_phi_data_face)
237 76884 : it.second.release();
238 :
239 70017 : for (auto & it : _ad_vector_grad_phi_data_face)
240 1418 : it.second.release();
241 :
242 68599 : _current_physical_points.release();
243 :
244 68599 : _coord.release();
245 68599 : _coord_neighbor.release();
246 68599 : _coord_msm.release();
247 :
248 68599 : _ad_JxW.release();
249 68599 : _ad_q_points.release();
250 68599 : _ad_JxW_face.release();
251 68599 : _ad_normals.release();
252 68599 : _ad_q_points_face.release();
253 68599 : _curvatures.release();
254 68599 : _ad_curvatures.release();
255 68599 : _ad_coord.release();
256 :
257 68599 : delete _qrule_msm;
258 137198 : }
259 :
260 : const MooseArray<Real> &
261 107 : Assembly::JxWNeighbor() const
262 : {
263 107 : _need_JxW_neighbor = true;
264 107 : return _current_JxW_neighbor;
265 : }
266 :
267 : void
268 464414 : Assembly::buildFE(FEType type) const
269 : {
270 464414 : if (!_building_helpers && type == _helper_type)
271 208869 : _user_added_fe_of_helper_type = true;
272 :
273 464414 : if (!_fe_shape_data[type])
274 103187 : _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 1859384 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
278 : {
279 1394970 : if (!_fe[dim][type])
280 303391 : _fe[dim][type] = FEGenericBase<Real>::build(dim, type).release();
281 :
282 1394970 : _fe[dim][type]->get_phi();
283 1394970 : _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 1394970 : _fe[dim][type]->get_xyz();
288 1394970 : if (_need_second_derivative.count(type))
289 72340 : _fe[dim][type]->get_d2phi();
290 : }
291 464414 : }
292 :
293 : void
294 439785 : Assembly::buildFaceFE(FEType type) const
295 : {
296 439785 : if (!_building_helpers && type == _helper_type)
297 203684 : _user_added_fe_face_of_helper_type = true;
298 :
299 439785 : if (!_fe_shape_data_face[type])
300 103187 : _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 1760808 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
304 : {
305 1321023 : if (!_fe_face[dim][type])
306 303391 : _fe_face[dim][type] = FEGenericBase<Real>::build(dim, type).release();
307 :
308 1321023 : _fe_face[dim][type]->get_phi();
309 1321023 : _fe_face[dim][type]->get_dphi();
310 1321023 : if (_need_second_derivative.count(type))
311 12838 : _fe_face[dim][type]->get_d2phi();
312 : }
313 439785 : }
314 :
315 : void
316 435549 : Assembly::buildNeighborFE(FEType type) const
317 : {
318 435549 : if (!_building_helpers && type == _helper_type)
319 203572 : _user_added_fe_neighbor_of_helper_type = true;
320 :
321 435549 : if (!_fe_shape_data_neighbor[type])
322 103187 : _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 1743851 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
326 : {
327 1308302 : if (!_fe_neighbor[dim][type])
328 303391 : _fe_neighbor[dim][type] = FEGenericBase<Real>::build(dim, type).release();
329 :
330 1308302 : _fe_neighbor[dim][type]->get_phi();
331 1308302 : _fe_neighbor[dim][type]->get_dphi();
332 1308302 : if (_need_second_derivative_neighbor.count(type))
333 117 : _fe_neighbor[dim][type]->get_d2phi();
334 : }
335 435549 : }
336 :
337 : void
338 435549 : Assembly::buildFaceNeighborFE(FEType type) const
339 : {
340 435549 : if (!_building_helpers && type == _helper_type)
341 203572 : _user_added_fe_face_neighbor_of_helper_type = true;
342 :
343 435549 : if (!_fe_shape_data_face_neighbor[type])
344 103187 : _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 1743851 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
348 : {
349 1308302 : if (!_fe_face_neighbor[dim][type])
350 303391 : _fe_face_neighbor[dim][type] = FEGenericBase<Real>::build(dim, type).release();
351 :
352 1308302 : _fe_face_neighbor[dim][type]->get_phi();
353 1308302 : _fe_face_neighbor[dim][type]->get_dphi();
354 1308302 : if (_need_second_derivative_neighbor.count(type))
355 117 : _fe_face_neighbor[dim][type]->get_d2phi();
356 : }
357 435549 : }
358 :
359 : void
360 769828 : Assembly::buildLowerDFE(FEType type) const
361 : {
362 769828 : if (!_building_helpers && type == _helper_type)
363 406604 : _user_added_fe_lower_of_helper_type = true;
364 :
365 769828 : if (!_fe_shape_data_lower[type])
366 98881 : _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 2331375 : for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
372 : {
373 1561547 : if (!_fe_lower[dim][type])
374 193541 : _fe_lower[dim][type] = FEGenericBase<Real>::build(dim, type).release();
375 :
376 1561547 : _fe_lower[dim][type]->get_phi();
377 1561547 : _fe_lower[dim][type]->get_dphi();
378 1561547 : if (_need_second_derivative.count(type))
379 0 : _fe_lower[dim][type]->get_d2phi();
380 : }
381 769828 : }
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 7480 : Assembly::buildVectorLowerDFE(FEType type) const
406 : {
407 7480 : if (!_vector_fe_shape_data_lower[type])
408 1458 : _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 7480 : unsigned int dim = ((type.family == LAGRANGE_VEC) || (type.family == MONOMIAL_VEC)) ? 0 : 2;
414 7480 : const auto ending_dim = cast_int<unsigned int>(_mesh_dimension - 1);
415 7480 : if (ending_dim < dim)
416 1552 : return;
417 17088 : for (; dim <= ending_dim; dim++)
418 : {
419 11160 : if (!_vector_fe_lower[dim][type])
420 2047 : _vector_fe_lower[dim][type] = FEVectorBase::build(dim, type).release();
421 :
422 11160 : _vector_fe_lower[dim][type]->get_phi();
423 11160 : _vector_fe_lower[dim][type]->get_dphi();
424 11160 : 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 582393 : Assembly::buildVectorFE(const FEType type) const
456 : {
457 582393 : if (!_vector_fe_shape_data[type])
458 1458 : _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 582393 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
463 3086 : type.family == L2_RAVIART_THOMAS)
464 579439 : min_dim = 2;
465 : else
466 2954 : 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 1657435 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
471 : {
472 1075042 : if (!_vector_fe[dim][type])
473 3691 : _vector_fe[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
474 :
475 1075042 : _vector_fe[dim][type]->get_phi();
476 1075042 : _vector_fe[dim][type]->get_dphi();
477 1075042 : if (_need_curl.count(type))
478 63886 : _vector_fe[dim][type]->get_curl_phi();
479 1075042 : if (_need_div.count(type))
480 1001408 : _vector_fe[dim][type]->get_div_phi();
481 1075042 : _vector_fe[dim][type]->get_xyz();
482 : }
483 582393 : }
484 :
485 : void
486 4493 : Assembly::buildVectorFaceFE(const FEType type) const
487 : {
488 4493 : if (!_vector_fe_shape_data_face[type])
489 1458 : _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 4493 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
494 2770 : type.family == L2_RAVIART_THOMAS)
495 1789 : min_dim = 2;
496 : else
497 2704 : 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 15250 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
502 : {
503 10757 : if (!_vector_fe_face[dim][type])
504 3691 : _vector_fe_face[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
505 :
506 10757 : _vector_fe_face[dim][type]->get_phi();
507 10757 : _vector_fe_face[dim][type]->get_dphi();
508 10757 : if (_need_curl.count(type))
509 281 : _vector_fe_face[dim][type]->get_curl_phi();
510 10757 : if (_need_face_div.count(type))
511 728 : _vector_fe_face[dim][type]->get_div_phi();
512 : }
513 4493 : }
514 :
515 : void
516 3740 : Assembly::buildVectorNeighborFE(const FEType type) const
517 : {
518 3740 : if (!_vector_fe_shape_data_neighbor[type])
519 1458 : _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 3740 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
524 2770 : type.family == L2_RAVIART_THOMAS)
525 1036 : min_dim = 2;
526 : else
527 2704 : 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 13488 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
532 : {
533 9748 : if (!_vector_fe_neighbor[dim][type])
534 3691 : _vector_fe_neighbor[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
535 :
536 9748 : _vector_fe_neighbor[dim][type]->get_phi();
537 9748 : _vector_fe_neighbor[dim][type]->get_dphi();
538 9748 : if (_need_curl.count(type))
539 0 : _vector_fe_neighbor[dim][type]->get_curl_phi();
540 9748 : if (_need_neighbor_div.count(type))
541 0 : _vector_fe_neighbor[dim][type]->get_div_phi();
542 : }
543 3740 : }
544 :
545 : void
546 4493 : Assembly::buildVectorFaceNeighborFE(const FEType type) const
547 : {
548 4493 : if (!_vector_fe_shape_data_face_neighbor[type])
549 1458 : _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 4493 : if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
554 2770 : type.family == L2_RAVIART_THOMAS)
555 1789 : min_dim = 2;
556 : else
557 2704 : 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 15250 : for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
562 : {
563 10757 : if (!_vector_fe_face_neighbor[dim][type])
564 3691 : _vector_fe_face_neighbor[dim][type] =
565 7382 : FEGenericBase<VectorValue<Real>>::build(dim, type).release();
566 :
567 10757 : _vector_fe_face_neighbor[dim][type]->get_phi();
568 10757 : _vector_fe_face_neighbor[dim][type]->get_dphi();
569 10757 : if (_need_curl.count(type))
570 281 : _vector_fe_face_neighbor[dim][type]->get_curl_phi();
571 10757 : if (_need_face_neighbor_div.count(type))
572 0 : _vector_fe_face_neighbor[dim][type]->get_div_phi();
573 : }
574 4493 : }
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 72055 : 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 72055 : auto & qvec = _qrules[block];
627 72055 : unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
628 72055 : if (qvec.size() != ndims)
629 72055 : qvec.resize(ndims);
630 :
631 281807 : for (unsigned int i = 0; i < qvec.size(); i++)
632 : {
633 209752 : int dim = i;
634 209752 : auto & q = qvec[dim];
635 209752 : q.vol = QBase::build(type, dim, volume_order);
636 209752 : q.vol->allow_rules_with_negative_weights = allow_negative_qweights;
637 209752 : q.face = QBase::build(type, dim - 1, face_order);
638 209752 : q.face->allow_rules_with_negative_weights = allow_negative_qweights;
639 209752 : q.fv_face = QBase::build(QMONOMIAL, dim - 1, CONSTANT);
640 209752 : q.fv_face->allow_rules_with_negative_weights = allow_negative_qweights;
641 209752 : q.neighbor = std::make_unique<ArbitraryQuadrature>(dim - 1, face_order);
642 209752 : q.neighbor->allow_rules_with_negative_weights = allow_negative_qweights;
643 209752 : q.arbitrary_vol = std::make_unique<ArbitraryQuadrature>(dim, order);
644 209752 : q.arbitrary_vol->allow_rules_with_negative_weights = allow_negative_qweights;
645 209752 : q.arbitrary_face = std::make_unique<ArbitraryQuadrature>(dim - 1, face_order);
646 209752 : q.arbitrary_face->allow_rules_with_negative_weights = allow_negative_qweights;
647 : }
648 :
649 72055 : delete _qrule_msm;
650 72055 : _custom_mortar_qrule = false;
651 72055 : _qrule_msm = QBase::build(type, _mesh_dimension - 1, face_order).release();
652 72055 : _qrule_msm->allow_rules_with_negative_weights = allow_negative_qweights;
653 72055 : _fe_msm->attach_quadrature_rule(_qrule_msm);
654 72055 : }
655 :
656 : void
657 227103 : Assembly::setVolumeQRule(QBase * qrule, unsigned int dim)
658 : {
659 227103 : _current_qrule = qrule;
660 :
661 227103 : if (qrule) // Don't set a NULL qrule
662 : {
663 547761 : for (auto & it : _fe[dim])
664 320658 : it.second->attach_quadrature_rule(qrule);
665 231216 : for (auto & it : _vector_fe[dim])
666 4113 : it.second->attach_quadrature_rule(qrule);
667 227103 : if (!_unique_fe_helper.empty())
668 : {
669 : mooseAssert(dim < _unique_fe_helper.size(), "We should not be indexing out of bounds");
670 219 : _unique_fe_helper[dim]->attach_quadrature_rule(qrule);
671 : }
672 : }
673 227103 : }
674 :
675 : void
676 1164478 : Assembly::setFaceQRule(QBase * qrule, unsigned int dim)
677 : {
678 1164478 : _current_qrule_face = qrule;
679 :
680 2912070 : for (auto & it : _fe_face[dim])
681 1747592 : it.second->attach_quadrature_rule(qrule);
682 1165037 : for (auto & it : _vector_fe_face[dim])
683 559 : it.second->attach_quadrature_rule(qrule);
684 1164478 : if (!_unique_fe_face_helper.empty())
685 : {
686 : mooseAssert(dim < _unique_fe_face_helper.size(), "We should not be indexing out of bounds");
687 197 : _unique_fe_face_helper[dim]->attach_quadrature_rule(qrule);
688 : }
689 1164478 : }
690 :
691 : void
692 564417 : Assembly::setLowerQRule(QBase * qrule, unsigned int dim)
693 : {
694 : // The lower-dimensional quadrature rule matches the face quadrature rule
695 564417 : setFaceQRule(qrule, dim);
696 :
697 564417 : _current_qrule_lower = qrule;
698 :
699 1406991 : for (auto & it : _fe_lower[dim])
700 842574 : it.second->attach_quadrature_rule(qrule);
701 564417 : for (auto & it : _vector_fe_lower[dim])
702 0 : it.second->attach_quadrature_rule(qrule);
703 564417 : 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 564417 : }
709 :
710 : void
711 20616347 : Assembly::setNeighborQRule(QBase * qrule, unsigned int dim)
712 : {
713 20616347 : _current_qrule_neighbor = qrule;
714 :
715 61333988 : for (auto & it : _fe_face_neighbor[dim])
716 40717641 : it.second->attach_quadrature_rule(qrule);
717 20642164 : for (auto & it : _vector_fe_face_neighbor[dim])
718 25817 : it.second->attach_quadrature_rule(qrule);
719 20616347 : 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 68724 : _unique_fe_face_neighbor_helper[dim]->attach_quadrature_rule(qrule);
724 : }
725 20616347 : }
726 :
727 : void
728 62946 : Assembly::clearCachedQRules()
729 : {
730 62946 : _current_qrule = nullptr;
731 62946 : _current_qrule_face = nullptr;
732 62946 : _current_qrule_lower = nullptr;
733 62946 : _current_qrule_neighbor = nullptr;
734 62946 : }
735 :
736 : void
737 18 : Assembly::setMortarQRule(Order order)
738 : {
739 18 : 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 18 : }
760 :
761 : void
762 385394707 : Assembly::reinitFE(const Elem * elem)
763 : {
764 385394707 : unsigned int dim = elem->dim();
765 :
766 862893924 : for (const auto & it : _fe[dim])
767 : {
768 477499241 : FEBase & fe = *it.second;
769 477499241 : const FEType & fe_type = it.first;
770 :
771 477499241 : _current_fe[fe_type] = &fe;
772 :
773 477499241 : FEShapeData & fesd = *_fe_shape_data[fe_type];
774 :
775 477499241 : fe.reinit(elem);
776 :
777 477499217 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_phi()));
778 477499217 : fesd._grad_phi.shallowCopy(
779 477499217 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_dphi()));
780 477499217 : if (_need_second_derivative.count(fe_type))
781 57575 : fesd._second_phi.shallowCopy(
782 57575 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_d2phi()));
783 : }
784 390250915 : for (const auto & it : _vector_fe[dim])
785 : {
786 4856232 : FEVectorBase & fe = *it.second;
787 4856232 : const FEType & fe_type = it.first;
788 :
789 4856232 : _current_vector_fe[fe_type] = &fe;
790 :
791 4856232 : VectorFEShapeData & fesd = *_vector_fe_shape_data[fe_type];
792 :
793 4856232 : fe.reinit(elem);
794 :
795 4856232 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_phi()));
796 4856232 : fesd._grad_phi.shallowCopy(
797 4856232 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_dphi()));
798 4856232 : 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 4856232 : if (_need_curl.count(fe_type))
802 2095820 : fesd._curl_phi.shallowCopy(
803 2095820 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_curl_phi()));
804 4856232 : if (_need_div.count(fe_type))
805 824370 : fesd._div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_div_phi()));
806 : }
807 385394683 : if (!_unique_fe_helper.empty())
808 : {
809 : mooseAssert(dim < _unique_fe_helper.size(), "We should be in bounds here");
810 764430 : _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 385394683 : _current_q_points.shallowCopy(
816 385394683 : const_cast<std::vector<Point> &>(_holder_fe_helper[dim]->get_xyz()));
817 385394683 : _current_JxW.shallowCopy(const_cast<std::vector<Real> &>(_holder_fe_helper[dim]->get_JxW()));
818 :
819 385394683 : if (_subproblem.haveADObjects())
820 : {
821 26894196 : auto n_qp = _current_qrule->n_points();
822 26894196 : resizeADMappingObjects(n_qp, dim);
823 26894196 : if (_displaced)
824 : {
825 529872 : const auto & qw = _current_qrule->get_weights();
826 2904422 : for (unsigned int qp = 0; qp != n_qp; qp++)
827 2374550 : computeSinglePointMapAD(elem, qw, qp, _holder_fe_helper[dim]);
828 : }
829 : else
830 : {
831 94174326 : for (unsigned qp = 0; qp < n_qp; ++qp)
832 67810002 : _ad_JxW[qp] = _current_JxW[qp];
833 26364324 : if (_calculate_xyz)
834 56976248 : for (unsigned qp = 0; qp < n_qp; ++qp)
835 44515734 : _ad_q_points[qp] = _current_q_points[qp];
836 : }
837 :
838 68363125 : for (const auto & it : _fe[dim])
839 : {
840 41468929 : FEBase & fe = *it.second;
841 41468929 : auto fe_type = it.first;
842 41468929 : auto num_shapes = FEInterface::n_shape_functions(fe_type, elem);
843 41468929 : auto & grad_phi = _ad_grad_phi_data[fe_type];
844 :
845 41468929 : grad_phi.resize(num_shapes);
846 167046946 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
847 125578017 : grad_phi[i].resize(n_qp);
848 :
849 41468929 : if (_displaced)
850 818574 : computeGradPhiAD(elem, n_qp, grad_phi, &fe);
851 : else
852 : {
853 40650355 : const auto & regular_grad_phi = _fe_shape_data[fe_type]->_grad_phi;
854 163404378 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
855 582567725 : for (unsigned qp = 0; qp < n_qp; ++qp)
856 459813702 : grad_phi[i][qp] = regular_grad_phi[i][qp];
857 : }
858 : }
859 29655532 : for (const auto & it : _vector_fe[dim])
860 : {
861 2761336 : FEVectorBase & fe = *it.second;
862 2761336 : auto fe_type = it.first;
863 2761336 : auto num_shapes = FEInterface::n_shape_functions(fe_type, elem);
864 2761336 : auto & grad_phi = _ad_vector_grad_phi_data[fe_type];
865 :
866 2761336 : grad_phi.resize(num_shapes);
867 16764904 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
868 14003568 : grad_phi[i].resize(n_qp);
869 :
870 2761336 : if (_displaced)
871 0 : computeGradPhiAD(elem, n_qp, grad_phi, &fe);
872 : else
873 : {
874 2761336 : const auto & regular_grad_phi = _vector_fe_shape_data[fe_type]->_grad_phi;
875 16764904 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
876 71678960 : for (unsigned qp = 0; qp < n_qp; ++qp)
877 57675392 : grad_phi[i][qp] = regular_grad_phi[i][qp];
878 : }
879 : }
880 : }
881 :
882 385394683 : auto n = numExtraElemIntegers();
883 389401269 : for (auto i : make_range(n))
884 4006586 : _extra_elem_ids[i] = _current_elem->get_extra_integer(i);
885 385394683 : _extra_elem_ids[n] = _current_elem->subdomain_id();
886 :
887 385394683 : if (_xfem != nullptr)
888 0 : modifyWeightsDueToXFEM(elem);
889 385394683 : }
890 :
891 : template <typename OutputType>
892 : void
893 943591 : 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 943591 : auto dim = elem->dim();
910 943591 : const auto & dphidxi = fe->get_dphidxi();
911 943591 : const auto & dphideta = fe->get_dphideta();
912 943591 : const auto & dphidzeta = fe->get_dphidzeta();
913 943591 : auto num_shapes = grad_phi.size();
914 :
915 943591 : 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 54074 : case 1:
926 : {
927 146121 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
928 296031 : for (unsigned qp = 0; qp < n_qp; ++qp)
929 : {
930 203984 : grad_phi[i][qp].slice(0) = dphidxi[i][qp] * _ad_dxidx_map[qp];
931 203984 : grad_phi[i][qp].slice(1) = dphidxi[i][qp] * _ad_dxidy_map[qp];
932 203984 : grad_phi[i][qp].slice(2) = dphidxi[i][qp] * _ad_dxidz_map[qp];
933 : }
934 54074 : break;
935 : }
936 :
937 889517 : case 2:
938 : {
939 4247110 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
940 19604835 : for (unsigned qp = 0; qp < n_qp; ++qp)
941 : {
942 32494484 : grad_phi[i][qp].slice(0) =
943 32494484 : dphidxi[i][qp] * _ad_dxidx_map[qp] + dphideta[i][qp] * _ad_detadx_map[qp];
944 32494484 : grad_phi[i][qp].slice(1) =
945 32494484 : dphidxi[i][qp] * _ad_dxidy_map[qp] + dphideta[i][qp] * _ad_detady_map[qp];
946 32494484 : grad_phi[i][qp].slice(2) =
947 32494484 : dphidxi[i][qp] * _ad_dxidz_map[qp] + dphideta[i][qp] * _ad_detadz_map[qp];
948 : }
949 889517 : 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 943591 : }
971 :
972 : void
973 29735285 : Assembly::resizeADMappingObjects(unsigned int n_qp, unsigned int dim)
974 : {
975 29735285 : _ad_dxyzdxi_map.resize(n_qp);
976 29735285 : _ad_dxidx_map.resize(n_qp);
977 29735285 : _ad_dxidy_map.resize(n_qp); // 1D element may live in 2D ...
978 29735285 : _ad_dxidz_map.resize(n_qp); // ... or 3D
979 :
980 29735285 : if (dim > 1)
981 : {
982 19020234 : _ad_dxyzdeta_map.resize(n_qp);
983 19020234 : _ad_detadx_map.resize(n_qp);
984 19020234 : _ad_detady_map.resize(n_qp);
985 19020234 : _ad_detadz_map.resize(n_qp);
986 :
987 19020234 : if (dim > 2)
988 : {
989 1124286 : _ad_dxyzdzeta_map.resize(n_qp);
990 1124286 : _ad_dzetadx_map.resize(n_qp);
991 1124286 : _ad_dzetady_map.resize(n_qp);
992 1124286 : _ad_dzetadz_map.resize(n_qp);
993 : }
994 : }
995 :
996 29735285 : _ad_jac.resize(n_qp);
997 29735285 : _ad_JxW.resize(n_qp);
998 29735285 : if (_calculate_xyz)
999 14797885 : _ad_q_points.resize(n_qp);
1000 29735285 : }
1001 :
1002 : void
1003 2479403 : 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 2479403 : auto dim = elem->dim();
1046 2479403 : const auto & elem_nodes = elem->get_nodes();
1047 2479403 : auto num_shapes = FEInterface::n_shape_functions(fe->get_fe_type(), elem);
1048 2479403 : const auto & phi_map = fe->get_fe_map().get_phi_map();
1049 2479403 : const auto & dphidxi_map = fe->get_fe_map().get_dphidxi_map();
1050 2479403 : const auto & dphideta_map = fe->get_fe_map().get_dphideta_map();
1051 2479403 : const auto & dphidzeta_map = fe->get_fe_map().get_dphidzeta_map();
1052 2479403 : const auto sys_num = _sys.number();
1053 : const bool do_derivatives =
1054 2479403 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
1055 :
1056 2479403 : 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 64315 : case 1:
1068 : {
1069 64315 : if (_calculate_xyz)
1070 10897 : _ad_q_points[p].zero();
1071 :
1072 64315 : _ad_dxyzdxi_map[p].zero();
1073 :
1074 204297 : for (std::size_t i = 0; i < num_shapes; i++)
1075 : {
1076 : libmesh_assert(elem_nodes[i]);
1077 139982 : const Node & node = *elem_nodes[i];
1078 139982 : libMesh::VectorValue<ADReal> elem_point = node;
1079 139982 : if (do_derivatives)
1080 59314 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1081 2417 : if (node.n_dofs(sys_num, disp_num))
1082 4834 : Moose::derivInsert(
1083 2417 : elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1084 :
1085 139982 : _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1086 :
1087 139982 : if (_calculate_xyz)
1088 28562 : _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1089 139982 : }
1090 :
1091 64315 : _ad_jac[p] = _ad_dxyzdxi_map[p].norm();
1092 :
1093 64315 : 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 64315 : const auto jacm2 = 1. / _ad_jac[p] / _ad_jac[p];
1108 64315 : _ad_dxidx_map[p] = jacm2 * _ad_dxyzdxi_map[p](0);
1109 64315 : _ad_dxidy_map[p] = jacm2 * _ad_dxyzdxi_map[p](1);
1110 64315 : _ad_dxidz_map[p] = jacm2 * _ad_dxyzdxi_map[p](2);
1111 :
1112 64315 : _ad_JxW[p] = _ad_jac[p] * qw[p];
1113 :
1114 64315 : break;
1115 64315 : }
1116 :
1117 2415088 : case 2:
1118 : {
1119 2415088 : if (_calculate_xyz)
1120 1641366 : _ad_q_points[p].zero();
1121 2415088 : _ad_dxyzdxi_map[p].zero();
1122 2415088 : _ad_dxyzdeta_map[p].zero();
1123 :
1124 14715816 : for (std::size_t i = 0; i < num_shapes; i++)
1125 : {
1126 : libmesh_assert(elem_nodes[i]);
1127 12300728 : const Node & node = *elem_nodes[i];
1128 12300728 : libMesh::VectorValue<ADReal> elem_point = node;
1129 12300728 : if (do_derivatives)
1130 1662013 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1131 125346 : if (node.n_dofs(sys_num, disp_num))
1132 250692 : Moose::derivInsert(
1133 125346 : elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1134 :
1135 12300728 : _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1136 12300728 : _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
1137 :
1138 12300728 : if (_calculate_xyz)
1139 9637714 : _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1140 12300728 : }
1141 :
1142 2415088 : const auto &dx_dxi = _ad_dxyzdxi_map[p](0), &dx_deta = _ad_dxyzdeta_map[p](0),
1143 2415088 : &dy_dxi = _ad_dxyzdxi_map[p](1), &dy_deta = _ad_dxyzdeta_map[p](1),
1144 2415088 : &dz_dxi = _ad_dxyzdxi_map[p](2), &dz_deta = _ad_dxyzdeta_map[p](2);
1145 :
1146 2415088 : const auto g11 = (dx_dxi * dx_dxi + dy_dxi * dy_dxi + dz_dxi * dz_dxi);
1147 :
1148 2415088 : const auto g12 = (dx_dxi * dx_deta + dy_dxi * dy_deta + dz_dxi * dz_deta);
1149 :
1150 2415088 : const auto & g21 = g12;
1151 :
1152 2415088 : const auto g22 = (dx_deta * dx_deta + dy_deta * dy_deta + dz_deta * dz_deta);
1153 :
1154 2415088 : auto det = (g11 * g22 - g12 * g21);
1155 :
1156 2415088 : 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 2415088 : else if (det.value() <= 0.)
1170 0 : det.value() = TOLERANCE * TOLERANCE;
1171 :
1172 2415088 : const auto inv_det = 1. / det;
1173 : using std::sqrt;
1174 2415088 : _ad_jac[p] = sqrt(det);
1175 :
1176 2415088 : _ad_JxW[p] = _ad_jac[p] * qw[p];
1177 :
1178 2415088 : const auto g11inv = g22 * inv_det;
1179 2415088 : const auto g12inv = -g12 * inv_det;
1180 2415088 : const auto g21inv = -g21 * inv_det;
1181 2415088 : const auto g22inv = g11 * inv_det;
1182 :
1183 2415088 : _ad_dxidx_map[p] = g11inv * dx_dxi + g12inv * dx_deta;
1184 2415088 : _ad_dxidy_map[p] = g11inv * dy_dxi + g12inv * dy_deta;
1185 2415088 : _ad_dxidz_map[p] = g11inv * dz_dxi + g12inv * dz_deta;
1186 :
1187 2415088 : _ad_detadx_map[p] = g21inv * dx_dxi + g22inv * dx_deta;
1188 2415088 : _ad_detady_map[p] = g21inv * dy_dxi + g22inv * dy_deta;
1189 2415088 : _ad_detadz_map[p] = g21inv * dz_dxi + g22inv * dz_deta;
1190 :
1191 2415088 : break;
1192 12075440 : }
1193 :
1194 0 : case 3:
1195 : {
1196 0 : if (_calculate_xyz)
1197 0 : _ad_q_points[p].zero();
1198 0 : _ad_dxyzdxi_map[p].zero();
1199 0 : _ad_dxyzdeta_map[p].zero();
1200 0 : _ad_dxyzdzeta_map[p].zero();
1201 :
1202 0 : for (std::size_t i = 0; i < num_shapes; i++)
1203 : {
1204 : libmesh_assert(elem_nodes[i]);
1205 0 : const Node & node = *elem_nodes[i];
1206 0 : libMesh::VectorValue<ADReal> elem_point = node;
1207 0 : if (do_derivatives)
1208 0 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1209 0 : if (node.n_dofs(sys_num, disp_num))
1210 0 : Moose::derivInsert(
1211 0 : elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1212 :
1213 0 : _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1214 0 : _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
1215 0 : _ad_dxyzdzeta_map[p].add_scaled(elem_point, dphidzeta_map[i][p]);
1216 :
1217 0 : if (_calculate_xyz)
1218 0 : _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1219 0 : }
1220 :
1221 0 : const auto &dx_dxi = _ad_dxyzdxi_map[p](0), &dy_dxi = _ad_dxyzdxi_map[p](1),
1222 0 : &dz_dxi = _ad_dxyzdxi_map[p](2), &dx_deta = _ad_dxyzdeta_map[p](0),
1223 0 : &dy_deta = _ad_dxyzdeta_map[p](1), &dz_deta = _ad_dxyzdeta_map[p](2),
1224 0 : &dx_dzeta = _ad_dxyzdzeta_map[p](0), &dy_dzeta = _ad_dxyzdzeta_map[p](1),
1225 0 : &dz_dzeta = _ad_dxyzdzeta_map[p](2);
1226 :
1227 0 : _ad_jac[p] = (dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) +
1228 0 : dy_dxi * (dz_deta * dx_dzeta - dx_deta * dz_dzeta) +
1229 0 : dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta));
1230 :
1231 0 : if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
1232 : {
1233 : static bool failing = false;
1234 0 : if (!failing)
1235 : {
1236 0 : failing = true;
1237 0 : elem->print_info(libMesh::err);
1238 0 : libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
1239 : << p << " in element " << elem->id());
1240 : }
1241 : else
1242 0 : return;
1243 : }
1244 :
1245 0 : _ad_JxW[p] = _ad_jac[p] * qw[p];
1246 :
1247 0 : const auto inv_jac = 1. / _ad_jac[p];
1248 :
1249 0 : _ad_dxidx_map[p] = (dy_deta * dz_dzeta - dz_deta * dy_dzeta) * inv_jac;
1250 0 : _ad_dxidy_map[p] = (dz_deta * dx_dzeta - dx_deta * dz_dzeta) * inv_jac;
1251 0 : _ad_dxidz_map[p] = (dx_deta * dy_dzeta - dy_deta * dx_dzeta) * inv_jac;
1252 :
1253 0 : _ad_detadx_map[p] = (dz_dxi * dy_dzeta - dy_dxi * dz_dzeta) * inv_jac;
1254 0 : _ad_detady_map[p] = (dx_dxi * dz_dzeta - dz_dxi * dx_dzeta) * inv_jac;
1255 0 : _ad_detadz_map[p] = (dy_dxi * dx_dzeta - dx_dxi * dy_dzeta) * inv_jac;
1256 :
1257 0 : _ad_dzetadx_map[p] = (dy_dxi * dz_deta - dz_dxi * dy_deta) * inv_jac;
1258 0 : _ad_dzetady_map[p] = (dz_dxi * dx_deta - dx_dxi * dz_deta) * inv_jac;
1259 0 : _ad_dzetadz_map[p] = (dx_dxi * dy_deta - dy_dxi * dx_deta) * inv_jac;
1260 :
1261 0 : break;
1262 0 : }
1263 :
1264 0 : default:
1265 0 : libmesh_error_msg("Invalid dim = " << dim);
1266 : }
1267 : }
1268 :
1269 : void
1270 9110361 : Assembly::reinitFEFace(const Elem * elem, unsigned int side)
1271 : {
1272 9110361 : unsigned int dim = elem->dim();
1273 :
1274 26373602 : for (const auto & it : _fe_face[dim])
1275 : {
1276 17263241 : FEBase & fe_face = *it.second;
1277 17263241 : const FEType & fe_type = it.first;
1278 17263241 : FEShapeData & fesd = *_fe_shape_data_face[fe_type];
1279 17263241 : fe_face.reinit(elem, side);
1280 17263241 : _current_fe_face[fe_type] = &fe_face;
1281 :
1282 17263241 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face.get_phi()));
1283 17263241 : fesd._grad_phi.shallowCopy(
1284 17263241 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_dphi()));
1285 17263241 : if (_need_second_derivative.count(fe_type))
1286 15136 : fesd._second_phi.shallowCopy(
1287 15136 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
1288 : }
1289 10609564 : for (const auto & it : _vector_fe_face[dim])
1290 : {
1291 1499203 : FEVectorBase & fe_face = *it.second;
1292 1499203 : const FEType & fe_type = it.first;
1293 :
1294 1499203 : _current_vector_fe_face[fe_type] = &fe_face;
1295 :
1296 1499203 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face[fe_type];
1297 :
1298 1499203 : fe_face.reinit(elem, side);
1299 :
1300 1499203 : fesd._phi.shallowCopy(
1301 1499203 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_phi()));
1302 1499203 : fesd._grad_phi.shallowCopy(
1303 1499203 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_dphi()));
1304 1499203 : if (_need_second_derivative.count(fe_type))
1305 0 : fesd._second_phi.shallowCopy(
1306 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
1307 1499203 : if (_need_curl.count(fe_type))
1308 418745 : fesd._curl_phi.shallowCopy(
1309 418745 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
1310 1499203 : if (_need_face_div.count(fe_type))
1311 52224 : fesd._div_phi.shallowCopy(
1312 52224 : const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
1313 : }
1314 9110361 : if (!_unique_fe_face_helper.empty())
1315 : {
1316 : mooseAssert(dim < _unique_fe_face_helper.size(), "We should be in bounds here");
1317 113670 : _unique_fe_face_helper[dim]->reinit(elem, side);
1318 : }
1319 :
1320 : // During that last loop the helper objects will have been reinitialized as well
1321 : // We need to dig out the q_points and JxW from it.
1322 9110361 : _current_q_points_face.shallowCopy(
1323 9110361 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_xyz()));
1324 9110361 : _current_JxW_face.shallowCopy(
1325 9110361 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_JxW()));
1326 9110361 : _current_normals.shallowCopy(
1327 9110361 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_normals()));
1328 :
1329 9110361 : _mapped_normals.resize(_current_normals.size(), Eigen::Map<RealDIMValue>(nullptr));
1330 30337888 : for (unsigned int i = 0; i < _current_normals.size(); i++)
1331 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
1332 21227527 : new (&_mapped_normals[i]) Eigen::Map<RealDIMValue>(const_cast<Real *>(&_current_normals[i](0)));
1333 :
1334 9110361 : if (_calculate_curvatures)
1335 416 : _curvatures.shallowCopy(
1336 416 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_curvatures()));
1337 :
1338 9110361 : computeADFace(*elem, side);
1339 :
1340 9110361 : if (_xfem != nullptr)
1341 0 : modifyFaceWeightsDueToXFEM(elem, side);
1342 :
1343 9110361 : auto n = numExtraElemIntegers();
1344 9254581 : for (auto i : make_range(n))
1345 144220 : _extra_elem_ids[i] = _current_elem->get_extra_integer(i);
1346 9110361 : _extra_elem_ids[n] = _current_elem->subdomain_id();
1347 9110361 : }
1348 :
1349 : void
1350 52499 : Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std::vector<Real> & qw)
1351 : {
1352 : // Important quantities calculated by this method:
1353 : // - _ad_JxW_face
1354 : // - _ad_q_points_face
1355 : // - _ad_normals
1356 : // - _ad_curvatures
1357 :
1358 52499 : const Elem & side_elem = _compute_face_map_side_elem_builder(elem, side);
1359 52499 : const auto dim = elem.dim();
1360 52499 : const auto n_qp = qw.size();
1361 52499 : const auto & dpsidxi_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsidxi();
1362 52499 : const auto & dpsideta_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsideta();
1363 52499 : const auto & psi_map = _holder_fe_face_helper[dim]->get_fe_map().get_psi();
1364 52499 : std::vector<std::vector<Real>> const * d2psidxi2_map = nullptr;
1365 52499 : std::vector<std::vector<Real>> const * d2psidxideta_map = nullptr;
1366 52499 : std::vector<std::vector<Real>> const * d2psideta2_map = nullptr;
1367 52499 : const auto sys_num = _sys.number();
1368 52499 : const bool do_derivatives = ADReal::do_derivatives && sys_num == _subproblem.currentNlSysNum();
1369 :
1370 52499 : if (_calculate_curvatures)
1371 : {
1372 0 : d2psidxi2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxi2();
1373 0 : d2psidxideta_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxideta();
1374 0 : d2psideta2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psideta2();
1375 : }
1376 :
1377 52499 : switch (dim)
1378 : {
1379 279 : case 1:
1380 : {
1381 279 : if (!n_qp)
1382 0 : break;
1383 :
1384 279 : if (side_elem.node_id(0) == elem.node_id(0))
1385 279 : _ad_normals[0] = Point(-1.);
1386 : else
1387 0 : _ad_normals[0] = Point(1.);
1388 :
1389 279 : VectorValue<ADReal> side_point;
1390 279 : if (_calculate_face_xyz)
1391 : {
1392 279 : const Node & node = side_elem.node_ref(0);
1393 279 : side_point = node;
1394 :
1395 279 : if (do_derivatives)
1396 126 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1397 126 : Moose::derivInsert(
1398 63 : side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1399 : }
1400 :
1401 558 : for (const auto p : make_range(n_qp))
1402 : {
1403 279 : if (_calculate_face_xyz)
1404 : {
1405 279 : _ad_q_points_face[p].zero();
1406 279 : _ad_q_points_face[p].add_scaled(side_point, psi_map[0][p]);
1407 : }
1408 :
1409 279 : _ad_normals[p] = _ad_normals[0];
1410 279 : _ad_JxW_face[p] = 1.0 * qw[p];
1411 : }
1412 :
1413 279 : break;
1414 279 : }
1415 :
1416 52220 : case 2:
1417 : {
1418 52220 : _ad_dxyzdxi_map.resize(n_qp);
1419 52220 : if (_calculate_curvatures)
1420 0 : _ad_d2xyzdxi2_map.resize(n_qp);
1421 :
1422 156794 : for (const auto p : make_range(n_qp))
1423 104574 : _ad_dxyzdxi_map[p].zero();
1424 52220 : if (_calculate_face_xyz)
1425 88704 : for (const auto p : make_range(n_qp))
1426 59136 : _ad_q_points_face[p].zero();
1427 52220 : if (_calculate_curvatures)
1428 0 : for (const auto p : make_range(n_qp))
1429 0 : _ad_d2xyzdxi2_map[p].zero();
1430 :
1431 : const auto n_mapping_shape_functions =
1432 52220 : FE<2, LAGRANGE>::n_dofs(&side_elem, side_elem.default_order());
1433 :
1434 181994 : for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1435 : {
1436 129774 : const Node & node = side_elem.node_ref(i);
1437 129774 : VectorValue<ADReal> side_point = node;
1438 :
1439 129774 : if (do_derivatives)
1440 40698 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1441 0 : Moose::derivInsert(
1442 0 : side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1443 :
1444 389724 : for (const auto p : make_range(n_qp))
1445 259950 : _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1446 129774 : if (_calculate_face_xyz)
1447 253008 : for (const auto p : make_range(n_qp))
1448 168672 : _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1449 129774 : if (_calculate_curvatures)
1450 0 : for (const auto p : make_range(n_qp))
1451 0 : _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1452 129774 : }
1453 :
1454 156794 : for (const auto p : make_range(n_qp))
1455 : {
1456 104574 : _ad_normals[p] =
1457 209148 : (VectorValue<ADReal>(_ad_dxyzdxi_map[p](1), -_ad_dxyzdxi_map[p](0), 0.)).unit();
1458 104574 : const auto the_jac = _ad_dxyzdxi_map[p].norm();
1459 104574 : _ad_JxW_face[p] = the_jac * qw[p];
1460 104574 : if (_calculate_curvatures)
1461 : {
1462 0 : const auto numerator = _ad_d2xyzdxi2_map[p] * _ad_normals[p];
1463 0 : const auto denominator = _ad_dxyzdxi_map[p].norm_sq();
1464 : libmesh_assert_not_equal_to(denominator, 0);
1465 0 : _ad_curvatures[p] = numerator / denominator;
1466 0 : }
1467 104574 : }
1468 :
1469 52220 : break;
1470 : }
1471 :
1472 0 : case 3:
1473 : {
1474 0 : _ad_dxyzdxi_map.resize(n_qp);
1475 0 : _ad_dxyzdeta_map.resize(n_qp);
1476 0 : if (_calculate_curvatures)
1477 : {
1478 0 : _ad_d2xyzdxi2_map.resize(n_qp);
1479 0 : _ad_d2xyzdxideta_map.resize(n_qp);
1480 0 : _ad_d2xyzdeta2_map.resize(n_qp);
1481 : }
1482 :
1483 0 : for (const auto p : make_range(n_qp))
1484 : {
1485 0 : _ad_dxyzdxi_map[p].zero();
1486 0 : _ad_dxyzdeta_map[p].zero();
1487 : }
1488 0 : if (_calculate_face_xyz)
1489 0 : for (const auto p : make_range(n_qp))
1490 0 : _ad_q_points_face[p].zero();
1491 0 : if (_calculate_curvatures)
1492 0 : for (const auto p : make_range(n_qp))
1493 : {
1494 0 : _ad_d2xyzdxi2_map[p].zero();
1495 0 : _ad_d2xyzdxideta_map[p].zero();
1496 0 : _ad_d2xyzdeta2_map[p].zero();
1497 : }
1498 :
1499 : const unsigned int n_mapping_shape_functions =
1500 0 : FE<3, LAGRANGE>::n_dofs(&side_elem, side_elem.default_order());
1501 :
1502 0 : for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1503 : {
1504 0 : const Node & node = side_elem.node_ref(i);
1505 0 : VectorValue<ADReal> side_point = node;
1506 :
1507 0 : if (do_derivatives)
1508 0 : for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1509 0 : Moose::derivInsert(
1510 0 : side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1511 :
1512 0 : for (const auto p : make_range(n_qp))
1513 : {
1514 0 : _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1515 0 : _ad_dxyzdeta_map[p].add_scaled(side_point, dpsideta_map[i][p]);
1516 : }
1517 0 : if (_calculate_face_xyz)
1518 0 : for (const auto p : make_range(n_qp))
1519 0 : _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1520 0 : if (_calculate_curvatures)
1521 0 : for (const auto p : make_range(n_qp))
1522 : {
1523 0 : _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1524 0 : _ad_d2xyzdxideta_map[p].add_scaled(side_point, (*d2psidxideta_map)[i][p]);
1525 0 : _ad_d2xyzdeta2_map[p].add_scaled(side_point, (*d2psideta2_map)[i][p]);
1526 : }
1527 0 : }
1528 :
1529 0 : for (const auto p : make_range(n_qp))
1530 : {
1531 0 : _ad_normals[p] = _ad_dxyzdxi_map[p].cross(_ad_dxyzdeta_map[p]).unit();
1532 :
1533 0 : const auto &dxdxi = _ad_dxyzdxi_map[p](0), &dxdeta = _ad_dxyzdeta_map[p](0),
1534 0 : &dydxi = _ad_dxyzdxi_map[p](1), &dydeta = _ad_dxyzdeta_map[p](1),
1535 0 : &dzdxi = _ad_dxyzdxi_map[p](2), &dzdeta = _ad_dxyzdeta_map[p](2);
1536 :
1537 0 : const auto g11 = (dxdxi * dxdxi + dydxi * dydxi + dzdxi * dzdxi);
1538 :
1539 0 : const auto g12 = (dxdxi * dxdeta + dydxi * dydeta + dzdxi * dzdeta);
1540 :
1541 0 : const auto & g21 = g12;
1542 :
1543 0 : const auto g22 = (dxdeta * dxdeta + dydeta * dydeta + dzdeta * dzdeta);
1544 :
1545 : using std::sqrt;
1546 0 : const auto the_jac = sqrt(g11 * g22 - g12 * g21);
1547 :
1548 0 : _ad_JxW_face[p] = the_jac * qw[p];
1549 :
1550 0 : if (_calculate_curvatures)
1551 : {
1552 0 : const auto L = -_ad_d2xyzdxi2_map[p] * _ad_normals[p];
1553 0 : const auto M = -_ad_d2xyzdxideta_map[p] * _ad_normals[p];
1554 0 : const auto N = -_ad_d2xyzdeta2_map[p] * _ad_normals[p];
1555 0 : const auto E = _ad_dxyzdxi_map[p].norm_sq();
1556 0 : const auto F = _ad_dxyzdxi_map[p] * _ad_dxyzdeta_map[p];
1557 0 : const auto G = _ad_dxyzdeta_map[p].norm_sq();
1558 :
1559 0 : const auto numerator = E * N - 2. * F * M + G * L;
1560 0 : const auto denominator = E * G - F * F;
1561 : libmesh_assert_not_equal_to(denominator, 0.);
1562 0 : _ad_curvatures[p] = 0.5 * numerator / denominator;
1563 0 : }
1564 0 : }
1565 :
1566 0 : break;
1567 : }
1568 :
1569 0 : default:
1570 0 : mooseError("Invalid dimension dim = ", dim);
1571 : }
1572 52499 : }
1573 :
1574 : void
1575 3920245 : Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1576 : {
1577 3920245 : unsigned int neighbor_dim = neighbor->dim();
1578 :
1579 : // reinit neighbor face
1580 11556638 : for (const auto & it : _fe_face_neighbor[neighbor_dim])
1581 : {
1582 7636393 : FEBase & fe_face_neighbor = *it.second;
1583 7636393 : FEType fe_type = it.first;
1584 7636393 : FEShapeData & fesd = *_fe_shape_data_face_neighbor[fe_type];
1585 :
1586 7636393 : fe_face_neighbor.reinit(neighbor, &reference_points);
1587 :
1588 7636393 : _current_fe_face_neighbor[fe_type] = &fe_face_neighbor;
1589 :
1590 7636393 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_phi()));
1591 7636393 : fesd._grad_phi.shallowCopy(
1592 7636393 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor.get_dphi()));
1593 7636393 : if (_need_second_derivative_neighbor.count(fe_type))
1594 8640 : fesd._second_phi.shallowCopy(
1595 8640 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
1596 : }
1597 3946062 : for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
1598 : {
1599 25817 : FEVectorBase & fe_face_neighbor = *it.second;
1600 25817 : const FEType & fe_type = it.first;
1601 :
1602 25817 : _current_vector_fe_face_neighbor[fe_type] = &fe_face_neighbor;
1603 :
1604 25817 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face_neighbor[fe_type];
1605 :
1606 25817 : fe_face_neighbor.reinit(neighbor, &reference_points);
1607 :
1608 25817 : fesd._phi.shallowCopy(
1609 25817 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor.get_phi()));
1610 25817 : fesd._grad_phi.shallowCopy(
1611 25817 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_dphi()));
1612 25817 : if (_need_second_derivative.count(fe_type))
1613 0 : fesd._second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
1614 0 : fe_face_neighbor.get_d2phi()));
1615 25817 : if (_need_curl.count(fe_type))
1616 105 : fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
1617 105 : fe_face_neighbor.get_curl_phi()));
1618 25817 : if (_need_face_neighbor_div.count(fe_type))
1619 0 : fesd._div_phi.shallowCopy(
1620 0 : const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
1621 : }
1622 3920245 : if (!_unique_fe_face_neighbor_helper.empty())
1623 : {
1624 : mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
1625 : "We should be in bounds here");
1626 68724 : _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
1627 : }
1628 :
1629 3920245 : _current_q_points_face_neighbor.shallowCopy(
1630 3920245 : const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
1631 3920245 : }
1632 :
1633 : void
1634 19880 : Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1635 : {
1636 19880 : unsigned int neighbor_dim = neighbor->dim();
1637 :
1638 : // reinit neighbor face
1639 39760 : for (const auto & it : _fe_neighbor[neighbor_dim])
1640 : {
1641 19880 : FEBase & fe_neighbor = *it.second;
1642 19880 : FEType fe_type = it.first;
1643 19880 : FEShapeData & fesd = *_fe_shape_data_neighbor[fe_type];
1644 :
1645 19880 : fe_neighbor.reinit(neighbor, &reference_points);
1646 :
1647 19880 : _current_fe_neighbor[fe_type] = &fe_neighbor;
1648 :
1649 19880 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_phi()));
1650 19880 : fesd._grad_phi.shallowCopy(
1651 19880 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_neighbor.get_dphi()));
1652 19880 : if (_need_second_derivative_neighbor.count(fe_type))
1653 0 : fesd._second_phi.shallowCopy(
1654 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_d2phi()));
1655 : }
1656 19880 : for (const auto & it : _vector_fe_neighbor[neighbor_dim])
1657 : {
1658 0 : FEVectorBase & fe_neighbor = *it.second;
1659 0 : const FEType & fe_type = it.first;
1660 :
1661 0 : _current_vector_fe_neighbor[fe_type] = &fe_neighbor;
1662 :
1663 0 : VectorFEShapeData & fesd = *_vector_fe_shape_data_neighbor[fe_type];
1664 :
1665 0 : fe_neighbor.reinit(neighbor, &reference_points);
1666 :
1667 0 : fesd._phi.shallowCopy(
1668 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_phi()));
1669 0 : fesd._grad_phi.shallowCopy(
1670 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_dphi()));
1671 0 : if (_need_second_derivative.count(fe_type))
1672 0 : fesd._second_phi.shallowCopy(
1673 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_neighbor.get_d2phi()));
1674 0 : if (_need_curl.count(fe_type))
1675 0 : fesd._curl_phi.shallowCopy(
1676 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_curl_phi()));
1677 0 : if (_need_neighbor_div.count(fe_type))
1678 0 : fesd._div_phi.shallowCopy(
1679 0 : const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_div_phi()));
1680 : }
1681 19880 : if (!_unique_fe_neighbor_helper.empty())
1682 : {
1683 : mooseAssert(neighbor_dim < _unique_fe_neighbor_helper.size(), "We should be in bounds here");
1684 0 : _unique_fe_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
1685 : }
1686 19880 : }
1687 :
1688 : void
1689 3940125 : Assembly::reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1690 : {
1691 3940125 : unsigned int neighbor_dim = neighbor->dim();
1692 : mooseAssert(_current_neighbor_subdomain_id == neighbor->subdomain_id(),
1693 : "Neighbor subdomain ID has not been correctly set");
1694 :
1695 : ArbitraryQuadrature * neighbor_rule =
1696 3940125 : qrules(neighbor_dim, _current_neighbor_subdomain_id).neighbor.get();
1697 3940125 : neighbor_rule->setPoints(reference_points);
1698 3940125 : setNeighborQRule(neighbor_rule, neighbor_dim);
1699 :
1700 3940125 : _current_neighbor_elem = neighbor;
1701 : mooseAssert(_current_neighbor_subdomain_id == _current_neighbor_elem->subdomain_id(),
1702 : "current neighbor subdomain has been set incorrectly");
1703 :
1704 : // Calculate the volume of the neighbor
1705 3940125 : if (_need_neighbor_elem_volume)
1706 : {
1707 2205452 : unsigned int dim = neighbor->dim();
1708 2205452 : FEBase & fe = *_holder_fe_neighbor_helper[dim];
1709 2205452 : QBase * qrule = qrules(dim).vol.get();
1710 :
1711 2205452 : fe.attach_quadrature_rule(qrule);
1712 2205452 : fe.reinit(neighbor);
1713 :
1714 2205452 : const std::vector<Real> & JxW = fe.get_JxW();
1715 2205452 : MooseArray<Point> q_points;
1716 2205452 : q_points.shallowCopy(const_cast<std::vector<Point> &>(fe.get_xyz()));
1717 :
1718 2205452 : setCoordinateTransformation(qrule, q_points, _coord_neighbor, _current_neighbor_subdomain_id);
1719 :
1720 2205452 : _current_neighbor_volume = 0.;
1721 17540623 : for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1722 15335171 : _current_neighbor_volume += JxW[qp] * _coord_neighbor[qp];
1723 2205452 : }
1724 :
1725 3940125 : auto n = numExtraElemIntegers();
1726 4021205 : for (auto i : make_range(n))
1727 81080 : _neighbor_extra_elem_ids[i] = _current_neighbor_elem->get_extra_integer(i);
1728 3940125 : _neighbor_extra_elem_ids[n] = _current_neighbor_elem->subdomain_id();
1729 3940125 : }
1730 :
1731 : template <typename Points, typename Coords>
1732 : void
1733 412622227 : Assembly::setCoordinateTransformation(const QBase * qrule,
1734 : const Points & q_points,
1735 : Coords & coord,
1736 : SubdomainID sub_id)
1737 : {
1738 :
1739 : mooseAssert(qrule, "The quadrature rule is null in Assembly::setCoordinateTransformation");
1740 412622227 : auto n_points = qrule->n_points();
1741 : mooseAssert(n_points == q_points.size(),
1742 : "The number of points in the quadrature rule doesn't match the number of passed-in "
1743 : "points in Assembly::setCoordinateTransformation");
1744 :
1745 : // Make sure to honor the name of this method and set the _coord_type member because users may
1746 : // make use of the const Moose::CoordinateSystem & coordTransformation() { return _coord_type; }
1747 : // API. MaterialBase for example uses it
1748 412622227 : _coord_type = _subproblem.getCoordSystem(sub_id);
1749 :
1750 412622227 : coord.resize(n_points);
1751 2220289070 : for (unsigned int qp = 0; qp < n_points; qp++)
1752 1807666843 : coordTransformFactor(_subproblem, sub_id, q_points[qp], coord[qp]);
1753 412622227 : }
1754 :
1755 : void
1756 385394683 : Assembly::computeCurrentElemVolume()
1757 : {
1758 385394683 : if (_current_elem_volume_computed)
1759 0 : return;
1760 :
1761 385394683 : setCoordinateTransformation(
1762 385394683 : _current_qrule, _current_q_points, _coord, _current_elem->subdomain_id());
1763 385394683 : if (_calculate_ad_coord)
1764 12786190 : setCoordinateTransformation(
1765 12786190 : _current_qrule, _ad_q_points, _ad_coord, _current_elem->subdomain_id());
1766 :
1767 385394683 : _current_elem_volume = 0.;
1768 2099969563 : for (unsigned int qp = 0; qp < _current_qrule->n_points(); qp++)
1769 1714574880 : _current_elem_volume += _current_JxW[qp] * _coord[qp];
1770 :
1771 385394683 : _current_elem_volume_computed = true;
1772 : }
1773 :
1774 : void
1775 9110361 : Assembly::computeCurrentFaceVolume()
1776 : {
1777 9110361 : if (_current_side_volume_computed)
1778 0 : return;
1779 :
1780 9110361 : setCoordinateTransformation(
1781 9110361 : _current_qrule_face, _current_q_points_face, _coord, _current_elem->subdomain_id());
1782 9110361 : if (_calculate_ad_coord)
1783 1994365 : setCoordinateTransformation(
1784 1994365 : _current_qrule_face, _ad_q_points_face, _ad_coord, _current_elem->subdomain_id());
1785 :
1786 9110361 : _current_side_volume = 0.;
1787 30337888 : for (unsigned int qp = 0; qp < _current_qrule_face->n_points(); qp++)
1788 21227527 : _current_side_volume += _current_JxW_face[qp] * _coord[qp];
1789 :
1790 9110361 : _current_side_volume_computed = true;
1791 : }
1792 :
1793 : void
1794 362544 : Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points)
1795 : {
1796 362544 : _current_elem = elem;
1797 362544 : _current_neighbor_elem = nullptr;
1798 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1799 : "current subdomain has been set incorrectly");
1800 362544 : _current_elem_volume_computed = false;
1801 :
1802 362544 : FEMap::inverse_map(elem->dim(), elem, physical_points, _temp_reference_points);
1803 :
1804 362544 : reinit(elem, _temp_reference_points);
1805 :
1806 : // Save off the physical points
1807 362544 : _current_physical_points = physical_points;
1808 362544 : }
1809 :
1810 : void
1811 385107527 : Assembly::setVolumeQRule(const Elem * const elem)
1812 : {
1813 385107527 : unsigned int elem_dimension = elem->dim();
1814 385107527 : _current_qrule_volume = qrules(elem_dimension).vol.get();
1815 : // Make sure the qrule is the right one
1816 385107527 : if (_current_qrule != _current_qrule_volume)
1817 196509 : setVolumeQRule(_current_qrule_volume, elem_dimension);
1818 385107527 : }
1819 :
1820 : void
1821 385032163 : Assembly::reinit(const Elem * elem)
1822 : {
1823 385032163 : _current_elem = elem;
1824 385032163 : _current_neighbor_elem = nullptr;
1825 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1826 : "current subdomain has been set incorrectly");
1827 385032163 : _current_elem_volume_computed = false;
1828 385032163 : setVolumeQRule(elem);
1829 385032163 : reinitFE(elem);
1830 :
1831 385032139 : computeCurrentElemVolume();
1832 385032139 : }
1833 :
1834 : void
1835 362544 : Assembly::reinit(const Elem * elem, const std::vector<Point> & reference_points)
1836 : {
1837 362544 : _current_elem = elem;
1838 362544 : _current_neighbor_elem = nullptr;
1839 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1840 : "current subdomain has been set incorrectly");
1841 362544 : _current_elem_volume_computed = false;
1842 :
1843 362544 : unsigned int elem_dimension = _current_elem->dim();
1844 :
1845 362544 : _current_qrule_arbitrary = qrules(elem_dimension).arbitrary_vol.get();
1846 :
1847 : // Make sure the qrule is the right one
1848 362544 : if (_current_qrule != _current_qrule_arbitrary)
1849 30594 : setVolumeQRule(_current_qrule_arbitrary, elem_dimension);
1850 :
1851 362544 : _current_qrule_arbitrary->setPoints(reference_points);
1852 :
1853 362544 : reinitFE(elem);
1854 :
1855 362544 : computeCurrentElemVolume();
1856 362544 : }
1857 :
1858 : void
1859 17673279 : Assembly::reinitFVFace(const FaceInfo & fi)
1860 : {
1861 17673279 : _current_elem = &fi.elem();
1862 17673279 : _current_neighbor_elem = fi.neighborPtr();
1863 17673279 : _current_side = fi.elemSideID();
1864 17673279 : _current_neighbor_side = fi.neighborSideID();
1865 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1866 : "current subdomain has been set incorrectly");
1867 :
1868 17673279 : _current_elem_volume_computed = false;
1869 17673279 : _current_side_volume_computed = false;
1870 :
1871 17673279 : prepareResidual();
1872 17673279 : prepareNeighbor();
1873 17673279 : prepareJacobianBlock();
1874 :
1875 17673279 : unsigned int dim = _current_elem->dim();
1876 17673279 : if (_current_qrule_face != qrules(dim).fv_face.get())
1877 : {
1878 4216 : setFaceQRule(qrules(dim).fv_face.get(), dim);
1879 : // The order of the element that is used for initing here doesn't matter since this will just
1880 : // be used for constant monomials (which only need a single integration point)
1881 4216 : if (dim == 3)
1882 11 : _current_qrule_face->init(QUAD4, /* p_level = */ 0, /* simple_type_only = */ true);
1883 : else
1884 4205 : _current_qrule_face->init(EDGE2, /* p_level = */ 0, /* simple_type_only = */ true);
1885 : }
1886 :
1887 17673279 : _current_side_elem = &_current_side_elem_builder(*_current_elem, _current_side);
1888 :
1889 : mooseAssert(_current_qrule_face->n_points() == 1,
1890 : "Our finite volume quadrature rule should always yield a single point");
1891 :
1892 : // We've initialized the reference points. Now we need to compute the physical location of the
1893 : // quadrature points. We do not do any FE initialization so we cannot simply copy over FE
1894 : // results like we do in reinitFEFace. Instead we handle the computation of the physical
1895 : // locations manually
1896 17673279 : _current_q_points_face.resize(1);
1897 17673279 : const auto & ref_points = _current_qrule_face->get_points();
1898 17673279 : const auto & ref_point = ref_points[0];
1899 17673279 : auto physical_point = FEMap::map(_current_side_elem->dim(), _current_side_elem, ref_point);
1900 17673279 : _current_q_points_face[0] = physical_point;
1901 :
1902 17673279 : if (_current_neighbor_elem)
1903 : {
1904 : mooseAssert(_current_neighbor_subdomain_id == _current_neighbor_elem->subdomain_id(),
1905 : "current neighbor subdomain has been set incorrectly");
1906 : // Now handle the neighbor qrule/qpoints
1907 : ArbitraryQuadrature * const neighbor_rule =
1908 16092122 : qrules(_current_neighbor_elem->dim(), _current_neighbor_subdomain_id).neighbor.get();
1909 : // Here we are setting a reference point that is correct for the neighbor *side* element. It
1910 : // would be wrong if this reference point is used for a volumetric FE reinit with the neighbor
1911 16092122 : neighbor_rule->setPoints(ref_points);
1912 16092122 : setNeighborQRule(neighbor_rule, _current_neighbor_elem->dim());
1913 16092122 : _current_q_points_face_neighbor.resize(1);
1914 16092122 : _current_q_points_face_neighbor[0] = std::move(physical_point);
1915 : }
1916 17673279 : }
1917 :
1918 : QBase *
1919 9110361 : Assembly::qruleFace(const Elem * elem, unsigned int side)
1920 : {
1921 18410035 : return qruleFaceHelper<QBase>(elem, side, [](QRules & q) { return q.face.get(); });
1922 : }
1923 :
1924 : ArbitraryQuadrature *
1925 584100 : Assembly::qruleArbitraryFace(const Elem * elem, unsigned int side)
1926 : {
1927 1168200 : return qruleFaceHelper<ArbitraryQuadrature>(
1928 1752300 : elem, side, [](QRules & q) { return q.arbitrary_face.get(); });
1929 : }
1930 :
1931 : void
1932 9110361 : Assembly::setFaceQRule(const Elem * const elem, const unsigned int side)
1933 : {
1934 9110361 : const auto elem_dimension = elem->dim();
1935 : //// Make sure the qrule is the right one
1936 9110361 : auto rule = qruleFace(elem, side);
1937 9110361 : if (_current_qrule_face != rule)
1938 11745 : setFaceQRule(rule, elem_dimension);
1939 9110361 : }
1940 :
1941 : void
1942 9110361 : Assembly::reinit(const Elem * const elem, const unsigned int side)
1943 : {
1944 9110361 : _current_elem = elem;
1945 9110361 : _current_neighbor_elem = nullptr;
1946 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1947 : "current subdomain has been set incorrectly");
1948 9110361 : _current_side = side;
1949 9110361 : _current_elem_volume_computed = false;
1950 9110361 : _current_side_volume_computed = false;
1951 :
1952 9110361 : _current_side_elem = &_current_side_elem_builder(*elem, side);
1953 :
1954 9110361 : setFaceQRule(elem, side);
1955 9110361 : reinitFEFace(elem, side);
1956 :
1957 9110361 : computeCurrentFaceVolume();
1958 9110361 : }
1959 :
1960 : void
1961 0 : Assembly::reinit(const Elem * elem, unsigned int side, const std::vector<Point> & reference_points)
1962 : {
1963 0 : _current_elem = elem;
1964 0 : _current_neighbor_elem = nullptr;
1965 : mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1966 : "current subdomain has been set incorrectly");
1967 0 : _current_side = side;
1968 0 : _current_elem_volume_computed = false;
1969 0 : _current_side_volume_computed = false;
1970 :
1971 0 : unsigned int elem_dimension = _current_elem->dim();
1972 :
1973 0 : _current_qrule_arbitrary_face = qruleArbitraryFace(elem, side);
1974 :
1975 : // Make sure the qrule is the right one
1976 0 : if (_current_qrule_face != _current_qrule_arbitrary_face)
1977 0 : setFaceQRule(_current_qrule_arbitrary_face, elem_dimension);
1978 :
1979 0 : _current_qrule_arbitrary->setPoints(reference_points);
1980 :
1981 0 : _current_side_elem = &_current_side_elem_builder(*elem, side);
1982 :
1983 0 : reinitFEFace(elem, side);
1984 :
1985 0 : computeCurrentFaceVolume();
1986 0 : }
1987 :
1988 : void
1989 100790072 : Assembly::reinit(const Node * node)
1990 : {
1991 100790072 : _current_node = node;
1992 100790072 : _current_neighbor_node = NULL;
1993 100790072 : }
1994 :
1995 : void
1996 3771719 : Assembly::reinitElemAndNeighbor(const Elem * elem,
1997 : unsigned int side,
1998 : const Elem * neighbor,
1999 : unsigned int neighbor_side,
2000 : const std::vector<Point> * neighbor_reference_points)
2001 : {
2002 3771719 : _current_neighbor_side = neighbor_side;
2003 :
2004 3771719 : reinit(elem, side);
2005 :
2006 3771719 : unsigned int neighbor_dim = neighbor->dim();
2007 :
2008 3771719 : if (neighbor_reference_points)
2009 64800 : _current_neighbor_ref_points = *neighbor_reference_points;
2010 : else
2011 3706919 : FEMap::inverse_map(
2012 7413838 : neighbor_dim, neighbor, _current_q_points_face.stdVector(), _current_neighbor_ref_points);
2013 :
2014 3771719 : _current_neighbor_side_elem = &_current_neighbor_side_elem_builder(*neighbor, neighbor_side);
2015 :
2016 3771719 : reinitFEFaceNeighbor(neighbor, _current_neighbor_ref_points);
2017 3771719 : reinitNeighbor(neighbor, _current_neighbor_ref_points);
2018 3771719 : }
2019 :
2020 : void
2021 584100 : Assembly::reinitElemFaceRef(const Elem * elem,
2022 : unsigned int elem_side,
2023 : Real tolerance,
2024 : const std::vector<Point> * const pts,
2025 : const std::vector<Real> * const weights)
2026 : {
2027 584100 : _current_elem = elem;
2028 :
2029 584100 : unsigned int elem_dim = elem->dim();
2030 :
2031 : // Attach the quadrature rules
2032 584100 : if (pts)
2033 : {
2034 584100 : auto face_rule = qruleArbitraryFace(elem, elem_side);
2035 584100 : face_rule->setPoints(*pts);
2036 584100 : setFaceQRule(face_rule, elem_dim);
2037 : }
2038 : else
2039 : {
2040 0 : auto rule = qruleFace(elem, elem_side);
2041 0 : if (_current_qrule_face != rule)
2042 0 : setFaceQRule(rule, elem_dim);
2043 : }
2044 :
2045 : // reinit face
2046 1461224 : for (const auto & it : _fe_face[elem_dim])
2047 : {
2048 877124 : FEBase & fe_face = *it.second;
2049 877124 : FEType fe_type = it.first;
2050 877124 : FEShapeData & fesd = *_fe_shape_data_face[fe_type];
2051 :
2052 877124 : fe_face.reinit(elem, elem_side, tolerance, pts, weights);
2053 :
2054 877124 : _current_fe_face[fe_type] = &fe_face;
2055 :
2056 877124 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face.get_phi()));
2057 877124 : fesd._grad_phi.shallowCopy(
2058 877124 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_face.get_dphi()));
2059 877124 : if (_need_second_derivative_neighbor.count(fe_type))
2060 0 : fesd._second_phi.shallowCopy(
2061 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
2062 : }
2063 584100 : for (const auto & it : _vector_fe_face[elem_dim])
2064 : {
2065 0 : FEVectorBase & fe_face = *it.second;
2066 0 : const FEType & fe_type = it.first;
2067 :
2068 0 : _current_vector_fe_face[fe_type] = &fe_face;
2069 :
2070 0 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face[fe_type];
2071 :
2072 0 : fe_face.reinit(elem, elem_side, tolerance, pts, weights);
2073 :
2074 0 : fesd._phi.shallowCopy(
2075 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_phi()));
2076 0 : fesd._grad_phi.shallowCopy(
2077 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_dphi()));
2078 0 : if (_need_second_derivative.count(fe_type))
2079 0 : fesd._second_phi.shallowCopy(
2080 0 : const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
2081 0 : if (_need_curl.count(fe_type))
2082 0 : fesd._curl_phi.shallowCopy(
2083 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
2084 0 : if (_need_face_div.count(fe_type))
2085 0 : fesd._div_phi.shallowCopy(
2086 0 : const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
2087 : }
2088 584100 : if (!_unique_fe_face_helper.empty())
2089 : {
2090 : mooseAssert(elem_dim < _unique_fe_face_helper.size(), "We should be in bounds here");
2091 0 : _unique_fe_face_helper[elem_dim]->reinit(elem, elem_side, tolerance, pts, weights);
2092 : }
2093 :
2094 : // During that last loop the helper objects will have been reinitialized
2095 584100 : _current_q_points_face.shallowCopy(
2096 584100 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_xyz()));
2097 584100 : _current_normals.shallowCopy(
2098 584100 : const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_normals()));
2099 584100 : _current_tangents.shallowCopy(const_cast<std::vector<std::vector<Point>> &>(
2100 584100 : _holder_fe_face_helper[elem_dim]->get_tangents()));
2101 : // Note that if the user did pass in points and not weights to this method, JxW will be garbage
2102 : // and should not be used
2103 584100 : _current_JxW_face.shallowCopy(
2104 584100 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_JxW()));
2105 584100 : if (_calculate_curvatures)
2106 0 : _curvatures.shallowCopy(
2107 0 : const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_curvatures()));
2108 :
2109 584100 : computeADFace(*elem, elem_side);
2110 584100 : }
2111 :
2112 : void
2113 9694461 : Assembly::computeADFace(const Elem & elem, const unsigned int side)
2114 : {
2115 9694461 : const auto dim = elem.dim();
2116 :
2117 9694461 : if (_subproblem.haveADObjects())
2118 : {
2119 2841089 : auto n_qp = _current_qrule_face->n_points();
2120 2841089 : resizeADMappingObjects(n_qp, dim);
2121 2841089 : _ad_normals.resize(n_qp);
2122 2841089 : _ad_JxW_face.resize(n_qp);
2123 2841089 : if (_calculate_face_xyz)
2124 2011695 : _ad_q_points_face.resize(n_qp);
2125 2841089 : if (_calculate_curvatures)
2126 416 : _ad_curvatures.resize(n_qp);
2127 :
2128 2841089 : if (_displaced)
2129 : {
2130 52499 : const auto & qw = _current_qrule_face->get_weights();
2131 52499 : computeFaceMap(elem, side, qw);
2132 52499 : const std::vector<Real> dummy_qw(n_qp, 1.);
2133 :
2134 157352 : for (unsigned int qp = 0; qp != n_qp; qp++)
2135 104853 : computeSinglePointMapAD(&elem, dummy_qw, qp, _holder_fe_face_helper[dim]);
2136 52499 : }
2137 : else
2138 : {
2139 10581383 : for (unsigned qp = 0; qp < n_qp; ++qp)
2140 : {
2141 7792793 : _ad_JxW_face[qp] = _current_JxW_face[qp];
2142 7792793 : _ad_normals[qp] = _current_normals[qp];
2143 : }
2144 2788590 : if (_calculate_face_xyz)
2145 6191912 : for (unsigned qp = 0; qp < n_qp; ++qp)
2146 4210064 : _ad_q_points_face[qp] = _current_q_points_face[qp];
2147 2788590 : if (_calculate_curvatures)
2148 832 : for (unsigned qp = 0; qp < n_qp; ++qp)
2149 416 : _ad_curvatures[qp] = _curvatures[qp];
2150 : }
2151 :
2152 7526774 : for (const auto & it : _fe_face[dim])
2153 : {
2154 4685685 : FEBase & fe = *it.second;
2155 4685685 : auto fe_type = it.first;
2156 4685685 : auto num_shapes = FEInterface::n_shape_functions(fe_type, &elem);
2157 4685685 : auto & grad_phi = _ad_grad_phi_data_face[fe_type];
2158 :
2159 4685685 : grad_phi.resize(num_shapes);
2160 36243049 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2161 31557364 : grad_phi[i].resize(n_qp);
2162 :
2163 4685685 : const auto & regular_grad_phi = _fe_shape_data_face[fe_type]->_grad_phi;
2164 :
2165 4685685 : if (_displaced)
2166 125017 : computeGradPhiAD(&elem, n_qp, grad_phi, &fe);
2167 : else
2168 35492386 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2169 145696693 : for (unsigned qp = 0; qp < n_qp; ++qp)
2170 114764975 : grad_phi[i][qp] = regular_grad_phi[i][qp];
2171 : }
2172 3269427 : for (const auto & it : _vector_fe_face[dim])
2173 : {
2174 428338 : FEVectorBase & fe = *it.second;
2175 428338 : auto fe_type = it.first;
2176 428338 : auto num_shapes = FEInterface::n_shape_functions(fe_type, &elem);
2177 428338 : auto & grad_phi = _ad_vector_grad_phi_data_face[fe_type];
2178 :
2179 428338 : grad_phi.resize(num_shapes);
2180 2255050 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2181 1826712 : grad_phi[i].resize(n_qp);
2182 :
2183 428338 : const auto & regular_grad_phi = _vector_fe_shape_data_face[fe_type]->_grad_phi;
2184 :
2185 428338 : if (_displaced)
2186 0 : computeGradPhiAD(&elem, n_qp, grad_phi, &fe);
2187 : else
2188 2255050 : for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2189 5589456 : for (unsigned qp = 0; qp < n_qp; ++qp)
2190 3762744 : grad_phi[i][qp] = regular_grad_phi[i][qp];
2191 : }
2192 : }
2193 9694461 : }
2194 :
2195 : void
2196 584100 : Assembly::reinitNeighborFaceRef(const Elem * neighbor,
2197 : unsigned int neighbor_side,
2198 : Real tolerance,
2199 : const std::vector<Point> * const pts,
2200 : const std::vector<Real> * const weights)
2201 : {
2202 584100 : _current_neighbor_elem = neighbor;
2203 :
2204 584100 : unsigned int neighbor_dim = neighbor->dim();
2205 :
2206 : ArbitraryQuadrature * neighbor_rule =
2207 584100 : qrules(neighbor_dim, neighbor->subdomain_id()).neighbor.get();
2208 584100 : neighbor_rule->setPoints(*pts);
2209 :
2210 : // Attach this quadrature rule to all the _fe_face_neighbor FE objects. This
2211 : // has to have garbage quadrature weights but that's ok because we never
2212 : // actually use the JxW coming from these FE reinit'd objects, e.g. we use the
2213 : // JxW coming from the element face reinit for DGKernels or we use the JxW
2214 : // coming from reinit of the mortar segment element in the case of mortar
2215 584100 : setNeighborQRule(neighbor_rule, neighbor_dim);
2216 :
2217 : // reinit neighbor face
2218 1461224 : for (const auto & it : _fe_face_neighbor[neighbor_dim])
2219 : {
2220 877124 : FEBase & fe_face_neighbor = *it.second;
2221 877124 : FEType fe_type = it.first;
2222 877124 : FEShapeData & fesd = *_fe_shape_data_face_neighbor[fe_type];
2223 :
2224 877124 : fe_face_neighbor.reinit(neighbor, neighbor_side, tolerance, pts, weights);
2225 :
2226 877124 : _current_fe_face_neighbor[fe_type] = &fe_face_neighbor;
2227 :
2228 877124 : fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_phi()));
2229 877124 : fesd._grad_phi.shallowCopy(
2230 877124 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor.get_dphi()));
2231 877124 : if (_need_second_derivative_neighbor.count(fe_type))
2232 0 : fesd._second_phi.shallowCopy(
2233 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
2234 : }
2235 584100 : for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
2236 : {
2237 0 : FEVectorBase & fe_face_neighbor = *it.second;
2238 0 : const FEType & fe_type = it.first;
2239 :
2240 0 : _current_vector_fe_face_neighbor[fe_type] = &fe_face_neighbor;
2241 :
2242 0 : VectorFEShapeData & fesd = *_vector_fe_shape_data_face_neighbor[fe_type];
2243 :
2244 0 : fe_face_neighbor.reinit(neighbor, neighbor_side, tolerance, pts, weights);
2245 :
2246 0 : fesd._phi.shallowCopy(
2247 0 : const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor.get_phi()));
2248 0 : fesd._grad_phi.shallowCopy(
2249 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_dphi()));
2250 0 : if (_need_second_derivative.count(fe_type))
2251 0 : fesd._second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
2252 0 : fe_face_neighbor.get_d2phi()));
2253 0 : if (_need_curl.count(fe_type))
2254 0 : fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
2255 0 : fe_face_neighbor.get_curl_phi()));
2256 0 : if (_need_face_neighbor_div.count(fe_type))
2257 0 : fesd._div_phi.shallowCopy(
2258 0 : const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
2259 : }
2260 584100 : if (!_unique_fe_face_neighbor_helper.empty())
2261 : {
2262 : mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
2263 : "We should be in bounds here");
2264 0 : _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(
2265 : neighbor, neighbor_side, tolerance, pts, weights);
2266 : }
2267 : // During that last loop the helper objects will have been reinitialized as well
2268 : // We need to dig out the q_points from it
2269 584100 : _current_q_points_face_neighbor.shallowCopy(
2270 584100 : const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
2271 584100 : }
2272 :
2273 : void
2274 4057 : Assembly::reinitDual(const Elem * elem,
2275 : const std::vector<Point> & pts,
2276 : const std::vector<Real> & JxW)
2277 : {
2278 4057 : const unsigned int elem_dim = elem->dim();
2279 : mooseAssert(elem_dim == _mesh_dimension - 1,
2280 : "Dual shape functions should only be computed on lower dimensional face elements");
2281 :
2282 8390 : for (const auto & it : _fe_lower[elem_dim])
2283 : {
2284 4333 : FEBase & fe_lower = *it.second;
2285 : // We use customized quadrature rule for integration along the mortar segment elements
2286 4333 : fe_lower.set_calculate_default_dual_coeff(false);
2287 4333 : fe_lower.reinit_dual_shape_coeffs(elem, pts, JxW);
2288 : }
2289 4057 : }
2290 :
2291 : void
2292 600003 : Assembly::reinitLowerDElem(const Elem * elem,
2293 : const std::vector<Point> * const pts,
2294 : const std::vector<Real> * const weights)
2295 : {
2296 600003 : _current_lower_d_elem = elem;
2297 :
2298 600003 : const unsigned int elem_dim = elem->dim();
2299 : mooseAssert(elem_dim < _mesh_dimension,
2300 : "The lower dimensional element should truly be a lower dimensional element");
2301 :
2302 600003 : if (pts)
2303 : {
2304 : // Lower rule matches the face rule for the higher dimensional element
2305 564004 : ArbitraryQuadrature * lower_rule = qrules(elem_dim + 1).arbitrary_face.get();
2306 :
2307 : // This also sets the quadrature weights to unity
2308 564004 : lower_rule->setPoints(*pts);
2309 :
2310 564004 : if (weights)
2311 0 : lower_rule->setWeights(*weights);
2312 :
2313 564004 : setLowerQRule(lower_rule, elem_dim);
2314 : }
2315 35999 : else if (_current_qrule_lower != qrules(elem_dim + 1).face.get())
2316 413 : setLowerQRule(qrules(elem_dim + 1).face.get(), elem_dim);
2317 :
2318 1527006 : for (const auto & it : _fe_lower[elem_dim])
2319 : {
2320 927003 : FEBase & fe_lower = *it.second;
2321 927003 : FEType fe_type = it.first;
2322 :
2323 927003 : fe_lower.reinit(elem);
2324 :
2325 927003 : if (FEShapeData * fesd = _fe_shape_data_lower[fe_type].get())
2326 : {
2327 927003 : fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower.get_phi()));
2328 927003 : fesd->_grad_phi.shallowCopy(
2329 927003 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower.get_dphi()));
2330 927003 : if (_need_second_derivative_neighbor.count(fe_type))
2331 0 : fesd->_second_phi.shallowCopy(
2332 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_d2phi()));
2333 : }
2334 :
2335 : // Dual shape functions need to be computed after primal basis being initialized
2336 927003 : if (FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type].get())
2337 : {
2338 12142 : fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower.get_dual_phi()));
2339 12142 : fesd->_grad_phi.shallowCopy(
2340 12142 : const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower.get_dual_dphi()));
2341 12142 : if (_need_second_derivative_neighbor.count(fe_type))
2342 0 : fesd->_second_phi.shallowCopy(
2343 0 : const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_dual_d2phi()));
2344 : }
2345 : }
2346 600003 : if (!_unique_fe_lower_helper.empty())
2347 : {
2348 : mooseAssert(elem_dim < _unique_fe_lower_helper.size(), "We should be in bounds here");
2349 0 : _unique_fe_lower_helper[elem_dim]->reinit(elem);
2350 : }
2351 :
2352 600003 : if (!_need_lower_d_elem_volume)
2353 190241 : return;
2354 :
2355 409762 : if (pts && !weights)
2356 : {
2357 : // We only have dummy weights so the JxWs computed during our FE reinits are meaningless and
2358 : // we cannot use them
2359 :
2360 813188 : if (_subproblem.getCoordSystem(elem->subdomain_id()) == Moose::CoordinateSystemType::COORD_XYZ)
2361 : // We are in a Cartesian coordinate system and we can just use the element volume method
2362 : // which has fast computation for certain element types
2363 406594 : _current_lower_d_elem_volume = elem->volume();
2364 : else
2365 : // We manually compute the volume taking the curvilinear coordinate transformations into
2366 : // account
2367 0 : _current_lower_d_elem_volume = elementVolume(elem);
2368 : }
2369 : else
2370 : {
2371 : // During that last loop the helper objects will have been reinitialized as well
2372 3168 : FEBase & helper_fe = *_holder_fe_lower_helper[elem_dim];
2373 3168 : const auto & physical_q_points = helper_fe.get_xyz();
2374 3168 : const auto & JxW = helper_fe.get_JxW();
2375 3168 : MooseArray<Real> coord;
2376 3168 : setCoordinateTransformation(
2377 3168 : _current_qrule_lower, physical_q_points, coord, elem->subdomain_id());
2378 3168 : _current_lower_d_elem_volume = 0;
2379 46944 : for (const auto qp : make_range(_current_qrule_lower->n_points()))
2380 43776 : _current_lower_d_elem_volume += JxW[qp] * coord[qp];
2381 3168 : }
2382 : }
2383 :
2384 : void
2385 564004 : Assembly::reinitNeighborLowerDElem(const Elem * elem)
2386 : {
2387 : mooseAssert(elem->dim() < _mesh_dimension,
2388 : "You should be calling reinitNeighborLowerDElem on a lower dimensional element");
2389 :
2390 564004 : _current_neighbor_lower_d_elem = elem;
2391 :
2392 564004 : if (!_need_neighbor_lower_d_elem_volume)
2393 157410 : return;
2394 :
2395 406594 : if (_subproblem.getCoordSystem(elem->subdomain_id()) == Moose::CoordinateSystemType::COORD_XYZ)
2396 : // We are in a Cartesian coordinate system and we can just use the element volume method which
2397 : // has fast computation for certain element types
2398 406594 : _current_neighbor_lower_d_elem_volume = elem->volume();
2399 : else
2400 : // We manually compute the volume taking the curvilinear coordinate transformations into
2401 : // account
2402 0 : _current_neighbor_lower_d_elem_volume = elementVolume(elem);
2403 : }
2404 :
2405 : void
2406 1128008 : Assembly::reinitMortarElem(const Elem * elem)
2407 : {
2408 : mooseAssert(elem->dim() == _mesh_dimension - 1,
2409 : "You should be calling reinitMortarElem on a lower dimensional element");
2410 :
2411 1128008 : _fe_msm->reinit(elem);
2412 1128008 : _msm_elem = elem;
2413 :
2414 1128008 : MooseArray<Point> array_q_points;
2415 1128008 : array_q_points.shallowCopy(const_cast<std::vector<Point> &>(_fe_msm->get_xyz()));
2416 1128008 : setCoordinateTransformation(_qrule_msm, array_q_points, _coord_msm, elem->subdomain_id());
2417 1128008 : }
2418 :
2419 : void
2420 148526 : Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
2421 : unsigned int neighbor_side,
2422 : const std::vector<Point> & physical_points)
2423 : {
2424 148526 : unsigned int neighbor_dim = neighbor->dim();
2425 148526 : FEMap::inverse_map(neighbor_dim, neighbor, physical_points, _current_neighbor_ref_points);
2426 :
2427 148526 : if (_need_JxW_neighbor)
2428 : {
2429 : mooseAssert(
2430 : physical_points.size() == 1,
2431 : "If reinitializing with more than one point, then I am dubious of your use case. Perhaps "
2432 : "you are performing a DG type method and you are reinitializing using points from the "
2433 : "element face. In such a case your neighbor JxW must have its index order 'match' the "
2434 : "element JxW index order, e.g. imagining a vertical 1D face with two quadrature points, "
2435 : "if "
2436 : "index 0 for elem JxW corresponds to the 'top' quadrature point, then index 0 for "
2437 : "neighbor "
2438 : "JxW must also correspond to the 'top' quadrature point. And libMesh/MOOSE has no way to "
2439 : "guarantee that with multiple quadrature points.");
2440 :
2441 97230 : _current_neighbor_side_elem = &_current_neighbor_side_elem_builder(*neighbor, neighbor_side);
2442 :
2443 : // With a single point our size-1 JxW should just be the element volume
2444 97230 : _current_JxW_neighbor.resize(1);
2445 97230 : _current_JxW_neighbor[0] = _current_neighbor_side_elem->volume();
2446 : }
2447 :
2448 148526 : reinitFEFaceNeighbor(neighbor, _current_neighbor_ref_points);
2449 148526 : reinitNeighbor(neighbor, _current_neighbor_ref_points);
2450 :
2451 : // Save off the physical points
2452 148526 : _current_physical_points = physical_points;
2453 148526 : }
2454 :
2455 : void
2456 19880 : Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
2457 : const std::vector<Point> & physical_points)
2458 : {
2459 19880 : unsigned int neighbor_dim = neighbor->dim();
2460 19880 : FEMap::inverse_map(neighbor_dim, neighbor, physical_points, _current_neighbor_ref_points);
2461 :
2462 19880 : reinitFENeighbor(neighbor, _current_neighbor_ref_points);
2463 19880 : reinitNeighbor(neighbor, _current_neighbor_ref_points);
2464 : // Save off the physical points
2465 19880 : _current_physical_points = physical_points;
2466 19880 : }
2467 :
2468 : void
2469 68960 : Assembly::init(const CouplingMatrix * cm)
2470 : {
2471 68960 : _cm = cm;
2472 :
2473 68960 : unsigned int n_vars = _sys.nVariables();
2474 :
2475 68960 : _cm_ss_entry.clear();
2476 68960 : _cm_sf_entry.clear();
2477 68960 : _cm_fs_entry.clear();
2478 68960 : _cm_ff_entry.clear();
2479 :
2480 68960 : auto & vars = _sys.getVariables(_tid);
2481 :
2482 68960 : _block_diagonal_matrix = true;
2483 136006 : for (auto & ivar : vars)
2484 : {
2485 67046 : auto i = ivar->number();
2486 67046 : if (i >= _component_block_diagonal.size())
2487 67046 : _component_block_diagonal.resize(i + 1, true);
2488 :
2489 67046 : auto ivar_start = _cm_ff_entry.size();
2490 135793 : for (unsigned int k = 0; k < ivar->count(); ++k)
2491 : {
2492 68747 : unsigned int iv = i + k;
2493 160387 : for (const auto & j : ConstCouplingRow(iv, *_cm))
2494 : {
2495 91640 : if (_sys.isScalarVariable(j))
2496 : {
2497 1151 : auto & jvar = _sys.getScalarVariable(_tid, j);
2498 1151 : _cm_fs_entry.push_back(std::make_pair(ivar, &jvar));
2499 1151 : _block_diagonal_matrix = false;
2500 : }
2501 : else
2502 : {
2503 90489 : auto & jvar = _sys.getVariable(_tid, j);
2504 90489 : auto pair = std::make_pair(ivar, &jvar);
2505 90489 : auto c = ivar_start;
2506 : // check if the pair has been pushed or not
2507 90489 : bool has_pair = false;
2508 117587 : for (; c < _cm_ff_entry.size(); ++c)
2509 32935 : if (_cm_ff_entry[c] == pair)
2510 : {
2511 5837 : has_pair = true;
2512 5837 : break;
2513 : }
2514 90489 : if (!has_pair)
2515 84652 : _cm_ff_entry.push_back(pair);
2516 : // only set having diagonal matrix to false when ivar and jvar numbers are different
2517 : // Note: for array variables, since we save the entire local Jacobian of all components,
2518 : // even there are couplings among components of the same array variable, we still
2519 : // do not set the flag to false.
2520 90489 : if (i != jvar.number())
2521 20036 : _block_diagonal_matrix = false;
2522 70453 : else if (iv != j)
2523 1706 : _component_block_diagonal[i] = false;
2524 : }
2525 : }
2526 : }
2527 : }
2528 :
2529 68960 : auto & scalar_vars = _sys.getScalarVariables(_tid);
2530 :
2531 70407 : for (auto & ivar : scalar_vars)
2532 : {
2533 1447 : auto i = ivar->number();
2534 1447 : if (i >= _component_block_diagonal.size())
2535 1383 : _component_block_diagonal.resize(i + 1, true);
2536 :
2537 4261 : for (const auto & j : ConstCouplingRow(i, *_cm))
2538 2814 : if (_sys.isScalarVariable(j))
2539 : {
2540 1663 : auto & jvar = _sys.getScalarVariable(_tid, j);
2541 1663 : _cm_ss_entry.push_back(std::make_pair(ivar, &jvar));
2542 : }
2543 : else
2544 : {
2545 1151 : auto & jvar = _sys.getVariable(_tid, j);
2546 1151 : _cm_sf_entry.push_back(std::make_pair(ivar, &jvar));
2547 : }
2548 : }
2549 :
2550 68960 : if (_block_diagonal_matrix && scalar_vars.size() != 0)
2551 409 : _block_diagonal_matrix = false;
2552 :
2553 68960 : auto num_vector_tags = _residual_vector_tags.size();
2554 :
2555 68960 : _sub_Re.resize(num_vector_tags);
2556 68960 : _sub_Rn.resize(num_vector_tags);
2557 68960 : _sub_Rl.resize(num_vector_tags);
2558 244924 : for (MooseIndex(_sub_Re) i = 0; i < _sub_Re.size(); i++)
2559 : {
2560 175964 : _sub_Re[i].resize(n_vars);
2561 175964 : _sub_Rn[i].resize(n_vars);
2562 175964 : _sub_Rl[i].resize(n_vars);
2563 : }
2564 :
2565 68960 : _cached_residual_values.resize(num_vector_tags);
2566 68960 : _cached_residual_rows.resize(num_vector_tags);
2567 :
2568 68960 : auto num_matrix_tags = _subproblem.numMatrixTags();
2569 :
2570 68960 : _cached_jacobian_values.resize(num_matrix_tags);
2571 68960 : _cached_jacobian_rows.resize(num_matrix_tags);
2572 68960 : _cached_jacobian_cols.resize(num_matrix_tags);
2573 :
2574 : // Element matrices
2575 68960 : _sub_Kee.resize(num_matrix_tags);
2576 68960 : _sub_Keg.resize(num_matrix_tags);
2577 68960 : _sub_Ken.resize(num_matrix_tags);
2578 68960 : _sub_Kne.resize(num_matrix_tags);
2579 68960 : _sub_Knn.resize(num_matrix_tags);
2580 68960 : _sub_Kll.resize(num_matrix_tags);
2581 68960 : _sub_Kle.resize(num_matrix_tags);
2582 68960 : _sub_Kln.resize(num_matrix_tags);
2583 68960 : _sub_Kel.resize(num_matrix_tags);
2584 68960 : _sub_Knl.resize(num_matrix_tags);
2585 :
2586 68960 : _jacobian_block_used.resize(num_matrix_tags);
2587 68960 : _jacobian_block_neighbor_used.resize(num_matrix_tags);
2588 68960 : _jacobian_block_lower_used.resize(num_matrix_tags);
2589 68960 : _jacobian_block_nonlocal_used.resize(num_matrix_tags);
2590 :
2591 209340 : for (MooseIndex(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
2592 : {
2593 140380 : _sub_Keg[tag].resize(n_vars);
2594 140380 : _sub_Ken[tag].resize(n_vars);
2595 140380 : _sub_Kne[tag].resize(n_vars);
2596 140380 : _sub_Knn[tag].resize(n_vars);
2597 140380 : _sub_Kee[tag].resize(n_vars);
2598 140380 : _sub_Kll[tag].resize(n_vars);
2599 140380 : _sub_Kle[tag].resize(n_vars);
2600 140380 : _sub_Kln[tag].resize(n_vars);
2601 140380 : _sub_Kel[tag].resize(n_vars);
2602 140380 : _sub_Knl[tag].resize(n_vars);
2603 :
2604 140380 : _jacobian_block_used[tag].resize(n_vars);
2605 140380 : _jacobian_block_neighbor_used[tag].resize(n_vars);
2606 140380 : _jacobian_block_lower_used[tag].resize(n_vars);
2607 140380 : _jacobian_block_nonlocal_used[tag].resize(n_vars);
2608 284495 : for (MooseIndex(n_vars) i = 0; i < n_vars; ++i)
2609 : {
2610 144115 : if (!_block_diagonal_matrix)
2611 : {
2612 32700 : _sub_Kee[tag][i].resize(n_vars);
2613 32700 : _sub_Keg[tag][i].resize(n_vars);
2614 32700 : _sub_Ken[tag][i].resize(n_vars);
2615 32700 : _sub_Kne[tag][i].resize(n_vars);
2616 32700 : _sub_Knn[tag][i].resize(n_vars);
2617 32700 : _sub_Kll[tag][i].resize(n_vars);
2618 32700 : _sub_Kle[tag][i].resize(n_vars);
2619 32700 : _sub_Kln[tag][i].resize(n_vars);
2620 32700 : _sub_Kel[tag][i].resize(n_vars);
2621 32700 : _sub_Knl[tag][i].resize(n_vars);
2622 :
2623 32700 : _jacobian_block_used[tag][i].resize(n_vars);
2624 32700 : _jacobian_block_neighbor_used[tag][i].resize(n_vars);
2625 32700 : _jacobian_block_lower_used[tag][i].resize(n_vars);
2626 32700 : _jacobian_block_nonlocal_used[tag][i].resize(n_vars);
2627 : }
2628 : else
2629 : {
2630 111415 : _sub_Kee[tag][i].resize(1);
2631 111415 : _sub_Keg[tag][i].resize(1);
2632 111415 : _sub_Ken[tag][i].resize(1);
2633 111415 : _sub_Kne[tag][i].resize(1);
2634 111415 : _sub_Knn[tag][i].resize(1);
2635 111415 : _sub_Kll[tag][i].resize(1);
2636 111415 : _sub_Kle[tag][i].resize(1);
2637 111415 : _sub_Kln[tag][i].resize(1);
2638 111415 : _sub_Kel[tag][i].resize(1);
2639 111415 : _sub_Knl[tag][i].resize(1);
2640 :
2641 111415 : _jacobian_block_used[tag][i].resize(1);
2642 111415 : _jacobian_block_neighbor_used[tag][i].resize(1);
2643 111415 : _jacobian_block_lower_used[tag][i].resize(1);
2644 111415 : _jacobian_block_nonlocal_used[tag][i].resize(1);
2645 : }
2646 : }
2647 : }
2648 68960 : }
2649 :
2650 : void
2651 63 : Assembly::initNonlocalCoupling()
2652 : {
2653 63 : _cm_nonlocal_entry.clear();
2654 :
2655 63 : auto & vars = _sys.getVariables(_tid);
2656 :
2657 189 : for (auto & ivar : vars)
2658 : {
2659 126 : auto i = ivar->number();
2660 126 : auto ivar_start = _cm_nonlocal_entry.size();
2661 252 : for (unsigned int k = 0; k < ivar->count(); ++k)
2662 : {
2663 126 : unsigned int iv = i + k;
2664 216 : for (const auto & j : ConstCouplingRow(iv, _nonlocal_cm))
2665 90 : if (!_sys.isScalarVariable(j))
2666 : {
2667 90 : auto & jvar = _sys.getVariable(_tid, j);
2668 90 : auto pair = std::make_pair(ivar, &jvar);
2669 90 : auto c = ivar_start;
2670 : // check if the pair has been pushed or not
2671 90 : bool has_pair = false;
2672 117 : for (; c < _cm_nonlocal_entry.size(); ++c)
2673 27 : if (_cm_nonlocal_entry[c] == pair)
2674 : {
2675 0 : has_pair = true;
2676 0 : break;
2677 : }
2678 90 : if (!has_pair)
2679 90 : _cm_nonlocal_entry.push_back(pair);
2680 : }
2681 : }
2682 : }
2683 63 : }
2684 :
2685 : void
2686 73225995 : Assembly::prepareJacobianBlock()
2687 : {
2688 179459825 : for (const auto & it : _cm_ff_entry)
2689 : {
2690 106233830 : MooseVariableFEBase & ivar = *(it.first);
2691 106233830 : MooseVariableFEBase & jvar = *(it.second);
2692 :
2693 106233830 : unsigned int vi = ivar.number();
2694 106233830 : unsigned int vj = jvar.number();
2695 :
2696 106233830 : const bool array_block_diagonal_purely_diagonal = vi == vj && _component_block_diagonal[vi];
2697 106233830 : auto num_cols = jvar.dofIndices().size();
2698 106233830 : if (array_block_diagonal_purely_diagonal)
2699 84043802 : num_cols /= jvar.count();
2700 :
2701 321343238 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2702 : {
2703 215109408 : jacobianBlock(vi, vj, LocalDataKey{}, tag).resize(ivar.dofIndices().size(), num_cols);
2704 215109408 : jacobianBlockUsed(tag, vi, vj, false);
2705 : }
2706 : }
2707 73225995 : }
2708 :
2709 : void
2710 403038943 : Assembly::prepareResidual()
2711 : {
2712 403038943 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2713 840914451 : for (const auto & var : vars)
2714 1641760533 : for (auto & tag_Re : _sub_Re)
2715 1203885025 : tag_Re[var->number()].resize(var->dofIndices().size());
2716 403038943 : }
2717 :
2718 : void
2719 550746 : Assembly::prepare()
2720 : {
2721 550746 : prepareJacobianBlock();
2722 550746 : prepareResidual();
2723 550746 : }
2724 :
2725 : void
2726 8824 : Assembly::prepareNonlocal()
2727 : {
2728 17862 : for (const auto & it : _cm_nonlocal_entry)
2729 : {
2730 9038 : MooseVariableFEBase & ivar = *(it.first);
2731 9038 : MooseVariableFEBase & jvar = *(it.second);
2732 :
2733 9038 : unsigned int vi = ivar.number();
2734 9038 : unsigned int vj = jvar.number();
2735 :
2736 9038 : const bool array_block_diagonal_purely_diagonal = vi == vj && _component_block_diagonal[vi];
2737 9038 : auto num_cols = jvar.allDofIndices().size();
2738 9038 : if (array_block_diagonal_purely_diagonal)
2739 334 : num_cols /= jvar.count();
2740 :
2741 27114 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2742 27114 : tag < _jacobian_block_nonlocal_used.size();
2743 : tag++)
2744 : {
2745 18076 : jacobianBlockNonlocal(vi, vj, LocalDataKey{}, tag).resize(ivar.dofIndices().size(), num_cols);
2746 18076 : jacobianBlockNonlocalUsed(tag, vi, vj, false);
2747 : }
2748 : }
2749 8824 : }
2750 :
2751 : void
2752 1338 : Assembly::prepareVariable(MooseVariableFEBase * var)
2753 : {
2754 5054 : for (const auto & it : _cm_ff_entry)
2755 : {
2756 3716 : MooseVariableFEBase & ivar = *(it.first);
2757 3716 : MooseVariableFEBase & jvar = *(it.second);
2758 :
2759 3716 : unsigned int vi = ivar.number();
2760 3716 : unsigned int vj = jvar.number();
2761 :
2762 3716 : const bool array_block_diagonal_purely_diagonal = vi == vj && _component_block_diagonal[vi];
2763 3716 : auto num_cols = jvar.dofIndices().size();
2764 3716 : if (array_block_diagonal_purely_diagonal)
2765 2444 : num_cols /= jvar.count();
2766 :
2767 3716 : if (vi == var->number() || vj == var->number())
2768 : {
2769 7830 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2770 : {
2771 5220 : jacobianBlock(vi, vj, LocalDataKey{}, tag).resize(ivar.dofIndices().size(), num_cols);
2772 5220 : jacobianBlockUsed(tag, vi, vj, false);
2773 : }
2774 : }
2775 : }
2776 :
2777 4650 : for (auto & tag_Re : _sub_Re)
2778 3312 : tag_Re[var->number()].resize(var->dofIndices().size());
2779 1338 : }
2780 :
2781 : void
2782 0 : Assembly::prepareVariableNonlocal(MooseVariableFEBase * var)
2783 : {
2784 0 : for (const auto & it : _cm_nonlocal_entry)
2785 : {
2786 0 : MooseVariableFEBase & ivar = *(it.first);
2787 0 : MooseVariableFEBase & jvar = *(it.second);
2788 :
2789 0 : unsigned int vi = ivar.number();
2790 0 : unsigned int vj = jvar.number();
2791 :
2792 0 : const bool array_block_diagonal_purely_diagonal = vi == vj && _component_block_diagonal[vi];
2793 0 : auto num_cols = jvar.dofIndices().size();
2794 0 : if (array_block_diagonal_purely_diagonal)
2795 0 : num_cols /= jvar.count();
2796 :
2797 0 : if (vi == var->number() || vj == var->number())
2798 : {
2799 0 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2800 0 : tag < _jacobian_block_nonlocal_used.size();
2801 : tag++)
2802 : {
2803 0 : jacobianBlockNonlocal(vi, vj, LocalDataKey{}, tag)
2804 0 : .resize(ivar.dofIndices().size(), num_cols);
2805 0 : jacobianBlockNonlocalUsed(tag, vi, vj);
2806 : }
2807 : }
2808 : }
2809 0 : }
2810 :
2811 : void
2812 22206913 : Assembly::prepareNeighbor()
2813 : {
2814 59793369 : for (const auto & it : _cm_ff_entry)
2815 : {
2816 37586456 : MooseVariableFEBase & ivar = *(it.first);
2817 37586456 : MooseVariableFEBase & jvar = *(it.second);
2818 :
2819 37586456 : unsigned int vi = ivar.number();
2820 37586456 : unsigned int vj = jvar.number();
2821 :
2822 37586456 : const bool array_block_diagonal_purely_diagonal = vi == vj && _component_block_diagonal[vi];
2823 37586456 : const auto dofs_divisor = array_block_diagonal_purely_diagonal ? jvar.count() : 1;
2824 :
2825 112766856 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
2826 112766856 : tag < _jacobian_block_neighbor_used.size();
2827 : tag++)
2828 : {
2829 75180400 : jacobianBlockNeighbor(Moose::ElementNeighbor, vi, vj, LocalDataKey{}, tag)
2830 75180400 : .resize(ivar.dofIndices().size(), jvar.dofIndicesNeighbor().size() / dofs_divisor);
2831 :
2832 75180400 : jacobianBlockNeighbor(Moose::NeighborElement, vi, vj, LocalDataKey{}, tag)
2833 75180400 : .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndices().size() / dofs_divisor);
2834 :
2835 75180400 : jacobianBlockNeighbor(Moose::NeighborNeighbor, vi, vj, LocalDataKey{}, tag)
2836 75180400 : .resize(ivar.dofIndicesNeighbor().size(),
2837 75180400 : jvar.dofIndicesNeighbor().size() / dofs_divisor);
2838 :
2839 75180400 : jacobianBlockNeighborUsed(tag, vi, vj, false);
2840 : }
2841 : }
2842 :
2843 22206913 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2844 48891739 : for (const auto & var : vars)
2845 93195727 : for (auto & tag_Rn : _sub_Rn)
2846 66510901 : tag_Rn[var->number()].resize(var->dofIndicesNeighbor().size());
2847 22206913 : }
2848 :
2849 : void
2850 600003 : Assembly::prepareLowerD()
2851 : {
2852 2802079 : for (const auto & it : _cm_ff_entry)
2853 : {
2854 2202076 : MooseVariableFEBase & ivar = *(it.first);
2855 2202076 : MooseVariableFEBase & jvar = *(it.second);
2856 :
2857 2202076 : unsigned int vi = ivar.number();
2858 2202076 : unsigned int vj = jvar.number();
2859 :
2860 2202076 : const bool array_block_diagonal_purely_diagonal = vi == vj && _component_block_diagonal[vi];
2861 2202076 : const auto dofs_divisor = array_block_diagonal_purely_diagonal ? jvar.count() : 1;
2862 :
2863 6606228 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
2864 : tag++)
2865 : {
2866 : // To cover all possible cases we should have 9 combinations below for every 2-permutation
2867 : // of Lower,Secondary,Primary. However, 4 cases will in general be covered by calls to
2868 : // prepare() and prepareNeighbor(). These calls will cover SecondarySecondary
2869 : // (ElementElement), SecondaryPrimary (ElementNeighbor), PrimarySecondary (NeighborElement),
2870 : // and PrimaryPrimary (NeighborNeighbor). With these covered we only need to prepare the 5
2871 : // remaining below
2872 :
2873 : // derivatives w.r.t. lower dimensional residuals
2874 4404152 : jacobianBlockMortar(Moose::LowerLower, vi, vj, LocalDataKey{}, tag)
2875 4404152 : .resize(ivar.dofIndicesLower().size(), jvar.dofIndicesLower().size() / dofs_divisor);
2876 :
2877 4404152 : jacobianBlockMortar(Moose::LowerSecondary, vi, vj, LocalDataKey{}, tag)
2878 4404152 : .resize(ivar.dofIndicesLower().size(), jvar.dofIndices().size() / dofs_divisor);
2879 :
2880 4404152 : jacobianBlockMortar(Moose::LowerPrimary, vi, vj, LocalDataKey{}, tag)
2881 4404152 : .resize(ivar.dofIndicesLower().size(), jvar.dofIndicesNeighbor().size() / dofs_divisor);
2882 :
2883 : // derivatives w.r.t. interior secondary residuals
2884 4404152 : jacobianBlockMortar(Moose::SecondaryLower, vi, vj, LocalDataKey{}, tag)
2885 4404152 : .resize(ivar.dofIndices().size(), jvar.dofIndicesLower().size() / dofs_divisor);
2886 :
2887 : // derivatives w.r.t. interior primary residuals
2888 4404152 : jacobianBlockMortar(Moose::PrimaryLower, vi, vj, LocalDataKey{}, tag)
2889 4404152 : .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndicesLower().size() / dofs_divisor);
2890 :
2891 4404152 : jacobianBlockLowerUsed(tag, vi, vj, false);
2892 : }
2893 : }
2894 :
2895 600003 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2896 1709955 : for (const auto & var : vars)
2897 3404396 : for (auto & tag_Rl : _sub_Rl)
2898 2294444 : tag_Rl[var->number()].resize(var->dofIndicesLower().size());
2899 600003 : }
2900 :
2901 : void
2902 0 : Assembly::prepareBlock(unsigned int ivar,
2903 : unsigned int jvar,
2904 : const std::vector<dof_id_type> & dof_indices)
2905 : {
2906 0 : const auto & iv = _sys.getVariable(_tid, ivar);
2907 0 : const auto & jv = _sys.getVariable(_tid, jvar);
2908 0 : const unsigned int ivn = iv.number();
2909 0 : const unsigned int jvn = jv.number();
2910 0 : const unsigned int icount = iv.count();
2911 0 : unsigned int jcount = jv.count();
2912 0 : if (ivn == jvn && _component_block_diagonal[ivn])
2913 0 : jcount = 1;
2914 :
2915 0 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2916 : {
2917 0 : jacobianBlock(ivn, jvn, LocalDataKey{}, tag)
2918 0 : .resize(dof_indices.size() * icount, dof_indices.size() * jcount);
2919 0 : jacobianBlockUsed(tag, ivn, jvn, false);
2920 : }
2921 :
2922 0 : for (auto & tag_Re : _sub_Re)
2923 0 : tag_Re[ivn].resize(dof_indices.size() * icount);
2924 0 : }
2925 :
2926 : void
2927 0 : Assembly::prepareBlockNonlocal(unsigned int ivar,
2928 : unsigned int jvar,
2929 : const std::vector<dof_id_type> & idof_indices,
2930 : const std::vector<dof_id_type> & jdof_indices)
2931 : {
2932 0 : const auto & iv = _sys.getVariable(_tid, ivar);
2933 0 : const auto & jv = _sys.getVariable(_tid, jvar);
2934 0 : const unsigned int ivn = iv.number();
2935 0 : const unsigned int jvn = jv.number();
2936 0 : const unsigned int icount = iv.count();
2937 0 : unsigned int jcount = jv.count();
2938 0 : if (ivn == jvn && _component_block_diagonal[ivn])
2939 0 : jcount = 1;
2940 :
2941 0 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2942 0 : tag < _jacobian_block_nonlocal_used.size();
2943 : tag++)
2944 : {
2945 0 : jacobianBlockNonlocal(ivn, jvn, LocalDataKey{}, tag)
2946 0 : .resize(idof_indices.size() * icount, jdof_indices.size() * jcount);
2947 :
2948 0 : jacobianBlockNonlocalUsed(tag, ivn, jvn, false);
2949 : }
2950 0 : }
2951 :
2952 : void
2953 8711897 : Assembly::prepareScalar()
2954 : {
2955 8711897 : const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
2956 8978216 : for (const auto & ivar : vars)
2957 : {
2958 266319 : auto idofs = ivar->dofIndices().size();
2959 :
2960 1052029 : for (auto & tag_Re : _sub_Re)
2961 785710 : tag_Re[ivar->number()].resize(idofs);
2962 :
2963 680436 : for (const auto & jvar : vars)
2964 : {
2965 414117 : auto jdofs = jvar->dofIndices().size();
2966 :
2967 1245629 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2968 : {
2969 831512 : jacobianBlock(ivar->number(), jvar->number(), LocalDataKey{}, tag).resize(idofs, jdofs);
2970 831512 : jacobianBlockUsed(tag, ivar->number(), jvar->number(), false);
2971 : }
2972 : }
2973 : }
2974 8711897 : }
2975 :
2976 : void
2977 185626 : Assembly::prepareOffDiagScalar()
2978 : {
2979 185626 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2980 185626 : const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
2981 :
2982 376665 : for (const auto & ivar : scalar_vars)
2983 : {
2984 191039 : auto idofs = ivar->dofIndices().size();
2985 :
2986 474044 : for (const auto & jvar : vars)
2987 : {
2988 283005 : auto jdofs = jvar->dofIndices().size() * jvar->count();
2989 849015 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2990 : {
2991 566010 : jacobianBlock(ivar->number(), jvar->number(), LocalDataKey{}, tag).resize(idofs, jdofs);
2992 566010 : jacobianBlockUsed(tag, ivar->number(), jvar->number(), false);
2993 :
2994 566010 : jacobianBlock(jvar->number(), ivar->number(), LocalDataKey{}, tag).resize(jdofs, idofs);
2995 566010 : jacobianBlockUsed(tag, jvar->number(), ivar->number(), false);
2996 : }
2997 : }
2998 : }
2999 185626 : }
3000 :
3001 : template <typename T>
3002 : void
3003 126180081 : Assembly::copyShapes(MooseVariableField<T> & v)
3004 : {
3005 126180081 : phi(v).shallowCopy(v.phi());
3006 126180081 : gradPhi(v).shallowCopy(v.gradPhi());
3007 126180081 : if (v.computingSecond())
3008 23852 : secondPhi(v).shallowCopy(v.secondPhi());
3009 126180081 : }
3010 :
3011 : void
3012 126180081 : Assembly::copyShapes(unsigned int var)
3013 : {
3014 126180081 : auto & v = _sys.getVariable(_tid, var);
3015 126180081 : if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3016 : {
3017 123799331 : auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3018 123799331 : copyShapes(v);
3019 : }
3020 2380750 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3021 : {
3022 562980 : auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3023 562980 : copyShapes(v);
3024 : }
3025 1817770 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3026 : {
3027 1817770 : auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3028 1817770 : copyShapes(v);
3029 1817770 : if (v.computingCurl())
3030 50508 : curlPhi(v).shallowCopy(v.curlPhi());
3031 1817770 : if (v.computingDiv())
3032 494208 : divPhi(v).shallowCopy(v.divPhi());
3033 : }
3034 : else
3035 0 : mooseError("Unsupported variable field type!");
3036 126180081 : }
3037 :
3038 : template <typename T>
3039 : void
3040 611444 : Assembly::copyFaceShapes(MooseVariableField<T> & v)
3041 : {
3042 611444 : phiFace(v).shallowCopy(v.phiFace());
3043 611444 : gradPhiFace(v).shallowCopy(v.gradPhiFace());
3044 611444 : if (v.computingSecond())
3045 4088 : secondPhiFace(v).shallowCopy(v.secondPhiFace());
3046 611444 : }
3047 :
3048 : void
3049 611444 : Assembly::copyFaceShapes(unsigned int var)
3050 : {
3051 611444 : auto & v = _sys.getVariable(_tid, var);
3052 611444 : if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3053 : {
3054 542102 : auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3055 542102 : copyFaceShapes(v);
3056 : }
3057 69342 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3058 : {
3059 12810 : auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3060 12810 : copyFaceShapes(v);
3061 : }
3062 56532 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3063 : {
3064 56532 : auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3065 56532 : copyFaceShapes(v);
3066 56532 : if (v.computingCurl())
3067 6380 : _vector_curl_phi_face.shallowCopy(v.curlPhi());
3068 56532 : if (v.computingDiv())
3069 26112 : _vector_div_phi_face.shallowCopy(v.divPhi());
3070 : }
3071 : else
3072 0 : mooseError("Unsupported variable field type!");
3073 611444 : }
3074 :
3075 : template <typename T>
3076 : void
3077 187886 : Assembly::copyNeighborShapes(MooseVariableField<T> & v)
3078 : {
3079 187886 : if (v.usesPhiNeighbor())
3080 : {
3081 187886 : phiFaceNeighbor(v).shallowCopy(v.phiFaceNeighbor());
3082 187886 : phiNeighbor(v).shallowCopy(v.phiNeighbor());
3083 : }
3084 187886 : if (v.usesGradPhiNeighbor())
3085 : {
3086 187886 : gradPhiFaceNeighbor(v).shallowCopy(v.gradPhiFaceNeighbor());
3087 187886 : gradPhiNeighbor(v).shallowCopy(v.gradPhiNeighbor());
3088 : }
3089 187886 : if (v.usesSecondPhiNeighbor())
3090 : {
3091 0 : secondPhiFaceNeighbor(v).shallowCopy(v.secondPhiFaceNeighbor());
3092 0 : secondPhiNeighbor(v).shallowCopy(v.secondPhiNeighbor());
3093 : }
3094 187886 : }
3095 :
3096 : void
3097 187886 : Assembly::copyNeighborShapes(unsigned int var)
3098 : {
3099 187886 : auto & v = _sys.getVariable(_tid, var);
3100 187886 : if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3101 : {
3102 184164 : auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3103 184164 : copyNeighborShapes(v);
3104 : }
3105 3722 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3106 : {
3107 3456 : auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3108 3456 : copyNeighborShapes(v);
3109 : }
3110 266 : else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3111 : {
3112 266 : auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3113 266 : copyNeighborShapes(v);
3114 : }
3115 : else
3116 0 : mooseError("Unsupported variable field type!");
3117 187886 : }
3118 :
3119 : DenseMatrix<Number> &
3120 226256578 : Assembly::jacobianBlockNeighbor(
3121 : Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
3122 : {
3123 226256578 : if (type == Moose::ElementElement)
3124 0 : jacobianBlockUsed(tag, ivar, jvar, true);
3125 : else
3126 226256578 : jacobianBlockNeighborUsed(tag, ivar, jvar, true);
3127 :
3128 226256578 : if (_block_diagonal_matrix)
3129 : {
3130 120082487 : switch (type)
3131 : {
3132 0 : default:
3133 : case Moose::ElementElement:
3134 0 : return _sub_Kee[tag][ivar][0];
3135 40028856 : case Moose::ElementNeighbor:
3136 40028856 : return _sub_Ken[tag][ivar][0];
3137 40022515 : case Moose::NeighborElement:
3138 40022515 : return _sub_Kne[tag][ivar][0];
3139 40031116 : case Moose::NeighborNeighbor:
3140 40031116 : return _sub_Knn[tag][ivar][0];
3141 : }
3142 : }
3143 : else
3144 : {
3145 106174091 : switch (type)
3146 : {
3147 0 : default:
3148 : case Moose::ElementElement:
3149 0 : return _sub_Kee[tag][ivar][jvar];
3150 35391593 : case Moose::ElementNeighbor:
3151 35391593 : return _sub_Ken[tag][ivar][jvar];
3152 35390905 : case Moose::NeighborElement:
3153 35390905 : return _sub_Kne[tag][ivar][jvar];
3154 35391593 : case Moose::NeighborNeighbor:
3155 35391593 : return _sub_Knn[tag][ivar][jvar];
3156 : }
3157 : }
3158 : }
3159 :
3160 : DenseMatrix<Number> &
3161 22761044 : Assembly::jacobianBlockMortar(Moose::ConstraintJacobianType type,
3162 : unsigned int ivar,
3163 : unsigned int jvar,
3164 : LocalDataKey,
3165 : TagID tag)
3166 : {
3167 22761044 : jacobianBlockLowerUsed(tag, ivar, jvar, true);
3168 22761044 : if (_block_diagonal_matrix)
3169 : {
3170 1897800 : switch (type)
3171 : {
3172 323576 : default:
3173 : case Moose::LowerLower:
3174 323576 : return _sub_Kll[tag][ivar][0];
3175 323192 : case Moose::LowerSecondary:
3176 323192 : return _sub_Kle[tag][ivar][0];
3177 323000 : case Moose::LowerPrimary:
3178 323000 : return _sub_Kln[tag][ivar][0];
3179 351376 : case Moose::SecondaryLower:
3180 351376 : return _sub_Kel[tag][ivar][0];
3181 56368 : case Moose::SecondarySecondary:
3182 56368 : return _sub_Kee[tag][ivar][0];
3183 56368 : case Moose::SecondaryPrimary:
3184 56368 : return _sub_Ken[tag][ivar][0];
3185 351184 : case Moose::PrimaryLower:
3186 351184 : return _sub_Knl[tag][ivar][0];
3187 56368 : case Moose::PrimarySecondary:
3188 56368 : return _sub_Kne[tag][ivar][0];
3189 56368 : case Moose::PrimaryPrimary:
3190 56368 : return _sub_Knn[tag][ivar][0];
3191 : }
3192 : }
3193 : else
3194 : {
3195 20863244 : switch (type)
3196 : {
3197 4154188 : default:
3198 : case Moose::LowerLower:
3199 4154188 : return _sub_Kll[tag][ivar][jvar];
3200 4153308 : case Moose::LowerSecondary:
3201 4153308 : return _sub_Kle[tag][ivar][jvar];
3202 4145740 : case Moose::LowerPrimary:
3203 4145740 : return _sub_Kln[tag][ivar][jvar];
3204 4155468 : case Moose::SecondaryLower:
3205 4155468 : return _sub_Kel[tag][ivar][jvar];
3206 26660 : case Moose::SecondarySecondary:
3207 26660 : return _sub_Kee[tag][ivar][jvar];
3208 26660 : case Moose::SecondaryPrimary:
3209 26660 : return _sub_Ken[tag][ivar][jvar];
3210 4147900 : case Moose::PrimaryLower:
3211 4147900 : return _sub_Knl[tag][ivar][jvar];
3212 26660 : case Moose::PrimarySecondary:
3213 26660 : return _sub_Kne[tag][ivar][jvar];
3214 26660 : case Moose::PrimaryPrimary:
3215 26660 : return _sub_Knn[tag][ivar][jvar];
3216 : }
3217 : }
3218 : }
3219 :
3220 : void
3221 970203856 : Assembly::processLocalResidual(DenseVector<Number> & res_block,
3222 : std::vector<dof_id_type> & dof_indices,
3223 : const std::vector<Real> & scaling_factor)
3224 : {
3225 : mooseAssert(res_block.size() == dof_indices.size(),
3226 : "The size of residual and degree of freedom container must be the same");
3227 :
3228 : // For an array variable, ndof is the number of dofs of the zero-th component and
3229 : // ntdof is the number of dofs of all components.
3230 : // For standard or vector variables, ndof will be the same as ntdof.
3231 970203856 : const auto ntdof = res_block.size();
3232 970203856 : const auto count = scaling_factor.size();
3233 970203856 : const auto ndof = ntdof / count;
3234 970203856 : if (count > 1)
3235 : {
3236 4392932 : unsigned int p = 0;
3237 13718346 : for (MooseIndex(count) j = 0; j < count; ++j)
3238 52453508 : for (MooseIndex(ndof) i = 0; i < ndof; ++i)
3239 43128094 : res_block(p++) *= scaling_factor[j];
3240 : }
3241 : else
3242 : {
3243 965810924 : if (scaling_factor[0] != 1.0)
3244 4268837 : res_block *= scaling_factor[0];
3245 : }
3246 :
3247 970203856 : _dof_map.constrain_element_vector(res_block, dof_indices, false);
3248 970203856 : }
3249 :
3250 : void
3251 13691991 : Assembly::addResidualBlock(NumericVector<Number> & residual,
3252 : DenseVector<Number> & res_block,
3253 : const std::vector<dof_id_type> & dof_indices,
3254 : const std::vector<Real> & scaling_factor)
3255 : {
3256 13691991 : if (dof_indices.size() > 0 && res_block.size())
3257 : {
3258 7007030 : _temp_dof_indices = dof_indices;
3259 7007030 : _tmp_Re = res_block;
3260 7007030 : processLocalResidual(_tmp_Re, _temp_dof_indices, scaling_factor);
3261 7007030 : residual.add_vector(_tmp_Re, _temp_dof_indices);
3262 : }
3263 13691991 : }
3264 :
3265 : void
3266 1032307164 : Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
3267 : std::vector<dof_id_type> & cached_residual_rows,
3268 : DenseVector<Number> & res_block,
3269 : const std::vector<dof_id_type> & dof_indices,
3270 : const std::vector<Real> & scaling_factor)
3271 : {
3272 1032307164 : if (dof_indices.size() > 0 && res_block.size())
3273 : {
3274 963196826 : _temp_dof_indices = dof_indices;
3275 963196826 : _tmp_Re = res_block;
3276 963196826 : processLocalResidual(_tmp_Re, _temp_dof_indices, scaling_factor);
3277 :
3278 4915160637 : for (MooseIndex(_tmp_Re) i = 0; i < _tmp_Re.size(); i++)
3279 : {
3280 3951963811 : cached_residual_values.push_back(_tmp_Re(i));
3281 3951963811 : cached_residual_rows.push_back(_temp_dof_indices[i]);
3282 : }
3283 : }
3284 :
3285 1032307164 : res_block.zero();
3286 1032307164 : }
3287 :
3288 : void
3289 0 : Assembly::setResidualBlock(NumericVector<Number> & residual,
3290 : DenseVector<Number> & res_block,
3291 : const std::vector<dof_id_type> & dof_indices,
3292 : const std::vector<Real> & scaling_factor)
3293 : {
3294 0 : if (dof_indices.size() > 0)
3295 : {
3296 0 : std::vector<dof_id_type> di(dof_indices);
3297 0 : _tmp_Re = res_block;
3298 0 : processLocalResidual(_tmp_Re, di, scaling_factor);
3299 0 : residual.insert(_tmp_Re, di);
3300 0 : }
3301 0 : }
3302 :
3303 : void
3304 786387 : Assembly::addResidual(const VectorTag & vector_tag)
3305 : {
3306 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3307 : "Non-residual tag in Assembly::addResidual");
3308 :
3309 786387 : auto & tag_Re = _sub_Re[vector_tag._type_id];
3310 786387 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3311 786387 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3312 2264954 : for (const auto & var : vars)
3313 1478567 : addResidualBlock(residual, tag_Re[var->number()], var->dofIndices(), var->arrayScalingFactor());
3314 786387 : }
3315 :
3316 : void
3317 271389 : Assembly::addResidual(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3318 : {
3319 1057776 : for (const auto & vector_tag : vector_tags)
3320 786387 : if (_sys.hasVector(vector_tag._id))
3321 786387 : addResidual(vector_tag);
3322 271389 : }
3323 :
3324 : void
3325 5011955 : Assembly::addResidualNeighbor(const VectorTag & vector_tag)
3326 : {
3327 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3328 : "Non-residual tag in Assembly::addResidualNeighbor");
3329 :
3330 5011955 : auto & tag_Rn = _sub_Rn[vector_tag._type_id];
3331 5011955 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3332 5011955 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3333 11036874 : for (const auto & var : vars)
3334 6024919 : addResidualBlock(
3335 6024919 : residual, tag_Rn[var->number()], var->dofIndicesNeighbor(), var->arrayScalingFactor());
3336 5011955 : }
3337 :
3338 : void
3339 2046393 : Assembly::addResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3340 : {
3341 7058348 : for (const auto & vector_tag : vector_tags)
3342 5011955 : if (_sys.hasVector(vector_tag._id))
3343 5011955 : addResidualNeighbor(vector_tag);
3344 2046393 : }
3345 :
3346 : void
3347 4985617 : Assembly::addResidualLower(const VectorTag & vector_tag)
3348 : {
3349 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3350 : "Non-residual tag in Assembly::addResidualLower");
3351 :
3352 4985617 : auto & tag_Rl = _sub_Rl[vector_tag._type_id];
3353 4985617 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3354 4985617 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3355 10977796 : for (const auto & var : vars)
3356 5992179 : addResidualBlock(
3357 5992179 : residual, tag_Rl[var->number()], var->dofIndicesLower(), var->arrayScalingFactor());
3358 4985617 : }
3359 :
3360 : void
3361 2028936 : Assembly::addResidualLower(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3362 : {
3363 7014553 : for (const auto & vector_tag : vector_tags)
3364 4985617 : if (_sys.hasVector(vector_tag._id))
3365 4985617 : addResidualLower(vector_tag);
3366 2028936 : }
3367 :
3368 : // private method, so no key required
3369 : void
3370 145140 : Assembly::addResidualScalar(const VectorTag & vector_tag)
3371 : {
3372 : mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3373 : "Non-residual tag in Assembly::addResidualScalar");
3374 :
3375 : // add the scalar variables residuals
3376 145140 : auto & tag_Re = _sub_Re[vector_tag._type_id];
3377 145140 : NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3378 145140 : const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
3379 341466 : for (const auto & var : vars)
3380 196326 : addResidualBlock(residual, tag_Re[var->number()], var->dofIndices(), var->arrayScalingFactor());
3381 145140 : }
3382 :
3383 : void
3384 48426 : Assembly::addResidualScalar(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3385 : {
3386 193566 : for (const auto & vector_tag : vector_tags)
3387 145140 : if (_sys.hasVector(vector_tag._id))
3388 145140 : addResidualScalar(vector_tag);
3389 48426 : }
3390 :
3391 : void
3392 321282928 : Assembly::cacheResidual(GlobalDataKey, const std::vector<VectorTag> & tags)
3393 : {
3394 321282928 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3395 673291306 : for (const auto & var : vars)
3396 1310646489 : for (const auto & vector_tag : tags)
3397 958638111 : if (_sys.hasVector(vector_tag._id))
3398 958638111 : cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3399 958638111 : _cached_residual_rows[vector_tag._type_id],
3400 958638111 : _sub_Re[vector_tag._type_id][var->number()],
3401 958638111 : var->dofIndices(),
3402 958638111 : var->arrayScalingFactor());
3403 321282928 : }
3404 :
3405 : // private method, so no key required
3406 : void
3407 70092390 : Assembly::cacheResidual(dof_id_type dof, Real value, TagID tag_id)
3408 : {
3409 70092390 : const VectorTag & tag = _subproblem.getVectorTag(tag_id);
3410 :
3411 70092390 : _cached_residual_values[tag._type_id].push_back(value);
3412 70092390 : _cached_residual_rows[tag._type_id].push_back(dof);
3413 70092390 : }
3414 :
3415 : // private method, so no key required
3416 : void
3417 69656283 : Assembly::cacheResidual(dof_id_type dof, Real value, const std::set<TagID> & tags)
3418 : {
3419 139748673 : for (auto & tag : tags)
3420 70092390 : cacheResidual(dof, value, tag);
3421 69656283 : }
3422 :
3423 : void
3424 0 : Assembly::cacheResidualNodes(const DenseVector<Number> & res,
3425 : const std::vector<dof_id_type> & dof_index,
3426 : LocalDataKey,
3427 : TagID tag)
3428 : {
3429 : // Add the residual value and dof_index to cached_residual_values and cached_residual_rows
3430 : // respectively.
3431 : // This is used by NodalConstraint.C to cache the residual calculated for primary and secondary
3432 : // node.
3433 0 : const VectorTag & vector_tag = _subproblem.getVectorTag(tag);
3434 0 : for (MooseIndex(dof_index) i = 0; i < dof_index.size(); ++i)
3435 : {
3436 0 : _cached_residual_values[vector_tag._type_id].push_back(res(i));
3437 0 : _cached_residual_rows[vector_tag._type_id].push_back(dof_index[i]);
3438 : }
3439 0 : }
3440 :
3441 : void
3442 25452135 : Assembly::cacheResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & tags)
3443 : {
3444 25452135 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3445 55149414 : for (const auto & var : vars)
3446 101998852 : for (const auto & vector_tag : tags)
3447 72301573 : if (_sys.hasVector(vector_tag._id))
3448 72301573 : cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3449 72301573 : _cached_residual_rows[vector_tag._type_id],
3450 72301573 : _sub_Rn[vector_tag._type_id][var->number()],
3451 72301573 : var->dofIndicesNeighbor(),
3452 72301573 : var->arrayScalingFactor());
3453 25452135 : }
3454 :
3455 : void
3456 373958 : Assembly::cacheResidualLower(GlobalDataKey, const std::vector<VectorTag> & tags)
3457 : {
3458 373958 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3459 1042142 : for (const auto & var : vars)
3460 2035664 : for (const auto & vector_tag : tags)
3461 1367480 : if (_sys.hasVector(vector_tag._id))
3462 1367480 : cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3463 1367480 : _cached_residual_rows[vector_tag._type_id],
3464 1367480 : _sub_Rl[vector_tag._type_id][var->number()],
3465 1367480 : var->dofIndicesLower(),
3466 1367480 : var->arrayScalingFactor());
3467 373958 : }
3468 :
3469 : void
3470 42340558 : Assembly::addCachedResiduals(GlobalDataKey, const std::vector<VectorTag> & tags)
3471 : {
3472 152825785 : for (const auto & vector_tag : tags)
3473 : {
3474 110485227 : if (!_sys.hasVector(vector_tag._id))
3475 : {
3476 0 : _cached_residual_values[vector_tag._type_id].clear();
3477 0 : _cached_residual_rows[vector_tag._type_id].clear();
3478 0 : continue;
3479 : }
3480 110485227 : addCachedResidualDirectly(_sys.getVector(vector_tag._id), GlobalDataKey{}, vector_tag);
3481 : }
3482 42340558 : }
3483 :
3484 : void
3485 11367 : Assembly::clearCachedResiduals(GlobalDataKey)
3486 : {
3487 44257 : for (const auto & vector_tag : _residual_vector_tags)
3488 32890 : clearCachedResiduals(vector_tag);
3489 11367 : }
3490 :
3491 : // private method, so no key required
3492 : void
3493 101832346 : Assembly::clearCachedResiduals(const VectorTag & vector_tag)
3494 : {
3495 101832346 : auto & values = _cached_residual_values[vector_tag._type_id];
3496 101832346 : auto & rows = _cached_residual_rows[vector_tag._type_id];
3497 :
3498 : mooseAssert(values.size() == rows.size(),
3499 : "Number of cached residuals and number of rows must match!");
3500 :
3501 : // Keep track of the largest size so we can use it to reserve and avoid
3502 : // as much dynamic allocation as possible
3503 101832346 : if (_max_cached_residuals < values.size())
3504 41760 : _max_cached_residuals = values.size();
3505 :
3506 : // Clear both vectors (keeps the capacity the same)
3507 101832346 : values.clear();
3508 101832346 : rows.clear();
3509 : // And then reserve: use 2 as a fudge factor to *really* avoid dynamic allocation!
3510 101832346 : values.reserve(_max_cached_residuals * 2);
3511 101832346 : rows.reserve(_max_cached_residuals * 2);
3512 101832346 : }
3513 :
3514 : void
3515 110506678 : Assembly::addCachedResidualDirectly(NumericVector<Number> & residual,
3516 : GlobalDataKey,
3517 : const VectorTag & vector_tag)
3518 : {
3519 110506678 : const auto & values = _cached_residual_values[vector_tag._type_id];
3520 110506678 : const auto & rows = _cached_residual_rows[vector_tag._type_id];
3521 :
3522 : mooseAssert(values.size() == rows.size(),
3523 : "Number of cached residuals and number of rows must match!");
3524 :
3525 110506678 : if (!values.empty())
3526 : {
3527 101799456 : residual.add_vector(values, rows);
3528 101799456 : clearCachedResiduals(vector_tag);
3529 : }
3530 110506678 : }
3531 :
3532 : void
3533 0 : Assembly::setResidual(NumericVector<Number> & residual, GlobalDataKey, const VectorTag & vector_tag)
3534 : {
3535 0 : auto & tag_Re = _sub_Re[vector_tag._type_id];
3536 0 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3537 0 : for (const auto & var : vars)
3538 0 : setResidualBlock(residual, tag_Re[var->number()], var->dofIndices(), var->arrayScalingFactor());
3539 0 : }
3540 :
3541 : void
3542 0 : Assembly::setResidualNeighbor(NumericVector<Number> & residual,
3543 : GlobalDataKey,
3544 : const VectorTag & vector_tag)
3545 : {
3546 0 : auto & tag_Rn = _sub_Rn[vector_tag._type_id];
3547 0 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3548 0 : for (const auto & var : vars)
3549 0 : setResidualBlock(
3550 0 : residual, tag_Rn[var->number()], var->dofIndicesNeighbor(), var->arrayScalingFactor());
3551 0 : }
3552 :
3553 : // private method, so no key required
3554 : void
3555 432160 : Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
3556 : DenseMatrix<Number> & jac_block,
3557 : const MooseVariableBase & ivar,
3558 : const MooseVariableBase & jvar,
3559 : const std::vector<dof_id_type> & idof_indices,
3560 : const std::vector<dof_id_type> & jdof_indices)
3561 : {
3562 432160 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3563 86478 : return;
3564 345682 : if (jac_block.n() == 0 || jac_block.m() == 0)
3565 0 : return;
3566 :
3567 345682 : const auto & scaling_factors = ivar.arrayScalingFactor();
3568 345682 : const unsigned int iv = ivar.number();
3569 345682 : const unsigned int jv = jvar.number();
3570 :
3571 705092 : for (unsigned int i = 0; i < ivar.count(); ++i)
3572 : {
3573 873841 : for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3574 : {
3575 514431 : if (jt < jv || jt >= jv + jvar.count())
3576 129485 : continue;
3577 384946 : unsigned int j = jt - jv;
3578 :
3579 384946 : auto di = ivar.componentDofIndices(idof_indices, i);
3580 384946 : auto dj = jvar.componentDofIndices(jdof_indices, j);
3581 384946 : auto indof = di.size();
3582 384946 : auto jndof = dj.size();
3583 :
3584 384946 : unsigned int jj = j;
3585 384946 : if (iv == jv && _component_block_diagonal[iv])
3586 : // here i must be equal to j
3587 309123 : jj = 0;
3588 :
3589 384946 : auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3590 384946 : if (scaling_factors[i] != 1.0)
3591 5216 : sub *= scaling_factors[i];
3592 :
3593 : // If we're computing the jacobian for automatically scaling variables we do not want
3594 : // to constrain the element matrix because it introduces 1s on the diagonal for the
3595 : // constrained dofs
3596 384946 : if (!_sys.computingScalingJacobian())
3597 384890 : _dof_map.constrain_element_matrix(sub, di, dj, false);
3598 :
3599 384946 : jacobian.add_matrix(sub, di, dj);
3600 384946 : }
3601 : }
3602 : }
3603 :
3604 : // private method, so no key required
3605 : void
3606 58335441 : Assembly::cacheJacobianBlock(const DenseMatrix<Number> & jac_block,
3607 : const MooseVariableBase & ivar,
3608 : const MooseVariableBase & jvar,
3609 : const std::vector<dof_id_type> & idof_indices,
3610 : const std::vector<dof_id_type> & jdof_indices,
3611 : TagID tag)
3612 : {
3613 58335441 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3614 390862 : return;
3615 57944579 : if (jac_block.n() == 0 || jac_block.m() == 0)
3616 0 : return;
3617 57944579 : if (!_sys.hasMatrix(tag))
3618 0 : return;
3619 :
3620 57944579 : auto & scaling_factors = ivar.arrayScalingFactor();
3621 57944579 : const unsigned int iv = ivar.number();
3622 57944579 : const unsigned int jv = jvar.number();
3623 :
3624 116135744 : for (unsigned int i = 0; i < ivar.count(); ++i)
3625 : {
3626 144310786 : for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3627 : {
3628 86119621 : if (jt < jv || jt >= jv + jvar.count())
3629 27828776 : continue;
3630 58290845 : unsigned int j = jt - jv;
3631 :
3632 58290845 : auto di = ivar.componentDofIndices(idof_indices, i);
3633 58290845 : auto dj = jvar.componentDofIndices(jdof_indices, j);
3634 58290845 : auto indof = di.size();
3635 58290845 : auto jndof = dj.size();
3636 :
3637 58290845 : unsigned int jj = j;
3638 58290845 : if (iv == jv && _component_block_diagonal[iv])
3639 : // here i must be equal to j
3640 48660161 : jj = 0;
3641 :
3642 58290845 : auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3643 58290845 : if (scaling_factors[i] != 1.0)
3644 1558874 : sub *= scaling_factors[i];
3645 :
3646 : // If we're computing the jacobian for automatically scaling variables we do not want
3647 : // to constrain the element matrix because it introduces 1s on the diagonal for the
3648 : // constrained dofs
3649 58290845 : if (!_sys.computingScalingJacobian())
3650 58094679 : _dof_map.constrain_element_matrix(sub, di, dj, false);
3651 :
3652 318449156 : for (MooseIndex(di) i = 0; i < di.size(); i++)
3653 1601078944 : for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3654 : {
3655 1340920633 : _cached_jacobian_values[tag].push_back(sub(i, j));
3656 1340920633 : _cached_jacobian_rows[tag].push_back(di[i]);
3657 1340920633 : _cached_jacobian_cols[tag].push_back(dj[j]);
3658 : }
3659 58290845 : }
3660 : }
3661 : }
3662 :
3663 : // private method, so no key required
3664 : void
3665 634 : Assembly::cacheJacobianBlockNonzero(const DenseMatrix<Number> & jac_block,
3666 : const MooseVariableBase & ivar,
3667 : const MooseVariableBase & jvar,
3668 : const std::vector<dof_id_type> & idof_indices,
3669 : const std::vector<dof_id_type> & jdof_indices,
3670 : TagID tag)
3671 : {
3672 634 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3673 0 : return;
3674 634 : if (jac_block.n() == 0 || jac_block.m() == 0)
3675 0 : return;
3676 634 : if (!_sys.hasMatrix(tag))
3677 0 : return;
3678 :
3679 634 : auto & scaling_factor = ivar.arrayScalingFactor();
3680 :
3681 1268 : for (unsigned int i = 0; i < ivar.count(); ++i)
3682 : {
3683 634 : unsigned int iv = ivar.number();
3684 1862 : for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3685 : {
3686 1228 : unsigned int jv = jvar.number();
3687 1228 : if (jt < jv || jt >= jv + jvar.count())
3688 594 : continue;
3689 634 : unsigned int j = jt - jv;
3690 :
3691 634 : auto di = ivar.componentDofIndices(idof_indices, i);
3692 634 : auto dj = jvar.componentDofIndices(jdof_indices, j);
3693 634 : auto indof = di.size();
3694 634 : auto jndof = dj.size();
3695 :
3696 634 : unsigned int jj = j;
3697 634 : if (iv == jv && _component_block_diagonal[iv])
3698 : // here i must be equal to j
3699 152 : jj = 0;
3700 :
3701 634 : auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3702 634 : if (scaling_factor[i] != 1.0)
3703 284 : sub *= scaling_factor[i];
3704 :
3705 634 : _dof_map.constrain_element_matrix(sub, di, dj, false);
3706 :
3707 3122 : for (MooseIndex(di) i = 0; i < di.size(); i++)
3708 195472 : for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3709 192984 : if (sub(i, j) != 0.0) // no storage allocated for unimplemented jacobian terms,
3710 : // maintaining maximum sparsity possible
3711 : {
3712 25528 : _cached_jacobian_values[tag].push_back(sub(i, j));
3713 25528 : _cached_jacobian_rows[tag].push_back(di[i]);
3714 25528 : _cached_jacobian_cols[tag].push_back(dj[j]);
3715 : }
3716 634 : }
3717 : }
3718 : }
3719 :
3720 : void
3721 2952760 : Assembly::cacheJacobianBlock(const DenseMatrix<Number> & jac_block,
3722 : const std::vector<dof_id_type> & idof_indices,
3723 : const std::vector<dof_id_type> & jdof_indices,
3724 : Real scaling_factor,
3725 : LocalDataKey,
3726 : const std::set<TagID> & tags)
3727 : {
3728 : const auto has_matrix =
3729 5905520 : std::any_of(tags.begin(), tags.end(), [this](const auto tag) { return _sys.hasMatrix(tag); });
3730 :
3731 : // Work on a reusable Assembly-owned copy so callers retain their local matrix. This also lets us
3732 : // apply constraints and scaling once before caching the same block to every requested matrix tag.
3733 2952760 : if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m() &&
3734 : has_matrix)
3735 : {
3736 2932066 : _row_indices.assign(idof_indices.begin(), idof_indices.end());
3737 2932066 : _column_indices.assign(jdof_indices.begin(), jdof_indices.end());
3738 2932066 : _element_matrix = jac_block;
3739 :
3740 : // If we're computing the jacobian for automatically scaling variables we do not want to
3741 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
3742 : // dofs
3743 2932066 : if (!_sys.computingScalingJacobian())
3744 2932066 : _dof_map.constrain_element_matrix(_element_matrix, _row_indices, _column_indices, false);
3745 :
3746 2932066 : if (scaling_factor != 1.0)
3747 6664 : _element_matrix *= scaling_factor;
3748 :
3749 17329961 : for (const auto i : index_range(_row_indices))
3750 95087812 : for (const auto j : index_range(_column_indices))
3751 80689917 : cacheJacobian(
3752 161379834 : _row_indices[i], _column_indices[j], _element_matrix(i, j), LocalDataKey{}, tags);
3753 : }
3754 2952760 : }
3755 :
3756 : Real
3757 0 : Assembly::elementVolume(const Elem * elem) const
3758 : {
3759 0 : FEType fe_type(elem->default_order(), LAGRANGE);
3760 0 : std::unique_ptr<FEBase> fe(FEBase::build(elem->dim(), fe_type));
3761 :
3762 : // references to the quadrature points and weights
3763 0 : const std::vector<Real> & JxW = fe->get_JxW();
3764 0 : const std::vector<Point> & q_points = fe->get_xyz();
3765 :
3766 : // The default quadrature rule should integrate the mass matrix,
3767 : // thus it should be plenty to compute the volume
3768 0 : QGauss qrule(elem->dim(), fe_type.default_quadrature_order());
3769 0 : fe->attach_quadrature_rule(&qrule);
3770 0 : fe->reinit(elem);
3771 :
3772 : // perform a sanity check to ensure that size of quad rule and size of q_points is
3773 : // identical
3774 : mooseAssert(qrule.n_points() == q_points.size(),
3775 : "The number of points in the quadrature rule doesn't match the number of passed-in "
3776 : "points in Assembly::setCoordinateTransformation");
3777 :
3778 : // compute the coordinate transformation
3779 0 : Real vol = 0;
3780 0 : for (unsigned int qp = 0; qp < qrule.n_points(); ++qp)
3781 : {
3782 : Real coord;
3783 0 : coordTransformFactor(_subproblem, elem->subdomain_id(), q_points[qp], coord);
3784 0 : vol += JxW[qp] * coord;
3785 : }
3786 0 : return vol;
3787 0 : }
3788 :
3789 : void
3790 21504 : Assembly::saveLocalADArray(std::vector<ADReal> & re,
3791 : unsigned int i,
3792 : unsigned int ntest,
3793 : const ADRealEigenVector & v) const
3794 : {
3795 64512 : for (unsigned int j = 0; j < v.size(); ++j, i += ntest)
3796 43008 : re[i] += v(j);
3797 21504 : }
3798 :
3799 : void
3800 3450273 : Assembly::addCachedJacobian(GlobalDataKey)
3801 : {
3802 : #ifndef NDEBUG
3803 : if (!_subproblem.checkNonlocalCouplingRequirement())
3804 : {
3805 : mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
3806 : "Error: Cached data sizes MUST be the same!");
3807 : for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3808 : mooseAssert(_cached_jacobian_rows[i].size() == _cached_jacobian_cols[i].size(),
3809 : "Error: Cached data sizes MUST be the same for a given tag!");
3810 : }
3811 : #endif
3812 :
3813 10422651 : for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3814 6972381 : if (_sys.hasMatrix(i))
3815 1560136186 : for (MooseIndex(_cached_jacobian_rows[i]) j = 0; j < _cached_jacobian_rows[i].size(); j++)
3816 3113340136 : _sys.getMatrix(i).add(_cached_jacobian_rows[i][j],
3817 1556670068 : _cached_jacobian_cols[i][j],
3818 1556670068 : _cached_jacobian_values[i][j]);
3819 :
3820 10422648 : for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3821 : {
3822 6972378 : if (!_sys.hasMatrix(i))
3823 3506260 : continue;
3824 :
3825 3466118 : if (_max_cached_jacobians < _cached_jacobian_values[i].size())
3826 42697 : _max_cached_jacobians = _cached_jacobian_values[i].size();
3827 :
3828 : // Try to be more efficient from now on
3829 : // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
3830 3466118 : _cached_jacobian_values[i].clear();
3831 3466118 : _cached_jacobian_values[i].reserve(_max_cached_jacobians * 2);
3832 :
3833 3466118 : _cached_jacobian_rows[i].clear();
3834 3466118 : _cached_jacobian_rows[i].reserve(_max_cached_jacobians * 2);
3835 :
3836 3466118 : _cached_jacobian_cols[i].clear();
3837 3466118 : _cached_jacobian_cols[i].reserve(_max_cached_jacobians * 2);
3838 : }
3839 3450270 : }
3840 :
3841 : inline void
3842 101576 : Assembly::addJacobianCoupledVarPair(const MooseVariableBase & ivar, const MooseVariableBase & jvar)
3843 : {
3844 101576 : auto i = ivar.number();
3845 101576 : auto j = jvar.number();
3846 305704 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
3847 204128 : if (jacobianBlockUsed(tag, i, j) && _sys.hasMatrix(tag))
3848 53074 : addJacobianBlock(_sys.getMatrix(tag),
3849 106148 : jacobianBlock(i, j, LocalDataKey{}, tag),
3850 : ivar,
3851 : jvar,
3852 53074 : ivar.dofIndices(),
3853 53074 : jvar.dofIndices());
3854 101576 : }
3855 :
3856 : void
3857 37615 : Assembly::addJacobian(GlobalDataKey)
3858 : {
3859 110943 : for (const auto & it : _cm_ff_entry)
3860 73328 : addJacobianCoupledVarPair(*it.first, *it.second);
3861 :
3862 37615 : for (const auto & it : _cm_sf_entry)
3863 0 : addJacobianCoupledVarPair(*it.first, *it.second);
3864 :
3865 37615 : for (const auto & it : _cm_fs_entry)
3866 0 : addJacobianCoupledVarPair(*it.first, *it.second);
3867 37615 : }
3868 :
3869 : void
3870 0 : Assembly::addJacobianNonlocal(GlobalDataKey)
3871 : {
3872 0 : for (const auto & it : _cm_nonlocal_entry)
3873 : {
3874 0 : auto ivar = it.first;
3875 0 : auto jvar = it.second;
3876 0 : auto i = ivar->number();
3877 0 : auto j = jvar->number();
3878 0 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
3879 0 : tag < _jacobian_block_nonlocal_used.size();
3880 : tag++)
3881 0 : if (jacobianBlockNonlocalUsed(tag, i, j) && _sys.hasMatrix(tag))
3882 0 : addJacobianBlock(_sys.getMatrix(tag),
3883 0 : jacobianBlockNonlocal(i, j, LocalDataKey{}, tag),
3884 : *ivar,
3885 : *jvar,
3886 0 : ivar->dofIndices(),
3887 : jvar->allDofIndices());
3888 : }
3889 0 : }
3890 :
3891 : void
3892 7896 : Assembly::addJacobianNeighbor(GlobalDataKey)
3893 : {
3894 37181 : for (const auto & it : _cm_ff_entry)
3895 : {
3896 29285 : auto ivar = it.first;
3897 29285 : auto jvar = it.second;
3898 29285 : auto i = ivar->number();
3899 29285 : auto j = jvar->number();
3900 88239 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3901 88239 : tag < _jacobian_block_neighbor_used.size();
3902 : tag++)
3903 58954 : if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
3904 : {
3905 21952 : addJacobianBlock(_sys.getMatrix(tag),
3906 21952 : jacobianBlockNeighbor(Moose::ElementNeighbor, i, j, LocalDataKey{}, tag),
3907 : *ivar,
3908 : *jvar,
3909 21952 : ivar->dofIndices(),
3910 21952 : jvar->dofIndicesNeighbor());
3911 :
3912 21952 : addJacobianBlock(_sys.getMatrix(tag),
3913 21952 : jacobianBlockNeighbor(Moose::NeighborElement, i, j, LocalDataKey{}, tag),
3914 : *ivar,
3915 : *jvar,
3916 21952 : ivar->dofIndicesNeighbor(),
3917 21952 : jvar->dofIndices());
3918 :
3919 21952 : addJacobianBlock(_sys.getMatrix(tag),
3920 43904 : jacobianBlockNeighbor(Moose::NeighborNeighbor, i, j, LocalDataKey{}, tag),
3921 : *ivar,
3922 : *jvar,
3923 21952 : ivar->dofIndicesNeighbor(),
3924 21952 : jvar->dofIndicesNeighbor());
3925 : }
3926 : }
3927 7896 : }
3928 :
3929 : void
3930 112176 : Assembly::addJacobianNeighborLowerD(GlobalDataKey)
3931 : {
3932 287463 : for (const auto & it : _cm_ff_entry)
3933 : {
3934 175287 : auto ivar = it.first;
3935 175287 : auto jvar = it.second;
3936 175287 : auto i = ivar->number();
3937 175287 : auto j = jvar->number();
3938 530949 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
3939 : tag++)
3940 355662 : if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
3941 : {
3942 7836 : addJacobianBlock(_sys.getMatrix(tag),
3943 7836 : jacobianBlockMortar(Moose::LowerLower, i, j, LocalDataKey{}, tag),
3944 : *ivar,
3945 : *jvar,
3946 7836 : ivar->dofIndicesLower(),
3947 7836 : jvar->dofIndicesLower());
3948 :
3949 7836 : addJacobianBlock(_sys.getMatrix(tag),
3950 7836 : jacobianBlockMortar(Moose::LowerSecondary, i, j, LocalDataKey{}, tag),
3951 : *ivar,
3952 : *jvar,
3953 7836 : ivar->dofIndicesLower(),
3954 7836 : jvar->dofIndicesNeighbor());
3955 :
3956 7836 : addJacobianBlock(_sys.getMatrix(tag),
3957 7836 : jacobianBlockMortar(Moose::LowerPrimary, i, j, LocalDataKey{}, tag),
3958 : *ivar,
3959 : *jvar,
3960 7836 : ivar->dofIndicesLower(),
3961 7836 : jvar->dofIndices());
3962 :
3963 7836 : addJacobianBlock(_sys.getMatrix(tag),
3964 7836 : jacobianBlockMortar(Moose::SecondaryLower, i, j, LocalDataKey{}, tag),
3965 : *ivar,
3966 : *jvar,
3967 7836 : ivar->dofIndicesNeighbor(),
3968 7836 : jvar->dofIndicesLower());
3969 :
3970 7836 : addJacobianBlock(_sys.getMatrix(tag),
3971 15672 : jacobianBlockMortar(Moose::PrimaryLower, i, j, LocalDataKey{}, tag),
3972 : *ivar,
3973 : *jvar,
3974 7836 : ivar->dofIndices(),
3975 7836 : jvar->dofIndicesLower());
3976 : }
3977 :
3978 530949 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3979 530949 : tag < _jacobian_block_neighbor_used.size();
3980 : tag++)
3981 355662 : if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
3982 : {
3983 84646 : addJacobianBlock(_sys.getMatrix(tag),
3984 84646 : jacobianBlockNeighbor(Moose::ElementNeighbor, i, j, LocalDataKey{}, tag),
3985 : *ivar,
3986 : *jvar,
3987 84646 : ivar->dofIndices(),
3988 84646 : jvar->dofIndicesNeighbor());
3989 :
3990 84646 : addJacobianBlock(_sys.getMatrix(tag),
3991 84646 : jacobianBlockNeighbor(Moose::NeighborElement, i, j, LocalDataKey{}, tag),
3992 : *ivar,
3993 : *jvar,
3994 84646 : ivar->dofIndicesNeighbor(),
3995 84646 : jvar->dofIndices());
3996 :
3997 84646 : addJacobianBlock(_sys.getMatrix(tag),
3998 169292 : jacobianBlockNeighbor(Moose::NeighborNeighbor, i, j, LocalDataKey{}, tag),
3999 : *ivar,
4000 : *jvar,
4001 84646 : ivar->dofIndicesNeighbor(),
4002 84646 : jvar->dofIndicesNeighbor());
4003 : }
4004 : }
4005 112176 : }
4006 :
4007 : void
4008 4888 : Assembly::addJacobianLowerD(GlobalDataKey)
4009 : {
4010 42320 : for (const auto & it : _cm_ff_entry)
4011 : {
4012 37432 : auto ivar = it.first;
4013 37432 : auto jvar = it.second;
4014 37432 : auto i = ivar->number();
4015 37432 : auto j = jvar->number();
4016 112296 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
4017 : tag++)
4018 74864 : if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
4019 : {
4020 6704 : addJacobianBlock(_sys.getMatrix(tag),
4021 6704 : jacobianBlockMortar(Moose::LowerLower, i, j, LocalDataKey{}, tag),
4022 : *ivar,
4023 : *jvar,
4024 6704 : ivar->dofIndicesLower(),
4025 6704 : jvar->dofIndicesLower());
4026 :
4027 6704 : addJacobianBlock(_sys.getMatrix(tag),
4028 6704 : jacobianBlockMortar(Moose::LowerSecondary, i, j, LocalDataKey{}, tag),
4029 : *ivar,
4030 : *jvar,
4031 6704 : ivar->dofIndicesLower(),
4032 6704 : jvar->dofIndices());
4033 :
4034 6704 : addJacobianBlock(_sys.getMatrix(tag),
4035 13408 : jacobianBlockMortar(Moose::SecondaryLower, i, j, LocalDataKey{}, tag),
4036 : *ivar,
4037 : *jvar,
4038 6704 : ivar->dofIndices(),
4039 6704 : jvar->dofIndicesLower());
4040 : }
4041 : }
4042 4888 : }
4043 :
4044 : void
4045 48176134 : Assembly::cacheJacobian(GlobalDataKey)
4046 : {
4047 118591179 : for (const auto & it : _cm_ff_entry)
4048 70415045 : cacheJacobianCoupledVarPair(*it.first, *it.second);
4049 :
4050 48434874 : for (const auto & it : _cm_fs_entry)
4051 258740 : cacheJacobianCoupledVarPair(*it.first, *it.second);
4052 :
4053 48434874 : for (const auto & it : _cm_sf_entry)
4054 258740 : cacheJacobianCoupledVarPair(*it.first, *it.second);
4055 48176134 : }
4056 :
4057 : // private method, so no key required
4058 : void
4059 70932525 : Assembly::cacheJacobianCoupledVarPair(const MooseVariableBase & ivar,
4060 : const MooseVariableBase & jvar)
4061 : {
4062 70932525 : auto i = ivar.number();
4063 70932525 : auto j = jvar.number();
4064 214786131 : for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
4065 143853606 : if (jacobianBlockUsed(tag, i, j) && _sys.hasMatrix(tag))
4066 57908442 : cacheJacobianBlock(jacobianBlock(i, j, LocalDataKey{}, tag),
4067 : ivar,
4068 : jvar,
4069 57908442 : ivar.dofIndices(),
4070 57908442 : jvar.dofIndices(),
4071 : tag);
4072 70932525 : }
4073 :
4074 : void
4075 4404 : Assembly::cacheJacobianNonlocal(GlobalDataKey)
4076 : {
4077 8912 : for (const auto & it : _cm_nonlocal_entry)
4078 : {
4079 4508 : auto ivar = it.first;
4080 4508 : auto jvar = it.second;
4081 4508 : auto i = ivar->number();
4082 4508 : auto j = jvar->number();
4083 13524 : for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
4084 13524 : tag < _jacobian_block_nonlocal_used.size();
4085 : tag++)
4086 9016 : if (jacobianBlockNonlocalUsed(tag, i, j) && _sys.hasMatrix(tag))
4087 1268 : cacheJacobianBlockNonzero(jacobianBlockNonlocal(i, j, LocalDataKey{}, tag),
4088 : *ivar,
4089 : *jvar,
4090 634 : ivar->dofIndices(),
4091 : jvar->allDofIndices(),
4092 : tag);
4093 : }
4094 4404 : }
4095 :
4096 : void
4097 9409 : Assembly::cacheJacobianNeighbor(GlobalDataKey)
4098 : {
4099 32642 : for (const auto & it : _cm_ff_entry)
4100 : {
4101 23233 : auto ivar = it.first;
4102 23233 : auto jvar = it.second;
4103 23233 : auto i = ivar->number();
4104 23233 : auto j = jvar->number();
4105 :
4106 69699 : for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
4107 69699 : tag < _jacobian_block_neighbor_used.size();
4108 : tag++)
4109 46466 : if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
4110 : {
4111 7717 : cacheJacobianBlock(jacobianBlockNeighbor(Moose::ElementNeighbor, i, j, LocalDataKey{}, tag),
4112 : *ivar,
4113 : *jvar,
4114 7717 : ivar->dofIndices(),
4115 7717 : jvar->dofIndicesNeighbor(),
4116 : tag);
4117 7717 : cacheJacobianBlock(jacobianBlockNeighbor(Moose::NeighborElement, i, j, LocalDataKey{}, tag),
4118 : *ivar,
4119 : *jvar,
4120 7717 : ivar->dofIndicesNeighbor(),
4121 7717 : jvar->dofIndices(),
4122 : tag);
4123 30868 : cacheJacobianBlock(
4124 7717 : jacobianBlockNeighbor(Moose::NeighborNeighbor, i, j, LocalDataKey{}, tag),
4125 : *ivar,
4126 : *jvar,
4127 7717 : ivar->dofIndicesNeighbor(),
4128 7717 : jvar->dofIndicesNeighbor(),
4129 : tag);
4130 : }
4131 : }
4132 9409 : }
4133 :
4134 : void
4135 191698 : Assembly::cacheJacobianMortar(GlobalDataKey)
4136 : {
4137 886858 : for (const auto & it : _cm_ff_entry)
4138 : {
4139 695160 : auto ivar = it.first;
4140 695160 : auto jvar = it.second;
4141 695160 : auto i = ivar->number();
4142 695160 : auto j = jvar->number();
4143 2085480 : for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
4144 : tag++)
4145 1390320 : if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
4146 : {
4147 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::LowerLower, i, j, LocalDataKey{}, tag),
4148 : *ivar,
4149 : *jvar,
4150 44872 : ivar->dofIndicesLower(),
4151 44872 : jvar->dofIndicesLower(),
4152 : tag);
4153 :
4154 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::LowerSecondary, i, j, LocalDataKey{}, tag),
4155 : *ivar,
4156 : *jvar,
4157 44872 : ivar->dofIndicesLower(),
4158 44872 : jvar->dofIndices(),
4159 : tag);
4160 :
4161 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::LowerPrimary, i, j, LocalDataKey{}, tag),
4162 : *ivar,
4163 : *jvar,
4164 44872 : ivar->dofIndicesLower(),
4165 44872 : jvar->dofIndicesNeighbor(),
4166 : tag);
4167 :
4168 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::SecondaryLower, i, j, LocalDataKey{}, tag),
4169 : *ivar,
4170 : *jvar,
4171 44872 : ivar->dofIndices(),
4172 44872 : jvar->dofIndicesLower(),
4173 : tag);
4174 :
4175 179488 : cacheJacobianBlock(
4176 44872 : jacobianBlockMortar(Moose::SecondarySecondary, i, j, LocalDataKey{}, tag),
4177 : *ivar,
4178 : *jvar,
4179 44872 : ivar->dofIndices(),
4180 44872 : jvar->dofIndices(),
4181 : tag);
4182 :
4183 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::SecondaryPrimary, i, j, LocalDataKey{}, tag),
4184 : *ivar,
4185 : *jvar,
4186 44872 : ivar->dofIndices(),
4187 44872 : jvar->dofIndicesNeighbor(),
4188 : tag);
4189 :
4190 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::PrimaryLower, i, j, LocalDataKey{}, tag),
4191 : *ivar,
4192 : *jvar,
4193 44872 : ivar->dofIndicesNeighbor(),
4194 44872 : jvar->dofIndicesLower(),
4195 : tag);
4196 :
4197 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::PrimarySecondary, i, j, LocalDataKey{}, tag),
4198 : *ivar,
4199 : *jvar,
4200 44872 : ivar->dofIndicesNeighbor(),
4201 44872 : jvar->dofIndices(),
4202 : tag);
4203 :
4204 44872 : cacheJacobianBlock(jacobianBlockMortar(Moose::PrimaryPrimary, i, j, LocalDataKey{}, tag),
4205 : *ivar,
4206 : *jvar,
4207 44872 : ivar->dofIndicesNeighbor(),
4208 44872 : jvar->dofIndicesNeighbor(),
4209 : tag);
4210 : }
4211 : }
4212 191698 : }
4213 :
4214 : void
4215 70832 : Assembly::addJacobianBlockTags(SparseMatrix<Number> & jacobian,
4216 : unsigned int ivar,
4217 : unsigned int jvar,
4218 : const DofMap & dof_map,
4219 : std::vector<dof_id_type> & dof_indices,
4220 : GlobalDataKey,
4221 : const std::set<TagID> & tags)
4222 : {
4223 184848 : for (auto tag : tags)
4224 114016 : addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices, GlobalDataKey{}, tag);
4225 70832 : }
4226 :
4227 : void
4228 114016 : Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
4229 : unsigned int ivar,
4230 : unsigned int jvar,
4231 : const DofMap & dof_map,
4232 : std::vector<dof_id_type> & dof_indices,
4233 : GlobalDataKey,
4234 : TagID tag)
4235 : {
4236 114016 : if (dof_indices.size() == 0)
4237 0 : return;
4238 114016 : if (!(*_cm)(ivar, jvar))
4239 0 : return;
4240 :
4241 114016 : auto & iv = _sys.getVariable(_tid, ivar);
4242 114016 : auto & jv = _sys.getVariable(_tid, jvar);
4243 114016 : auto & scaling_factor = iv.arrayScalingFactor();
4244 :
4245 114016 : const unsigned int ivn = iv.number();
4246 114016 : const unsigned int jvn = jv.number();
4247 114016 : auto & ke = jacobianBlock(ivn, jvn, LocalDataKey{}, tag);
4248 :
4249 : // It is guaranteed by design iv.number <= ivar since iv is obtained
4250 : // through SystemBase::getVariable with ivar.
4251 : // Most of times ivar will just be equal to iv.number except for array variables,
4252 : // where ivar could be a number for a component of an array variable but calling
4253 : // getVariable will return the array variable that has the number of the 0th component.
4254 : // It is the same for jvar.
4255 114016 : const unsigned int i = ivar - ivn;
4256 114016 : const unsigned int j = jvar - jvn;
4257 :
4258 : // DoF indices are independently given
4259 114016 : auto di = dof_indices;
4260 114016 : auto dj = dof_indices;
4261 :
4262 114016 : auto indof = di.size();
4263 114016 : auto jndof = dj.size();
4264 :
4265 114016 : unsigned int jj = j;
4266 114016 : if (ivar == jvar && _component_block_diagonal[ivn])
4267 109568 : jj = 0;
4268 :
4269 114016 : auto sub = ke.sub_matrix(i * indof, indof, jj * jndof, jndof);
4270 : // If we're computing the jacobian for automatically scaling variables we do not want to
4271 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4272 : // dofs
4273 114016 : if (!_sys.computingScalingJacobian())
4274 114016 : dof_map.constrain_element_matrix(sub, di, dj, false);
4275 :
4276 114016 : if (scaling_factor[i] != 1.0)
4277 0 : sub *= scaling_factor[i];
4278 :
4279 114016 : jacobian.add_matrix(sub, di, dj);
4280 114016 : }
4281 :
4282 : void
4283 0 : Assembly::addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
4284 : const unsigned int ivar,
4285 : const unsigned int jvar,
4286 : const DofMap & dof_map,
4287 : const std::vector<dof_id_type> & idof_indices,
4288 : const std::vector<dof_id_type> & jdof_indices,
4289 : GlobalDataKey,
4290 : const TagID tag)
4291 : {
4292 0 : if (idof_indices.size() == 0 || jdof_indices.size() == 0)
4293 0 : return;
4294 0 : if (jacobian.n() == 0 || jacobian.m() == 0)
4295 0 : return;
4296 0 : if (!(*_cm)(ivar, jvar))
4297 0 : return;
4298 :
4299 0 : auto & iv = _sys.getVariable(_tid, ivar);
4300 0 : auto & jv = _sys.getVariable(_tid, jvar);
4301 0 : auto & scaling_factor = iv.arrayScalingFactor();
4302 :
4303 0 : const unsigned int ivn = iv.number();
4304 0 : const unsigned int jvn = jv.number();
4305 0 : auto & keg = jacobianBlockNonlocal(ivn, jvn, LocalDataKey{}, tag);
4306 :
4307 : // It is guaranteed by design iv.number <= ivar since iv is obtained
4308 : // through SystemBase::getVariable with ivar.
4309 : // Most of times ivar will just be equal to iv.number except for array variables,
4310 : // where ivar could be a number for a component of an array variable but calling
4311 : // getVariable will return the array variable that has the number of the 0th component.
4312 : // It is the same for jvar.
4313 0 : const unsigned int i = ivar - ivn;
4314 0 : const unsigned int j = jvar - jvn;
4315 :
4316 : // DoF indices are independently given
4317 0 : auto di = idof_indices;
4318 0 : auto dj = jdof_indices;
4319 :
4320 0 : auto indof = di.size();
4321 0 : auto jndof = dj.size();
4322 :
4323 0 : unsigned int jj = j;
4324 0 : if (ivar == jvar && _component_block_diagonal[ivn])
4325 0 : jj = 0;
4326 :
4327 0 : auto sub = keg.sub_matrix(i * indof, indof, jj * jndof, jndof);
4328 : // If we're computing the jacobian for automatically scaling variables we do not want to
4329 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4330 : // dofs
4331 0 : if (!_sys.computingScalingJacobian())
4332 0 : dof_map.constrain_element_matrix(sub, di, dj, false);
4333 :
4334 0 : if (scaling_factor[i] != 1.0)
4335 0 : sub *= scaling_factor[i];
4336 :
4337 0 : jacobian.add_matrix(sub, di, dj);
4338 0 : }
4339 :
4340 : void
4341 0 : Assembly::addJacobianBlockNonlocalTags(SparseMatrix<Number> & jacobian,
4342 : const unsigned int ivar,
4343 : const unsigned int jvar,
4344 : const DofMap & dof_map,
4345 : const std::vector<dof_id_type> & idof_indices,
4346 : const std::vector<dof_id_type> & jdof_indices,
4347 : GlobalDataKey,
4348 : const std::set<TagID> & tags)
4349 : {
4350 0 : for (auto tag : tags)
4351 0 : addJacobianBlockNonlocal(
4352 0 : jacobian, ivar, jvar, dof_map, idof_indices, jdof_indices, GlobalDataKey{}, tag);
4353 0 : }
4354 :
4355 : void
4356 1536 : Assembly::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
4357 : const unsigned int ivar,
4358 : const unsigned int jvar,
4359 : const DofMap & dof_map,
4360 : std::vector<dof_id_type> & dof_indices,
4361 : std::vector<dof_id_type> & neighbor_dof_indices,
4362 : GlobalDataKey,
4363 : const TagID tag)
4364 : {
4365 1536 : if (dof_indices.size() == 0 && neighbor_dof_indices.size() == 0)
4366 0 : return;
4367 1536 : if (!(*_cm)(ivar, jvar))
4368 0 : return;
4369 :
4370 1536 : auto & iv = _sys.getVariable(_tid, ivar);
4371 1536 : auto & jv = _sys.getVariable(_tid, jvar);
4372 1536 : auto & scaling_factor = iv.arrayScalingFactor();
4373 :
4374 1536 : const unsigned int ivn = iv.number();
4375 1536 : const unsigned int jvn = jv.number();
4376 1536 : auto & ken = jacobianBlockNeighbor(Moose::ElementNeighbor, ivn, jvn, LocalDataKey{}, tag);
4377 1536 : auto & kne = jacobianBlockNeighbor(Moose::NeighborElement, ivn, jvn, LocalDataKey{}, tag);
4378 1536 : auto & knn = jacobianBlockNeighbor(Moose::NeighborNeighbor, ivn, jvn, LocalDataKey{}, tag);
4379 :
4380 : // It is guaranteed by design iv.number <= ivar since iv is obtained
4381 : // through SystemBase::getVariable with ivar.
4382 : // Most of times ivar will just be equal to iv.number except for array variables,
4383 : // where ivar could be a number for a component of an array variable but calling
4384 : // getVariable will return the array variable that has the number of the 0th component.
4385 : // It is the same for jvar.
4386 1536 : const unsigned int i = ivar - ivn;
4387 1536 : const unsigned int j = jvar - jvn;
4388 : // DoF indices are independently given
4389 1536 : auto dc = dof_indices;
4390 1536 : auto dn = neighbor_dof_indices;
4391 1536 : auto cndof = dc.size();
4392 1536 : auto nndof = dn.size();
4393 :
4394 1536 : unsigned int jj = j;
4395 1536 : if (ivar == jvar && _component_block_diagonal[ivn])
4396 1536 : jj = 0;
4397 :
4398 1536 : auto suben = ken.sub_matrix(i * cndof, cndof, jj * nndof, nndof);
4399 1536 : auto subne = kne.sub_matrix(i * nndof, nndof, jj * cndof, cndof);
4400 1536 : auto subnn = knn.sub_matrix(i * nndof, nndof, jj * nndof, nndof);
4401 :
4402 : // If we're computing the jacobian for automatically scaling variables we do not want to
4403 : // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4404 : // dofs
4405 1536 : if (!_sys.computingScalingJacobian())
4406 : {
4407 1536 : dof_map.constrain_element_matrix(suben, dc, dn, false);
4408 1536 : dof_map.constrain_element_matrix(subne, dn, dc, false);
4409 1536 : dof_map.constrain_element_matrix(subnn, dn, dn, false);
4410 : }
4411 :
4412 1536 : if (scaling_factor[i] != 1.0)
4413 : {
4414 0 : suben *= scaling_factor[i];
4415 0 : subne *= scaling_factor[i];
4416 0 : subnn *= scaling_factor[i];
4417 : }
4418 :
4419 1536 : jacobian.add_matrix(suben, dc, dn);
4420 1536 : jacobian.add_matrix(subne, dn, dc);
4421 1536 : jacobian.add_matrix(subnn, dn, dn);
4422 1536 : }
4423 :
4424 : void
4425 768 : Assembly::addJacobianNeighborTags(SparseMatrix<Number> & jacobian,
4426 : const unsigned int ivar,
4427 : const unsigned int jvar,
4428 : const DofMap & dof_map,
4429 : std::vector<dof_id_type> & dof_indices,
4430 : std::vector<dof_id_type> & neighbor_dof_indices,
4431 : GlobalDataKey,
4432 : const std::set<TagID> & tags)
4433 : {
4434 2304 : for (const auto tag : tags)
4435 1536 : addJacobianNeighbor(
4436 3072 : jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices, GlobalDataKey{}, tag);
4437 768 : }
4438 :
4439 : void
4440 11603 : Assembly::addJacobianScalar(GlobalDataKey)
4441 : {
4442 26924 : for (const auto & it : _cm_ss_entry)
4443 15321 : addJacobianCoupledVarPair(*it.first, *it.second);
4444 11603 : }
4445 :
4446 : void
4447 30002 : Assembly::addJacobianOffDiagScalar(unsigned int ivar, GlobalDataKey)
4448 : {
4449 30002 : const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
4450 30002 : MooseVariableScalar & var_i = _sys.getScalarVariable(_tid, ivar);
4451 42929 : for (const auto & var_j : vars)
4452 12927 : addJacobianCoupledVarPair(var_i, *var_j);
4453 30002 : }
4454 :
4455 : void
4456 223982165 : Assembly::cacheJacobian(
4457 : numeric_index_type i, numeric_index_type j, Real value, LocalDataKey, TagID tag)
4458 : {
4459 223982165 : _cached_jacobian_rows[tag].push_back(i);
4460 223982165 : _cached_jacobian_cols[tag].push_back(j);
4461 223982165 : _cached_jacobian_values[tag].push_back(value);
4462 223982165 : }
4463 :
4464 : void
4465 223928421 : Assembly::cacheJacobian(numeric_index_type i,
4466 : numeric_index_type j,
4467 : Real value,
4468 : LocalDataKey,
4469 : const std::set<TagID> & tags)
4470 : {
4471 463141763 : for (auto tag : tags)
4472 239213342 : if (_sys.hasMatrix(tag))
4473 223982165 : cacheJacobian(i, j, value, LocalDataKey{}, tag);
4474 223928421 : }
4475 :
4476 : void
4477 397598 : Assembly::setCachedJacobian(GlobalDataKey)
4478 : {
4479 1204621 : for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4480 807023 : if (_sys.hasMatrix(tag))
4481 : {
4482 : // First zero the rows (including the diagonals) to prepare for
4483 : // setting the cached values.
4484 398999 : _sys.getMatrix(tag).zero_rows(_cached_jacobian_rows[tag], 0.0);
4485 :
4486 : // TODO: Use SparseMatrix::set_values() for efficiency
4487 8649349 : for (MooseIndex(_cached_jacobian_values) i = 0; i < _cached_jacobian_values[tag].size(); ++i)
4488 16500700 : _sys.getMatrix(tag).set(_cached_jacobian_rows[tag][i],
4489 8250350 : _cached_jacobian_cols[tag][i],
4490 8250350 : _cached_jacobian_values[tag][i]);
4491 : }
4492 :
4493 397598 : clearCachedJacobian();
4494 397598 : }
4495 :
4496 : void
4497 0 : Assembly::zeroCachedJacobian(GlobalDataKey)
4498 : {
4499 0 : for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4500 0 : if (_sys.hasMatrix(tag))
4501 0 : _sys.getMatrix(tag).zero_rows(_cached_jacobian_rows[tag], 0.0);
4502 :
4503 0 : clearCachedJacobian();
4504 0 : }
4505 :
4506 : void
4507 397598 : Assembly::clearCachedJacobian()
4508 : {
4509 1204621 : for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4510 : {
4511 807023 : _cached_jacobian_rows[tag].clear();
4512 807023 : _cached_jacobian_cols[tag].clear();
4513 807023 : _cached_jacobian_values[tag].clear();
4514 : }
4515 397598 : }
4516 :
4517 : void
4518 0 : Assembly::modifyWeightsDueToXFEM(const Elem * elem)
4519 : {
4520 : mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
4521 :
4522 0 : if (_current_qrule == _current_qrule_arbitrary)
4523 0 : return;
4524 :
4525 0 : MooseArray<Real> xfem_weight_multipliers;
4526 0 : if (_xfem->getXFEMWeights(xfem_weight_multipliers, elem, _current_qrule, _current_q_points))
4527 : {
4528 : mooseAssert(xfem_weight_multipliers.size() == _current_JxW.size(),
4529 : "Size of weight multipliers in xfem doesn't match number of quadrature points");
4530 0 : for (unsigned i = 0; i < xfem_weight_multipliers.size(); i++)
4531 0 : _current_JxW[i] = _current_JxW[i] * xfem_weight_multipliers[i];
4532 :
4533 0 : xfem_weight_multipliers.release();
4534 : }
4535 0 : }
4536 :
4537 : void
4538 0 : Assembly::modifyFaceWeightsDueToXFEM(const Elem * elem, unsigned int side)
4539 : {
4540 : mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
4541 :
4542 0 : if (_current_qrule_face == _current_qrule_arbitrary)
4543 0 : return;
4544 :
4545 0 : MooseArray<Real> xfem_face_weight_multipliers;
4546 0 : if (_xfem->getXFEMFaceWeights(
4547 0 : xfem_face_weight_multipliers, elem, _current_qrule_face, _current_q_points_face, side))
4548 : {
4549 : mooseAssert(xfem_face_weight_multipliers.size() == _current_JxW_face.size(),
4550 : "Size of weight multipliers in xfem doesn't match number of quadrature points");
4551 0 : for (unsigned i = 0; i < xfem_face_weight_multipliers.size(); i++)
4552 0 : _current_JxW_face[i] = _current_JxW_face[i] * xfem_face_weight_multipliers[i];
4553 :
4554 0 : xfem_face_weight_multipliers.release();
4555 : }
4556 0 : }
4557 :
4558 : void
4559 450 : Assembly::hasScalingVector()
4560 : {
4561 900 : _scaling_vector = &_sys.getVector("scaling_factors");
4562 450 : }
4563 :
4564 : void
4565 0 : Assembly::modifyArbitraryWeights(const std::vector<Real> & weights)
4566 : {
4567 : mooseAssert(_current_qrule == _current_qrule_arbitrary, "Rule should be arbitrary");
4568 : mooseAssert(weights.size() == _current_physical_points.size(), "Size mismatch");
4569 :
4570 0 : for (MooseIndex(weights.size()) i = 0; i < weights.size(); ++i)
4571 0 : _current_JxW[i] = weights[i];
4572 0 : }
4573 :
4574 : template <>
4575 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4576 1870 : Assembly::fePhi<VectorValue<Real>>(FEType type) const
4577 : {
4578 1870 : buildVectorFE(type);
4579 1870 : return _vector_fe_shape_data[type]->_phi;
4580 : }
4581 :
4582 : template <>
4583 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4584 1870 : Assembly::feGradPhi<VectorValue<Real>>(FEType type) const
4585 : {
4586 1870 : buildVectorFE(type);
4587 1870 : return _vector_fe_shape_data[type]->_grad_phi;
4588 : }
4589 :
4590 : template <>
4591 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4592 0 : Assembly::feSecondPhi<VectorValue<Real>>(FEType type) const
4593 : {
4594 0 : _need_second_derivative.insert(type);
4595 0 : buildVectorFE(type);
4596 0 : return _vector_fe_shape_data[type]->_second_phi;
4597 : }
4598 :
4599 : template <>
4600 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4601 3740 : Assembly::fePhiLower<VectorValue<Real>>(FEType type) const
4602 : {
4603 3740 : buildVectorLowerDFE(type);
4604 3740 : return _vector_fe_shape_data_lower[type]->_phi;
4605 : }
4606 :
4607 : template <>
4608 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4609 0 : Assembly::feDualPhiLower<VectorValue<Real>>(FEType type) const
4610 : {
4611 0 : buildVectorDualLowerDFE(type);
4612 0 : return _vector_fe_shape_data_dual_lower[type]->_phi;
4613 : }
4614 :
4615 : template <>
4616 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4617 3740 : Assembly::feGradPhiLower<VectorValue<Real>>(FEType type) const
4618 : {
4619 3740 : buildVectorLowerDFE(type);
4620 3740 : return _vector_fe_shape_data_lower[type]->_grad_phi;
4621 : }
4622 :
4623 : template <>
4624 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4625 0 : Assembly::feGradDualPhiLower<VectorValue<Real>>(FEType type) const
4626 : {
4627 0 : buildVectorDualLowerDFE(type);
4628 0 : return _vector_fe_shape_data_dual_lower[type]->_grad_phi;
4629 : }
4630 :
4631 : template <>
4632 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4633 1870 : Assembly::fePhiFace<VectorValue<Real>>(FEType type) const
4634 : {
4635 1870 : buildVectorFaceFE(type);
4636 1870 : return _vector_fe_shape_data_face[type]->_phi;
4637 : }
4638 :
4639 : template <>
4640 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4641 1870 : Assembly::feGradPhiFace<VectorValue<Real>>(FEType type) const
4642 : {
4643 1870 : buildVectorFaceFE(type);
4644 1870 : return _vector_fe_shape_data_face[type]->_grad_phi;
4645 : }
4646 :
4647 : template <>
4648 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4649 0 : Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type) const
4650 : {
4651 0 : _need_second_derivative.insert(type);
4652 0 : buildVectorFaceFE(type);
4653 :
4654 : // If we're building for a face we probably need to build for a
4655 : // neighbor while _need_second_derivative is set;
4656 : // onInterface/reinitNeighbor/etc don't distinguish
4657 0 : buildVectorFaceNeighborFE(type);
4658 :
4659 0 : return _vector_fe_shape_data_face[type]->_second_phi;
4660 : }
4661 :
4662 : template <>
4663 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4664 1870 : Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type) const
4665 : {
4666 1870 : buildVectorNeighborFE(type);
4667 1870 : return _vector_fe_shape_data_neighbor[type]->_phi;
4668 : }
4669 :
4670 : template <>
4671 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4672 1870 : Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type) const
4673 : {
4674 1870 : buildVectorNeighborFE(type);
4675 1870 : return _vector_fe_shape_data_neighbor[type]->_grad_phi;
4676 : }
4677 :
4678 : template <>
4679 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4680 0 : Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type) const
4681 : {
4682 0 : _need_second_derivative_neighbor.insert(type);
4683 0 : buildVectorNeighborFE(type);
4684 0 : return _vector_fe_shape_data_neighbor[type]->_second_phi;
4685 : }
4686 :
4687 : template <>
4688 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
4689 1870 : Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4690 : {
4691 1870 : buildVectorFaceNeighborFE(type);
4692 1870 : return _vector_fe_shape_data_face_neighbor[type]->_phi;
4693 : }
4694 :
4695 : template <>
4696 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
4697 1870 : Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4698 : {
4699 1870 : buildVectorFaceNeighborFE(type);
4700 1870 : return _vector_fe_shape_data_face_neighbor[type]->_grad_phi;
4701 : }
4702 :
4703 : template <>
4704 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
4705 0 : Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4706 : {
4707 0 : _need_second_derivative_neighbor.insert(type);
4708 0 : buildVectorFaceNeighborFE(type);
4709 0 : return _vector_fe_shape_data_face_neighbor[type]->_second_phi;
4710 : }
4711 :
4712 : template <>
4713 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4714 57237 : Assembly::feCurlPhi<VectorValue<Real>>(FEType type) const
4715 : {
4716 57237 : _need_curl.insert(type);
4717 57237 : buildVectorFE(type);
4718 57237 : return _vector_fe_shape_data[type]->_curl_phi;
4719 : }
4720 :
4721 : template <>
4722 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4723 207 : Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type) const
4724 : {
4725 207 : _need_curl.insert(type);
4726 207 : buildVectorFaceFE(type);
4727 :
4728 : // If we're building for a face we probably need to build for a
4729 : // neighbor while _need_curl is set;
4730 : // onInterface/reinitNeighbor/etc don't distinguish
4731 207 : buildVectorFaceNeighborFE(type);
4732 :
4733 207 : return _vector_fe_shape_data_face[type]->_curl_phi;
4734 : }
4735 :
4736 : template <>
4737 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4738 0 : Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type) const
4739 : {
4740 0 : _need_curl.insert(type);
4741 0 : buildVectorNeighborFE(type);
4742 0 : return _vector_fe_shape_data_neighbor[type]->_curl_phi;
4743 : }
4744 :
4745 : template <>
4746 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
4747 0 : Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4748 : {
4749 0 : _need_curl.insert(type);
4750 0 : buildVectorFaceNeighborFE(type);
4751 :
4752 0 : return _vector_fe_shape_data_face_neighbor[type]->_curl_phi;
4753 : }
4754 :
4755 : template <>
4756 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4757 521416 : Assembly::feDivPhi<VectorValue<Real>>(FEType type) const
4758 : {
4759 521416 : _need_div.insert(type);
4760 521416 : buildVectorFE(type);
4761 521416 : return _vector_fe_shape_data[type]->_div_phi;
4762 : }
4763 :
4764 : template <>
4765 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4766 546 : Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const
4767 : {
4768 546 : _need_face_div.insert(type);
4769 546 : buildVectorFaceFE(type);
4770 :
4771 : // If we're building for a face we probably need to build for a
4772 : // neighbor while _need_face_div is set;
4773 : // onInterface/reinitNeighbor/etc don't distinguish
4774 546 : buildVectorFaceNeighborFE(type);
4775 :
4776 546 : return _vector_fe_shape_data_face[type]->_div_phi;
4777 : }
4778 :
4779 : template <>
4780 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4781 0 : Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const
4782 : {
4783 0 : _need_neighbor_div.insert(type);
4784 0 : buildVectorNeighborFE(type);
4785 0 : return _vector_fe_shape_data_neighbor[type]->_div_phi;
4786 : }
4787 :
4788 : template <>
4789 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
4790 0 : Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4791 : {
4792 0 : _need_face_neighbor_div.insert(type);
4793 0 : buildVectorFaceNeighborFE(type);
4794 0 : return _vector_fe_shape_data_face_neighbor[type]->_div_phi;
4795 : }
4796 :
4797 : const MooseArray<ADReal> &
4798 26 : Assembly::adCurvatures() const
4799 : {
4800 26 : _calculate_curvatures = true;
4801 26 : const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
4802 26 : const FEType helper_type(helper_order, LAGRANGE);
4803 : // Must prerequest the second derivatives. Sadly because there is only one
4804 : // _need_second_derivative map for both volumetric and face FE objects we must request both here
4805 26 : feSecondPhi<Real>(helper_type);
4806 26 : feSecondPhiFace<Real>(helper_type);
4807 26 : return _ad_curvatures;
4808 : }
4809 :
4810 : void
4811 73112 : Assembly::helpersRequestData()
4812 : {
4813 286316 : for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
4814 : {
4815 213204 : _holder_fe_helper[dim]->get_phi();
4816 213204 : _holder_fe_helper[dim]->get_dphi();
4817 213204 : _holder_fe_helper[dim]->get_xyz();
4818 213204 : _holder_fe_helper[dim]->get_JxW();
4819 :
4820 213204 : _holder_fe_face_helper[dim]->get_phi();
4821 213204 : _holder_fe_face_helper[dim]->get_dphi();
4822 213204 : _holder_fe_face_helper[dim]->get_xyz();
4823 213204 : _holder_fe_face_helper[dim]->get_JxW();
4824 213204 : _holder_fe_face_helper[dim]->get_normals();
4825 :
4826 213204 : _holder_fe_face_neighbor_helper[dim]->get_xyz();
4827 213204 : _holder_fe_face_neighbor_helper[dim]->get_JxW();
4828 213204 : _holder_fe_face_neighbor_helper[dim]->get_normals();
4829 :
4830 213204 : _holder_fe_neighbor_helper[dim]->get_xyz();
4831 213204 : _holder_fe_neighbor_helper[dim]->get_JxW();
4832 : }
4833 :
4834 213204 : for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
4835 : {
4836 : // We need these computations in order to compute correct lower-d element volumes in
4837 : // curvilinear coordinates
4838 140092 : _holder_fe_lower_helper[dim]->get_xyz();
4839 140092 : _holder_fe_lower_helper[dim]->get_JxW();
4840 : }
4841 73112 : }
4842 :
4843 : void
4844 260 : Assembly::havePRefinement(const std::unordered_set<FEFamily> & disable_families)
4845 : {
4846 260 : if (_have_p_refinement)
4847 : // Already performed tasks for p-refinement
4848 0 : return;
4849 :
4850 260 : const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
4851 260 : const FEType helper_type(helper_order, LAGRANGE);
4852 : auto process_fe =
4853 2600 : [&disable_families](const unsigned int num_dimensionalities, auto & fe_container)
4854 : {
4855 2600 : if (!disable_families.empty())
4856 8034 : for (const auto dim : make_range(num_dimensionalities))
4857 : {
4858 5954 : auto fe_container_it = fe_container.find(dim);
4859 5954 : if (fe_container_it != fe_container.end())
4860 11897 : for (auto & [fe_type, fe_ptr] : fe_container_it->second)
4861 8008 : if (disable_families.count(fe_type.family))
4862 2847 : fe_ptr->add_p_level_in_reinit(false);
4863 : }
4864 2860 : };
4865 1300 : auto process_fe_and_helpers = [process_fe, &helper_type](auto & unique_helper_container,
4866 : auto & helper_container,
4867 : const unsigned int num_dimensionalities,
4868 : const bool user_added_helper_type,
4869 : auto & fe_container)
4870 : {
4871 1300 : unique_helper_container.resize(num_dimensionalities);
4872 5005 : for (const auto dim : make_range(num_dimensionalities))
4873 : {
4874 3705 : auto & unique_helper = unique_helper_container[dim];
4875 3705 : unique_helper = FEGenericBase<Real>::build(dim, helper_type);
4876 : // don't participate in p-refinement
4877 3705 : unique_helper->add_p_level_in_reinit(false);
4878 3705 : helper_container[dim] = unique_helper.get();
4879 :
4880 : // If the user did not request the helper type then we should erase it from our FE container
4881 : // so that they're not penalized (in the "we should be able to do p-refinement sense") for
4882 : // our perhaps silly helpers
4883 3705 : if (!user_added_helper_type)
4884 : {
4885 2730 : auto & fe_container_dim = libmesh_map_find(fe_container, dim);
4886 2730 : auto fe_it = fe_container_dim.find(helper_type);
4887 : mooseAssert(fe_it != fe_container_dim.end(), "We should have the helper type");
4888 2730 : delete fe_it->second;
4889 2730 : fe_container_dim.erase(fe_it);
4890 : }
4891 : }
4892 :
4893 1300 : process_fe(num_dimensionalities, fe_container);
4894 1300 : };
4895 :
4896 : // Handle scalar field families
4897 260 : process_fe_and_helpers(_unique_fe_helper,
4898 260 : _holder_fe_helper,
4899 260 : _mesh_dimension + 1,
4900 260 : _user_added_fe_of_helper_type,
4901 260 : _fe);
4902 260 : process_fe_and_helpers(_unique_fe_face_helper,
4903 260 : _holder_fe_face_helper,
4904 260 : _mesh_dimension + 1,
4905 260 : _user_added_fe_face_of_helper_type,
4906 260 : _fe_face);
4907 260 : process_fe_and_helpers(_unique_fe_face_neighbor_helper,
4908 260 : _holder_fe_face_neighbor_helper,
4909 260 : _mesh_dimension + 1,
4910 260 : _user_added_fe_face_neighbor_of_helper_type,
4911 260 : _fe_face_neighbor);
4912 260 : process_fe_and_helpers(_unique_fe_neighbor_helper,
4913 260 : _holder_fe_neighbor_helper,
4914 260 : _mesh_dimension + 1,
4915 260 : _user_added_fe_neighbor_of_helper_type,
4916 260 : _fe_neighbor);
4917 260 : process_fe_and_helpers(_unique_fe_lower_helper,
4918 260 : _holder_fe_lower_helper,
4919 : _mesh_dimension,
4920 260 : _user_added_fe_lower_of_helper_type,
4921 260 : _fe_lower);
4922 : // Handle vector field families
4923 260 : process_fe(_mesh_dimension + 1, _vector_fe);
4924 260 : process_fe(_mesh_dimension + 1, _vector_fe_face);
4925 260 : process_fe(_mesh_dimension + 1, _vector_fe_neighbor);
4926 260 : process_fe(_mesh_dimension + 1, _vector_fe_face_neighbor);
4927 260 : process_fe(_mesh_dimension, _vector_fe_lower);
4928 :
4929 260 : helpersRequestData();
4930 :
4931 260 : _have_p_refinement = true;
4932 : }
4933 :
4934 : template void coordTransformFactor<Point, Real>(const SubProblem & s,
4935 : SubdomainID sub_id,
4936 : const Point & point,
4937 : Real & factor,
4938 : SubdomainID neighbor_sub_id);
4939 : template void coordTransformFactor<ADPoint, ADReal>(const SubProblem & s,
4940 : SubdomainID sub_id,
4941 : const ADPoint & point,
4942 : ADReal & factor,
4943 : SubdomainID neighbor_sub_id);
4944 : template void coordTransformFactor<Point, Real>(const MooseMesh & mesh,
4945 : SubdomainID sub_id,
4946 : const Point & point,
4947 : Real & factor,
4948 : SubdomainID neighbor_sub_id);
4949 : template void coordTransformFactor<ADPoint, ADReal>(const MooseMesh & mesh,
4950 : SubdomainID sub_id,
4951 : const ADPoint & point,
4952 : ADReal & factor,
4953 : SubdomainID neighbor_sub_id);
4954 :
4955 : template <>
4956 : const MooseArray<Moose::GenericType<Point, false>> &
4957 34 : Assembly::genericQPoints<false>() const
4958 : {
4959 34 : return qPoints();
4960 : }
4961 :
4962 : template <>
4963 : const MooseArray<Moose::GenericType<Point, true>> &
4964 10 : Assembly::genericQPoints<true>() const
4965 : {
4966 10 : return adQPoints();
4967 : }
|
