Reference documentation for deal.II version GIT relicensing-761-gcc3619e707 2024-05-25 17:40: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_base.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 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
22
24
26
27#include <deal.II/fe/fe.h>
29#include <deal.II/fe/mapping.h>
30
33
34#include <deal.II/lac/vector.h>
35
36#include <boost/container/small_vector.hpp>
37
38#include <iomanip>
39#include <memory>
40#include <type_traits>
41
43
44
45namespace internal
46{
47 template <int dim, int spacedim>
48 inline std::vector<unsigned int>
50 {
51 std::vector<unsigned int> shape_function_to_row_table(
53 unsigned int row = 0;
54 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
55 {
56 // loop over all components that are nonzero for this particular
57 // shape function. if a component is zero then we leave the
58 // value in the table unchanged (at the invalid value)
59 // otherwise it is mapped to the next free entry
60 unsigned int nth_nonzero_component = 0;
61 for (unsigned int c = 0; c < fe.n_components(); ++c)
62 if (fe.get_nonzero_components(i)[c] == true)
63 {
64 shape_function_to_row_table[i * fe.n_components() + c] =
65 row + nth_nonzero_component;
66 ++nth_nonzero_component;
67 }
68 row += fe.n_nonzero_components(i);
69 }
70
71 return shape_function_to_row_table;
72 }
73
74 namespace
75 {
76 // Check to see if a DoF value is zero, implying that subsequent operations
77 // with the value have no effect.
78 template <typename Number, typename T = void>
79 struct CheckForZero
80 {
81 static bool
82 value(const Number &value)
83 {
85 }
86 };
87
88 // For auto-differentiable numbers, the fact that a DoF value is zero
89 // does not imply that its derivatives are zero as well. So we
90 // can't filter by value for these number types.
91 // Note that we also want to avoid actually checking the value itself,
92 // since some AD numbers are not contextually convertible to booleans.
93 template <typename Number>
94 struct CheckForZero<
95 Number,
96 std::enable_if_t<Differentiation::AD::is_ad_number<Number>::value>>
97 {
98 static bool
99 value(const Number & /*value*/)
100 {
101 return false;
102 }
103 };
104 } // namespace
105} // namespace internal
106
107/* ------------ FEValuesBase<dim,spacedim>::CellIteratorWrapper ----------- */
108
109
110template <int dim, int spacedim>
115
116
117
118template <int dim, int spacedim>
123
124
125
126template <int dim, int spacedim>
131
132
133
134template <int dim, int spacedim>
135bool
137{
138 return cell.has_value();
139}
140
141
142
143template <int dim, int spacedim>
145operator typename Triangulation<dim, spacedim>::cell_iterator() const
146{
147 Assert(is_initialized(), ExcNotReinited());
148
149 // We can always convert to a tria iterator, regardless of which of
150 // the three types of cell we store.
151 return std::visit(
152 [](auto &cell_iterator) ->
154 return cell_iterator;
155 },
156 cell.value());
157}
158
159
160
161template <int dim, int spacedim>
164{
165 Assert(is_initialized(), ExcNotReinited());
166
167 switch (cell.value().index())
168 {
169 case 1:
170 return std::get<1>(cell.value())->get_dof_handler().n_dofs();
171 case 2:
172 return std::get<2>(cell.value())->get_dof_handler().n_dofs();
173 default:
174 Assert(false, ExcNeedsDoFHandler());
176 }
177}
178
179
180
181template <int dim, int spacedim>
182template <typename Number>
183void
185 const ReadVector<Number> &in,
186 Vector<Number> &out) const
187{
188 Assert(is_initialized(), ExcNotReinited());
189
190 switch (cell.value().index())
191 {
192 case 1:
193 std::get<1>(cell.value())->get_interpolated_dof_values(in, out);
194 break;
195
196 case 2:
197 std::get<2>(cell.value())->get_interpolated_dof_values(in, out);
198 break;
199
200 default:
201 Assert(false, ExcNeedsDoFHandler());
202 break;
203 }
204}
205
206
207
208/*------------------------------- FEValuesBase ---------------------------*/
209
210
211template <int dim, int spacedim>
213 const unsigned int n_q_points,
214 const unsigned int dofs_per_cell,
215 const UpdateFlags flags,
218 : n_quadrature_points(n_q_points)
219 , max_n_quadrature_points(n_q_points)
221 , mapping(&mapping, typeid(*this).name())
222 , fe(&fe, typeid(*this).name())
223 , cell_similarity(CellSimilarity::Similarity::none)
226{
227 Assert(n_q_points > 0,
228 ExcMessage("There is nothing useful you can do with an FEValues "
229 "object when using a quadrature formula with zero "
230 "quadrature points!"));
231 this->update_flags = flags;
232}
233
234
235
236template <int dim, int spacedim>
242
243
244
245template <int dim, int spacedim>
246void
252
253
254
255namespace internal
256{
257 // put shape function part of get_function_xxx methods into separate
258 // internal functions. this allows us to reuse the same code for several
259 // functions (e.g. both the versions with and without indices) as well as
260 // the same code for gradients and Hessians. Moreover, this speeds up
261 // compilation and reduces the size of the final file since all the
262 // different global vectors get channeled through the same code.
263
264 template <typename Number, typename Number2>
265 void
266 do_function_values(const ArrayView<Number2> &dof_values,
267 const ::Table<2, double> &shape_values,
268 std::vector<Number> &values)
269 {
270 // scalar finite elements, so shape_values.size() == dofs_per_cell
271 const unsigned int dofs_per_cell = shape_values.n_rows();
272 const unsigned int n_quadrature_points = values.size();
273
274 // initialize with zero
275 std::fill_n(values.begin(),
276 n_quadrature_points,
278
279 // add up contributions of trial functions. note that here we deal with
280 // scalar finite elements, so no need to check for non-primitivity of
281 // shape functions. in order to increase the speed of this function, we
282 // directly access the data in the shape_values array, and increment
283 // pointers for accessing the data. this saves some lookup time and
284 // indexing. moreover, the order of the loops is such that we can access
285 // the shape_values data stored contiguously
286 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
287 {
288 const Number2 value = dof_values[shape_func];
289 // For auto-differentiable numbers, the fact that a DoF value is zero
290 // does not imply that its derivatives are zero as well. So we
291 // can't filter by value for these number types.
294 continue;
295
296 const double *shape_value_ptr = &shape_values(shape_func, 0);
297 for (unsigned int point = 0; point < n_quadrature_points; ++point)
298 values[point] += value * (*shape_value_ptr++);
299 }
300 }
301
302
303
304 template <int dim, int spacedim, typename VectorType>
305 void
306 do_function_values(
308 const ::Table<2, double> &shape_values,
310 const std::vector<unsigned int> &shape_function_to_row_table,
312 const bool quadrature_points_fastest = false,
313 const unsigned int component_multiple = 1)
314 {
315 using Number = typename VectorType::value_type;
316 // initialize with zero
317 for (unsigned int i = 0; i < values.size(); ++i)
318 std::fill_n(values[i].begin(),
319 values[i].size(),
320 typename VectorType::value_type());
321
322 // see if there the current cell has DoFs at all, and if not
323 // then there is nothing else to do.
324 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
325 if (dofs_per_cell == 0)
326 return;
327
328 const unsigned int n_quadrature_points =
329 quadrature_points_fastest ? values[0].size() : values.size();
330 const unsigned int n_components = fe.n_components();
331
332 // Assert that we can write all components into the result vectors
333 const unsigned result_components = n_components * component_multiple;
334 (void)result_components;
335 if (quadrature_points_fastest)
336 {
337 AssertDimension(values.size(), result_components);
338 for (unsigned int i = 0; i < values.size(); ++i)
339 AssertDimension(values[i].size(), n_quadrature_points);
340 }
341 else
342 {
343 AssertDimension(values.size(), n_quadrature_points);
344 for (unsigned int i = 0; i < values.size(); ++i)
345 AssertDimension(values[i].size(), result_components);
346 }
347
348 // add up contributions of trial functions. now check whether the shape
349 // function is primitive or not. if it is, then set its only non-zero
350 // component, otherwise loop over components
351 for (unsigned int mc = 0; mc < component_multiple; ++mc)
352 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
353 ++shape_func)
354 {
355 const Number &value = dof_values[shape_func + mc * dofs_per_cell];
356 // For auto-differentiable numbers, the fact that a DoF value is zero
357 // does not imply that its derivatives are zero as well. So we
358 // can't filter by value for these number types.
359 if (::internal::CheckForZero<Number>::value(value) == true)
360 continue;
361
362 if (fe.is_primitive(shape_func))
363 {
364 const unsigned int comp =
365 fe.system_to_component_index(shape_func).first +
366 mc * n_components;
367 const unsigned int row =
368 shape_function_to_row_table[shape_func * n_components + comp];
369
370 const double *shape_value_ptr = &shape_values(row, 0);
371
372 if (quadrature_points_fastest)
373 {
374 VectorType &values_comp = values[comp];
375 for (unsigned int point = 0; point < n_quadrature_points;
376 ++point)
377 values_comp[point] += value * (*shape_value_ptr++);
378 }
379 else
380 for (unsigned int point = 0; point < n_quadrature_points;
381 ++point)
382 values[point][comp] += value * (*shape_value_ptr++);
383 }
384 else
385 for (unsigned int c = 0; c < n_components; ++c)
386 {
387 if (fe.get_nonzero_components(shape_func)[c] == false)
388 continue;
389
390 const unsigned int row =
391 shape_function_to_row_table[shape_func * n_components + c];
392
393 const double *shape_value_ptr = &shape_values(row, 0);
394 const unsigned int comp = c + mc * n_components;
395
396 if (quadrature_points_fastest)
397 {
398 VectorType &values_comp = values[comp];
399 for (unsigned int point = 0; point < n_quadrature_points;
400 ++point)
401 values_comp[point] += value * (*shape_value_ptr++);
402 }
403 else
404 for (unsigned int point = 0; point < n_quadrature_points;
405 ++point)
406 values[point][comp] += value * (*shape_value_ptr++);
407 }
408 }
409 }
410
411
412
413 // use the same implementation for gradients and Hessians, distinguish them
414 // by the rank of the tensors
415 template <int order, int spacedim, typename Number>
416 void
417 do_function_derivatives(
418 const ArrayView<Number> &dof_values,
419 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
420 std::vector<Tensor<order, spacedim, Number>> &derivatives)
421 {
422 const unsigned int dofs_per_cell = shape_derivatives.size()[0];
423 const unsigned int n_quadrature_points = derivatives.size();
424
425 // initialize with zero
426 std::fill_n(derivatives.begin(),
427 n_quadrature_points,
429
430 // add up contributions of trial functions. note that here we deal with
431 // scalar finite elements, so no need to check for non-primitivity of
432 // shape functions. in order to increase the speed of this function, we
433 // directly access the data in the shape_gradients/hessians array, and
434 // increment pointers for accessing the data. this saves some lookup time
435 // and indexing. moreover, the order of the loops is such that we can
436 // access the shape_gradients/hessians data stored contiguously
437 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
438 {
439 const Number &value = dof_values[shape_func];
440 // For auto-differentiable numbers, the fact that a DoF value is zero
441 // does not imply that its derivatives are zero as well. So we
442 // can't filter by value for these number types.
443 if (::internal::CheckForZero<Number>::value(value) == true)
444 continue;
445
446 const Tensor<order, spacedim> *shape_derivative_ptr =
447 &shape_derivatives[shape_func][0];
448 for (unsigned int point = 0; point < n_quadrature_points; ++point)
449 derivatives[point] += value * (*shape_derivative_ptr++);
450 }
451 }
452
453
454
455 template <int order, int dim, int spacedim, typename Number>
456 void
457 do_function_derivatives(
458 const ArrayView<Number> &dof_values,
459 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
461 const std::vector<unsigned int> &shape_function_to_row_table,
462 ArrayView<std::vector<Tensor<order, spacedim, Number>>> derivatives,
463 const bool quadrature_points_fastest = false,
464 const unsigned int component_multiple = 1)
465 {
466 // initialize with zero
467 for (unsigned int i = 0; i < derivatives.size(); ++i)
468 std::fill_n(derivatives[i].begin(),
469 derivatives[i].size(),
471
472 // see if there the current cell has DoFs at all, and if not
473 // then there is nothing else to do.
474 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
475 if (dofs_per_cell == 0)
476 return;
477
478
479 const unsigned int n_quadrature_points =
480 quadrature_points_fastest ? derivatives[0].size() : derivatives.size();
481 const unsigned int n_components = fe.n_components();
482
483 // Assert that we can write all components into the result vectors
484 const unsigned result_components = n_components * component_multiple;
485 (void)result_components;
486 if (quadrature_points_fastest)
488 AssertDimension(derivatives.size(), result_components);
489 for (unsigned int i = 0; i < derivatives.size(); ++i)
490 AssertDimension(derivatives[i].size(), n_quadrature_points);
491 }
492 else
493 {
494 AssertDimension(derivatives.size(), n_quadrature_points);
495 for (unsigned int i = 0; i < derivatives.size(); ++i)
496 AssertDimension(derivatives[i].size(), result_components);
497 }
498
499 // add up contributions of trial functions. now check whether the shape
500 // function is primitive or not. if it is, then set its only non-zero
501 // component, otherwise loop over components
502 for (unsigned int mc = 0; mc < component_multiple; ++mc)
503 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
504 ++shape_func)
506 const Number &value = dof_values[shape_func + mc * dofs_per_cell];
507 // For auto-differentiable numbers, the fact that a DoF value is zero
508 // does not imply that its derivatives are zero as well. So we
509 // can't filter by value for these number types.
510 if (::internal::CheckForZero<Number>::value(value) == true)
511 continue;
512
513 if (fe.is_primitive(shape_func))
514 {
515 const unsigned int comp =
516 fe.system_to_component_index(shape_func).first +
517 mc * n_components;
518 const unsigned int row =
519 shape_function_to_row_table[shape_func * n_components + comp];
520
521 const Tensor<order, spacedim> *shape_derivative_ptr =
522 &shape_derivatives[row][0];
523
524 if (quadrature_points_fastest)
525 for (unsigned int point = 0; point < n_quadrature_points;
526 ++point)
527 derivatives[comp][point] += value * (*shape_derivative_ptr++);
528 else
529 for (unsigned int point = 0; point < n_quadrature_points;
530 ++point)
531 derivatives[point][comp] += value * (*shape_derivative_ptr++);
532 }
533 else
534 for (unsigned int c = 0; c < n_components; ++c)
535 {
536 if (fe.get_nonzero_components(shape_func)[c] == false)
537 continue;
538
539 const unsigned int row =
540 shape_function_to_row_table[shape_func * n_components + c];
541
542 const Tensor<order, spacedim> *shape_derivative_ptr =
543 &shape_derivatives[row][0];
544 const unsigned int comp = c + mc * n_components;
545
546 if (quadrature_points_fastest)
547 for (unsigned int point = 0; point < n_quadrature_points;
548 ++point)
549 derivatives[comp][point] +=
550 value * (*shape_derivative_ptr++);
551 else
552 for (unsigned int point = 0; point < n_quadrature_points;
553 ++point)
554 derivatives[point][comp] +=
555 value * (*shape_derivative_ptr++);
556 }
557 }
558 }
559
560
561
562 template <int spacedim, typename Number, typename Number2>
563 void
564 do_function_laplacians(
565 const ArrayView<Number2> &dof_values,
566 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
567 std::vector<Number> &laplacians)
568 {
569 const unsigned int dofs_per_cell = shape_hessians.size()[0];
570 const unsigned int n_quadrature_points = laplacians.size();
571
572 // initialize with zero
573 std::fill_n(laplacians.begin(),
574 n_quadrature_points,
576
577 // add up contributions of trial functions. note that here we deal with
578 // scalar finite elements and also note that the Laplacian is
579 // the trace of the Hessian.
580 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
581 {
582 const Number2 value = dof_values[shape_func];
583 // For auto-differentiable numbers, the fact that a DoF value is zero
584 // does not imply that its derivatives are zero as well. So we
585 // can't filter by value for these number types.
588 continue;
589
590 const Tensor<2, spacedim> *shape_hessian_ptr =
591 &shape_hessians[shape_func][0];
592 for (unsigned int point = 0; point < n_quadrature_points; ++point)
593 laplacians[point] += value * trace(*shape_hessian_ptr++);
594 }
595 }
596
597
598
599 template <int dim, int spacedim, typename VectorType, typename Number>
600 void
601 do_function_laplacians(
602 const ArrayView<Number> &dof_values,
603 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
605 const std::vector<unsigned int> &shape_function_to_row_table,
606 std::vector<VectorType> &laplacians,
607 const bool quadrature_points_fastest = false,
608 const unsigned int component_multiple = 1)
609 {
610 // initialize with zero
611 for (unsigned int i = 0; i < laplacians.size(); ++i)
612 std::fill_n(laplacians[i].begin(),
613 laplacians[i].size(),
614 typename VectorType::value_type());
615
616 // see if there the current cell has DoFs at all, and if not
617 // then there is nothing else to do.
618 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
619 if (dofs_per_cell == 0)
620 return;
621
622
623 const unsigned int n_quadrature_points = laplacians.size();
624 const unsigned int n_components = fe.n_components();
625
626 // Assert that we can write all components into the result vectors
627 const unsigned result_components = n_components * component_multiple;
628 (void)result_components;
629 if (quadrature_points_fastest)
630 {
631 AssertDimension(laplacians.size(), result_components);
632 for (unsigned int i = 0; i < laplacians.size(); ++i)
633 AssertDimension(laplacians[i].size(), n_quadrature_points);
634 }
635 else
636 {
637 AssertDimension(laplacians.size(), n_quadrature_points);
638 for (unsigned int i = 0; i < laplacians.size(); ++i)
639 AssertDimension(laplacians[i].size(), result_components);
640 }
641
642 // add up contributions of trial functions. now check whether the shape
643 // function is primitive or not. if it is, then set its only non-zero
644 // component, otherwise loop over components
645 for (unsigned int mc = 0; mc < component_multiple; ++mc)
646 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
647 ++shape_func)
648 {
649 const Number &value = dof_values[shape_func + mc * dofs_per_cell];
650 // For auto-differentiable numbers, the fact that a DoF value is zero
651 // does not imply that its derivatives are zero as well. So we
652 // can't filter by value for these number types.
653 if (::internal::CheckForZero<Number>::value(value) == true)
654 continue;
655
656 if (fe.is_primitive(shape_func))
657 {
658 const unsigned int comp =
659 fe.system_to_component_index(shape_func).first +
660 mc * n_components;
661 const unsigned int row =
662 shape_function_to_row_table[shape_func * n_components + comp];
663
664 const Tensor<2, spacedim> *shape_hessian_ptr =
665 &shape_hessians[row][0];
666 if (quadrature_points_fastest)
667 {
668 VectorType &laplacians_comp = laplacians[comp];
669 for (unsigned int point = 0; point < n_quadrature_points;
670 ++point)
671 laplacians_comp[point] +=
672 value * trace(*shape_hessian_ptr++);
673 }
674 else
675 for (unsigned int point = 0; point < n_quadrature_points;
676 ++point)
677 laplacians[point][comp] +=
678 value * trace(*shape_hessian_ptr++);
679 }
680 else
681 for (unsigned int c = 0; c < n_components; ++c)
682 {
683 if (fe.get_nonzero_components(shape_func)[c] == false)
684 continue;
685
686 const unsigned int row =
687 shape_function_to_row_table[shape_func * n_components + c];
688
689 const Tensor<2, spacedim> *shape_hessian_ptr =
690 &shape_hessians[row][0];
691 const unsigned int comp = c + mc * n_components;
692
693 if (quadrature_points_fastest)
694 {
695 VectorType &laplacians_comp = laplacians[comp];
696 for (unsigned int point = 0; point < n_quadrature_points;
697 ++point)
698 laplacians_comp[point] +=
699 value * trace(*shape_hessian_ptr++);
700 }
701 else
702 for (unsigned int point = 0; point < n_quadrature_points;
703 ++point)
704 laplacians[point][comp] +=
705 value * trace(*shape_hessian_ptr++);
706 }
707 }
708 }
709} // namespace internal
710
711
713template <int dim, int spacedim>
714template <typename Number>
715void
717 const ReadVector<Number> &fe_function,
718 std::vector<Number> &values) const
719{
721 ExcAccessToUninitializedField("update_values"));
725
726 // get function values of dofs on this cell
727 Vector<Number> dof_values(dofs_per_cell);
728 present_cell.get_interpolated_dof_values(fe_function, dof_values);
730 dof_values.end()),
731 this->finite_element_output.shape_values,
732 values);
733}
734
735
736
737template <int dim, int spacedim>
738template <typename Number>
739void
741 const ReadVector<Number> &fe_function,
743 std::vector<Number> &values) const
744{
746 ExcAccessToUninitializedField("update_values"));
749
750 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
751 auto view = make_array_view(dof_values.begin(), dof_values.end());
752 fe_function.extract_subvector_to(indices, view);
755 values);
756}
757
758
759
760template <int dim, int spacedim>
761template <typename Number>
762void
764 const ReadVector<Number> &fe_function,
765 std::vector<Vector<Number>> &values) const
766{
768
770 ExcAccessToUninitializedField("update_values"));
772
773 // get function values of dofs on this cell
774 Vector<Number> dof_values(dofs_per_cell);
775 present_cell.get_interpolated_dof_values(fe_function, dof_values);
777 make_array_view(dof_values.begin(), dof_values.end()),
778 this->finite_element_output.shape_values,
779 *fe,
780 this->finite_element_output.shape_function_to_row_table,
781 make_array_view(values.begin(), values.end()));
782}
783
784
785
786template <int dim, int spacedim>
787template <typename Number>
788void
790 const ReadVector<Number> &fe_function,
792 std::vector<Vector<Number>> &values) const
793{
794 // Size of indices must be a multiple of dofs_per_cell such that an integer
795 // number of function values is generated in each point.
796 Assert(indices.size() % dofs_per_cell == 0,
797 ExcNotMultiple(indices.size(), dofs_per_cell));
799 ExcAccessToUninitializedField("update_values"));
800
801 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
802 auto view = make_array_view(dof_values.begin(), dof_values.end());
803 fe_function.extract_subvector_to(indices, view);
805 view,
807 *fe,
808 this->finite_element_output.shape_function_to_row_table,
809 make_array_view(values.begin(), values.end()),
810 false,
811 indices.size() / dofs_per_cell);
812}
813
814
815
816template <int dim, int spacedim>
817template <typename Number>
818void
820 const ReadVector<Number> &fe_function,
822 ArrayView<std::vector<Number>> values,
823 const bool quadrature_points_fastest) const
824{
825 Assert(this->update_flags & update_values,
826 ExcAccessToUninitializedField("update_values"));
827
828 // Size of indices must be a multiple of dofs_per_cell such that an integer
829 // number of function values is generated in each point.
830 Assert(indices.size() % dofs_per_cell == 0,
831 ExcNotMultiple(indices.size(), dofs_per_cell));
832
833 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
834 auto view = make_array_view(dof_values.begin(), dof_values.end());
835 fe_function.extract_subvector_to(indices, view);
837 view,
839 *fe,
840 this->finite_element_output.shape_function_to_row_table,
841 make_array_view(values.begin(), values.end()),
842 quadrature_points_fastest,
843 indices.size() / dofs_per_cell);
844}
845
846
847
848template <int dim, int spacedim>
849template <typename Number>
850void
852 const ReadVector<Number> &fe_function,
853 std::vector<Tensor<1, spacedim, Number>> &gradients) const
854{
856 ExcAccessToUninitializedField("update_gradients"));
860
861 // get function values of dofs on this cell
862 Vector<Number> dof_values(dofs_per_cell);
863 present_cell.get_interpolated_dof_values(fe_function, dof_values);
865 dof_values.end()),
866 this->finite_element_output.shape_gradients,
867 gradients);
868}
869
870
871
872template <int dim, int spacedim>
873template <typename Number>
874void
876 const ReadVector<Number> &fe_function,
878 std::vector<Tensor<1, spacedim, Number>> &gradients) const
879{
880 Assert(this->update_flags & update_gradients,
881 ExcAccessToUninitializedField("update_gradients"));
884
885 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
886 auto view = make_array_view(dof_values.begin(), dof_values.end());
887 fe_function.extract_subvector_to(indices, view);
890 gradients);
891}
892
893
894
895template <int dim, int spacedim>
896template <typename Number>
897void
899 const ReadVector<Number> &fe_function,
900 std::vector<std::vector<Tensor<1, spacedim, Number>>> &gradients) const
902 Assert(this->update_flags & update_gradients,
903 ExcAccessToUninitializedField("update_gradients"));
906
907 // get function values of dofs on this cell
908 Vector<Number> dof_values(dofs_per_cell);
909 present_cell.get_interpolated_dof_values(fe_function, dof_values);
911 make_array_view(dof_values.begin(), dof_values.end()),
912 this->finite_element_output.shape_gradients,
913 *fe,
914 this->finite_element_output.shape_function_to_row_table,
915 make_array_view(gradients.begin(), gradients.end()));
917
918
919
920template <int dim, int spacedim>
921template <typename Number>
922void
924 const ReadVector<Number> &fe_function,
926 ArrayView<std::vector<Tensor<1, spacedim, Number>>> gradients,
927 const bool quadrature_points_fastest) const
928{
929 // Size of indices must be a multiple of dofs_per_cell such that an integer
930 // number of function values is generated in each point.
931 Assert(indices.size() % dofs_per_cell == 0,
932 ExcNotMultiple(indices.size(), dofs_per_cell));
933 Assert(this->update_flags & update_gradients,
934 ExcAccessToUninitializedField("update_gradients"));
935
936 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
937 auto view = make_array_view(dof_values.begin(), dof_values.end());
938 fe_function.extract_subvector_to(indices, view);
940 view,
942 *fe,
943 this->finite_element_output.shape_function_to_row_table,
944 make_array_view(gradients.begin(), gradients.end()),
945 quadrature_points_fastest,
946 indices.size() / dofs_per_cell);
947}
948
949
950
951template <int dim, int spacedim>
952template <typename Number>
953void
955 const ReadVector<Number> &fe_function,
956 std::vector<Tensor<2, spacedim, Number>> &hessians) const
957{
960 ExcAccessToUninitializedField("update_hessians"));
963
964 // get function values of dofs on this cell
965 Vector<Number> dof_values(dofs_per_cell);
966 present_cell.get_interpolated_dof_values(fe_function, dof_values);
968 dof_values.end()),
969 this->finite_element_output.shape_hessians,
970 hessians);
971}
972
973
974
975template <int dim, int spacedim>
976template <typename Number>
977void
979 const ReadVector<Number> &fe_function,
981 std::vector<Tensor<2, spacedim, Number>> &hessians) const
982{
984 ExcAccessToUninitializedField("update_hessians"));
987
988 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
989 auto view = make_array_view(dof_values.begin(), dof_values.end());
990 fe_function.extract_subvector_to(indices, view);
993 hessians);
994}
995
996
997
998template <int dim, int spacedim>
999template <typename Number>
1000void
1002 const ReadVector<Number> &fe_function,
1003 std::vector<std::vector<Tensor<2, spacedim, Number>>> &hessians,
1004 const bool quadrature_points_fastest) const
1005{
1007 ExcAccessToUninitializedField("update_hessians"));
1010
1011 // get function values of dofs on this cell
1012 Vector<Number> dof_values(dofs_per_cell);
1013 present_cell.get_interpolated_dof_values(fe_function, dof_values);
1015 make_array_view(dof_values.begin(), dof_values.end()),
1016 this->finite_element_output.shape_hessians,
1017 *fe,
1018 this->finite_element_output.shape_function_to_row_table,
1019 make_array_view(hessians.begin(), hessians.end()),
1020 quadrature_points_fastest);
1022
1023
1024
1025template <int dim, int spacedim>
1026template <typename Number>
1027void
1029 const ReadVector<Number> &fe_function,
1031 ArrayView<std::vector<Tensor<2, spacedim, Number>>> hessians,
1032 const bool quadrature_points_fastest) const
1033{
1034 Assert(this->update_flags & update_hessians,
1035 ExcAccessToUninitializedField("update_hessians"));
1036 Assert(indices.size() % dofs_per_cell == 0,
1037 ExcNotMultiple(indices.size(), dofs_per_cell));
1038
1039 boost::container::small_vector<Number, 200> dof_values(indices.size());
1040 auto view = make_array_view(dof_values.begin(), dof_values.end());
1041 fe_function.extract_subvector_to(indices, view);
1043 view,
1045 *fe,
1046 this->finite_element_output.shape_function_to_row_table,
1047 make_array_view(hessians.begin(), hessians.end()),
1048 quadrature_points_fastest,
1049 indices.size() / dofs_per_cell);
1050}
1051
1052
1053
1054template <int dim, int spacedim>
1055template <typename Number>
1056void
1058 const ReadVector<Number> &fe_function,
1059 std::vector<Number> &laplacians) const
1060{
1061 Assert(this->update_flags & update_hessians,
1062 ExcAccessToUninitializedField("update_hessians"));
1066
1067 // get function values of dofs on this cell
1068 Vector<Number> dof_values(dofs_per_cell);
1069 present_cell.get_interpolated_dof_values(fe_function, dof_values);
1071 dof_values.end()),
1072 this->finite_element_output.shape_hessians,
1073 laplacians);
1074}
1075
1076
1077
1078template <int dim, int spacedim>
1079template <typename Number>
1080void
1082 const ReadVector<Number> &fe_function,
1084 std::vector<Number> &laplacians) const
1085{
1086 Assert(this->update_flags & update_hessians,
1087 ExcAccessToUninitializedField("update_hessians"));
1090
1091 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
1092 auto view = make_array_view(dof_values.begin(), dof_values.end());
1093 fe_function.extract_subvector_to(indices, view);
1096 laplacians);
1097}
1098
1099
1100
1101template <int dim, int spacedim>
1102template <typename Number>
1103void
1105 const ReadVector<Number> &fe_function,
1106 std::vector<Vector<Number>> &laplacians) const
1107{
1109 Assert(this->update_flags & update_hessians,
1110 ExcAccessToUninitializedField("update_hessians"));
1112
1113 // get function values of dofs on this cell
1114 Vector<Number> dof_values(dofs_per_cell);
1115 present_cell.get_interpolated_dof_values(fe_function, dof_values);
1117 make_array_view(dof_values.begin(), dof_values.end()),
1118 this->finite_element_output.shape_hessians,
1119 *fe,
1120 this->finite_element_output.shape_function_to_row_table,
1121 laplacians);
1122}
1123
1124
1125
1126template <int dim, int spacedim>
1127template <typename Number>
1128void
1130 const ReadVector<Number> &fe_function,
1132 std::vector<Vector<Number>> &laplacians) const
1133{
1134 // Size of indices must be a multiple of dofs_per_cell such that an integer
1135 // number of function values is generated in each point.
1136 Assert(indices.size() % dofs_per_cell == 0,
1137 ExcNotMultiple(indices.size(), dofs_per_cell));
1138 Assert(this->update_flags & update_hessians,
1139 ExcAccessToUninitializedField("update_hessians"));
1140
1141 boost::container::small_vector<Number, 200> dof_values(indices.size());
1142 auto view = make_array_view(dof_values.begin(), dof_values.end());
1143 fe_function.extract_subvector_to(indices, view);
1145 view,
1147 *fe,
1148 this->finite_element_output.shape_function_to_row_table,
1149 laplacians,
1150 false,
1151 indices.size() / dofs_per_cell);
1152}
1153
1154
1155
1156template <int dim, int spacedim>
1157template <typename Number>
1158void
1160 const ReadVector<Number> &fe_function,
1162 std::vector<std::vector<Number>> &laplacians,
1163 const bool quadrature_points_fastest) const
1164{
1165 Assert(indices.size() % dofs_per_cell == 0,
1166 ExcNotMultiple(indices.size(), dofs_per_cell));
1167 Assert(this->update_flags & update_hessians,
1168 ExcAccessToUninitializedField("update_hessians"));
1169
1170 boost::container::small_vector<Number, 200> dof_values(indices.size());
1171 auto view = make_array_view(dof_values.begin(), dof_values.end());
1172 fe_function.extract_subvector_to(indices, view);
1174 view,
1176 *fe,
1177 this->finite_element_output.shape_function_to_row_table,
1178 laplacians,
1179 quadrature_points_fastest,
1180 indices.size() / dofs_per_cell);
1181}
1182
1183
1184
1185template <int dim, int spacedim>
1186template <typename Number>
1187void
1189 const ReadVector<Number> &fe_function,
1190 std::vector<Tensor<3, spacedim, Number>> &third_derivatives) const
1191{
1194 ExcAccessToUninitializedField("update_3rd_derivatives"));
1197
1198 // get function values of dofs on this cell
1199 Vector<Number> dof_values(dofs_per_cell);
1200 present_cell.get_interpolated_dof_values(fe_function, dof_values);
1202 make_array_view(dof_values.begin(), dof_values.end()),
1203 this->finite_element_output.shape_3rd_derivatives,
1204 third_derivatives);
1205}
1206
1207
1208
1209template <int dim, int spacedim>
1210template <typename Number>
1211void
1213 const ReadVector<Number> &fe_function,
1215 std::vector<Tensor<3, spacedim, Number>> &third_derivatives) const
1216{
1218 ExcAccessToUninitializedField("update_3rd_derivatives"));
1221
1222 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
1223 auto view = make_array_view(dof_values.begin(), dof_values.end());
1224 fe_function.extract_subvector_to(indices, view);
1226 view, this->finite_element_output.shape_3rd_derivatives, third_derivatives);
1227}
1228
1229
1230
1231template <int dim, int spacedim>
1232template <typename Number>
1233void
1235 const ReadVector<Number> &fe_function,
1236 std::vector<std::vector<Tensor<3, spacedim, Number>>> &third_derivatives,
1237 const bool quadrature_points_fastest) const
1238{
1240 ExcAccessToUninitializedField("update_3rd_derivatives"));
1243
1244 // get function values of dofs on this cell
1245 Vector<Number> dof_values(dofs_per_cell);
1246 present_cell.get_interpolated_dof_values(fe_function, dof_values);
1248 make_array_view(dof_values.begin(), dof_values.end()),
1249 this->finite_element_output.shape_3rd_derivatives,
1250 *fe,
1251 this->finite_element_output.shape_function_to_row_table,
1252 make_array_view(third_derivatives.begin(), third_derivatives.end()),
1253 quadrature_points_fastest);
1254}
1255
1256
1257
1258template <int dim, int spacedim>
1259template <typename Number>
1260void
1262 const ReadVector<Number> &fe_function,
1264 ArrayView<std::vector<Tensor<3, spacedim, Number>>> third_derivatives,
1265 const bool quadrature_points_fastest) const
1266{
1268 ExcAccessToUninitializedField("update_3rd_derivatives"));
1269 Assert(indices.size() % dofs_per_cell == 0,
1270 ExcNotMultiple(indices.size(), dofs_per_cell));
1271
1272 boost::container::small_vector<Number, 200> dof_values(indices.size());
1273 auto view = make_array_view(dof_values.begin(), dof_values.end());
1274 fe_function.extract_subvector_to(indices, view);
1276 view,
1278 *fe,
1279 this->finite_element_output.shape_function_to_row_table,
1280 make_array_view(third_derivatives.begin(), third_derivatives.end()),
1281 quadrature_points_fastest,
1282 indices.size() / dofs_per_cell);
1283}
1284
1285
1286
1287template <int dim, int spacedim>
1293
1294
1295
1296template <int dim, int spacedim>
1297const std::vector<Tensor<1, spacedim>> &
1306
1307
1308
1309template <int dim, int spacedim>
1310std::size_t
1327
1328
1329
1330template <int dim, int spacedim>
1333 const UpdateFlags update_flags) const
1334{
1335 // first find out which objects need to be recomputed on each
1336 // cell we visit. this we have to ask the finite element and mapping.
1337 // elements are first since they might require update in mapping
1338 //
1339 // there is no need to iterate since mappings will never require
1340 // the finite element to compute something for them
1342 flags |= mapping->requires_update_flags(flags);
1343
1344 return flags;
1345}
1346
1347
1348
1349template <int dim, int spacedim>
1350void
1352{
1353 // if there is no present cell, then we shouldn't be
1354 // connected via a signal to a triangulation
1356
1357 // so delete the present cell and
1358 // disconnect from the signal we have with
1359 // it
1360 tria_listener_refinement.disconnect();
1361 tria_listener_mesh_transform.disconnect();
1362 present_cell = {};
1363}
1364
1365
1366
1367template <int dim, int spacedim>
1368void
1371{
1373 {
1374 if (&cell->get_triangulation() !=
1376 .
1378 ->get_triangulation())
1379 {
1380 // the triangulations for the previous cell and the current cell
1381 // do not match. disconnect from the previous triangulation and
1382 // connect to the current one; also invalidate the previous
1383 // cell because we shouldn't be comparing cells from different
1384 // triangulations
1387 cell->get_triangulation().signals.any_change.connect(
1388 [this]() { this->invalidate_present_cell(); });
1390 cell->get_triangulation().signals.mesh_movement.connect(
1391 [this]() { this->invalidate_present_cell(); });
1392 }
1393 }
1394 else
1395 {
1396 // if this FEValues has never been set to any cell at all, then
1397 // at least subscribe to the triangulation to get notified of
1398 // changes
1400 cell->get_triangulation().signals.post_refinement.connect(
1401 [this]() { this->invalidate_present_cell(); });
1403 cell->get_triangulation().signals.mesh_movement.connect(
1404 [this]() { this->invalidate_present_cell(); });
1405 }
1406}
1407
1408
1409
1410template <int dim, int spacedim>
1411inline void
1414{
1416 {
1418 return;
1419 }
1420
1421 // case that there has not been any cell before
1422 if (this->present_cell.is_initialized() == false)
1424 else
1425 // in MappingQ, data can have been modified during the last call. Then, we
1426 // can't use that data on the new cell.
1429 else
1431 (cell->is_translation_of(
1432 static_cast<
1434 this->present_cell)) ?
1437
1438 if ((dim == spacedim - 1) && (cell_similarity == CellSimilarity::translation))
1439 {
1440 if (static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
1441 &>(this->present_cell)
1442 ->direction_flag() != cell->direction_flag())
1444 }
1445 // TODO: here, one could implement other checks for similarity, e.g. for
1446 // children of a parallelogram.
1447}
1448
1449
1450
1451template <int dim, int spacedim>
1457
1458
1459
1460template <int dim, int spacedim>
1462
1463
1464
1465template <int dim, int spacedim>
1467
1468/*-------------------------- Explicit Instantiations -------------------------*/
1469
1470
1471#include "fe_values_base.inst"
1472
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
std::size_t size() const
Definition array_view.h:684
typename LevelSelector::cell_iterator level_cell_iterator
types::global_dof_index n_dofs_for_dof_handler() const
void get_interpolated_dof_values(const ReadVector< Number > &in, Vector< Number > &out) const
CellSimilarity::Similarity cell_similarity
::internal::FEValuesViews::Cache< dim, spacedim > fe_values_views_cache
Triangulation< dim, spacedim >::cell_iterator get_cell() const
boost::signals2::connection tria_listener_mesh_transform
void get_function_values(const ReadVector< Number > &fe_function, std::vector< Number > &values) const
FEValuesBase(const unsigned int n_q_points, const unsigned int dofs_per_cell, const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe)
internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > mapping_output
virtual ~FEValuesBase() override
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
const unsigned int dofs_per_cell
void check_cell_similarity(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
CellIteratorWrapper present_cell
UpdateFlags update_flags
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
void get_function_hessians(const ReadVector< Number > &fe_function, std::vector< Tensor< 2, spacedim, Number > > &hessians) const
void always_allow_check_for_cell_similarity(const bool allow)
void get_function_laplacians(const ReadVector< Number > &fe_function, std::vector< Number > &laplacians) const
const unsigned int n_quadrature_points
CellSimilarity::Similarity get_cell_similarity() const
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
std::size_t memory_consumption() const
void get_function_third_derivatives(const ReadVector< Number > &fe_function, std::vector< Tensor< 3, spacedim, Number > > &third_derivatives) const
UpdateFlags compute_update_flags(const UpdateFlags update_flags) const
std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > mapping_data
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
boost::signals2::connection tria_listener_refinement
void invalidate_present_cell()
bool check_for_cell_similarity_allowed
::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > finite_element_output
std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > fe_data
void maybe_invalidate_previous_present_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
const unsigned int max_n_quadrature_points
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const =0
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
unsigned int n_nonzero_components(const unsigned int i) const
Abstract base class for mapping classes.
Definition mapping.h:318
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const =0
virtual size_type size() const =0
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, ArrayView< Number > &elements) const =0
iterator end()
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcAccessToUninitializedField(std::string arg1)
static ::ExceptionBase & ExcNotMultiple(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotReinited()
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition tria.h:1550
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_gradients
Shape function gradients.
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void do_function_laplacians(const ArrayView< Number2 > &dof_values, const ::Table< 2, Tensor< 2, spacedim > > &shape_hessians, std::vector< Number > &laplacians)
void do_function_derivatives(const ArrayView< Number > &dof_values, const ::Table< 2, Tensor< order, spacedim > > &shape_derivatives, std::vector< Tensor< order, spacedim, Number > > &derivatives)
std::vector< unsigned int > make_shape_function_to_row_table(const FiniteElement< dim, spacedim > &fe)
Definition fe_values.cc:49
void do_function_values(const ArrayView< Number2 > &dof_values, const ::Table< 2, double > &shape_values, std::vector< Number > &values)
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252
STL namespace.
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:702
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)