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