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