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