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