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 495537 : 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 991074 : _fe_type(var.feType()),
39 495537 : _var_num(var.number()),
40 495537 : _assembly(_subproblem.assembly(_tid, var.kind() == Moose::VAR_SOLVER ? sys.number() : 0)),
41 495537 : _element_type(element_type),
42 495537 : _ad_zero(0),
43 495537 : _need_ad_u_dot(false),
44 495537 : _need_ad_u_dotdot(false),
45 495537 : _need_second(false),
46 495537 : _need_second_old(false),
47 495537 : _need_second_older(false),
48 495537 : _need_second_previous_nl(false),
49 495537 : _need_curl(false),
50 495537 : _need_curl_old(false),
51 495537 : _need_curl_older(false),
52 495537 : _need_div(false),
53 495537 : _need_div_old(false),
54 495537 : _need_div_older(false),
55 495537 : _need_ad(false),
56 495537 : _need_ad_u(false),
57 495537 : _need_ad_grad_u(false),
58 495537 : _need_ad_grad_u_dot(false),
59 495537 : _need_ad_second_u(false),
60 495537 : _need_ad_curl_u(false),
61 495537 : _has_dof_indices(false),
62 495537 : _qrule(qrule_in),
63 495537 : _qrule_face(qrule_face_in),
64 495537 : _use_dual(var.useDual()),
65 495537 : _second_phi_assembly_method(nullptr),
66 495537 : _second_phi_face_assembly_method(nullptr),
67 495537 : _curl_phi_assembly_method(nullptr),
68 495537 : _curl_phi_face_assembly_method(nullptr),
69 495537 : _div_phi_assembly_method(nullptr),
70 495537 : _div_phi_face_assembly_method(nullptr),
71 495537 : _ad_grad_phi_assembly_method(nullptr),
72 495537 : _ad_grad_phi_face_assembly_method(nullptr),
73 495537 : _time_integrator(nullptr),
74 495537 : _node(node),
75 495537 : _elem(elem),
76 495537 : _displaced(dynamic_cast<const DisplacedSystem *>(&_sys) ? true : false),
77 2477685 : _current_side(_assembly.side())
78 : {
79 495537 : _continuity = FEInterface::get_continuity(_fe_type);
80 :
81 495537 : _is_nodal = (_continuity == C_ZERO || _continuity == C_ONE);
82 :
83 495537 : _time_integrator = _sys.queryTimeIntegrator(_var_num);
84 :
85 : // Initialize AD zero with zero derivatives
86 495537 : const auto old_do = ADReal::do_derivatives;
87 495537 : ADReal::do_derivatives = true;
88 495537 : _ad_zero = 0.;
89 495537 : ADReal::do_derivatives = old_do;
90 :
91 495537 : switch (_element_type)
92 : {
93 165179 : case Moose::ElementType::Element:
94 : {
95 165179 : _phi_assembly_method = &Assembly::fePhi<OutputShape>;
96 165179 : _phi_face_assembly_method = &Assembly::fePhiFace<OutputShape>;
97 165179 : _grad_phi_assembly_method = &Assembly::feGradPhi<OutputShape>;
98 165179 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFace<OutputShape>;
99 165179 : _second_phi_assembly_method = &Assembly::feSecondPhi<OutputShape>;
100 165179 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFace<OutputShape>;
101 165179 : _curl_phi_assembly_method = &Assembly::feCurlPhi<OutputShape>;
102 165179 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFace<OutputShape>;
103 165179 : _div_phi_assembly_method = &Assembly::feDivPhi<OutputShape>;
104 165179 : _div_phi_face_assembly_method = &Assembly::feDivPhiFace<OutputShape>;
105 165179 : _ad_grad_phi_assembly_method = &Assembly::feADGradPhi<OutputShape>;
106 165179 : _ad_grad_phi_face_assembly_method = &Assembly::feADGradPhiFace<OutputShape>;
107 :
108 165179 : _ad_grad_phi = &_ad_grad_phi_assembly_method(_assembly, _fe_type);
109 165179 : _ad_grad_phi_face = &_ad_grad_phi_face_assembly_method(_assembly, _fe_type);
110 165179 : break;
111 : }
112 165179 : case Moose::ElementType::Neighbor:
113 : {
114 165179 : _phi_assembly_method = &Assembly::fePhiNeighbor<OutputShape>;
115 165179 : _phi_face_assembly_method = &Assembly::fePhiFaceNeighbor<OutputShape>;
116 165179 : _grad_phi_assembly_method = &Assembly::feGradPhiNeighbor<OutputShape>;
117 165179 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFaceNeighbor<OutputShape>;
118 165179 : _second_phi_assembly_method = &Assembly::feSecondPhiNeighbor<OutputShape>;
119 165179 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFaceNeighbor<OutputShape>;
120 165179 : _curl_phi_assembly_method = &Assembly::feCurlPhiNeighbor<OutputShape>;
121 165179 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFaceNeighbor<OutputShape>;
122 165179 : _div_phi_assembly_method = &Assembly::feDivPhiNeighbor<OutputShape>;
123 165179 : _div_phi_face_assembly_method = &Assembly::feDivPhiFaceNeighbor<OutputShape>;
124 :
125 165179 : _ad_grad_phi = nullptr;
126 165179 : _ad_grad_phi_face = nullptr;
127 165179 : break;
128 : }
129 165179 : case Moose::ElementType::Lower:
130 : {
131 165179 : if (_use_dual)
132 : {
133 135 : _phi_assembly_method = &Assembly::feDualPhiLower<OutputType>;
134 135 : _phi_face_assembly_method = &Assembly::feDualPhiLower<OutputType>; // Place holder
135 135 : _grad_phi_assembly_method = &Assembly::feGradDualPhiLower<OutputType>;
136 135 : _grad_phi_face_assembly_method = &Assembly::feGradDualPhiLower<OutputType>; // Place holder
137 : }
138 : else
139 : {
140 165044 : _phi_assembly_method = &Assembly::fePhiLower<OutputType>;
141 165044 : _phi_face_assembly_method = &Assembly::fePhiLower<OutputType>; // Place holder
142 165044 : _grad_phi_assembly_method = &Assembly::feGradPhiLower<OutputType>;
143 165044 : _grad_phi_face_assembly_method = &Assembly::feGradPhiLower<OutputType>; // Place holder
144 : }
145 :
146 165179 : _ad_grad_phi = nullptr;
147 165179 : _ad_grad_phi_face = nullptr;
148 165179 : break;
149 : }
150 : }
151 495537 : _phi = &_phi_assembly_method(_assembly, _fe_type);
152 495537 : _phi_face = &_phi_face_assembly_method(_assembly, _fe_type);
153 495537 : _grad_phi = &_grad_phi_assembly_method(_assembly, _fe_type);
154 495537 : _grad_phi_face = &_grad_phi_face_assembly_method(_assembly, _fe_type);
155 495537 : }
156 :
157 : template <typename OutputType>
158 : void
159 705589435 : MooseVariableData<OutputType>::setGeometry(Moose::GeometryType gm_type)
160 : {
161 705589435 : switch (gm_type)
162 : {
163 662981011 : case Moose::Volume:
164 : {
165 662981011 : _current_qrule = _qrule;
166 662981011 : _current_phi = _phi;
167 662981011 : _current_grad_phi = _grad_phi;
168 662981011 : _current_second_phi = _second_phi;
169 662981011 : _current_curl_phi = _curl_phi;
170 662981011 : _current_div_phi = _div_phi;
171 662981011 : _current_ad_grad_phi = _ad_grad_phi;
172 662981011 : break;
173 : }
174 42608424 : case Moose::Face:
175 : {
176 42608424 : _current_qrule = _qrule_face;
177 42608424 : _current_phi = _phi_face;
178 42608424 : _current_grad_phi = _grad_phi_face;
179 42608424 : _current_second_phi = _second_phi_face;
180 42608424 : _current_curl_phi = _curl_phi_face;
181 42608424 : _current_div_phi = _div_phi_face;
182 42608424 : _current_ad_grad_phi = _ad_grad_phi_face;
183 42608424 : break;
184 : }
185 : }
186 705589435 : }
187 :
188 : template <typename OutputType>
189 : const typename MooseVariableData<OutputType>::FieldVariableValue &
190 15965 : MooseVariableData<OutputType>::uDot() const
191 : {
192 15965 : if (_sys.solutionUDot())
193 : {
194 15965 : _need_u_dot = true;
195 15965 : 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 204 : MooseVariableData<OutputType>::uDotDot() const
205 : {
206 204 : if (_sys.solutionUDotDot())
207 : {
208 204 : _need_u_dotdot = true;
209 204 : 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 37 : MooseVariableData<OutputType>::gradSlnDot() const
250 : {
251 37 : if (_sys.solutionUDot())
252 : {
253 37 : _need_grad_dot = true;
254 37 : 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 187 : MooseVariableData<OutputType>::secondSln(Moose::SolutionState state) const
279 : {
280 187 : secondPhi();
281 187 : secondPhiFace();
282 187 : switch (state)
283 : {
284 187 : case Moose::Current:
285 : {
286 187 : _need_second = true;
287 187 : 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 159 : MooseVariableData<OutputType>::curlSln(Moose::SolutionState state) const
318 : {
319 159 : curlPhi();
320 159 : curlPhiFace();
321 159 : switch (state)
322 : {
323 159 : case Moose::Current:
324 : {
325 159 : _need_curl = true;
326 159 : 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 312 : MooseVariableData<OutputType>::divSln(Moose::SolutionState state) const
349 : {
350 312 : divPhi();
351 312 : divPhiFace();
352 312 : switch (state)
353 : {
354 312 : case Moose::Current:
355 : {
356 312 : _need_div = true;
357 312 : 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 26886 : MooseVariableData<OutputType>::secondPhi() const
380 : {
381 26886 : _second_phi = &_second_phi_assembly_method(_assembly, _fe_type);
382 26886 : return *_second_phi;
383 : }
384 :
385 : template <typename OutputType>
386 : const typename MooseVariableData<OutputType>::FieldVariablePhiSecond &
387 5408 : MooseVariableData<OutputType>::secondPhiFace() const
388 : {
389 5408 : _second_phi_face = &_second_phi_face_assembly_method(_assembly, _fe_type);
390 5408 : return *_second_phi_face;
391 : }
392 :
393 : template <typename OutputType>
394 : const typename MooseVariableData<OutputType>::FieldVariablePhiCurl &
395 57245 : MooseVariableData<OutputType>::curlPhi() const
396 : {
397 57245 : _curl_phi = &_curl_phi_assembly_method(_assembly, _fe_type);
398 57245 : return *_curl_phi;
399 : }
400 :
401 : template <typename OutputType>
402 : const typename MooseVariableData<OutputType>::FieldVariablePhiCurl &
403 209 : MooseVariableData<OutputType>::curlPhiFace() const
404 : {
405 209 : _curl_phi_face = &_curl_phi_face_assembly_method(_assembly, _fe_type);
406 209 : return *_curl_phi_face;
407 : }
408 :
409 : template <typename OutputType>
410 : const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
411 469068 : MooseVariableData<OutputType>::divPhi() const
412 : {
413 469068 : _div_phi = &_div_phi_assembly_method(_assembly, _fe_type);
414 469068 : return *_div_phi;
415 : }
416 :
417 : template <typename OutputType>
418 : const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
419 312 : MooseVariableData<OutputType>::divPhiFace() const
420 : {
421 312 : _div_phi_face = &_div_phi_face_assembly_method(_assembly, _fe_type);
422 312 : return *_div_phi_face;
423 : }
424 :
425 : template <typename OutputType>
426 : template <bool monomial>
427 : void
428 704878519 : MooseVariableData<OutputType>::computeValuesInternal()
429 : {
430 704878519 : const auto num_dofs = _dof_indices.size();
431 704878519 : if (num_dofs > 0)
432 615556963 : fetchDoFValues();
433 :
434 704878519 : const bool is_transient = _subproblem.isTransient();
435 704878519 : const auto nqp = _current_qrule->n_points();
436 704878519 : 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 4039217 : if (_qrule == _current_qrule)
442 : {
443 3825363 : _mapped_grad_phi.resize(num_dofs);
444 18115121 : for (const auto i : make_range(num_dofs))
445 : {
446 14289758 : _mapped_grad_phi[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
447 95693161 : for (const auto qp : make_range(nqp))
448 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
449 162806806 : new (&_mapped_grad_phi[i][qp])
450 81403403 : Eigen::Map<RealDIMValue>(const_cast<Real *>(&(*_current_grad_phi)[i][qp](0)));
451 : }
452 : }
453 : else
454 : {
455 213854 : _mapped_grad_phi_face.resize(num_dofs);
456 854266 : for (const auto i : make_range(num_dofs))
457 : {
458 640412 : _mapped_grad_phi_face[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
459 2261916 : for (const auto qp : make_range(nqp))
460 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
461 3243008 : new (&_mapped_grad_phi_face[i][qp])
462 1621504 : 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 2067126397 : const auto fill = [this, nqp, num_dofs](auto & dest, const auto & phi, const auto & dof_values)
478 : {
479 : if constexpr (monomial)
480 158395425 : libmesh_ignore(num_dofs);
481 :
482 : // Deduce OutputType
483 1362247878 : constexpr bool is_real = std::is_same_v<OutputType, Real>;
484 1362247878 : constexpr bool is_real_vector = std::is_same_v<OutputType, RealVectorValue>;
485 1362247878 : 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 1353979693 : 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 1362247878 : constexpr bool is_value = std::is_same_v<dest_array_type, OutputType>;
495 1362247878 : constexpr bool is_gradient = std::is_same_v<dest_array_type, OutputGradient>;
496 1362247878 : constexpr bool is_second = std::is_same_v<dest_array_type, OutputSecond>;
497 1362247878 : 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 18023583441 : const auto set_zero = [this, &dest](const auto qp)
503 : {
504 : if constexpr (!is_eigen)
505 5514928353 : libmesh_ignore(this);
506 :
507 : if constexpr (is_real || is_real_vector)
508 5514928353 : dest[qp] = 0;
509 : else if constexpr (is_eigen)
510 : {
511 : if constexpr (is_value)
512 27141450 : dest[qp].setZero(this->_count);
513 : else if constexpr (is_gradient)
514 11708718 : 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 >10138*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 15792145361 : dest[qp] += phi[i][qp] * dof_values[i];
531 : else if constexpr (is_gradient || is_second)
532 9197620972 : 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 265652364 : for (const auto d : make_range(Moose::dim))
541 199239273 : 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 1362247878 : 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 158395425 : set_zero(0);
563 158395425 : accumulate(0, 0);
564 708341111 : for (unsigned int qp = 1; qp < nqp; ++qp)
565 549945686 : dest[qp] = dest[0];
566 : }
567 : // Non monomial case
568 : else
569 : {
570 6599235549 : for (const auto qp : make_range(nqp))
571 5395383096 : set_zero(qp);
572 6249472193 : for (const auto i : make_range(num_dofs))
573 29943403739 : for (const auto qp : make_range(nqp))
574 24897783999 : accumulate(i, qp);
575 : }
576 : };
577 :
578 : // Curl
579 704878519 : if (_need_curl)
580 71020 : fill(_curl_u, *_current_curl_phi, _vector_tags_dof_u[_solution_tag]);
581 704878519 : if (_need_curl_old)
582 0 : fill(_curl_u_old, *_current_curl_phi, _vector_tags_dof_u[_old_solution_tag]);
583 :
584 : // Div
585 704878519 : if (_need_div)
586 352128 : fill(_div_u, *_current_div_phi, _vector_tags_dof_u[_solution_tag]);
587 704878519 : 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 704878519 : if (_need_second)
592 83053 : fill(_second_u, *_current_second_phi, _vector_tags_dof_u[_solution_tag]);
593 704878519 : 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 1418722938 : for (auto tag : _required_vector_tags)
599 : {
600 713844419 : if (_need_vector_tag_u[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed())
601 711383905 : fill(_vector_tag_u[tag], *_current_phi, _vector_tags_dof_u[tag]);
602 713844419 : if (_need_vector_tag_grad[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed())
603 406690873 : fill(_vector_tag_grad[tag], *_current_grad_phi, _vector_tags_dof_u[tag]);
604 : }
605 :
606 : // Matrix tags
607 704878903 : for (auto tag : active_coupleable_matrix_tags)
608 384 : if (_need_matrix_tag_u[tag])
609 128 : fill(_matrix_tag_u[tag], *_current_phi, _matrix_tags_dof_u[tag]);
610 :
611 : // Derivatives and old values
612 704878519 : if (is_transient)
613 : {
614 526179330 : if (_need_second_old)
615 0 : fill(_second_u_old, *_current_second_phi, _vector_tags_dof_u[_old_solution_tag]);
616 526179330 : if (_need_second_older)
617 0 : fill(_second_u_older, *_current_second_phi, _vector_tags_dof_u[_older_solution_tag]);
618 526179330 : if (_need_u_dot)
619 243535210 : fill(_u_dot, *_current_phi, _dof_values_dot);
620 526179330 : if (_need_u_dotdot)
621 6206 : fill(_u_dotdot, *_current_phi, _dof_values_dotdot);
622 526179330 : if (_need_u_dot_old)
623 0 : fill(_u_dot_old, *_current_phi, _dof_values_dot_old);
624 526179330 : if (_need_u_dotdot_old)
625 0 : fill(_u_dotdot_old, *_current_phi, _dof_values_dotdot_old);
626 :
627 526179330 : if (_need_du_dot_du)
628 : {
629 243657082 : _du_dot_du.resize(nqp);
630 1297579381 : for (const auto qp : make_range(nqp))
631 1053922299 : _du_dot_du[qp] = 0.;
632 1298837156 : for (const auto i : make_range(num_dofs))
633 6048738460 : for (const auto qp : make_range(nqp))
634 4993558386 : _du_dot_du[qp] = _dof_du_dot_du[i];
635 : }
636 526179330 : if (_need_du_dotdot_du)
637 : {
638 5062 : _du_dotdot_du.resize(nqp);
639 24140 : for (const auto qp : make_range(nqp))
640 19078 : _du_dotdot_du[qp] = 0.;
641 42828 : for (const auto i : make_range(num_dofs))
642 187660 : for (const auto qp : make_range(nqp))
643 149894 : _du_dotdot_du[qp] = _dof_du_dotdot_du[i];
644 : }
645 :
646 526179330 : if (_need_grad_dot)
647 125355 : fill(_grad_u_dot, *_current_grad_phi, _dof_values_dot);
648 526179330 : 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 700839302 : if (_need_ad)
655 25227938 : computeAD(num_dofs, nqp);
656 704878519 : }
657 :
658 : template <typename OutputType>
659 : void
660 556415243 : MooseVariableData<OutputType>::computeValues()
661 : {
662 556415243 : computeValuesInternal</* monomial = */ false>();
663 556415243 : }
664 :
665 : template <typename OutputType>
666 : void
667 149247540 : MooseVariableData<OutputType>::computeMonomialValues()
668 : {
669 149247540 : if (_dof_indices.size() == 0)
670 710916 : return;
671 :
672 : // Monomial optimizations are not appropriate after p-refinement
673 148536624 : if (_elem->p_level())
674 73348 : computeValues();
675 : else
676 148463276 : computeValuesInternal</* monomial = */ true>();
677 : }
678 :
679 : template <typename OutputType>
680 : void
681 25227938 : MooseVariableData<OutputType>::computeAD(const unsigned int num_dofs, const unsigned int nqp)
682 : {
683 25227938 : const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys);
684 :
685 25227938 : _ad_dof_values.resize(num_dofs);
686 139090169 : for (const auto i : make_range(num_dofs))
687 113862231 : _ad_dof_values[i] = (*_sys.currentSolution())(_dof_indices[i]);
688 : // NOTE! You have to do this AFTER setting the value!
689 25227938 : if (do_derivatives)
690 29033695 : for (const auto i : make_range(num_dofs))
691 24198936 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
692 :
693 25227938 : if (_need_ad_u)
694 : {
695 23903797 : _ad_u.resize(nqp);
696 104141594 : for (const auto qp : make_range(nqp))
697 80237797 : _ad_u[qp] = _ad_zero;
698 :
699 133501630 : for (const auto i : make_range(num_dofs))
700 533194000 : for (const auto qp : make_range(nqp))
701 423596167 : _ad_u[qp] += _ad_dof_values[i] * (*_current_phi)[i][qp];
702 : }
703 :
704 25227938 : if (_need_ad_grad_u)
705 : {
706 18614789 : _ad_grad_u.resize(nqp);
707 84159449 : for (const auto qp : make_range(nqp))
708 65544660 : _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 18614789 : if (_displaced && _current_ad_grad_phi)
714 2952192 : for (const auto i : make_range(num_dofs))
715 16241820 : for (const auto qp : make_range(nqp))
716 14361702 : _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_ad_grad_phi)[i][qp];
717 : else
718 93658364 : for (const auto i : make_range(num_dofs))
719 398310926 : for (const auto qp : make_range(nqp))
720 322731314 : _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_grad_phi)[i][qp];
721 : }
722 :
723 25227938 : 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 25227938 : if (_need_ad_curl_u)
737 : {
738 2858120 : _ad_curl_u.resize(nqp);
739 13362760 : for (const auto qp : make_range(nqp))
740 10504640 : _ad_curl_u[qp] = _ad_zero;
741 :
742 14359400 : for (const auto i : make_range(num_dofs))
743 53955040 : 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 42453760 : _ad_curl_u[qp] += _ad_dof_values[i] * (*_current_curl_phi)[i][qp];
750 : }
751 : }
752 :
753 25227938 : bool is_transient = _subproblem.isTransient();
754 25227938 : if (is_transient)
755 : {
756 10331616 : if (_need_ad_u_dot)
757 : {
758 2359413 : _ad_dofs_dot.resize(num_dofs);
759 2359413 : if (_need_ad_u_dotdot)
760 394 : _ad_dofs_dotdot.resize(num_dofs);
761 2359413 : _ad_u_dot.resize(nqp);
762 10198771 : for (const auto qp : make_range(nqp))
763 7839358 : _ad_u_dot[qp] = _ad_zero;
764 :
765 2359413 : if (_time_integrator && _time_integrator->dt())
766 : {
767 12021220 : for (const auto i : make_range(num_dofs))
768 9680378 : _ad_dofs_dot[i] = _ad_dof_values[i];
769 12021220 : for (const auto i : make_range(num_dofs))
770 19359634 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
771 9680378 : _dof_indices[i],
772 9680378 : _need_ad_u_dotdot ? _ad_dofs_dotdot[i]
773 : : _ad_real_dummy);
774 :
775 12021220 : for (const auto i : make_range(num_dofs))
776 42326164 : for (const auto qp : make_range(nqp))
777 32645786 : _ad_u_dot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dot[i];
778 : }
779 18571 : 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 10331616 : if (_need_ad_u_dotdot)
789 : {
790 394 : _ad_u_dotdot.resize(nqp);
791 1580 : for (const auto qp : make_range(nqp))
792 1186 : _ad_u_dotdot[qp] = _ad_zero;
793 :
794 394 : if (_time_integrator && _time_integrator->dt())
795 1476 : for (const auto i : make_range(num_dofs))
796 5316 : for (const auto qp : make_range(nqp))
797 4194 : _ad_u_dotdot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dotdot[i];
798 40 : 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 10331616 : 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 25227938 : }
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 7623827 : MooseVariableData<OutputType>::setDofValue(const OutputData & value, unsigned int index)
841 : {
842 7623827 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
843 7623827 : dof_values[index] = value;
844 7623827 : _has_dof_values = true;
845 :
846 7623827 : auto & u = _vector_tag_u[_solution_tag];
847 7623827 : const auto nqps = u.size();
848 7623827 : const auto ndofs = dof_values.size();
849 40712458 : for (const auto qp : make_range(nqps))
850 33088631 : u[qp] *= 0.;
851 40712458 : for (const auto qp : make_range(nqps))
852 282386514 : for (const auto i : make_range(ndofs))
853 249297883 : u[qp] += (*_phi)[i][qp] * dof_values[i];
854 7623827 : }
855 :
856 : template <typename OutputType>
857 : void
858 573480 : MooseVariableData<OutputType>::setDofValues(const DenseVector<OutputData> & values)
859 : {
860 573480 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
861 3399546 : for (unsigned int i = 0; i < values.size(); i++)
862 2826066 : dof_values[i] = values(i);
863 :
864 573480 : _has_dof_values = true;
865 :
866 573480 : auto & u = _vector_tag_u[_solution_tag];
867 573480 : const auto nqps = u.size();
868 573480 : const auto ndofs = dof_values.size();
869 4921392 : for (const auto qp : make_range(nqps))
870 4347912 : u[qp] *= 0.;
871 4921392 : for (const auto qp : make_range(nqps))
872 31890188 : for (const auto i : make_range(ndofs))
873 27542276 : u[qp] += (*_phi)[i][qp] * dof_values[i];
874 573480 : }
875 :
876 : template <typename OutputType>
877 : void
878 54525825 : MooseVariableData<OutputType>::insertNodalValue(NumericVector<Number> & residual,
879 : const OutputData & v)
880 : {
881 54525825 : residual.set(_nodal_dof_index, v);
882 54525825 : }
883 :
884 : template <>
885 : void
886 275872 : MooseVariableData<RealEigenVector>::insertNodalValue(NumericVector<Number> & residual,
887 : const RealEigenVector & v)
888 : {
889 983456 : for (const auto j : make_range(_count))
890 707584 : residual.set(_nodal_dof_index + j, v(j));
891 275872 : }
892 :
893 : template <typename OutputType>
894 : typename MooseVariableData<OutputType>::OutputData
895 104554 : 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 104554 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
907 :
908 104554 : switch (state)
909 : {
910 1220 : case Moose::Current:
911 1220 : return (*_sys.currentSolution())(dof);
912 :
913 103334 : case Moose::Old:
914 103334 : 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 70888 : MooseVariableData<OutputType>::getElementalValue(const Elem * elem,
967 : Moose::SolutionState state,
968 : unsigned int idx) const
969 : {
970 70888 : std::vector<dof_id_type> dof_indices;
971 70888 : _dof_map.dof_indices(elem, dof_indices, _var_num);
972 :
973 70888 : switch (state)
974 : {
975 69928 : case Moose::Current:
976 69928 : return (*_sys.currentSolution())(dof_indices[idx]);
977 :
978 480 : case Moose::Old:
979 480 : return _sys.solutionOld()(dof_indices[idx]);
980 :
981 480 : case Moose::Older:
982 480 : return _sys.solutionOlder()(dof_indices[idx]);
983 :
984 0 : default:
985 0 : mooseError("PreviousNL not currently supported for getElementalValue");
986 : }
987 70888 : }
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 18860 : MooseVariableData<OutputType>::getDofIndices(const Elem * elem,
1028 : std::vector<dof_id_type> & dof_indices) const
1029 : {
1030 18860 : _dof_map.dof_indices(elem, dof_indices, _var_num);
1031 18860 : }
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 2720 : MooseVariableData<RealEigenVector>::addSolution(NumericVector<Number> & sol,
1044 : const DenseVector<Number> & v) const
1045 : {
1046 2720 : unsigned int p = 0;
1047 2720 : const auto num_dofs = _dof_indices.size();
1048 8160 : for (const auto j : make_range(_count))
1049 : {
1050 5440 : unsigned int inc = (isNodal() ? j : j * num_dofs);
1051 27200 : for (const auto i : make_range(num_dofs))
1052 21760 : sol.add(_dof_indices[i] + inc, v(p++));
1053 : }
1054 2720 : }
1055 :
1056 : template <typename OutputType>
1057 : const typename MooseVariableData<OutputType>::DoFValue &
1058 305 : MooseVariableData<OutputType>::dofValuesDot() const
1059 : {
1060 305 : if (_sys.solutionUDot())
1061 : {
1062 305 : _need_dof_values_dot = true;
1063 305 : 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 122 : MooseVariableData<OutputType>::dofValuesDuDotDu() const
1118 : {
1119 122 : _need_dof_du_dot_du = true;
1120 122 : 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 1774 : MooseVariableData<OutputType>::computeIncrementAtQps(const NumericVector<Number> & increment_vec)
1134 : {
1135 1774 : unsigned int nqp = _qrule->n_points();
1136 :
1137 1774 : _increment.resize(nqp);
1138 : // Compute the increment at each quadrature point
1139 1774 : unsigned int num_dofs = _dof_indices.size();
1140 15108 : for (const auto qp : make_range(nqp))
1141 : {
1142 13334 : _increment[qp] = 0.;
1143 129924 : for (const auto i : make_range(num_dofs))
1144 116590 : _increment[qp] += (*_phi)[i][qp] * increment_vec(_dof_indices[i]);
1145 : }
1146 1774 : }
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 7440 : MooseVariableData<OutputType>::computeIncrementAtNode(const NumericVector<Number> & increment_vec)
1185 : {
1186 7440 : if (!isNodal())
1187 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1188 :
1189 7440 : _increment.resize(1);
1190 :
1191 : // Compute the increment for the current DOF
1192 7440 : _increment[0] = increment_vec(_dof_indices[0]);
1193 7440 : }
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 270795437 : MooseVariableData<OutputType>::computeNodalValues()
1314 : {
1315 270795437 : if (_has_dof_indices)
1316 : {
1317 270795227 : fetchDoFValues();
1318 270795227 : assignNodalValue();
1319 :
1320 270795227 : if (_need_ad)
1321 3114893 : fetchADNodalValues();
1322 : }
1323 : else
1324 210 : zeroSizeDofValues();
1325 270795437 : }
1326 :
1327 : template <typename OutputType>
1328 : void
1329 3114893 : MooseVariableData<OutputType>::fetchADNodalValues()
1330 : {
1331 3114893 : auto n = _dof_indices.size();
1332 : libmesh_assert(n);
1333 3114893 : _ad_dof_values.resize(n);
1334 :
1335 3114893 : if (_need_ad_u_dot)
1336 566899 : _ad_dofs_dot.resize(n);
1337 3114893 : if (_need_ad_u_dotdot)
1338 1216 : _ad_dofs_dotdot.resize(n);
1339 :
1340 3114893 : const bool do_derivatives =
1341 3114893 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
1342 :
1343 6396639 : for (decltype(n) i = 0; i < n; ++i)
1344 : {
1345 3281746 : _ad_dof_values[i] = _vector_tags_dof_u[_solution_tag][i];
1346 3281746 : if (do_derivatives)
1347 1341842 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
1348 3281746 : assignADNodalValue(_ad_dof_values[i], i);
1349 :
1350 3281746 : if (_need_ad_u_dot)
1351 : {
1352 574359 : if (_time_integrator && _time_integrator->dt())
1353 : {
1354 574359 : _ad_dofs_dot[i] = _ad_dof_values[i];
1355 1147502 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
1356 574359 : _dof_indices[i],
1357 574359 : _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 3114893 : }
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 2979111 : MooseVariableData<Real>::assignADNodalValue(const ADReal & value, const unsigned int &)
1380 : {
1381 2979111 : _ad_nodal_value = value;
1382 2979111 : }
1383 :
1384 : template <>
1385 : void
1386 302635 : MooseVariableData<RealVectorValue>::assignADNodalValue(const ADReal & value,
1387 : const unsigned int & component)
1388 : {
1389 302635 : _ad_nodal_value(component) = value;
1390 302635 : }
1391 :
1392 : template <typename OutputType>
1393 : void
1394 1998668 : MooseVariableData<OutputType>::prepareIC()
1395 : {
1396 1998668 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1397 1998668 : _vector_tags_dof_u[_solution_tag].resize(_dof_indices.size());
1398 :
1399 1998668 : unsigned int nqp = _qrule->n_points();
1400 1998668 : _vector_tag_u[_solution_tag].resize(nqp);
1401 1998668 : }
1402 :
1403 : template <typename OutputType>
1404 : void
1405 481567578 : MooseVariableData<OutputType>::prepare()
1406 : {
1407 481567578 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1408 481567578 : _has_dof_values = false;
1409 :
1410 : // FIXME: remove this when the Richard's module is migrated to use the new NodalCoupleable
1411 : // interface.
1412 481567578 : _has_dof_indices = _dof_indices.size();
1413 481567578 : }
1414 :
1415 : template <typename OutputType>
1416 : void
1417 307731586 : MooseVariableData<OutputType>::reinitNode()
1418 : {
1419 307731586 : if (std::size_t n_dofs = _node->n_dofs(_sys.number(), _var_num))
1420 : {
1421 270789042 : _dof_indices.resize(n_dofs);
1422 542598685 : for (std::size_t i = 0; i < n_dofs; ++i)
1423 271809643 : _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 270789042 : _nodal_dof_index = _dof_indices[0];
1427 270789042 : _has_dof_indices = true;
1428 : }
1429 : else
1430 : {
1431 36942544 : _dof_indices.clear(); // Clear these so Assembly doesn't think there's dofs here
1432 36942544 : _has_dof_indices = false;
1433 : }
1434 307731586 : }
1435 :
1436 : template <typename OutputType>
1437 : void
1438 148512161 : MooseVariableData<OutputType>::reinitAux()
1439 : {
1440 : /* FIXME: this method is only for elemental auxiliary variables, so
1441 : * we may want to rename it */
1442 148512161 : if (_elem)
1443 : {
1444 147542702 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1445 147542702 : 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 146687106 : _nodal_dof_index = _elem->dof_number(_sys.number(), _var_num, 0);
1449 :
1450 146687106 : fetchDoFValues();
1451 :
1452 146687106 : const auto num_dofs = _dof_indices.size();
1453 1207449423 : for (auto & dof_u : _vector_tags_dof_u)
1454 1060762317 : dof_u.resize(num_dofs);
1455 :
1456 440072014 : for (auto & dof_u : _matrix_tags_dof_u)
1457 293384908 : dof_u.resize(num_dofs);
1458 :
1459 146687106 : _has_dof_indices = true;
1460 : }
1461 : else
1462 855596 : _has_dof_indices = false;
1463 : }
1464 : else
1465 969459 : _has_dof_indices = false;
1466 148512161 : }
1467 :
1468 : template <typename OutputType>
1469 : void
1470 6990 : MooseVariableData<OutputType>::reinitNodes(const std::vector<dof_id_type> & nodes)
1471 : {
1472 6990 : _dof_indices.clear();
1473 37684 : for (const auto & node_id : nodes)
1474 : {
1475 30694 : auto && nd = _subproblem.mesh().getMesh().query_node_ptr(node_id);
1476 30694 : if (nd && (_subproblem.mesh().isSemiLocal(const_cast<Node *>(nd))))
1477 : {
1478 30628 : if (nd->n_dofs(_sys.number(), _var_num) > 0)
1479 : {
1480 29617 : dof_id_type dof = nd->dof_number(_sys.number(), _var_num, 0);
1481 29617 : _dof_indices.push_back(dof);
1482 : }
1483 : }
1484 : }
1485 :
1486 6990 : if (!_dof_indices.empty())
1487 6780 : _has_dof_indices = true;
1488 : else
1489 210 : _has_dof_indices = false;
1490 6990 : }
1491 :
1492 : template class MooseVariableData<Real>;
1493 : template class MooseVariableData<RealVectorValue>;
1494 : template class MooseVariableData<RealEigenVector>;
|
