deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_values_views_internal.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2023 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/array_view.h>
14#include <deal.II/base/numbers.h>
15
16#include <deal.II/differentiation/ad.h>
17
18#include <deal.II/fe/fe_values_views_internal.h>
19
20#ifdef DEAL_II_WITH_ADOLC
21# include <adolc/adouble.h>
22# include <adolc/adtl.h>
23#endif
24
25#include <type_traits>
26
27DEAL_II_NAMESPACE_OPEN
28
29namespace FEValuesViews
30{
31 namespace internal
32 {
33 namespace
34 {
35 // Check to see if a DoF value is zero, implying that subsequent
36 // operations with the value have no effect.
37 template <typename Number, typename T = void>
38 struct CheckForZero
39 {
40 static bool
41 value(const Number &value)
42 {
43 return value == ::internal::NumberType<Number>::value(0.0);
44 }
45 };
46
47 // For auto-differentiable numbers, the fact that a DoF value is zero
48 // does not imply that its derivatives are zero as well. So we
49 // can't filter by value for these number types.
50 // Note that we also want to avoid actually checking the value itself,
51 // since some AD numbers are not contextually convertible to booleans.
52 template <typename Number>
53 struct CheckForZero<
54 Number,
55 std::enable_if_t<Differentiation::AD::is_ad_number<Number>::value>>
56 {
57 static bool
58 value(const Number & /*value*/)
59 {
60 return false;
61 }
62 };
63 } // namespace
64
65 template <int dim, int spacedim, typename Number>
66 void
67 do_function_values(
68 const ArrayView<const Number> &dof_values,
69 const Table<2, double> &shape_values,
70 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
71 &shape_function_data,
72 std::vector<typename ProductType<Number, double>::type> &values)
73 {
74 const unsigned int dofs_per_cell = dof_values.size();
75 const unsigned int n_quadrature_points = values.size();
76
77 std::fill(values.begin(),
78 values.end(),
79 ::internal::NumberType<Number>::value(0.0));
80
81 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
82 ++shape_function)
83 if (shape_function_data[shape_function]
84 .is_nonzero_shape_function_component)
85 {
86 const Number &value = dof_values[shape_function];
87 // For auto-differentiable numbers, the fact that a DoF value is
88 // zero does not imply that its derivatives are zero as well. So we
89 // can't filter by value for these number types.
90 if (CheckForZero<Number>::value(value) == true)
91 continue;
92
93 const double *shape_value_ptr =
94 &shape_values(shape_function_data[shape_function].row_index, 0);
95 for (unsigned int q_point = 0; q_point < n_quadrature_points;
96 ++q_point, ++shape_value_ptr)
97 values[q_point] += value * (*shape_value_ptr);
98 }
99 }
100
101
102
103 template <int order, int dim, int spacedim, typename Number>
104 void
105 do_function_derivatives(
106 const ArrayView<const Number> &dof_values,
107 const Table<2, ::Tensor<order, spacedim>> &shape_derivatives,
108 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
109 &shape_function_data,
110 std::vector<
111 typename ProductType<Number, ::Tensor<order, spacedim>>::type>
112 &derivatives)
113 {
114 const unsigned int dofs_per_cell = dof_values.size();
115 const unsigned int n_quadrature_points = derivatives.size();
116
117 std::fill(
118 derivatives.begin(),
119 derivatives.end(),
120 typename ProductType<Number, ::Tensor<order, spacedim>>::type());
121
122 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
123 ++shape_function)
124 if (shape_function_data[shape_function]
125 .is_nonzero_shape_function_component)
126 {
127 const Number &value = dof_values[shape_function];
128 // For auto-differentiable numbers, the fact that a DoF value is
129 // zero does not imply that its derivatives are zero as well. So we
130 // can't filter by value for these number types.
131 if (CheckForZero<Number>::value(value) == true)
132 continue;
133
134 const ::Tensor<order, spacedim> *shape_derivative_ptr =
135 &shape_derivatives[shape_function_data[shape_function].row_index]
136 [0];
137 for (unsigned int q_point = 0; q_point < n_quadrature_points;
138 ++q_point)
139 derivatives[q_point] += value * (*shape_derivative_ptr++);
140 }
141 }
142
143
144
145 template <int dim, int spacedim, typename Number>
146 void
147 do_function_laplacians(
148 const ArrayView<const Number> &dof_values,
149 const Table<2, ::Tensor<2, spacedim>> &shape_hessians,
150 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
151 &shape_function_data,
152 std::vector<typename Scalar<dim, spacedim>::
153 template solution_laplacian_type<Number>> &laplacians)
154 {
155 const unsigned int dofs_per_cell = dof_values.size();
156 const unsigned int n_quadrature_points = laplacians.size();
157
158 std::fill(
159 laplacians.begin(),
160 laplacians.end(),
161 typename Scalar<dim,
162 spacedim>::template solution_laplacian_type<Number>());
163
164 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
165 ++shape_function)
166 if (shape_function_data[shape_function]
167 .is_nonzero_shape_function_component)
168 {
169 const Number &value = dof_values[shape_function];
170 // For auto-differentiable numbers, the fact that a DoF value is
171 // zero does not imply that its derivatives are zero as well. So we
172 // can't filter by value for these number types.
173 if (CheckForZero<Number>::value(value) == true)
174 continue;
175
176 const ::Tensor<2, spacedim> *shape_hessian_ptr =
177 &shape_hessians[shape_function_data[shape_function].row_index][0];
178 for (unsigned int q_point = 0; q_point < n_quadrature_points;
179 ++q_point)
180 laplacians[q_point] += value * trace(*shape_hessian_ptr++);
181 }
182 }
183
184
185
186 // ----------------------------- vector part ---------------------------
187
188 template <int dim, int spacedim, typename Number>
189 void
190 do_function_values(
191 const ArrayView<const Number> &dof_values,
192 const Table<2, double> &shape_values,
193 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
194 &shape_function_data,
195 std::vector<
196 typename ProductType<Number, ::Tensor<1, spacedim>>::type>
197 &values)
198 {
199 const unsigned int dofs_per_cell = dof_values.size();
200 const unsigned int n_quadrature_points = values.size();
201
202 std::fill(
203 values.begin(),
204 values.end(),
205 typename ProductType<Number, ::Tensor<1, spacedim>>::type());
206
207 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
208 ++shape_function)
209 {
210 const int snc =
211 shape_function_data[shape_function].single_nonzero_component;
212
213 if (snc == -2)
214 // shape function is zero for the selected components
215 continue;
216
217 const Number &value = dof_values[shape_function];
218 // For auto-differentiable numbers, the fact that a DoF value is zero
219 // does not imply that its derivatives are zero as well. So we
220 // can't filter by value for these number types.
221 if (CheckForZero<Number>::value(value) == true)
222 continue;
223
224 if (snc != -1)
225 {
226 const unsigned int comp = shape_function_data[shape_function]
227 .single_nonzero_component_index;
228 const double *shape_value_ptr = &shape_values(snc, 0);
229 for (unsigned int q_point = 0; q_point < n_quadrature_points;
230 ++q_point, ++shape_value_ptr)
231 values[q_point][comp] += value * (*shape_value_ptr);
232 }
233 else
234 for (unsigned int d = 0; d < spacedim; ++d)
235 if (shape_function_data[shape_function]
236 .is_nonzero_shape_function_component[d])
237 {
238 const double *shape_value_ptr = &shape_values(
239 shape_function_data[shape_function].row_index[d], 0);
240 for (unsigned int q_point = 0; q_point < n_quadrature_points;
241 ++q_point, ++shape_value_ptr)
242 values[q_point][d] += value * (*shape_value_ptr);
243 }
244 }
245 }
246
247
248
249 template <int order, int dim, int spacedim, typename Number>
250 void
251 do_function_derivatives(
252 const ArrayView<const Number> &dof_values,
253 const Table<2, ::Tensor<order, spacedim>> &shape_derivatives,
254 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
255 &shape_function_data,
256 std::vector<
257 typename ProductType<Number, ::Tensor<order + 1, spacedim>>::type>
258 &derivatives)
259 {
260 const unsigned int dofs_per_cell = dof_values.size();
261 const unsigned int n_quadrature_points = derivatives.size();
262
263 std::fill(
264 derivatives.begin(),
265 derivatives.end(),
266 typename ProductType<Number,
267 ::Tensor<order + 1, spacedim>>::type());
268
269 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
270 ++shape_function)
271 {
272 const int snc =
273 shape_function_data[shape_function].single_nonzero_component;
274
275 if (snc == -2)
276 // shape function is zero for the selected components
277 continue;
278
279 const Number &value = dof_values[shape_function];
280 // For auto-differentiable numbers, the fact that a DoF value is zero
281 // does not imply that its derivatives are zero as well. So we
282 // can't filter by value for these number types.
283 if (CheckForZero<Number>::value(value) == true)
284 continue;
285
286 if (snc != -1)
287 {
288 const unsigned int comp = shape_function_data[shape_function]
289 .single_nonzero_component_index;
290 const ::Tensor<order, spacedim> *shape_derivative_ptr =
291 &shape_derivatives[snc][0];
292 for (unsigned int q_point = 0; q_point < n_quadrature_points;
293 ++q_point)
294 derivatives[q_point][comp] += value * (*shape_derivative_ptr++);
295 }
296 else
297 for (unsigned int d = 0; d < spacedim; ++d)
298 if (shape_function_data[shape_function]
299 .is_nonzero_shape_function_component[d])
300 {
301 const ::Tensor<order, spacedim> *shape_derivative_ptr =
302 &shape_derivatives[shape_function_data[shape_function]
303 .row_index[d]][0];
304 for (unsigned int q_point = 0; q_point < n_quadrature_points;
305 ++q_point)
306 derivatives[q_point][d] +=
307 value * (*shape_derivative_ptr++);
308 }
309 }
310 }
311
312
313
314 template <int dim, int spacedim, typename Number>
315 void
316 do_function_symmetric_gradients(
317 const ArrayView<const Number> &dof_values,
318 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
319 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
320 &shape_function_data,
321 std::vector<
322 typename ProductType<Number,
323 ::SymmetricTensor<2, spacedim>>::type>
324 &symmetric_gradients)
325 {
326 const unsigned int dofs_per_cell = dof_values.size();
327 const unsigned int n_quadrature_points = symmetric_gradients.size();
328
329 std::fill(
330 symmetric_gradients.begin(),
331 symmetric_gradients.end(),
332 typename ProductType<Number,
333 ::SymmetricTensor<2, spacedim>>::type());
334
335 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
336 ++shape_function)
337 {
338 const int snc =
339 shape_function_data[shape_function].single_nonzero_component;
340
341 if (snc == -2)
342 // shape function is zero for the selected components
343 continue;
344
345 const Number &value = dof_values[shape_function];
346 // For auto-differentiable numbers, the fact that a DoF value is zero
347 // does not imply that its derivatives are zero as well. So we
348 // can't filter by value for these number types.
349 if (CheckForZero<Number>::value(value) == true)
350 continue;
351
352 if (snc != -1)
353 {
354 const unsigned int comp = shape_function_data[shape_function]
355 .single_nonzero_component_index;
356 const ::Tensor<1, spacedim> *shape_gradient_ptr =
357 &shape_gradients[snc][0];
358 for (unsigned int q_point = 0; q_point < n_quadrature_points;
359 ++q_point)
360 {
361 for (unsigned int d = 0; d < dim; ++d)
362 symmetric_gradients[q_point][comp][d] +=
363 0.5 * value * (*shape_gradient_ptr)[d];
364 symmetric_gradients[q_point][comp][comp] +=
365 0.5 * value * (*shape_gradient_ptr++)[comp];
366 }
367 }
368 else
369 for (unsigned int q_point = 0; q_point < n_quadrature_points;
370 ++q_point)
371 {
372 typename ProductType<Number, ::Tensor<2, spacedim>>::type
373 grad;
374 for (unsigned int d = 0; d < spacedim; ++d)
375 if (shape_function_data[shape_function]
376 .is_nonzero_shape_function_component[d])
377 grad[d] =
378 value *
379 shape_gradients[shape_function_data[shape_function]
380 .row_index[d]][q_point];
381 symmetric_gradients[q_point] += symmetrize(grad);
382 }
383 }
384 }
385
386
387
388 template <int dim, int spacedim, typename Number>
389 void
390 do_function_divergences(
391 const ArrayView<const Number> &dof_values,
392 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
393 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
394 &shape_function_data,
395 std::vector<typename Vector<dim, spacedim>::
396 template solution_divergence_type<Number>> &divergences)
397 {
398 const unsigned int dofs_per_cell = dof_values.size();
399 const unsigned int n_quadrature_points = divergences.size();
400
401 std::fill(
402 divergences.begin(),
403 divergences.end(),
404 typename Vector<dim,
405 spacedim>::template solution_divergence_type<Number>());
406
407 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
408 ++shape_function)
409 {
410 const int snc =
411 shape_function_data[shape_function].single_nonzero_component;
412
413 if (snc == -2)
414 // shape function is zero for the selected components
415 continue;
416
417 const Number &value = dof_values[shape_function];
418 // For auto-differentiable numbers, the fact that a DoF value is zero
419 // does not imply that its derivatives are zero as well. So we
420 // can't filter by value for these number types.
421 if (CheckForZero<Number>::value(value) == true)
422 continue;
423
424 if (snc != -1)
425 {
426 const unsigned int comp = shape_function_data[shape_function]
427 .single_nonzero_component_index;
428 const ::Tensor<1, spacedim> *shape_gradient_ptr =
429 &shape_gradients[snc][0];
430 for (unsigned int q_point = 0; q_point < n_quadrature_points;
431 ++q_point)
432 divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
433 }
434 else
435 for (unsigned int d = 0; d < spacedim; ++d)
436 if (shape_function_data[shape_function]
437 .is_nonzero_shape_function_component[d])
438 {
439 const ::Tensor<1, spacedim> *shape_gradient_ptr =
440 &shape_gradients[shape_function_data[shape_function]
441 .row_index[d]][0];
442 for (unsigned int q_point = 0; q_point < n_quadrature_points;
443 ++q_point)
444 divergences[q_point] += value * (*shape_gradient_ptr++)[d];
445 }
446 }
447 }
448
449
450
451 template <int dim, int spacedim, typename Number>
452 void
453 do_function_curls(
454 const ArrayView<const Number> &dof_values,
455 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
456 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
457 &shape_function_data,
458 std::vector<typename ProductType<
459 Number,
460 typename ::internal::CurlType<spacedim>>::type> &curls)
461 {
462 const unsigned int dofs_per_cell = dof_values.size();
463 const unsigned int n_quadrature_points = curls.size();
464
465 std::fill(curls.begin(),
466 curls.end(),
467 typename ProductType<
468 Number,
469 typename ::internal::CurlType<spacedim>>::type());
470
471 if constexpr (spacedim == 1)
472 {
473 Assert(false,
474 ExcMessage(
475 "Computing the curl in 1d is not a useful operation"));
476 }
477 else if constexpr (spacedim == 2)
478 {
479 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
480 ++shape_function)
481 {
482 const int snc =
483 shape_function_data[shape_function].single_nonzero_component;
484
485 if (snc == -2)
486 // shape function is zero for the selected components
487 continue;
488
489 const Number &value = dof_values[shape_function];
490 // For auto-differentiable numbers, the fact that a DoF value
491 // is zero does not imply that its derivatives are zero as
492 // well. So we can't filter by value for these number types.
493 if (CheckForZero<Number>::value(value) == true)
494 continue;
495
496 if (snc != -1)
497 {
498 const ::Tensor<1, spacedim> *shape_gradient_ptr =
499 &shape_gradients[snc][0];
500
501 Assert(shape_function_data[shape_function]
502 .single_nonzero_component >= 0,
503 ExcInternalError());
504 // we're in 2d, so the formula for the curl is simple:
505 if (shape_function_data[shape_function]
506 .single_nonzero_component_index == 0)
507 for (unsigned int q_point = 0;
508 q_point < n_quadrature_points;
509 ++q_point)
510 curls[q_point] -= value * (*shape_gradient_ptr++)[1];
511 else
512 for (unsigned int q_point = 0;
513 q_point < n_quadrature_points;
514 ++q_point)
515 curls[q_point] += value * (*shape_gradient_ptr++)[0];
516 }
517 else
518 // we have multiple non-zero components in the shape
519 // functions. not all of them must necessarily be within the
520 // 2-component window this FEValuesViews::Vector object
521 // considers, however.
522 {
523 if (shape_function_data[shape_function]
524 .is_nonzero_shape_function_component[0])
525 {
526 const ::Tensor<1, spacedim> *shape_gradient_ptr =
527 &shape_gradients[shape_function_data[shape_function]
528 .row_index[0]][0];
529
530 for (unsigned int q_point = 0;
531 q_point < n_quadrature_points;
532 ++q_point)
533 curls[q_point] -= value * (*shape_gradient_ptr++)[1];
534 }
535
536 if (shape_function_data[shape_function]
537 .is_nonzero_shape_function_component[1])
538 {
539 const ::Tensor<1, spacedim> *shape_gradient_ptr =
540 &shape_gradients[shape_function_data[shape_function]
541 .row_index[1]][0];
542
543 for (unsigned int q_point = 0;
544 q_point < n_quadrature_points;
545 ++q_point)
546 curls[q_point] += value * (*shape_gradient_ptr++)[0];
547 }
548 }
549 }
550 }
551 else if constexpr (spacedim == 3)
552 {
553 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
554 ++shape_function)
555 {
556 const int snc =
557 shape_function_data[shape_function].single_nonzero_component;
558
559 if (snc == -2)
560 // shape function is zero for the selected components
561 continue;
562
563 const Number &value = dof_values[shape_function];
564 // For auto-differentiable numbers, the fact that a DoF value
565 // is zero does not imply that its derivatives are zero as
566 // well. So we can't filter by value for these number types.
567 if (CheckForZero<Number>::value(value) == true)
568 continue;
569
570 if (snc != -1)
571 {
572 const ::Tensor<1, spacedim> *shape_gradient_ptr =
573 &shape_gradients[snc][0];
574
575 switch (shape_function_data[shape_function]
576 .single_nonzero_component_index)
577 {
578 case 0:
579 {
580 for (unsigned int q_point = 0;
581 q_point < n_quadrature_points;
582 ++q_point)
583 {
584 curls[q_point][1] +=
585 value * (*shape_gradient_ptr)[2];
586 curls[q_point][2] -=
587 value * (*shape_gradient_ptr++)[1];
588 }
589
590 break;
591 }
592
593 case 1:
594 {
595 for (unsigned int q_point = 0;
596 q_point < n_quadrature_points;
597 ++q_point)
598 {
599 curls[q_point][0] -=
600 value * (*shape_gradient_ptr)[2];
601 curls[q_point][2] +=
602 value * (*shape_gradient_ptr++)[0];
603 }
604
605 break;
606 }
607
608 case 2:
609 {
610 for (unsigned int q_point = 0;
611 q_point < n_quadrature_points;
612 ++q_point)
613 {
614 curls[q_point][0] +=
615 value * (*shape_gradient_ptr)[1];
616 curls[q_point][1] -=
617 value * (*shape_gradient_ptr++)[0];
618 }
619 break;
620 }
621
622 default:
623 DEAL_II_ASSERT_UNREACHABLE();
624 }
625 }
626 else
627 // we have multiple non-zero components in the shape
628 // functions. not all of them must necessarily be within the
629 // 3-component window this FEValuesViews::Vector object
630 // considers, however.
631 {
632 if (shape_function_data[shape_function]
633 .is_nonzero_shape_function_component[0])
634 {
635 const ::Tensor<1, spacedim> *shape_gradient_ptr =
636 &shape_gradients[shape_function_data[shape_function]
637 .row_index[0]][0];
638
639 for (unsigned int q_point = 0;
640 q_point < n_quadrature_points;
641 ++q_point)
642 {
643 curls[q_point][1] += value * (*shape_gradient_ptr)[2];
644 curls[q_point][2] -=
645 value * (*shape_gradient_ptr++)[1];
646 }
647 }
648
649 if (shape_function_data[shape_function]
650 .is_nonzero_shape_function_component[1])
651 {
652 const ::Tensor<1, spacedim> *shape_gradient_ptr =
653 &shape_gradients[shape_function_data[shape_function]
654 .row_index[1]][0];
655
656 for (unsigned int q_point = 0;
657 q_point < n_quadrature_points;
658 ++q_point)
659 {
660 curls[q_point][0] -= value * (*shape_gradient_ptr)[2];
661 curls[q_point][2] +=
662 value * (*shape_gradient_ptr++)[0];
663 }
664 }
665
666 if (shape_function_data[shape_function]
667 .is_nonzero_shape_function_component[2])
668 {
669 const ::Tensor<1, spacedim> *shape_gradient_ptr =
670 &shape_gradients[shape_function_data[shape_function]
671 .row_index[2]][0];
672
673 for (unsigned int q_point = 0;
674 q_point < n_quadrature_points;
675 ++q_point)
676 {
677 curls[q_point][0] += value * (*shape_gradient_ptr)[1];
678 curls[q_point][1] -=
679 value * (*shape_gradient_ptr++)[0];
680 }
681 }
682 }
683 }
684 }
685 else
686 DEAL_II_ASSERT_UNREACHABLE();
687 }
688
689
690
691 template <int dim, int spacedim, typename Number>
692 void
693 do_function_laplacians(
694 const ArrayView<const Number> &dof_values,
695 const Table<2, ::Tensor<2, spacedim>> &shape_hessians,
696 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
697 &shape_function_data,
698 std::vector<typename Vector<dim, spacedim>::
699 template solution_laplacian_type<Number>> &laplacians)
700 {
701 const unsigned int dofs_per_cell = dof_values.size();
702 const unsigned int n_quadrature_points = laplacians.size();
703
704 std::fill(
705 laplacians.begin(),
706 laplacians.end(),
707 typename Vector<dim,
708 spacedim>::template solution_laplacian_type<Number>());
709
710 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
711 ++shape_function)
712 {
713 const int snc =
714 shape_function_data[shape_function].single_nonzero_component;
715
716 if (snc == -2)
717 // shape function is zero for the selected components
718 continue;
719
720 const Number &value = dof_values[shape_function];
721 // For auto-differentiable numbers, the fact that a DoF value is zero
722 // does not imply that its derivatives are zero as well. So we
723 // can't filter by value for these number types.
724 if (CheckForZero<Number>::value(value) == true)
725 continue;
726
727 if (snc != -1)
728 {
729 const unsigned int comp = shape_function_data[shape_function]
730 .single_nonzero_component_index;
731 const ::Tensor<2, spacedim> *shape_hessian_ptr =
732 &shape_hessians[snc][0];
733 for (unsigned int q_point = 0; q_point < n_quadrature_points;
734 ++q_point)
735 laplacians[q_point][comp] +=
736 value * trace(*shape_hessian_ptr++);
737 }
738 else
739 for (unsigned int d = 0; d < spacedim; ++d)
740 if (shape_function_data[shape_function]
741 .is_nonzero_shape_function_component[d])
742 {
743 const ::Tensor<2, spacedim> *shape_hessian_ptr =
744 &shape_hessians[shape_function_data[shape_function]
745 .row_index[d]][0];
746 for (unsigned int q_point = 0; q_point < n_quadrature_points;
747 ++q_point)
748 laplacians[q_point][d] +=
749 value * trace(*shape_hessian_ptr++);
750 }
751 }
752 }
753
754
755
756 // ---------------------- symmetric tensor part ------------------------
757
758 template <int dim, int spacedim, typename Number>
759 void
760 do_function_values(
761 const ArrayView<const Number> &dof_values,
762 const ::Table<2, double> &shape_values,
763 const std::vector<
764 typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
765 &shape_function_data,
766 std::vector<
767 typename ProductType<Number,
768 ::SymmetricTensor<2, spacedim>>::type>
769 &values)
770 {
771 const unsigned int dofs_per_cell = dof_values.size();
772 const unsigned int n_quadrature_points = values.size();
773
774 std::fill(
775 values.begin(),
776 values.end(),
777 typename ProductType<Number,
778 ::SymmetricTensor<2, spacedim>>::type());
779
780 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
781 ++shape_function)
782 {
783 const int snc =
784 shape_function_data[shape_function].single_nonzero_component;
785
786 if (snc == -2)
787 // shape function is zero for the selected components
788 continue;
789
790 const Number &value = dof_values[shape_function];
791 // For auto-differentiable numbers, the fact that a DoF value is zero
792 // does not imply that its derivatives are zero as well. So we
793 // can't filter by value for these number types.
794 if (CheckForZero<Number>::value(value) == true)
795 continue;
796
797 if (snc != -1)
798 {
799 const TableIndices<2> comp = ::
800 SymmetricTensor<2, spacedim>::unrolled_to_component_indices(
801 shape_function_data[shape_function]
802 .single_nonzero_component_index);
803 const double *shape_value_ptr = &shape_values(snc, 0);
804 for (unsigned int q_point = 0; q_point < n_quadrature_points;
805 ++q_point, ++shape_value_ptr)
806 values[q_point][comp] += value * (*shape_value_ptr);
807 }
808 else
809 for (unsigned int d = 0;
810 d <
811 ::SymmetricTensor<2, spacedim>::n_independent_components;
812 ++d)
813 if (shape_function_data[shape_function]
814 .is_nonzero_shape_function_component[d])
815 {
816 const TableIndices<2> comp =
817 ::SymmetricTensor<2, spacedim>::
818 unrolled_to_component_indices(d);
819 const double *shape_value_ptr = &shape_values(
820 shape_function_data[shape_function].row_index[d], 0);
821 for (unsigned int q_point = 0; q_point < n_quadrature_points;
822 ++q_point, ++shape_value_ptr)
823 values[q_point][comp] += value * (*shape_value_ptr);
824 }
825 }
826 }
827
828
829
830 template <int dim, int spacedim, typename Number>
831 void
832 do_function_divergences(
833 const ArrayView<const Number> &dof_values,
834 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
835 const std::vector<
836 typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
837 &shape_function_data,
838 std::vector<typename SymmetricTensor<2, dim, spacedim>::
839 template solution_divergence_type<Number>> &divergences)
840 {
841 const unsigned int dofs_per_cell = dof_values.size();
842 const unsigned int n_quadrature_points = divergences.size();
843
844 std::fill(divergences.begin(),
845 divergences.end(),
846 typename SymmetricTensor<2, dim, spacedim>::
847 template solution_divergence_type<Number>());
848
849 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
850 ++shape_function)
851 {
852 const int snc =
853 shape_function_data[shape_function].single_nonzero_component;
854
855 if (snc == -2)
856 // shape function is zero for the selected components
857 continue;
858
859 const Number &value = dof_values[shape_function];
860 // For auto-differentiable numbers, the fact that a DoF value is zero
861 // does not imply that its derivatives are zero as well. So we
862 // can't filter by value for these number types.
863 if (CheckForZero<Number>::value(value) == true)
864 continue;
865
866 if (snc != -1)
867 {
868 const unsigned int comp = shape_function_data[shape_function]
869 .single_nonzero_component_index;
870
871 const ::Tensor<1, spacedim> *shape_gradient_ptr =
872 &shape_gradients[snc][0];
873
874 const unsigned int ii = ::SymmetricTensor<2, spacedim>::
875 unrolled_to_component_indices(comp)[0];
876 const unsigned int jj = ::SymmetricTensor<2, spacedim>::
877 unrolled_to_component_indices(comp)[1];
878
879 for (unsigned int q_point = 0; q_point < n_quadrature_points;
880 ++q_point, ++shape_gradient_ptr)
881 {
882 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
883
884 if (ii != jj)
885 divergences[q_point][jj] +=
886 value * (*shape_gradient_ptr)[ii];
887 }
888 }
889 else
890 {
891 for (unsigned int d = 0;
892 d <
893 ::SymmetricTensor<2,
894 spacedim>::n_independent_components;
895 ++d)
896 if (shape_function_data[shape_function]
897 .is_nonzero_shape_function_component[d])
898 {
899 DEAL_II_NOT_IMPLEMENTED();
900
901 // the following implementation needs to be looked over -- I
902 // think it can't be right, because we are in a case where
903 // there is no single nonzero component
904 //
905 // the following is not implemented! we need to consider the
906 // interplay between multiple non-zero entries in shape
907 // function and the representation as a symmetric
908 // second-order tensor
909 const unsigned int comp =
910 shape_function_data[shape_function]
911 .single_nonzero_component_index;
912
913 const ::Tensor<1, spacedim> *shape_gradient_ptr =
914 &shape_gradients[shape_function_data[shape_function]
915 .row_index[d]][0];
916 for (unsigned int q_point = 0;
917 q_point < n_quadrature_points;
918 ++q_point, ++shape_gradient_ptr)
919 {
920 for (unsigned int j = 0; j < spacedim;
921 ++j, ++shape_gradient_ptr)
922 {
923 const unsigned int vector_component =
924 ::SymmetricTensor<2, spacedim>::
925 component_to_unrolled_index(
926 TableIndices<2>(comp, j));
927 divergences[q_point][vector_component] +=
928 value * (*shape_gradient_ptr)[j];
929 }
930 }
931 }
932 }
933 }
934 }
935
936 // ---------------------- non-symmetric tensor part ------------------------
937
938 template <int dim, int spacedim, typename Number>
939 void
940 do_function_values(
941 const ArrayView<const Number> &dof_values,
942 const ::Table<2, double> &shape_values,
943 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
944 &shape_function_data,
945 std::vector<
946 typename ProductType<Number, ::Tensor<2, spacedim>>::type>
947 &values)
948 {
949 const unsigned int dofs_per_cell = dof_values.size();
950 const unsigned int n_quadrature_points = values.size();
951
952 std::fill(
953 values.begin(),
954 values.end(),
955 typename ProductType<Number, ::Tensor<2, spacedim>>::type());
956
957 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
958 ++shape_function)
959 {
960 const int snc =
961 shape_function_data[shape_function].single_nonzero_component;
962
963 if (snc == -2)
964 // shape function is zero for the selected components
965 continue;
966
967 const Number &value = dof_values[shape_function];
968 // For auto-differentiable numbers, the fact that a DoF value is zero
969 // does not imply that its derivatives are zero as well. So we
970 // can't filter by value for these number types.
971 if (CheckForZero<Number>::value(value) == true)
972 continue;
973
974 if (snc != -1)
975 {
976 const unsigned int comp = shape_function_data[shape_function]
977 .single_nonzero_component_index;
978
979 const TableIndices<2> indices =
980 ::Tensor<2, spacedim>::unrolled_to_component_indices(
981 comp);
982
983 const double *shape_value_ptr = &shape_values(snc, 0);
984 for (unsigned int q_point = 0; q_point < n_quadrature_points;
985 ++q_point, ++shape_value_ptr)
986 values[q_point][indices] += value * (*shape_value_ptr);
987 }
988 else
989 for (unsigned int d = 0; d < dim * dim; ++d)
990 if (shape_function_data[shape_function]
991 .is_nonzero_shape_function_component[d])
992 {
993 const TableIndices<2> indices =
994 ::Tensor<2, spacedim>::unrolled_to_component_indices(
995 d);
996
997 const double *shape_value_ptr = &shape_values(
998 shape_function_data[shape_function].row_index[d], 0);
999 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1000 ++q_point, ++shape_value_ptr)
1001 values[q_point][indices] += value * (*shape_value_ptr);
1002 }
1003 }
1004 }
1005
1006
1007
1008 template <int dim, int spacedim, typename Number>
1009 void
1010 do_function_divergences(
1011 const ArrayView<const Number> &dof_values,
1012 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1013 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1014 &shape_function_data,
1015 std::vector<typename Tensor<2, dim, spacedim>::
1016 template solution_divergence_type<Number>> &divergences)
1017 {
1018 const unsigned int dofs_per_cell = dof_values.size();
1019 const unsigned int n_quadrature_points = divergences.size();
1020
1021 std::fill(
1022 divergences.begin(),
1023 divergences.end(),
1024 typename Tensor<2, dim, spacedim>::template solution_divergence_type<
1025 Number>());
1026
1027 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1028 ++shape_function)
1029 {
1030 const int snc =
1031 shape_function_data[shape_function].single_nonzero_component;
1032
1033 if (snc == -2)
1034 // shape function is zero for the selected components
1035 continue;
1036
1037 const Number &value = dof_values[shape_function];
1038 // For auto-differentiable numbers, the fact that a DoF value is zero
1039 // does not imply that its derivatives are zero as well. So we
1040 // can't filter by value for these number types.
1041 if (CheckForZero<Number>::value(value) == true)
1042 continue;
1043
1044 if (snc != -1)
1045 {
1046 const unsigned int comp = shape_function_data[shape_function]
1047 .single_nonzero_component_index;
1048
1049 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1050 &shape_gradients[snc][0];
1051
1052 const TableIndices<2> indices =
1053 ::Tensor<2, spacedim>::unrolled_to_component_indices(
1054 comp);
1055 const unsigned int ii = indices[0];
1056 const unsigned int jj = indices[1];
1057
1058 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1059 ++q_point, ++shape_gradient_ptr)
1060 {
1061 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1062 }
1063 }
1064 else
1065 {
1066 for (unsigned int d = 0; d < dim * dim; ++d)
1067 if (shape_function_data[shape_function]
1068 .is_nonzero_shape_function_component[d])
1069 {
1070 DEAL_II_NOT_IMPLEMENTED();
1071 }
1072 }
1073 }
1074 }
1075
1076
1077
1078 template <int dim, int spacedim, typename Number>
1079 void
1080 do_function_gradients(
1081 const ArrayView<const Number> &dof_values,
1082 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1083 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1084 &shape_function_data,
1085 std::vector<typename Tensor<2, dim, spacedim>::
1086 template solution_gradient_type<Number>> &gradients)
1087 {
1088 const unsigned int dofs_per_cell = dof_values.size();
1089 const unsigned int n_quadrature_points = gradients.size();
1090
1091 std::fill(
1092 gradients.begin(),
1093 gradients.end(),
1094 typename Tensor<2, dim, spacedim>::template solution_gradient_type<
1095 Number>());
1096
1097 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1098 ++shape_function)
1099 {
1100 const int snc =
1101 shape_function_data[shape_function].single_nonzero_component;
1102
1103 if (snc == -2)
1104 // shape function is zero for the selected components
1105 continue;
1106
1107 const Number &value = dof_values[shape_function];
1108 // For auto-differentiable numbers, the fact that a DoF value is zero
1109 // does not imply that its derivatives are zero as well. So we
1110 // can't filter by value for these number types.
1111 if (CheckForZero<Number>::value(value) == true)
1112 continue;
1113
1114 if (snc != -1)
1115 {
1116 const unsigned int comp = shape_function_data[shape_function]
1117 .single_nonzero_component_index;
1118
1119 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1120 &shape_gradients[snc][0];
1121
1122 const TableIndices<2> indices =
1123 ::Tensor<2, spacedim>::unrolled_to_component_indices(
1124 comp);
1125 const unsigned int ii = indices[0];
1126 const unsigned int jj = indices[1];
1127
1128 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1129 ++q_point, ++shape_gradient_ptr)
1130 {
1131 gradients[q_point][ii][jj] += value * (*shape_gradient_ptr);
1132 }
1133 }
1134 else
1135 {
1136 for (unsigned int d = 0; d < dim * dim; ++d)
1137 if (shape_function_data[shape_function]
1138 .is_nonzero_shape_function_component[d])
1139 {
1140 DEAL_II_NOT_IMPLEMENTED();
1141 }
1142 }
1143 }
1144 }
1145 } // end of namespace internal
1146} // namespace FEValuesViews
1147
1148
1149
1150/*------------------------------- Explicit Instantiations -------------*/
1151
1152#include "fe/fe_values_views_internal.inst"
1153
1154DEAL_II_NAMESPACE_CLOSE
array_view.h
ArrayView::size
std::size_t size() const
Definition array_view.h:737
SymmetricTensor::component_to_unrolled_index
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensor::unrolled_to_component_indices
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Tensor::unrolled_to_component_indices
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exception_macros.h:781
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exception_macros.h:705
fe_values_views_internal.h
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
FEValuesViews::internal::do_function_divergences
void do_function_divergences(const ArrayView< const Number > &dof_values, const Table< 2, ::Tensor< 1, spacedim > > &shape_gradients, const std::vector< typename Vector< dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename Vector< dim, spacedim >::template solution_divergence_type< Number > > &divergences)
Definition fe_values_views_internal.cc:390
FEValuesViews::internal::do_function_symmetric_gradients
void do_function_symmetric_gradients(const ArrayView< const Number > &dof_values, const Table< 2, ::Tensor< 1, spacedim > > &shape_gradients, const std::vector< typename Vector< dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename ProductType< Number, ::SymmetricTensor< 2, spacedim > >::type > &symmetric_gradients)
Definition fe_values_views_internal.cc:316
FEValuesViews::internal::do_function_derivatives
void do_function_derivatives(const ArrayView< const Number > &dof_values, const Table< 2, ::Tensor< order, spacedim > > &shape_derivatives, const std::vector< typename Scalar< dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename ProductType< Number, ::Tensor< order, spacedim > >::type > &derivatives)
Definition fe_values_views_internal.cc:105
FEValuesViews::internal::do_function_gradients
void do_function_gradients(const ArrayView< const Number > &dof_values, const Table< 2, ::Tensor< 1, spacedim > > &shape_gradients, const std::vector< typename Tensor< 2, dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename Tensor< 2, dim, spacedim >::template solution_gradient_type< Number > > &gradients)
Definition fe_values_views_internal.cc:1080
FEValuesViews::internal::do_function_values
void do_function_values(const ArrayView< const Number > &dof_values, const Table< 2, double > &shape_values, const std::vector< typename Scalar< dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename ProductType< Number, double >::type > &values)
Definition fe_values_views_internal.cc:67
FEValuesViews::internal::do_function_curls
void do_function_curls(const ArrayView< const Number > &dof_values, const Table< 2, ::Tensor< 1, spacedim > > &shape_gradients, const std::vector< typename Vector< dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename ProductType< Number, ::internal::CurlType< spacedim > >::type > &curls)
Definition fe_values_views_internal.cc:453
FEValuesViews::internal::do_function_laplacians
void do_function_laplacians(const ArrayView< const Number > &dof_values, const Table< 2, ::Tensor< 2, spacedim > > &shape_hessians, const std::vector< typename Scalar< dim, spacedim >::ShapeFunctionData > &shape_function_data, std::vector< typename Scalar< dim, spacedim >::template solution_laplacian_type< Number > > &laplacians)
Definition fe_values_views_internal.cc:147
PlasticityLab::Number
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
Definition ALE_Finite_Strain_Plasticity.h:8733
std
STL namespace.
ProductType::type
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
Definition template_constraints.h:451
internal::NumberType::value
static constexpr const T & value(const T &t)
Definition numbers.h:662
symmetrize
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition symmetric_tensor.h:3572
trace
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
Generated by doxygen 1.9.8