Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "MooseVariableData.h"
11 : #include "MooseVariableField.h"
12 : #include "Assembly.h"
13 : #include "MooseError.h"
14 : #include "DisplacedSystem.h"
15 : #include "TimeIntegrator.h"
16 : #include "MooseVariableFE.h"
17 : #include "MooseTypes.h"
18 :
19 : #include "libmesh/quadrature.h"
20 : #include "libmesh/fe_base.h"
21 : #include "libmesh/system.h"
22 : #include "libmesh/type_n_tensor.h"
23 : #include "libmesh/fe_interface.h"
24 :
25 : using namespace libMesh;
26 :
27 : template <typename OutputType>
28 530991 : MooseVariableData<OutputType>::MooseVariableData(const MooseVariableField<OutputType> & var,
29 : SystemBase & sys,
30 : THREAD_ID tid,
31 : Moose::ElementType element_type,
32 : const QBase * const & qrule_in,
33 : const QBase * const & qrule_face_in,
34 : const Node * const & node,
35 : const Elem * const & elem)
36 :
37 : : MooseVariableDataBase<OutputType>(var, sys, tid),
38 1061982 : _fe_type(var.feType()),
39 530991 : _var_num(var.number()),
40 530991 : _assembly(_subproblem.assembly(_tid, var.kind() == Moose::VAR_SOLVER ? sys.number() : 0)),
41 530991 : _element_type(element_type),
42 530991 : _ad_zero(0),
43 530991 : _need_ad_u_dot(false),
44 530991 : _need_ad_u_dotdot(false),
45 530991 : _need_second(false),
46 530991 : _need_second_old(false),
47 530991 : _need_second_older(false),
48 530991 : _need_second_previous_nl(false),
49 530991 : _need_curl(false),
50 530991 : _need_curl_old(false),
51 530991 : _need_curl_older(false),
52 530991 : _need_div(false),
53 530991 : _need_div_old(false),
54 530991 : _need_div_older(false),
55 530991 : _need_ad(false),
56 530991 : _need_ad_u(false),
57 530991 : _need_ad_grad_u(false),
58 530991 : _need_ad_grad_u_dot(false),
59 530991 : _need_ad_second_u(false),
60 530991 : _need_ad_curl_u(false),
61 530991 : _has_dof_indices(false),
62 530991 : _qrule(qrule_in),
63 530991 : _qrule_face(qrule_face_in),
64 530991 : _use_dual(var.useDual()),
65 530991 : _second_phi_assembly_method(nullptr),
66 530991 : _second_phi_face_assembly_method(nullptr),
67 530991 : _curl_phi_assembly_method(nullptr),
68 530991 : _curl_phi_face_assembly_method(nullptr),
69 530991 : _div_phi_assembly_method(nullptr),
70 530991 : _div_phi_face_assembly_method(nullptr),
71 530991 : _ad_grad_phi_assembly_method(nullptr),
72 530991 : _ad_grad_phi_face_assembly_method(nullptr),
73 530991 : _time_integrator(nullptr),
74 530991 : _node(node),
75 530991 : _elem(elem),
76 530991 : _displaced(dynamic_cast<const DisplacedSystem *>(&_sys) ? true : false),
77 2654955 : _current_side(_assembly.side())
78 : {
79 530991 : _continuity = FEInterface::get_continuity(_fe_type);
80 :
81 530991 : _is_nodal = (_continuity == C_ZERO || _continuity == C_ONE);
82 :
83 530991 : _time_integrator = _sys.queryTimeIntegrator(_var_num);
84 :
85 : // Initialize AD zero with zero derivatives
86 530991 : const auto old_do = ADReal::do_derivatives;
87 530991 : ADReal::do_derivatives = true;
88 530991 : _ad_zero = 0.;
89 530991 : ADReal::do_derivatives = old_do;
90 :
91 530991 : switch (_element_type)
92 : {
93 176997 : case Moose::ElementType::Element:
94 : {
95 176997 : _phi_assembly_method = &Assembly::fePhi<OutputShape>;
96 176997 : _phi_face_assembly_method = &Assembly::fePhiFace<OutputShape>;
97 176997 : _grad_phi_assembly_method = &Assembly::feGradPhi<OutputShape>;
98 176997 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFace<OutputShape>;
99 176997 : _second_phi_assembly_method = &Assembly::feSecondPhi<OutputShape>;
100 176997 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFace<OutputShape>;
101 176997 : _curl_phi_assembly_method = &Assembly::feCurlPhi<OutputShape>;
102 176997 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFace<OutputShape>;
103 176997 : _div_phi_assembly_method = &Assembly::feDivPhi<OutputShape>;
104 176997 : _div_phi_face_assembly_method = &Assembly::feDivPhiFace<OutputShape>;
105 176997 : _ad_grad_phi_assembly_method = &Assembly::feADGradPhi<OutputShape>;
106 176997 : _ad_grad_phi_face_assembly_method = &Assembly::feADGradPhiFace<OutputShape>;
107 :
108 176997 : _ad_grad_phi = &_ad_grad_phi_assembly_method(_assembly, _fe_type);
109 176997 : _ad_grad_phi_face = &_ad_grad_phi_face_assembly_method(_assembly, _fe_type);
110 176997 : break;
111 : }
112 176997 : case Moose::ElementType::Neighbor:
113 : {
114 176997 : _phi_assembly_method = &Assembly::fePhiNeighbor<OutputShape>;
115 176997 : _phi_face_assembly_method = &Assembly::fePhiFaceNeighbor<OutputShape>;
116 176997 : _grad_phi_assembly_method = &Assembly::feGradPhiNeighbor<OutputShape>;
117 176997 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFaceNeighbor<OutputShape>;
118 176997 : _second_phi_assembly_method = &Assembly::feSecondPhiNeighbor<OutputShape>;
119 176997 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFaceNeighbor<OutputShape>;
120 176997 : _curl_phi_assembly_method = &Assembly::feCurlPhiNeighbor<OutputShape>;
121 176997 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFaceNeighbor<OutputShape>;
122 176997 : _div_phi_assembly_method = &Assembly::feDivPhiNeighbor<OutputShape>;
123 176997 : _div_phi_face_assembly_method = &Assembly::feDivPhiFaceNeighbor<OutputShape>;
124 :
125 176997 : _ad_grad_phi = nullptr;
126 176997 : _ad_grad_phi_face = nullptr;
127 176997 : break;
128 : }
129 176997 : case Moose::ElementType::Lower:
130 : {
131 176997 : if (_use_dual)
132 : {
133 146 : _phi_assembly_method = &Assembly::feDualPhiLower<OutputType>;
134 146 : _phi_face_assembly_method = &Assembly::feDualPhiLower<OutputType>; // Place holder
135 146 : _grad_phi_assembly_method = &Assembly::feGradDualPhiLower<OutputType>;
136 146 : _grad_phi_face_assembly_method = &Assembly::feGradDualPhiLower<OutputType>; // Place holder
137 : }
138 : else
139 : {
140 176851 : _phi_assembly_method = &Assembly::fePhiLower<OutputType>;
141 176851 : _phi_face_assembly_method = &Assembly::fePhiLower<OutputType>; // Place holder
142 176851 : _grad_phi_assembly_method = &Assembly::feGradPhiLower<OutputType>;
143 176851 : _grad_phi_face_assembly_method = &Assembly::feGradPhiLower<OutputType>; // Place holder
144 : }
145 :
146 176997 : _ad_grad_phi = nullptr;
147 176997 : _ad_grad_phi_face = nullptr;
148 176997 : break;
149 : }
150 : }
151 530991 : _phi = &_phi_assembly_method(_assembly, _fe_type);
152 530991 : _phi_face = &_phi_face_assembly_method(_assembly, _fe_type);
153 530991 : _grad_phi = &_grad_phi_assembly_method(_assembly, _fe_type);
154 530991 : _grad_phi_face = &_grad_phi_face_assembly_method(_assembly, _fe_type);
155 530991 : }
156 :
157 : template <typename OutputType>
158 : void
159 787811377 : MooseVariableData<OutputType>::setGeometry(Moose::GeometryType gm_type)
160 : {
161 787811377 : switch (gm_type)
162 : {
163 742546044 : case Moose::Volume:
164 : {
165 742546044 : _current_qrule = _qrule;
166 742546044 : _current_phi = _phi;
167 742546044 : _current_grad_phi = _grad_phi;
168 742546044 : _current_second_phi = _second_phi;
169 742546044 : _current_curl_phi = _curl_phi;
170 742546044 : _current_div_phi = _div_phi;
171 742546044 : _current_ad_grad_phi = _ad_grad_phi;
172 742546044 : break;
173 : }
174 45265333 : case Moose::Face:
175 : {
176 45265333 : _current_qrule = _qrule_face;
177 45265333 : _current_phi = _phi_face;
178 45265333 : _current_grad_phi = _grad_phi_face;
179 45265333 : _current_second_phi = _second_phi_face;
180 45265333 : _current_curl_phi = _curl_phi_face;
181 45265333 : _current_div_phi = _div_phi_face;
182 45265333 : _current_ad_grad_phi = _ad_grad_phi_face;
183 45265333 : break;
184 : }
185 : }
186 787811377 : }
187 :
188 : template <typename OutputType>
189 : const typename MooseVariableData<OutputType>::FieldVariableValue &
190 17145 : MooseVariableData<OutputType>::uDot() const
191 : {
192 17145 : if (_sys.solutionUDot())
193 : {
194 17145 : _need_u_dot = true;
195 17145 : return _u_dot;
196 : }
197 : else
198 0 : mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set "
199 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
200 : }
201 :
202 : template <typename OutputType>
203 : const typename MooseVariableData<OutputType>::FieldVariableValue &
204 219 : MooseVariableData<OutputType>::uDotDot() const
205 : {
206 219 : if (_sys.solutionUDotDot())
207 : {
208 219 : _need_u_dotdot = true;
209 219 : return _u_dotdot;
210 : }
211 : else
212 0 : mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. "
213 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
214 : "`u_dotdot`.");
215 : }
216 :
217 : template <typename OutputType>
218 : const typename MooseVariableData<OutputType>::FieldVariableValue &
219 0 : MooseVariableData<OutputType>::uDotOld() const
220 : {
221 0 : if (_sys.solutionUDotOld())
222 : {
223 0 : _need_u_dot_old = true;
224 0 : return _u_dot_old;
225 : }
226 : else
227 0 : mooseError("MooseVariableFE: Old time derivative of solution (`u_dot_old`) is not stored. "
228 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
229 : "`u_dot_old`.");
230 : }
231 :
232 : template <typename OutputType>
233 : const typename MooseVariableData<OutputType>::FieldVariableValue &
234 0 : MooseVariableData<OutputType>::uDotDotOld() const
235 : {
236 0 : if (_sys.solutionUDotDotOld())
237 : {
238 0 : _need_u_dotdot_old = true;
239 0 : return _u_dotdot_old;
240 : }
241 : else
242 0 : mooseError("MooseVariableFE: Old second time derivative of solution (`u_dotdot_old`) is not "
243 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
244 : "requesting `u_dotdot_old`");
245 : }
246 :
247 : template <typename OutputType>
248 : const typename MooseVariableData<OutputType>::FieldVariableGradient &
249 39 : MooseVariableData<OutputType>::gradSlnDot() const
250 : {
251 39 : if (_sys.solutionUDot())
252 : {
253 39 : _need_grad_dot = true;
254 39 : return _grad_u_dot;
255 : }
256 : else
257 0 : mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set "
258 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
259 : }
260 :
261 : template <typename OutputType>
262 : const typename MooseVariableData<OutputType>::FieldVariableGradient &
263 0 : MooseVariableData<OutputType>::gradSlnDotDot() const
264 : {
265 0 : if (_sys.solutionUDotDot())
266 : {
267 0 : _need_grad_dotdot = true;
268 0 : return _grad_u_dotdot;
269 : }
270 : else
271 0 : mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. "
272 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
273 : "`u_dotdot`.");
274 : }
275 :
276 : template <typename OutputType>
277 : const typename MooseVariableData<OutputType>::FieldVariableSecond &
278 201 : MooseVariableData<OutputType>::secondSln(Moose::SolutionState state) const
279 : {
280 201 : secondPhi();
281 201 : secondPhiFace();
282 201 : switch (state)
283 : {
284 201 : case Moose::Current:
285 : {
286 201 : _need_second = true;
287 201 : return _second_u;
288 : }
289 :
290 0 : case Moose::Old:
291 : {
292 0 : _need_second_old = true;
293 0 : return _second_u_old;
294 : }
295 :
296 0 : case Moose::Older:
297 : {
298 0 : _need_second_older = true;
299 0 : return _second_u_older;
300 : }
301 :
302 0 : case Moose::PreviousNL:
303 : {
304 0 : _need_second_previous_nl = true;
305 0 : return _second_u_previous_nl;
306 : }
307 :
308 0 : default:
309 : // We should never get here but gcc requires that we have a default. See
310 : // htpps://stackoverflow.com/questions/18680378/after-defining-case-for-all-enum-values-compiler-still-says-control-reaches-e
311 0 : mooseError("Unknown SolutionState!");
312 : }
313 : }
314 :
315 : template <typename OutputType>
316 : const typename MooseVariableData<OutputType>::FieldVariableCurl &
317 171 : MooseVariableData<OutputType>::curlSln(Moose::SolutionState state) const
318 : {
319 171 : curlPhi();
320 171 : curlPhiFace();
321 171 : switch (state)
322 : {
323 171 : case Moose::Current:
324 : {
325 171 : _need_curl = true;
326 171 : return _curl_u;
327 : }
328 :
329 0 : case Moose::Old:
330 : {
331 0 : _need_curl_old = true;
332 0 : return _curl_u_old;
333 : }
334 :
335 0 : case Moose::Older:
336 : {
337 0 : _need_curl_older = true;
338 0 : return _curl_u_older;
339 : }
340 :
341 0 : default:
342 0 : mooseError("We don't currently support curl from the previous non-linear iteration");
343 : }
344 : }
345 :
346 : template <typename OutputType>
347 : const typename MooseVariableData<OutputType>::FieldVariableDivergence &
348 336 : MooseVariableData<OutputType>::divSln(Moose::SolutionState state) const
349 : {
350 336 : divPhi();
351 336 : divPhiFace();
352 336 : switch (state)
353 : {
354 336 : case Moose::Current:
355 : {
356 336 : _need_div = true;
357 336 : return _div_u;
358 : }
359 :
360 0 : case Moose::Old:
361 : {
362 0 : _need_div_old = true;
363 0 : return _div_u_old;
364 : }
365 :
366 0 : case Moose::Older:
367 : {
368 0 : _need_div_older = true;
369 0 : return _div_u_older;
370 : }
371 :
372 0 : default:
373 0 : mooseError("We don't currently support divergence from the previous non-linear iteration");
374 : }
375 : }
376 :
377 : template <typename OutputType>
378 : const typename MooseVariableData<OutputType>::FieldVariablePhiSecond &
379 29260 : MooseVariableData<OutputType>::secondPhi() const
380 : {
381 29260 : _second_phi = &_second_phi_assembly_method(_assembly, _fe_type);
382 29260 : return *_second_phi;
383 : }
384 :
385 : template <typename OutputType>
386 : const typename MooseVariableData<OutputType>::FieldVariablePhiSecond &
387 6304 : MooseVariableData<OutputType>::secondPhiFace() const
388 : {
389 6304 : _second_phi_face = &_second_phi_face_assembly_method(_assembly, _fe_type);
390 6304 : return *_second_phi_face;
391 : }
392 :
393 : template <typename OutputType>
394 : const typename MooseVariableData<OutputType>::FieldVariablePhiCurl &
395 64382 : MooseVariableData<OutputType>::curlPhi() const
396 : {
397 64382 : _curl_phi = &_curl_phi_assembly_method(_assembly, _fe_type);
398 64382 : return *_curl_phi;
399 : }
400 :
401 : template <typename OutputType>
402 : const typename MooseVariableData<OutputType>::FieldVariablePhiCurl &
403 224 : MooseVariableData<OutputType>::curlPhiFace() const
404 : {
405 224 : _curl_phi_face = &_curl_phi_face_assembly_method(_assembly, _fe_type);
406 224 : return *_curl_phi_face;
407 : }
408 :
409 : template <typename OutputType>
410 : const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
411 527632 : MooseVariableData<OutputType>::divPhi() const
412 : {
413 527632 : _div_phi = &_div_phi_assembly_method(_assembly, _fe_type);
414 527632 : return *_div_phi;
415 : }
416 :
417 : template <typename OutputType>
418 : const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
419 336 : MooseVariableData<OutputType>::divPhiFace() const
420 : {
421 336 : _div_phi_face = &_div_phi_face_assembly_method(_assembly, _fe_type);
422 336 : return *_div_phi_face;
423 : }
424 :
425 : template <typename OutputType>
426 : template <bool monomial>
427 : void
428 787033688 : MooseVariableData<OutputType>::computeValuesInternal()
429 : {
430 787033688 : const auto num_dofs = _dof_indices.size();
431 787033688 : if (num_dofs > 0)
432 685711470 : fetchDoFValues();
433 :
434 787033688 : const bool is_transient = _subproblem.isTransient();
435 787033688 : const auto nqp = _current_qrule->n_points();
436 787033688 : auto && active_coupleable_matrix_tags = _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid);
437 :
438 : // Map grad_phi using Eigen so that we can perform array operations easier
439 : if constexpr (std::is_same_v<OutputType, RealEigenVector>)
440 : {
441 4580352 : if (_qrule == _current_qrule)
442 : {
443 4336258 : _mapped_grad_phi.resize(num_dofs);
444 20341835 : for (const auto i : make_range(num_dofs))
445 : {
446 16005577 : _mapped_grad_phi[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
447 106141802 : for (const auto qp : make_range(nqp))
448 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
449 180272450 : new (&_mapped_grad_phi[i][qp])
450 90136225 : Eigen::Map<RealDIMValue>(const_cast<Real *>(&(*_current_grad_phi)[i][qp](0)));
451 : }
452 : }
453 : else
454 : {
455 244094 : _mapped_grad_phi_face.resize(num_dofs);
456 973106 : for (const auto i : make_range(num_dofs))
457 : {
458 729012 : _mapped_grad_phi_face[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
459 2569476 : for (const auto qp : make_range(nqp))
460 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
461 3680928 : new (&_mapped_grad_phi_face[i][qp])
462 1840464 : Eigen::Map<RealDIMValue>(const_cast<Real *>(&(*_current_grad_phi)[i][qp](0)));
463 : }
464 : }
465 : }
466 :
467 : mooseAssert(
468 : !(_need_second || _need_second_old || _need_second_older || _need_second_previous_nl) ||
469 : _current_second_phi,
470 : "We're requiring a second calculation but have not set a second shape function!");
471 : mooseAssert(!(_need_curl || _need_curl_old) || _current_curl_phi,
472 : "We're requiring a curl calculation but have not set a curl shape function!");
473 : mooseAssert(!(_need_div || _need_div_old) || _current_div_phi,
474 : "We're requiring a divergence calculation but have not set a div shape function!");
475 :
476 : // Helper for filling values
477 2309745372 : const auto fill = [this, nqp, num_dofs](auto & dest, const auto & phi, const auto & dof_values)
478 : {
479 : if constexpr (monomial)
480 178776732 : libmesh_ignore(num_dofs);
481 :
482 : // Deduce OutputType
483 1522711684 : constexpr bool is_real = std::is_same_v<OutputType, Real>;
484 1522711684 : constexpr bool is_real_vector = std::is_same_v<OutputType, RealVectorValue>;
485 1522711684 : constexpr bool is_eigen = std::is_same_v<OutputType, RealEigenVector>;
486 : static_assert(is_real || is_real_vector || is_eigen, "Unsupported type");
487 :
488 : // this is only used in the RealEigenVector case to get this->_count
489 : if constexpr (!is_eigen)
490 1513377002 : libmesh_ignore(this);
491 :
492 : // Deduce type of value within dest MooseArray
493 : using dest_array_type = typename std::remove_reference_t<decltype(dest)>::value_type;
494 1522711684 : constexpr bool is_value = std::is_same_v<dest_array_type, OutputType>;
495 1522711684 : constexpr bool is_gradient = std::is_same_v<dest_array_type, OutputGradient>;
496 1522711684 : constexpr bool is_second = std::is_same_v<dest_array_type, OutputSecond>;
497 1522711684 : constexpr bool is_divergence = std::is_same_v<dest_array_type, OutputDivergence>;
498 : static_assert(is_value || is_gradient || is_second || is_divergence,
499 : "Unsupported destination array type");
500 :
501 : // Sets a value to zero at a quadrature point
502 20210914264 : const auto set_zero = [this, &dest](const auto qp)
503 : {
504 : if constexpr (!is_eigen)
505 6186007530 : libmesh_ignore(this);
506 :
507 : if constexpr (is_real || is_real_vector)
508 6186007530 : dest[qp] = 0;
509 : else if constexpr (is_eigen)
510 : {
511 : if constexpr (is_value)
512 30379213 : dest[qp].setZero(this->_count);
513 : else if constexpr (is_gradient)
514 13014117 : dest[qp].setZero(this->_count, LIBMESH_DIM);
515 : else if constexpr (is_second)
516 0 : dest[qp].setZero(this->_count, LIBMESH_DIM * LIBMESH_DIM);
517 : else
518 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
519 : }
520 : else
521 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
522 : };
523 :
524 : // Accumulates a value
525 >11385*10^7 : const auto accumulate = [&dest, &phi, &dof_values](const auto i, const auto qp)
526 : {
527 : if constexpr (is_real || is_real_vector || (is_eigen && is_value))
528 : {
529 : if constexpr (is_value || is_divergence)
530 17738237230 : dest[qp] += phi[i][qp] * dof_values[i];
531 : else if constexpr (is_gradient || is_second)
532 10327503635 : dest[qp].add_scaled(phi[i][qp], dof_values[i]);
533 : else
534 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
535 : }
536 : else if constexpr (is_eigen)
537 : {
538 : if constexpr (is_gradient)
539 : {
540 291557732 : for (const auto d : make_range(Moose::dim))
541 218668299 : dest[qp].col(d) += phi[i][qp](d) * dof_values[i];
542 : }
543 : else if constexpr (is_second)
544 : {
545 0 : for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1)
546 0 : for (const auto d2 : make_range(Moose::dim))
547 0 : dest[qp].col(d++) += phi[i][qp](d1, d2) * dof_values[i];
548 : }
549 : else
550 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
551 : }
552 : else
553 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
554 : };
555 :
556 1522711684 : dest.resize(nqp);
557 :
558 : // Monomial case, accumulate dest[0] and set dest[>0] to dest[0]
559 : if constexpr (monomial)
560 : {
561 : mooseAssert(num_dofs == 1, "Should have only one dof");
562 178776732 : set_zero(0);
563 178776732 : accumulate(0, 0);
564 807389550 : for (unsigned int qp = 1; qp < nqp; ++qp)
565 628612818 : dest[qp] = dest[0];
566 : }
567 : // Non monomial case
568 : else
569 : {
570 7394559080 : for (const auto qp : make_range(nqp))
571 6050624128 : set_zero(qp);
572 6981298999 : for (const auto i : make_range(num_dofs))
573 33597217613 : for (const auto qp : make_range(nqp))
574 27959853566 : accumulate(i, qp);
575 : }
576 : };
577 :
578 : // Curl
579 787033688 : if (_need_curl)
580 79901 : fill(_curl_u, *_current_curl_phi, _vector_tags_dof_u[_solution_tag]);
581 787033688 : if (_need_curl_old)
582 0 : fill(_curl_u_old, *_current_curl_phi, _vector_tags_dof_u[_old_solution_tag]);
583 :
584 : // Div
585 787033688 : if (_need_div)
586 396144 : fill(_div_u, *_current_div_phi, _vector_tags_dof_u[_solution_tag]);
587 787033688 : if (is_transient && _need_div_old)
588 0 : fill(_div_u_old, *_current_div_phi, _vector_tags_dof_u[_old_solution_tag]);
589 :
590 : // Second
591 787033688 : if (_need_second)
592 92492 : fill(_second_u, *_current_second_phi, _vector_tags_dof_u[_solution_tag]);
593 787033688 : if (_need_second_previous_nl)
594 0 : fill(
595 0 : _second_u_previous_nl, *_current_second_phi, _vector_tags_dof_u[_previous_nl_solution_tag]);
596 :
597 : // Vector tags
598 1584195720 : for (auto tag : _required_vector_tags)
599 : {
600 797162032 : if (_need_vector_tag_u[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed())
601 794387882 : fill(_vector_tag_u[tag], *_current_phi, _vector_tags_dof_u[tag]);
602 797162032 : if (_need_vector_tag_grad[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed())
603 453915862 : fill(_vector_tag_grad[tag], *_current_grad_phi, _vector_tags_dof_u[tag]);
604 : }
605 :
606 : // Matrix tags
607 787034120 : for (auto tag : active_coupleable_matrix_tags)
608 432 : if (_need_matrix_tag_u[tag])
609 144 : fill(_matrix_tag_u[tag], *_current_phi, _matrix_tags_dof_u[tag]);
610 :
611 : // Derivatives and old values
612 787033688 : if (is_transient)
613 : {
614 591968398 : if (_need_second_old)
615 0 : fill(_second_u_old, *_current_second_phi, _vector_tags_dof_u[_old_solution_tag]);
616 591968398 : if (_need_second_older)
617 0 : fill(_second_u_older, *_current_second_phi, _vector_tags_dof_u[_older_solution_tag]);
618 591968398 : if (_need_u_dot)
619 273690781 : fill(_u_dot, *_current_phi, _dof_values_dot);
620 591968398 : if (_need_u_dotdot)
621 7145 : fill(_u_dotdot, *_current_phi, _dof_values_dotdot);
622 591968398 : if (_need_u_dot_old)
623 0 : fill(_u_dot_old, *_current_phi, _dof_values_dot_old);
624 591968398 : if (_need_u_dotdot_old)
625 0 : fill(_u_dotdot_old, *_current_phi, _dof_values_dotdot_old);
626 :
627 591968398 : if (_need_du_dot_du)
628 : {
629 273828089 : _du_dot_du.resize(nqp);
630 1463298721 : for (const auto qp : make_range(nqp))
631 1189470632 : _du_dot_du[qp] = 0.;
632 1464815557 : for (const auto i : make_range(num_dofs))
633 6860169170 : for (const auto qp : make_range(nqp))
634 5669181702 : _du_dot_du[qp] = _dof_du_dot_du[i];
635 : }
636 591968398 : if (_need_du_dotdot_du)
637 : {
638 5822 : _du_dotdot_du.resize(nqp);
639 27772 : for (const auto qp : make_range(nqp))
640 21950 : _du_dotdot_du[qp] = 0.;
641 49276 : for (const auto i : make_range(num_dofs))
642 215932 : for (const auto qp : make_range(nqp))
643 172478 : _du_dotdot_du[qp] = _dof_du_dotdot_du[i];
644 : }
645 :
646 591968398 : if (_need_grad_dot)
647 141333 : fill(_grad_u_dot, *_current_grad_phi, _dof_values_dot);
648 591968398 : if (_need_grad_dotdot)
649 0 : fill(_grad_u_dotdot, *_current_grad_phi, _dof_values_dotdot);
650 : }
651 :
652 : // Automatic differentiation, not for eigen
653 : if constexpr (!std::is_same_v<OutputType, RealEigenVector>)
654 782453336 : if (_need_ad)
655 26243703 : computeAD(num_dofs, nqp);
656 787033688 : }
657 :
658 : template <typename OutputType>
659 : void
660 619501374 : MooseVariableData<OutputType>::computeValues()
661 : {
662 619501374 : computeValuesInternal</* monomial = */ false>();
663 619501374 : }
664 :
665 : template <typename OutputType>
666 : void
667 168391805 : MooseVariableData<OutputType>::computeMonomialValues()
668 : {
669 168391805 : if (_dof_indices.size() == 0)
670 777689 : return;
671 :
672 : // Monomial optimizations are not appropriate after p-refinement
673 167614116 : if (_elem->p_level())
674 81802 : computeValues();
675 : else
676 167532314 : computeValuesInternal</* monomial = */ true>();
677 : }
678 :
679 : template <typename OutputType>
680 : void
681 26243703 : MooseVariableData<OutputType>::computeAD(const unsigned int num_dofs, const unsigned int nqp)
682 : {
683 26243703 : const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys);
684 :
685 26243703 : _ad_dof_values.resize(num_dofs);
686 144511926 : for (const auto i : make_range(num_dofs))
687 118268223 : _ad_dof_values[i] = (*_sys.currentSolution())(_dof_indices[i]);
688 : // NOTE! You have to do this AFTER setting the value!
689 26243703 : if (do_derivatives)
690 30578185 : for (const auto i : make_range(num_dofs))
691 25459516 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
692 :
693 26243703 : if (_need_ad_u)
694 : {
695 24756613 : _ad_u.resize(nqp);
696 108243899 : for (const auto qp : make_range(nqp))
697 83487286 : _ad_u[qp] = _ad_zero;
698 :
699 138238926 : for (const auto i : make_range(num_dofs))
700 554258719 : for (const auto qp : make_range(nqp))
701 440776406 : _ad_u[qp] += _ad_dof_values[i] * (*_current_phi)[i][qp];
702 : }
703 :
704 26243703 : if (_need_ad_grad_u)
705 : {
706 19372059 : _ad_grad_u.resize(nqp);
707 87872160 : for (const auto qp : make_range(nqp))
708 68500101 : _ad_grad_u[qp] = _ad_zero;
709 :
710 : // The latter check here is for handling the fact that we have not yet implemented
711 : // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that
712 : // situation we need to default to using the non-ad grad_phi
713 19372059 : if (_displaced && _current_ad_grad_phi)
714 2998844 : for (const auto i : make_range(num_dofs))
715 16381446 : for (const auto qp : make_range(nqp))
716 14476760 : _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_ad_grad_phi)[i][qp];
717 : else
718 97753320 : for (const auto i : make_range(num_dofs))
719 417259985 : for (const auto qp : make_range(nqp))
720 338331645 : _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_grad_phi)[i][qp];
721 : }
722 :
723 26243703 : if (_need_ad_second_u)
724 : {
725 0 : _ad_second_u.resize(nqp);
726 0 : for (const auto qp : make_range(nqp))
727 0 : _ad_second_u[qp] = _ad_zero;
728 :
729 0 : for (const auto i : make_range(num_dofs))
730 0 : for (const auto qp : make_range(nqp))
731 : // Note that this will not carry any derivatives with respect to displacements because
732 : // those calculations have not yet been implemented in Assembly
733 0 : _ad_second_u[qp] += _ad_dof_values[i] * (*_current_second_phi)[i][qp];
734 : }
735 :
736 26243703 : if (_need_ad_curl_u)
737 : {
738 2860555 : _ad_curl_u.resize(nqp);
739 13377115 : for (const auto qp : make_range(nqp))
740 10516560 : _ad_curl_u[qp] = _ad_zero;
741 :
742 14380175 : for (const auto i : make_range(num_dofs))
743 54075460 : for (const auto qp : make_range(nqp))
744 : {
745 : mooseAssert(_current_curl_phi,
746 : "We're requiring a curl calculation but have not set a curl shape function!");
747 :
748 : // Note that the current version of _ad_curl_u is not yet implemented for mesh displacement
749 42555840 : _ad_curl_u[qp] += _ad_dof_values[i] * (*_current_curl_phi)[i][qp];
750 : }
751 : }
752 :
753 26243703 : bool is_transient = _subproblem.isTransient();
754 26243703 : if (is_transient)
755 : {
756 10929951 : if (_need_ad_u_dot)
757 : {
758 2490231 : _ad_dofs_dot.resize(num_dofs);
759 2490231 : if (_need_ad_u_dotdot)
760 449 : _ad_dofs_dotdot.resize(num_dofs);
761 2490231 : _ad_u_dot.resize(nqp);
762 10805745 : for (const auto qp : make_range(nqp))
763 8315514 : _ad_u_dot[qp] = _ad_zero;
764 :
765 2490231 : if (_time_integrator && _time_integrator->dt())
766 : {
767 12686220 : for (const auto i : make_range(num_dofs))
768 10216650 : _ad_dofs_dot[i] = _ad_dof_values[i];
769 12686220 : for (const auto i : make_range(num_dofs))
770 20432014 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
771 10216650 : _dof_indices[i],
772 10216650 : _need_ad_u_dotdot ? _ad_dofs_dotdot[i]
773 : : _ad_real_dummy);
774 :
775 12686220 : for (const auto i : make_range(num_dofs))
776 44862272 : for (const auto qp : make_range(nqp))
777 34645622 : _ad_u_dot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dot[i];
778 : }
779 20661 : else if (!_time_integrator)
780 : {
781 0 : for (const auto i : make_range(num_dofs))
782 0 : _ad_dofs_dot[i] = _dof_values_dot[i];
783 0 : for (const auto qp : make_range(nqp))
784 0 : _ad_u_dot[qp] = _u_dot[qp];
785 : }
786 : }
787 :
788 10929951 : if (_need_ad_u_dotdot)
789 : {
790 449 : _ad_u_dotdot.resize(nqp);
791 1807 : for (const auto qp : make_range(nqp))
792 1358 : _ad_u_dotdot[qp] = _ad_zero;
793 :
794 449 : if (_time_integrator && _time_integrator->dt())
795 1690 : for (const auto i : make_range(num_dofs))
796 6100 : for (const auto qp : make_range(nqp))
797 4814 : _ad_u_dotdot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dotdot[i];
798 45 : else if (!_time_integrator)
799 0 : for (const auto qp : make_range(nqp))
800 0 : _ad_u_dotdot[qp] = _u_dotdot[qp];
801 : }
802 :
803 10929951 : if (_need_ad_grad_u_dot)
804 : {
805 0 : _ad_grad_u_dot.resize(nqp);
806 0 : for (const auto qp : make_range(nqp))
807 0 : _ad_grad_u_dot[qp] = _ad_zero;
808 :
809 0 : if (_time_integrator && _time_integrator->dt())
810 : {
811 : // The latter check here is for handling the fact that we have not yet implemented
812 : // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that
813 : // situation we need to default to using the non-ad grad_phi
814 0 : if (_displaced && _current_ad_grad_phi)
815 0 : for (const auto i : make_range(num_dofs))
816 0 : for (const auto qp : make_range(nqp))
817 0 : _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_ad_grad_phi)[i][qp];
818 : else
819 0 : for (const auto i : make_range(num_dofs))
820 0 : for (const auto qp : make_range(nqp))
821 0 : _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_grad_phi)[i][qp];
822 : }
823 0 : else if (!_time_integrator)
824 0 : for (const auto qp : make_range(nqp))
825 0 : _ad_grad_u_dot[qp] = _grad_u_dot[qp];
826 : }
827 : }
828 26243703 : }
829 :
830 : template <>
831 : void
832 0 : MooseVariableData<RealEigenVector>::computeAD(const unsigned int /*num_dofs*/,
833 : const unsigned int /*nqp*/)
834 : {
835 0 : mooseError("AD for array variable has not been implemented");
836 : }
837 :
838 : template <typename OutputType>
839 : void
840 8515175 : MooseVariableData<OutputType>::setDofValue(const OutputData & value, unsigned int index)
841 : {
842 8515175 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
843 8515175 : dof_values[index] = value;
844 8515175 : _has_dof_values = true;
845 :
846 8515175 : auto & u = _vector_tag_u[_solution_tag];
847 8515175 : const auto nqps = u.size();
848 8515175 : const auto ndofs = dof_values.size();
849 45454857 : for (const auto qp : make_range(nqps))
850 36939682 : u[qp] *= 0.;
851 45454857 : for (const auto qp : make_range(nqps))
852 315055554 : for (const auto i : make_range(ndofs))
853 278115872 : u[qp] += (*_phi)[i][qp] * dof_values[i];
854 8515175 : }
855 :
856 : template <typename OutputType>
857 : void
858 651223 : MooseVariableData<OutputType>::setDofValues(const DenseVector<OutputData> & values)
859 : {
860 651223 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
861 3862132 : for (unsigned int i = 0; i < values.size(); i++)
862 3210909 : dof_values[i] = values(i);
863 :
864 651223 : _has_dof_values = true;
865 :
866 651223 : auto & u = _vector_tag_u[_solution_tag];
867 651223 : const auto nqps = u.size();
868 651223 : const auto ndofs = dof_values.size();
869 5572180 : for (const auto qp : make_range(nqps))
870 4920957 : u[qp] *= 0.;
871 5572180 : for (const auto qp : make_range(nqps))
872 36061755 : for (const auto i : make_range(ndofs))
873 31140798 : u[qp] += (*_phi)[i][qp] * dof_values[i];
874 651223 : }
875 :
876 : template <typename OutputType>
877 : void
878 60828997 : MooseVariableData<OutputType>::insertNodalValue(NumericVector<Number> & residual,
879 : const OutputData & v)
880 : {
881 60828997 : residual.set(_nodal_dof_index, v);
882 60828997 : }
883 :
884 : template <>
885 : void
886 313701 : MooseVariableData<RealEigenVector>::insertNodalValue(NumericVector<Number> & residual,
887 : const RealEigenVector & v)
888 : {
889 1116623 : for (const auto j : make_range(_count))
890 802922 : residual.set(_nodal_dof_index + j, v(j));
891 313701 : }
892 :
893 : template <typename OutputType>
894 : typename MooseVariableData<OutputType>::OutputData
895 121377 : MooseVariableData<OutputType>::getNodalValue(const Node & node, Moose::SolutionState state) const
896 : {
897 : mooseAssert(_subproblem.mesh().isSemiLocal(const_cast<Node *>(&node)), "Node is not Semilocal");
898 :
899 : // Make sure that the node has DOFs
900 : /* Note, this is a reproduction of an assert within libMesh::Node::dof_number, this is done to
901 : * produce a better error (see misc/check_error.node_value_off_block) */
902 : mooseAssert(node.n_dofs(_sys.number(), _var_num) > 0,
903 : "Node " << node.id() << " does not contain any dofs for the "
904 : << _sys.system().variable_name(_var_num) << " variable");
905 :
906 121377 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
907 :
908 121377 : switch (state)
909 : {
910 1345 : case Moose::Current:
911 1345 : return (*_sys.currentSolution())(dof);
912 :
913 120032 : case Moose::Old:
914 120032 : return _sys.solutionOld()(dof);
915 :
916 0 : case Moose::Older:
917 0 : return _sys.solutionOlder()(dof);
918 :
919 0 : default:
920 0 : mooseError("PreviousNL not currently supported for getNodalValue");
921 : }
922 : }
923 :
924 : template <>
925 : RealEigenVector
926 0 : MooseVariableData<RealEigenVector>::getNodalValue(const Node & node,
927 : Moose::SolutionState state) const
928 : {
929 : mooseAssert(_subproblem.mesh().isSemiLocal(const_cast<Node *>(&node)), "Node is not Semilocal");
930 :
931 : // Make sure that the node has DOFs
932 : /* Note, this is a reproduction of an assert within libMesh::Node::dof_number, this is done to
933 : * produce a better error (see misc/check_error.node_value_off_block) */
934 : mooseAssert(node.n_dofs(_sys.number(), _var_num) > 0,
935 : "Node " << node.id() << " does not contain any dofs for the "
936 : << _sys.system().variable_name(_var_num) << " variable");
937 :
938 0 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
939 :
940 0 : RealEigenVector v(_count);
941 0 : switch (state)
942 : {
943 0 : case Moose::Current:
944 0 : for (unsigned int i = 0; i < _count; ++i)
945 0 : v(i) = (*_sys.currentSolution())(dof++);
946 0 : break;
947 :
948 0 : case Moose::Old:
949 0 : for (unsigned int i = 0; i < _count; ++i)
950 0 : v(i) = _sys.solutionOld()(dof++);
951 0 : break;
952 :
953 0 : case Moose::Older:
954 0 : for (unsigned int i = 0; i < _count; ++i)
955 0 : v(i) = _sys.solutionOlder()(dof++);
956 0 : break;
957 :
958 0 : default:
959 0 : mooseError("PreviousNL not currently supported for getNodalValue");
960 : }
961 0 : return v;
962 0 : }
963 :
964 : template <typename OutputType>
965 : typename MooseVariableData<OutputType>::OutputData
966 79650 : MooseVariableData<OutputType>::getElementalValue(const Elem * elem,
967 : Moose::SolutionState state,
968 : unsigned int idx) const
969 : {
970 79650 : std::vector<dof_id_type> dof_indices;
971 79650 : _dof_map.dof_indices(elem, dof_indices, _var_num);
972 :
973 79650 : switch (state)
974 : {
975 78546 : case Moose::Current:
976 78546 : return (*_sys.currentSolution())(dof_indices[idx]);
977 :
978 552 : case Moose::Old:
979 552 : return _sys.solutionOld()(dof_indices[idx]);
980 :
981 552 : case Moose::Older:
982 552 : return _sys.solutionOlder()(dof_indices[idx]);
983 :
984 0 : default:
985 0 : mooseError("PreviousNL not currently supported for getElementalValue");
986 : }
987 79650 : }
988 :
989 : template <>
990 : RealEigenVector
991 0 : MooseVariableData<RealEigenVector>::getElementalValue(const Elem * elem,
992 : Moose::SolutionState state,
993 : unsigned int idx) const
994 : {
995 0 : std::vector<dof_id_type> dof_indices;
996 0 : _dof_map.dof_indices(elem, dof_indices, _var_num);
997 :
998 0 : dof_id_type dof = dof_indices[idx];
999 :
1000 0 : RealEigenVector v(_count);
1001 :
1002 0 : switch (state)
1003 : {
1004 0 : case Moose::Current:
1005 0 : for (unsigned int i = 0; i < _count; ++i)
1006 0 : v(i) = (*_sys.currentSolution())(dof++);
1007 0 : break;
1008 :
1009 0 : case Moose::Old:
1010 0 : for (unsigned int i = 0; i < _count; ++i)
1011 0 : v(i) = _sys.solutionOld()(dof++);
1012 0 : break;
1013 :
1014 0 : case Moose::Older:
1015 0 : for (unsigned int i = 0; i < _count; ++i)
1016 0 : v(i) = _sys.solutionOlder()(dof++);
1017 0 : break;
1018 :
1019 0 : default:
1020 0 : mooseError("PreviousNL not currently supported for getElementalValue");
1021 : }
1022 0 : return v;
1023 0 : }
1024 :
1025 : template <typename OutputType>
1026 : void
1027 21514 : MooseVariableData<OutputType>::getDofIndices(const Elem * elem,
1028 : std::vector<dof_id_type> & dof_indices) const
1029 : {
1030 21514 : _dof_map.dof_indices(elem, dof_indices, _var_num);
1031 21514 : }
1032 :
1033 : template <typename OutputType>
1034 : void
1035 0 : MooseVariableData<OutputType>::addSolution(NumericVector<Number> & sol,
1036 : const DenseVector<Number> & v) const
1037 : {
1038 0 : sol.add_vector(v, _dof_indices);
1039 0 : }
1040 :
1041 : template <>
1042 : void
1043 3060 : MooseVariableData<RealEigenVector>::addSolution(NumericVector<Number> & sol,
1044 : const DenseVector<Number> & v) const
1045 : {
1046 3060 : unsigned int p = 0;
1047 3060 : const auto num_dofs = _dof_indices.size();
1048 9180 : for (const auto j : make_range(_count))
1049 : {
1050 6120 : unsigned int inc = (isNodal() ? j : j * num_dofs);
1051 30600 : for (const auto i : make_range(num_dofs))
1052 24480 : sol.add(_dof_indices[i] + inc, v(p++));
1053 : }
1054 3060 : }
1055 :
1056 : template <typename OutputType>
1057 : const typename MooseVariableData<OutputType>::DoFValue &
1058 327 : MooseVariableData<OutputType>::dofValuesDot() const
1059 : {
1060 327 : if (_sys.solutionUDot())
1061 : {
1062 327 : _need_dof_values_dot = true;
1063 327 : return _dof_values_dot;
1064 : }
1065 : else
1066 0 : mooseError("MooseVariableData: Time derivative of solution (`u_dot`) is not stored. Please set "
1067 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
1068 : }
1069 :
1070 : template <typename OutputType>
1071 : const typename MooseVariableData<OutputType>::DoFValue &
1072 4 : MooseVariableData<OutputType>::dofValuesDotDot() const
1073 : {
1074 4 : if (_sys.solutionUDotDot())
1075 : {
1076 4 : _need_dof_values_dotdot = true;
1077 4 : return _dof_values_dotdot;
1078 : }
1079 : else
1080 0 : mooseError("MooseVariableData: Second time derivative of solution (`u_dotdot`) is not stored. "
1081 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
1082 : "`u_dotdot`.");
1083 : }
1084 :
1085 : template <typename OutputType>
1086 : const typename MooseVariableData<OutputType>::DoFValue &
1087 0 : MooseVariableData<OutputType>::dofValuesDotOld() const
1088 : {
1089 0 : if (_sys.solutionUDotOld())
1090 : {
1091 0 : _need_dof_values_dot_old = true;
1092 0 : return _dof_values_dot_old;
1093 : }
1094 : else
1095 0 : mooseError("MooseVariableData: Old time derivative of solution (`u_dot_old`) is not stored. "
1096 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
1097 : "`u_dot_old`.");
1098 : }
1099 :
1100 : template <typename OutputType>
1101 : const typename MooseVariableData<OutputType>::DoFValue &
1102 0 : MooseVariableData<OutputType>::dofValuesDotDotOld() const
1103 : {
1104 0 : if (_sys.solutionUDotDotOld())
1105 : {
1106 0 : _need_dof_values_dotdot_old = true;
1107 0 : return _dof_values_dotdot_old;
1108 : }
1109 : else
1110 0 : mooseError("MooseVariableData: Old second time derivative of solution (`u_dotdot_old`) is not "
1111 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
1112 : "requesting `u_dotdot_old`.");
1113 : }
1114 :
1115 : template <typename OutputType>
1116 : const MooseArray<Number> &
1117 131 : MooseVariableData<OutputType>::dofValuesDuDotDu() const
1118 : {
1119 131 : _need_dof_du_dot_du = true;
1120 131 : return _dof_du_dot_du;
1121 : }
1122 :
1123 : template <typename OutputType>
1124 : const MooseArray<Number> &
1125 0 : MooseVariableData<OutputType>::dofValuesDuDotDotDu() const
1126 : {
1127 0 : _need_dof_du_dotdot_du = true;
1128 0 : return _dof_du_dotdot_du;
1129 : }
1130 :
1131 : template <typename OutputType>
1132 : void
1133 1776 : MooseVariableData<OutputType>::computeIncrementAtQps(const NumericVector<Number> & increment_vec)
1134 : {
1135 1776 : unsigned int nqp = _qrule->n_points();
1136 :
1137 1776 : _increment.resize(nqp);
1138 : // Compute the increment at each quadrature point
1139 1776 : unsigned int num_dofs = _dof_indices.size();
1140 15142 : for (const auto qp : make_range(nqp))
1141 : {
1142 13366 : _increment[qp] = 0.;
1143 130356 : for (const auto i : make_range(num_dofs))
1144 116990 : _increment[qp] += (*_phi)[i][qp] * increment_vec(_dof_indices[i]);
1145 : }
1146 1776 : }
1147 :
1148 : template <>
1149 : void
1150 0 : MooseVariableData<RealEigenVector>::computeIncrementAtQps(
1151 : const NumericVector<Number> & increment_vec)
1152 : {
1153 0 : unsigned int nqp = _qrule->n_points();
1154 :
1155 0 : _increment.resize(nqp);
1156 : // Compute the increment at each quadrature point
1157 0 : unsigned int num_dofs = _dof_indices.size();
1158 0 : if (isNodal())
1159 : {
1160 0 : for (const auto qp : make_range(nqp))
1161 : {
1162 0 : for (const auto i : make_range(num_dofs))
1163 0 : for (const auto j : make_range(_count))
1164 0 : _increment[qp](j) += (*_phi)[i][qp] * increment_vec(_dof_indices[i] + j);
1165 : }
1166 : }
1167 : else
1168 : {
1169 0 : for (const auto qp : make_range(nqp))
1170 : {
1171 0 : unsigned int n = 0;
1172 0 : for (const auto j : make_range(_count))
1173 0 : for (const auto i : make_range(num_dofs))
1174 : {
1175 0 : _increment[qp](j) += (*_phi)[i][qp] * increment_vec(_dof_indices[i] + n);
1176 0 : n += num_dofs;
1177 : }
1178 : }
1179 : }
1180 0 : }
1181 :
1182 : template <typename OutputType>
1183 : void
1184 7441 : MooseVariableData<OutputType>::computeIncrementAtNode(const NumericVector<Number> & increment_vec)
1185 : {
1186 7441 : if (!isNodal())
1187 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1188 :
1189 7441 : _increment.resize(1);
1190 :
1191 : // Compute the increment for the current DOF
1192 7441 : _increment[0] = increment_vec(_dof_indices[0]);
1193 7441 : }
1194 :
1195 : template <>
1196 : void
1197 0 : MooseVariableData<RealEigenVector>::computeIncrementAtNode(
1198 : const NumericVector<Number> & increment_vec)
1199 : {
1200 0 : if (!isNodal())
1201 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1202 :
1203 0 : _increment.resize(1);
1204 :
1205 : // Compute the increment for the current DOF
1206 0 : if (isNodal())
1207 0 : for (unsigned int j = 0; j < _count; j++)
1208 0 : _increment[0](j) = increment_vec(_dof_indices[0] + j);
1209 : else
1210 : {
1211 0 : unsigned int n = 0;
1212 0 : const auto n_dof_indices = _dof_indices.size();
1213 0 : for (const auto j : make_range(_count))
1214 : {
1215 0 : _increment[0](j) = increment_vec(_dof_indices[0] + n);
1216 0 : n += n_dof_indices;
1217 : }
1218 : }
1219 0 : }
1220 :
1221 : template <typename OutputType>
1222 : const OutputType &
1223 0 : MooseVariableData<OutputType>::nodalValueDot() const
1224 : {
1225 0 : if (isNodal())
1226 : {
1227 0 : if (_sys.solutionUDot())
1228 : {
1229 0 : _need_dof_values_dot = true;
1230 0 : return _nodal_value_dot;
1231 : }
1232 : else
1233 0 : mooseError(
1234 : "MooseVariableData: Time derivative of solution (`u_dot`) is not stored. Please set "
1235 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
1236 : }
1237 : else
1238 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1239 0 : var().name(),
1240 : "' is not nodal.");
1241 : }
1242 :
1243 : template <typename OutputType>
1244 : const OutputType &
1245 0 : MooseVariableData<OutputType>::nodalValueDotDot() const
1246 : {
1247 0 : if (isNodal())
1248 : {
1249 0 : if (_sys.solutionUDotDot())
1250 : {
1251 0 : _need_dof_values_dotdot = true;
1252 0 : return _nodal_value_dotdot;
1253 : }
1254 : else
1255 0 : mooseError(
1256 : "MooseVariableData: Second time derivative of solution (`u_dotdot`) is not stored. "
1257 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
1258 : "`u_dotdot`.");
1259 : }
1260 : else
1261 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1262 0 : var().name(),
1263 : "' is not nodal.");
1264 : }
1265 :
1266 : template <typename OutputType>
1267 : const OutputType &
1268 0 : MooseVariableData<OutputType>::nodalValueDotOld() const
1269 : {
1270 0 : if (isNodal())
1271 : {
1272 0 : if (_sys.solutionUDotOld())
1273 : {
1274 0 : _need_dof_values_dot_old = true;
1275 0 : return _nodal_value_dot_old;
1276 : }
1277 : else
1278 0 : mooseError("MooseVariableData: Old time derivative of solution (`u_dot_old`) is not stored. "
1279 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
1280 : "`u_dot_old`.");
1281 : }
1282 : else
1283 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1284 0 : var().name(),
1285 : "' is not nodal.");
1286 : }
1287 :
1288 : template <typename OutputType>
1289 : const OutputType &
1290 0 : MooseVariableData<OutputType>::nodalValueDotDotOld() const
1291 : {
1292 0 : if (isNodal())
1293 : {
1294 0 : if (_sys.solutionUDotDotOld())
1295 : {
1296 0 : _need_dof_values_dotdot_old = true;
1297 0 : return _nodal_value_dotdot_old;
1298 : }
1299 : else
1300 0 : mooseError(
1301 : "MooseVariableData: Old second time derivative of solution (`u_dotdot_old`) is not "
1302 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
1303 : "requesting `u_dotdot_old`.");
1304 : }
1305 : else
1306 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1307 0 : var().name(),
1308 : "' is not nodal.");
1309 : }
1310 :
1311 : template <typename OutputType>
1312 : void
1313 303977080 : MooseVariableData<OutputType>::computeNodalValues()
1314 : {
1315 303977080 : if (_has_dof_indices)
1316 : {
1317 303976840 : fetchDoFValues();
1318 303976840 : assignNodalValue();
1319 :
1320 303976840 : if (_need_ad)
1321 3446940 : fetchADNodalValues();
1322 : }
1323 : else
1324 240 : zeroSizeDofValues();
1325 303977080 : }
1326 :
1327 : template <typename OutputType>
1328 : void
1329 3446940 : MooseVariableData<OutputType>::fetchADNodalValues()
1330 : {
1331 3446940 : auto n = _dof_indices.size();
1332 : libmesh_assert(n);
1333 3446940 : _ad_dof_values.resize(n);
1334 :
1335 3446940 : if (_need_ad_u_dot)
1336 617628 : _ad_dofs_dot.resize(n);
1337 3446940 : if (_need_ad_u_dotdot)
1338 1396 : _ad_dofs_dotdot.resize(n);
1339 :
1340 3446940 : const bool do_derivatives =
1341 3446940 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
1342 :
1343 7073456 : for (decltype(n) i = 0; i < n; ++i)
1344 : {
1345 3626516 : _ad_dof_values[i] = _vector_tags_dof_u[_solution_tag][i];
1346 3626516 : if (do_derivatives)
1347 1513350 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
1348 3626516 : assignADNodalValue(_ad_dof_values[i], i);
1349 :
1350 3626516 : if (_need_ad_u_dot)
1351 : {
1352 626012 : if (_time_integrator && _time_integrator->dt())
1353 : {
1354 626012 : _ad_dofs_dot[i] = _ad_dof_values[i];
1355 1250628 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
1356 626012 : _dof_indices[i],
1357 626012 : _need_ad_u_dotdot ? _ad_dofs_dotdot[i]
1358 : : _ad_real_dummy);
1359 : }
1360 : // Executing something with a time derivative at initial should not put a NaN
1361 0 : else if (_time_integrator && !_time_integrator->dt())
1362 0 : _ad_dofs_dot[i] = 0.;
1363 : else
1364 0 : mooseError("AD nodal time derivatives not implemented for variables without a time "
1365 : "integrator (auxiliary variables)");
1366 : }
1367 : }
1368 3446940 : }
1369 :
1370 : template <>
1371 : void
1372 0 : MooseVariableData<RealEigenVector>::fetchADNodalValues()
1373 : {
1374 0 : mooseError("I do not know how to support AD with array variables");
1375 : }
1376 :
1377 : template <>
1378 : void
1379 3298751 : MooseVariableData<Real>::assignADNodalValue(const ADReal & value, const unsigned int &)
1380 : {
1381 3298751 : _ad_nodal_value = value;
1382 3298751 : }
1383 :
1384 : template <>
1385 : void
1386 327765 : MooseVariableData<RealVectorValue>::assignADNodalValue(const ADReal & value,
1387 : const unsigned int & component)
1388 : {
1389 327765 : _ad_nodal_value(component) = value;
1390 327765 : }
1391 :
1392 : template <typename OutputType>
1393 : void
1394 2230440 : MooseVariableData<OutputType>::prepareIC()
1395 : {
1396 2230440 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1397 2230440 : _vector_tags_dof_u[_solution_tag].resize(_dof_indices.size());
1398 :
1399 2230440 : unsigned int nqp = _qrule->n_points();
1400 2230440 : _vector_tag_u[_solution_tag].resize(nqp);
1401 2230440 : }
1402 :
1403 : template <typename OutputType>
1404 : void
1405 535686926 : MooseVariableData<OutputType>::prepare()
1406 : {
1407 535686926 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1408 535686926 : _has_dof_values = false;
1409 :
1410 : // FIXME: remove this when the Richard's module is migrated to use the new NodalCoupleable
1411 : // interface.
1412 535686926 : _has_dof_indices = _dof_indices.size();
1413 535686926 : }
1414 :
1415 : template <typename OutputType>
1416 : void
1417 345912904 : MooseVariableData<OutputType>::reinitNode()
1418 : {
1419 345912904 : if (std::size_t n_dofs = _node->n_dofs(_sys.number(), _var_num))
1420 : {
1421 303969889 : _dof_indices.resize(n_dofs);
1422 609044450 : for (std::size_t i = 0; i < n_dofs; ++i)
1423 305074561 : _dof_indices[i] = _node->dof_number(_sys.number(), _var_num, i);
1424 : // For standard variables. _nodal_dof_index is retrieved by nodalDofIndex() which is used in
1425 : // NodalBC for example
1426 303969889 : _nodal_dof_index = _dof_indices[0];
1427 303969889 : _has_dof_indices = true;
1428 : }
1429 : else
1430 : {
1431 41943015 : _dof_indices.clear(); // Clear these so Assembly doesn't think there's dofs here
1432 41943015 : _has_dof_indices = false;
1433 : }
1434 345912904 : }
1435 :
1436 : template <typename OutputType>
1437 : void
1438 167757020 : MooseVariableData<OutputType>::reinitAux()
1439 : {
1440 : /* FIXME: this method is only for elemental auxiliary variables, so
1441 : * we may want to rename it */
1442 167757020 : if (_elem)
1443 : {
1444 166695848 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1445 166695848 : if (_elem->n_dofs(_sys.number(), _var_num) > 0)
1446 : {
1447 : // FIXME: check if the following is equivalent with '_nodal_dof_index = _dof_indices[0];'?
1448 165739244 : _nodal_dof_index = _elem->dof_number(_sys.number(), _var_num, 0);
1449 :
1450 165739244 : fetchDoFValues();
1451 :
1452 165739244 : const auto num_dofs = _dof_indices.size();
1453 1364947458 : for (auto & dof_u : _vector_tags_dof_u)
1454 1199208214 : dof_u.resize(num_dofs);
1455 :
1456 497229780 : for (auto & dof_u : _matrix_tags_dof_u)
1457 331490536 : dof_u.resize(num_dofs);
1458 :
1459 165739244 : _has_dof_indices = true;
1460 : }
1461 : else
1462 956604 : _has_dof_indices = false;
1463 : }
1464 : else
1465 1061172 : _has_dof_indices = false;
1466 167757020 : }
1467 :
1468 : template <typename OutputType>
1469 : void
1470 7866 : MooseVariableData<OutputType>::reinitNodes(const std::vector<dof_id_type> & nodes)
1471 : {
1472 7866 : _dof_indices.clear();
1473 42615 : for (const auto & node_id : nodes)
1474 : {
1475 34749 : auto && nd = _subproblem.mesh().getMesh().query_node_ptr(node_id);
1476 34749 : if (nd && (_subproblem.mesh().isSemiLocal(const_cast<Node *>(nd))))
1477 : {
1478 34683 : if (nd->n_dofs(_sys.number(), _var_num) > 0)
1479 : {
1480 33492 : dof_id_type dof = nd->dof_number(_sys.number(), _var_num, 0);
1481 33492 : _dof_indices.push_back(dof);
1482 : }
1483 : }
1484 : }
1485 :
1486 7866 : if (!_dof_indices.empty())
1487 7626 : _has_dof_indices = true;
1488 : else
1489 240 : _has_dof_indices = false;
1490 7866 : }
1491 :
1492 : template class MooseVariableData<Real>;
1493 : template class MooseVariableData<RealVectorValue>;
1494 : template class MooseVariableData<RealEigenVector>;
|
