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