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