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