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 543225 : 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 1086450 : _fe_type(var.feType()),
39 543225 : _var_num(var.number()),
40 543225 : _assembly(_subproblem.assembly(_tid, var.kind() == Moose::VAR_SOLVER ? sys.number() : 0)),
41 543225 : _element_type(element_type),
42 543225 : _ad_zero(0),
43 543225 : _need_ad_u_dot(false),
44 543225 : _need_ad_u_dotdot(false),
45 543225 : _need_second(false),
46 543225 : _need_second_old(false),
47 543225 : _need_second_older(false),
48 543225 : _need_second_previous_nl(false),
49 543225 : _need_curl(false),
50 543225 : _need_curl_old(false),
51 543225 : _need_curl_older(false),
52 543225 : _need_div(false),
53 543225 : _need_div_old(false),
54 543225 : _need_div_older(false),
55 543225 : _need_ad(false),
56 543225 : _need_ad_u(false),
57 543225 : _need_ad_grad_u(false),
58 543225 : _need_ad_grad_u_dot(false),
59 543225 : _need_ad_second_u(false),
60 543225 : _need_ad_curl_u(false),
61 543225 : _has_dof_indices(false),
62 543225 : _qrule(qrule_in),
63 543225 : _qrule_face(qrule_face_in),
64 543225 : _use_dual(var.useDual()),
65 543225 : _second_phi_assembly_method(nullptr),
66 543225 : _second_phi_face_assembly_method(nullptr),
67 543225 : _curl_phi_assembly_method(nullptr),
68 543225 : _curl_phi_face_assembly_method(nullptr),
69 543225 : _div_phi_assembly_method(nullptr),
70 543225 : _div_phi_face_assembly_method(nullptr),
71 543225 : _ad_grad_phi_assembly_method(nullptr),
72 543225 : _ad_grad_phi_face_assembly_method(nullptr),
73 543225 : _time_integrator(nullptr),
74 543225 : _node(node),
75 543225 : _elem(elem),
76 543225 : _displaced(dynamic_cast<const DisplacedSystem *>(&_sys) ? true : false),
77 2716125 : _current_side(_assembly.side())
78 : {
79 543225 : _continuity = FEInterface::get_continuity(_fe_type);
80 :
81 543225 : _is_nodal = (_continuity == C_ZERO || _continuity == C_ONE);
82 :
83 543225 : _time_integrator = _sys.queryTimeIntegrator(_var_num);
84 :
85 : // Initialize AD zero with zero derivatives
86 543225 : const auto old_do = ADReal::do_derivatives;
87 543225 : ADReal::do_derivatives = true;
88 543225 : _ad_zero = 0.;
89 543225 : ADReal::do_derivatives = old_do;
90 :
91 543225 : switch (_element_type)
92 : {
93 181075 : case Moose::ElementType::Element:
94 : {
95 181075 : _phi_assembly_method = &Assembly::fePhi<OutputShape>;
96 181075 : _phi_face_assembly_method = &Assembly::fePhiFace<OutputShape>;
97 181075 : _grad_phi_assembly_method = &Assembly::feGradPhi<OutputShape>;
98 181075 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFace<OutputShape>;
99 181075 : _second_phi_assembly_method = &Assembly::feSecondPhi<OutputShape>;
100 181075 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFace<OutputShape>;
101 181075 : _curl_phi_assembly_method = &Assembly::feCurlPhi<OutputShape>;
102 181075 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFace<OutputShape>;
103 181075 : _div_phi_assembly_method = &Assembly::feDivPhi<OutputShape>;
104 181075 : _div_phi_face_assembly_method = &Assembly::feDivPhiFace<OutputShape>;
105 181075 : _ad_grad_phi_assembly_method = &Assembly::feADGradPhi<OutputShape>;
106 181075 : _ad_grad_phi_face_assembly_method = &Assembly::feADGradPhiFace<OutputShape>;
107 :
108 181075 : _ad_grad_phi = &_ad_grad_phi_assembly_method(_assembly, _fe_type);
109 181075 : _ad_grad_phi_face = &_ad_grad_phi_face_assembly_method(_assembly, _fe_type);
110 181075 : break;
111 : }
112 181075 : case Moose::ElementType::Neighbor:
113 : {
114 181075 : _phi_assembly_method = &Assembly::fePhiNeighbor<OutputShape>;
115 181075 : _phi_face_assembly_method = &Assembly::fePhiFaceNeighbor<OutputShape>;
116 181075 : _grad_phi_assembly_method = &Assembly::feGradPhiNeighbor<OutputShape>;
117 181075 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFaceNeighbor<OutputShape>;
118 181075 : _second_phi_assembly_method = &Assembly::feSecondPhiNeighbor<OutputShape>;
119 181075 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFaceNeighbor<OutputShape>;
120 181075 : _curl_phi_assembly_method = &Assembly::feCurlPhiNeighbor<OutputShape>;
121 181075 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFaceNeighbor<OutputShape>;
122 181075 : _div_phi_assembly_method = &Assembly::feDivPhiNeighbor<OutputShape>;
123 181075 : _div_phi_face_assembly_method = &Assembly::feDivPhiFaceNeighbor<OutputShape>;
124 :
125 181075 : _ad_grad_phi = nullptr;
126 181075 : _ad_grad_phi_face = nullptr;
127 181075 : break;
128 : }
129 181075 : case Moose::ElementType::Lower:
130 : {
131 181075 : 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 180929 : _phi_assembly_method = &Assembly::fePhiLower<OutputType>;
141 180929 : _phi_face_assembly_method = &Assembly::fePhiLower<OutputType>; // Place holder
142 180929 : _grad_phi_assembly_method = &Assembly::feGradPhiLower<OutputType>;
143 180929 : _grad_phi_face_assembly_method = &Assembly::feGradPhiLower<OutputType>; // Place holder
144 : }
145 :
146 181075 : _ad_grad_phi = nullptr;
147 181075 : _ad_grad_phi_face = nullptr;
148 181075 : break;
149 : }
150 : }
151 543225 : _phi = &_phi_assembly_method(_assembly, _fe_type);
152 543225 : _phi_face = &_phi_face_assembly_method(_assembly, _fe_type);
153 543225 : _grad_phi = &_grad_phi_assembly_method(_assembly, _fe_type);
154 543225 : _grad_phi_face = &_grad_phi_face_assembly_method(_assembly, _fe_type);
155 543225 : }
156 :
157 : template <typename OutputType>
158 : void
159 793899748 : MooseVariableData<OutputType>::setGeometry(Moose::GeometryType gm_type)
160 : {
161 793899748 : switch (gm_type)
162 : {
163 748010118 : case Moose::Volume:
164 : {
165 748010118 : _current_qrule = _qrule;
166 748010118 : _current_phi = _phi;
167 748010118 : _current_grad_phi = _grad_phi;
168 748010118 : _current_second_phi = _second_phi;
169 748010118 : _current_curl_phi = _curl_phi;
170 748010118 : _current_div_phi = _div_phi;
171 748010118 : _current_ad_grad_phi = _ad_grad_phi;
172 748010118 : break;
173 : }
174 45889630 : case Moose::Face:
175 : {
176 45889630 : _current_qrule = _qrule_face;
177 45889630 : _current_phi = _phi_face;
178 45889630 : _current_grad_phi = _grad_phi_face;
179 45889630 : _current_second_phi = _second_phi_face;
180 45889630 : _current_curl_phi = _curl_phi_face;
181 45889630 : _current_div_phi = _div_phi_face;
182 45889630 : _current_ad_grad_phi = _ad_grad_phi_face;
183 45889630 : break;
184 : }
185 : }
186 793899748 : }
187 :
188 : template <typename OutputType>
189 : const typename MooseVariableData<OutputType>::FieldVariableValue &
190 17361 : MooseVariableData<OutputType>::uDot() const
191 : {
192 17361 : if (_sys.solutionUDot())
193 : {
194 17361 : _need_u_dot = true;
195 17361 : 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 793095923 : MooseVariableData<OutputType>::computeValuesInternal()
429 : {
430 793095923 : const auto num_dofs = _dof_indices.size();
431 793095923 : if (num_dofs > 0)
432 689460369 : fetchDoFValues();
433 :
434 793095923 : const bool is_transient = _subproblem.isTransient();
435 793095923 : const auto nqp = _current_qrule->n_points();
436 793095923 : 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 4563548 : if (_qrule == _current_qrule)
442 : {
443 4328546 : _mapped_grad_phi.resize(num_dofs);
444 20317451 : for (const auto i : make_range(num_dofs))
445 : {
446 15988905 : _mapped_grad_phi[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
447 106056274 : for (const auto qp : make_range(nqp))
448 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
449 180134738 : new (&_mapped_grad_phi[i][qp])
450 90067369 : Eigen::Map<RealDIMValue>(const_cast<Real *>(&(*_current_grad_phi)[i][qp](0)));
451 : }
452 : }
453 : else
454 : {
455 235002 : _mapped_grad_phi_face.resize(num_dofs);
456 928630 : for (const auto i : make_range(num_dofs))
457 : {
458 693628 : _mapped_grad_phi_face[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
459 2401164 : for (const auto qp : make_range(nqp))
460 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
461 3415072 : new (&_mapped_grad_phi_face[i][qp])
462 1707536 : 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 2327098332 : const auto fill = [this, nqp, num_dofs](auto & dest, const auto & phi, const auto & dof_values)
478 : {
479 : if constexpr (monomial)
480 178970714 : libmesh_ignore(num_dofs);
481 :
482 : // Deduce OutputType
483 1534002409 : constexpr bool is_real = std::is_same_v<OutputType, Real>;
484 1534002409 : constexpr bool is_real_vector = std::is_same_v<OutputType, RealVectorValue>;
485 1534002409 : 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 1524691979 : 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 1534002409 : constexpr bool is_value = std::is_same_v<dest_array_type, OutputType>;
495 1534002409 : constexpr bool is_gradient = std::is_same_v<dest_array_type, OutputGradient>;
496 1534002409 : constexpr bool is_second = std::is_same_v<dest_array_type, OutputSecond>;
497 1534002409 : 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 14085730427 : const auto set_zero = [this, &dest](const auto qp)
503 : {
504 : if constexpr (!is_eigen)
505 6232564727 : libmesh_ignore(this);
506 :
507 : if constexpr (is_real || is_real_vector)
508 6232564727 : dest[qp] = 0;
509 : else if constexpr (is_eigen)
510 : {
511 : if constexpr (is_value)
512 30313885 : dest[qp].setZero(this->_count);
513 : else if constexpr (is_gradient)
514 12985397 : 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 58016304170 : 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 17823873277 : dest[qp] += phi[i][qp] * dof_values[i];
531 : else if constexpr (is_gradient || is_second)
532 10380919691 : 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 290863300 : for (const auto d : make_range(Moose::dim))
541 218147475 : 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 1534002409 : 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 178970714 : set_zero(0);
563 178970714 : accumulate(0, 0);
564 808175498 : for (unsigned int qp = 1; qp < nqp; ++qp)
565 629204784 : dest[qp] = dest[0];
566 : }
567 : // Non monomial case
568 : else
569 : {
570 7451924990 : for (const auto qp : make_range(nqp))
571 6096893295 : set_zero(qp);
572 7020975571 : for (const auto i : make_range(num_dofs))
573 33764481955 : for (const auto qp : make_range(nqp))
574 28098538079 : accumulate(i, qp);
575 : }
576 : };
577 :
578 : // Curl
579 793095923 : if (_need_curl)
580 79901 : fill(_curl_u, *_current_curl_phi, _vector_tags_dof_u[_solution_tag]);
581 793095923 : if (is_transient && _need_curl_old)
582 0 : fill(_curl_u_old, *_current_curl_phi, _vector_tags_dof_u[_old_solution_tag]);
583 :
584 : // Div
585 793095923 : if (_need_div)
586 396144 : fill(_div_u, *_current_div_phi, _vector_tags_dof_u[_solution_tag]);
587 793095923 : 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 793095923 : if (_need_second)
592 92492 : fill(_second_u, *_current_second_phi, _vector_tags_dof_u[_solution_tag]);
593 793095923 : 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 1596373646 : for (auto tag : _required_vector_tags)
599 : {
600 803277723 : if (_need_vector_tag_u[tag] && _sys.hasVector(tag))
601 : {
602 : mooseAssert(_sys.getVector(tag).closed(), "Vector should be closed");
603 800450065 : fill(_vector_tag_u[tag], *_current_phi, _vector_tags_dof_u[tag]);
604 : }
605 803277723 : if (_need_vector_tag_grad[tag] && _sys.hasVector(tag))
606 : {
607 : mooseAssert(_sys.getVector(tag).closed(), "Vector should be closed");
608 457640020 : fill(_vector_tag_grad[tag], *_current_grad_phi, _vector_tags_dof_u[tag]);
609 : }
610 : }
611 :
612 : // Matrix tags
613 793096355 : for (auto tag : active_coupleable_matrix_tags)
614 432 : if (_need_matrix_tag_u[tag])
615 144 : fill(_matrix_tag_u[tag], *_current_phi, _matrix_tags_dof_u[tag]);
616 :
617 : // Derivatives and old values
618 793095923 : if (is_transient)
619 : {
620 596718567 : if (_need_second_old)
621 0 : fill(_second_u_old, *_current_second_phi, _vector_tags_dof_u[_old_solution_tag]);
622 596718567 : if (_need_second_older)
623 0 : fill(_second_u_older, *_current_second_phi, _vector_tags_dof_u[_older_solution_tag]);
624 596718567 : if (_need_u_dot)
625 275195165 : fill(_u_dot, *_current_phi, _dof_values_dot);
626 596718567 : if (_need_u_dotdot)
627 7145 : fill(_u_dotdot, *_current_phi, _dof_values_dotdot);
628 596718567 : if (_need_u_dot_old)
629 0 : fill(_u_dot_old, *_current_phi, _dof_values_dot_old);
630 596718567 : if (_need_u_dotdot_old)
631 0 : fill(_u_dotdot_old, *_current_phi, _dof_values_dotdot_old);
632 :
633 596718567 : if (_need_du_dot_du)
634 : {
635 275332473 : _du_dot_du.resize(nqp);
636 1470955998 : for (const auto i : make_range(num_dofs))
637 6888199300 : for (const auto qp : make_range(nqp))
638 5692575775 : _du_dot_du[qp] = _dof_du_dot_du[i];
639 : }
640 596718567 : if (_need_du_dotdot_du)
641 : {
642 5822 : _du_dotdot_du.resize(nqp);
643 49276 : for (const auto i : make_range(num_dofs))
644 215932 : for (const auto qp : make_range(nqp))
645 172478 : _du_dotdot_du[qp] = _dof_du_dotdot_du[i];
646 : }
647 :
648 596718567 : if (_need_grad_dot)
649 141333 : fill(_grad_u_dot, *_current_grad_phi, _dof_values_dot);
650 596718567 : if (_need_grad_dotdot)
651 0 : fill(_grad_u_dotdot, *_current_grad_phi, _dof_values_dotdot);
652 : }
653 :
654 : // Automatic differentiation, not for eigen
655 : if constexpr (!std::is_same_v<OutputType, RealEigenVector>)
656 788532375 : if (_need_ad)
657 28392888 : computeAD(num_dofs, nqp);
658 793095923 : }
659 :
660 : template <typename OutputType>
661 : void
662 625383684 : MooseVariableData<OutputType>::computeValues()
663 : {
664 625383684 : computeValuesInternal</* monomial = */ false>();
665 625383684 : }
666 :
667 : template <typename OutputType>
668 : void
669 168597866 : MooseVariableData<OutputType>::computeMonomialValues()
670 : {
671 168597866 : if (_dof_indices.size() == 0)
672 803825 : return;
673 :
674 : // Monomial optimizations are not appropriate after p-refinement
675 167794041 : if (_elem->p_level())
676 81802 : computeValues();
677 : else
678 167712239 : computeValuesInternal</* monomial = */ true>();
679 : }
680 :
681 : template <typename OutputType>
682 : void
683 28392888 : MooseVariableData<OutputType>::computeAD(const unsigned int num_dofs, const unsigned int nqp)
684 : {
685 28392888 : const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys);
686 :
687 28392888 : _ad_dof_values.resize(num_dofs);
688 152107456 : for (const auto i : make_range(num_dofs))
689 123714568 : _ad_dof_values[i] = _vector_tags_dof_u[_solution_tag][i];
690 : // NOTE! You have to do this AFTER setting the value!
691 28392888 : if (do_derivatives)
692 31872362 : for (const auto i : make_range(num_dofs))
693 26397896 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
694 :
695 28392888 : if (_need_ad_u)
696 : {
697 26908828 : _ad_u.resize(nqp);
698 115377874 : for (const auto qp : make_range(nqp))
699 88469046 : _ad_u[qp] = _ad_zero;
700 :
701 145843566 : for (const auto i : make_range(num_dofs))
702 573087747 : for (const auto qp : make_range(nqp))
703 454153009 : _ad_u[qp] += _ad_dof_values[i] * (*_current_phi)[i][qp];
704 : }
705 :
706 28392888 : if (_need_ad_grad_u)
707 : {
708 21513814 : _ad_grad_u.resize(nqp);
709 94969617 : for (const auto qp : make_range(nqp))
710 73455803 : _ad_grad_u[qp] = _ad_zero;
711 :
712 : // The latter check here is for handling the fact that we have not yet implemented
713 : // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that
714 : // situation we need to default to using the non-ad grad_phi
715 21513814 : if (_displaced && _current_ad_grad_phi)
716 3029251 : for (const auto i : make_range(num_dofs))
717 16544328 : for (const auto qp : make_range(nqp))
718 14621539 : _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_ad_grad_phi)[i][qp];
719 : else
720 105257657 : for (const auto i : make_range(num_dofs))
721 435716851 : for (const auto qp : make_range(nqp))
722 351419777 : _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_grad_phi)[i][qp];
723 : }
724 :
725 28392888 : if (_need_ad_second_u)
726 : {
727 0 : _ad_second_u.resize(nqp);
728 0 : for (const auto qp : make_range(nqp))
729 0 : _ad_second_u[qp] = _ad_zero;
730 :
731 0 : for (const auto i : make_range(num_dofs))
732 0 : for (const auto qp : make_range(nqp))
733 : // Note that this will not carry any derivatives with respect to displacements because
734 : // those calculations have not yet been implemented in Assembly
735 0 : _ad_second_u[qp] += _ad_dof_values[i] * (*_current_second_phi)[i][qp];
736 : }
737 :
738 28392888 : if (_need_ad_curl_u)
739 : {
740 2860555 : _ad_curl_u.resize(nqp);
741 13377115 : for (const auto qp : make_range(nqp))
742 10516560 : _ad_curl_u[qp] = _ad_zero;
743 :
744 14380175 : for (const auto i : make_range(num_dofs))
745 54075460 : for (const auto qp : make_range(nqp))
746 : {
747 : mooseAssert(_current_curl_phi,
748 : "We're requiring a curl calculation but have not set a curl shape function!");
749 :
750 : // Note that the current version of _ad_curl_u is not yet implemented for mesh displacement
751 42555840 : _ad_curl_u[qp] += _ad_dof_values[i] * (*_current_curl_phi)[i][qp];
752 : }
753 : }
754 :
755 28392888 : bool is_transient = _subproblem.isTransient();
756 28392888 : if (is_transient)
757 : {
758 12321742 : if (_need_ad_u_dot)
759 : {
760 2556451 : _ad_dofs_dot.resize(num_dofs);
761 2556451 : if (_need_ad_u_dotdot)
762 449 : _ad_dofs_dotdot.resize(num_dofs);
763 2556451 : _ad_u_dot.resize(nqp);
764 :
765 2556451 : if (_time_integrator && _time_integrator->dt())
766 : {
767 13008833 : for (const auto i : make_range(num_dofs))
768 10474304 : _ad_dofs_dot[i] = _ad_dof_values[i];
769 13008833 : for (const auto i : make_range(num_dofs))
770 20947322 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
771 10474304 : _dof_indices[i],
772 10474304 : _need_ad_u_dotdot ? _ad_dofs_dotdot[i]
773 : : _ad_real_dummy);
774 :
775 11021901 : for (const auto qp : make_range(nqp))
776 8487372 : _ad_u_dot[qp] = _ad_zero;
777 13008833 : for (const auto i : make_range(num_dofs))
778 46114330 : for (const auto qp : make_range(nqp))
779 35640026 : _ad_u_dot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dot[i];
780 : }
781 : // We are too early in the setup to have a time integrator, so we are not really using the
782 : // AD-derivatives. We set the AD value of the derivatives to the nonAD value
783 21922 : else if (!_time_integrator)
784 : {
785 0 : for (const auto i : make_range(num_dofs))
786 0 : _ad_dofs_dot[i] = _dof_values_dot[i];
787 0 : for (const auto qp : make_range(nqp))
788 0 : _ad_u_dot[qp] = _u_dot[qp];
789 : }
790 : }
791 :
792 12321742 : if (_need_ad_u_dotdot)
793 : {
794 449 : _ad_u_dotdot.resize(nqp);
795 1807 : for (const auto qp : make_range(nqp))
796 1358 : _ad_u_dotdot[qp] = _ad_zero;
797 :
798 449 : if (_time_integrator && _time_integrator->dt())
799 1690 : for (const auto i : make_range(num_dofs))
800 6100 : for (const auto qp : make_range(nqp))
801 4814 : _ad_u_dotdot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dotdot[i];
802 45 : else if (!_time_integrator)
803 0 : for (const auto qp : make_range(nqp))
804 0 : _ad_u_dotdot[qp] = _u_dotdot[qp];
805 : }
806 :
807 12321742 : if (_need_ad_grad_u_dot)
808 : {
809 0 : _ad_grad_u_dot.resize(nqp);
810 :
811 0 : if (_time_integrator && _time_integrator->dt())
812 : {
813 0 : for (const auto qp : make_range(nqp))
814 0 : _ad_grad_u_dot[qp] = _ad_zero;
815 :
816 : // The latter check here is for handling the fact that we have not yet implemented
817 : // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that
818 : // situation we need to default to using the non-ad grad_phi
819 0 : if (_displaced && _current_ad_grad_phi)
820 0 : for (const auto i : make_range(num_dofs))
821 0 : for (const auto qp : make_range(nqp))
822 0 : _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_ad_grad_phi)[i][qp];
823 : else
824 0 : for (const auto i : make_range(num_dofs))
825 0 : for (const auto qp : make_range(nqp))
826 0 : _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_grad_phi)[i][qp];
827 : }
828 0 : else if (!_time_integrator)
829 0 : for (const auto qp : make_range(nqp))
830 0 : _ad_grad_u_dot[qp] = _grad_u_dot[qp];
831 : }
832 : }
833 28392888 : }
834 :
835 : template <>
836 : void
837 0 : MooseVariableData<RealEigenVector>::computeAD(const unsigned int /*num_dofs*/,
838 : const unsigned int /*nqp*/)
839 : {
840 0 : mooseError("AD for array variable has not been implemented");
841 : }
842 :
843 : template <typename OutputType>
844 : void
845 8654464 : MooseVariableData<OutputType>::setDofValue(const OutputData & value, unsigned int index)
846 : {
847 8654464 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
848 8654464 : dof_values[index] = value;
849 8654464 : _has_dof_values = true;
850 :
851 8654464 : auto & u = _vector_tag_u[_solution_tag];
852 8654464 : const auto nqps = u.size();
853 8654464 : const auto ndofs = dof_values.size();
854 45983266 : for (const auto qp : make_range(nqps))
855 37328802 : u[qp] *= 0.;
856 45983266 : for (const auto qp : make_range(nqps))
857 317274646 : for (const auto i : make_range(ndofs))
858 279945844 : u[qp] += (*_phi)[i][qp] * dof_values[i];
859 8654464 : }
860 :
861 : template <typename OutputType>
862 : void
863 651227 : MooseVariableData<OutputType>::setDofValues(const DenseVector<OutputData> & values)
864 : {
865 651227 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
866 3862142 : for (unsigned int i = 0; i < values.size(); i++)
867 3210915 : dof_values[i] = values(i);
868 :
869 651227 : _has_dof_values = true;
870 :
871 651227 : auto & u = _vector_tag_u[_solution_tag];
872 651227 : const auto nqps = u.size();
873 651227 : const auto ndofs = dof_values.size();
874 5572192 : for (const auto qp : make_range(nqps))
875 4920965 : u[qp] *= 0.;
876 5572192 : for (const auto qp : make_range(nqps))
877 36061775 : for (const auto i : make_range(ndofs))
878 31140810 : u[qp] += (*_phi)[i][qp] * dof_values[i];
879 651227 : }
880 :
881 : template <typename OutputType>
882 : void
883 61405632 : MooseVariableData<OutputType>::insertNodalValue(NumericVector<Number> & residual,
884 : const OutputData & v)
885 : {
886 61405632 : residual.set(_nodal_dof_index, v);
887 61405632 : }
888 :
889 : template <>
890 : void
891 313701 : MooseVariableData<RealEigenVector>::insertNodalValue(NumericVector<Number> & residual,
892 : const RealEigenVector & v)
893 : {
894 1116623 : for (const auto j : make_range(_count))
895 802922 : residual.set(_nodal_dof_index + j, v(j));
896 313701 : }
897 :
898 : template <typename OutputType>
899 : typename MooseVariableData<OutputType>::OutputData
900 121377 : MooseVariableData<OutputType>::getNodalValue(const Node & node, Moose::SolutionState state) const
901 : {
902 : mooseAssert(_subproblem.mesh().isSemiLocal(const_cast<Node *>(&node)), "Node is not Semilocal");
903 :
904 : // Make sure that the node has DOFs
905 : /* Note, this is a reproduction of an assert within libMesh::Node::dof_number, this is done to
906 : * produce a better error (see misc/check_error.node_value_off_block) */
907 : mooseAssert(node.n_dofs(_sys.number(), _var_num) > 0,
908 : "Node " << node.id() << " does not contain any dofs for the "
909 : << _sys.system().variable_name(_var_num) << " variable");
910 :
911 121377 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
912 :
913 121377 : switch (state)
914 : {
915 1345 : case Moose::Current:
916 1345 : return (*_sys.currentSolution())(dof);
917 :
918 120032 : case Moose::Old:
919 120032 : return _sys.solutionOld()(dof);
920 :
921 0 : case Moose::Older:
922 0 : return _sys.solutionOlder()(dof);
923 :
924 0 : default:
925 0 : mooseError("PreviousNL not currently supported for getNodalValue");
926 : }
927 : }
928 :
929 : template <>
930 : RealEigenVector
931 0 : MooseVariableData<RealEigenVector>::getNodalValue(const Node & node,
932 : Moose::SolutionState state) const
933 : {
934 : mooseAssert(_subproblem.mesh().isSemiLocal(const_cast<Node *>(&node)), "Node is not Semilocal");
935 :
936 : // Make sure that the node has DOFs
937 : /* Note, this is a reproduction of an assert within libMesh::Node::dof_number, this is done to
938 : * produce a better error (see misc/check_error.node_value_off_block) */
939 : mooseAssert(node.n_dofs(_sys.number(), _var_num) > 0,
940 : "Node " << node.id() << " does not contain any dofs for the "
941 : << _sys.system().variable_name(_var_num) << " variable");
942 :
943 0 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
944 :
945 0 : RealEigenVector v(_count);
946 0 : switch (state)
947 : {
948 0 : case Moose::Current:
949 0 : for (unsigned int i = 0; i < _count; ++i)
950 0 : v(i) = (*_sys.currentSolution())(dof++);
951 0 : break;
952 :
953 0 : case Moose::Old:
954 0 : for (unsigned int i = 0; i < _count; ++i)
955 0 : v(i) = _sys.solutionOld()(dof++);
956 0 : break;
957 :
958 0 : case Moose::Older:
959 0 : for (unsigned int i = 0; i < _count; ++i)
960 0 : v(i) = _sys.solutionOlder()(dof++);
961 0 : break;
962 :
963 0 : default:
964 0 : mooseError("PreviousNL not currently supported for getNodalValue");
965 : }
966 0 : return v;
967 0 : }
968 :
969 : template <typename OutputType>
970 : typename MooseVariableData<OutputType>::OutputData
971 79650 : MooseVariableData<OutputType>::getElementalValue(const Elem * elem,
972 : Moose::SolutionState state,
973 : unsigned int idx) const
974 : {
975 79650 : std::vector<dof_id_type> dof_indices;
976 79650 : _dof_map.dof_indices(elem, dof_indices, _var_num);
977 :
978 79650 : switch (state)
979 : {
980 78546 : case Moose::Current:
981 78546 : return (*_sys.currentSolution())(dof_indices[idx]);
982 :
983 552 : case Moose::Old:
984 552 : return _sys.solutionOld()(dof_indices[idx]);
985 :
986 552 : case Moose::Older:
987 552 : return _sys.solutionOlder()(dof_indices[idx]);
988 :
989 0 : default:
990 0 : mooseError("PreviousNL not currently supported for getElementalValue");
991 : }
992 79650 : }
993 :
994 : template <>
995 : RealEigenVector
996 0 : MooseVariableData<RealEigenVector>::getElementalValue(const Elem * elem,
997 : Moose::SolutionState state,
998 : unsigned int idx) const
999 : {
1000 0 : std::vector<dof_id_type> dof_indices;
1001 0 : _dof_map.dof_indices(elem, dof_indices, _var_num);
1002 :
1003 0 : dof_id_type dof = dof_indices[idx];
1004 :
1005 0 : RealEigenVector v(_count);
1006 :
1007 0 : switch (state)
1008 : {
1009 0 : case Moose::Current:
1010 0 : for (unsigned int i = 0; i < _count; ++i)
1011 0 : v(i) = (*_sys.currentSolution())(dof++);
1012 0 : break;
1013 :
1014 0 : case Moose::Old:
1015 0 : for (unsigned int i = 0; i < _count; ++i)
1016 0 : v(i) = _sys.solutionOld()(dof++);
1017 0 : break;
1018 :
1019 0 : case Moose::Older:
1020 0 : for (unsigned int i = 0; i < _count; ++i)
1021 0 : v(i) = _sys.solutionOlder()(dof++);
1022 0 : break;
1023 :
1024 0 : default:
1025 0 : mooseError("PreviousNL not currently supported for getElementalValue");
1026 : }
1027 0 : return v;
1028 0 : }
1029 :
1030 : template <typename OutputType>
1031 : void
1032 21514 : MooseVariableData<OutputType>::getDofIndices(const Elem * elem,
1033 : std::vector<dof_id_type> & dof_indices) const
1034 : {
1035 21514 : _dof_map.dof_indices(elem, dof_indices, _var_num);
1036 21514 : }
1037 :
1038 : template <typename OutputType>
1039 : void
1040 0 : MooseVariableData<OutputType>::addSolution(NumericVector<Number> & sol,
1041 : const DenseVector<Number> & v) const
1042 : {
1043 0 : sol.add_vector(v, _dof_indices);
1044 0 : }
1045 :
1046 : template <>
1047 : void
1048 3060 : MooseVariableData<RealEigenVector>::addSolution(NumericVector<Number> & sol,
1049 : const DenseVector<Number> & v) const
1050 : {
1051 3060 : unsigned int p = 0;
1052 3060 : const auto num_dofs = _dof_indices.size();
1053 9180 : for (const auto j : make_range(_count))
1054 : {
1055 6120 : unsigned int inc = (isNodal() ? j : j * num_dofs);
1056 30600 : for (const auto i : make_range(num_dofs))
1057 24480 : sol.add(_dof_indices[i] + inc, v(p++));
1058 : }
1059 3060 : }
1060 :
1061 : template <typename OutputType>
1062 : const typename MooseVariableData<OutputType>::DoFValue &
1063 354 : MooseVariableData<OutputType>::dofValuesDot() const
1064 : {
1065 354 : if (_sys.solutionUDot())
1066 : {
1067 354 : _need_dof_values_dot = true;
1068 354 : return _dof_values_dot;
1069 : }
1070 : else
1071 0 : mooseError("MooseVariableData: Time derivative of solution (`u_dot`) is not stored. Please set "
1072 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
1073 : }
1074 :
1075 : template <typename OutputType>
1076 : const typename MooseVariableData<OutputType>::DoFValue &
1077 4 : MooseVariableData<OutputType>::dofValuesDotDot() const
1078 : {
1079 4 : if (_sys.solutionUDotDot())
1080 : {
1081 4 : _need_dof_values_dotdot = true;
1082 4 : return _dof_values_dotdot;
1083 : }
1084 : else
1085 0 : mooseError("MooseVariableData: Second time derivative of solution (`u_dotdot`) is not stored. "
1086 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
1087 : "`u_dotdot`.");
1088 : }
1089 :
1090 : template <typename OutputType>
1091 : const typename MooseVariableData<OutputType>::DoFValue &
1092 0 : MooseVariableData<OutputType>::dofValuesDotOld() const
1093 : {
1094 0 : if (_sys.solutionUDotOld())
1095 : {
1096 0 : _need_dof_values_dot_old = true;
1097 0 : return _dof_values_dot_old;
1098 : }
1099 : else
1100 0 : mooseError("MooseVariableData: Old time derivative of solution (`u_dot_old`) is not stored. "
1101 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
1102 : "`u_dot_old`.");
1103 : }
1104 :
1105 : template <typename OutputType>
1106 : const typename MooseVariableData<OutputType>::DoFValue &
1107 0 : MooseVariableData<OutputType>::dofValuesDotDotOld() const
1108 : {
1109 0 : if (_sys.solutionUDotDotOld())
1110 : {
1111 0 : _need_dof_values_dotdot_old = true;
1112 0 : return _dof_values_dotdot_old;
1113 : }
1114 : else
1115 0 : mooseError("MooseVariableData: Old second time derivative of solution (`u_dotdot_old`) is not "
1116 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
1117 : "requesting `u_dotdot_old`.");
1118 : }
1119 :
1120 : template <typename OutputType>
1121 : const MooseArray<Number> &
1122 156 : MooseVariableData<OutputType>::dofValuesDuDotDu() const
1123 : {
1124 156 : _need_dof_du_dot_du = true;
1125 156 : return _dof_du_dot_du;
1126 : }
1127 :
1128 : template <typename OutputType>
1129 : const MooseArray<Number> &
1130 0 : MooseVariableData<OutputType>::dofValuesDuDotDotDu() const
1131 : {
1132 0 : _need_dof_du_dotdot_du = true;
1133 0 : return _dof_du_dotdot_du;
1134 : }
1135 :
1136 : template <typename OutputType>
1137 : void
1138 1776 : MooseVariableData<OutputType>::computeIncrementAtQps(const NumericVector<Number> & increment_vec)
1139 : {
1140 1776 : unsigned int nqp = _qrule->n_points();
1141 :
1142 1776 : _increment.resize(nqp);
1143 : // Compute the increment at each quadrature point
1144 1776 : unsigned int num_dofs = _dof_indices.size();
1145 15142 : for (const auto qp : make_range(nqp))
1146 : {
1147 13366 : _increment[qp] = 0.;
1148 130356 : for (const auto i : make_range(num_dofs))
1149 116990 : _increment[qp] += (*_phi)[i][qp] * increment_vec(_dof_indices[i]);
1150 : }
1151 1776 : }
1152 :
1153 : template <>
1154 : void
1155 0 : MooseVariableData<RealEigenVector>::computeIncrementAtQps(
1156 : const NumericVector<Number> & increment_vec)
1157 : {
1158 0 : unsigned int nqp = _qrule->n_points();
1159 :
1160 0 : _increment.resize(nqp);
1161 : // Compute the increment at each quadrature point
1162 0 : unsigned int num_dofs = _dof_indices.size();
1163 0 : if (isNodal())
1164 : {
1165 0 : for (const auto qp : make_range(nqp))
1166 : {
1167 0 : for (const auto i : make_range(num_dofs))
1168 0 : for (const auto j : make_range(_count))
1169 0 : _increment[qp](j) += (*_phi)[i][qp] * increment_vec(_dof_indices[i] + j);
1170 : }
1171 : }
1172 : else
1173 : {
1174 0 : for (const auto qp : make_range(nqp))
1175 : {
1176 0 : unsigned int n = 0;
1177 0 : for (const auto j : make_range(_count))
1178 0 : for (const auto i : make_range(num_dofs))
1179 : {
1180 0 : _increment[qp](j) += (*_phi)[i][qp] * increment_vec(_dof_indices[i] + n);
1181 0 : n += num_dofs;
1182 : }
1183 : }
1184 : }
1185 0 : }
1186 :
1187 : template <typename OutputType>
1188 : void
1189 7441 : MooseVariableData<OutputType>::computeIncrementAtNode(const NumericVector<Number> & increment_vec)
1190 : {
1191 7441 : if (!isNodal())
1192 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1193 :
1194 7441 : _increment.resize(1);
1195 :
1196 : // Compute the increment for the current DOF
1197 7441 : _increment[0] = increment_vec(_dof_indices[0]);
1198 7441 : }
1199 :
1200 : template <>
1201 : void
1202 0 : MooseVariableData<RealEigenVector>::computeIncrementAtNode(
1203 : const NumericVector<Number> & increment_vec)
1204 : {
1205 0 : if (!isNodal())
1206 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1207 :
1208 0 : _increment.resize(1);
1209 :
1210 : // Compute the increment for the current DOF
1211 0 : if (isNodal())
1212 0 : for (unsigned int j = 0; j < _count; j++)
1213 0 : _increment[0](j) = increment_vec(_dof_indices[0] + j);
1214 : else
1215 : {
1216 0 : unsigned int n = 0;
1217 0 : const auto n_dof_indices = _dof_indices.size();
1218 0 : for (const auto j : make_range(_count))
1219 : {
1220 0 : _increment[0](j) = increment_vec(_dof_indices[0] + n);
1221 0 : n += n_dof_indices;
1222 : }
1223 : }
1224 0 : }
1225 :
1226 : template <typename OutputType>
1227 : const OutputType &
1228 0 : MooseVariableData<OutputType>::nodalValueDot() const
1229 : {
1230 0 : if (isNodal())
1231 : {
1232 0 : if (_sys.solutionUDot())
1233 : {
1234 0 : _need_dof_values_dot = true;
1235 0 : return _nodal_value_dot;
1236 : }
1237 : else
1238 0 : mooseError(
1239 : "MooseVariableData: Time derivative of solution (`u_dot`) is not stored. Please set "
1240 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
1241 : }
1242 : else
1243 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1244 0 : var().name(),
1245 : "' is not nodal.");
1246 : }
1247 :
1248 : template <typename OutputType>
1249 : const OutputType &
1250 0 : MooseVariableData<OutputType>::nodalValueDotDot() const
1251 : {
1252 0 : if (isNodal())
1253 : {
1254 0 : if (_sys.solutionUDotDot())
1255 : {
1256 0 : _need_dof_values_dotdot = true;
1257 0 : return _nodal_value_dotdot;
1258 : }
1259 : else
1260 0 : mooseError(
1261 : "MooseVariableData: Second time derivative of solution (`u_dotdot`) is not stored. "
1262 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
1263 : "`u_dotdot`.");
1264 : }
1265 : else
1266 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1267 0 : var().name(),
1268 : "' is not nodal.");
1269 : }
1270 :
1271 : template <typename OutputType>
1272 : const OutputType &
1273 0 : MooseVariableData<OutputType>::nodalValueDotOld() const
1274 : {
1275 0 : if (isNodal())
1276 : {
1277 0 : if (_sys.solutionUDotOld())
1278 : {
1279 0 : _need_dof_values_dot_old = true;
1280 0 : return _nodal_value_dot_old;
1281 : }
1282 : else
1283 0 : mooseError("MooseVariableData: Old time derivative of solution (`u_dot_old`) is not stored. "
1284 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
1285 : "`u_dot_old`.");
1286 : }
1287 : else
1288 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1289 0 : var().name(),
1290 : "' is not nodal.");
1291 : }
1292 :
1293 : template <typename OutputType>
1294 : const OutputType &
1295 0 : MooseVariableData<OutputType>::nodalValueDotDotOld() const
1296 : {
1297 0 : if (isNodal())
1298 : {
1299 0 : if (_sys.solutionUDotDotOld())
1300 : {
1301 0 : _need_dof_values_dotdot_old = true;
1302 0 : return _nodal_value_dotdot_old;
1303 : }
1304 : else
1305 0 : mooseError(
1306 : "MooseVariableData: Old second time derivative of solution (`u_dotdot_old`) is not "
1307 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
1308 : "requesting `u_dotdot_old`.");
1309 : }
1310 : else
1311 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1312 0 : var().name(),
1313 : "' is not nodal.");
1314 : }
1315 :
1316 : template <typename OutputType>
1317 : void
1318 276007858 : MooseVariableData<OutputType>::computeNodalValues()
1319 : {
1320 276007858 : if (_has_dof_indices)
1321 : {
1322 276007618 : fetchDoFValues();
1323 276007618 : assignNodalValue();
1324 :
1325 276007618 : if (_need_ad)
1326 3423463 : fetchADNodalValues();
1327 : }
1328 : else
1329 240 : zeroSizeDofValues();
1330 276007858 : }
1331 :
1332 : template <typename OutputType>
1333 : void
1334 3423463 : MooseVariableData<OutputType>::fetchADNodalValues()
1335 : {
1336 3423463 : auto n = _dof_indices.size();
1337 : libmesh_assert(n);
1338 3423463 : _ad_dof_values.resize(n);
1339 :
1340 3423463 : if (_need_ad_u_dot)
1341 581910 : _ad_dofs_dot.resize(n);
1342 3423463 : if (_need_ad_u_dotdot)
1343 956 : _ad_dofs_dotdot.resize(n);
1344 :
1345 3423463 : const bool do_derivatives =
1346 3423463 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
1347 :
1348 6992750 : for (decltype(n) i = 0; i < n; ++i)
1349 : {
1350 3569287 : _ad_dof_values[i] = _vector_tags_dof_u[_solution_tag][i];
1351 3569287 : if (do_derivatives)
1352 1436026 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
1353 3569287 : assignADNodalValue(_ad_dof_values[i], i);
1354 :
1355 3569287 : if (_need_ad_u_dot)
1356 : {
1357 581910 : if (_time_integrator && _time_integrator->dt())
1358 : {
1359 581910 : _ad_dofs_dot[i] = _ad_dof_values[i];
1360 1162864 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
1361 581910 : _dof_indices[i],
1362 581910 : _need_ad_u_dotdot ? _ad_dofs_dotdot[i]
1363 : : _ad_real_dummy);
1364 : }
1365 : // Executing something with a time derivative at initial should not put a NaN
1366 0 : else if (_time_integrator && !_time_integrator->dt())
1367 0 : _ad_dofs_dot[i] = 0.;
1368 : else
1369 0 : mooseError("AD nodal time derivatives not implemented for variables without a time "
1370 : "integrator (auxiliary variables)");
1371 : }
1372 : }
1373 3423463 : }
1374 :
1375 : template <>
1376 : void
1377 0 : MooseVariableData<RealEigenVector>::fetchADNodalValues()
1378 : {
1379 0 : mooseError("I do not know how to support AD with array variables");
1380 : }
1381 :
1382 : template <>
1383 : void
1384 3275098 : MooseVariableData<Real>::assignADNodalValue(const ADReal & value, const unsigned int &)
1385 : {
1386 3275098 : _ad_nodal_value = value;
1387 3275098 : }
1388 :
1389 : template <>
1390 : void
1391 294189 : MooseVariableData<RealVectorValue>::assignADNodalValue(const ADReal & value,
1392 : const unsigned int & component)
1393 : {
1394 294189 : _ad_nodal_value(component) = value;
1395 294189 : }
1396 :
1397 : template <typename OutputType>
1398 : void
1399 2277553 : MooseVariableData<OutputType>::prepareIC()
1400 : {
1401 2277553 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1402 2277553 : _vector_tags_dof_u[_solution_tag].resize(_dof_indices.size());
1403 :
1404 2277553 : unsigned int nqp = _qrule->n_points();
1405 2277553 : _vector_tag_u[_solution_tag].resize(nqp);
1406 2277553 : }
1407 :
1408 : template <typename OutputType>
1409 : void
1410 539683531 : MooseVariableData<OutputType>::prepare()
1411 : {
1412 539683531 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1413 539683531 : _has_dof_values = false;
1414 :
1415 : // FIXME: remove this when the Richard's module is migrated to use the new NodalCoupleable
1416 : // interface.
1417 539683531 : _has_dof_indices = _dof_indices.size();
1418 539683531 : }
1419 :
1420 : template <typename OutputType>
1421 : void
1422 314997342 : MooseVariableData<OutputType>::reinitNode()
1423 : {
1424 314997342 : if (std::size_t n_dofs = _node->n_dofs(_sys.number(), _var_num))
1425 : {
1426 276000667 : _dof_indices.resize(n_dofs);
1427 552967232 : for (std::size_t i = 0; i < n_dofs; ++i)
1428 276966565 : _dof_indices[i] = _node->dof_number(_sys.number(), _var_num, i);
1429 : // For standard variables. _nodal_dof_index is retrieved by nodalDofIndex() which is used in
1430 : // NodalBC for example
1431 276000667 : _nodal_dof_index = _dof_indices[0];
1432 276000667 : _has_dof_indices = true;
1433 : }
1434 : else
1435 : {
1436 38996675 : _dof_indices.clear(); // Clear these so Assembly doesn't think there's dofs here
1437 38996675 : _has_dof_indices = false;
1438 : }
1439 314997342 : }
1440 :
1441 : template <typename OutputType>
1442 : void
1443 167971304 : MooseVariableData<OutputType>::reinitAux()
1444 : {
1445 : /* FIXME: this method is only for elemental auxiliary variables, so
1446 : * we may want to rename it */
1447 167971304 : if (_elem)
1448 : {
1449 166897363 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1450 166897363 : if (_elem->n_dofs(_sys.number(), _var_num) > 0)
1451 : {
1452 165919712 : _nodal_dof_index = _dof_indices[0];
1453 :
1454 165919712 : fetchDoFValues();
1455 :
1456 165919712 : const auto num_dofs = _dof_indices.size();
1457 1200500910 : for (auto & dof_u : _vector_tags_dof_u)
1458 1034581198 : dof_u.resize(num_dofs);
1459 :
1460 388129798 : for (auto & dof_u : _matrix_tags_dof_u)
1461 222210086 : dof_u.resize(num_dofs);
1462 :
1463 165919712 : _has_dof_indices = true;
1464 : }
1465 : else
1466 977651 : _has_dof_indices = false;
1467 : }
1468 : else
1469 1073941 : _has_dof_indices = false;
1470 167971304 : }
1471 :
1472 : template <typename OutputType>
1473 : void
1474 7866 : MooseVariableData<OutputType>::reinitNodes(const std::vector<dof_id_type> & nodes)
1475 : {
1476 7866 : _dof_indices.clear();
1477 42615 : for (const auto & node_id : nodes)
1478 : {
1479 34749 : auto && nd = _subproblem.mesh().getMesh().query_node_ptr(node_id);
1480 34749 : if (nd && (_subproblem.mesh().isSemiLocal(const_cast<Node *>(nd))))
1481 : {
1482 34683 : if (nd->n_dofs(_sys.number(), _var_num) > 0)
1483 : {
1484 33492 : dof_id_type dof = nd->dof_number(_sys.number(), _var_num, 0);
1485 33492 : _dof_indices.push_back(dof);
1486 : }
1487 : }
1488 : }
1489 :
1490 7866 : if (!_dof_indices.empty())
1491 7626 : _has_dof_indices = true;
1492 : else
1493 240 : _has_dof_indices = false;
1494 7866 : }
1495 :
1496 : template class MooseVariableData<Real>;
1497 : template class MooseVariableData<RealVectorValue>;
1498 : template class MooseVariableData<RealEigenVector>;
|
