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 : #pragma once
11 :
12 : #include "DenseMatrix.h"
13 : #include "MooseArray.h"
14 : #include "MooseTypes.h"
15 : #include "MooseVariableFE.h"
16 : #include "MoosePassKey.h"
17 : #include "ArbitraryQuadrature.h"
18 :
19 : #include "libmesh/dense_vector.h"
20 : #include "libmesh/enum_quadrature_type.h"
21 : #include "libmesh/fe_type.h"
22 : #include "libmesh/point.h"
23 : #include "libmesh/fe_base.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/elem_side_builder.h"
26 :
27 : #include <unordered_map>
28 :
29 : // libMesh forward declarations
30 : namespace libMesh
31 : {
32 : class DofMap;
33 : class CouplingMatrix;
34 : class Elem;
35 : template <typename>
36 : class VectorValue;
37 : typedef VectorValue<Real> RealVectorValue;
38 : template <typename T>
39 : class FEGenericBase;
40 : typedef FEGenericBase<Real> FEBase;
41 : typedef FEGenericBase<VectorValue<Real>> FEVectorBase;
42 : class Node;
43 : template <typename T>
44 : class NumericVector;
45 : template <typename T>
46 : class SparseMatrix;
47 : class StaticCondensation;
48 : }
49 :
50 : // MOOSE Forward Declares
51 : class FaceInfo;
52 : class MooseMesh;
53 : class ArbitraryQuadrature;
54 : class SystemBase;
55 : class MooseVariableFieldBase;
56 : class MooseVariableBase;
57 : template <typename>
58 : class MooseVariableFE;
59 : class MooseVariableScalar;
60 : typedef MooseVariableFE<Real> MooseVariable;
61 : typedef MooseVariableFE<RealVectorValue> VectorMooseVariable;
62 : typedef MooseVariableFE<RealEigenVector> ArrayMooseVariable;
63 : class XFEMInterface;
64 : class SubProblem;
65 : class NodeFaceConstraint;
66 :
67 : #ifdef MOOSE_KOKKOS_ENABLED
68 : namespace Moose::Kokkos
69 : {
70 : class Assembly;
71 : }
72 : #endif
73 :
74 : // Assembly.h does not import Moose.h nor libMeshReducedNamespace.h
75 : using libMesh::FEBase;
76 : using libMesh::FEFamily;
77 : using libMesh::FEType;
78 : using libMesh::FEVectorBase;
79 : using libMesh::LAGRANGE_VEC;
80 : using libMesh::Order;
81 : using libMesh::QuadratureType;
82 :
83 : /// Computes a conversion multiplier for use when computing integraals for the
84 : /// current coordinate system type. This allows us to handle cases where we use RZ,
85 : /// spherical, or other non-cartesian coordinate systems. The factor returned
86 : /// by this function should generally be multiplied against all integration
87 : /// terms. Note that the computed factor is particular to a specific point on
88 : /// the mesh. The result is stored in the factor argument. point is the point
89 : /// at which to compute the factor. point and factor can be either Point and
90 : /// Real or ADPoint and ADReal.
91 : template <typename P, typename C>
92 : void coordTransformFactor(const SubProblem & s,
93 : SubdomainID sub_id,
94 : const P & point,
95 : C & factor,
96 : SubdomainID neighbor_sub_id = libMesh::Elem::invalid_subdomain_id);
97 :
98 : template <typename P, typename C>
99 : void coordTransformFactor(const MooseMesh & mesh,
100 : SubdomainID sub_id,
101 : const P & point,
102 : C & factor,
103 : SubdomainID neighbor_sub_id = libMesh::Elem::invalid_subdomain_id);
104 :
105 : /**
106 : * Keeps track of stuff related to assembling
107 : *
108 : */
109 : class Assembly
110 : {
111 : public:
112 : Assembly(SystemBase & sys, THREAD_ID tid);
113 : virtual ~Assembly();
114 :
115 : /**
116 : * Workaround for C++ compilers thinking they can't just cast a
117 : * const-reference-to-pointer to const-reference-to-const-pointer
118 : */
119 : template <typename T>
120 69876158 : static const T * const & constify_ref(T * const & inref)
121 : {
122 69876158 : const T * const * ptr = &inref;
123 69876158 : return *ptr;
124 : }
125 :
126 : /**
127 : * Get a reference to a pointer that will contain the current volume FE.
128 : * @param type The type of FE
129 : * @param dim The dimension of the current volume
130 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
131 : */
132 4591 : const FEBase * const & getFE(FEType type, unsigned int dim) const
133 : {
134 4591 : buildFE(type);
135 4591 : return constify_ref(_fe[dim][type]);
136 : }
137 :
138 : /**
139 : * Get a reference to a pointer that will contain the current 'neighbor' FE.
140 : * @param type The type of FE
141 : * @param dim The dimension of the current volume
142 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
143 : */
144 : const FEBase * const & getFENeighbor(FEType type, unsigned int dim) const
145 : {
146 : buildNeighborFE(type);
147 : return constify_ref(_fe_neighbor[dim][type]);
148 : }
149 :
150 : /**
151 : * Get a reference to a pointer that will contain the current "face" FE.
152 : * @param type The type of FE
153 : * @param dim The dimension of the current face
154 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
155 : */
156 : const FEBase * const & getFEFace(FEType type, unsigned int dim) const
157 : {
158 : buildFaceFE(type);
159 : return constify_ref(_fe_face[dim][type]);
160 : }
161 :
162 : /**
163 : * Get a reference to a pointer that will contain the current "neighbor" FE.
164 : * @param type The type of FE
165 : * @param dim The dimension of the neighbor face
166 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
167 : */
168 : const FEBase * const & getFEFaceNeighbor(FEType type, unsigned int dim) const
169 : {
170 : buildFaceNeighborFE(type);
171 : return constify_ref(_fe_face_neighbor[dim][type]);
172 : }
173 :
174 : /**
175 : * Get a reference to a pointer that will contain the current volume FEVector.
176 : * @param type The type of FEVector
177 : * @param dim The dimension of the current volume
178 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
179 : */
180 : const FEVectorBase * const & getVectorFE(FEType type, unsigned int dim) const
181 : {
182 : buildVectorFE(type);
183 : return constify_ref(_vector_fe[dim][type]);
184 : }
185 :
186 : /**
187 : * GetVector a reference to a pointer that will contain the current 'neighbor' FE.
188 : * @param type The type of FE
189 : * @param dim The dimension of the current volume
190 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
191 : */
192 : const FEVectorBase * const & getVectorFENeighbor(FEType type, unsigned int dim) const
193 : {
194 : buildVectorNeighborFE(type);
195 : return constify_ref(_vector_fe_neighbor[dim][type]);
196 : }
197 :
198 : /**
199 : * GetVector a reference to a pointer that will contain the current "face" FE.
200 : * @param type The type of FE
201 : * @param dim The dimension of the current face
202 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
203 : */
204 : const FEVectorBase * const & getVectorFEFace(FEType type, unsigned int dim) const
205 : {
206 : buildVectorFaceFE(type);
207 : return constify_ref(_vector_fe_face[dim][type]);
208 : }
209 :
210 : /**
211 : * GetVector a reference to a pointer that will contain the current "neighbor" FE.
212 : * @param type The type of FE
213 : * @param dim The dimension of the neighbor face
214 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
215 : */
216 : const FEVectorBase * const & getVectorFEFaceNeighbor(FEType type, unsigned int dim) const
217 : {
218 : buildVectorFaceNeighborFE(type);
219 : return constify_ref(_vector_fe_face_neighbor[dim][type]);
220 : }
221 :
222 : #ifdef MOOSE_KOKKOS_ENABLED
223 : /**
224 : * Key structure for APIs manipulating internal shape and quadrature data. Developers in blessed
225 : * classes may create keys using simple curly braces \p {} or may be more explicit and use \p
226 : * Assembly::InternalDataKey{}
227 : */
228 : using InternalDataKey = Moose::PassKey<Moose::Kokkos::Assembly>;
229 : #endif
230 :
231 : /**
232 : * Returns the reference to the current quadrature being used
233 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
234 : */
235 33454905 : const libMesh::QBase * const & qRule() const { return constify_ref(_current_qrule); }
236 :
237 : /**
238 : * Returns the reference to the current quadrature being used
239 : * @return A _reference_ to the pointer. Make sure to store this as a reference!
240 : */
241 79116 : libMesh::QBase * const & writeableQRule() { return _current_qrule; }
242 :
243 : #ifdef MOOSE_KOKKOS_ENABLED
244 : /**
245 : * Returns the pointer to the quadrature of specified block and dimension
246 : * @return A pointer.
247 : */
248 5892 : libMesh::QBase * writeableQRule(unsigned int dim, SubdomainID block, InternalDataKey)
249 : {
250 5892 : return qrules(dim, block).vol.get();
251 : }
252 : #endif
253 :
254 : /**
255 : * Returns the reference to the quadrature points
256 : * @return A _reference_. Make sure to store this as a reference!
257 : */
258 195532 : const MooseArray<Point> & qPoints() const { return _current_q_points; }
259 :
260 : /**
261 : * Returns the reference to the mortar segment element quadrature points
262 : * @return A _reference_. Make sure to store this as a reference!
263 : */
264 1224 : const std::vector<Point> & qPointsMortar() const { return _fe_msm->get_xyz(); }
265 :
266 : /**
267 : * The current points in physical space where we have reinited through reinitAtPhysical()
268 : * @return A _reference_. Make sure to store this as a reference!
269 : */
270 954 : const MooseArray<Point> & physicalPoints() const { return _current_physical_points; }
271 :
272 : /**
273 : * Returns the reference to the transformed jacobian weights
274 : * @return A _reference_. Make sure to store this as a reference!
275 : */
276 192474 : const MooseArray<Real> & JxW() const { return _current_JxW; }
277 :
278 4752 : const MooseArray<ADReal> & adJxW() const { return _ad_JxW; }
279 :
280 1704 : const MooseArray<ADReal> & adJxWFace() const { return _ad_JxW_face; }
281 :
282 : const MooseArray<ADReal> & adCurvatures() const;
283 :
284 : /**
285 : * Returns the reference to the coordinate transformation coefficients
286 : * @return A _reference_. Make sure to store this as a reference!
287 : */
288 247628 : const MooseArray<Real> & coordTransformation() const { return _coord; }
289 :
290 : /**
291 : * Returns the reference to the coordinate transformation coefficients on the mortar segment mesh
292 : * @return A _reference_. Make sure to store this as a reference!
293 : */
294 1326 : const MooseArray<Real> & mortarCoordTransformation() const { return _coord_msm; }
295 :
296 : /**
297 : * Returns the reference to the AD version of the coordinate transformation coefficients
298 : * @return A _reference_. Make sure to store this as a reference!
299 : */
300 6456 : const MooseArray<ADReal> & adCoordTransformation() const
301 : {
302 : // Coord values for non-cartesian coordinate systems are functions of the locations of the
303 : // quadrature points in physical space. We also have no way of knowing whether this was called
304 : // from a volumetric or face object so we should set both volumetric and face xyz to true
305 6456 : _calculate_xyz = true;
306 6456 : _calculate_face_xyz = true;
307 :
308 6456 : _calculate_ad_coord = true;
309 6456 : return _ad_coord;
310 : }
311 :
312 : /**
313 : * Get the coordinate system type
314 : * @return A reference to the coordinate system type
315 : */
316 149595 : const Moose::CoordinateSystemType & coordSystem() const { return _coord_type; }
317 :
318 : /**
319 : * Returns the reference to the current quadrature being used on a current face
320 : * @return A _reference_. Make sure to store this as a reference!
321 : */
322 35091326 : const libMesh::QBase * const & qRuleFace() const { return constify_ref(_current_qrule_face); }
323 :
324 : /**
325 : * Returns the reference to the current quadrature being used on a current face
326 : * @return A _reference_. Make sure to store this as a reference!
327 : */
328 132 : libMesh::QBase * const & writeableQRuleFace() { return _current_qrule_face; }
329 :
330 : #ifdef MOOSE_KOKKOS_ENABLED
331 : /**
332 : * Returns the pointer to the quadrature used on a face of specified block and dimension
333 : * @return A pointer.
334 : */
335 5892 : libMesh::QBase * writeableQRuleFace(unsigned int dim, SubdomainID block, InternalDataKey)
336 : {
337 5892 : return qrules(dim, block).face.get();
338 : }
339 : #endif
340 :
341 : /**
342 : * Returns the reference to the current quadrature being used
343 : * @return A _reference_. Make sure to store this as a reference!
344 : */
345 66543 : const MooseArray<Point> & qPointsFace() const { return _current_q_points_face; }
346 :
347 : /**
348 : * Returns the reference to the transformed jacobian weights on a current face
349 : * @return A _reference_. Make sure to store this as a reference!
350 : */
351 55922 : const MooseArray<Real> & JxWFace() const { return _current_JxW_face; }
352 :
353 : /**
354 : * Returns the array of normals for quadrature points on a current side
355 : * @return A _reference_. Make sure to store this as a reference!
356 : */
357 56970 : const MooseArray<Point> & normals() const { return _current_normals; }
358 :
359 : /***
360 : * Returns the array of normals for quadrature points on a current side
361 : */
362 178 : const std::vector<Eigen::Map<RealDIMValue>> & mappedNormals() const { return _mapped_normals; }
363 :
364 : /**
365 : * Returns the array of tangents for quadrature points on a current side
366 : * @return A _reference_. Make sure to store this as a reference!
367 : */
368 1224 : const MooseArray<std::vector<Point>> & tangents() const { return _current_tangents; }
369 :
370 : /**
371 : * Number of extra element integers Assembly tracked
372 : */
373 397342004 : unsigned int numExtraElemIntegers() const { return _extra_elem_ids.size() - 1; }
374 :
375 : /**
376 : * Returns an integer ID of the current element given the index associated with the integer
377 : */
378 508 : const dof_id_type & extraElemID(unsigned int id) const
379 : {
380 : mooseAssert(id < _extra_elem_ids.size(), "An invalid extra element integer id");
381 508 : return _extra_elem_ids[id];
382 : }
383 :
384 : /**
385 : * Returns an integer ID of the current element given the index associated with the integer
386 : */
387 52 : const dof_id_type & extraElemIDNeighbor(unsigned int id) const
388 : {
389 : mooseAssert(id < _neighbor_extra_elem_ids.size(), "An invalid extra element integer id");
390 52 : return _neighbor_extra_elem_ids[id];
391 : }
392 :
393 1583 : const MooseArray<ADPoint> & adNormals() const { return _ad_normals; }
394 :
395 4870 : const MooseArray<ADPoint> & adQPoints() const
396 : {
397 4870 : _calculate_xyz = true;
398 4870 : return _ad_q_points;
399 : }
400 :
401 1583 : const MooseArray<ADPoint> & adQPointsFace() const
402 : {
403 1583 : _calculate_face_xyz = true;
404 1583 : return _ad_q_points_face;
405 : }
406 :
407 : template <bool is_ad>
408 : const MooseArray<Moose::GenericType<Point, is_ad>> & genericQPoints() const;
409 :
410 : /**
411 : * Return the current element
412 : * @return A _reference_. Make sure to store this as a reference!
413 : */
414 392334585 : const Elem * const & elem() const { return _current_elem; }
415 :
416 : /**
417 : * Return the current subdomain ID
418 : */
419 60032 : const SubdomainID & currentSubdomainID() const { return _current_subdomain_id; }
420 :
421 : /**
422 : * set the current subdomain ID
423 : */
424 416661275 : void setCurrentSubdomainID(SubdomainID i) { _current_subdomain_id = i; }
425 :
426 : /**
427 : * Return the current boundary ID
428 : */
429 85125 : const BoundaryID & currentBoundaryID() const { return _current_boundary_id; }
430 :
431 : /**
432 : * set the current boundary ID
433 : */
434 131567919 : void setCurrentBoundaryID(BoundaryID i) { _current_boundary_id = i; }
435 :
436 : /**
437 : * Returns the reference to the current element volume
438 : * @return A _reference_. Make sure to store this as a reference!
439 : */
440 8828497 : const Real & elemVolume() const { return _current_elem_volume; }
441 :
442 : /**
443 : * Returns the current side
444 : * @return A _reference_. Make sure to store this as a reference!
445 : */
446 11067859 : const unsigned int & side() const { return _current_side; }
447 :
448 : /**
449 : * Returns the current neighboring side
450 : * @return A _reference_. Make sure to store this as a reference!
451 : */
452 2134408 : const unsigned int & neighborSide() const { return _current_neighbor_side; }
453 :
454 : /**
455 : * Returns the side element
456 : * @return A _reference_. Make sure to store this as a reference!
457 : */
458 19734 : const Elem * const & sideElem() const { return _current_side_elem; }
459 :
460 : /**
461 : * Returns the reference to the volume of current side element
462 : * @return A _reference_. Make sure to store this as a reference!
463 : */
464 89328 : const Real & sideElemVolume() const { return _current_side_volume; }
465 :
466 : /**
467 : * Return the neighbor element
468 : * @return A _reference_. Make sure to store this as a reference!
469 : */
470 6837285 : const Elem * const & neighbor() const { return _current_neighbor_elem; }
471 :
472 : /**
473 : * Return the lower dimensional element
474 : * @return A _reference_. Make sure to store this as a reference!
475 : */
476 245953 : const Elem * const & lowerDElem() const { return _current_lower_d_elem; }
477 :
478 : /**
479 : * Return the neighboring lower dimensional element
480 : * @return A _reference_. Make sure to store this as a reference!
481 : */
482 1326 : const Elem * const & neighborLowerDElem() const { return _current_neighbor_lower_d_elem; }
483 :
484 : /*
485 : * @return The current lower-dimensional element volume
486 : */
487 : const Real & lowerDElemVolume() const;
488 :
489 : /*
490 : * @return The current neighbor lower-dimensional element volume
491 : */
492 : const Real & neighborLowerDElemVolume() const;
493 :
494 : /**
495 : * Return the current subdomain ID
496 : */
497 10559 : const SubdomainID & currentNeighborSubdomainID() const { return _current_neighbor_subdomain_id; }
498 :
499 : /**
500 : * set the current subdomain ID
501 : */
502 1501013526 : void setCurrentNeighborSubdomainID(SubdomainID i) { _current_neighbor_subdomain_id = i; }
503 :
504 : /**
505 : * Returns the reference to the current neighbor volume
506 : * @return A _reference_. Make sure to store this as a reference!
507 : */
508 3580 : const Real & neighborVolume()
509 : {
510 3580 : _need_neighbor_elem_volume = true;
511 3580 : return _current_neighbor_volume;
512 : }
513 :
514 : /**
515 : * Returns the reference to the current quadrature being used on a current neighbor
516 : * @return A _reference_. Make sure to store this as a reference!
517 : */
518 1311015 : const libMesh::QBase * const & qRuleNeighbor() const
519 : {
520 1311015 : return constify_ref(_current_qrule_neighbor);
521 : }
522 :
523 : /**
524 : * Returns the reference to the current quadrature being used on a current neighbor
525 : * @return A _reference_. Make sure to store this as a reference!
526 : */
527 : libMesh::QBase * const & writeableQRuleNeighbor() { return _current_qrule_neighbor; }
528 :
529 : /**
530 : * Returns the reference to the transformed jacobian weights on a current face
531 : * @return A _reference_. Make sure to store this as a reference!
532 : */
533 : const MooseArray<Real> & JxWNeighbor() const;
534 :
535 : /**
536 : * Returns the reference to the current quadrature points being used on the neighbor face
537 : * @return A _reference_. Make sure to store this as a reference!
538 : */
539 11885 : const MooseArray<Point> & qPointsFaceNeighbor() const { return _current_q_points_face_neighbor; }
540 :
541 : /**
542 : * Returns the reference to the node
543 : * @return A _reference_. Make sure to store this as a reference!
544 : */
545 752788 : const Node * const & node() const { return _current_node; }
546 :
547 : /**
548 : * Returns the reference to the neighboring node
549 : * @return A _reference_. Make sure to store this as a reference!
550 : */
551 174071 : const Node * const & nodeNeighbor() const { return _current_neighbor_node; }
552 :
553 : /**
554 : * Creates block-specific volume, face and arbitrary qrules based on the
555 : * orders and the flag of whether or not to allow negative qweights passed in.
556 : * Any quadrature rules specified using this function override those created
557 : * via in the non-block-specific/global createQRules function. order is used
558 : * for arbitrary volume quadrature rules, while volume_order and face_order
559 : * are for elem and face quadrature respectively.
560 : */
561 : void createQRules(QuadratureType type,
562 : Order order,
563 : Order volume_order,
564 : Order face_order,
565 : SubdomainID block,
566 : bool allow_negative_qweights = true);
567 :
568 : /**
569 : * Increases the element/volume quadrature order for the specified mesh
570 : * block if and only if the current volume quadrature order is lower. This
571 : * works exactly like the bumpAllQRuleOrder function, except it only
572 : * affects the volume quadrature rule (not face quadrature).
573 : */
574 : void bumpVolumeQRuleOrder(Order volume_order, SubdomainID block);
575 :
576 : /**
577 : * Increases the element/volume and face/area quadrature orders for the specified mesh
578 : * block if and only if the current volume or face quadrature order is lower. This
579 : * can only cause the quadrature level to increase. If order is
580 : * lower than or equal to the current volume+face quadrature rule order,
581 : * then nothing is done (i.e. this function is idempotent).
582 : */
583 : void bumpAllQRuleOrder(Order order, SubdomainID block);
584 :
585 : /**
586 : * Set the qrule to be used for volume integration.
587 : *
588 : * Note: This is normally set internally, only use if you know what you are doing!
589 : *
590 : * @param qrule The qrule you want to set
591 : * @param dim The spatial dimension of the qrule
592 : */
593 : void setVolumeQRule(libMesh::QBase * qrule, unsigned int dim);
594 :
595 : /**
596 : * Set the qrule to be used for face integration.
597 : *
598 : * Note: This is normally set internally, only use if you know what you are doing!
599 : *
600 : * @param qrule The qrule you want to set
601 : * @param dim The spatial dimension of the qrule
602 : */
603 : void setFaceQRule(libMesh::QBase * qrule, unsigned int dim);
604 :
605 : /**
606 : * Specifies a custom qrule for integration on mortar segment mesh
607 : *
608 : * Used to properly integrate QUAD face elements using quadrature on TRI mortar segment elements.
609 : * For example, to exactly integrate a FIRST order QUAD element, SECOND order quadrature on TRI
610 : * mortar segments is needed.
611 : */
612 : void setMortarQRule(Order order);
613 :
614 : /**
615 : * Indicates that dual shape functions are used for mortar constraint
616 : */
617 124 : void activateDual() { _need_dual = true; }
618 :
619 : /**
620 : * Indicates whether dual shape functions are used (computation is now repeated on each element
621 : * so expense of computing dual shape functions is no longer trivial)
622 : */
623 312289 : bool needDual() const { return _need_dual; }
624 :
625 : /**
626 : * Set the cached quadrature rules to nullptr
627 : */
628 : void clearCachedQRules();
629 :
630 : private:
631 : /**
632 : * Set the qrule to be used for lower dimensional integration.
633 : *
634 : * @param qrule The qrule you want to set
635 : * @param dim The spatial dimension of the qrule
636 : */
637 : void setLowerQRule(libMesh::QBase * qrule, unsigned int dim);
638 :
639 : public:
640 : /**
641 : * Set the qrule to be used for neighbor integration.
642 : *
643 : * Note: This is normally set internally, only use if you know what you are doing!
644 : *
645 : * @param qrule The qrule you want to set
646 : * @param dim The spatial dimension of the qrule
647 : */
648 : void setNeighborQRule(libMesh::QBase * qrule, unsigned int dim);
649 :
650 : /**
651 : * Reinitialize objects (JxW, q_points, ...) for an elements
652 : *
653 : * @param elem The element we want to reinitialize on
654 : */
655 : void reinit(const Elem * elem);
656 :
657 : /**
658 : * Set the volumetric quadrature rule based on the provided element
659 : */
660 : void setVolumeQRule(const Elem * elem);
661 :
662 : /**
663 : * Reinitialize FE data for the given element on the given side, optionally
664 : * with a given set of reference points
665 : */
666 : void reinitElemFaceRef(const Elem * elem,
667 : unsigned int elem_side,
668 : Real tolerance,
669 : const std::vector<Point> * const pts = nullptr,
670 : const std::vector<Real> * const weights = nullptr);
671 :
672 : /**
673 : * Reinitialize FE data for the given neighbor_element on the given side with a given set of
674 : * reference points
675 : */
676 : void reinitNeighborFaceRef(const Elem * neighbor_elem,
677 : unsigned int neighbor_side,
678 : Real tolerance,
679 : const std::vector<Point> * const pts,
680 : const std::vector<Real> * const weights = nullptr);
681 :
682 : /**
683 : * Reintialize dual basis coefficients based on a customized quadrature rule
684 : */
685 : void reinitDual(const Elem * elem, const std::vector<Point> & pts, const std::vector<Real> & JxW);
686 :
687 : /**
688 : * Reinitialize FE data for a lower dimenesional element with a given set of reference points
689 : */
690 : void reinitLowerDElem(const Elem * elem,
691 : const std::vector<Point> * const pts = nullptr,
692 : const std::vector<Real> * const weights = nullptr);
693 :
694 : /**
695 : * reinitialize a neighboring lower dimensional element
696 : */
697 : void reinitNeighborLowerDElem(const Elem * elem);
698 :
699 : /**
700 : * reinitialize a mortar segment mesh element in order to get a proper JxW
701 : */
702 : void reinitMortarElem(const Elem * elem);
703 :
704 : /**
705 : * Returns a reference to JxW for mortar segment elements
706 : */
707 14321 : const std::vector<Real> & jxWMortar() const { return *_JxW_msm; }
708 :
709 : /**
710 : * Returns a reference to the quadrature rule for the mortar segments
711 : */
712 14321 : const libMesh::QBase * const & qRuleMortar() const { return constify_ref(_qrule_msm); }
713 :
714 : private:
715 : /**
716 : * compute AD things on an element face
717 : */
718 : void computeADFace(const Elem & elem, const unsigned int side);
719 :
720 : public:
721 : /**
722 : * Reinitialize the assembly data at specific physical point in the given element.
723 : */
724 : void reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points);
725 :
726 : /**
727 : * Reinitialize the assembly data at specific points in the reference element.
728 : */
729 : void reinit(const Elem * elem, const std::vector<Point> & reference_points);
730 :
731 : /**
732 : * Set the face quadrature rule based on the provided element and side
733 : */
734 : void setFaceQRule(const Elem * const elem, const unsigned int side);
735 :
736 : /**
737 : * Reinitialize the assembly data on an side of an element
738 : */
739 : void reinit(const Elem * elem, unsigned int side);
740 :
741 : /**
742 : * Reinitialize the assembly data on the side of a element at the custom reference points
743 : */
744 : void reinit(const Elem * elem, unsigned int side, const std::vector<Point> & reference_points);
745 :
746 : void reinitFVFace(const FaceInfo & fi);
747 :
748 : /**
749 : * Reinitialize an element and its neighbor along a particular side.
750 : *
751 : * @param elem Element being reinitialized
752 : * @param side Side of the element
753 : * @param neighbor Neighbor facing the element on the side 'side'
754 : * @param neighbor_side The side id on the neighboring element.
755 : * @param neighbor_reference_points Optional argument specifying the neighbor reference points. If
756 : * not passed, then neighbor reference points will be determined by doing an inverse map based on
757 : * the physical location of the \p elem quadrature points
758 : */
759 : void reinitElemAndNeighbor(const Elem * elem,
760 : unsigned int side,
761 : const Elem * neighbor,
762 : unsigned int neighbor_side,
763 : const std::vector<Point> * neighbor_reference_points = nullptr);
764 :
765 : /**
766 : * Reinitializes the neighbor at the physical coordinates on neighbor side given.
767 : */
768 : void reinitNeighborAtPhysical(const Elem * neighbor,
769 : unsigned int neighbor_side,
770 : const std::vector<Point> & physical_points);
771 :
772 : /**
773 : * Reinitializes the neighbor at the physical coordinates within element given.
774 : */
775 : void reinitNeighborAtPhysical(const Elem * neighbor, const std::vector<Point> & physical_points);
776 :
777 : /**
778 : * Reinitializes the neighbor side using reference coordinates.
779 : */
780 : void reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);
781 :
782 : /**
783 : * Reinitialize assembly data for a node
784 : */
785 : void reinit(const Node * node);
786 :
787 : /**
788 : * Initialize the Assembly object and set the CouplingMatrix for use throughout.
789 : */
790 : void init(const libMesh::CouplingMatrix * cm);
791 :
792 : /// Create pair of variables requiring nonlocal jacobian contributions
793 : void initNonlocalCoupling();
794 :
795 : /// Sizes and zeroes the Jacobian blocks used for the current element
796 : void prepareJacobianBlock();
797 :
798 : /// Sizes and zeroes the residual for the current element
799 : void prepareResidual();
800 :
801 : void prepare();
802 : void prepareNonlocal();
803 :
804 : /**
805 : * Used for preparing the dense residual and jacobian blocks for one particular variable.
806 : *
807 : * @param var The variable that needs to have its datastructures prepared
808 : */
809 : void prepareVariable(MooseVariableFieldBase * var);
810 : void prepareVariableNonlocal(MooseVariableFieldBase * var);
811 : void prepareNeighbor();
812 :
813 : /**
814 : * Prepare the Jacobians and residuals for a lower dimensional element. This method may be called
815 : * when performing mortar finite element simulations
816 : */
817 : void prepareLowerD();
818 :
819 : void prepareBlock(unsigned int ivar, unsigned jvar, const std::vector<dof_id_type> & dof_indices);
820 : void prepareBlockNonlocal(unsigned int ivar,
821 : unsigned jvar,
822 : const std::vector<dof_id_type> & idof_indices,
823 : const std::vector<dof_id_type> & jdof_indices);
824 : void prepareScalar();
825 : void prepareOffDiagScalar();
826 :
827 : template <typename T>
828 : void copyShapes(MooseVariableField<T> & v);
829 : void copyShapes(unsigned int var);
830 :
831 : template <typename T>
832 : void copyFaceShapes(MooseVariableField<T> & v);
833 : void copyFaceShapes(unsigned int var);
834 :
835 : template <typename T>
836 : void copyNeighborShapes(MooseVariableField<T> & v);
837 : void copyNeighborShapes(unsigned int var);
838 :
839 : /**
840 : * Key structure for APIs manipulating global vectors/matrices. Developers in blessed classes may
841 : * create keys using simple curly braces \p {} or may be more explicit and use \p
842 : * Assembly::GlobalDataKey{}
843 : */
844 : class GlobalDataKey
845 : {
846 : // Blessed classes
847 : friend class Assembly;
848 : friend class SubProblem;
849 : friend class FEProblemBase;
850 : friend class DisplacedProblem;
851 : friend class ComputeMortarFunctor;
852 : friend class NonlinearSystemBase;
853 554806119 : GlobalDataKey() {}
854 : GlobalDataKey(const GlobalDataKey &) {}
855 : };
856 :
857 : /**
858 : * Key structure for APIs adding/caching local element residuals/Jacobians. Developers in blessed
859 : * classes may create keys using simple curly braces \p {} or may be more explicit and use \p
860 : * Assembly::LocalDataKey{}
861 : */
862 : class LocalDataKey
863 : {
864 : // Blessed classes
865 : friend class Assembly;
866 : friend class TaggingInterface;
867 1638324738 : LocalDataKey() {}
868 : LocalDataKey(const LocalDataKey &) {}
869 : };
870 :
871 : /**
872 : * Add local residuals of all field variables for a set of tags onto the global residual vectors
873 : * associated with the tags.
874 : */
875 : void addResidual(GlobalDataKey, const std::vector<VectorTag> & vector_tags);
876 : /**
877 : * Add local neighbor residuals of all field variables for a set of tags onto the global residual
878 : * vectors associated with the tags.
879 : */
880 : void addResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & vector_tags);
881 : /**
882 : * Add local neighbor residuals of all field variables for a set of tags onto the global residual
883 : * vectors associated with the tags.
884 : */
885 : void addResidualLower(GlobalDataKey, const std::vector<VectorTag> & vector_tags);
886 :
887 : /**
888 : * Add residuals of all scalar variables for a set of tags onto the global residual vectors
889 : * associated with the tags.
890 : */
891 : void addResidualScalar(GlobalDataKey, const std::vector<VectorTag> & vector_tags);
892 :
893 : /**
894 : * Takes the values that are currently in _sub_Re of all field variables and appends them to
895 : * the cached values.
896 : */
897 : void cacheResidual(GlobalDataKey, const std::vector<VectorTag> & tags);
898 :
899 : /**
900 : * Takes the values that are currently in _sub_Rn of all field variables and appends them to
901 : * the cached values.
902 : */
903 : void cacheResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & tags);
904 :
905 : /**
906 : * Takes the values that are currently in _sub_Rl and appends them to the cached values.
907 : */
908 : void cacheResidualLower(GlobalDataKey, const std::vector<VectorTag> & tags);
909 :
910 : /**
911 : * Pushes all cached residuals to the global residual vectors associated with each tag.
912 : *
913 : * Note that this will also clear the cache.
914 : */
915 : void addCachedResiduals(GlobalDataKey, const std::vector<VectorTag> & tags);
916 :
917 : /**
918 : * Clears all of the residuals in _cached_residual_rows and _cached_residual_values
919 : *
920 : * This method is designed specifically for use after calling
921 : * FEProblemBase::addCachedResidualDirectly() and DisplacedProblem::addCachedResidualDirectly() to
922 : * ensure that we don't have any extra residuals hanging around that we didn't have the vectors
923 : * for
924 : */
925 : void clearCachedResiduals(GlobalDataKey);
926 :
927 : /**
928 : * Adds the values that have been cached by calling cacheResidual(), cacheResidualNeighbor(),
929 : * and/or cacheResidualLower() to a user-defined residual (that is, not necessarily the vector
930 : * that vector_tag points to)
931 : *
932 : * Note that this will also clear the cache.
933 : */
934 : void addCachedResidualDirectly(NumericVector<Number> & residual,
935 : GlobalDataKey,
936 : const VectorTag & vector_tag);
937 :
938 : /**
939 : * Sets local residuals of all field variables to the global residual vector for a tag.
940 : */
941 : void setResidual(NumericVector<Number> & residual, GlobalDataKey, const VectorTag & vector_tag);
942 :
943 : /**
944 : * Sets local neighbor residuals of all field variables to the global residual vector for a tag.
945 : */
946 : void setResidualNeighbor(NumericVector<Number> & residual,
947 : GlobalDataKey,
948 : const VectorTag & vector_tag);
949 :
950 : /**
951 : * Adds all local Jacobian to the global Jacobian matrices.
952 : */
953 : void addJacobian(GlobalDataKey);
954 :
955 : /**
956 : * Adds non-local Jacobian to the global Jacobian matrices.
957 : */
958 : void addJacobianNonlocal(GlobalDataKey);
959 :
960 : /**
961 : * Add ElementNeighbor, NeighborElement, and NeighborNeighbor portions of the Jacobian for compute
962 : * objects like DGKernels
963 : */
964 : void addJacobianNeighbor(GlobalDataKey);
965 :
966 : /**
967 : * Add Jacobians for pairs of scalar variables into the global Jacobian matrices.
968 : */
969 : void addJacobianScalar(GlobalDataKey);
970 :
971 : /**
972 : * Add Jacobians for a scalar variables with all other field variables into the global Jacobian
973 : * matrices.
974 : */
975 : void addJacobianOffDiagScalar(unsigned int ivar, GlobalDataKey);
976 :
977 : /**
978 : * Adds element matrix for ivar rows and jvar columns to the global Jacobian matrix.
979 : */
980 : void addJacobianBlock(libMesh::SparseMatrix<Number> & jacobian,
981 : unsigned int ivar,
982 : unsigned int jvar,
983 : const libMesh::DofMap & dof_map,
984 : std::vector<dof_id_type> & dof_indices,
985 : GlobalDataKey,
986 : TagID tag);
987 :
988 : /**
989 : * Add element matrix for ivar rows and jvar columns to the global Jacobian matrix for given
990 : * tags.
991 : */
992 : void addJacobianBlockTags(libMesh::SparseMatrix<Number> & jacobian,
993 : unsigned int ivar,
994 : unsigned int jvar,
995 : const libMesh::DofMap & dof_map,
996 : std::vector<dof_id_type> & dof_indices,
997 : GlobalDataKey,
998 : const std::set<TagID> & tags);
999 :
1000 : /**
1001 : * Adds non-local element matrix for ivar rows and jvar columns to the global Jacobian matrix.
1002 : */
1003 : void addJacobianBlockNonlocal(libMesh::SparseMatrix<Number> & jacobian,
1004 : unsigned int ivar,
1005 : unsigned int jvar,
1006 : const libMesh::DofMap & dof_map,
1007 : const std::vector<dof_id_type> & idof_indices,
1008 : const std::vector<dof_id_type> & jdof_indices,
1009 : GlobalDataKey,
1010 : TagID tag);
1011 :
1012 : /**
1013 : * Adds non-local element matrix for ivar rows and jvar columns to the global Jacobian matrix.
1014 : */
1015 : void addJacobianBlockNonlocalTags(libMesh::SparseMatrix<Number> & jacobian,
1016 : unsigned int ivar,
1017 : unsigned int jvar,
1018 : const libMesh::DofMap & dof_map,
1019 : const std::vector<dof_id_type> & idof_indices,
1020 : const std::vector<dof_id_type> & jdof_indices,
1021 : GlobalDataKey,
1022 : const std::set<TagID> & tags);
1023 :
1024 : /**
1025 : * Add *all* portions of the Jacobian except PrimaryPrimary, e.g. LowerLower, LowerSecondary,
1026 : * LowerPrimary, SecondaryLower, SecondarySecondary, SecondaryPrimary, PrimaryLower,
1027 : * PrimarySecondary, for mortar-like objects. Primary indicates the interior parent element on the
1028 : * primary side of the mortar interface. Secondary indicates the neighbor of the interior parent
1029 : * element. Lower denotes the lower-dimensional element living on the primary side of the mortar
1030 : * interface.
1031 : */
1032 : void addJacobianNeighborLowerD(GlobalDataKey);
1033 :
1034 : /**
1035 : * Add portions of the Jacobian of LowerLower, LowerSecondary, and SecondaryLower for
1036 : * boundary conditions. Secondary indicates the boundary element. Lower denotes the
1037 : * lower-dimensional element living on the boundary side.
1038 : */
1039 : void addJacobianLowerD(GlobalDataKey);
1040 :
1041 : /**
1042 : * Cache *all* portions of the Jacobian, e.g. LowerLower, LowerSecondary, LowerPrimary,
1043 : * SecondaryLower, SecondarySecondary, SecondaryPrimary, PrimaryLower, PrimarySecondary,
1044 : * PrimaryPrimary for mortar-like objects. Primary indicates the interior parent element on the
1045 : * primary side of the mortar interface. Secondary indicates the interior parent element on the
1046 : * secondary side of the interface. Lower denotes the lower-dimensional element living on the
1047 : * secondary side of the mortar interface; it's the boundary face of the \p Secondary element.
1048 : */
1049 : void cacheJacobianMortar(GlobalDataKey);
1050 :
1051 : /**
1052 : * Adds three neighboring element matrices for ivar rows and jvar columns to the global Jacobian
1053 : * matrix.
1054 : */
1055 : void addJacobianNeighbor(libMesh::SparseMatrix<Number> & jacobian,
1056 : unsigned int ivar,
1057 : unsigned int jvar,
1058 : const libMesh::DofMap & dof_map,
1059 : std::vector<dof_id_type> & dof_indices,
1060 : std::vector<dof_id_type> & neighbor_dof_indices,
1061 : GlobalDataKey,
1062 : TagID tag);
1063 :
1064 : /**
1065 : * Adds three neighboring element matrices for ivar rows and jvar columns to the global Jacobian
1066 : * matrix.
1067 : */
1068 : void addJacobianNeighborTags(libMesh::SparseMatrix<Number> & jacobian,
1069 : unsigned int ivar,
1070 : unsigned int jvar,
1071 : const libMesh::DofMap & dof_map,
1072 : std::vector<dof_id_type> & dof_indices,
1073 : std::vector<dof_id_type> & neighbor_dof_indices,
1074 : GlobalDataKey,
1075 : const std::set<TagID> & tags);
1076 :
1077 : /**
1078 : * Takes the values that are currently in _sub_Kee and appends them to the cached values.
1079 : */
1080 : void cacheJacobian(GlobalDataKey);
1081 :
1082 : /**
1083 : * Takes the values that are currently in _sub_Keg and appends them to the cached values.
1084 : */
1085 : void cacheJacobianNonlocal(GlobalDataKey);
1086 :
1087 : /**
1088 : * Takes the values that are currently in the neighbor Dense Matrices and appends them to the
1089 : * cached values.
1090 : */
1091 : void cacheJacobianNeighbor(GlobalDataKey);
1092 :
1093 : /**
1094 : * Adds the values that have been cached by calling cacheJacobian() and or cacheJacobianNeighbor()
1095 : * to the jacobian matrix.
1096 : *
1097 : * Note that this will also clear the cache.
1098 : */
1099 : void addCachedJacobian(GlobalDataKey);
1100 :
1101 : /**
1102 : * Sets previously-cached Jacobian values via SparseMatrix::set() calls.
1103 : */
1104 : void setCachedJacobian(GlobalDataKey);
1105 :
1106 : /**
1107 : * Zero out previously-cached Jacobian rows.
1108 : */
1109 : void zeroCachedJacobian(GlobalDataKey);
1110 :
1111 : /**
1112 : * Get local residual block for a variable and a tag. Only blessed framework classes may call this
1113 : * API by creating the requisiste \p LocalDataKey class
1114 : */
1115 667794784 : DenseVector<Number> & residualBlock(unsigned int var_num, LocalDataKey, TagID tag_id)
1116 : {
1117 667794784 : return _sub_Re[tag_id][var_num];
1118 : }
1119 :
1120 : /**
1121 : * Get local neighbor residual block for a variable and a tag. Only blessed framework classes may
1122 : * call this API by creating the requisiste \p LocalDataKey class
1123 : */
1124 15448450 : DenseVector<Number> & residualBlockNeighbor(unsigned int var_num, LocalDataKey, TagID tag_id)
1125 : {
1126 15448450 : return _sub_Rn[tag_id][var_num];
1127 : }
1128 :
1129 : /**
1130 : * Get residual block for lower. Only blessed framework classes may call this API by creating the
1131 : * requisiste \p LocalDataKey class
1132 : */
1133 100546 : DenseVector<Number> & residualBlockLower(unsigned int var_num, LocalDataKey, TagID tag_id)
1134 : {
1135 100546 : return _sub_Rl[tag_id][var_num];
1136 : }
1137 :
1138 : /**
1139 : * Get local Jacobian block for a pair of variables and a tag. Only blessed framework classes may
1140 : * call this API by creating the requisiste \p LocalDataKey class
1141 : */
1142 438324716 : DenseMatrix<Number> & jacobianBlock(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
1143 : {
1144 438324716 : jacobianBlockUsed(tag, ivar, jvar, true);
1145 438324716 : return _sub_Kee[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
1146 : }
1147 :
1148 : /**
1149 : * Get local Jacobian block from non-local contribution for a pair of variables and a tag. Only
1150 : * blessed framework classes may call this API by creating the requisiste \p LocalDataKey class
1151 : */
1152 : DenseMatrix<Number> &
1153 19344 : jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
1154 : {
1155 19344 : jacobianBlockNonlocalUsed(tag, ivar, jvar, true);
1156 19344 : return _sub_Keg[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
1157 : }
1158 :
1159 : /**
1160 : * Get local Jacobian block of a DG Jacobian type for a pair of variables and a tag. Only blessed
1161 : * framework classes may call this API by creating the requisiste \p LocalDataKey class
1162 : */
1163 : DenseMatrix<Number> & jacobianBlockNeighbor(
1164 : Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag);
1165 :
1166 : /**
1167 : * Returns the jacobian block for the given mortar Jacobian type. This jacobian block can involve
1168 : * degrees of freedom from the secondary side interior parent, the primary side
1169 : * interior parent, or the lower-dimensional element (located on the secondary
1170 : * side). Only blessed framework classes may call this API by creating the requisiste \p
1171 : * LocalDataKey class
1172 : */
1173 : DenseMatrix<Number> & jacobianBlockMortar(Moose::ConstraintJacobianType type,
1174 : unsigned int ivar,
1175 : unsigned int jvar,
1176 : LocalDataKey,
1177 : TagID tag);
1178 :
1179 : /**
1180 : * Lets an external class cache residual at a set of nodes. Only blessed framework classes may
1181 : * call this API by creating the requisiste \p LocalDataKey class
1182 : */
1183 : void cacheResidualNodes(const DenseVector<Number> & res,
1184 : const std::vector<dof_id_type> & dof_index,
1185 : LocalDataKey,
1186 : TagID tag);
1187 :
1188 : /**
1189 : * Caches the Jacobian entry 'value', to eventually be
1190 : * added/set in the (i,j) location of the matrix.
1191 : *
1192 : * We use numeric_index_type for the index arrays (rather than
1193 : * dof_id_type) since that is what the SparseMatrix interface uses,
1194 : * but at the time of this writing, those two types are equivalent.
1195 : *
1196 : * Only blessed framework classes may call this API by creating the requisiste \p LocalDataKey
1197 : * class
1198 : */
1199 : void
1200 : cacheJacobian(numeric_index_type i, numeric_index_type j, Real value, LocalDataKey, TagID tag);
1201 :
1202 : /**
1203 : * Caches the Jacobian entry 'value', to eventually be
1204 : * added/set in the (i,j) location of the matrices in corresponding to \p tags.
1205 : *
1206 : * We use numeric_index_type for the index arrays (rather than
1207 : * dof_id_type) since that is what the SparseMatrix interface uses,
1208 : * but at the time of this writing, those two types are equivalent.
1209 : *
1210 : * Only blessed framework classes may call this API by creating the requisiste \p LocalDataKey
1211 : * class
1212 : */
1213 : void cacheJacobian(numeric_index_type i,
1214 : numeric_index_type j,
1215 : Real value,
1216 : LocalDataKey,
1217 : const std::set<TagID> & tags);
1218 :
1219 : /**
1220 : * Cache a local Jacobian block with the provided rows (\p idof_indices) and columns (\p
1221 : * jdof_indices) for eventual accumulation into the global matrix specified by \p tag. The \p
1222 : * scaling_factor will be applied before caching. Only blessed framework classes may call this API
1223 : * by creating the requisiste \p LocalDataKey class
1224 : */
1225 : void cacheJacobianBlock(DenseMatrix<Number> & jac_block,
1226 : const std::vector<dof_id_type> & idof_indices,
1227 : const std::vector<dof_id_type> & jdof_indices,
1228 : Real scaling_factor,
1229 : LocalDataKey,
1230 : TagID tag);
1231 :
1232 : /**
1233 : * Process the supplied residual values. This is a mirror of of the non-templated version of \p
1234 : * addResiduals except that it's meant for \emph only processing residuals (and not their
1235 : * derivatives/Jacobian). We supply this API such that residual objects that leverage the AD
1236 : * version of this method when computing the Jacobian (or residual + Jacobian) can mirror the same
1237 : * behavior when doing pure residual evaluations, such as when evaluating linear residuals during
1238 : * (P)JFNK. This method will call \p constrain_element_vector on the supplied residuals. Only
1239 : * blessed framework classes may call this API by creating the requisiste \p LocalDataKey class
1240 : */
1241 : template <typename Residuals, typename Indices>
1242 : void cacheResiduals(const Residuals & residuals,
1243 : const Indices & row_indices,
1244 : Real scaling_factor,
1245 : LocalDataKey,
1246 : const std::set<TagID> & vector_tags);
1247 :
1248 : /**
1249 : * Process the \p derivatives() data of a vector of \p ADReals. This
1250 : * method simply caches the derivative values for the corresponding column indices for the
1251 : * provided \p matrix_tags. Note that this overload will call \p DofMap::constrain_element_matrix.
1252 : * Only blessed framework classes may call this API by creating the requisiste \p LocalDataKey
1253 : * class
1254 : */
1255 : template <typename Residuals, typename Indices>
1256 : void cacheJacobian(const Residuals & residuals,
1257 : const Indices & row_indices,
1258 : Real scaling_factor,
1259 : LocalDataKey,
1260 : const std::set<TagID> & matrix_tags);
1261 :
1262 : /**
1263 : * Process the supplied residual values. This is a mirror of of the non-templated version of \p
1264 : * addResiduals except that it's meant for \emph only processing residuals (and not their
1265 : * derivatives/Jacobian). We supply this API such that residual objects that leverage the AD
1266 : * version of this method when computing the Jacobian (or residual + Jacobian) can mirror the same
1267 : * behavior when doing pure residual evaluations, such as when evaluating linear residuals during
1268 : * (P)JFNK. This method will \emph not call \p constrain_element_vector on the supplied residuals.
1269 : * Only blessed framework classes may call this API by creating the requisiste \p LocalDataKey
1270 : * class
1271 : */
1272 : template <typename Residuals, typename Indices>
1273 : void cacheResidualsWithoutConstraints(const Residuals & residuals,
1274 : const Indices & row_indices,
1275 : Real scaling_factor,
1276 : LocalDataKey,
1277 : const std::set<TagID> & vector_tags);
1278 :
1279 : /**
1280 : * Process the \p derivatives() data of a vector of \p ADReals. This
1281 : * method simply caches the derivative values for the corresponding column indices for the
1282 : * provided \p matrix_tags. Note that this overload will \emph not call \p
1283 : * DofMap::constrain_element_matrix. Only blessed framework classes may call this API by creating
1284 : * the requisiste \p LocalDataKey class
1285 : */
1286 : template <typename Residuals, typename Indices>
1287 : void cacheJacobianWithoutConstraints(const Residuals & residuals,
1288 : const Indices & row_indices,
1289 : Real scaling_factor,
1290 : LocalDataKey,
1291 : const std::set<TagID> & matrix_tags);
1292 :
1293 15184614 : std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> & couplingEntries()
1294 : {
1295 15184614 : return _cm_ff_entry;
1296 : }
1297 : const std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> &
1298 37577 : couplingEntries() const
1299 : {
1300 37577 : return _cm_ff_entry;
1301 : }
1302 : std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> &
1303 4162 : nonlocalCouplingEntries()
1304 : {
1305 4162 : return _cm_nonlocal_entry;
1306 : }
1307 : const std::vector<std::pair<MooseVariableFieldBase *, MooseVariableScalar *>> &
1308 : fieldScalarCouplingEntries() const
1309 : {
1310 : return _cm_fs_entry;
1311 : }
1312 : const std::vector<std::pair<MooseVariableScalar *, MooseVariableFieldBase *>> &
1313 3092 : scalarFieldCouplingEntries() const
1314 : {
1315 3092 : return _cm_sf_entry;
1316 : }
1317 :
1318 : // Read-only references
1319 44 : const VariablePhiValue & phi() const { return _phi; }
1320 : template <typename T>
1321 4539 : const ADTemplateVariablePhiGradient<T> & adGradPhi(const MooseVariableFE<T> & v) const
1322 : {
1323 4539 : return _ad_grad_phi_data.at(v.feType());
1324 : }
1325 : const VariablePhiValue & phi(const MooseVariableField<Real> &) const { return _phi; }
1326 44 : const VariablePhiGradient & gradPhi() const { return _grad_phi; }
1327 : const VariablePhiGradient & gradPhi(const MooseVariableField<Real> &) const { return _grad_phi; }
1328 : const VariablePhiSecond & secondPhi() const { return _second_phi; }
1329 : const VariablePhiSecond & secondPhi(const MooseVariableField<Real> &) const
1330 : {
1331 : return _second_phi;
1332 : }
1333 :
1334 56 : const VariablePhiValue & phiFace() const { return _phi_face; }
1335 : const VariablePhiValue & phiFace(const MooseVariableField<Real> &) const { return _phi_face; }
1336 56 : const VariablePhiGradient & gradPhiFace() const { return _grad_phi_face; }
1337 : const VariablePhiGradient & gradPhiFace(const MooseVariableField<Real> &) const
1338 : {
1339 : return _grad_phi_face;
1340 : }
1341 : const VariablePhiSecond & secondPhiFace(const MooseVariableField<Real> &) const
1342 : {
1343 : return _second_phi_face;
1344 : }
1345 :
1346 : const VariablePhiValue & phiNeighbor(const MooseVariableField<Real> &) const
1347 : {
1348 : return _phi_neighbor;
1349 : }
1350 : const VariablePhiGradient & gradPhiNeighbor(const MooseVariableField<Real> &) const
1351 : {
1352 : return _grad_phi_neighbor;
1353 : }
1354 : const VariablePhiSecond & secondPhiNeighbor(const MooseVariableField<Real> &) const
1355 : {
1356 : return _second_phi_neighbor;
1357 : }
1358 :
1359 : const VariablePhiValue & phiFaceNeighbor(const MooseVariableField<Real> &) const
1360 : {
1361 : return _phi_face_neighbor;
1362 : }
1363 : const VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<Real> &) const
1364 : {
1365 : return _grad_phi_face_neighbor;
1366 : }
1367 : const VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<Real> &) const
1368 : {
1369 : return _second_phi_face_neighbor;
1370 : }
1371 :
1372 : const VectorVariablePhiValue & phi(const MooseVariableField<RealVectorValue> &) const
1373 : {
1374 : return _vector_phi;
1375 : }
1376 : const VectorVariablePhiGradient & gradPhi(const MooseVariableField<RealVectorValue> &) const
1377 : {
1378 : return _vector_grad_phi;
1379 : }
1380 : const VectorVariablePhiSecond & secondPhi(const MooseVariableField<RealVectorValue> &) const
1381 : {
1382 : return _vector_second_phi;
1383 : }
1384 : const VectorVariablePhiCurl & curlPhi(const MooseVariableField<RealVectorValue> &) const
1385 : {
1386 : return _vector_curl_phi;
1387 : }
1388 : const VectorVariablePhiDivergence & divPhi(const MooseVariableField<RealVectorValue> &) const
1389 : {
1390 : return _vector_div_phi;
1391 : }
1392 :
1393 : const VectorVariablePhiValue & phiFace(const MooseVariableField<RealVectorValue> &) const
1394 : {
1395 : return _vector_phi_face;
1396 : }
1397 : const VectorVariablePhiGradient & gradPhiFace(const MooseVariableField<RealVectorValue> &) const
1398 : {
1399 : return _vector_grad_phi_face;
1400 : }
1401 : const VectorVariablePhiSecond & secondPhiFace(const MooseVariableField<RealVectorValue> &) const
1402 : {
1403 : return _vector_second_phi_face;
1404 : }
1405 : const VectorVariablePhiCurl & curlPhiFace(const MooseVariableField<RealVectorValue> &) const
1406 : {
1407 : return _vector_curl_phi_face;
1408 : }
1409 : const VectorVariablePhiDivergence & divPhiFace(const MooseVariableField<RealVectorValue> &) const
1410 : {
1411 : return _vector_div_phi_face;
1412 : }
1413 :
1414 : const VectorVariablePhiValue & phiNeighbor(const MooseVariableField<RealVectorValue> &) const
1415 : {
1416 : return _vector_phi_neighbor;
1417 : }
1418 : const VectorVariablePhiGradient &
1419 : gradPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
1420 : {
1421 : return _vector_grad_phi_neighbor;
1422 : }
1423 : const VectorVariablePhiSecond &
1424 : secondPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
1425 : {
1426 : return _vector_second_phi_neighbor;
1427 : }
1428 : const VectorVariablePhiCurl & curlPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
1429 : {
1430 : return _vector_curl_phi_neighbor;
1431 : }
1432 : const VectorVariablePhiDivergence &
1433 : divPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
1434 : {
1435 : return _vector_div_phi_neighbor;
1436 : }
1437 :
1438 : const VectorVariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
1439 : {
1440 : return _vector_phi_face_neighbor;
1441 : }
1442 : const VectorVariablePhiGradient &
1443 : gradPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
1444 : {
1445 : return _vector_grad_phi_face_neighbor;
1446 : }
1447 : const VectorVariablePhiSecond &
1448 : secondPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
1449 : {
1450 : return _vector_second_phi_face_neighbor;
1451 : }
1452 : const VectorVariablePhiCurl &
1453 : curlPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
1454 : {
1455 : return _vector_curl_phi_face_neighbor;
1456 : }
1457 : const VectorVariablePhiDivergence &
1458 : divPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
1459 : {
1460 : return _vector_div_phi_face_neighbor;
1461 : }
1462 :
1463 : // Writeable references
1464 123480166 : VariablePhiValue & phi(const MooseVariableField<Real> &) { return _phi; }
1465 123478483 : VariablePhiGradient & gradPhi(const MooseVariableField<Real> &) { return _grad_phi; }
1466 24309 : VariablePhiSecond & secondPhi(const MooseVariableField<Real> &) { return _second_phi; }
1467 :
1468 520769 : VariablePhiValue & phiFace(const MooseVariableField<Real> &) { return _phi_face; }
1469 520671 : VariablePhiGradient & gradPhiFace(const MooseVariableField<Real> &) { return _grad_phi_face; }
1470 4248 : VariablePhiSecond & secondPhiFace(const MooseVariableField<Real> &) { return _second_phi_face; }
1471 :
1472 185122 : VariablePhiValue & phiNeighbor(const MooseVariableField<Real> &) { return _phi_neighbor; }
1473 185039 : VariablePhiGradient & gradPhiNeighbor(const MooseVariableField<Real> &)
1474 : {
1475 185039 : return _grad_phi_neighbor;
1476 : }
1477 0 : VariablePhiSecond & secondPhiNeighbor(const MooseVariableField<Real> &)
1478 : {
1479 0 : return _second_phi_neighbor;
1480 : }
1481 :
1482 187115 : VariablePhiValue & phiFaceNeighbor(const MooseVariableField<Real> &)
1483 : {
1484 187115 : return _phi_face_neighbor;
1485 : }
1486 186981 : VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<Real> &)
1487 : {
1488 186981 : return _grad_phi_face_neighbor;
1489 : }
1490 0 : VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<Real> &)
1491 : {
1492 0 : return _second_phi_face_neighbor;
1493 : }
1494 :
1495 : // Writeable references with vector variable
1496 1818404 : VectorVariablePhiValue & phi(const MooseVariableField<RealVectorValue> &) { return _vector_phi; }
1497 1818325 : VectorVariablePhiGradient & gradPhi(const MooseVariableField<RealVectorValue> &)
1498 : {
1499 1818325 : return _vector_grad_phi;
1500 : }
1501 0 : VectorVariablePhiSecond & secondPhi(const MooseVariableField<RealVectorValue> &)
1502 : {
1503 0 : return _vector_second_phi;
1504 : }
1505 50650 : VectorVariablePhiCurl & curlPhi(const MooseVariableField<RealVectorValue> &)
1506 : {
1507 50650 : return _vector_curl_phi;
1508 : }
1509 494442 : VectorVariablePhiDivergence & divPhi(const MooseVariableField<RealVectorValue> &)
1510 : {
1511 494442 : return _vector_div_phi;
1512 : }
1513 :
1514 56778 : VectorVariablePhiValue & phiFace(const MooseVariableField<RealVectorValue> &)
1515 : {
1516 56778 : return _vector_phi_face;
1517 : }
1518 56394 : VectorVariablePhiGradient & gradPhiFace(const MooseVariableField<RealVectorValue> &)
1519 : {
1520 56394 : return _vector_grad_phi_face;
1521 : }
1522 0 : VectorVariablePhiSecond & secondPhiFace(const MooseVariableField<RealVectorValue> &)
1523 : {
1524 0 : return _vector_second_phi_face;
1525 : }
1526 : VectorVariablePhiCurl & curlPhiFace(const MooseVariableField<RealVectorValue> &)
1527 : {
1528 : return _vector_curl_phi_face;
1529 : }
1530 : VectorVariablePhiDivergence & divPhiFace(const MooseVariableField<RealVectorValue> &)
1531 : {
1532 : return _vector_div_phi_face;
1533 : }
1534 :
1535 266 : VectorVariablePhiValue & phiNeighbor(const MooseVariableField<RealVectorValue> &)
1536 : {
1537 266 : return _vector_phi_neighbor;
1538 : }
1539 266 : VectorVariablePhiGradient & gradPhiNeighbor(const MooseVariableField<RealVectorValue> &)
1540 : {
1541 266 : return _vector_grad_phi_neighbor;
1542 : }
1543 0 : VectorVariablePhiSecond & secondPhiNeighbor(const MooseVariableField<RealVectorValue> &)
1544 : {
1545 0 : return _vector_second_phi_neighbor;
1546 : }
1547 : VectorVariablePhiCurl & curlPhiNeighbor(const MooseVariableField<RealVectorValue> &)
1548 : {
1549 : return _vector_curl_phi_neighbor;
1550 : }
1551 : VectorVariablePhiDivergence & divPhiNeighbor(const MooseVariableField<RealVectorValue> &)
1552 : {
1553 : return _vector_div_phi_neighbor;
1554 : }
1555 311 : VectorVariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
1556 : {
1557 311 : return _vector_phi_face_neighbor;
1558 : }
1559 288 : VectorVariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
1560 : {
1561 288 : return _vector_grad_phi_face_neighbor;
1562 : }
1563 0 : VectorVariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
1564 : {
1565 0 : return _vector_second_phi_face_neighbor;
1566 : }
1567 : VectorVariablePhiCurl & curlPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
1568 : {
1569 : return _vector_curl_phi_face_neighbor;
1570 : }
1571 : VectorVariablePhiDivergence & divPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
1572 : {
1573 : return _vector_div_phi_face_neighbor;
1574 : }
1575 :
1576 : // Writeable references with array variable
1577 1177340 : VariablePhiValue & phi(const MooseVariableField<RealEigenVector> &) { return _phi; }
1578 1177340 : VariablePhiGradient & gradPhi(const MooseVariableField<RealEigenVector> &) { return _grad_phi; }
1579 0 : VariablePhiSecond & secondPhi(const MooseVariableField<RealEigenVector> &) { return _second_phi; }
1580 :
1581 13084 : VariablePhiValue & phiFace(const MooseVariableField<RealEigenVector> &) { return _phi_face; }
1582 12810 : VariablePhiGradient & gradPhiFace(const MooseVariableField<RealEigenVector> &)
1583 : {
1584 12810 : return _grad_phi_face;
1585 : }
1586 0 : VariablePhiSecond & secondPhiFace(const MooseVariableField<RealEigenVector> &)
1587 : {
1588 0 : return _second_phi_face;
1589 : }
1590 :
1591 3456 : VariablePhiValue & phiNeighbor(const MooseVariableField<RealEigenVector> &)
1592 : {
1593 3456 : return _phi_neighbor;
1594 : }
1595 3456 : VariablePhiGradient & gradPhiNeighbor(const MooseVariableField<RealEigenVector> &)
1596 : {
1597 3456 : return _grad_phi_neighbor;
1598 : }
1599 0 : VariablePhiSecond & secondPhiNeighbor(const MooseVariableField<RealEigenVector> &)
1600 : {
1601 0 : return _second_phi_neighbor;
1602 : }
1603 :
1604 3456 : VariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealEigenVector> &)
1605 : {
1606 3456 : return _phi_face_neighbor;
1607 : }
1608 3456 : VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<RealEigenVector> &)
1609 : {
1610 3456 : return _grad_phi_face_neighbor;
1611 : }
1612 0 : VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<RealEigenVector> &)
1613 : {
1614 0 : return _second_phi_face_neighbor;
1615 : }
1616 :
1617 : template <typename OutputType>
1618 179166 : const typename OutputTools<OutputType>::VariablePhiValue & fePhi(FEType type) const
1619 : {
1620 179166 : buildFE(type);
1621 179166 : return _fe_shape_data[type]->_phi;
1622 : }
1623 :
1624 : template <typename OutputType>
1625 179361 : const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhi(FEType type) const
1626 : {
1627 179361 : buildFE(type);
1628 179361 : return _fe_shape_data[type]->_grad_phi;
1629 : }
1630 :
1631 : template <typename OutputType>
1632 172517 : const ADTemplateVariablePhiGradient<OutputType> & feADGradPhi(FEType type) const
1633 : {
1634 172517 : return _ad_grad_phi_data[type];
1635 : }
1636 :
1637 : template <typename OutputType>
1638 24509 : const typename OutputTools<OutputType>::VariablePhiSecond & feSecondPhi(FEType type) const
1639 : {
1640 24509 : _need_second_derivative.insert(type);
1641 24509 : buildFE(type);
1642 24509 : return _fe_shape_data[type]->_second_phi;
1643 : }
1644 :
1645 : template <typename OutputType>
1646 : const typename OutputTools<OutputType>::VariablePhiValue & fePhiLower(FEType type) const;
1647 :
1648 : template <typename OutputType>
1649 : const typename OutputTools<OutputType>::VariablePhiValue & feDualPhiLower(FEType type) const;
1650 :
1651 : template <typename OutputType>
1652 : const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhiLower(FEType type) const;
1653 :
1654 : template <typename OutputType>
1655 : const typename OutputTools<OutputType>::VariablePhiGradient &
1656 : feGradDualPhiLower(FEType type) const;
1657 :
1658 : template <typename OutputType>
1659 179166 : const typename OutputTools<OutputType>::VariablePhiValue & fePhiFace(FEType type) const
1660 : {
1661 179166 : buildFaceFE(type);
1662 179166 : return _fe_shape_data_face[type]->_phi;
1663 : }
1664 :
1665 : template <typename OutputType>
1666 179166 : const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhiFace(FEType type) const
1667 : {
1668 179166 : buildFaceFE(type);
1669 179166 : return _fe_shape_data_face[type]->_grad_phi;
1670 : }
1671 :
1672 : template <typename OutputType>
1673 172517 : const ADTemplateVariablePhiGradient<OutputType> & feADGradPhiFace(FEType type) const
1674 : {
1675 172517 : return _ad_grad_phi_data_face[type];
1676 : }
1677 :
1678 : template <typename OutputType>
1679 4435 : const typename OutputTools<OutputType>::VariablePhiSecond & feSecondPhiFace(FEType type) const
1680 : {
1681 4435 : _need_second_derivative.insert(type);
1682 4435 : buildFaceFE(type);
1683 4435 : return _fe_shape_data_face[type]->_second_phi;
1684 : }
1685 :
1686 : template <typename OutputType>
1687 179166 : const typename OutputTools<OutputType>::VariablePhiValue & fePhiNeighbor(FEType type) const
1688 : {
1689 179166 : buildNeighborFE(type);
1690 179166 : return _fe_shape_data_neighbor[type]->_phi;
1691 : }
1692 :
1693 : template <typename OutputType>
1694 179166 : const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhiNeighbor(FEType type) const
1695 : {
1696 179166 : buildNeighborFE(type);
1697 179166 : return _fe_shape_data_neighbor[type]->_grad_phi;
1698 : }
1699 :
1700 : template <typename OutputType>
1701 39 : const typename OutputTools<OutputType>::VariablePhiSecond & feSecondPhiNeighbor(FEType type) const
1702 : {
1703 39 : _need_second_derivative_neighbor.insert(type);
1704 39 : buildNeighborFE(type);
1705 39 : return _fe_shape_data_neighbor[type]->_second_phi;
1706 : }
1707 :
1708 : template <typename OutputType>
1709 179166 : const typename OutputTools<OutputType>::VariablePhiValue & fePhiFaceNeighbor(FEType type) const
1710 : {
1711 179166 : buildFaceNeighborFE(type);
1712 179166 : return _fe_shape_data_face_neighbor[type]->_phi;
1713 : }
1714 :
1715 : template <typename OutputType>
1716 : const typename OutputTools<OutputType>::VariablePhiGradient &
1717 179166 : feGradPhiFaceNeighbor(FEType type) const
1718 : {
1719 179166 : buildFaceNeighborFE(type);
1720 179166 : return _fe_shape_data_face_neighbor[type]->_grad_phi;
1721 : }
1722 :
1723 : template <typename OutputType>
1724 : const typename OutputTools<OutputType>::VariablePhiSecond &
1725 39 : feSecondPhiFaceNeighbor(FEType type) const
1726 : {
1727 39 : _need_second_derivative_neighbor.insert(type);
1728 39 : buildFaceNeighborFE(type);
1729 39 : return _fe_shape_data_face_neighbor[type]->_second_phi;
1730 : }
1731 :
1732 : template <typename OutputType>
1733 0 : const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhi(FEType type) const
1734 : {
1735 0 : _need_curl.insert(type);
1736 0 : buildFE(type);
1737 0 : return _fe_shape_data[type]->_curl_phi;
1738 : }
1739 :
1740 : template <typename OutputType>
1741 0 : const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhiFace(FEType type) const
1742 : {
1743 0 : _need_curl.insert(type);
1744 0 : buildFaceFE(type);
1745 0 : return _fe_shape_data_face[type]->_curl_phi;
1746 : }
1747 :
1748 : template <typename OutputType>
1749 0 : const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhiNeighbor(FEType type) const
1750 : {
1751 0 : _need_curl.insert(type);
1752 0 : buildNeighborFE(type);
1753 0 : return _fe_shape_data_neighbor[type]->_curl_phi;
1754 : }
1755 :
1756 : template <typename OutputType>
1757 0 : const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhiFaceNeighbor(FEType type) const
1758 : {
1759 0 : _need_curl.insert(type);
1760 0 : buildFaceNeighborFE(type);
1761 0 : return _fe_shape_data_face_neighbor[type]->_curl_phi;
1762 : }
1763 :
1764 : template <typename OutputType>
1765 0 : const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhi(FEType type) const
1766 : {
1767 0 : buildFE(type);
1768 0 : return _fe_shape_data[type]->_div_phi;
1769 : }
1770 :
1771 : template <typename OutputType>
1772 0 : const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhiFace(FEType type) const
1773 : {
1774 0 : buildFaceFE(type);
1775 0 : return _fe_shape_data_face[type]->_div_phi;
1776 : }
1777 :
1778 : template <typename OutputType>
1779 : const typename OutputTools<OutputType>::VariablePhiDivergence &
1780 0 : feDivPhiNeighbor(FEType type) const
1781 : {
1782 0 : buildNeighborFE(type);
1783 0 : return _fe_shape_data_neighbor[type]->_div_phi;
1784 : }
1785 :
1786 : template <typename OutputType>
1787 : const typename OutputTools<OutputType>::VariablePhiDivergence &
1788 0 : feDivPhiFaceNeighbor(FEType type) const
1789 : {
1790 0 : buildFaceNeighborFE(type);
1791 0 : return _fe_shape_data_face_neighbor[type]->_div_phi;
1792 : }
1793 :
1794 : /// On-demand computation of volume element accounting for RZ/RSpherical
1795 : Real elementVolume(const Elem * elem) const;
1796 :
1797 : /**
1798 : * Set the pointer to the XFEM controller object
1799 : */
1800 0 : void setXFEM(std::shared_ptr<XFEMInterface> xfem) { _xfem = xfem; }
1801 :
1802 : /**
1803 : * Assign the displacement numbers and directions
1804 : */
1805 : void assignDisplacements(
1806 : std::vector<std::pair<unsigned int, unsigned short>> && disp_numbers_and_directions);
1807 :
1808 : /**
1809 : * Helper function for assembling residual contriubutions on local
1810 : * quadrature points for an array kernel, bc, etc.
1811 : * @param re The local residual
1812 : * @param i The local test function index
1813 : * @param ntest The number of test functions
1814 : * @param v The residual contribution on the current qp
1815 : */
1816 112248352 : void saveLocalArrayResidual(DenseVector<Number> & re,
1817 : unsigned int i,
1818 : unsigned int ntest,
1819 : const RealEigenVector & v) const
1820 : {
1821 371932320 : for (unsigned int j = 0; j < v.size(); ++j, i += ntest)
1822 259683968 : re(i) += v(j);
1823 112248352 : }
1824 :
1825 : void saveLocalADArray(std::vector<ADReal> & re,
1826 : unsigned int i,
1827 : unsigned int ntest,
1828 : const ADRealEigenVector & v) const;
1829 :
1830 : /**
1831 : * Helper function for assembling diagonal Jacobian contriubutions on local
1832 : * quadrature points for an array kernel, bc, etc.
1833 : * @param ke The local Jacobian
1834 : * @param i The local test function index
1835 : * @param ntest The number of test functions
1836 : * @param j The local shape function index
1837 : * @param nphi The number of shape functions
1838 : * @param v The diagonal Jacobian contribution on the current qp
1839 : */
1840 31244128 : void saveDiagLocalArrayJacobian(DenseMatrix<Number> & ke,
1841 : unsigned int i,
1842 : unsigned int ntest,
1843 : unsigned int j,
1844 : unsigned int nphi,
1845 : unsigned int ivar,
1846 : const RealEigenVector & v) const
1847 : {
1848 31244128 : unsigned int pace = (_component_block_diagonal[ivar] ? 0 : nphi);
1849 97464864 : for (unsigned int k = 0; k < v.size(); ++k, i += ntest, j += pace)
1850 66220736 : ke(i, j) += v(k);
1851 31244128 : }
1852 :
1853 : /**
1854 : * Helper function for assembling full Jacobian contriubutions on local
1855 : * quadrature points for an array kernel, bc, etc.
1856 : * @param ke The local Jacobian
1857 : * @param i The local test function index
1858 : * @param ntest The number of test functions
1859 : * @param j The local shape function index
1860 : * @param nphi The number of shape functions
1861 : * @param ivar The array variable index
1862 : * @param jvar The contributing variable index
1863 : * @param v The full Jacobian contribution from a variable on the current qp
1864 : */
1865 52730848 : void saveFullLocalArrayJacobian(DenseMatrix<Number> & ke,
1866 : unsigned int i,
1867 : unsigned int ntest,
1868 : unsigned int j,
1869 : unsigned int nphi,
1870 : unsigned int ivar,
1871 : unsigned int jvar,
1872 : const RealEigenMatrix & v) const
1873 : {
1874 52730848 : if (ivar == jvar && _component_block_diagonal[ivar])
1875 : {
1876 148032 : for (unsigned int k = 0; k < v.rows(); ++k, i += ntest)
1877 95520 : ke(i, j) += v(k, k);
1878 : }
1879 : else
1880 : {
1881 52678336 : const unsigned int saved_j = j;
1882 236755008 : for (unsigned int k = 0; k < v.rows(); ++k, i += ntest)
1883 : {
1884 184076672 : j = saved_j;
1885 866815104 : for (unsigned int l = 0; l < v.cols(); ++l, j += nphi)
1886 682738432 : ke(i, j) += v(k, l);
1887 : }
1888 : }
1889 52730848 : }
1890 :
1891 2399 : DenseVector<Real> getJacobianDiagonal(DenseMatrix<Number> & ke)
1892 : {
1893 2399 : unsigned int rows = ke.m();
1894 2399 : unsigned int cols = ke.n();
1895 2399 : DenseVector<Real> diag(rows);
1896 16347 : for (unsigned int i = 0; i < rows; i++)
1897 : // % operation is needed to account for cases of no component coupling of array variables
1898 13948 : diag(i) = ke(i, i % cols);
1899 2399 : return diag;
1900 : }
1901 :
1902 : /**
1903 : * Attaches the current elem/volume quadrature rule to the given fe. The
1904 : * current subdomain (as set via setCurrentSubdomainID is used to determine
1905 : * the correct rule. The attached quadrature rule is also returned.
1906 : */
1907 : inline const libMesh::QBase * attachQRuleElem(unsigned int dim, FEBase & fe)
1908 : {
1909 : auto qrule = qrules(dim).vol.get();
1910 : fe.attach_quadrature_rule(qrule);
1911 : return qrule;
1912 : }
1913 :
1914 : /**
1915 : * Attaches the current face/area quadrature rule to the given fe. The
1916 : * current subdomain (as set via setCurrentSubdomainID is used to determine
1917 : * the correct rule. The attached quadrature rule is also returned.
1918 : */
1919 : inline const libMesh::QBase * attachQRuleFace(unsigned int dim, FEBase & fe)
1920 : {
1921 : auto qrule = qrules(dim).face.get();
1922 : fe.attach_quadrature_rule(qrule);
1923 : return qrule;
1924 : }
1925 :
1926 : /**
1927 : * signals this object that a vector containing variable scaling factors should be used when
1928 : * doing residual and matrix assembly
1929 : */
1930 : void hasScalingVector();
1931 :
1932 : /**
1933 : * Modify the weights when using the arbitrary quadrature rule. The intention is to use this when
1934 : * you wish to supply your own quadrature after calling reinit at physical points.
1935 : *
1936 : * You should only use this if the arbitrary quadrature is the current quadrature rule!
1937 : *
1938 : * @param weights The weights to fill into _current_JxW
1939 : */
1940 : void modifyArbitraryWeights(const std::vector<Real> & weights);
1941 :
1942 : /**
1943 : * @return whether we are computing a residual
1944 : */
1945 30271857 : bool computingResidual() const { return _computing_residual; }
1946 :
1947 : /**
1948 : * @return whether we are computing a Jacobian
1949 : */
1950 26443601 : bool computingJacobian() const { return _computing_jacobian; }
1951 :
1952 : /**
1953 : * @return whether we are computing a residual and a Jacobian simultaneously
1954 : */
1955 : bool computingResidualAndJacobian() const { return _computing_residual_and_jacobian; }
1956 :
1957 : /**
1958 : * @return The current mortar segment element
1959 : */
1960 1326 : const Elem * const & msmElem() const { return _msm_elem; }
1961 :
1962 : /**
1963 : * Indicate that we have p-refinement. This method will perform the following tasks:
1964 : * - Disable p-refinement as requested by the user with \p disable_p_refinement_for_families
1965 : * -.Disable p-refinement of Lagrange helper types that we use for getting things like the
1966 : * physical locations of quadrature points and JxW. (Don't worry, we still use the element
1967 : * p-level when initializing the quadrature rule attached to the Lagrange helper so the number
1968 : * of quadrature points reflects the element p-level)
1969 : * @param disable_p_refinement_for_families Families that we should disable p-refinement for
1970 : */
1971 : void havePRefinement(const std::unordered_set<FEFamily> & disable_p_refinement_for_families);
1972 :
1973 : /**
1974 : * Set the current lower dimensional element. This may be null
1975 : */
1976 : void setCurrentLowerDElem(const Elem * const lower_d_elem);
1977 :
1978 : private:
1979 : /**
1980 : * Just an internal helper function to reinit the volume FE objects.
1981 : *
1982 : * @param elem The element we are using to reinit
1983 : */
1984 : void reinitFE(const Elem * elem);
1985 :
1986 : /**
1987 : * Just an internal helper function to reinit the face FE objects.
1988 : *
1989 : * @param elem The element we are using to reinit
1990 : * @param side The side of the element we are reiniting on
1991 : */
1992 : void reinitFEFace(const Elem * elem, unsigned int side);
1993 :
1994 : void computeFaceMap(const Elem & elem, const unsigned int side, const std::vector<Real> & qw);
1995 :
1996 : void reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);
1997 :
1998 : void reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);
1999 :
2000 : template <typename Points, typename Coords>
2001 : void setCoordinateTransformation(const libMesh::QBase * qrule,
2002 : const Points & q_points,
2003 : Coords & coord,
2004 : SubdomainID sub_id);
2005 :
2006 : void computeCurrentElemVolume();
2007 :
2008 : void computeCurrentFaceVolume();
2009 :
2010 : void computeCurrentNeighborVolume();
2011 :
2012 : /**
2013 : * Update the integration weights for XFEM partial elements.
2014 : * This only affects the weights if XFEM is used and if the element is cut.
2015 : * @param elem The element for which the weights are adjusted
2016 : */
2017 : void modifyWeightsDueToXFEM(const Elem * elem);
2018 :
2019 : /**
2020 : * Update the face integration weights for XFEM partial elements.
2021 : * This only affects the weights if XFEM is used and if the element is cut.
2022 : * @param elem The element for which the weights are adjusted
2023 : * @param side The side of element for which the weights are adjusted
2024 : */
2025 : void modifyFaceWeightsDueToXFEM(const Elem * elem, unsigned int side = 0);
2026 :
2027 : /**
2028 : * compute gradient of phi possibly with derivative information with respect to nonlinear
2029 : * displacement variables
2030 : */
2031 : template <typename OutputType>
2032 : void computeGradPhiAD(const Elem * elem,
2033 : unsigned int n_qp,
2034 : ADTemplateVariablePhiGradient<OutputType> & grad_phi,
2035 : libMesh::FEGenericBase<OutputType> * fe);
2036 :
2037 : /**
2038 : * resize any objects that contribute to automatic differentiation-related mapping calculations
2039 : */
2040 : void resizeADMappingObjects(unsigned int n_qp, unsigned int dim);
2041 :
2042 : /**
2043 : * compute the finite element reference-physical mapping quantities (such as JxW) with possible
2044 : * dependence on nonlinear displacement variables at a single quadrature point
2045 : */
2046 : void
2047 : computeSinglePointMapAD(const Elem * elem, const std::vector<Real> & qw, unsigned p, FEBase * fe);
2048 :
2049 : /**
2050 : * Add local residuals of all field variables for a tag onto the tag's residual vector
2051 : */
2052 : void addResidual(const VectorTag & vector_tag);
2053 : /**
2054 : * Add local neighbor residuals of all field variables for a tag onto the tag's residual vector
2055 : */
2056 : void addResidualNeighbor(const VectorTag & vector_tag);
2057 : /**
2058 : * Add local lower-dimensional block residuals of all field variables for a tag onto the tag's
2059 : * residual vector
2060 : */
2061 : void addResidualLower(const VectorTag & vector_tag);
2062 : /**
2063 : * Add residuals of all scalar variables for a tag onto the tag's residual vector
2064 : */
2065 : void addResidualScalar(const VectorTag & vector_tag);
2066 :
2067 : /**
2068 : * Clears all of the cached residuals for a specific vector tag
2069 : */
2070 : void clearCachedResiduals(const VectorTag & vector_tag);
2071 :
2072 : /**
2073 : * Cache individual residual contributions. These will ultimately get added to the residual when
2074 : * addCachedResidual() is called.
2075 : *
2076 : * @param dof The degree of freedom to add the residual contribution to
2077 : * @param value The value of the residual contribution.
2078 : * @param TagID the contribution should go to this tagged residual
2079 : */
2080 : void cacheResidual(dof_id_type dof, Real value, TagID tag_id);
2081 :
2082 : /**
2083 : * Cache individual residual contributions. These will ultimately get added to the residual when
2084 : * addCachedResidual() is called.
2085 : *
2086 : * @param dof The degree of freedom to add the residual contribution to
2087 : * @param value The value of the residual contribution.
2088 : * @param tags the contribution should go to all these tags
2089 : */
2090 : void cacheResidual(dof_id_type dof, Real value, const std::set<TagID> & tags);
2091 :
2092 : /**
2093 : * Appling scaling, constraints to the local residual block and populate the full DoF indices
2094 : * for array variable.
2095 : */
2096 : void processLocalResidual(DenseVector<Number> & res_block,
2097 : std::vector<dof_id_type> & dof_indices,
2098 : const std::vector<Real> & scaling_factor);
2099 :
2100 : /**
2101 : * Add a local residual block to a global residual vector with proper scaling.
2102 : */
2103 : void addResidualBlock(NumericVector<Number> & residual,
2104 : DenseVector<Number> & res_block,
2105 : const std::vector<dof_id_type> & dof_indices,
2106 : const std::vector<Real> & scaling_factor);
2107 :
2108 : /**
2109 : * Push a local residual block with proper scaling into cache.
2110 : */
2111 : void cacheResidualBlock(std::vector<Real> & cached_residual_values,
2112 : std::vector<dof_id_type> & cached_residual_rows,
2113 : DenseVector<Number> & res_block,
2114 : const std::vector<dof_id_type> & dof_indices,
2115 : const std::vector<Real> & scaling_factor);
2116 :
2117 : /**
2118 : * Set a local residual block to a global residual vector with proper scaling.
2119 : */
2120 : void setResidualBlock(NumericVector<Number> & residual,
2121 : DenseVector<Number> & res_block,
2122 : const std::vector<dof_id_type> & dof_indices,
2123 : const std::vector<Real> & scaling_factor);
2124 :
2125 : /**
2126 : * Add a local Jacobian block to a global Jacobian with proper scaling.
2127 : */
2128 : void addJacobianBlock(libMesh::SparseMatrix<Number> & jacobian,
2129 : DenseMatrix<Number> & jac_block,
2130 : const MooseVariableBase & ivar,
2131 : const MooseVariableBase & jvar,
2132 : const std::vector<dof_id_type> & idof_indices,
2133 : const std::vector<dof_id_type> & jdof_indices);
2134 :
2135 : /**
2136 : * Push a local Jacobian block with proper scaling into cache for a certain tag.
2137 : */
2138 : void cacheJacobianBlock(DenseMatrix<Number> & jac_block,
2139 : const MooseVariableBase & ivar,
2140 : const MooseVariableBase & jvar,
2141 : const std::vector<dof_id_type> & idof_indices,
2142 : const std::vector<dof_id_type> & jdof_indices,
2143 : TagID tag);
2144 :
2145 : /**
2146 : * Push non-zeros of a local Jacobian block with proper scaling into cache for a certain tag.
2147 : */
2148 : void cacheJacobianBlockNonzero(DenseMatrix<Number> & jac_block,
2149 : const MooseVariableBase & ivar,
2150 : const MooseVariableBase & jvar,
2151 : const std::vector<dof_id_type> & idof_indices,
2152 : const std::vector<dof_id_type> & jdof_indices,
2153 : TagID tag);
2154 :
2155 : /**
2156 : * Adds element matrices for ivar rows and jvar columns to the global Jacobian matrices.
2157 : */
2158 : void addJacobianCoupledVarPair(const MooseVariableBase & ivar, const MooseVariableBase & jvar);
2159 :
2160 : /**
2161 : * Caches element matrix for ivar rows and jvar columns
2162 : */
2163 : void cacheJacobianCoupledVarPair(const MooseVariableBase & ivar, const MooseVariableBase & jvar);
2164 :
2165 : /**
2166 : * Clear any currently cached jacobians
2167 : *
2168 : * This is automatically called by setCachedJacobian
2169 : */
2170 : void clearCachedJacobian();
2171 :
2172 : /**
2173 : * Build FEs with a type
2174 : * @param type The type of FE
2175 : */
2176 : void buildFE(FEType type) const;
2177 :
2178 : /**
2179 : * Build FEs for a face with a type
2180 : * @param type The type of FE
2181 : */
2182 : void buildFaceFE(FEType type) const;
2183 :
2184 : /**
2185 : * Build FEs for a neighbor with a type
2186 : * @param type The type of FE
2187 : */
2188 : void buildNeighborFE(FEType type) const;
2189 :
2190 : /**
2191 : * Build FEs for a neighbor face with a type
2192 : * @param type The type of FE
2193 : */
2194 : void buildFaceNeighborFE(FEType type) const;
2195 :
2196 : /**
2197 : * Build FEs for a lower dimensional element with a type
2198 : * @param type The type of FE
2199 : */
2200 : void buildLowerDFE(FEType type) const;
2201 :
2202 : void buildLowerDDualFE(FEType type) const;
2203 :
2204 : /**
2205 : * Build Vector FEs with a type
2206 : * @param type The type of FE
2207 : */
2208 : void buildVectorFE(FEType type) const;
2209 :
2210 : /**
2211 : * Build Vector FEs for a face with a type
2212 : * @param type The type of FE
2213 : */
2214 : void buildVectorFaceFE(FEType type) const;
2215 :
2216 : /**
2217 : * Build Vector FEs for a neighbor with a type
2218 : * @param type The type of FE
2219 : */
2220 : void buildVectorNeighborFE(FEType type) const;
2221 :
2222 : /**
2223 : * Build Vector FEs for a neighbor face with a type
2224 : * @param type The type of FE
2225 : */
2226 : void buildVectorFaceNeighborFE(FEType type) const;
2227 :
2228 : /**
2229 : * Build Vector FEs for a lower dimensional element with a type
2230 : * @param type The type of FE
2231 : */
2232 : void buildVectorLowerDFE(FEType type) const;
2233 : void buildVectorDualLowerDFE(FEType type) const;
2234 :
2235 : /**
2236 : * Sets whether or not Jacobian coupling between \p ivar and \p jvar is used
2237 : * to the value \p used
2238 : */
2239 654115478 : void jacobianBlockUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
2240 : {
2241 654115478 : _jacobian_block_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
2242 654115478 : }
2243 :
2244 : /**
2245 : * Return a flag to indicate if a particular coupling Jacobian block
2246 : * between \p ivar and \p jvar is used
2247 : */
2248 143694582 : char jacobianBlockUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
2249 : {
2250 143694582 : return _jacobian_block_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
2251 : }
2252 :
2253 : /**
2254 : * Sets whether or not neighbor Jacobian coupling between \p ivar and \p jvar is used
2255 : * to the value \p used
2256 : */
2257 291396769 : void jacobianBlockNeighborUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
2258 : {
2259 291396769 : _jacobian_block_neighbor_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
2260 291396769 : }
2261 :
2262 : /**
2263 : * Return a flag to indicate if a particular coupling neighbor Jacobian block
2264 : * between \p ivar and \p jvar is used
2265 : */
2266 462592 : char jacobianBlockNeighborUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
2267 : {
2268 462592 : return _jacobian_block_neighbor_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
2269 : }
2270 :
2271 : /**
2272 : * Sets whether or not lower Jacobian coupling between \p ivar and \p jvar is used
2273 : * to the value \p used
2274 : */
2275 12316684 : void jacobianBlockLowerUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
2276 : {
2277 12316684 : _jacobian_block_lower_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
2278 12316684 : }
2279 :
2280 : /**
2281 : * Return a flag to indicate if a particular coupling lower Jacobian block
2282 : * between \p ivar and \p jvar is used
2283 : */
2284 948462 : char jacobianBlockLowerUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
2285 : {
2286 948462 : return _jacobian_block_lower_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
2287 : }
2288 :
2289 : /**
2290 : * Sets whether or not nonlocal Jacobian coupling between \p ivar and \p jvar is used
2291 : * to the value \p used
2292 : */
2293 37420 : void jacobianBlockNonlocalUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
2294 : {
2295 37420 : _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
2296 37420 : }
2297 :
2298 : /**
2299 : * Return a flag to indicate if a particular coupling nonlocal Jacobian block
2300 : * between \p ivar and \p jvar is used
2301 : */
2302 9016 : char jacobianBlockNonlocalUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
2303 : {
2304 9016 : return _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
2305 : }
2306 :
2307 : /**
2308 : * request phi, dphi, xyz, JxW, etc. data through the FE helper functions
2309 : */
2310 : void helpersRequestData();
2311 :
2312 : SystemBase & _sys;
2313 : SubProblem & _subproblem;
2314 :
2315 : const bool _displaced;
2316 :
2317 : /// Coupling matrices
2318 : const libMesh::CouplingMatrix * _cm;
2319 : const libMesh::CouplingMatrix & _nonlocal_cm;
2320 :
2321 : /// Whether we are currently computing the residual
2322 : const bool & _computing_residual;
2323 :
2324 : /// Whether we are currently computing the Jacobian
2325 : const bool & _computing_jacobian;
2326 :
2327 : /// Whether we are currently computing the residual and Jacobian
2328 : const bool & _computing_residual_and_jacobian;
2329 :
2330 : /// Entries in the coupling matrix for field variables
2331 : std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> _cm_ff_entry;
2332 : /// Entries in the coupling matrix for field variables vs scalar variables
2333 : std::vector<std::pair<MooseVariableFieldBase *, MooseVariableScalar *>> _cm_fs_entry;
2334 : /// Entries in the coupling matrix for scalar variables vs field variables
2335 : std::vector<std::pair<MooseVariableScalar *, MooseVariableFieldBase *>> _cm_sf_entry;
2336 : /// Entries in the coupling matrix for scalar variables
2337 : std::vector<std::pair<MooseVariableScalar *, MooseVariableScalar *>> _cm_ss_entry;
2338 : /// Entries in the coupling matrix for field variables for nonlocal calculations
2339 : std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> _cm_nonlocal_entry;
2340 : /// Flag that indicates if the jacobian block was used
2341 : std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_used;
2342 : std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_nonlocal_used;
2343 : /// Flag that indicates if the jacobian block for neighbor was used
2344 : std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_neighbor_used;
2345 : /// Flag that indicates if the jacobian block for the lower dimensional element was used
2346 : std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_lower_used;
2347 : /// DOF map
2348 : const libMesh::DofMap & _dof_map;
2349 : /// Thread number (id)
2350 : THREAD_ID _tid;
2351 :
2352 : MooseMesh & _mesh;
2353 :
2354 : unsigned int _mesh_dimension;
2355 :
2356 : /// The finite element type of the FE helper classes. The helper class gives us data like JxW, the
2357 : /// physical quadrature point locations, etc.
2358 : const FEType _helper_type;
2359 :
2360 : /// Whether user code requested a \p FEType the same as our \p _helper_type
2361 : mutable bool _user_added_fe_of_helper_type;
2362 : mutable bool _user_added_fe_face_of_helper_type;
2363 : mutable bool _user_added_fe_face_neighbor_of_helper_type;
2364 : mutable bool _user_added_fe_neighbor_of_helper_type;
2365 : mutable bool _user_added_fe_lower_of_helper_type;
2366 :
2367 : /// Containers for holding unique FE helper types if we are doing p-refinement. If we are not
2368 : /// doing p-refinement then the helper data is owned by the \p _fe data members
2369 : std::vector<std::unique_ptr<FEBase>> _unique_fe_helper;
2370 : std::vector<std::unique_ptr<FEBase>> _unique_fe_face_helper;
2371 : std::vector<std::unique_ptr<FEBase>> _unique_fe_face_neighbor_helper;
2372 : std::vector<std::unique_ptr<FEBase>> _unique_fe_neighbor_helper;
2373 : std::vector<std::unique_ptr<FEBase>> _unique_fe_lower_helper;
2374 :
2375 : /// Whether we are currently building the FE classes for the helpers
2376 : bool _building_helpers;
2377 :
2378 : /// The XFEM controller
2379 : std::shared_ptr<XFEMInterface> _xfem;
2380 :
2381 : /// The "volume" fe object that matches the current elem
2382 : std::map<FEType, FEBase *> _current_fe;
2383 : /// The "face" fe object that matches the current elem
2384 : std::map<FEType, FEBase *> _current_fe_face;
2385 : /// The "neighbor" fe object that matches the current elem
2386 : std::map<FEType, FEBase *> _current_fe_neighbor;
2387 : /// The "neighbor face" fe object that matches the current elem
2388 : std::map<FEType, FEBase *> _current_fe_face_neighbor;
2389 :
2390 : /// The "volume" vector fe object that matches the current elem
2391 : std::map<FEType, FEVectorBase *> _current_vector_fe;
2392 : /// The "face" vector fe object that matches the current elem
2393 : std::map<FEType, FEVectorBase *> _current_vector_fe_face;
2394 : /// The "neighbor" vector fe object that matches the current elem
2395 : std::map<FEType, FEVectorBase *> _current_vector_fe_neighbor;
2396 : /// The "neighbor face" vector fe object that matches the current elem
2397 : std::map<FEType, FEVectorBase *> _current_vector_fe_face_neighbor;
2398 :
2399 : /**** Volume Stuff ****/
2400 :
2401 : /// Each dimension's actual fe objects indexed on type
2402 : mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe;
2403 : /// Each dimension's actual vector fe objects indexed on type
2404 : mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe;
2405 : /// Each dimension's helper objects
2406 : std::map<unsigned int, FEBase *> _holder_fe_helper;
2407 : /// The current helper object for transforming coordinates
2408 : FEBase * _current_fe_helper;
2409 : /// The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac kernels)
2410 : libMesh::QBase * _current_qrule;
2411 : /// The current volumetric quadrature for the element
2412 : libMesh::QBase * _current_qrule_volume;
2413 : /// The current arbitrary quadrature rule used within the element interior
2414 : ArbitraryQuadrature * _current_qrule_arbitrary;
2415 : /// The current arbitrary quadrature rule used on the element face
2416 : ArbitraryQuadrature * _current_qrule_arbitrary_face;
2417 : /// The current list of quadrature points
2418 : MooseArray<Point> _current_q_points;
2419 : /// The current list of transformed jacobian weights
2420 : MooseArray<Real> _current_JxW;
2421 : /// The coordinate system
2422 : Moose::CoordinateSystemType _coord_type;
2423 : /// The current coordinate transformation coefficients
2424 : MooseArray<Real> _coord;
2425 : /// The AD version of the current coordinate transformation coefficients
2426 : MooseArray<ADReal> _ad_coord;
2427 :
2428 : /// Data structure for tracking/grouping a set of quadrature rules for a
2429 : /// particular dimensionality of mesh element.
2430 : struct QRules
2431 : {
2432 204588 : QRules()
2433 204588 : : vol(nullptr),
2434 204588 : face(nullptr),
2435 204588 : arbitrary_vol(nullptr),
2436 204588 : arbitrary_face(nullptr),
2437 409176 : neighbor(nullptr)
2438 : {
2439 204588 : }
2440 :
2441 : /// volume/elem (meshdim) quadrature rule
2442 : std::unique_ptr<libMesh::QBase> vol;
2443 : /// area/face (meshdim-1) quadrature rule
2444 : std::unique_ptr<libMesh::QBase> face;
2445 : /// finite volume face/flux quadrature rule (meshdim-1)
2446 : std::unique_ptr<libMesh::QBase> fv_face;
2447 : /// volume/elem (meshdim) custom points quadrature rule
2448 : std::unique_ptr<ArbitraryQuadrature> arbitrary_vol;
2449 : /// area/face (meshdim-1) custom points quadrature rule
2450 : std::unique_ptr<ArbitraryQuadrature> arbitrary_face;
2451 : /// area/face (meshdim-1) custom points quadrature rule for DG
2452 : std::unique_ptr<ArbitraryQuadrature> neighbor;
2453 : };
2454 :
2455 : /// Holds quadrature rules for each dimension. These are created up front
2456 : /// at the start of the simulation and reused/referenced for the remainder of
2457 : /// the sim. This data structure should generally be read/accessed via the
2458 : /// qrules() function.
2459 : std::unordered_map<SubdomainID, std::vector<QRules>> _qrules;
2460 :
2461 : /// This is an abstraction over the internal qrules function. This is
2462 : /// necessary for faces because (nodes of) faces can exists in more than one
2463 : /// subdomain. When this is the case, we need to use the quadrature rule from
2464 : /// the subdomain that has the highest specified quadrature order. So when
2465 : /// you need to access a face quadrature rule, you should retrieve it via this
2466 : /// function.
2467 : libMesh::QBase * qruleFace(const Elem * elem, unsigned int side);
2468 : ArbitraryQuadrature * qruleArbitraryFace(const Elem * elem, unsigned int side);
2469 :
2470 : template <typename T>
2471 9287907 : T * qruleFaceHelper(const Elem * elem, unsigned int side, std::function<T *(QRules &)> rule_fn)
2472 : {
2473 9287907 : auto dim = elem->dim();
2474 9287907 : auto neighbor = elem->neighbor_ptr(side);
2475 9287907 : auto q = rule_fn(qrules(dim, elem->subdomain_id()));
2476 9287907 : if (!neighbor)
2477 4136890 : return q;
2478 :
2479 : // find the maximum face quadrature order for all blocks the face is in
2480 5151017 : auto neighbor_block = neighbor->subdomain_id();
2481 5151017 : if (neighbor_block == elem->subdomain_id())
2482 4961763 : return q;
2483 :
2484 189254 : auto q_neighbor = rule_fn(qrules(dim, neighbor_block));
2485 189254 : if (q->get_order() > q_neighbor->get_order())
2486 343 : return q;
2487 188911 : return q_neighbor;
2488 : }
2489 :
2490 404644994 : inline QRules & qrules(unsigned int dim) { return qrules(dim, _current_subdomain_id); }
2491 :
2492 : /// This is a helper function for accessing quadrature rules for a
2493 : /// particular dimensionality of element. All access to quadrature rules in
2494 : /// Assembly should be done via this accessor function.
2495 434443102 : inline QRules & qrules(unsigned int dim, SubdomainID block)
2496 : {
2497 434443102 : if (_qrules.find(block) == _qrules.end())
2498 : {
2499 : mooseAssert(_qrules.find(Moose::ANY_BLOCK_ID) != _qrules.end(),
2500 : "missing quadrature rules for specified block");
2501 : mooseAssert(_qrules[Moose::ANY_BLOCK_ID].size() > dim,
2502 : "quadrature rules not sized property for dimension");
2503 434441355 : return _qrules[Moose::ANY_BLOCK_ID][dim];
2504 : }
2505 : mooseAssert(_qrules.find(block) != _qrules.end(),
2506 : "missing quadrature rules for specified block");
2507 : mooseAssert(_qrules[block].size() > dim, "quadrature rules not sized property for dimension");
2508 1747 : return _qrules[block][dim];
2509 : }
2510 :
2511 : /**** Face Stuff ****/
2512 :
2513 : /// types of finite elements
2514 : mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_face;
2515 : /// types of vector finite elements
2516 : mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face;
2517 : /// Each dimension's helper objects
2518 : std::map<unsigned int, FEBase *> _holder_fe_face_helper;
2519 : /// helper object for transforming coordinates
2520 : FEBase * _current_fe_face_helper;
2521 : /// quadrature rule used on faces
2522 : libMesh::QBase * _current_qrule_face;
2523 : /// The current arbitrary quadrature rule used on element faces
2524 : ArbitraryQuadrature * _current_qface_arbitrary;
2525 : /// The current quadrature points on a face
2526 : MooseArray<Point> _current_q_points_face;
2527 : /// The current transformed jacobian weights on a face
2528 : MooseArray<Real> _current_JxW_face;
2529 : /// The current Normal vectors at the quadrature points.
2530 : MooseArray<Point> _current_normals;
2531 : /// Mapped normals
2532 : std::vector<Eigen::Map<RealDIMValue>> _mapped_normals;
2533 : /// The current tangent vectors at the quadrature points
2534 : MooseArray<std::vector<Point>> _current_tangents;
2535 :
2536 : /// Extra element IDs
2537 : std::vector<dof_id_type> _extra_elem_ids;
2538 : /// Extra element IDs of neighbor
2539 : std::vector<dof_id_type> _neighbor_extra_elem_ids;
2540 : /// Holds pointers to the dimension's normal vectors
2541 : std::map<unsigned int, const std::vector<Point> *> _holder_normals;
2542 :
2543 : /**** Neighbor Stuff ****/
2544 :
2545 : /// types of finite elements
2546 : mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_neighbor;
2547 : mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_face_neighbor;
2548 : mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_neighbor;
2549 : mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face_neighbor;
2550 :
2551 : /// Each dimension's helper objects
2552 : std::map<unsigned int, FEBase *> _holder_fe_neighbor_helper;
2553 : std::map<unsigned int, FEBase *> _holder_fe_face_neighbor_helper;
2554 :
2555 : /// FE objects for lower dimensional elements
2556 : mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_lower;
2557 : /// Vector FE objects for lower dimensional elements
2558 : mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_lower;
2559 : /// helper object for transforming coordinates for lower dimensional element quadrature points
2560 : std::map<unsigned int, FEBase *> _holder_fe_lower_helper;
2561 :
2562 : /// quadrature rule used on neighbors
2563 : libMesh::QBase * _current_qrule_neighbor;
2564 : /// The current quadrature points on the neighbor face
2565 : MooseArray<Point> _current_q_points_face_neighbor;
2566 : /// Flag to indicate that JxW_neighbor is needed
2567 : mutable bool _need_JxW_neighbor;
2568 : /// The current transformed jacobian weights on a neighbor's face
2569 : MooseArray<Real> _current_JxW_neighbor;
2570 : /// The current coordinate transformation coefficients
2571 : MooseArray<Real> _coord_neighbor;
2572 : /// The coordinate transformation coefficients evaluated on the quadrature points of the mortar
2573 : /// segment mesh
2574 : MooseArray<Real> _coord_msm;
2575 :
2576 : /********** mortar stuff *************/
2577 :
2578 : /// A JxW for working on mortar segement elements
2579 : const std::vector<Real> * _JxW_msm;
2580 : /// A FE object for working on mortar segement elements
2581 : std::unique_ptr<FEBase> _fe_msm;
2582 : /// A qrule object for working on mortar segement elements. This needs to be a
2583 : /// raw pointer because we need to be able to return a reference to it because
2584 : /// we will be constructing other objects that need the qrule before the qrule
2585 : /// is actually created
2586 : libMesh::QBase * _qrule_msm;
2587 : /// Flag specifying whether a custom quadrature rule has been specified for mortar segment mesh
2588 : bool _custom_mortar_qrule;
2589 :
2590 : /// quadrature rule used on lower dimensional elements. This should always be
2591 : /// the same as the face qrule
2592 : libMesh::QBase * _current_qrule_lower;
2593 :
2594 : protected:
2595 : /// The current "element" we are currently on.
2596 : const Elem * _current_elem;
2597 : /// The current subdomain ID
2598 : SubdomainID _current_subdomain_id;
2599 : /// The current boundary ID
2600 : BoundaryID _current_boundary_id;
2601 : /// Volume of the current element
2602 : Real _current_elem_volume;
2603 : /// The current side of the selected element (valid only when working with sides)
2604 : unsigned int _current_side;
2605 : /// The current "element" making up the side we are currently on.
2606 : const Elem * _current_side_elem;
2607 : /// Volume of the current side element
2608 : Real _current_side_volume;
2609 : /// The current neighbor "element"
2610 : const Elem * _current_neighbor_elem;
2611 : /// The current neighbor subdomain ID
2612 : SubdomainID _current_neighbor_subdomain_id;
2613 : /// The current side of the selected neighboring element (valid only when working with sides)
2614 : unsigned int _current_neighbor_side;
2615 : /// The current side element of the ncurrent neighbor element
2616 : const Elem * _current_neighbor_side_elem;
2617 : /// true is apps need to compute neighbor element volume
2618 : mutable bool _need_neighbor_elem_volume;
2619 : /// Volume of the current neighbor
2620 : Real _current_neighbor_volume;
2621 : /// The current node we are working with
2622 : const Node * _current_node;
2623 : /// The current neighboring node we are working with
2624 : const Node * _current_neighbor_node;
2625 : /// Boolean to indicate whether current element volumes has been computed
2626 : bool _current_elem_volume_computed;
2627 : /// Boolean to indicate whether current element side volumes has been computed
2628 : bool _current_side_volume_computed;
2629 :
2630 : /// The current lower dimensional element
2631 : const Elem * _current_lower_d_elem;
2632 : /// The current neighboring lower dimensional element
2633 : const Elem * _current_neighbor_lower_d_elem;
2634 : /// Whether we need to compute the lower dimensional element volume
2635 : mutable bool _need_lower_d_elem_volume;
2636 : /// The current lower dimensional element volume
2637 : Real _current_lower_d_elem_volume;
2638 : /// Whether we need to compute the neighboring lower dimensional element volume
2639 : mutable bool _need_neighbor_lower_d_elem_volume;
2640 : /// The current neighboring lower dimensional element volume
2641 : Real _current_neighbor_lower_d_elem_volume;
2642 : /// Whether dual shape functions need to be computed for mortar constraints
2643 : bool _need_dual;
2644 :
2645 : /// This will be filled up with the physical points passed into reinitAtPhysical() if it is called. Invalid at all other times.
2646 : MooseArray<Point> _current_physical_points;
2647 :
2648 : /*
2649 : * Residual contributions <tag_index, ivar>
2650 : *
2651 : * tag_index is the index into _residual_vector_tags, that is, _sub_Re[0] corresponds to the tag
2652 : * with TagID _residual_vector_tags[0]._id
2653 : *
2654 : * When ivar corresponds to an array variable, the dense vector is in size of ndof * count,
2655 : * where count is the number of components of the array variable. The local residual is ordered
2656 : * as (r_i,j, i = 1,...,ndof; j = 1,...,count).
2657 : *
2658 : * Dense vectors for variables (ivar+i, i = 1,...,count) are empty.
2659 : */
2660 : std::vector<std::vector<DenseVector<Number>>> _sub_Re;
2661 : std::vector<std::vector<DenseVector<Number>>> _sub_Rn;
2662 : /// residual contributions for each variable from the lower dimensional element
2663 : std::vector<std::vector<DenseVector<Number>>> _sub_Rl;
2664 :
2665 : /// auxiliary vector for scaling residuals (optimization to avoid expensive construction/destruction)
2666 : DenseVector<Number> _tmp_Re;
2667 :
2668 : /*
2669 : * Jacobian contributions <Tag, ivar, jvar>
2670 : * When ivar corresponds to an array variable, the number of rows of the dense matrix is in size
2671 : * of indof * icount, where icount is the number of components of ivar. When jvar corresponds to
2672 : * an array variable, the number of columns of the dense matrix is in size of jndof * jcount,
2673 : * where jcount is the number of components of jvar. The local residual is ordered as
2674 : * (K_(i,j,k,l), k=1,...,jndof; l = 1,...,jcout; i = 1,...,indof; j = 1,...,icount).
2675 : *
2676 : * Dense matrices for variables (ivar+i, i = 1,...,icount) or (jvar+j, j = 1,...,jcount) are
2677 : * empty.
2678 : */
2679 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kee;
2680 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Keg;
2681 :
2682 : /// jacobian contributions from the element and neighbor <Tag, ivar, jvar>
2683 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Ken;
2684 : /// jacobian contributions from the neighbor and element <Tag, ivar, jvar>
2685 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kne;
2686 : /// jacobian contributions from the neighbor <Tag, ivar, jvar>
2687 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Knn;
2688 : /// dlower/dlower
2689 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kll;
2690 : /// dlower/dsecondary (or dlower/delement)
2691 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kle;
2692 : /// dlower/dprimary (or dlower/dneighbor)
2693 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kln;
2694 : /// dsecondary/dlower (or delement/dlower)
2695 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kel;
2696 : /// dprimary/dlower (or dneighbor/dlower)
2697 : std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Knl;
2698 :
2699 : /// auxiliary matrix for scaling jacobians (optimization to avoid expensive construction/destruction)
2700 : DenseMatrix<Number> _tmp_Ke;
2701 :
2702 : // Shape function values, gradients. second derivatives
2703 : VariablePhiValue _phi;
2704 : VariablePhiGradient _grad_phi;
2705 : VariablePhiSecond _second_phi;
2706 :
2707 : VariablePhiValue _phi_face;
2708 : VariablePhiGradient _grad_phi_face;
2709 : VariablePhiSecond _second_phi_face;
2710 :
2711 : VariablePhiValue _phi_neighbor;
2712 : VariablePhiGradient _grad_phi_neighbor;
2713 : VariablePhiSecond _second_phi_neighbor;
2714 :
2715 : VariablePhiValue _phi_face_neighbor;
2716 : VariablePhiGradient _grad_phi_face_neighbor;
2717 : VariablePhiSecond _second_phi_face_neighbor;
2718 :
2719 : // Shape function values, gradients, second derivatives
2720 : VectorVariablePhiValue _vector_phi;
2721 : VectorVariablePhiGradient _vector_grad_phi;
2722 : VectorVariablePhiSecond _vector_second_phi;
2723 : VectorVariablePhiCurl _vector_curl_phi;
2724 : VectorVariablePhiDivergence _vector_div_phi;
2725 :
2726 : VectorVariablePhiValue _vector_phi_face;
2727 : VectorVariablePhiGradient _vector_grad_phi_face;
2728 : VectorVariablePhiSecond _vector_second_phi_face;
2729 : VectorVariablePhiCurl _vector_curl_phi_face;
2730 : VectorVariablePhiDivergence _vector_div_phi_face;
2731 :
2732 : VectorVariablePhiValue _vector_phi_neighbor;
2733 : VectorVariablePhiGradient _vector_grad_phi_neighbor;
2734 : VectorVariablePhiSecond _vector_second_phi_neighbor;
2735 : VectorVariablePhiCurl _vector_curl_phi_neighbor;
2736 : VectorVariablePhiDivergence _vector_div_phi_neighbor;
2737 :
2738 : VectorVariablePhiValue _vector_phi_face_neighbor;
2739 : VectorVariablePhiGradient _vector_grad_phi_face_neighbor;
2740 : VectorVariablePhiSecond _vector_second_phi_face_neighbor;
2741 : VectorVariablePhiCurl _vector_curl_phi_face_neighbor;
2742 : VectorVariablePhiDivergence _vector_div_phi_face_neighbor;
2743 :
2744 : class FEShapeData
2745 : {
2746 : public:
2747 : VariablePhiValue _phi;
2748 : VariablePhiGradient _grad_phi;
2749 : VariablePhiSecond _second_phi;
2750 : VariablePhiCurl _curl_phi;
2751 : VariablePhiDivergence _div_phi;
2752 : };
2753 :
2754 : class VectorFEShapeData
2755 : {
2756 : public:
2757 : VectorVariablePhiValue _phi;
2758 : VectorVariablePhiGradient _grad_phi;
2759 : VectorVariablePhiSecond _second_phi;
2760 : VectorVariablePhiCurl _curl_phi;
2761 : VectorVariablePhiDivergence _div_phi;
2762 : };
2763 :
2764 : /// Shape function values, gradients, second derivatives for each FE type
2765 : mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data;
2766 : mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_face;
2767 : mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_neighbor;
2768 : mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_face_neighbor;
2769 : mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_lower;
2770 : mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_dual_lower;
2771 :
2772 : /// Shape function values, gradients, second derivatives for each vector FE type
2773 : mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data;
2774 : mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_face;
2775 : mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_neighbor;
2776 : mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_face_neighbor;
2777 : mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_lower;
2778 : mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_dual_lower;
2779 :
2780 : mutable std::map<FEType, ADTemplateVariablePhiGradient<Real>> _ad_grad_phi_data;
2781 : mutable std::map<FEType, ADTemplateVariablePhiGradient<RealVectorValue>> _ad_vector_grad_phi_data;
2782 : mutable std::map<FEType, ADTemplateVariablePhiGradient<Real>> _ad_grad_phi_data_face;
2783 : mutable std::map<FEType, ADTemplateVariablePhiGradient<RealVectorValue>>
2784 : _ad_vector_grad_phi_data_face;
2785 :
2786 : /**
2787 : * The residual vector tags that Assembly could possibly contribute to.
2788 : *
2789 : * The following variables are all indexed with this vector (i.e., index 0 in the following
2790 : * vectors corresponds to the tag with TagID _residual_vector_tags[0]._id):
2791 : * _sub_Re, _sub_Rn, _sub_Rl, _cached_residual_rows, _cached_residual_values,
2792 : *
2793 : * This index is also available in VectorTag::_type_id
2794 : */
2795 : const std::vector<VectorTag> & _residual_vector_tags;
2796 :
2797 : /// Values cached by calling cacheResidual() (the first vector is for TIME vs NONTIME)
2798 : std::vector<std::vector<Real>> _cached_residual_values;
2799 :
2800 : /// Where the cached values should go (the first vector is for TIME vs NONTIME)
2801 : std::vector<std::vector<dof_id_type>> _cached_residual_rows;
2802 :
2803 : unsigned int _max_cached_residuals;
2804 :
2805 : /// Values cached by calling cacheJacobian()
2806 : std::vector<std::vector<Real>> _cached_jacobian_values;
2807 : /// Row where the corresponding cached value should go
2808 : std::vector<std::vector<dof_id_type>> _cached_jacobian_rows;
2809 : /// Column where the corresponding cached value should go
2810 : std::vector<std::vector<dof_id_type>> _cached_jacobian_cols;
2811 :
2812 : unsigned int _max_cached_jacobians;
2813 :
2814 : /// Will be true if our preconditioning matrix is a block-diagonal matrix. Which means that we can take some shortcuts.
2815 : bool _block_diagonal_matrix;
2816 : /// An flag array Indiced by variable index to show if there is no component-wise
2817 : /// coupling for the variable.
2818 : std::vector<bool> _component_block_diagonal;
2819 :
2820 : /// Temporary work vector to keep from reallocating it
2821 : std::vector<dof_id_type> _temp_dof_indices;
2822 :
2823 : /// Temporary work data for reinitAtPhysical()
2824 : std::vector<Point> _temp_reference_points;
2825 :
2826 : /// AD quantities
2827 : std::vector<VectorValue<ADReal>> _ad_dxyzdxi_map;
2828 : std::vector<VectorValue<ADReal>> _ad_dxyzdeta_map;
2829 : std::vector<VectorValue<ADReal>> _ad_dxyzdzeta_map;
2830 : std::vector<VectorValue<ADReal>> _ad_d2xyzdxi2_map;
2831 : std::vector<VectorValue<ADReal>> _ad_d2xyzdxideta_map;
2832 : std::vector<VectorValue<ADReal>> _ad_d2xyzdeta2_map;
2833 : std::vector<ADReal> _ad_jac;
2834 : MooseArray<ADReal> _ad_JxW;
2835 : MooseArray<VectorValue<ADReal>> _ad_q_points;
2836 : std::vector<ADReal> _ad_dxidx_map;
2837 : std::vector<ADReal> _ad_dxidy_map;
2838 : std::vector<ADReal> _ad_dxidz_map;
2839 : std::vector<ADReal> _ad_detadx_map;
2840 : std::vector<ADReal> _ad_detady_map;
2841 : std::vector<ADReal> _ad_detadz_map;
2842 : std::vector<ADReal> _ad_dzetadx_map;
2843 : std::vector<ADReal> _ad_dzetady_map;
2844 : std::vector<ADReal> _ad_dzetadz_map;
2845 :
2846 : MooseArray<ADReal> _ad_JxW_face;
2847 : MooseArray<VectorValue<ADReal>> _ad_normals;
2848 : MooseArray<VectorValue<ADReal>> _ad_q_points_face;
2849 : MooseArray<Real> _curvatures;
2850 : MooseArray<ADReal> _ad_curvatures;
2851 :
2852 : /**
2853 : * Container of displacement numbers and directions
2854 : */
2855 : std::vector<std::pair<unsigned int, unsigned short>> _disp_numbers_and_directions;
2856 :
2857 : mutable bool _calculate_xyz;
2858 : mutable bool _calculate_face_xyz;
2859 : mutable bool _calculate_curvatures;
2860 :
2861 : /// Whether to calculate coord with AD. This will only be set to \p true if a consumer calls
2862 : /// adCoordTransformation()
2863 : mutable bool _calculate_ad_coord;
2864 :
2865 : mutable std::set<FEType> _need_second_derivative;
2866 : mutable std::set<FEType> _need_second_derivative_neighbor;
2867 : mutable std::set<FEType> _need_curl;
2868 : mutable std::set<FEType> _need_div;
2869 : mutable std::set<FEType> _need_face_div;
2870 : mutable std::set<FEType> _need_neighbor_div;
2871 : mutable std::set<FEType> _need_face_neighbor_div;
2872 :
2873 : /// The map from global index to variable scaling factor
2874 : const NumericVector<Real> * _scaling_vector = nullptr;
2875 :
2876 : /// In place side element builder for _current_side_elem
2877 : libMesh::ElemSideBuilder _current_side_elem_builder;
2878 : /// In place side element builder for _current_neighbor_side_elem
2879 : libMesh::ElemSideBuilder _current_neighbor_side_elem_builder;
2880 : /// In place side element builder for computeFaceMap()
2881 : libMesh::ElemSideBuilder _compute_face_map_side_elem_builder;
2882 :
2883 : const Elem * _msm_elem = nullptr;
2884 :
2885 : /// A working vector to avoid repeated heap allocations when caching residuals that must have
2886 : /// libMesh-level constraints (hanging nodes, periodic bcs) applied to them. This stores local
2887 : /// residual values
2888 : DenseVector<Number> _element_vector;
2889 :
2890 : /// A working matrix to avoid repeated heap allocations when caching Jacobians that must have
2891 : /// libMesh-level constraints (hanging nodes, periodic bcs) applied to them. This stores local
2892 : /// Jacobian values
2893 : DenseMatrix<Number> _element_matrix;
2894 :
2895 : /// Working vectors to avoid repeated heap allocations when caching residuals/Jacobians that must
2896 : /// have libMesh-level constraints (hanging nodes, periodic bcs) applied to them. These are for
2897 : /// storing the dof indices
2898 : std::vector<dof_id_type> _row_indices, _column_indices;
2899 :
2900 : /// Whether we have ever conducted p-refinement
2901 : bool _have_p_refinement;
2902 :
2903 : /// The current reference points on the neighbor element
2904 : std::vector<Point> _current_neighbor_ref_points;
2905 : };
2906 :
2907 : template <typename OutputType>
2908 : const typename OutputTools<OutputType>::VariablePhiValue &
2909 344764 : Assembly::fePhiLower(FEType type) const
2910 : {
2911 344764 : buildLowerDFE(type);
2912 344764 : return _fe_shape_data_lower[type]->_phi;
2913 : }
2914 :
2915 : template <typename OutputType>
2916 : const typename OutputTools<OutputType>::VariablePhiValue &
2917 270 : Assembly::feDualPhiLower(FEType type) const
2918 : {
2919 270 : buildLowerDDualFE(type);
2920 270 : return _fe_shape_data_dual_lower[type]->_phi;
2921 : }
2922 :
2923 : template <typename OutputType>
2924 : const typename OutputTools<OutputType>::VariablePhiGradient &
2925 344764 : Assembly::feGradPhiLower(FEType type) const
2926 : {
2927 344764 : buildLowerDFE(type);
2928 344764 : return _fe_shape_data_lower[type]->_grad_phi;
2929 : }
2930 :
2931 : template <typename OutputType>
2932 : const typename OutputTools<OutputType>::VariablePhiGradient &
2933 270 : Assembly::feGradDualPhiLower(FEType type) const
2934 : {
2935 270 : buildLowerDDualFE(type);
2936 270 : return _fe_shape_data_dual_lower[type]->_grad_phi;
2937 : }
2938 :
2939 : template <>
2940 : inline const ADTemplateVariablePhiGradient<RealVectorValue> &
2941 1554 : Assembly::feADGradPhi<RealVectorValue>(FEType type) const
2942 : {
2943 1554 : return _ad_vector_grad_phi_data[type];
2944 : }
2945 :
2946 : template <>
2947 : inline const ADTemplateVariablePhiGradient<RealVectorValue> &
2948 1554 : Assembly::feADGradPhiFace<RealVectorValue>(FEType type) const
2949 : {
2950 1554 : return _ad_vector_grad_phi_data_face[type];
2951 : }
2952 :
2953 : template <>
2954 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
2955 : Assembly::fePhi<VectorValue<Real>>(FEType type) const;
2956 :
2957 : template <>
2958 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
2959 : Assembly::feGradPhi<VectorValue<Real>>(FEType type) const;
2960 :
2961 : template <>
2962 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
2963 : Assembly::feSecondPhi<VectorValue<Real>>(FEType type) const;
2964 :
2965 : template <>
2966 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
2967 : Assembly::fePhiLower<VectorValue<Real>>(FEType type) const;
2968 :
2969 : template <>
2970 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
2971 : Assembly::feDualPhiLower<VectorValue<Real>>(FEType type) const;
2972 :
2973 : template <>
2974 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
2975 : Assembly::feGradPhiLower<VectorValue<Real>>(FEType type) const;
2976 :
2977 : template <>
2978 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
2979 : Assembly::feGradDualPhiLower<VectorValue<Real>>(FEType type) const;
2980 :
2981 : template <>
2982 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
2983 : Assembly::fePhiFace<VectorValue<Real>>(FEType type) const;
2984 :
2985 : template <>
2986 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
2987 : Assembly::feGradPhiFace<VectorValue<Real>>(FEType type) const;
2988 :
2989 : template <>
2990 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
2991 : Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type) const;
2992 :
2993 : template <>
2994 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
2995 : Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type) const;
2996 :
2997 : template <>
2998 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
2999 : Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type) const;
3000 :
3001 : template <>
3002 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
3003 : Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type) const;
3004 :
3005 : template <>
3006 : const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
3007 : Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
3008 :
3009 : template <>
3010 : const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
3011 : Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
3012 :
3013 : template <>
3014 : const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
3015 : Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
3016 :
3017 : template <>
3018 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
3019 : Assembly::feCurlPhi<VectorValue<Real>>(FEType type) const;
3020 :
3021 : template <>
3022 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
3023 : Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type) const;
3024 :
3025 : template <>
3026 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
3027 : Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type) const;
3028 :
3029 : template <>
3030 : const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
3031 : Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
3032 :
3033 : template <>
3034 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
3035 : Assembly::feDivPhi<VectorValue<Real>>(FEType type) const;
3036 :
3037 : template <>
3038 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
3039 : Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const;
3040 :
3041 : template <>
3042 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
3043 : Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const;
3044 :
3045 : template <>
3046 : const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
3047 : Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;
3048 :
3049 : template <>
3050 : inline const ADTemplateVariablePhiGradient<RealVectorValue> &
3051 200 : Assembly::adGradPhi<RealVectorValue>(const MooseVariableFE<RealVectorValue> & v) const
3052 : {
3053 200 : return _ad_vector_grad_phi_data.at(v.feType());
3054 : }
3055 :
3056 : template <typename Residuals, typename Indices>
3057 : void
3058 27185210 : Assembly::cacheResiduals(const Residuals & residuals,
3059 : const Indices & input_row_indices,
3060 : const Real scaling_factor,
3061 : LocalDataKey,
3062 : const std::set<TagID> & vector_tags)
3063 : {
3064 : mooseAssert(residuals.size() == input_row_indices.size(),
3065 : "The number of residuals should match the number of dof indices");
3066 : mooseAssert(residuals.size() >= 1, "Why you calling me with no residuals?");
3067 :
3068 27185210 : if (!computingResidual() || vector_tags.empty())
3069 8857061 : return;
3070 :
3071 18328149 : if (residuals.size() == 1)
3072 : {
3073 : // No constraining is required. (This is likely a finite volume computation if we only have a
3074 : // single dof)
3075 2913738 : cacheResidualsWithoutConstraints(
3076 2913738 : residuals, input_row_indices, scaling_factor, LocalDataKey{}, vector_tags);
3077 2913738 : return;
3078 : }
3079 :
3080 : // Need to make a copy because we might modify this in constrain_element_vector
3081 15414411 : _row_indices.assign(input_row_indices.begin(), input_row_indices.end());
3082 :
3083 15414411 : _element_vector.resize(_row_indices.size());
3084 82092521 : for (const auto i : index_range(_row_indices))
3085 66678110 : _element_vector(i) = MetaPhysicL::raw_value(residuals[i]) * scaling_factor;
3086 :
3087 : // At time of writing, this method doesn't do anything with the asymmetric_constraint_rows
3088 : // argument, but we set it to false to be consistent with processLocalResidual
3089 15414411 : _dof_map.constrain_element_vector(
3090 15414411 : _element_vector, _row_indices, /*asymmetric_constraint_rows=*/false);
3091 :
3092 82094813 : for (const auto i : index_range(_row_indices))
3093 66680402 : cacheResidual(_row_indices[i], _element_vector(i), vector_tags);
3094 : }
3095 :
3096 : template <typename Residuals, typename Indices>
3097 : void
3098 3073900 : Assembly::cacheResidualsWithoutConstraints(const Residuals & residuals,
3099 : const Indices & row_indices,
3100 : const Real scaling_factor,
3101 : LocalDataKey,
3102 : const std::set<TagID> & vector_tags)
3103 : {
3104 : mooseAssert(residuals.size() == row_indices.size(),
3105 : "The number of residuals should match the number of dof indices");
3106 : mooseAssert(residuals.size() >= 1, "Why you calling me with no residuals?");
3107 :
3108 3073900 : if (computingResidual() && !vector_tags.empty())
3109 5861652 : for (const auto i : index_range(row_indices))
3110 2939370 : cacheResidual(
3111 2939370 : row_indices[i], MetaPhysicL::raw_value(residuals[i]) * scaling_factor, vector_tags);
3112 3073900 : }
3113 :
3114 : template <typename Residuals, typename Indices>
3115 : void
3116 14365542 : Assembly::cacheJacobian(const Residuals & residuals,
3117 : const Indices & input_row_indices,
3118 : const Real scaling_factor,
3119 : LocalDataKey,
3120 : const std::set<TagID> & matrix_tags)
3121 : {
3122 14365542 : if (!computingJacobian() || matrix_tags.empty())
3123 0 : return;
3124 :
3125 14365542 : if (residuals.size() == 1)
3126 : {
3127 : // No constraining is required. (This is likely a finite volume computation if we only have a
3128 : // single dof)
3129 11905342 : cacheJacobianWithoutConstraints(
3130 11905342 : residuals, input_row_indices, scaling_factor, LocalDataKey{}, matrix_tags);
3131 11905342 : return;
3132 : }
3133 :
3134 2460200 : const auto & compare_dofs = residuals[0].derivatives().nude_indices();
3135 : #ifndef NDEBUG
3136 : auto compare_dofs_set = std::set<dof_id_type>(compare_dofs.begin(), compare_dofs.end());
3137 :
3138 : for (const auto i : make_range(decltype(residuals.size())(1), residuals.size()))
3139 : {
3140 : const auto & residual = residuals[i];
3141 : auto current_dofs_set = std::set<dof_id_type>(residual.derivatives().nude_indices().begin(),
3142 : residual.derivatives().nude_indices().end());
3143 : mooseAssert(compare_dofs_set == current_dofs_set,
3144 : "We're going to see whether the dof sets are the same. IIRC the degree of freedom "
3145 : "dependence (as indicated by the dof index set held by the ADReal) has to be the "
3146 : "same for every residual passed to this method otherwise constrain_element_matrix "
3147 : "will not work.");
3148 : }
3149 : #endif
3150 2460200 : _column_indices.assign(compare_dofs.begin(), compare_dofs.end());
3151 :
3152 : // If there's no derivatives then there is nothing to do. Moreover, if we pass zero size column
3153 : // indices to constrain_element_matrix then we will potentially get errors out of BLAS
3154 2460200 : if (!_column_indices.size())
3155 194421 : return;
3156 :
3157 : // Need to make a copy because we might modify this in constrain_element_matrix
3158 2265779 : _row_indices.assign(input_row_indices.begin(), input_row_indices.end());
3159 :
3160 2265779 : _element_matrix.resize(_row_indices.size(), _column_indices.size());
3161 11370329 : for (const auto i : index_range(_row_indices))
3162 : {
3163 9104550 : const auto & sparse_derivatives = residuals[i].derivatives();
3164 :
3165 64857340 : for (const auto j : index_range(_column_indices))
3166 55752790 : _element_matrix(i, j) = sparse_derivatives[_column_indices[j]] * scaling_factor;
3167 : }
3168 :
3169 2265779 : _dof_map.constrain_element_matrix(_element_matrix, _row_indices, _column_indices);
3170 :
3171 11370817 : for (const auto i : index_range(_row_indices))
3172 64867372 : for (const auto j : index_range(_column_indices))
3173 55762334 : cacheJacobian(_row_indices[i], _column_indices[j], _element_matrix(i, j), {}, matrix_tags);
3174 : }
3175 :
3176 : template <typename Residuals, typename Indices>
3177 : void
3178 12065312 : Assembly::cacheJacobianWithoutConstraints(const Residuals & residuals,
3179 : const Indices & row_indices,
3180 : const Real scaling_factor,
3181 : LocalDataKey,
3182 : const std::set<TagID> & matrix_tags)
3183 : {
3184 : mooseAssert(residuals.size() == row_indices.size(),
3185 : "The number of residuals should match the number of dof indices");
3186 : mooseAssert(residuals.size() >= 1, "Why you calling me with no residuals?");
3187 :
3188 12065312 : if (!computingJacobian() || matrix_tags.empty())
3189 0 : return;
3190 :
3191 25697578 : for (const auto i : index_range(row_indices))
3192 : {
3193 13632266 : const auto row_index = row_indices[i];
3194 :
3195 13632266 : const auto & sparse_derivatives = residuals[i].derivatives();
3196 13632266 : const auto & column_indices = sparse_derivatives.nude_indices();
3197 13632266 : const auto & raw_derivatives = sparse_derivatives.nude_data();
3198 :
3199 62922106 : for (std::size_t j = 0; j < column_indices.size(); ++j)
3200 98579680 : cacheJacobian(
3201 49289840 : row_index, column_indices[j], raw_derivatives[j] * scaling_factor, {}, matrix_tags);
3202 : }
3203 : }
3204 :
3205 : inline const Real &
3206 434 : Assembly::lowerDElemVolume() const
3207 : {
3208 434 : _need_lower_d_elem_volume = true;
3209 434 : return _current_lower_d_elem_volume;
3210 : }
3211 :
3212 : inline const Real &
3213 408 : Assembly::neighborLowerDElemVolume() const
3214 : {
3215 408 : _need_neighbor_lower_d_elem_volume = true;
3216 408 : return _current_neighbor_lower_d_elem_volume;
3217 : }
3218 :
3219 : inline void
3220 2221 : Assembly::assignDisplacements(
3221 : std::vector<std::pair<unsigned int, unsigned short>> && disp_numbers_and_directions)
3222 : {
3223 2221 : _disp_numbers_and_directions = std::move(disp_numbers_and_directions);
3224 2221 : }
3225 :
3226 : inline void
3227 139453 : Assembly::setCurrentLowerDElem(const Elem * const lower_d_elem)
3228 : {
3229 139453 : _current_lower_d_elem = lower_d_elem;
3230 139453 : }
|
