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