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 522213 : MooseVariableData<OutputType>::MooseVariableData(const MooseVariableFE<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 522213 : _var(var),
39 1044426 : _fe_type(var.feType()),
40 522213 : _var_num(var.number()),
41 522213 : _assembly(_subproblem.assembly(_tid, var.kind() == Moose::VAR_SOLVER ? sys.number() : 0)),
42 522213 : _element_type(element_type),
43 522213 : _ad_zero(0),
44 522213 : _has_dof_indices(false),
45 522213 : _qrule(qrule_in),
46 522213 : _qrule_face(qrule_face_in),
47 522213 : _use_dual(var.useDual()),
48 522213 : _second_phi_assembly_method(nullptr),
49 522213 : _second_phi_face_assembly_method(nullptr),
50 522213 : _curl_phi_assembly_method(nullptr),
51 522213 : _curl_phi_face_assembly_method(nullptr),
52 522213 : _div_phi_assembly_method(nullptr),
53 522213 : _div_phi_face_assembly_method(nullptr),
54 522213 : _ad_grad_phi_assembly_method(nullptr),
55 522213 : _ad_grad_phi_face_assembly_method(nullptr),
56 522213 : _time_integrator(nullptr),
57 522213 : _node(node),
58 522213 : _elem(elem),
59 522213 : _displaced(dynamic_cast<const DisplacedSystem *>(&_sys) ? true : false),
60 2611065 : _current_side(_assembly.side())
61 : {
62 522213 : _continuity = FEInterface::get_continuity(_fe_type);
63 :
64 522213 : _is_nodal = (_continuity == C_ZERO || _continuity == C_ONE);
65 :
66 522213 : _time_integrator = _sys.queryTimeIntegrator(_var_num);
67 :
68 : // Initialize AD zero with zero derivatives
69 522213 : const auto old_do = ADReal::do_derivatives;
70 522213 : ADReal::do_derivatives = true;
71 522213 : _ad_zero = 0.;
72 522213 : ADReal::do_derivatives = old_do;
73 :
74 522213 : switch (_element_type)
75 : {
76 174071 : case Moose::ElementType::Element:
77 : {
78 174071 : _phi_assembly_method = &Assembly::fePhi<OutputShape>;
79 174071 : _phi_face_assembly_method = &Assembly::fePhiFace<OutputShape>;
80 174071 : _grad_phi_assembly_method = &Assembly::feGradPhi<OutputShape>;
81 174071 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFace<OutputShape>;
82 174071 : _second_phi_assembly_method = &Assembly::feSecondPhi<OutputShape>;
83 174071 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFace<OutputShape>;
84 174071 : _curl_phi_assembly_method = &Assembly::feCurlPhi<OutputShape>;
85 174071 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFace<OutputShape>;
86 174071 : _div_phi_assembly_method = &Assembly::feDivPhi<OutputShape>;
87 174071 : _div_phi_face_assembly_method = &Assembly::feDivPhiFace<OutputShape>;
88 174071 : _ad_grad_phi_assembly_method = &Assembly::feADGradPhi<OutputShape>;
89 174071 : _ad_grad_phi_face_assembly_method = &Assembly::feADGradPhiFace<OutputShape>;
90 :
91 174071 : _ad_grad_phi = &_ad_grad_phi_assembly_method(_assembly, _fe_type);
92 174071 : _ad_grad_phi_face = &_ad_grad_phi_face_assembly_method(_assembly, _fe_type);
93 174071 : break;
94 : }
95 174071 : case Moose::ElementType::Neighbor:
96 : {
97 174071 : _phi_assembly_method = &Assembly::fePhiNeighbor<OutputShape>;
98 174071 : _phi_face_assembly_method = &Assembly::fePhiFaceNeighbor<OutputShape>;
99 174071 : _grad_phi_assembly_method = &Assembly::feGradPhiNeighbor<OutputShape>;
100 174071 : _grad_phi_face_assembly_method = &Assembly::feGradPhiFaceNeighbor<OutputShape>;
101 174071 : _second_phi_assembly_method = &Assembly::feSecondPhiNeighbor<OutputShape>;
102 174071 : _second_phi_face_assembly_method = &Assembly::feSecondPhiFaceNeighbor<OutputShape>;
103 174071 : _curl_phi_assembly_method = &Assembly::feCurlPhiNeighbor<OutputShape>;
104 174071 : _curl_phi_face_assembly_method = &Assembly::feCurlPhiFaceNeighbor<OutputShape>;
105 174071 : _div_phi_assembly_method = &Assembly::feDivPhiNeighbor<OutputShape>;
106 174071 : _div_phi_face_assembly_method = &Assembly::feDivPhiFaceNeighbor<OutputShape>;
107 :
108 174071 : _ad_grad_phi = nullptr;
109 174071 : _ad_grad_phi_face = nullptr;
110 174071 : break;
111 : }
112 174071 : case Moose::ElementType::Lower:
113 : {
114 174071 : if (_use_dual)
115 : {
116 135 : _phi_assembly_method = &Assembly::feDualPhiLower<OutputType>;
117 135 : _phi_face_assembly_method = &Assembly::feDualPhiLower<OutputType>; // Place holder
118 135 : _grad_phi_assembly_method = &Assembly::feGradDualPhiLower<OutputType>;
119 135 : _grad_phi_face_assembly_method = &Assembly::feGradDualPhiLower<OutputType>; // Place holder
120 : }
121 : else
122 : {
123 173936 : _phi_assembly_method = &Assembly::fePhiLower<OutputType>;
124 173936 : _phi_face_assembly_method = &Assembly::fePhiLower<OutputType>; // Place holder
125 173936 : _grad_phi_assembly_method = &Assembly::feGradPhiLower<OutputType>;
126 173936 : _grad_phi_face_assembly_method = &Assembly::feGradPhiLower<OutputType>; // Place holder
127 : }
128 :
129 174071 : _ad_grad_phi = nullptr;
130 174071 : _ad_grad_phi_face = nullptr;
131 174071 : break;
132 : }
133 : }
134 522213 : _phi = &_phi_assembly_method(_assembly, _fe_type);
135 522213 : _phi_face = &_phi_face_assembly_method(_assembly, _fe_type);
136 522213 : _grad_phi = &_grad_phi_assembly_method(_assembly, _fe_type);
137 522213 : _grad_phi_face = &_grad_phi_face_assembly_method(_assembly, _fe_type);
138 522213 : }
139 :
140 : template <typename OutputType>
141 : void
142 696979845 : MooseVariableData<OutputType>::setGeometry(Moose::GeometryType gm_type)
143 : {
144 696979845 : switch (gm_type)
145 : {
146 662714234 : case Moose::Volume:
147 : {
148 662714234 : _current_qrule = _qrule;
149 662714234 : _current_phi = _phi;
150 662714234 : _current_grad_phi = _grad_phi;
151 662714234 : _current_second_phi = _second_phi;
152 662714234 : _current_curl_phi = _curl_phi;
153 662714234 : _current_div_phi = _div_phi;
154 662714234 : _current_ad_grad_phi = _ad_grad_phi;
155 662714234 : break;
156 : }
157 34265611 : case Moose::Face:
158 : {
159 34265611 : _current_qrule = _qrule_face;
160 34265611 : _current_phi = _phi_face;
161 34265611 : _current_grad_phi = _grad_phi_face;
162 34265611 : _current_second_phi = _second_phi_face;
163 34265611 : _current_curl_phi = _curl_phi_face;
164 34265611 : _current_div_phi = _div_phi_face;
165 34265611 : _current_ad_grad_phi = _ad_grad_phi_face;
166 34265611 : break;
167 : }
168 : }
169 696979845 : }
170 :
171 : template <typename OutputType>
172 : const typename MooseVariableData<OutputType>::FieldVariableValue &
173 15810 : MooseVariableData<OutputType>::uDot() const
174 : {
175 15810 : if (_sys.solutionUDot())
176 : {
177 15810 : _need_u_dot = true;
178 15810 : return _u_dot;
179 : }
180 : else
181 0 : mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set "
182 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
183 : }
184 :
185 : template <typename OutputType>
186 : const typename MooseVariableData<OutputType>::FieldVariableValue &
187 202 : MooseVariableData<OutputType>::uDotDot() const
188 : {
189 202 : if (_sys.solutionUDotDot())
190 : {
191 202 : _need_u_dotdot = true;
192 202 : return _u_dotdot;
193 : }
194 : else
195 0 : mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. "
196 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
197 : "`u_dotdot`.");
198 : }
199 :
200 : template <typename OutputType>
201 : const typename MooseVariableData<OutputType>::FieldVariableValue &
202 0 : MooseVariableData<OutputType>::uDotOld() const
203 : {
204 0 : if (_sys.solutionUDotOld())
205 : {
206 0 : _need_u_dot_old = true;
207 0 : return _u_dot_old;
208 : }
209 : else
210 0 : mooseError("MooseVariableFE: Old time derivative of solution (`u_dot_old`) is not stored. "
211 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
212 : "`u_dot_old`.");
213 : }
214 :
215 : template <typename OutputType>
216 : const typename MooseVariableData<OutputType>::FieldVariableValue &
217 0 : MooseVariableData<OutputType>::uDotDotOld() const
218 : {
219 0 : if (_sys.solutionUDotDotOld())
220 : {
221 0 : _need_u_dotdot_old = true;
222 0 : return _u_dotdot_old;
223 : }
224 : else
225 0 : mooseError("MooseVariableFE: Old second time derivative of solution (`u_dotdot_old`) is not "
226 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
227 : "requesting `u_dotdot_old`");
228 : }
229 :
230 : template <typename OutputType>
231 : const typename MooseVariableData<OutputType>::FieldVariableGradient &
232 36 : MooseVariableData<OutputType>::gradSlnDot() const
233 : {
234 36 : if (_sys.solutionUDot())
235 : {
236 36 : _need_grad_dot = true;
237 36 : return _grad_u_dot;
238 : }
239 : else
240 0 : mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set "
241 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
242 : }
243 :
244 : template <typename OutputType>
245 : const typename MooseVariableData<OutputType>::FieldVariableGradient &
246 0 : MooseVariableData<OutputType>::gradSlnDotDot() const
247 : {
248 0 : if (_sys.solutionUDotDot())
249 : {
250 0 : _need_grad_dotdot = true;
251 0 : return _grad_u_dotdot;
252 : }
253 : else
254 0 : mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. "
255 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
256 : "`u_dotdot`.");
257 : }
258 :
259 : template <typename OutputType>
260 : const typename MooseVariableData<OutputType>::FieldVariableSecond &
261 187 : MooseVariableData<OutputType>::secondSln(Moose::SolutionState state) const
262 : {
263 187 : secondPhi();
264 187 : secondPhiFace();
265 187 : switch (state)
266 : {
267 187 : case Moose::Current:
268 : {
269 187 : _need_second = true;
270 187 : return _second_u;
271 : }
272 :
273 0 : case Moose::Old:
274 : {
275 0 : _need_second_old = true;
276 0 : return _second_u_old;
277 : }
278 :
279 0 : case Moose::Older:
280 : {
281 0 : _need_second_older = true;
282 0 : return _second_u_older;
283 : }
284 :
285 0 : case Moose::PreviousNL:
286 : {
287 0 : _need_second_previous_nl = true;
288 0 : return _second_u_previous_nl;
289 : }
290 :
291 0 : default:
292 : // We should never get here but gcc requires that we have a default. See
293 : // htpps://stackoverflow.com/questions/18680378/after-defining-case-for-all-enum-values-compiler-still-says-control-reaches-e
294 0 : mooseError("Unknown SolutionState!");
295 : }
296 : }
297 :
298 : template <typename OutputType>
299 : const typename MooseVariableData<OutputType>::FieldVariableCurl &
300 158 : MooseVariableData<OutputType>::curlSln(Moose::SolutionState state) const
301 : {
302 158 : curlPhi();
303 158 : curlPhiFace();
304 158 : switch (state)
305 : {
306 158 : case Moose::Current:
307 : {
308 158 : _need_curl = true;
309 158 : return _curl_u;
310 : }
311 :
312 0 : case Moose::Old:
313 : {
314 0 : _need_curl_old = true;
315 0 : return _curl_u_old;
316 : }
317 :
318 0 : case Moose::Older:
319 : {
320 0 : _need_curl_older = true;
321 0 : return _curl_u_older;
322 : }
323 :
324 0 : default:
325 0 : mooseError("We don't currently support curl from the previous non-linear iteration");
326 : }
327 : }
328 :
329 : template <typename OutputType>
330 : const typename MooseVariableData<OutputType>::FieldVariableDivergence &
331 546 : MooseVariableData<OutputType>::divSln(Moose::SolutionState state) const
332 : {
333 546 : divPhi();
334 546 : divPhiFace();
335 546 : switch (state)
336 : {
337 546 : case Moose::Current:
338 : {
339 546 : _need_div = true;
340 546 : return _div_u;
341 : }
342 :
343 0 : case Moose::Old:
344 : {
345 0 : _need_div_old = true;
346 0 : return _div_u_old;
347 : }
348 :
349 0 : case Moose::Older:
350 : {
351 0 : _need_div_older = true;
352 0 : return _div_u_older;
353 : }
354 :
355 0 : default:
356 0 : mooseError("We don't currently support divergence from the previous non-linear iteration");
357 : }
358 : }
359 :
360 : template <typename OutputType>
361 : const typename MooseVariableData<OutputType>::FieldVariablePhiSecond &
362 24522 : MooseVariableData<OutputType>::secondPhi() const
363 : {
364 24522 : _second_phi = &_second_phi_assembly_method(_assembly, _fe_type);
365 24522 : return *_second_phi;
366 : }
367 :
368 : template <typename OutputType>
369 : const typename MooseVariableData<OutputType>::FieldVariablePhiSecond &
370 4448 : MooseVariableData<OutputType>::secondPhiFace() const
371 : {
372 4448 : _second_phi_face = &_second_phi_face_assembly_method(_assembly, _fe_type);
373 4448 : return *_second_phi_face;
374 : }
375 :
376 : template <typename OutputType>
377 : const typename MooseVariableData<OutputType>::FieldVariablePhiCurl &
378 57237 : MooseVariableData<OutputType>::curlPhi() const
379 : {
380 57237 : _curl_phi = &_curl_phi_assembly_method(_assembly, _fe_type);
381 57237 : return *_curl_phi;
382 : }
383 :
384 : template <typename OutputType>
385 : const typename MooseVariableData<OutputType>::FieldVariablePhiCurl &
386 207 : MooseVariableData<OutputType>::curlPhiFace() const
387 : {
388 207 : _curl_phi_face = &_curl_phi_face_assembly_method(_assembly, _fe_type);
389 207 : return *_curl_phi_face;
390 : }
391 :
392 : template <typename OutputType>
393 : const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
394 521416 : MooseVariableData<OutputType>::divPhi() const
395 : {
396 521416 : _div_phi = &_div_phi_assembly_method(_assembly, _fe_type);
397 521416 : return *_div_phi;
398 : }
399 :
400 : template <typename OutputType>
401 : const typename MooseVariableData<OutputType>::FieldVariablePhiDivergence &
402 546 : MooseVariableData<OutputType>::divPhiFace() const
403 : {
404 546 : _div_phi_face = &_div_phi_face_assembly_method(_assembly, _fe_type);
405 546 : return *_div_phi_face;
406 : }
407 :
408 : template <typename OutputType>
409 : template <bool constant_monomial,
410 : typename DestinationType,
411 : typename ShapeType,
412 : typename DofValuesType>
413 : void
414 1378930204 : MooseVariableData<OutputType>::fill(DestinationType & dest,
415 : const ShapeType & phi,
416 : const DofValuesType & dof_values,
417 : const unsigned int nqp,
418 : const std::size_t num_shapes)
419 : {
420 : if constexpr (constant_monomial)
421 160958727 : libmesh_ignore(num_shapes);
422 :
423 : // Deduce OutputType
424 1378930204 : constexpr bool is_real = std::is_same_v<OutputType, Real>;
425 1378930204 : constexpr bool is_real_vector = std::is_same_v<OutputType, RealVectorValue>;
426 1378930204 : constexpr bool is_eigen = std::is_same_v<OutputType, RealEigenVector>;
427 : static_assert(is_real || is_real_vector || is_eigen, "Unsupported type");
428 :
429 : // this is only used in the RealEigenVector case to get this->_count
430 : if constexpr (!is_eigen)
431 1370145650 : libmesh_ignore(this);
432 :
433 : // Deduce type of value within dest MooseArray
434 : using dest_array_type = typename std::remove_reference_t<decltype(dest)>::value_type;
435 1378930204 : constexpr bool is_value =
436 : std::is_same_v<dest_array_type, OutputType> ||
437 : std::is_same_v<dest_array_type, typename Moose::ADType<OutputType>::type>;
438 1378930204 : constexpr bool is_gradient =
439 : std::is_same_v<dest_array_type, OutputGradient> ||
440 : std::is_same_v<dest_array_type, typename Moose::ADType<OutputGradient>::type>;
441 1378930204 : constexpr bool is_second =
442 : std::is_same_v<dest_array_type, OutputSecond> ||
443 : std::is_same_v<dest_array_type, typename Moose::ADType<OutputSecond>::type>;
444 1378930204 : constexpr bool is_divergence =
445 : std::is_same_v<dest_array_type, OutputDivergence> ||
446 : std::is_same_v<dest_array_type, typename Moose::ADType<OutputDivergence>::type>;
447 : static_assert(is_value || is_gradient || is_second || is_divergence,
448 : "Unsupported destination array type");
449 :
450 : // Sets a value to zero at a quadrature point
451 12666852444 : const auto set_zero = [this, &dest](const auto qp)
452 : {
453 : if constexpr (!is_eigen)
454 5603533126 : libmesh_ignore(this);
455 :
456 : if constexpr (is_real || is_real_vector)
457 5603533126 : dest[qp] = 0;
458 : else if constexpr (is_eigen)
459 : {
460 : if constexpr (is_value)
461 28244734 : dest[qp].setZero(this->_count);
462 : else if constexpr (is_gradient)
463 12183260 : dest[qp].setZero(this->_count, LIBMESH_DIM);
464 : else if constexpr (is_second)
465 0 : dest[qp].setZero(this->_count, LIBMESH_DIM * LIBMESH_DIM);
466 : else
467 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
468 : }
469 : else
470 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
471 : };
472 :
473 : // Accumulates a value
474 52461597413 : const auto accumulate = [&dest, &phi, &dof_values](const auto i, const auto qp)
475 : {
476 : if constexpr (is_real || is_real_vector || (is_eigen && is_value))
477 : {
478 : if constexpr (is_value || is_divergence)
479 16021101980 : dest[qp] += phi[i][qp] * dof_values[i];
480 : else if constexpr (is_gradient || is_second)
481 9486784721 : dest[qp].add_scaled(phi[i][qp], dof_values[i]);
482 : else
483 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
484 : }
485 : else if constexpr (is_eigen)
486 : {
487 : if constexpr (is_gradient)
488 : {
489 267575228 : for (const auto d : make_range(Moose::dim))
490 200681421 : dest[qp].col(d) += phi[i][qp](d) * dof_values[i];
491 : }
492 : else if constexpr (is_second)
493 : {
494 0 : for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1)
495 0 : for (const auto d2 : make_range(Moose::dim))
496 0 : dest[qp].col(d++) += phi[i][qp](d1, d2) * dof_values[i];
497 : }
498 : else
499 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
500 : }
501 : else
502 : static_assert(Moose::always_false<OutputType, dest_array_type>, "Unsupported type");
503 : };
504 :
505 1378930204 : dest.resize(nqp);
506 :
507 : // Monomial case, accumulate dest[0] and set dest[>0] to dest[0]
508 : if constexpr (constant_monomial)
509 : {
510 : mooseAssert(num_shapes == 1, "Should have only one shape function for a constant monomial");
511 160958727 : set_zero(0);
512 160958727 : accumulate(0, 0);
513 728898509 : for (unsigned int qp = 1; qp < nqp; ++qp)
514 567939782 : dest[qp] = dest[0];
515 : }
516 : // Non constant monomial case
517 : else
518 : {
519 6700973870 : for (const auto qp : make_range(nqp))
520 5483002393 : set_zero(qp);
521 6303654117 : for (const auto i : make_range(num_shapes))
522 30499504421 : for (const auto qp : make_range(nqp))
523 25413821781 : accumulate(i, qp);
524 : }
525 1378930204 : }
526 :
527 : template <typename OutputType>
528 : template <bool constant_monomial>
529 : void
530 696338852 : MooseVariableData<OutputType>::computeValuesInternal()
531 : {
532 696338852 : const auto num_dofs = _dof_indices.size();
533 696338852 : const auto num_shapes = num_dofs / _count;
534 :
535 696338852 : if (num_dofs > 0)
536 603420182 : fetchDofValues();
537 :
538 696338852 : const bool is_transient = _subproblem.isTransient();
539 696338852 : const auto nqp = _current_qrule->n_points();
540 : const auto & active_coupleable_matrix_tags =
541 696338852 : _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid);
542 :
543 : // Map grad_phi using Eigen so that we can perform array operations easier
544 : if constexpr (std::is_same_v<OutputType, RealEigenVector>)
545 : {
546 4226956 : if (_qrule == _current_qrule)
547 : {
548 4019062 : _mapped_grad_phi.resize(num_shapes);
549 18881751 : for (const auto i : make_range(num_shapes))
550 : {
551 14862689 : _mapped_grad_phi[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
552 97256488 : for (const auto qp : make_range(nqp))
553 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
554 164787598 : new (&_mapped_grad_phi[i][qp])
555 82393799 : Eigen::Map<RealDIMValue>(const_cast<Real *>(&(*_current_grad_phi)[i][qp](0)));
556 : }
557 : }
558 : else
559 : {
560 207894 : _mapped_grad_phi_face.resize(num_shapes);
561 820342 : for (const auto i : make_range(num_shapes))
562 : {
563 612448 : _mapped_grad_phi_face[i].resize(nqp, Eigen::Map<RealDIMValue>(nullptr));
564 2121984 : for (const auto qp : make_range(nqp))
565 : // Note: this does NOT do any allocation. It is "reconstructing" the object in place
566 3019072 : new (&_mapped_grad_phi_face[i][qp])
567 1509536 : Eigen::Map<RealDIMValue>(const_cast<Real *>(&(*_current_grad_phi)[i][qp](0)));
568 : }
569 : }
570 : }
571 :
572 : mooseAssert(
573 : !(_need_second || _need_second_old || _need_second_older || _need_second_previous_nl) ||
574 : _current_second_phi,
575 : "We're requiring a second calculation but have not set a second shape function!");
576 : mooseAssert(!(_need_curl || _need_curl_old) || _current_curl_phi,
577 : "We're requiring a curl calculation but have not set a curl shape function!");
578 : mooseAssert(!(_need_div || _need_div_old) || _current_div_phi,
579 : "We're requiring a divergence calculation but have not set a div shape function!");
580 :
581 : // Curl
582 696338852 : if (_need_curl)
583 71006 : fill<constant_monomial>(
584 71006 : _curl_u, *_current_curl_phi, _vector_tags_dof_u[_solution_tag], nqp, num_shapes);
585 696338852 : if (is_transient && _need_curl_old)
586 0 : fill<constant_monomial>(
587 0 : _curl_u_old, *_current_curl_phi, _vector_tags_dof_u[_old_solution_tag], nqp, num_shapes);
588 :
589 : // Div
590 696338852 : if (_need_div)
591 426624 : fill<constant_monomial>(
592 426624 : _div_u, *_current_div_phi, _vector_tags_dof_u[_solution_tag], nqp, num_shapes);
593 696338852 : if (is_transient && _need_div_old)
594 0 : fill<constant_monomial>(
595 0 : _div_u_old, *_current_div_phi, _vector_tags_dof_u[_old_solution_tag], nqp, num_shapes);
596 :
597 : // Second
598 696338852 : if (_need_second)
599 80167 : fill<constant_monomial>(
600 80167 : _second_u, *_current_second_phi, _vector_tags_dof_u[_solution_tag], nqp, num_shapes);
601 696338852 : if (_need_second_previous_nl)
602 0 : fill<constant_monomial>(_second_u_previous_nl,
603 0 : *_current_second_phi,
604 0 : _vector_tags_dof_u[_previous_nl_solution_tag],
605 : nqp,
606 : num_shapes);
607 :
608 : // Vector tags
609 1399232438 : for (auto tag : _required_vector_tags)
610 : {
611 702893586 : if (_need_vector_tag_u[tag] && _sys.hasVector(tag))
612 : {
613 : mooseAssert(_sys.getVector(tag).closed(), "Vector should be closed");
614 700361386 : fill<constant_monomial>(
615 700361386 : _vector_tag_u[tag], *_current_phi, _vector_tags_dof_u[tag], nqp, num_shapes);
616 : }
617 702893586 : if (_need_vector_tag_grad[tag] && _sys.hasVector(tag))
618 : {
619 : mooseAssert(_sys.getVector(tag).closed(), "Vector should be closed");
620 398456747 : fill<constant_monomial>(
621 398456747 : _vector_tag_grad[tag], *_current_grad_phi, _vector_tags_dof_u[tag], nqp, num_shapes);
622 : }
623 : }
624 :
625 : // Matrix tags
626 696339236 : for (auto tag : active_coupleable_matrix_tags)
627 384 : if (_need_matrix_tag_u[tag])
628 128 : fill<constant_monomial>(
629 128 : _matrix_tag_u[tag], *_current_phi, _matrix_tags_dof_u[tag], nqp, num_shapes);
630 :
631 : // Derivatives and old values
632 696338852 : if (is_transient)
633 : {
634 529073097 : if (_need_second_old)
635 0 : fill<constant_monomial>(_second_u_old,
636 0 : *_current_second_phi,
637 0 : _vector_tags_dof_u[_old_solution_tag],
638 : nqp,
639 : num_shapes);
640 529073097 : if (_need_second_older)
641 0 : fill<constant_monomial>(_second_u_older,
642 0 : *_current_second_phi,
643 0 : _vector_tags_dof_u[_older_solution_tag],
644 : nqp,
645 : num_shapes);
646 529073097 : if (_need_u_dot)
647 238972611 : fill<constant_monomial>(_u_dot, *_current_phi, _dof_values_dot, nqp, num_shapes);
648 529073097 : if (_need_u_dotdot)
649 6378 : fill<constant_monomial>(_u_dotdot, *_current_phi, _dof_values_dotdot, nqp, num_shapes);
650 529073097 : if (_need_u_dot_old)
651 0 : fill<constant_monomial>(_u_dot_old, *_current_phi, _dof_values_dot_old, nqp, num_shapes);
652 529073097 : if (_need_u_dotdot_old)
653 0 : fill<constant_monomial>(
654 0 : _u_dotdot_old, *_current_phi, _dof_values_dotdot_old, nqp, num_shapes);
655 :
656 529073097 : if (_need_du_dot_du)
657 : {
658 239094856 : _du_dot_du.resize(nqp);
659 1280824765 : for (const auto i : make_range(num_shapes))
660 6028656765 : for (const auto qp : make_range(nqp))
661 4986926856 : _du_dot_du[qp] = _dof_du_dot_du[i];
662 : }
663 529073097 : if (_need_du_dotdot_du)
664 : {
665 5198 : _du_dotdot_du.resize(nqp);
666 43996 : for (const auto i : make_range(num_shapes))
667 192796 : for (const auto qp : make_range(nqp))
668 153998 : _du_dotdot_du[qp] = _dof_du_dotdot_du[i];
669 : }
670 :
671 529073097 : if (_need_grad_dot)
672 125502 : fill<constant_monomial>(_grad_u_dot, *_current_grad_phi, _dof_values_dot, nqp, num_shapes);
673 529073097 : if (_need_grad_dotdot)
674 0 : fill<constant_monomial>(
675 0 : _grad_u_dotdot, *_current_grad_phi, _dof_values_dotdot, nqp, num_shapes);
676 : }
677 :
678 696338852 : if (_need_ad)
679 20471932 : computeAD<constant_monomial>(num_dofs, nqp);
680 696338852 : }
681 :
682 : template <typename OutputType>
683 : void
684 545089169 : MooseVariableData<OutputType>::computeValues()
685 : {
686 545089169 : computeValuesInternal</* constant_monomial = */ false>();
687 545089169 : }
688 :
689 : template <typename OutputType>
690 : void
691 152027104 : MooseVariableData<OutputType>::computeConstantMonomialValues()
692 : {
693 152027104 : if (_dof_indices.size() == 0)
694 640993 : return;
695 :
696 : // Monomial optimizations are not appropriate after p-refinement
697 151386111 : if (_elem->p_level())
698 136428 : computeValues();
699 : else
700 151249683 : computeValuesInternal</* constant_monomial = */ true>();
701 : }
702 :
703 : template <typename OutputType>
704 : void
705 23466328 : MooseVariableData<OutputType>::fetchADDofValues()
706 : {
707 23466328 : const auto num_dofs = _dof_indices.size();
708 :
709 23466328 : const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys);
710 :
711 23466328 : _ad_dof_values.resize(num_dofs);
712 110168890 : for (const auto i : make_range(num_dofs))
713 86702562 : _ad_dof_values[i] = _vector_tags_dof_u[_solution_tag][i];
714 : // NOTE! You have to do this AFTER setting the value!
715 23466328 : if (do_derivatives)
716 18130212 : for (const auto i : make_range(num_dofs))
717 13988217 : Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.);
718 :
719 23466328 : const bool is_transient = _subproblem.isTransient();
720 23466328 : if (is_transient && _need_ad_u_dot)
721 : {
722 2740880 : _ad_dofs_dot.resize(num_dofs);
723 2740880 : if (_need_ad_u_dotdot)
724 1252 : _ad_dofs_dotdot.resize(num_dofs);
725 :
726 2740880 : if (_time_integrator)
727 : {
728 2740880 : if (_time_integrator->dt())
729 : {
730 12365778 : for (const auto i : make_range(num_dofs))
731 9644250 : _ad_dofs_dot[i] = _ad_dof_values[i];
732 12365778 : for (const auto i : make_range(num_dofs))
733 19286502 : _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i],
734 9644250 : _dof_indices[i],
735 9644250 : _need_ad_u_dotdot ? _ad_dofs_dotdot[i]
736 : : _ad_real_dummy);
737 : }
738 : else
739 : // Executing something with a time derivative at initial should not put a NaN
740 95724 : for (const auto i : make_range(num_dofs))
741 : {
742 76372 : _ad_dofs_dot[i] = 0.;
743 76372 : if (_need_ad_u_dotdot)
744 64 : _ad_dofs_dotdot[i] = 0;
745 : }
746 : }
747 : // We are too early in the setup to have a time integrator, so we are not really using the
748 : // AD-derivatives. We set the AD value of the derivatives to the nonAD value
749 : else
750 0 : for (const auto i : make_range(num_dofs))
751 : {
752 0 : _ad_dofs_dot[i] = _dof_values_dot[i];
753 0 : if (_need_ad_u_dotdot)
754 0 : _ad_dofs_dotdot[i] = _dof_values_dotdot[i];
755 : }
756 : }
757 23466328 : }
758 :
759 : template <>
760 : void
761 38798 : MooseVariableData<RealEigenVector>::fetchADDofValues()
762 : {
763 38798 : const auto num_dofs = _dof_indices.size();
764 :
765 38798 : const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys);
766 38798 : const auto n_test = num_dofs / _count;
767 : mooseAssert(num_dofs == _count * n_test,
768 : "Our assertions around number of dofs, test functions, and count are incorrect");
769 :
770 38798 : _ad_dof_values.resize(n_test);
771 : // Test is outer, count is inner
772 147032 : for (const auto i : make_range(n_test))
773 : {
774 108234 : _ad_dof_values[i].resize(_count);
775 324702 : for (const auto j : make_range(_count))
776 : {
777 216468 : auto & dual_number = _ad_dof_values[i](j);
778 216468 : const auto global_dof_index = _dof_indices[j * n_test + i];
779 216468 : dual_number = (*_sys.currentSolution())(global_dof_index);
780 : // NOTE! You have to do this AFTER setting the value!
781 216468 : if (do_derivatives)
782 16172 : Moose::derivInsert(dual_number.derivatives(), global_dof_index, 1.);
783 : }
784 : }
785 38798 : }
786 :
787 : template <typename OutputType>
788 : template <bool constant_monomial>
789 : void
790 20471932 : MooseVariableData<OutputType>::computeAD(const unsigned int num_dofs, const unsigned int nqp)
791 : {
792 20471932 : fetchADDofValues();
793 20471932 : const auto n_test = num_dofs / _count;
794 :
795 : // Values
796 20471932 : if (_need_ad_u)
797 19150413 : fill<constant_monomial>(_ad_u, *_current_phi, _ad_dof_values, nqp, n_test);
798 : // Grad
799 20471932 : if (_need_ad_grad_u)
800 : {
801 : // The latter check here is for handling the fact that we have not yet implemented
802 : // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that
803 : // situation we need to default to using the non-ad grad_phi
804 16596918 : if (_displaced && _current_ad_grad_phi)
805 472706 : fill<constant_monomial>(_ad_grad_u, *_current_ad_grad_phi, _ad_dof_values, nqp, n_test);
806 : else
807 16124212 : fill<constant_monomial>(_ad_grad_u, *_current_grad_phi, _ad_dof_values, nqp, n_test);
808 : }
809 : // Second
810 : if constexpr (std::is_same_v<OutputType, Real>)
811 17279736 : if (_need_ad_second_u)
812 0 : fill<constant_monomial>(_ad_second_u, *_current_second_phi, _ad_dof_values, nqp, n_test);
813 : // Curl
814 20471932 : if (_need_ad_curl_u)
815 2452600 : fill<constant_monomial>(_ad_curl_u, *_current_curl_phi, _ad_dof_values, nqp, n_test);
816 :
817 20471932 : const bool is_transient = _subproblem.isTransient();
818 20471932 : if (is_transient)
819 : {
820 10824251 : if (_need_ad_u_dot)
821 : {
822 2229324 : if (_time_integrator)
823 2229324 : fill<constant_monomial>(_ad_u_dot, *_current_phi, _ad_dofs_dot, nqp, n_test);
824 : // We are too early in the setup to have a time integrator, so we are not really using the
825 : // AD-derivatives. We set the AD value of the derivatives to the nonAD value
826 : else
827 : {
828 0 : _ad_u_dot.resize(nqp);
829 0 : for (const auto qp : make_range(nqp))
830 0 : _ad_u_dot[qp] = _u_dot[qp];
831 : }
832 : }
833 :
834 10824251 : if (_need_ad_u_dotdot)
835 : {
836 400 : if (_time_integrator)
837 400 : fill<constant_monomial>(_ad_u_dotdot, *_current_phi, _ad_dofs_dotdot, nqp, n_test);
838 : else
839 : {
840 0 : _ad_u_dotdot.resize(nqp);
841 0 : for (const auto qp : make_range(nqp))
842 0 : _ad_u_dotdot[qp] = _u_dotdot[qp];
843 : }
844 : }
845 :
846 10824251 : if (_need_ad_grad_u_dot)
847 : {
848 0 : if (_time_integrator)
849 : {
850 : // The latter check here is for handling the fact that we have not yet implemented
851 : // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that
852 : // situation we need to default to using the non-ad grad_phi
853 0 : if (_displaced && _current_ad_grad_phi)
854 0 : fill<constant_monomial>(_ad_grad_u_dot, *_current_ad_grad_phi, _ad_dofs_dot, nqp, n_test);
855 : else
856 0 : fill<constant_monomial>(_ad_grad_u_dot, *_current_grad_phi, _ad_dofs_dot, nqp, n_test);
857 : }
858 : else
859 : {
860 0 : _ad_grad_u_dot.resize(nqp);
861 0 : for (const auto qp : make_range(nqp))
862 0 : _ad_grad_u_dot[qp] = _grad_u_dot[qp];
863 : }
864 : }
865 : }
866 20471932 : }
867 :
868 : template <typename OutputType>
869 : void
870 8421231 : MooseVariableData<OutputType>::setDofValue(const DofValue & value, unsigned int index)
871 : {
872 8421231 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
873 8421231 : dof_values[index] = value;
874 8421231 : _has_dof_values = true;
875 :
876 8421231 : auto & u = _vector_tag_u[_solution_tag];
877 8421231 : const auto nqps = u.size();
878 8421231 : const auto ndofs = dof_values.size();
879 42752917 : for (const auto qp : make_range(nqps))
880 34331686 : u[qp] *= 0.;
881 42752917 : for (const auto qp : make_range(nqps))
882 292452170 : for (const auto i : make_range(ndofs))
883 258120484 : u[qp] += (*_phi)[i][qp] * dof_values[i];
884 8421231 : }
885 :
886 : template <typename OutputType>
887 : void
888 582887 : MooseVariableData<OutputType>::setDofValues(const DenseVector<DofValue> & values)
889 : {
890 582887 : auto & dof_values = _vector_tags_dof_u[_solution_tag];
891 3453696 : for (unsigned int i = 0; i < values.size(); i++)
892 2870809 : dof_values[i] = values(i);
893 :
894 582887 : _has_dof_values = true;
895 :
896 582887 : auto & u = _vector_tag_u[_solution_tag];
897 582887 : const auto nqps = u.size();
898 582887 : const auto ndofs = dof_values.size();
899 5007718 : for (const auto qp : make_range(nqps))
900 4424831 : u[qp] *= 0.;
901 5007718 : for (const auto qp : make_range(nqps))
902 32308741 : for (const auto i : make_range(ndofs))
903 27883910 : u[qp] += (*_phi)[i][qp] * dof_values[i];
904 582887 : }
905 :
906 : template <typename OutputType>
907 : void
908 55397720 : MooseVariableData<OutputType>::insertNodalValue(NumericVector<Number> & residual,
909 : const DofValue & v)
910 : {
911 55397720 : residual.set(_nodal_dof_index, v);
912 55397720 : }
913 :
914 : template <>
915 : void
916 293942 : MooseVariableData<RealEigenVector>::insertNodalValue(NumericVector<Number> & residual,
917 : const RealEigenVector & v)
918 : {
919 1039426 : for (const auto j : make_range(_count))
920 745484 : residual.set(_nodal_dof_index + j, v(j));
921 293942 : }
922 :
923 : template <typename OutputType>
924 : typename MooseVariableData<OutputType>::DofValue
925 110099 : MooseVariableData<OutputType>::getNodalValue(const Node & node, Moose::SolutionState state) const
926 : {
927 : mooseAssert(_subproblem.mesh().isSemiLocal(const_cast<Node *>(&node)), "Node is not Semilocal");
928 :
929 : // Make sure that the node has DOFs
930 : /* Note, this is a reproduction of an assert within libMesh::Node::dof_number, this is done to
931 : * produce a better error (see misc/check_error.node_value_off_block) */
932 : mooseAssert(node.n_dofs(_sys.number(), _var_num) > 0,
933 : "Node " << node.id() << " does not contain any dofs for the "
934 : << _sys.system().variable_name(_var_num) << " variable");
935 :
936 110099 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
937 :
938 110099 : switch (state)
939 : {
940 1199 : case Moose::Current:
941 1199 : return (*_sys.currentSolution())(dof);
942 :
943 108900 : case Moose::Old:
944 108900 : return _sys.solutionOld()(dof);
945 :
946 0 : case Moose::Older:
947 0 : return _sys.solutionOlder()(dof);
948 :
949 0 : default:
950 0 : mooseError("PreviousNL not currently supported for getNodalValue");
951 : }
952 : }
953 :
954 : template <>
955 : RealEigenVector
956 0 : MooseVariableData<RealEigenVector>::getNodalValue(const Node & node,
957 : Moose::SolutionState state) const
958 : {
959 : mooseAssert(_subproblem.mesh().isSemiLocal(const_cast<Node *>(&node)), "Node is not Semilocal");
960 :
961 : // Make sure that the node has DOFs
962 : /* Note, this is a reproduction of an assert within libMesh::Node::dof_number, this is done to
963 : * produce a better error (see misc/check_error.node_value_off_block) */
964 : mooseAssert(node.n_dofs(_sys.number(), _var_num) > 0,
965 : "Node " << node.id() << " does not contain any dofs for the "
966 : << _sys.system().variable_name(_var_num) << " variable");
967 :
968 0 : dof_id_type dof = node.dof_number(_sys.number(), _var_num, 0);
969 :
970 0 : RealEigenVector v(_count);
971 0 : switch (state)
972 : {
973 0 : case Moose::Current:
974 0 : for (unsigned int i = 0; i < _count; ++i)
975 0 : v(i) = (*_sys.currentSolution())(dof++);
976 0 : break;
977 :
978 0 : case Moose::Old:
979 0 : for (unsigned int i = 0; i < _count; ++i)
980 0 : v(i) = _sys.solutionOld()(dof++);
981 0 : break;
982 :
983 0 : case Moose::Older:
984 0 : for (unsigned int i = 0; i < _count; ++i)
985 0 : v(i) = _sys.solutionOlder()(dof++);
986 0 : break;
987 :
988 0 : default:
989 0 : mooseError("PreviousNL not currently supported for getNodalValue");
990 : }
991 0 : return v;
992 0 : }
993 :
994 : template <typename OutputType>
995 : typename MooseVariableData<OutputType>::DofValue
996 70787 : MooseVariableData<OutputType>::getElementalValue(const Elem * elem,
997 : const Moose::SolutionState state,
998 : const unsigned int idx) const
999 : {
1000 70787 : static thread_local std::vector<dof_id_type> dof_indices;
1001 70787 : _dof_map.dof_indices(elem, dof_indices, _var_num);
1002 :
1003 70787 : switch (state)
1004 : {
1005 69803 : case Moose::Current:
1006 69803 : return (*_sys.currentSolution())(dof_indices[idx]);
1007 :
1008 492 : case Moose::Old:
1009 492 : return _sys.solutionOld()(dof_indices[idx]);
1010 :
1011 492 : case Moose::Older:
1012 492 : return _sys.solutionOlder()(dof_indices[idx]);
1013 :
1014 0 : default:
1015 0 : mooseError("PreviousNL not currently supported for getElementalValue");
1016 : }
1017 : }
1018 :
1019 : template <>
1020 : RealEigenVector
1021 0 : MooseVariableData<RealEigenVector>::getElementalValue(const Elem * elem,
1022 : const Moose::SolutionState state,
1023 : const unsigned int idx) const
1024 : {
1025 0 : mooseDeprecated(
1026 : "getElementalValue has a really bad API name. It is retrieving a value from the solution "
1027 : "vector for a particular dof index. Generally speaking it has absolutely no equivalence to "
1028 : "an 'elemental' value, which most people would consider to be something like an element "
1029 : "average value");
1030 :
1031 0 : static thread_local std::vector<dof_id_type> dof_indices;
1032 0 : _dof_map.array_dof_indices(elem, dof_indices, _var_num);
1033 : mooseAssert(dof_indices.size() % _count == 0,
1034 : "The number of array dof indices should divide cleanly by the variable count");
1035 0 : const auto num_shapes = dof_indices.size() / _count;
1036 :
1037 0 : RealEigenVector v(_count);
1038 :
1039 0 : switch (state)
1040 : {
1041 0 : case Moose::Current:
1042 0 : for (unsigned int i = 0; i < _count; ++i)
1043 0 : v(i) = (*_sys.currentSolution())(dof_indices[i * num_shapes + idx]);
1044 0 : break;
1045 :
1046 0 : case Moose::Old:
1047 0 : for (unsigned int i = 0; i < _count; ++i)
1048 0 : v(i) = _sys.solutionOld()(dof_indices[i * num_shapes + idx]);
1049 0 : break;
1050 :
1051 0 : case Moose::Older:
1052 0 : for (unsigned int i = 0; i < _count; ++i)
1053 0 : v(i) = _sys.solutionOlder()(dof_indices[i * num_shapes + idx]);
1054 0 : break;
1055 :
1056 0 : default:
1057 0 : mooseError("PreviousNL not currently supported for getElementalValue");
1058 : }
1059 0 : return v;
1060 0 : }
1061 :
1062 : template <typename OutputType>
1063 : void
1064 18960 : MooseVariableData<OutputType>::getDofIndices(const Elem * elem,
1065 : std::vector<dof_id_type> & dof_indices) const
1066 : {
1067 : if constexpr (std::is_same<OutputType, RealEigenVector>::value)
1068 1928 : _dof_map.array_dof_indices(elem, dof_indices, _var_num);
1069 : else
1070 17032 : _dof_map.dof_indices(elem, dof_indices, _var_num);
1071 18960 : }
1072 :
1073 : template <typename OutputType>
1074 : void
1075 2720 : MooseVariableData<OutputType>::addSolution(NumericVector<Number> & sol,
1076 : const DenseVector<Number> & v) const
1077 : {
1078 2720 : sol.add_vector(v, _dof_indices);
1079 2720 : }
1080 :
1081 : template <typename OutputType>
1082 : const typename MooseVariableData<OutputType>::DofValues &
1083 348 : MooseVariableData<OutputType>::dofValuesDot() const
1084 : {
1085 348 : if (_sys.solutionUDot())
1086 : {
1087 348 : _need_dof_values_dot = true;
1088 348 : return _dof_values_dot;
1089 : }
1090 : else
1091 0 : mooseError("MooseVariableData: Time derivative of solution (`u_dot`) is not stored. Please set "
1092 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
1093 : }
1094 :
1095 : template <typename OutputType>
1096 : const typename MooseVariableData<OutputType>::DofValues &
1097 3 : MooseVariableData<OutputType>::dofValuesDotDot() const
1098 : {
1099 3 : if (_sys.solutionUDotDot())
1100 : {
1101 3 : _need_dof_values_dotdot = true;
1102 3 : return _dof_values_dotdot;
1103 : }
1104 : else
1105 0 : mooseError("MooseVariableData: Second time derivative of solution (`u_dotdot`) is not stored. "
1106 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
1107 : "`u_dotdot`.");
1108 : }
1109 :
1110 : template <typename OutputType>
1111 : const typename MooseVariableData<OutputType>::DofValues &
1112 0 : MooseVariableData<OutputType>::dofValuesDotOld() const
1113 : {
1114 0 : if (_sys.solutionUDotOld())
1115 : {
1116 0 : _need_dof_values_dot_old = true;
1117 0 : return _dof_values_dot_old;
1118 : }
1119 : else
1120 0 : mooseError("MooseVariableData: Old time derivative of solution (`u_dot_old`) is not stored. "
1121 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
1122 : "`u_dot_old`.");
1123 : }
1124 :
1125 : template <typename OutputType>
1126 : const typename MooseVariableData<OutputType>::DofValues &
1127 0 : MooseVariableData<OutputType>::dofValuesDotDotOld() const
1128 : {
1129 0 : if (_sys.solutionUDotDotOld())
1130 : {
1131 0 : _need_dof_values_dotdot_old = true;
1132 0 : return _dof_values_dotdot_old;
1133 : }
1134 : else
1135 0 : mooseError("MooseVariableData: Old second time derivative of solution (`u_dotdot_old`) is not "
1136 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
1137 : "requesting `u_dotdot_old`.");
1138 : }
1139 :
1140 : template <typename OutputType>
1141 : const MooseArray<Number> &
1142 165 : MooseVariableData<OutputType>::dofValuesDuDotDu() const
1143 : {
1144 165 : _need_dof_du_dot_du = true;
1145 165 : return _dof_du_dot_du;
1146 : }
1147 :
1148 : template <typename OutputType>
1149 : const MooseArray<Number> &
1150 0 : MooseVariableData<OutputType>::dofValuesDuDotDotDu() const
1151 : {
1152 0 : _need_dof_du_dotdot_du = true;
1153 0 : return _dof_du_dotdot_du;
1154 : }
1155 :
1156 : template <typename OutputType>
1157 : void
1158 1563 : MooseVariableData<OutputType>::computeIncrementAtQps(const NumericVector<Number> & increment_vec)
1159 : {
1160 1563 : unsigned int nqp = _qrule->n_points();
1161 :
1162 1563 : _increment.resize(nqp);
1163 : // Compute the increment at each quadrature point
1164 1563 : unsigned int num_dofs = _dof_indices.size();
1165 13082 : for (const auto qp : make_range(nqp))
1166 : {
1167 11519 : _increment[qp] = 0.;
1168 112278 : for (const auto i : make_range(num_dofs))
1169 100759 : _increment[qp] += (*_phi)[i][qp] * increment_vec(_dof_indices[i]);
1170 : }
1171 1563 : }
1172 :
1173 : template <>
1174 : void
1175 0 : MooseVariableData<RealEigenVector>::computeIncrementAtQps(
1176 : const NumericVector<Number> & increment_vec)
1177 : {
1178 0 : unsigned int nqp = _qrule->n_points();
1179 :
1180 0 : _increment.resize(nqp);
1181 : // Compute the increment at each quadrature point
1182 0 : unsigned int num_dofs = _dof_indices.size();
1183 0 : if (isNodal())
1184 : {
1185 0 : for (const auto qp : make_range(nqp))
1186 : {
1187 0 : for (const auto i : make_range(num_dofs))
1188 0 : for (const auto j : make_range(_count))
1189 0 : _increment[qp](j) += (*_phi)[i][qp] * increment_vec(_dof_indices[i] + j);
1190 : }
1191 : }
1192 : else
1193 : {
1194 0 : for (const auto qp : make_range(nqp))
1195 : {
1196 0 : unsigned int n = 0;
1197 0 : for (const auto j : make_range(_count))
1198 0 : for (const auto i : make_range(num_dofs))
1199 : {
1200 0 : _increment[qp](j) += (*_phi)[i][qp] * increment_vec(_dof_indices[i] + n);
1201 0 : n += num_dofs;
1202 : }
1203 : }
1204 : }
1205 0 : }
1206 :
1207 : template <typename OutputType>
1208 : void
1209 7525 : MooseVariableData<OutputType>::computeIncrementAtNode(const NumericVector<Number> & increment_vec)
1210 : {
1211 7525 : if (!isNodal())
1212 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1213 :
1214 7525 : _increment.resize(1);
1215 :
1216 : // Compute the increment for the current DOF
1217 7525 : _increment[0] = increment_vec(_dof_indices[0]);
1218 7525 : }
1219 :
1220 : template <>
1221 : void
1222 0 : MooseVariableData<RealEigenVector>::computeIncrementAtNode(
1223 : const NumericVector<Number> & increment_vec)
1224 : {
1225 0 : if (!isNodal())
1226 0 : mooseError("computeIncrementAtNode can only be called for nodal variables");
1227 :
1228 0 : _increment.resize(1);
1229 :
1230 : // Compute the increment for the current DOF
1231 0 : if (isNodal())
1232 0 : for (unsigned int j = 0; j < _count; j++)
1233 0 : _increment[0](j) = increment_vec(_dof_indices[0] + j);
1234 : else
1235 : {
1236 0 : unsigned int n = 0;
1237 0 : const auto n_dof_indices = _dof_indices.size();
1238 0 : for (const auto j : make_range(_count))
1239 : {
1240 0 : _increment[0](j) = increment_vec(_dof_indices[0] + n);
1241 0 : n += n_dof_indices;
1242 : }
1243 : }
1244 0 : }
1245 :
1246 : template <typename OutputType>
1247 : const OutputType &
1248 0 : MooseVariableData<OutputType>::nodalValueDot() const
1249 : {
1250 0 : if (isNodal())
1251 : {
1252 0 : if (_sys.solutionUDot())
1253 : {
1254 0 : _need_dof_values_dot = true;
1255 0 : return _nodal_value_dot;
1256 : }
1257 : else
1258 0 : mooseError(
1259 : "MooseVariableData: Time derivative of solution (`u_dot`) is not stored. Please set "
1260 : "uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
1261 : }
1262 : else
1263 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1264 0 : var().name(),
1265 : "' is not nodal.");
1266 : }
1267 :
1268 : template <typename OutputType>
1269 : const OutputType &
1270 0 : MooseVariableData<OutputType>::nodalValueDotDot() const
1271 : {
1272 0 : if (isNodal())
1273 : {
1274 0 : if (_sys.solutionUDotDot())
1275 : {
1276 0 : _need_dof_values_dotdot = true;
1277 0 : return _nodal_value_dotdot;
1278 : }
1279 : else
1280 0 : mooseError(
1281 : "MooseVariableData: Second time derivative of solution (`u_dotdot`) is not stored. "
1282 : "Please set uDotDotRequested() to true in FEProblemBase before requesting "
1283 : "`u_dotdot`.");
1284 : }
1285 : else
1286 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1287 0 : var().name(),
1288 : "' is not nodal.");
1289 : }
1290 :
1291 : template <typename OutputType>
1292 : const OutputType &
1293 0 : MooseVariableData<OutputType>::nodalValueDotOld() const
1294 : {
1295 0 : if (isNodal())
1296 : {
1297 0 : if (_sys.solutionUDotOld())
1298 : {
1299 0 : _need_dof_values_dot_old = true;
1300 0 : return _nodal_value_dot_old;
1301 : }
1302 : else
1303 0 : mooseError("MooseVariableData: Old time derivative of solution (`u_dot_old`) is not stored. "
1304 : "Please set uDotOldRequested() to true in FEProblemBase before requesting "
1305 : "`u_dot_old`.");
1306 : }
1307 : else
1308 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1309 0 : var().name(),
1310 : "' is not nodal.");
1311 : }
1312 :
1313 : template <typename OutputType>
1314 : const OutputType &
1315 0 : MooseVariableData<OutputType>::nodalValueDotDotOld() const
1316 : {
1317 0 : if (isNodal())
1318 : {
1319 0 : if (_sys.solutionUDotDotOld())
1320 : {
1321 0 : _need_dof_values_dotdot_old = true;
1322 0 : return _nodal_value_dotdot_old;
1323 : }
1324 : else
1325 0 : mooseError(
1326 : "MooseVariableData: Old second time derivative of solution (`u_dotdot_old`) is not "
1327 : "stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
1328 : "requesting `u_dotdot_old`.");
1329 : }
1330 : else
1331 0 : mooseError("Nodal values can be requested only on nodal variables, variable '",
1332 0 : var().name(),
1333 : "' is not nodal.");
1334 : }
1335 :
1336 : template <typename OutputType>
1337 : void
1338 246348834 : MooseVariableData<OutputType>::computeNodalValues()
1339 : {
1340 246348834 : if (_has_dof_indices)
1341 : {
1342 246348630 : fetchDofValues();
1343 246348630 : assignNodalValue();
1344 :
1345 246348630 : if (_need_ad)
1346 : {
1347 3033194 : fetchADDofValues();
1348 3033194 : assignADNodalValue();
1349 : }
1350 : }
1351 : else
1352 204 : zeroSizeDofValues();
1353 246348834 : }
1354 :
1355 : template <typename OutputType>
1356 : void
1357 2903180 : MooseVariableData<OutputType>::assignADNodalValue()
1358 : {
1359 : mooseAssert(_ad_dof_values.size(), "The AD dof values container must have size greater than 0");
1360 2903180 : _ad_nodal_value = _ad_dof_values[0];
1361 2903180 : }
1362 :
1363 : template <>
1364 : void
1365 130014 : MooseVariableData<RealVectorValue>::assignADNodalValue()
1366 : {
1367 130014 : const auto num_dofs = _dof_indices.size();
1368 : mooseAssert(_ad_dof_values.size() == num_dofs,
1369 : "Our dof values container size should match the dof indices container size");
1370 387858 : for (const auto i : make_range(num_dofs))
1371 257844 : _ad_nodal_value(i) = _ad_dof_values[i];
1372 130014 : }
1373 :
1374 : template <typename OutputType>
1375 : void
1376 2207079 : MooseVariableData<OutputType>::prepareIC()
1377 : {
1378 : if constexpr (std::is_same<RealEigenVector, OutputType>::value)
1379 55374 : _dof_map.array_dof_indices(_elem, _dof_indices, _var_num);
1380 : else
1381 2151705 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1382 :
1383 : mooseAssert(_dof_indices.size() % _count == 0,
1384 : "The number of dof indices should divide cleanly by the variable count");
1385 2207079 : const auto num_shapes = _dof_indices.size() / _count;
1386 2207079 : _vector_tags_dof_u[_solution_tag].resize(num_shapes);
1387 :
1388 2207079 : unsigned int nqp = _qrule->n_points();
1389 2207079 : _vector_tag_u[_solution_tag].resize(nqp);
1390 2207079 : }
1391 :
1392 : template <typename OutputType>
1393 : void
1394 473876275 : MooseVariableData<OutputType>::prepare()
1395 : {
1396 : if constexpr (std::is_same<OutputType, RealEigenVector>::value)
1397 3506105 : _dof_map.array_dof_indices(_elem, _dof_indices, _var_num);
1398 : else
1399 470370170 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1400 :
1401 473876275 : _has_dof_values = false;
1402 :
1403 : // FIXME: remove this when the Richard's module is migrated to use the new NodalCoupleable
1404 : // interface.
1405 473876275 : _has_dof_indices = _dof_indices.size();
1406 473876275 : }
1407 :
1408 : template <typename OutputType>
1409 : void
1410 282198903 : MooseVariableData<OutputType>::reinitNode()
1411 : {
1412 : if constexpr (std::is_same<OutputType, RealEigenVector>::value)
1413 2304330 : _dof_map.array_dof_indices(_node, _dof_indices, _var_num);
1414 : else
1415 279894573 : _dof_map.dof_indices(_node, _dof_indices, _var_num);
1416 :
1417 282198903 : const auto n_dofs = _dof_indices.size();
1418 282198903 : if (n_dofs)
1419 : {
1420 : // For standard variables. _nodal_dof_index is retrieved by nodalDofIndex() which is used in
1421 : // NodalBC for example
1422 246342821 : _nodal_dof_index = _dof_indices[0];
1423 246342821 : _has_dof_indices = true;
1424 : }
1425 : else
1426 35856082 : _has_dof_indices = false;
1427 282198903 : }
1428 :
1429 : template <typename OutputType>
1430 : void
1431 152208408 : MooseVariableData<OutputType>::reinitAux()
1432 : {
1433 : /* FIXME: this method is only for elemental auxiliary variables, so
1434 : * we may want to rename it */
1435 152208408 : if (_elem)
1436 : {
1437 : if constexpr (std::is_same<RealEigenVector, OutputType>::value)
1438 79714 : _dof_map.array_dof_indices(_elem, _dof_indices, _var_num);
1439 : else
1440 151182461 : _dof_map.dof_indices(_elem, _dof_indices, _var_num);
1441 151262175 : if (_elem->n_dofs(_sys.number(), _var_num) > 0)
1442 : {
1443 : // FIXME: Setting _nodal_dof_index inside of a method that apparently is only for "elemental"
1444 : // aux variables seems absolutely absurd
1445 150389846 : _nodal_dof_index = _dof_indices[0];
1446 :
1447 150389846 : fetchDofValues();
1448 :
1449 : mooseAssert(_dof_indices.size() % _count == 0,
1450 : "The number of dof indices should be cleanly divisible by the variable count");
1451 150389846 : const auto num_shapes = _dof_indices.size() / _count;
1452 1088686722 : for (auto & dof_u : _vector_tags_dof_u)
1453 938296876 : dof_u.resize(num_shapes);
1454 :
1455 347817072 : for (auto & dof_u : _matrix_tags_dof_u)
1456 197427226 : dof_u.resize(num_shapes);
1457 :
1458 150389846 : _has_dof_indices = true;
1459 : }
1460 : else
1461 872329 : _has_dof_indices = false;
1462 : }
1463 : else
1464 946233 : _has_dof_indices = false;
1465 152208408 : }
1466 :
1467 : template <typename OutputType>
1468 : void
1469 6612 : MooseVariableData<OutputType>::reinitNodes(const std::vector<dof_id_type> & nodes)
1470 : {
1471 6612 : _dof_indices.clear();
1472 32917 : for (const auto & node_id : nodes)
1473 : {
1474 26305 : auto && nd = _subproblem.mesh().getMesh().query_node_ptr(node_id);
1475 26305 : if (nd && (_subproblem.mesh().isSemiLocal(const_cast<Node *>(nd))))
1476 : {
1477 26239 : if (nd->n_dofs(_sys.number(), _var_num) > 0)
1478 : {
1479 : if constexpr (std::is_same<RealEigenVector, OutputType>::value)
1480 : {
1481 96 : static thread_local std::vector<dof_id_type> dof_indices;
1482 96 : _dof_map.array_dof_indices(nd, dof_indices, _var_num);
1483 288 : for (const auto dof : dof_indices)
1484 192 : _dof_indices.push_back(dof);
1485 : }
1486 : else
1487 : {
1488 25105 : dof_id_type dof = nd->dof_number(_sys.number(), _var_num, 0);
1489 25105 : _dof_indices.push_back(dof);
1490 : }
1491 : }
1492 : }
1493 : }
1494 :
1495 6612 : if (!_dof_indices.empty())
1496 6408 : _has_dof_indices = true;
1497 : else
1498 204 : _has_dof_indices = false;
1499 6612 : }
1500 :
1501 : template class MooseVariableData<Real>;
1502 : template class MooseVariableData<RealVectorValue>;
1503 : template class MooseVariableData<RealEigenVector>;
|
