deal.II version GIT relicensing-2171-g98af52917d 2024-11-23 14:40:00+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
15#include <deal.II/base/array_view.h>
16#include <deal.II/base/numbers.h>
17
18#include <deal.II/differentiation/ad.h>
19
20#include <deal.II/fe/fe_values_views_internal.h>
21
22#include <type_traits>
23
24DEAL_II_NAMESPACE_OPEN
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 {
40 return value == ::internal::NumberType<Number>::value(0.0);
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
64 do_function_values(
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(),
76 ::internal::NumberType<Number>::value(0.0));
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, ++shape_value_ptr)
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
102 do_function_derivatives(
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
144 do_function_laplacians(
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
187 do_function_values(
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, ++shape_value_ptr)
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, ++shape_value_ptr)
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
248 do_function_derivatives(
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,
264 ::Tensor<order + 1, spacedim>>::type());
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
313 do_function_symmetric_gradients(
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,
320 ::SymmetricTensor<2, spacedim>>::type>
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,
330 ::SymmetricTensor<2, spacedim>>::type());
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 {
369 typename ProductType<Number, ::Tensor<2, spacedim>>::type
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
387 do_function_divergences(
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
450 do_function_curls(
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,
473 ExcMessage(
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,
505 ExcInternalError());
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:
634 DEAL_II_ASSERT_UNREACHABLE();
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
710 do_function_laplacians(
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
777 do_function_values(
778 const ArrayView<const Number> &dof_values,
779 const ::Table<2, double> &shape_values,
780 const std::vector<
781 typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
782 &shape_function_data,
783 std::vector<
784 typename ProductType<Number,
785 ::SymmetricTensor<2, spacedim>>::type>
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,
795 ::SymmetricTensor<2, spacedim>>::type());
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 = ::
817 SymmetricTensor<2, spacedim>::unrolled_to_component_indices(
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, ++shape_value_ptr)
823 values[q_point][comp] += value * (*shape_value_ptr);
824 }
825 else
826 for (unsigned int d = 0;
827 d <
828 ::SymmetricTensor<2, spacedim>::n_independent_components;
829 ++d)
830 if (shape_function_data[shape_function]
831 .is_nonzero_shape_function_component[d])
832 {
833 const TableIndices<2> comp =
834 ::SymmetricTensor<2, spacedim>::
835 unrolled_to_component_indices(d);
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, ++shape_value_ptr)
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
849 do_function_divergences(
850 const ArrayView<const Number> &dof_values,
851 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
852 const std::vector<
853 typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
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(),
863 typename SymmetricTensor<2, dim, spacedim>::
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>::
892 unrolled_to_component_indices(comp)[0];
893 const unsigned int jj = ::SymmetricTensor<2, spacedim>::
894 unrolled_to_component_indices(comp)[1];
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 <
910 ::SymmetricTensor<2,
911 spacedim>::n_independent_components;
912 ++d)
913 if (shape_function_data[shape_function]
914 .is_nonzero_shape_function_component[d])
915 {
916 DEAL_II_NOT_IMPLEMENTED();
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;
938 ++j, ++shape_gradient_ptr)
939 {
940 const unsigned int vector_component =
941 ::SymmetricTensor<2, spacedim>::
942 component_to_unrolled_index(
943 TableIndices<2>(comp, j));
944 divergences[q_point][vector_component] +=
945 value * (*shape_gradient_ptr)[j];
946 }
947 }
948 }
949 }
950 }
951 }
952
953 // ---------------------- non-symmetric tensor part ------------------------
954
955 template <int dim, int spacedim, typename Number>
956 void
957 do_function_values(
958 const ArrayView<const Number> &dof_values,
959 const ::Table<2, double> &shape_values,
960 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
961 &shape_function_data,
962 std::vector<
963 typename ProductType<Number, ::Tensor<2, spacedim>>::type>
964 &values)
965 {
966 const unsigned int dofs_per_cell = dof_values.size();
967 const unsigned int n_quadrature_points = values.size();
968
969 std::fill(
970 values.begin(),
971 values.end(),
972 typename ProductType<Number, ::Tensor<2, spacedim>>::type());
973
974 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
975 ++shape_function)
976 {
977 const int snc =
978 shape_function_data[shape_function].single_nonzero_component;
979
980 if (snc == -2)
981 // shape function is zero for the selected components
982 continue;
983
984 const Number &value = dof_values[shape_function];
985 // For auto-differentiable numbers, the fact that a DoF value is zero
986 // does not imply that its derivatives are zero as well. So we
987 // can't filter by value for these number types.
988 if (CheckForZero<Number>::value(value) == true)
989 continue;
990
991 if (snc != -1)
992 {
993 const unsigned int comp = shape_function_data[shape_function]
994 .single_nonzero_component_index;
995
996 const TableIndices<2> indices =
997 ::Tensor<2, spacedim>::unrolled_to_component_indices(
998 comp);
999
1000 const double *shape_value_ptr = &shape_values(snc, 0);
1001 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1002 ++q_point, ++shape_value_ptr)
1003 values[q_point][indices] += value * (*shape_value_ptr);
1004 }
1005 else
1006 for (unsigned int d = 0; d < dim * dim; ++d)
1007 if (shape_function_data[shape_function]
1008 .is_nonzero_shape_function_component[d])
1009 {
1010 const TableIndices<2> indices =
1011 ::Tensor<2, spacedim>::unrolled_to_component_indices(
1012 d);
1013
1014 const double *shape_value_ptr = &shape_values(
1015 shape_function_data[shape_function].row_index[d], 0);
1016 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1017 ++q_point, ++shape_value_ptr)
1018 values[q_point][indices] += value * (*shape_value_ptr);
1019 }
1020 }
1021 }
1022
1023
1024
1025 template <int dim, int spacedim, typename Number>
1026 void
1027 do_function_divergences(
1028 const ArrayView<const Number> &dof_values,
1029 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1030 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1031 &shape_function_data,
1032 std::vector<typename Tensor<2, dim, spacedim>::
1033 template solution_divergence_type<Number>> &divergences)
1034 {
1035 const unsigned int dofs_per_cell = dof_values.size();
1036 const unsigned int n_quadrature_points = divergences.size();
1037
1038 std::fill(
1039 divergences.begin(),
1040 divergences.end(),
1041 typename Tensor<2, dim, spacedim>::template solution_divergence_type<
1042 Number>());
1043
1044 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1045 ++shape_function)
1046 {
1047 const int snc =
1048 shape_function_data[shape_function].single_nonzero_component;
1049
1050 if (snc == -2)
1051 // shape function is zero for the selected components
1052 continue;
1053
1054 const Number &value = dof_values[shape_function];
1055 // For auto-differentiable numbers, the fact that a DoF value is zero
1056 // does not imply that its derivatives are zero as well. So we
1057 // can't filter by value for these number types.
1058 if (CheckForZero<Number>::value(value) == true)
1059 continue;
1060
1061 if (snc != -1)
1062 {
1063 const unsigned int comp = shape_function_data[shape_function]
1064 .single_nonzero_component_index;
1065
1066 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1067 &shape_gradients[snc][0];
1068
1069 const TableIndices<2> indices =
1070 ::Tensor<2, spacedim>::unrolled_to_component_indices(
1071 comp);
1072 const unsigned int ii = indices[0];
1073 const unsigned int jj = indices[1];
1074
1075 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1076 ++q_point, ++shape_gradient_ptr)
1077 {
1078 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1079 }
1080 }
1081 else
1082 {
1083 for (unsigned int d = 0; d < dim * dim; ++d)
1084 if (shape_function_data[shape_function]
1085 .is_nonzero_shape_function_component[d])
1086 {
1087 DEAL_II_NOT_IMPLEMENTED();
1088 }
1089 }
1090 }
1091 }
1092
1093
1094
1095 template <int dim, int spacedim, typename Number>
1096 void
1097 do_function_gradients(
1098 const ArrayView<const Number> &dof_values,
1099 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1100 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1101 &shape_function_data,
1102 std::vector<typename Tensor<2, dim, spacedim>::
1103 template solution_gradient_type<Number>> &gradients)
1104 {
1105 const unsigned int dofs_per_cell = dof_values.size();
1106 const unsigned int n_quadrature_points = gradients.size();
1107
1108 std::fill(
1109 gradients.begin(),
1110 gradients.end(),
1111 typename Tensor<2, dim, spacedim>::template solution_gradient_type<
1112 Number>());
1113
1114 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1115 ++shape_function)
1116 {
1117 const int snc =
1118 shape_function_data[shape_function].single_nonzero_component;
1119
1120 if (snc == -2)
1121 // shape function is zero for the selected components
1122 continue;
1123
1124 const Number &value = dof_values[shape_function];
1125 // For auto-differentiable numbers, the fact that a DoF value is zero
1126 // does not imply that its derivatives are zero as well. So we
1127 // can't filter by value for these number types.
1128 if (CheckForZero<Number>::value(value) == true)
1129 continue;
1130
1131 if (snc != -1)
1132 {
1133 const unsigned int comp = shape_function_data[shape_function]
1134 .single_nonzero_component_index;
1135
1136 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1137 &shape_gradients[snc][0];
1138
1139 const TableIndices<2> indices =
1140 ::Tensor<2, spacedim>::unrolled_to_component_indices(
1141 comp);
1142 const unsigned int ii = indices[0];
1143 const unsigned int jj = indices[1];
1144
1145 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1146 ++q_point, ++shape_gradient_ptr)
1147 {
1148 gradients[q_point][ii][jj] += value * (*shape_gradient_ptr);
1149 }
1150 }
1151 else
1152 {
1153 for (unsigned int d = 0; d < dim * dim; ++d)
1154 if (shape_function_data[shape_function]
1155 .is_nonzero_shape_function_component[d])
1156 {
1157 DEAL_II_NOT_IMPLEMENTED();
1158 }
1159 }
1160 }
1161 }
1162 } // end of namespace internal
1163} // namespace FEValuesViews
1164
1165
1166
1167/*------------------------------- Explicit Instantiations -------------*/
1168
1169#include "fe_values_views_internal.inst"
1170
1171DEAL_II_NAMESPACE_CLOSE
array_view.h
ArrayView::size
std::size_t size() const
Definition array_view.h:684
SymmetricTensor::component_to_unrolled_index
static DEAL_II_HOST constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensor::unrolled_to_component_indices
static DEAL_II_HOST 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:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
fe_values_views_internal.h
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exceptions.h:1938
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
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:387
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:313
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:102
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, typename ::internal::CurlType< spacedim >::type >::type > &curls)
Definition fe_values_views_internal.cc:450
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:1097
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:64
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:144
std
STL namespace.
ProductType::type
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
Definition template_constraints.h:462
internal::NumberType::value
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:697
symmetrize
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition symmetric_tensor.h:3448
trace
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
Generated by doxygen 1.9.8