Reference documentation for deal.II version 9.5.0
\(\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.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
23
25
27
28#include <deal.II/fe/fe.h>
30#include <deal.II/fe/mapping.h>
31
34
45#include <deal.II/lac/vector.h>
47
48#include <boost/container/small_vector.hpp>
49
50#include <iomanip>
51#include <memory>
52#include <type_traits>
53
55
56
57namespace internal
58{
59 template <class VectorType>
60 typename VectorType::value_type inline get_vector_element(
61 const VectorType & vector,
62 const types::global_dof_index cell_number)
63 {
64 return internal::ElementAccess<VectorType>::get(vector, cell_number);
65 }
66
67
68
70 const IndexSet & is,
71 const types::global_dof_index cell_number)
72 {
73 return (is.is_element(cell_number) ? 1 : 0);
74 }
75
76
77
78 template <int dim, int spacedim>
79 inline std::vector<unsigned int>
81 {
82 std::vector<unsigned int> shape_function_to_row_table(
84 unsigned int row = 0;
85 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
86 {
87 // loop over all components that are nonzero for this particular
88 // shape function. if a component is zero then we leave the
89 // value in the table unchanged (at the invalid value)
90 // otherwise it is mapped to the next free entry
91 unsigned int nth_nonzero_component = 0;
92 for (unsigned int c = 0; c < fe.n_components(); ++c)
93 if (fe.get_nonzero_components(i)[c] == true)
94 {
95 shape_function_to_row_table[i * fe.n_components() + c] =
96 row + nth_nonzero_component;
97 ++nth_nonzero_component;
98 }
99 row += fe.n_nonzero_components(i);
100 }
101
102 return shape_function_to_row_table;
103 }
104
105 namespace
106 {
107 // Check to see if a DoF value is zero, implying that subsequent operations
108 // with the value have no effect.
109 template <typename Number, typename T = void>
110 struct CheckForZero
111 {
112 static bool
113 value(const Number &value)
114 {
116 }
117 };
118
119 // For auto-differentiable numbers, the fact that a DoF value is zero
120 // does not imply that its derivatives are zero as well. So we
121 // can't filter by value for these number types.
122 // Note that we also want to avoid actually checking the value itself,
123 // since some AD numbers are not contextually convertible to booleans.
124 template <typename Number>
125 struct CheckForZero<
126 Number,
127 std::enable_if_t<Differentiation::AD::is_ad_number<Number>::value>>
128 {
129 static bool
130 value(const Number & /*value*/)
131 {
132 return false;
133 }
134 };
135 } // namespace
136} // namespace internal
137
138
139
140namespace FEValuesViews
141{
142 template <int dim, int spacedim>
144 const unsigned int component)
145 : fe_values(&fe_values)
146 , component(component)
147 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
148 {
149 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
151
152 // TODO: we'd like to use the fields with the same name as these
153 // variables from FEValuesBase, but they aren't initialized yet
154 // at the time we get here, so re-create it all
155 const std::vector<unsigned int> shape_function_to_row_table =
157
158 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
159 {
160 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
161
162 if (is_primitive == true)
163 shape_function_data[i].is_nonzero_shape_function_component =
164 (component == fe.system_to_component_index(i).first);
165 else
166 shape_function_data[i].is_nonzero_shape_function_component =
167 (fe.get_nonzero_components(i)[component] == true);
168
169 if (shape_function_data[i].is_nonzero_shape_function_component == true)
170 shape_function_data[i].row_index =
171 shape_function_to_row_table[i * fe.n_components() + component];
172 else
174 }
175 }
176
177
178
179 template <int dim, int spacedim>
181 : fe_values(nullptr)
182 , component(numbers::invalid_unsigned_int)
183 {}
184
185
186
187 template <int dim, int spacedim>
189 const unsigned int first_vector_component)
190 : fe_values(&fe_values)
191 , first_vector_component(first_vector_component)
192 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
193 {
194 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
196
197 // TODO: we'd like to use the fields with the same name as these
198 // variables from FEValuesBase, but they aren't initialized yet
199 // at the time we get here, so re-create it all
200 const std::vector<unsigned int> shape_function_to_row_table =
202
203 for (unsigned int d = 0; d < spacedim; ++d)
204 {
205 const unsigned int component = first_vector_component + d;
206
207 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
208 {
209 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
210
211 if (is_primitive == true)
212 shape_function_data[i].is_nonzero_shape_function_component[d] =
213 (component == fe.system_to_component_index(i).first);
214 else
215 shape_function_data[i].is_nonzero_shape_function_component[d] =
216 (fe.get_nonzero_components(i)[component] == true);
217
218 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
219 true)
220 shape_function_data[i].row_index[d] =
221 shape_function_to_row_table[i * fe.n_components() + component];
222 else
223 shape_function_data[i].row_index[d] =
225 }
226 }
227
228 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
229 {
230 unsigned int n_nonzero_components = 0;
231 for (unsigned int d = 0; d < spacedim; ++d)
232 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
233 true)
234 ++n_nonzero_components;
235
236 if (n_nonzero_components == 0)
237 shape_function_data[i].single_nonzero_component = -2;
238 else if (n_nonzero_components > 1)
239 shape_function_data[i].single_nonzero_component = -1;
240 else
241 {
242 for (unsigned int d = 0; d < spacedim; ++d)
244 .is_nonzero_shape_function_component[d] == true)
245 {
246 shape_function_data[i].single_nonzero_component =
247 shape_function_data[i].row_index[d];
248 shape_function_data[i].single_nonzero_component_index = d;
249 break;
250 }
251 }
252 }
253 }
254
255
256
257 template <int dim, int spacedim>
259 : fe_values(nullptr)
260 , first_vector_component(numbers::invalid_unsigned_int)
261 {}
262
263
264
265 template <int dim, int spacedim>
267 const FEValuesBase<dim, spacedim> &fe_values,
268 const unsigned int first_tensor_component)
269 : fe_values(&fe_values)
270 , first_tensor_component(first_tensor_component)
271 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
272 {
273 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
274 Assert(first_tensor_component + (dim * dim + dim) / 2 - 1 <
275 fe.n_components(),
277 first_tensor_component +
279 0,
280 fe.n_components()));
281 // TODO: we'd like to use the fields with the same name as these
282 // variables from FEValuesBase, but they aren't initialized yet
283 // at the time we get here, so re-create it all
284 const std::vector<unsigned int> shape_function_to_row_table =
286
287 for (unsigned int d = 0;
288 d < ::SymmetricTensor<2, dim>::n_independent_components;
289 ++d)
290 {
291 const unsigned int component = first_tensor_component + d;
292
293 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
294 {
295 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
296
297 if (is_primitive == true)
298 shape_function_data[i].is_nonzero_shape_function_component[d] =
299 (component == fe.system_to_component_index(i).first);
300 else
301 shape_function_data[i].is_nonzero_shape_function_component[d] =
302 (fe.get_nonzero_components(i)[component] == true);
303
304 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
305 true)
306 shape_function_data[i].row_index[d] =
307 shape_function_to_row_table[i * fe.n_components() + component];
308 else
309 shape_function_data[i].row_index[d] =
311 }
312 }
313
314 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
315 {
316 unsigned int n_nonzero_components = 0;
317 for (unsigned int d = 0;
318 d < ::SymmetricTensor<2, dim>::n_independent_components;
319 ++d)
320 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
321 true)
322 ++n_nonzero_components;
323
324 if (n_nonzero_components == 0)
325 shape_function_data[i].single_nonzero_component = -2;
326 else if (n_nonzero_components > 1)
327 shape_function_data[i].single_nonzero_component = -1;
328 else
329 {
330 for (unsigned int d = 0;
331 d < ::SymmetricTensor<2, dim>::n_independent_components;
332 ++d)
333 if (shape_function_data[i]
334 .is_nonzero_shape_function_component[d] == true)
335 {
336 shape_function_data[i].single_nonzero_component =
337 shape_function_data[i].row_index[d];
338 shape_function_data[i].single_nonzero_component_index = d;
339 break;
340 }
341 }
342 }
343 }
344
345
346
347 template <int dim, int spacedim>
349 : fe_values(nullptr)
350 , first_tensor_component(numbers::invalid_unsigned_int)
351 {}
352
353
354
355 template <int dim, int spacedim>
357 const unsigned int first_tensor_component)
358 : fe_values(&fe_values)
359 , first_tensor_component(first_tensor_component)
360 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
361 {
362 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
363 AssertIndexRange(first_tensor_component + dim * dim - 1, fe.n_components());
364 // TODO: we'd like to use the fields with the same name as these
365 // variables from FEValuesBase, but they aren't initialized yet
366 // at the time we get here, so re-create it all
367 const std::vector<unsigned int> shape_function_to_row_table =
369
370 for (unsigned int d = 0; d < dim * dim; ++d)
371 {
372 const unsigned int component = first_tensor_component + d;
373
374 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
375 {
376 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
377
378 if (is_primitive == true)
379 shape_function_data[i].is_nonzero_shape_function_component[d] =
380 (component == fe.system_to_component_index(i).first);
381 else
382 shape_function_data[i].is_nonzero_shape_function_component[d] =
383 (fe.get_nonzero_components(i)[component] == true);
384
385 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
386 true)
387 shape_function_data[i].row_index[d] =
388 shape_function_to_row_table[i * fe.n_components() + component];
389 else
390 shape_function_data[i].row_index[d] =
392 }
393 }
394
395 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
396 {
397 unsigned int n_nonzero_components = 0;
398 for (unsigned int d = 0; d < dim * dim; ++d)
399 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
400 true)
401 ++n_nonzero_components;
402
403 if (n_nonzero_components == 0)
404 shape_function_data[i].single_nonzero_component = -2;
405 else if (n_nonzero_components > 1)
406 shape_function_data[i].single_nonzero_component = -1;
407 else
408 {
409 for (unsigned int d = 0; d < dim * dim; ++d)
410 if (shape_function_data[i]
411 .is_nonzero_shape_function_component[d] == true)
413 shape_function_data[i].single_nonzero_component =
414 shape_function_data[i].row_index[d];
415 shape_function_data[i].single_nonzero_component_index = d;
416 break;
417 }
418 }
419 }
420 }
422
423
424 template <int dim, int spacedim>
426 : fe_values(nullptr)
427 , first_tensor_component(numbers::invalid_unsigned_int)
428 {}
429
430
431
432 namespace internal
433 {
434 // Given values of degrees of freedom, evaluate the
435 // values/gradients/... at quadrature points
436
437 // ------------------------- scalar functions --------------------------
438 template <int dim, int spacedim, typename Number>
439 void
441 const ArrayView<Number> &dof_values,
442 const Table<2, double> & shape_values,
443 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
444 &shape_function_data,
445 std::vector<typename ProductType<Number, double>::type> &values)
446 {
447 const unsigned int dofs_per_cell = dof_values.size();
448 const unsigned int n_quadrature_points = values.size();
449
450 std::fill(values.begin(),
451 values.end(),
453
454 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
455 ++shape_function)
456 if (shape_function_data[shape_function]
457 .is_nonzero_shape_function_component)
458 {
459 const Number &value = dof_values[shape_function];
460 // For auto-differentiable numbers, the fact that a DoF value is
461 // zero does not imply that its derivatives are zero as well. So we
462 // can't filter by value for these number types.
463 if (::internal::CheckForZero<Number>::value(value) == true)
464 continue;
466 const double *shape_value_ptr =
467 &shape_values(shape_function_data[shape_function].row_index, 0);
468 for (unsigned int q_point = 0; q_point < n_quadrature_points;
469 ++q_point)
470 values[q_point] += value * (*shape_value_ptr++);
471 }
472 }
473
474
475
476 // same code for gradient and Hessian, template argument 'order' to give
477 // the order of the derivative (= rank of gradient/Hessian tensor)
478 template <int order, int dim, int spacedim, typename Number>
479 void
481 const ArrayView<Number> & dof_values,
482 const Table<2, ::Tensor<order, spacedim>> &shape_derivatives,
483 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
484 &shape_function_data,
485 std::vector<
486 typename ProductType<Number, ::Tensor<order, spacedim>>::type>
487 &derivatives)
488 {
489 const unsigned int dofs_per_cell = dof_values.size();
490 const unsigned int n_quadrature_points = derivatives.size();
491
492 std::fill(
493 derivatives.begin(),
494 derivatives.end(),
495 typename ProductType<Number, ::Tensor<order, spacedim>>::type());
496
497 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
498 ++shape_function)
499 if (shape_function_data[shape_function]
500 .is_nonzero_shape_function_component)
501 {
502 const Number &value = dof_values[shape_function];
503 // For auto-differentiable numbers, the fact that a DoF value is
504 // zero does not imply that its derivatives are zero as well. So we
505 // can't filter by value for these number types.
506 if (::internal::CheckForZero<Number>::value(value) == true)
507 continue;
508
509 const ::Tensor<order, spacedim> *shape_derivative_ptr =
510 &shape_derivatives[shape_function_data[shape_function].row_index]
511 [0];
512 for (unsigned int q_point = 0; q_point < n_quadrature_points;
513 ++q_point)
514 derivatives[q_point] += value * (*shape_derivative_ptr++);
515 }
516 }
517
518
519
520 template <int dim, int spacedim, typename Number>
521 void
523 const ArrayView<Number> & dof_values,
524 const Table<2, ::Tensor<2, spacedim>> &shape_hessians,
525 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
526 &shape_function_data,
527 std::vector<typename Scalar<dim, spacedim>::
528 template solution_laplacian_type<Number>> &laplacians)
529 {
530 const unsigned int dofs_per_cell = dof_values.size();
531 const unsigned int n_quadrature_points = laplacians.size();
532
533 std::fill(
534 laplacians.begin(),
535 laplacians.end(),
536 typename Scalar<dim,
537 spacedim>::template solution_laplacian_type<Number>());
538
539 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
540 ++shape_function)
541 if (shape_function_data[shape_function]
542 .is_nonzero_shape_function_component)
543 {
544 const Number &value = dof_values[shape_function];
545 // For auto-differentiable numbers, the fact that a DoF value is
546 // zero does not imply that its derivatives are zero as well. So we
547 // can't filter by value for these number types.
548 if (::internal::CheckForZero<Number>::value(value) == true)
549 continue;
550
551 const ::Tensor<2, spacedim> *shape_hessian_ptr =
552 &shape_hessians[shape_function_data[shape_function].row_index][0];
553 for (unsigned int q_point = 0; q_point < n_quadrature_points;
554 ++q_point)
555 laplacians[q_point] += value * trace(*shape_hessian_ptr++);
556 }
557 }
558
559
560
561 // ----------------------------- vector part ---------------------------
562
563 template <int dim, int spacedim, typename Number>
564 void
566 const ArrayView<Number> &dof_values,
567 const Table<2, double> & shape_values,
568 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
569 &shape_function_data,
570 std::vector<
571 typename ProductType<Number, ::Tensor<1, spacedim>>::type>
572 &values)
573 {
574 const unsigned int dofs_per_cell = dof_values.size();
575 const unsigned int n_quadrature_points = values.size();
576
577 std::fill(
578 values.begin(),
579 values.end(),
580 typename ProductType<Number, ::Tensor<1, spacedim>>::type());
581
582 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
583 ++shape_function)
584 {
585 const int snc =
586 shape_function_data[shape_function].single_nonzero_component;
587
588 if (snc == -2)
589 // shape function is zero for the selected components
590 continue;
591
592 const Number &value = dof_values[shape_function];
593 // For auto-differentiable numbers, the fact that a DoF value is zero
594 // does not imply that its derivatives are zero as well. So we
595 // can't filter by value for these number types.
596 if (::internal::CheckForZero<Number>::value(value) == true)
597 continue;
598
599 if (snc != -1)
600 {
601 const unsigned int comp = shape_function_data[shape_function]
602 .single_nonzero_component_index;
603 const double *shape_value_ptr = &shape_values(snc, 0);
604 for (unsigned int q_point = 0; q_point < n_quadrature_points;
605 ++q_point)
606 values[q_point][comp] += value * (*shape_value_ptr++);
607 }
608 else
609 for (unsigned int d = 0; d < spacedim; ++d)
610 if (shape_function_data[shape_function]
611 .is_nonzero_shape_function_component[d])
612 {
613 const double *shape_value_ptr = &shape_values(
614 shape_function_data[shape_function].row_index[d], 0);
615 for (unsigned int q_point = 0; q_point < n_quadrature_points;
616 ++q_point)
617 values[q_point][d] += value * (*shape_value_ptr++);
619 }
620 }
621
622
623
624 template <int order, int dim, int spacedim, typename Number>
625 void
627 const ArrayView<Number> & dof_values,
628 const Table<2, ::Tensor<order, spacedim>> &shape_derivatives,
629 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
630 &shape_function_data,
631 std::vector<
632 typename ProductType<Number, ::Tensor<order + 1, spacedim>>::type>
633 &derivatives)
634 {
635 const unsigned int dofs_per_cell = dof_values.size();
636 const unsigned int n_quadrature_points = derivatives.size();
637
638 std::fill(
639 derivatives.begin(),
640 derivatives.end(),
641 typename ProductType<Number,
643
644 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
645 ++shape_function)
646 {
647 const int snc =
648 shape_function_data[shape_function].single_nonzero_component;
649
650 if (snc == -2)
651 // shape function is zero for the selected components
652 continue;
653
654 const Number &value = dof_values[shape_function];
655 // For auto-differentiable numbers, the fact that a DoF value is zero
656 // does not imply that its derivatives are zero as well. So we
657 // can't filter by value for these number types.
658 if (::internal::CheckForZero<Number>::value(value) == true)
659 continue;
660
661 if (snc != -1)
662 {
663 const unsigned int comp = shape_function_data[shape_function]
664 .single_nonzero_component_index;
665 const ::Tensor<order, spacedim> *shape_derivative_ptr =
666 &shape_derivatives[snc][0];
667 for (unsigned int q_point = 0; q_point < n_quadrature_points;
668 ++q_point)
669 derivatives[q_point][comp] += value * (*shape_derivative_ptr++);
670 }
671 else
672 for (unsigned int d = 0; d < spacedim; ++d)
673 if (shape_function_data[shape_function]
674 .is_nonzero_shape_function_component[d])
675 {
676 const ::Tensor<order, spacedim> *shape_derivative_ptr =
677 &shape_derivatives[shape_function_data[shape_function]
678 .row_index[d]][0];
679 for (unsigned int q_point = 0; q_point < n_quadrature_points;
680 ++q_point)
681 derivatives[q_point][d] +=
682 value * (*shape_derivative_ptr++);
683 }
684 }
685 }
686
687
688
689 template <int dim, int spacedim, typename Number>
690 void
692 const ArrayView<Number> & dof_values,
693 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
694 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
695 &shape_function_data,
696 std::vector<
697 typename ProductType<Number,
699 &symmetric_gradients)
700 {
701 const unsigned int dofs_per_cell = dof_values.size();
702 const unsigned int n_quadrature_points = symmetric_gradients.size();
703
704 std::fill(
705 symmetric_gradients.begin(),
706 symmetric_gradients.end(),
707 typename ProductType<Number,
709
710 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
711 ++shape_function)
712 {
713 const int snc =
714 shape_function_data[shape_function].single_nonzero_component;
715
716 if (snc == -2)
717 // shape function is zero for the selected components
718 continue;
719
720 const Number &value = dof_values[shape_function];
721 // For auto-differentiable numbers, the fact that a DoF value is zero
722 // does not imply that its derivatives are zero as well. So we
723 // can't filter by value for these number types.
724 if (::internal::CheckForZero<Number>::value(value) == true)
725 continue;
726
727 if (snc != -1)
728 {
729 const unsigned int comp = shape_function_data[shape_function]
730 .single_nonzero_component_index;
731 const ::Tensor<1, spacedim> *shape_gradient_ptr =
732 &shape_gradients[snc][0];
733 for (unsigned int q_point = 0; q_point < n_quadrature_points;
734 ++q_point)
735 symmetric_gradients[q_point] +=
737 symmetrize_single_row(comp, *shape_gradient_ptr++));
738 }
739 else
740 for (unsigned int q_point = 0; q_point < n_quadrature_points;
741 ++q_point)
742 {
744 grad;
745 for (unsigned int d = 0; d < spacedim; ++d)
746 if (shape_function_data[shape_function]
747 .is_nonzero_shape_function_component[d])
748 grad[d] =
749 value *
750 shape_gradients[shape_function_data[shape_function]
751 .row_index[d]][q_point];
752 symmetric_gradients[q_point] += symmetrize(grad);
753 }
754 }
755 }
756
757
758
759 template <int dim, int spacedim, typename Number>
760 void
762 const ArrayView<Number> & dof_values,
763 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
764 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
765 &shape_function_data,
766 std::vector<typename Vector<dim, spacedim>::
767 template solution_divergence_type<Number>> &divergences)
768 {
769 const unsigned int dofs_per_cell = dof_values.size();
770 const unsigned int n_quadrature_points = divergences.size();
771
772 std::fill(
773 divergences.begin(),
774 divergences.end(),
775 typename Vector<dim,
776 spacedim>::template solution_divergence_type<Number>());
777
778 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
779 ++shape_function)
780 {
781 const int snc =
782 shape_function_data[shape_function].single_nonzero_component;
783
784 if (snc == -2)
785 // shape function is zero for the selected components
786 continue;
787
788 const Number &value = dof_values[shape_function];
789 // For auto-differentiable numbers, the fact that a DoF value is zero
790 // does not imply that its derivatives are zero as well. So we
791 // can't filter by value for these number types.
792 if (::internal::CheckForZero<Number>::value(value) == true)
793 continue;
794
795 if (snc != -1)
796 {
797 const unsigned int comp = shape_function_data[shape_function]
798 .single_nonzero_component_index;
799 const ::Tensor<1, spacedim> *shape_gradient_ptr =
800 &shape_gradients[snc][0];
801 for (unsigned int q_point = 0; q_point < n_quadrature_points;
802 ++q_point)
803 divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
804 }
805 else
806 for (unsigned int d = 0; d < spacedim; ++d)
807 if (shape_function_data[shape_function]
808 .is_nonzero_shape_function_component[d])
809 {
810 const ::Tensor<1, spacedim> *shape_gradient_ptr =
811 &shape_gradients[shape_function_data[shape_function]
812 .row_index[d]][0];
813 for (unsigned int q_point = 0; q_point < n_quadrature_points;
814 ++q_point)
815 divergences[q_point] += value * (*shape_gradient_ptr++)[d];
816 }
817 }
818 }
819
820
821
822 template <int dim, int spacedim, typename Number>
823 void
825 const ArrayView<Number> & dof_values,
826 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
827 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
828 &shape_function_data,
829 std::vector<typename ProductType<
830 Number,
831 typename ::internal::CurlType<spacedim>::type>::type> &curls)
832 {
833 const unsigned int dofs_per_cell = dof_values.size();
834 const unsigned int n_quadrature_points = curls.size();
835
836 std::fill(curls.begin(),
837 curls.end(),
838 typename ProductType<
839 Number,
840 typename ::internal::CurlType<spacedim>::type>::type());
841
842 switch (spacedim)
843 {
844 case 1:
845 {
846 Assert(false,
848 "Computing the curl in 1d is not a useful operation"));
849 break;
850 }
851
852 case 2:
853 {
854 for (unsigned int shape_function = 0;
855 shape_function < dofs_per_cell;
856 ++shape_function)
857 {
858 const int snc = shape_function_data[shape_function]
859 .single_nonzero_component;
860
861 if (snc == -2)
862 // shape function is zero for the selected components
863 continue;
864
865 const Number &value = dof_values[shape_function];
866 // For auto-differentiable numbers, the fact that a DoF value
867 // is zero does not imply that its derivatives are zero as
868 // well. So we can't filter by value for these number types.
869 if (::internal::CheckForZero<Number>::value(value) ==
870 true)
871 continue;
872
873 if (snc != -1)
874 {
875 const ::Tensor<1, spacedim> *shape_gradient_ptr =
876 &shape_gradients[snc][0];
877
878 Assert(shape_function_data[shape_function]
879 .single_nonzero_component >= 0,
881 // we're in 2d, so the formula for the curl is simple:
882 if (shape_function_data[shape_function]
883 .single_nonzero_component_index == 0)
884 for (unsigned int q_point = 0;
885 q_point < n_quadrature_points;
886 ++q_point)
887 curls[q_point][0] -=
888 value * (*shape_gradient_ptr++)[1];
889 else
890 for (unsigned int q_point = 0;
891 q_point < n_quadrature_points;
892 ++q_point)
893 curls[q_point][0] +=
894 value * (*shape_gradient_ptr++)[0];
895 }
896 else
897 // we have multiple non-zero components in the shape
898 // functions. not all of them must necessarily be within the
899 // 2-component window this FEValuesViews::Vector object
900 // considers, however.
901 {
902 if (shape_function_data[shape_function]
903 .is_nonzero_shape_function_component[0])
904 {
905 const ::Tensor<1,
906 spacedim> *shape_gradient_ptr =
907 &shape_gradients[shape_function_data[shape_function]
908 .row_index[0]][0];
909
910 for (unsigned int q_point = 0;
911 q_point < n_quadrature_points;
912 ++q_point)
913 curls[q_point][0] -=
914 value * (*shape_gradient_ptr++)[1];
915 }
916
917 if (shape_function_data[shape_function]
918 .is_nonzero_shape_function_component[1])
919 {
920 const ::Tensor<1,
921 spacedim> *shape_gradient_ptr =
922 &shape_gradients[shape_function_data[shape_function]
923 .row_index[1]][0];
924
925 for (unsigned int q_point = 0;
926 q_point < n_quadrature_points;
927 ++q_point)
928 curls[q_point][0] +=
929 value * (*shape_gradient_ptr++)[0];
931 }
932 }
933 break;
934 }
935
936 case 3:
937 {
938 for (unsigned int shape_function = 0;
939 shape_function < dofs_per_cell;
940 ++shape_function)
941 {
942 const int snc = shape_function_data[shape_function]
943 .single_nonzero_component;
944
945 if (snc == -2)
946 // shape function is zero for the selected components
947 continue;
948
949 const Number &value = dof_values[shape_function];
950 // For auto-differentiable numbers, the fact that a DoF value
951 // is zero does not imply that its derivatives are zero as
952 // well. So we can't filter by value for these number types.
953 if (::internal::CheckForZero<Number>::value(value) ==
954 true)
955 continue;
956
957 if (snc != -1)
958 {
959 const ::Tensor<1, spacedim> *shape_gradient_ptr =
960 &shape_gradients[snc][0];
961
962 switch (shape_function_data[shape_function]
963 .single_nonzero_component_index)
964 {
965 case 0:
966 {
967 for (unsigned int q_point = 0;
968 q_point < n_quadrature_points;
969 ++q_point)
970 {
971 curls[q_point][1] +=
972 value * (*shape_gradient_ptr)[2];
973 curls[q_point][2] -=
974 value * (*shape_gradient_ptr++)[1];
975 }
976
977 break;
978 }
979
980 case 1:
981 {
982 for (unsigned int q_point = 0;
983 q_point < n_quadrature_points;
984 ++q_point)
985 {
986 curls[q_point][0] -=
987 value * (*shape_gradient_ptr)[2];
988 curls[q_point][2] +=
989 value * (*shape_gradient_ptr++)[0];
990 }
991
992 break;
993 }
994
995 case 2:
996 {
997 for (unsigned int q_point = 0;
998 q_point < n_quadrature_points;
999 ++q_point)
1000 {
1001 curls[q_point][0] +=
1002 value * (*shape_gradient_ptr)[1];
1003 curls[q_point][1] -=
1004 value * (*shape_gradient_ptr++)[0];
1005 }
1006 break;
1007 }
1008
1009 default:
1010 Assert(false, ExcInternalError());
1011 }
1012 }
1013
1014 else
1015 // we have multiple non-zero components in the shape
1016 // functions. not all of them must necessarily be within the
1017 // 3-component window this FEValuesViews::Vector object
1018 // considers, however.
1019 {
1020 if (shape_function_data[shape_function]
1021 .is_nonzero_shape_function_component[0])
1022 {
1023 const ::Tensor<1,
1024 spacedim> *shape_gradient_ptr =
1025 &shape_gradients[shape_function_data[shape_function]
1026 .row_index[0]][0];
1027
1028 for (unsigned int q_point = 0;
1029 q_point < n_quadrature_points;
1030 ++q_point)
1031 {
1032 curls[q_point][1] +=
1033 value * (*shape_gradient_ptr)[2];
1034 curls[q_point][2] -=
1035 value * (*shape_gradient_ptr++)[1];
1036 }
1037 }
1038
1039 if (shape_function_data[shape_function]
1040 .is_nonzero_shape_function_component[1])
1041 {
1042 const ::Tensor<1,
1043 spacedim> *shape_gradient_ptr =
1044 &shape_gradients[shape_function_data[shape_function]
1045 .row_index[1]][0];
1046
1047 for (unsigned int q_point = 0;
1048 q_point < n_quadrature_points;
1049 ++q_point)
1050 {
1051 curls[q_point][0] -=
1052 value * (*shape_gradient_ptr)[2];
1053 curls[q_point][2] +=
1054 value * (*shape_gradient_ptr++)[0];
1055 }
1056 }
1057
1058 if (shape_function_data[shape_function]
1059 .is_nonzero_shape_function_component[2])
1060 {
1061 const ::Tensor<1,
1062 spacedim> *shape_gradient_ptr =
1063 &shape_gradients[shape_function_data[shape_function]
1064 .row_index[2]][0];
1065
1066 for (unsigned int q_point = 0;
1067 q_point < n_quadrature_points;
1068 ++q_point)
1069 {
1070 curls[q_point][0] +=
1071 value * (*shape_gradient_ptr)[1];
1072 curls[q_point][1] -=
1073 value * (*shape_gradient_ptr++)[0];
1074 }
1075 }
1076 }
1077 }
1078 }
1079 }
1080 }
1081
1082
1083
1084 template <int dim, int spacedim, typename Number>
1085 void
1087 const ArrayView<Number> & dof_values,
1088 const Table<2, ::Tensor<2, spacedim>> &shape_hessians,
1089 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
1090 &shape_function_data,
1091 std::vector<typename Vector<dim, spacedim>::
1092 template solution_laplacian_type<Number>> &laplacians)
1093 {
1094 const unsigned int dofs_per_cell = dof_values.size();
1095 const unsigned int n_quadrature_points = laplacians.size();
1096
1097 std::fill(
1098 laplacians.begin(),
1099 laplacians.end(),
1100 typename Vector<dim,
1101 spacedim>::template solution_laplacian_type<Number>());
1102
1103 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1104 ++shape_function)
1105 {
1106 const int snc =
1107 shape_function_data[shape_function].single_nonzero_component;
1108
1109 if (snc == -2)
1110 // shape function is zero for the selected components
1111 continue;
1112
1113 const Number &value = dof_values[shape_function];
1114 // For auto-differentiable numbers, the fact that a DoF value is zero
1115 // does not imply that its derivatives are zero as well. So we
1116 // can't filter by value for these number types.
1117 if (::internal::CheckForZero<Number>::value(value) == true)
1118 continue;
1119
1120 if (snc != -1)
1121 {
1122 const unsigned int comp = shape_function_data[shape_function]
1123 .single_nonzero_component_index;
1124 const ::Tensor<2, spacedim> *shape_hessian_ptr =
1125 &shape_hessians[snc][0];
1126 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1127 ++q_point)
1128 laplacians[q_point][comp] +=
1129 value * trace(*shape_hessian_ptr++);
1130 }
1131 else
1132 for (unsigned int d = 0; d < spacedim; ++d)
1133 if (shape_function_data[shape_function]
1134 .is_nonzero_shape_function_component[d])
1135 {
1136 const ::Tensor<2, spacedim> *shape_hessian_ptr =
1137 &shape_hessians[shape_function_data[shape_function]
1138 .row_index[d]][0];
1139 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1140 ++q_point)
1141 laplacians[q_point][d] +=
1142 value * trace(*shape_hessian_ptr++);
1143 }
1144 }
1145 }
1146
1147
1148
1149 // ---------------------- symmetric tensor part ------------------------
1150
1151 template <int dim, int spacedim, typename Number>
1152 void
1154 const ArrayView<Number> & dof_values,
1155 const ::Table<2, double> &shape_values,
1156 const std::vector<
1158 &shape_function_data,
1159 std::vector<
1160 typename ProductType<Number,
1162 &values)
1163 {
1164 const unsigned int dofs_per_cell = dof_values.size();
1165 const unsigned int n_quadrature_points = values.size();
1166
1167 std::fill(
1168 values.begin(),
1169 values.end(),
1170 typename ProductType<Number,
1172
1173 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1174 ++shape_function)
1175 {
1176 const int snc =
1177 shape_function_data[shape_function].single_nonzero_component;
1179 if (snc == -2)
1180 // shape function is zero for the selected components
1181 continue;
1182
1183 const Number &value = dof_values[shape_function];
1184 // For auto-differentiable numbers, the fact that a DoF value is zero
1185 // does not imply that its derivatives are zero as well. So we
1186 // can't filter by value for these number types.
1187 if (::internal::CheckForZero<Number>::value(value) == true)
1188 continue;
1189
1190 if (snc != -1)
1192 const TableIndices<2> comp = ::
1194 shape_function_data[shape_function]
1195 .single_nonzero_component_index);
1196 const double *shape_value_ptr = &shape_values(snc, 0);
1197 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1198 ++q_point)
1199 values[q_point][comp] += value * (*shape_value_ptr++);
1200 }
1201 else
1202 for (unsigned int d = 0;
1203 d <
1205 ++d)
1206 if (shape_function_data[shape_function]
1207 .is_nonzero_shape_function_component[d])
1208 {
1209 const TableIndices<2> comp =
1212 const double *shape_value_ptr = &shape_values(
1213 shape_function_data[shape_function].row_index[d], 0);
1214 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1215 ++q_point)
1216 values[q_point][comp] += value * (*shape_value_ptr++);
1217 }
1218 }
1219 }
1220
1222
1223 template <int dim, int spacedim, typename Number>
1224 void
1226 const ArrayView<Number> & dof_values,
1227 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1228 const std::vector<
1230 &shape_function_data,
1231 std::vector<typename SymmetricTensor<2, dim, spacedim>::
1232 template solution_divergence_type<Number>> &divergences)
1233 {
1234 const unsigned int dofs_per_cell = dof_values.size();
1235 const unsigned int n_quadrature_points = divergences.size();
1236
1237 std::fill(divergences.begin(),
1238 divergences.end(),
1240 template solution_divergence_type<Number>());
1241
1242 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1243 ++shape_function)
1244 {
1245 const int snc =
1246 shape_function_data[shape_function].single_nonzero_component;
1247
1248 if (snc == -2)
1249 // shape function is zero for the selected components
1250 continue;
1251
1252 const Number &value = dof_values[shape_function];
1253 // For auto-differentiable numbers, the fact that a DoF value is zero
1254 // does not imply that its derivatives are zero as well. So we
1255 // can't filter by value for these number types.
1256 if (::internal::CheckForZero<Number>::value(value) == true)
1257 continue;
1258
1259 if (snc != -1)
1260 {
1261 const unsigned int comp = shape_function_data[shape_function]
1262 .single_nonzero_component_index;
1263
1264 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1265 &shape_gradients[snc][0];
1266
1267 const unsigned int ii = ::SymmetricTensor<2, spacedim>::
1269 const unsigned int jj = ::SymmetricTensor<2, spacedim>::
1271
1272 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1273 ++q_point, ++shape_gradient_ptr)
1275 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1276
1277 if (ii != jj)
1278 divergences[q_point][jj] +=
1279 value * (*shape_gradient_ptr)[ii];
1280 }
1281 }
1282 else
1283 {
1284 for (unsigned int d = 0;
1285 d <
1287 spacedim>::n_independent_components;
1288 ++d)
1289 if (shape_function_data[shape_function]
1290 .is_nonzero_shape_function_component[d])
1291 {
1292 Assert(false, ExcNotImplemented());
1293
1294 // the following implementation needs to be looked over -- I
1295 // think it can't be right, because we are in a case where
1296 // there is no single nonzero component
1297 //
1298 // the following is not implemented! we need to consider the
1299 // interplay between multiple non-zero entries in shape
1300 // function and the representation as a symmetric
1301 // second-order tensor
1302 const unsigned int comp =
1303 shape_function_data[shape_function]
1304 .single_nonzero_component_index;
1305
1306 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1307 &shape_gradients[shape_function_data[shape_function]
1308 .row_index[d]][0];
1309 for (unsigned int q_point = 0;
1310 q_point < n_quadrature_points;
1311 ++q_point, ++shape_gradient_ptr)
1313 for (unsigned int j = 0; j < spacedim; ++j)
1314 {
1315 const unsigned int vector_component =
1318 TableIndices<2>(comp, j));
1319 divergences[q_point][vector_component] +=
1320 value * (*shape_gradient_ptr++)[j];
1321 }
1322 }
1323 }
1324 }
1325 }
1326 }
1327
1328 // ---------------------- non-symmetric tensor part ------------------------
1329
1330 template <int dim, int spacedim, typename Number>
1331 void
1333 const ArrayView<Number> & dof_values,
1334 const ::Table<2, double> &shape_values,
1335 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1336 &shape_function_data,
1337 std::vector<
1338 typename ProductType<Number, ::Tensor<2, spacedim>>::type>
1339 &values)
1340 {
1341 const unsigned int dofs_per_cell = dof_values.size();
1342 const unsigned int n_quadrature_points = values.size();
1343
1344 std::fill(
1345 values.begin(),
1346 values.end(),
1347 typename ProductType<Number, ::Tensor<2, spacedim>>::type());
1348
1349 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1350 ++shape_function)
1351 {
1352 const int snc =
1353 shape_function_data[shape_function].single_nonzero_component;
1354
1355 if (snc == -2)
1356 // shape function is zero for the selected components
1357 continue;
1358
1359 const Number &value = dof_values[shape_function];
1360 // For auto-differentiable numbers, the fact that a DoF value is zero
1361 // does not imply that its derivatives are zero as well. So we
1362 // can't filter by value for these number types.
1363 if (::internal::CheckForZero<Number>::value(value) == true)
1364 continue;
1365
1366 if (snc != -1)
1367 {
1368 const unsigned int comp = shape_function_data[shape_function]
1369 .single_nonzero_component_index;
1370
1371 const TableIndices<2> indices =
1373 comp);
1375 const double *shape_value_ptr = &shape_values(snc, 0);
1376 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1377 ++q_point)
1378 values[q_point][indices] += value * (*shape_value_ptr++);
1379 }
1380 else
1381 for (unsigned int d = 0; d < dim * dim; ++d)
1382 if (shape_function_data[shape_function]
1383 .is_nonzero_shape_function_component[d])
1384 {
1385 const TableIndices<2> indices =
1387 d);
1388
1389 const double *shape_value_ptr = &shape_values(
1390 shape_function_data[shape_function].row_index[d], 0);
1391 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1392 ++q_point)
1393 values[q_point][indices] += value * (*shape_value_ptr++);
1394 }
1395 }
1396 }
1397
1398
1399
1400 template <int dim, int spacedim, typename Number>
1401 void
1403 const ArrayView<Number> & dof_values,
1404 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1405 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1406 &shape_function_data,
1407 std::vector<typename Tensor<2, dim, spacedim>::
1408 template solution_divergence_type<Number>> &divergences)
1409 {
1410 const unsigned int dofs_per_cell = dof_values.size();
1411 const unsigned int n_quadrature_points = divergences.size();
1413 std::fill(
1414 divergences.begin(),
1415 divergences.end(),
1417 Number>());
1418
1419 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1420 ++shape_function)
1421 {
1422 const int snc =
1423 shape_function_data[shape_function].single_nonzero_component;
1424
1425 if (snc == -2)
1426 // shape function is zero for the selected components
1427 continue;
1428
1429 const Number &value = dof_values[shape_function];
1430 // For auto-differentiable numbers, the fact that a DoF value is zero
1431 // does not imply that its derivatives are zero as well. So we
1432 // can't filter by value for these number types.
1433 if (::internal::CheckForZero<Number>::value(value) == true)
1434 continue;
1435
1436 if (snc != -1)
1437 {
1438 const unsigned int comp = shape_function_data[shape_function]
1439 .single_nonzero_component_index;
1440
1441 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1442 &shape_gradients[snc][0];
1443
1444 const TableIndices<2> indices =
1446 comp);
1447 const unsigned int ii = indices[0];
1448 const unsigned int jj = indices[1];
1449
1450 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1451 ++q_point, ++shape_gradient_ptr)
1452 {
1453 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1454 }
1455 }
1456 else
1457 {
1458 for (unsigned int d = 0; d < dim * dim; ++d)
1459 if (shape_function_data[shape_function]
1460 .is_nonzero_shape_function_component[d])
1461 {
1462 Assert(false, ExcNotImplemented());
1463 }
1464 }
1465 }
1466 }
1467
1468
1469
1470 template <int dim, int spacedim, typename Number>
1471 void
1473 const ArrayView<Number> & dof_values,
1474 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1475 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1476 &shape_function_data,
1477 std::vector<typename Tensor<2, dim, spacedim>::
1478 template solution_gradient_type<Number>> &gradients)
1479 {
1480 const unsigned int dofs_per_cell = dof_values.size();
1481 const unsigned int n_quadrature_points = gradients.size();
1482
1483 std::fill(
1484 gradients.begin(),
1485 gradients.end(),
1486 typename Tensor<2, dim, spacedim>::template solution_gradient_type<
1487 Number>());
1488
1489 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1490 ++shape_function)
1491 {
1492 const int snc =
1493 shape_function_data[shape_function].single_nonzero_component;
1494
1495 if (snc == -2)
1496 // shape function is zero for the selected components
1497 continue;
1498
1499 const Number &value = dof_values[shape_function];
1500 // For auto-differentiable numbers, the fact that a DoF value is zero
1501 // does not imply that its derivatives are zero as well. So we
1502 // can't filter by value for these number types.
1503 if (::internal::CheckForZero<Number>::value(value) == true)
1504 continue;
1505
1506 if (snc != -1)
1507 {
1508 const unsigned int comp = shape_function_data[shape_function]
1509 .single_nonzero_component_index;
1510
1511 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1512 &shape_gradients[snc][0];
1513
1514 const TableIndices<2> indices =
1516 comp);
1517 const unsigned int ii = indices[0];
1518 const unsigned int jj = indices[1];
1519
1520 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1521 ++q_point, ++shape_gradient_ptr)
1522 {
1523 gradients[q_point][ii][jj] += value * (*shape_gradient_ptr);
1524 }
1525 }
1526 else
1527 {
1528 for (unsigned int d = 0; d < dim * dim; ++d)
1529 if (shape_function_data[shape_function]
1530 .is_nonzero_shape_function_component[d])
1531 {
1532 Assert(false, ExcNotImplemented());
1533 }
1534 }
1535 }
1536 }
1537
1538 } // end of namespace internal
1539
1540
1541
1542 template <int dim, int spacedim>
1543 template <class InputVector>
1544 void
1546 const InputVector &fe_function,
1548 const
1549 {
1550 Assert(fe_values->update_flags & update_values,
1552 "update_values")));
1553 Assert(fe_values->present_cell.is_initialized(),
1555 AssertDimension(fe_function.size(),
1557
1558 // get function values of dofs on this cell and call internal worker
1559 // function
1561 fe_values->dofs_per_cell);
1562 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1563 dof_values);
1564 internal::do_function_values<dim, spacedim>(
1565 make_array_view(dof_values.begin(), dof_values.end()),
1567 shape_function_data,
1568 values);
1569 }
1570
1571
1572
1573 template <int dim, int spacedim>
1574 template <class InputVector>
1575 void
1577 const InputVector &dof_values,
1579 const
1580 {
1581 Assert(fe_values->update_flags & update_values,
1583 "update_values")));
1584 Assert(fe_values->present_cell.is_initialized(),
1586 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1587
1588 internal::do_function_values<dim, spacedim>(
1589 make_array_view(dof_values.begin(), dof_values.end()),
1591 shape_function_data,
1592 values);
1593 }
1594
1595
1596
1597 template <int dim, int spacedim>
1598 template <class InputVector>
1599 void
1601 const InputVector &fe_function,
1603 &gradients) const
1604 {
1607 "update_gradients")));
1608 Assert(fe_values->present_cell.is_initialized(),
1610 AssertDimension(fe_function.size(),
1612
1613 // get function values of dofs on this cell
1615 fe_values->dofs_per_cell);
1616 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1617 dof_values);
1618 internal::do_function_derivatives<1, dim, spacedim>(
1619 make_array_view(dof_values.begin(), dof_values.end()),
1621 shape_function_data,
1622 gradients);
1623 }
1624
1625
1626
1627 template <int dim, int spacedim>
1628 template <class InputVector>
1629 void
1631 const InputVector &dof_values,
1633 &gradients) const
1634 {
1637 "update_gradients")));
1638 Assert(fe_values->present_cell.is_initialized(),
1640 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1641
1642 internal::do_function_derivatives<1, dim, spacedim>(
1643 make_array_view(dof_values.begin(), dof_values.end()),
1645 shape_function_data,
1646 gradients);
1647 }
1648
1649
1650
1651 template <int dim, int spacedim>
1652 template <class InputVector>
1653 void
1655 const InputVector &fe_function,
1657 &hessians) const
1658 {
1659 Assert(fe_values->update_flags & update_hessians,
1661 "update_hessians")));
1662 Assert(fe_values->present_cell.is_initialized(),
1664 AssertDimension(fe_function.size(),
1666
1667 // get function values of dofs on this cell
1669 fe_values->dofs_per_cell);
1670 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1671 dof_values);
1672 internal::do_function_derivatives<2, dim, spacedim>(
1673 make_array_view(dof_values.begin(), dof_values.end()),
1675 shape_function_data,
1676 hessians);
1677 }
1678
1679
1680
1681 template <int dim, int spacedim>
1682 template <class InputVector>
1683 void
1685 const InputVector &dof_values,
1687 &hessians) const
1688 {
1689 Assert(fe_values->update_flags & update_hessians,
1691 "update_hessians")));
1692 Assert(fe_values->present_cell.is_initialized(),
1694 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1695
1696 internal::do_function_derivatives<2, dim, spacedim>(
1697 make_array_view(dof_values.begin(), dof_values.end()),
1699 shape_function_data,
1700 hessians);
1701 }
1702
1703
1704
1705 template <int dim, int spacedim>
1706 template <class InputVector>
1707 void
1709 const InputVector &fe_function,
1711 &laplacians) const
1712 {
1713 Assert(fe_values->update_flags & update_hessians,
1715 "update_hessians")));
1716 Assert(fe_values->present_cell.is_initialized(),
1718 AssertDimension(fe_function.size(),
1720
1721 // get function values of dofs on this cell
1723 fe_values->dofs_per_cell);
1724 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1725 dof_values);
1726 internal::do_function_laplacians<dim, spacedim>(
1727 make_array_view(dof_values.begin(), dof_values.end()),
1729 shape_function_data,
1730 laplacians);
1731 }
1732
1733
1734
1735 template <int dim, int spacedim>
1736 template <class InputVector>
1737 void
1739 const InputVector &dof_values,
1741 &laplacians) const
1742 {
1743 Assert(fe_values->update_flags & update_hessians,
1745 "update_hessians")));
1746 Assert(fe_values->present_cell.is_initialized(),
1748 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1749
1750 internal::do_function_laplacians<dim, spacedim>(
1751 make_array_view(dof_values.begin(), dof_values.end()),
1753 shape_function_data,
1754 laplacians);
1755 }
1756
1757
1758
1759 template <int dim, int spacedim>
1760 template <class InputVector>
1761 void
1763 const InputVector &fe_function,
1764 std::vector<
1766 &third_derivatives) const
1767 {
1770 "update_3rd_derivatives")));
1771 Assert(fe_values->present_cell.is_initialized(),
1773 AssertDimension(fe_function.size(),
1775
1776 // get function values of dofs on this cell
1778 fe_values->dofs_per_cell);
1779 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1780 dof_values);
1781 internal::do_function_derivatives<3, dim, spacedim>(
1782 make_array_view(dof_values.begin(), dof_values.end()),
1784 shape_function_data,
1785 third_derivatives);
1786 }
1787
1788
1789
1790 template <int dim, int spacedim>
1791 template <class InputVector>
1792 void
1794 const InputVector &dof_values,
1795 std::vector<
1797 &third_derivatives) const
1798 {
1801 "update_3rd_derivatives")));
1802 Assert(fe_values->present_cell.is_initialized(),
1804 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1805
1806 internal::do_function_derivatives<3, dim, spacedim>(
1807 make_array_view(dof_values.begin(), dof_values.end()),
1809 shape_function_data,
1810 third_derivatives);
1811 }
1812
1813
1814
1815 template <int dim, int spacedim>
1816 template <class InputVector>
1817 void
1819 const InputVector &fe_function,
1821 const
1822 {
1823 Assert(fe_values->update_flags & update_values,
1825 "update_values")));
1826 Assert(fe_values->present_cell.is_initialized(),
1828 AssertDimension(fe_function.size(),
1830
1831 // get function values of dofs on this cell
1833 fe_values->dofs_per_cell);
1834 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1835 dof_values);
1836 internal::do_function_values<dim, spacedim>(
1837 make_array_view(dof_values.begin(), dof_values.end()),
1839 shape_function_data,
1840 values);
1841 }
1842
1843
1844
1845 template <int dim, int spacedim>
1846 template <class InputVector>
1847 void
1849 const InputVector &dof_values,
1851 const
1852 {
1853 Assert(fe_values->update_flags & update_values,
1855 "update_values")));
1856 Assert(fe_values->present_cell.is_initialized(),
1858 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1859
1860 internal::do_function_values<dim, spacedim>(
1861 make_array_view(dof_values.begin(), dof_values.end()),
1863 shape_function_data,
1864 values);
1865 }
1866
1867
1868
1869 template <int dim, int spacedim>
1870 template <class InputVector>
1871 void
1873 const InputVector &fe_function,
1875 &gradients) const
1876 {
1879 "update_gradients")));
1880 Assert(fe_values->present_cell.is_initialized(),
1882 AssertDimension(fe_function.size(),
1884
1885 // get function values of dofs on this cell
1887 fe_values->dofs_per_cell);
1888 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1889 dof_values);
1890 internal::do_function_derivatives<1, dim, spacedim>(
1891 make_array_view(dof_values.begin(), dof_values.end()),
1893 shape_function_data,
1894 gradients);
1895 }
1896
1897
1898
1899 template <int dim, int spacedim>
1900 template <class InputVector>
1901 void
1903 const InputVector &dof_values,
1905 &gradients) const
1906 {
1909 "update_gradients")));
1910 Assert(fe_values->present_cell.is_initialized(),
1912 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1913
1914 internal::do_function_derivatives<1, dim, spacedim>(
1915 make_array_view(dof_values.begin(), dof_values.end()),
1917 shape_function_data,
1918 gradients);
1919 }
1920
1921
1922
1923 template <int dim, int spacedim>
1924 template <class InputVector>
1925 void
1927 const InputVector &fe_function,
1928 std::vector<
1930 &symmetric_gradients) const
1931 {
1934 "update_gradients")));
1935 Assert(fe_values->present_cell.is_initialized(),
1937 AssertDimension(fe_function.size(),
1939
1940 // get function values of dofs on this cell
1942 fe_values->dofs_per_cell);
1943 fe_values->present_cell.get_interpolated_dof_values(fe_function,
1944 dof_values);
1945 internal::do_function_symmetric_gradients<dim, spacedim>(
1946 make_array_view(dof_values.begin(), dof_values.end()),
1948 shape_function_data,
1949 symmetric_gradients);
1950 }
1951
1952
1953
1954 template <int dim, int spacedim>
1955 template <class InputVector>
1956 void
1958 const InputVector &dof_values,
1959 std::vector<
1961 &symmetric_gradients) const
1962 {
1965 "update_gradients")));
1966 Assert(fe_values->present_cell.is_initialized(),
1968 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1969
1970 internal::do_function_symmetric_gradients<dim, spacedim>(
1971 make_array_view(dof_values.begin(), dof_values.end()),
1973 shape_function_data,
1974 symmetric_gradients);
1975 }
1976
1977
1978
1979 template <int dim, int spacedim>
1980 template <class InputVector>
1981 void
1983 const InputVector &fe_function,
1985 &divergences) const
1986 {
1989 "update_gradients")));
1990 Assert(fe_values->present_cell.is_initialized(),
1992 AssertDimension(fe_function.size(),
1994
1995 // get function values of dofs
1996 // on this cell
1998 fe_values->dofs_per_cell);
1999 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2000 dof_values);
2001 internal::do_function_divergences<dim, spacedim>(
2002 make_array_view(dof_values.begin(), dof_values.end()),
2004 shape_function_data,
2005 divergences);
2006 }
2007
2008
2009
2010 template <int dim, int spacedim>
2011 template <class InputVector>
2012 void
2014 const InputVector &dof_values,
2016 &divergences) const
2017 {
2020 "update_gradients")));
2021 Assert(fe_values->present_cell.is_initialized(),
2023 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2024
2025 internal::do_function_divergences<dim, spacedim>(
2026 make_array_view(dof_values.begin(), dof_values.end()),
2028 shape_function_data,
2029 divergences);
2030 }
2031
2032
2033
2034 template <int dim, int spacedim>
2035 template <class InputVector>
2036 void
2038 const InputVector &fe_function,
2040 const
2041 {
2044 "update_gradients")));
2045 Assert(fe_values->present_cell.is_initialized(),
2046 ExcMessage("FEValues object is not reinited to any cell"));
2047 AssertDimension(fe_function.size(),
2049
2050 // get function values of dofs on this cell
2052 fe_values->dofs_per_cell);
2053 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2054 dof_values);
2055 internal::do_function_curls<dim, spacedim>(
2056 make_array_view(dof_values.begin(), dof_values.end()),
2058 shape_function_data,
2059 curls);
2060 }
2061
2062
2063
2064 template <int dim, int spacedim>
2065 template <class InputVector>
2066 void
2068 const InputVector &dof_values,
2070 const
2071 {
2074 "update_gradients")));
2075 Assert(fe_values->present_cell.is_initialized(),
2076 ExcMessage("FEValues object is not reinited to any cell"));
2077 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2078
2079 internal::do_function_curls<dim, spacedim>(
2080 make_array_view(dof_values.begin(), dof_values.end()),
2082 shape_function_data,
2083 curls);
2084 }
2085
2086
2087
2088 template <int dim, int spacedim>
2089 template <class InputVector>
2090 void
2092 const InputVector &fe_function,
2094 &hessians) const
2095 {
2096 Assert(fe_values->update_flags & update_hessians,
2098 "update_hessians")));
2099 Assert(fe_values->present_cell.is_initialized(),
2101 AssertDimension(fe_function.size(),
2103
2104 // get function values of dofs on this cell
2106 fe_values->dofs_per_cell);
2107 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2108 dof_values);
2109 internal::do_function_derivatives<2, dim, spacedim>(
2110 make_array_view(dof_values.begin(), dof_values.end()),
2112 shape_function_data,
2113 hessians);
2114 }
2115
2116
2117
2118 template <int dim, int spacedim>
2119 template <class InputVector>
2120 void
2122 const InputVector &dof_values,
2124 &hessians) const
2125 {
2126 Assert(fe_values->update_flags & update_hessians,
2128 "update_hessians")));
2129 Assert(fe_values->present_cell.is_initialized(),
2131 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2132
2133 internal::do_function_derivatives<2, dim, spacedim>(
2134 make_array_view(dof_values.begin(), dof_values.end()),
2136 shape_function_data,
2137 hessians);
2138 }
2139
2140
2141
2142 template <int dim, int spacedim>
2143 template <class InputVector>
2144 void
2146 const InputVector &fe_function,
2148 &laplacians) const
2149 {
2150 Assert(fe_values->update_flags & update_hessians,
2152 "update_hessians")));
2153 Assert(laplacians.size() == fe_values->n_quadrature_points,
2154 ExcDimensionMismatch(laplacians.size(),
2155 fe_values->n_quadrature_points));
2156 Assert(fe_values->present_cell.is_initialized(),
2158 Assert(
2159 fe_function.size() == fe_values->present_cell.n_dofs_for_dof_handler(),
2160 ExcDimensionMismatch(fe_function.size(),
2161 fe_values->present_cell.n_dofs_for_dof_handler()));
2162
2163 // get function values of dofs on this cell
2165 fe_values->dofs_per_cell);
2166 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2167 dof_values);
2168 internal::do_function_laplacians<dim, spacedim>(
2169 make_array_view(dof_values.begin(), dof_values.end()),
2171 shape_function_data,
2172 laplacians);
2173 }
2174
2175
2176
2177 template <int dim, int spacedim>
2178 template <class InputVector>
2179 void
2181 const InputVector &dof_values,
2183 &laplacians) const
2184 {
2185 Assert(fe_values->update_flags & update_hessians,
2187 "update_hessians")));
2188 Assert(laplacians.size() == fe_values->n_quadrature_points,
2189 ExcDimensionMismatch(laplacians.size(),
2190 fe_values->n_quadrature_points));
2191 Assert(fe_values->present_cell.is_initialized(),
2193 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2194
2195 internal::do_function_laplacians<dim, spacedim>(
2196 make_array_view(dof_values.begin(), dof_values.end()),
2198 shape_function_data,
2199 laplacians);
2200 }
2201
2202
2203
2204 template <int dim, int spacedim>
2205 template <class InputVector>
2206 void
2208 const InputVector &fe_function,
2209 std::vector<
2211 &third_derivatives) const
2212 {
2215 "update_3rd_derivatives")));
2216 Assert(fe_values->present_cell.is_initialized(),
2218 AssertDimension(fe_function.size(),
2220
2221 // get function values of dofs on this cell
2223 fe_values->dofs_per_cell);
2224 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2225 dof_values);
2226 internal::do_function_derivatives<3, dim, spacedim>(
2227 make_array_view(dof_values.begin(), dof_values.end()),
2229 shape_function_data,
2230 third_derivatives);
2231 }
2232
2233
2234
2235 template <int dim, int spacedim>
2236 template <class InputVector>
2237 void
2239 const InputVector &dof_values,
2240 std::vector<
2242 &third_derivatives) const
2243 {
2246 "update_3rd_derivatives")));
2247 Assert(fe_values->present_cell.is_initialized(),
2249 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2250
2251 internal::do_function_derivatives<3, dim, spacedim>(
2252 make_array_view(dof_values.begin(), dof_values.end()),
2254 shape_function_data,
2255 third_derivatives);
2256 }
2257
2258
2259
2260 template <int dim, int spacedim>
2261 template <class InputVector>
2262 void
2264 const InputVector &fe_function,
2266 const
2267 {
2268 Assert(fe_values->update_flags & update_values,
2270 "update_values")));
2271 Assert(fe_values->present_cell.is_initialized(),
2273 AssertDimension(fe_function.size(),
2275
2276 // get function values of dofs on this cell
2278 fe_values->dofs_per_cell);
2279 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2280 dof_values);
2281 internal::do_function_values<dim, spacedim>(
2282 make_array_view(dof_values.begin(), dof_values.end()),
2284 shape_function_data,
2285 values);
2286 }
2287
2288
2289
2290 template <int dim, int spacedim>
2291 template <class InputVector>
2292 void
2294 const InputVector &dof_values,
2296 const
2297 {
2298 Assert(fe_values->update_flags & update_values,
2300 "update_values")));
2301 Assert(fe_values->present_cell.is_initialized(),
2303 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2304
2305 internal::do_function_values<dim, spacedim>(
2306 make_array_view(dof_values.begin(), dof_values.end()),
2308 shape_function_data,
2309 values);
2310 }
2311
2312
2313
2314 template <int dim, int spacedim>
2315 template <class InputVector>
2316 void
2318 const InputVector &fe_function,
2320 &divergences) const
2321 {
2324 "update_gradients")));
2325 Assert(fe_values->present_cell.is_initialized(),
2327 AssertDimension(fe_function.size(),
2329
2330 // get function values of dofs
2331 // on this cell
2333 fe_values->dofs_per_cell);
2334 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2335 dof_values);
2336 internal::do_function_divergences<dim, spacedim>(
2337 make_array_view(dof_values.begin(), dof_values.end()),
2339 shape_function_data,
2340 divergences);
2341 }
2342
2343
2344
2345 template <int dim, int spacedim>
2346 template <class InputVector>
2347 void
2350 const InputVector &dof_values,
2352 &divergences) const
2353 {
2356 "update_gradients")));
2357 Assert(fe_values->present_cell.is_initialized(),
2359 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2360
2361 internal::do_function_divergences<dim, spacedim>(
2362 make_array_view(dof_values.begin(), dof_values.end()),
2364 shape_function_data,
2365 divergences);
2366 }
2367
2368
2369
2370 template <int dim, int spacedim>
2371 template <class InputVector>
2372 void
2374 const InputVector &fe_function,
2376 const
2377 {
2378 Assert(fe_values->update_flags & update_values,
2380 "update_values")));
2381 Assert(fe_values->present_cell.is_initialized(),
2383 AssertDimension(fe_function.size(),
2385
2386 // get function values of dofs on this cell
2388 fe_values->dofs_per_cell);
2389 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2390 dof_values);
2391 internal::do_function_values<dim, spacedim>(
2392 make_array_view(dof_values.begin(), dof_values.end()),
2394 shape_function_data,
2395 values);
2396 }
2397
2398
2399
2400 template <int dim, int spacedim>
2401 template <class InputVector>
2402 void
2404 const InputVector &dof_values,
2406 const
2407 {
2408 Assert(fe_values->update_flags & update_values,
2410 "update_values")));
2411 Assert(fe_values->present_cell.is_initialized(),
2413 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2414
2415 internal::do_function_values<dim, spacedim>(
2416 make_array_view(dof_values.begin(), dof_values.end()),
2418 shape_function_data,
2419 values);
2420 }
2421
2422
2423
2424 template <int dim, int spacedim>
2425 template <class InputVector>
2426 void
2428 const InputVector &fe_function,
2430 &divergences) const
2431 {
2434 "update_gradients")));
2435 Assert(fe_values->present_cell.is_initialized(),
2437 AssertDimension(fe_function.size(),
2439
2440 // get function values of dofs
2441 // on this cell
2443 fe_values->dofs_per_cell);
2444 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2445 dof_values);
2446 internal::do_function_divergences<dim, spacedim>(
2447 make_array_view(dof_values.begin(), dof_values.end()),
2449 shape_function_data,
2450 divergences);
2451 }
2452
2453
2454
2455 template <int dim, int spacedim>
2456 template <class InputVector>
2457 void
2459 const InputVector &dof_values,
2461 &divergences) const
2462 {
2465 "update_gradients")));
2466 Assert(fe_values->present_cell.is_initialized(),
2468 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2469
2470 internal::do_function_divergences<dim, spacedim>(
2471 make_array_view(dof_values.begin(), dof_values.end()),
2473 shape_function_data,
2474 divergences);
2475 }
2476
2477
2478
2479 template <int dim, int spacedim>
2480 template <class InputVector>
2481 void
2483 const InputVector &fe_function,
2484 std::vector<solution_gradient_type<typename InputVector::value_type>>
2485 &gradients) const
2486 {
2489 "update_gradients")));
2490 Assert(fe_values->present_cell.is_initialized(),
2492 AssertDimension(fe_function.size(),
2494
2495 // get function values of dofs
2496 // on this cell
2498 fe_values->dofs_per_cell);
2499 fe_values->present_cell.get_interpolated_dof_values(fe_function,
2500 dof_values);
2501 internal::do_function_gradients<dim, spacedim>(
2502 make_array_view(dof_values.begin(), dof_values.end()),
2504 shape_function_data,
2505 gradients);
2506 }
2507
2508
2509
2510 template <int dim, int spacedim>
2511 template <class InputVector>
2512 void
2514 const InputVector &dof_values,
2516 &gradients) const
2517 {
2520 "update_gradients")));
2521 Assert(fe_values->present_cell.is_initialized(),
2523 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2524
2525 internal::do_function_gradients<dim, spacedim>(
2526 make_array_view(dof_values.begin(), dof_values.end()),
2528 shape_function_data,
2529 gradients);
2530 }
2531
2532} // namespace FEValuesViews
2533
2534
2535namespace internal
2536{
2537 namespace FEValuesViews
2538 {
2539 template <int dim, int spacedim>
2541 {
2542 const FiniteElement<dim, spacedim> &fe = fe_values.get_fe();
2543
2544 const unsigned int n_scalars = fe.n_components();
2545 scalars.reserve(n_scalars);
2546 for (unsigned int component = 0; component < n_scalars; ++component)
2547 scalars.emplace_back(fe_values, component);
2548
2549 // compute number of vectors that we can fit into this finite element.
2550 // note that this is based on the dimensionality 'dim' of the manifold,
2551 // not 'spacedim' of the output vector
2552 const unsigned int n_vectors =
2555 1 :
2556 0);
2557 vectors.reserve(n_vectors);
2558 for (unsigned int component = 0; component < n_vectors; ++component)
2559 vectors.emplace_back(fe_values, component);
2560
2561 // compute number of symmetric tensors in the same way as above
2562 const unsigned int n_symmetric_second_order_tensors =
2563 (fe.n_components() >=
2565 fe.n_components() -
2567 0);
2568 symmetric_second_order_tensors.reserve(n_symmetric_second_order_tensors);
2569 for (unsigned int component = 0;
2570 component < n_symmetric_second_order_tensors;
2571 ++component)
2572 symmetric_second_order_tensors.emplace_back(fe_values, component);
2573
2574
2575 // compute number of symmetric tensors in the same way as above
2576 const unsigned int n_second_order_tensors =
2579 1 :
2580 0);
2581 second_order_tensors.reserve(n_second_order_tensors);
2582 for (unsigned int component = 0; component < n_second_order_tensors;
2583 ++component)
2584 second_order_tensors.emplace_back(fe_values, component);
2585 }
2586 } // namespace FEValuesViews
2587} // namespace internal
2588
2589
2590/* ---------------- FEValuesBase<dim,spacedim>::CellIteratorContainer ---------
2591 */
2592
2593template <int dim, int spacedim>
2595 : initialized(false)
2596 , cell(typename Triangulation<dim, spacedim>::cell_iterator(nullptr, -1, -1))
2597 , dof_handler(nullptr)
2598 , level_dof_access(false)
2599{}
2600
2601
2602
2603template <int dim, int spacedim>
2606 : initialized(true)
2607 , cell(cell)
2608 , dof_handler(nullptr)
2609 , level_dof_access(false)
2610{}
2611
2612
2613
2614template <int dim, int spacedim>
2615bool
2617{
2618 return initialized;
2619}
2620
2621
2622
2623template <int dim, int spacedim>
2625operator typename Triangulation<dim, spacedim>::cell_iterator() const
2626{
2627 Assert(is_initialized(), ExcNotReinited());
2628
2629 return cell;
2630}
2631
2632
2633
2634template <int dim, int spacedim>
2637 const
2638{
2639 Assert(is_initialized(), ExcNotReinited());
2640 Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
2641
2642 return dof_handler->n_dofs();
2643}
2644
2645
2646
2647template <int dim, int spacedim>
2648template <typename VectorType>
2649void
2651 const VectorType & in,
2653{
2654 Assert(is_initialized(), ExcNotReinited());
2655 Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
2656
2657 if (level_dof_access)
2658 DoFCellAccessor<dim, spacedim, true>(&cell->get_triangulation(),
2659 cell->level(),
2660 cell->index(),
2661 dof_handler)
2663 else
2664 DoFCellAccessor<dim, spacedim, false>(&cell->get_triangulation(),
2665 cell->level(),
2666 cell->index(),
2667 dof_handler)
2669}
2670
2671
2672
2673template <int dim, int spacedim>
2674void
2676 const IndexSet & in,
2678{
2679 Assert(is_initialized(), ExcNotReinited());
2680 Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
2681 Assert(level_dof_access == false, ExcNotImplemented());
2682
2684 &cell->get_triangulation(), cell->level(), cell->index(), dof_handler);
2685
2686 std::vector<types::global_dof_index> dof_indices(
2687 cell_dofs.get_fe().n_dofs_per_cell());
2688 cell_dofs.get_dof_indices(dof_indices);
2689
2690 for (unsigned int i = 0; i < cell_dofs.get_fe().n_dofs_per_cell(); ++i)
2691 out[i] = (in.is_element(dof_indices[i]) ? 1 : 0);
2692}
2693
2694
2695
2696namespace internal
2697{
2698 namespace FEValuesImplementation
2699 {
2700 template <int dim, int spacedim>
2701 void
2703 const unsigned int n_quadrature_points,
2705 const UpdateFlags flags)
2706 {
2707 // initialize the table mapping from shape function number to
2708 // the rows in the tables storing the data by shape function and
2709 // nonzero component
2710 this->shape_function_to_row_table =
2712
2713 // count the total number of non-zero components accumulated
2714 // over all shape functions
2715 unsigned int n_nonzero_shape_components = 0;
2716 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
2717 n_nonzero_shape_components += fe.n_nonzero_components(i);
2718 Assert(n_nonzero_shape_components >= fe.n_dofs_per_cell(),
2720
2721 // with the number of rows now known, initialize those fields
2722 // that we will need to their correct size
2723 if (flags & update_values)
2724 {
2725 this->shape_values.reinit(n_nonzero_shape_components,
2726 n_quadrature_points);
2727 this->shape_values.fill(numbers::signaling_nan<double>());
2728 }
2729
2730 if (flags & update_gradients)
2731 {
2732 this->shape_gradients.reinit(n_nonzero_shape_components,
2733 n_quadrature_points);
2734 this->shape_gradients.fill(
2736 }
2737
2738 if (flags & update_hessians)
2739 {
2740 this->shape_hessians.reinit(n_nonzero_shape_components,
2741 n_quadrature_points);
2742 this->shape_hessians.fill(
2744 }
2745
2746 if (flags & update_3rd_derivatives)
2747 {
2748 this->shape_3rd_derivatives.reinit(n_nonzero_shape_components,
2749 n_quadrature_points);
2750 this->shape_3rd_derivatives.fill(
2752 }
2753 }
2754
2755
2756
2757 template <int dim, int spacedim>
2758 std::size_t
2760 {
2761 return (
2763 MemoryConsumption::memory_consumption(shape_gradients) +
2765 MemoryConsumption::memory_consumption(shape_3rd_derivatives) +
2766 MemoryConsumption::memory_consumption(shape_function_to_row_table));
2767 }
2768 } // namespace FEValuesImplementation
2769} // namespace internal
2770
2771
2772
2773/*------------------------------- FEValuesBase ---------------------------*/
2774
2775
2776template <int dim, int spacedim>
2778 const unsigned int n_q_points,
2779 const unsigned int dofs_per_cell,
2780 const UpdateFlags flags,
2783 : n_quadrature_points(n_q_points)
2784 , max_n_quadrature_points(n_q_points)
2786 , mapping(&mapping, typeid(*this).name())
2787 , fe(&fe, typeid(*this).name())
2788 , cell_similarity(CellSimilarity::Similarity::none)
2789 , fe_values_views_cache(*this)
2790{
2791 Assert(n_q_points > 0,
2792 ExcMessage("There is nothing useful you can do with an FEValues "
2793 "object when using a quadrature formula with zero "
2794 "quadrature points!"));
2795 this->update_flags = flags;
2796}
2797
2798
2799
2800template <int dim, int spacedim>
2802{
2803 tria_listener_refinement.disconnect();
2804 tria_listener_mesh_transform.disconnect();
2805}
2806
2807
2808
2809namespace internal
2810{
2811 // put shape function part of get_function_xxx methods into separate
2812 // internal functions. this allows us to reuse the same code for several
2813 // functions (e.g. both the versions with and without indices) as well as
2814 // the same code for gradients and Hessians. Moreover, this speeds up
2815 // compilation and reduces the size of the final file since all the
2816 // different global vectors get channeled through the same code.
2817
2818 template <typename Number, typename Number2>
2819 void
2821 const ::Table<2, double> &shape_values,
2822 std::vector<Number> & values)
2823 {
2824 // scalar finite elements, so shape_values.size() == dofs_per_cell
2825 const unsigned int dofs_per_cell = shape_values.n_rows();
2826 const unsigned int n_quadrature_points = values.size();
2827
2828 // initialize with zero
2829 std::fill_n(values.begin(),
2830 n_quadrature_points,
2832
2833 // add up contributions of trial functions. note that here we deal with
2834 // scalar finite elements, so no need to check for non-primitivity of
2835 // shape functions. in order to increase the speed of this function, we
2836 // directly access the data in the shape_values array, and increment
2837 // pointers for accessing the data. this saves some lookup time and
2838 // indexing. moreover, the order of the loops is such that we can access
2839 // the shape_values data stored contiguously
2840 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
2841 {
2842 const Number2 value = dof_values[shape_func];
2843 // For auto-differentiable numbers, the fact that a DoF value is zero
2844 // does not imply that its derivatives are zero as well. So we
2845 // can't filter by value for these number types.
2848 continue;
2849
2850 const double *shape_value_ptr = &shape_values(shape_func, 0);
2851 for (unsigned int point = 0; point < n_quadrature_points; ++point)
2852 values[point] += value * (*shape_value_ptr++);
2853 }
2854 }
2856
2857
2858 template <int dim, int spacedim, typename VectorType>
2859 void
2862 const ::Table<2, double> & shape_values,
2864 const std::vector<unsigned int> &shape_function_to_row_table,
2865 ArrayView<VectorType> values,
2866 const bool quadrature_points_fastest = false,
2867 const unsigned int component_multiple = 1)
2868 {
2869 using Number = typename VectorType::value_type;
2870 // initialize with zero
2871 for (unsigned int i = 0; i < values.size(); ++i)
2872 std::fill_n(values[i].begin(),
2873 values[i].size(),
2874 typename VectorType::value_type());
2875
2876 // see if there the current cell has DoFs at all, and if not
2877 // then there is nothing else to do.
2878 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
2879 if (dofs_per_cell == 0)
2880 return;
2881
2882 const unsigned int n_quadrature_points =
2883 quadrature_points_fastest ? values[0].size() : values.size();
2884 const unsigned int n_components = fe.n_components();
2885
2886 // Assert that we can write all components into the result vectors
2887 const unsigned result_components = n_components * component_multiple;
2888 (void)result_components;
2889 if (quadrature_points_fastest)
2890 {
2891 AssertDimension(values.size(), result_components);
2892 for (unsigned int i = 0; i < values.size(); ++i)
2893 AssertDimension(values[i].size(), n_quadrature_points);
2894 }
2895 else
2896 {
2897 AssertDimension(values.size(), n_quadrature_points);
2898 for (unsigned int i = 0; i < values.size(); ++i)
2899 AssertDimension(values[i].size(), result_components);
2900 }
2901
2902 // add up contributions of trial functions. now check whether the shape
2903 // function is primitive or not. if it is, then set its only non-zero
2904 // component, otherwise loop over components
2905 for (unsigned int mc = 0; mc < component_multiple; ++mc)
2906 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
2907 ++shape_func)
2908 {
2909 const Number &value = dof_values[shape_func + mc * dofs_per_cell];
2910 // For auto-differentiable numbers, the fact that a DoF value is zero
2911 // does not imply that its derivatives are zero as well. So we
2912 // can't filter by value for these number types.
2913 if (::internal::CheckForZero<Number>::value(value) == true)
2914 continue;
2915
2916 if (fe.is_primitive(shape_func))
2917 {
2918 const unsigned int comp =
2919 fe.system_to_component_index(shape_func).first +
2920 mc * n_components;
2921 const unsigned int row =
2922 shape_function_to_row_table[shape_func * n_components + comp];
2923
2924 const double *shape_value_ptr = &shape_values(row, 0);
2925
2926 if (quadrature_points_fastest)
2927 {
2928 VectorType &values_comp = values[comp];
2929 for (unsigned int point = 0; point < n_quadrature_points;
2930 ++point)
2931 values_comp[point] += value * (*shape_value_ptr++);
2933 else
2934 for (unsigned int point = 0; point < n_quadrature_points;
2935 ++point)
2936 values[point][comp] += value * (*shape_value_ptr++);
2937 }
2938 else
2939 for (unsigned int c = 0; c < n_components; ++c)
2940 {
2941 if (fe.get_nonzero_components(shape_func)[c] == false)
2942 continue;
2943
2944 const unsigned int row =
2945 shape_function_to_row_table[shape_func * n_components + c];
2946
2947 const double * shape_value_ptr = &shape_values(row, 0);
2948 const unsigned int comp = c + mc * n_components;
2949
2950 if (quadrature_points_fastest)
2951 {
2952 VectorType &values_comp = values[comp];
2953 for (unsigned int point = 0; point < n_quadrature_points;
2954 ++point)
2955 values_comp[point] += value * (*shape_value_ptr++);
2956 }
2957 else
2958 for (unsigned int point = 0; point < n_quadrature_points;
2959 ++point)
2960 values[point][comp] += value * (*shape_value_ptr++);
2961 }
2962 }
2963 }
2965
2966
2967 // use the same implementation for gradients and Hessians, distinguish them
2968 // by the rank of the tensors
2969 template <int order, int spacedim, typename Number>
2970 void
2972 const ArrayView<Number> & dof_values,
2973 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
2974 std::vector<Tensor<order, spacedim, Number>> & derivatives)
2975 {
2976 const unsigned int dofs_per_cell = shape_derivatives.size()[0];
2977 const unsigned int n_quadrature_points = derivatives.size();
2978
2979 // initialize with zero
2980 std::fill_n(derivatives.begin(),
2981 n_quadrature_points,
2983
2984 // add up contributions of trial functions. note that here we deal with
2985 // scalar finite elements, so no need to check for non-primitivity of
2986 // shape functions. in order to increase the speed of this function, we
2987 // directly access the data in the shape_gradients/hessians array, and
2988 // increment pointers for accessing the data. this saves some lookup time
2989 // and indexing. moreover, the order of the loops is such that we can
2990 // access the shape_gradients/hessians data stored contiguously
2991 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
2992 {
2993 const Number &value = dof_values[shape_func];
2994 // For auto-differentiable numbers, the fact that a DoF value is zero
2995 // does not imply that its derivatives are zero as well. So we
2996 // can't filter by value for these number types.
2997 if (::internal::CheckForZero<Number>::value(value) == true)
2998 continue;
2999
3000 const Tensor<order, spacedim> *shape_derivative_ptr =
3001 &shape_derivatives[shape_func][0];
3002 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3003 derivatives[point] += value * (*shape_derivative_ptr++);
3004 }
3005 }
3006
3007
3008
3009 template <int order, int dim, int spacedim, typename Number>
3010 void
3012 const ArrayView<Number> & dof_values,
3013 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
3015 const std::vector<unsigned int> &shape_function_to_row_table,
3016 ArrayView<std::vector<Tensor<order, spacedim, Number>>> derivatives,
3017 const bool quadrature_points_fastest = false,
3018 const unsigned int component_multiple = 1)
3019 {
3020 // initialize with zero
3021 for (unsigned int i = 0; i < derivatives.size(); ++i)
3022 std::fill_n(derivatives[i].begin(),
3023 derivatives[i].size(),
3025
3026 // see if there the current cell has DoFs at all, and if not
3027 // then there is nothing else to do.
3028 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3029 if (dofs_per_cell == 0)
3030 return;
3031
3032
3033 const unsigned int n_quadrature_points =
3034 quadrature_points_fastest ? derivatives[0].size() : derivatives.size();
3035 const unsigned int n_components = fe.n_components();
3036
3037 // Assert that we can write all components into the result vectors
3038 const unsigned result_components = n_components * component_multiple;
3039 (void)result_components;
3040 if (quadrature_points_fastest)
3042 AssertDimension(derivatives.size(), result_components);
3043 for (unsigned int i = 0; i < derivatives.size(); ++i)
3044 AssertDimension(derivatives[i].size(), n_quadrature_points);
3045 }
3046 else
3047 {
3048 AssertDimension(derivatives.size(), n_quadrature_points);
3049 for (unsigned int i = 0; i < derivatives.size(); ++i)
3050 AssertDimension(derivatives[i].size(), result_components);
3051 }
3052
3053 // add up contributions of trial functions. now check whether the shape
3054 // function is primitive or not. if it is, then set its only non-zero
3055 // component, otherwise loop over components
3056 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3057 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3058 ++shape_func)
3059 {
3060 const Number &value = dof_values[shape_func + mc * dofs_per_cell];
3061 // For auto-differentiable numbers, the fact that a DoF value is zero
3062 // does not imply that its derivatives are zero as well. So we
3063 // can't filter by value for these number types.
3064 if (::internal::CheckForZero<Number>::value(value) == true)
3065 continue;
3066
3067 if (fe.is_primitive(shape_func))
3068 {
3069 const unsigned int comp =
3070 fe.system_to_component_index(shape_func).first +
3071 mc * n_components;
3072 const unsigned int row =
3073 shape_function_to_row_table[shape_func * n_components + comp];
3075 const Tensor<order, spacedim> *shape_derivative_ptr =
3076 &shape_derivatives[row][0];
3077
3078 if (quadrature_points_fastest)
3079 for (unsigned int point = 0; point < n_quadrature_points;
3080 ++point)
3081 derivatives[comp][point] += value * (*shape_derivative_ptr++);
3082 else
3083 for (unsigned int point = 0; point < n_quadrature_points;
3084 ++point)
3085 derivatives[point][comp] += value * (*shape_derivative_ptr++);
3086 }
3087 else
3088 for (unsigned int c = 0; c < n_components; ++c)
3089 {
3090 if (fe.get_nonzero_components(shape_func)[c] == false)
3091 continue;
3092
3093 const unsigned int row =
3094 shape_function_to_row_table[shape_func * n_components + c];
3095
3096 const Tensor<order, spacedim> *shape_derivative_ptr =
3097 &shape_derivatives[row][0];
3098 const unsigned int comp = c + mc * n_components;
3099
3100 if (quadrature_points_fastest)
3101 for (unsigned int point = 0; point < n_quadrature_points;
3102 ++point)
3103 derivatives[comp][point] +=
3104 value * (*shape_derivative_ptr++);
3105 else
3106 for (unsigned int point = 0; point < n_quadrature_points;
3107 ++point)
3108 derivatives[point][comp] +=
3109 value * (*shape_derivative_ptr++);
3110 }
3111 }
3112 }
3113
3114
3115
3116 template <int spacedim, typename Number, typename Number2>
3117 void
3119 const ArrayView<Number2> & dof_values,
3120 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
3121 std::vector<Number> & laplacians)
3122 {
3123 const unsigned int dofs_per_cell = shape_hessians.size()[0];
3124 const unsigned int n_quadrature_points = laplacians.size();
3125
3126 // initialize with zero
3127 std::fill_n(laplacians.begin(),
3128 n_quadrature_points,
3130
3131 // add up contributions of trial functions. note that here we deal with
3132 // scalar finite elements and also note that the Laplacian is
3133 // the trace of the Hessian.
3134 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
3135 {
3136 const Number2 value = dof_values[shape_func];
3137 // For auto-differentiable numbers, the fact that a DoF value is zero
3138 // does not imply that its derivatives are zero as well. So we
3139 // can't filter by value for these number types.
3142 continue;
3143
3144 const Tensor<2, spacedim> *shape_hessian_ptr =
3145 &shape_hessians[shape_func][0];
3146 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3147 laplacians[point] += value * trace(*shape_hessian_ptr++);
3148 }
3150
3151
3152
3153 template <int dim, int spacedim, typename VectorType, typename Number>
3154 void
3156 const ArrayView<Number> & dof_values,
3157 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
3159 const std::vector<unsigned int> & shape_function_to_row_table,
3160 std::vector<VectorType> & laplacians,
3161 const bool quadrature_points_fastest = false,
3162 const unsigned int component_multiple = 1)
3164 // initialize with zero
3165 for (unsigned int i = 0; i < laplacians.size(); ++i)
3166 std::fill_n(laplacians[i].begin(),
3167 laplacians[i].size(),
3168 typename VectorType::value_type());
3169
3170 // see if there the current cell has DoFs at all, and if not
3171 // then there is nothing else to do.
3172 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3173 if (dofs_per_cell == 0)
3174 return;
3175
3176
3177 const unsigned int n_quadrature_points = laplacians.size();
3178 const unsigned int n_components = fe.n_components();
3179
3180 // Assert that we can write all components into the result vectors
3181 const unsigned result_components = n_components * component_multiple;
3182 (void)result_components;
3183 if (quadrature_points_fastest)
3184 {
3185 AssertDimension(laplacians.size(), result_components);
3186 for (unsigned int i = 0; i < laplacians.size(); ++i)
3187 AssertDimension(laplacians[i].size(), n_quadrature_points);
3188 }
3189 else
3190 {
3191 AssertDimension(laplacians.size(), n_quadrature_points);
3192 for (unsigned int i = 0; i < laplacians.size(); ++i)
3193 AssertDimension(laplacians[i].size(), result_components);
3194 }
3195
3196 // add up contributions of trial functions. now check whether the shape
3197 // function is primitive or not. if it is, then set its only non-zero
3198 // component, otherwise loop over components
3199 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3200 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3201 ++shape_func)
3202 {
3203 const Number &value = dof_values[shape_func + mc * dofs_per_cell];
3204 // For auto-differentiable numbers, the fact that a DoF value is zero
3205 // does not imply that its derivatives are zero as well. So we
3206 // can't filter by value for these number types.
3207 if (::internal::CheckForZero<Number>::value(value) == true)
3208 continue;
3209
3210 if (fe.is_primitive(shape_func))
3211 {
3212 const unsigned int comp =
3213 fe.system_to_component_index(shape_func).first +
3214 mc * n_components;
3215 const unsigned int row =
3216 shape_function_to_row_table[shape_func * n_components + comp];
3217
3218 const Tensor<2, spacedim> *shape_hessian_ptr =
3219 &shape_hessians[row][0];
3220 if (quadrature_points_fastest)
3221 {
3222 VectorType &laplacians_comp = laplacians[comp];
3223 for (unsigned int point = 0; point < n_quadrature_points;
3224 ++point)
3225 laplacians_comp[point] +=
3226 value * trace(*shape_hessian_ptr++);
3227 }
3228 else
3229 for (unsigned int point = 0; point < n_quadrature_points;
3230 ++point)
3231 laplacians[point][comp] +=
3232 value * trace(*shape_hessian_ptr++);
3233 }
3234 else
3235 for (unsigned int c = 0; c < n_components; ++c)
3236 {
3237 if (fe.get_nonzero_components(shape_func)[c] == false)
3238 continue;
3239
3240 const unsigned int row =
3241 shape_function_to_row_table[shape_func * n_components + c];
3242
3243 const Tensor<2, spacedim> *shape_hessian_ptr =
3244 &shape_hessians[row][0];
3245 const unsigned int comp = c + mc * n_components;
3246
3247 if (quadrature_points_fastest)
3248 {
3249 VectorType &laplacians_comp = laplacians[comp];
3250 for (unsigned int point = 0; point < n_quadrature_points;
3251 ++point)
3252 laplacians_comp[point] +=
3253 value * trace(*shape_hessian_ptr++);
3254 }
3255 else
3256 for (unsigned int point = 0; point < n_quadrature_points;
3257 ++point)
3258 laplacians[point][comp] +=
3259 value * trace(*shape_hessian_ptr++);
3260 }
3261 }
3262 }
3263} // namespace internal
3264
3265
3266
3267template <int dim, int spacedim>
3268template <class InputVector>
3269void
3271 const InputVector & fe_function,
3272 std::vector<typename InputVector::value_type> &values) const
3273{
3274 using Number = typename InputVector::value_type;
3275 Assert(this->update_flags & update_values,
3276 ExcAccessToUninitializedField("update_values"));
3277 AssertDimension(fe->n_components(), 1);
3278 Assert(present_cell.is_initialized(), ExcNotReinited());
3279 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3280
3281 // get function values of dofs on this cell
3282 Vector<Number> dof_values(dofs_per_cell);
3283 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3285 dof_values.end()),
3286 this->finite_element_output.shape_values,
3287 values);
3288}
3289
3290
3291
3292template <int dim, int spacedim>
3293template <class InputVector>
3294void
3296 const InputVector & fe_function,
3298 std::vector<typename InputVector::value_type> & values) const
3299{
3300 using Number = typename InputVector::value_type;
3301 Assert(this->update_flags & update_values,
3302 ExcAccessToUninitializedField("update_values"));
3303 AssertDimension(fe->n_components(), 1);
3304 AssertDimension(indices.size(), dofs_per_cell);
3305
3306 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3307 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3308 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3310 dof_values.end()),
3311 this->finite_element_output.shape_values,
3312 values);
3313}
3314
3315
3316
3317template <int dim, int spacedim>
3318template <class InputVector>
3319void
3321 const InputVector & fe_function,
3322 std::vector<Vector<typename InputVector::value_type>> &values) const
3323{
3324 using Number = typename InputVector::value_type;
3325 Assert(present_cell.is_initialized(), ExcNotReinited());
3326
3327 Assert(this->update_flags & update_values,
3328 ExcAccessToUninitializedField("update_values"));
3329 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3330
3331 // get function values of dofs on this cell
3332 Vector<Number> dof_values(dofs_per_cell);
3333 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3335 make_array_view(dof_values.begin(), dof_values.end()),
3336 this->finite_element_output.shape_values,
3337 *fe,
3338 this->finite_element_output.shape_function_to_row_table,
3339 make_array_view(values.begin(), values.end()));
3340}
3341
3342
3343
3344template <int dim, int spacedim>
3345template <class InputVector>
3346void
3348 const InputVector & fe_function,
3350 std::vector<Vector<typename InputVector::value_type>> &values) const
3351{
3352 using Number = typename InputVector::value_type;
3353 // Size of indices must be a multiple of dofs_per_cell such that an integer
3354 // number of function values is generated in each point.
3355 Assert(indices.size() % dofs_per_cell == 0,
3356 ExcNotMultiple(indices.size(), dofs_per_cell));
3357 Assert(this->update_flags & update_values,
3358 ExcAccessToUninitializedField("update_values"));
3359
3360 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3361 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3362 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3364 make_array_view(dof_values.begin(), dof_values.end()),
3365 this->finite_element_output.shape_values,
3366 *fe,
3367 this->finite_element_output.shape_function_to_row_table,
3368 make_array_view(values.begin(), values.end()),
3369 false,
3370 indices.size() / dofs_per_cell);
3371}
3372
3373
3374
3375template <int dim, int spacedim>
3376template <class InputVector>
3377void
3379 const InputVector & fe_function,
3381 ArrayView<std::vector<typename InputVector::value_type>> values,
3382 const bool quadrature_points_fastest) const
3383{
3384 using Number = typename InputVector::value_type;
3385 Assert(this->update_flags & update_values,
3386 ExcAccessToUninitializedField("update_values"));
3387
3388 // Size of indices must be a multiple of dofs_per_cell such that an integer
3389 // number of function values is generated in each point.
3390 Assert(indices.size() % dofs_per_cell == 0,
3391 ExcNotMultiple(indices.size(), dofs_per_cell));
3392
3393 boost::container::small_vector<Number, 200> dof_values(indices.size());
3394 for (unsigned int i = 0; i < indices.size(); ++i)
3395 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3397 make_array_view(dof_values.begin(), dof_values.end()),
3398 this->finite_element_output.shape_values,
3399 *fe,
3400 this->finite_element_output.shape_function_to_row_table,
3401 make_array_view(values.begin(), values.end()),
3402 quadrature_points_fastest,
3403 indices.size() / dofs_per_cell);
3404}
3405
3406
3407
3408template <int dim, int spacedim>
3409template <class InputVector>
3410void
3412 const InputVector &fe_function,
3414 const
3415{
3416 using Number = typename InputVector::value_type;
3417 Assert(this->update_flags & update_gradients,
3418 ExcAccessToUninitializedField("update_gradients"));
3419 AssertDimension(fe->n_components(), 1);
3420 Assert(present_cell.is_initialized(), ExcNotReinited());
3421 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3422
3423 // get function values of dofs on this cell
3424 Vector<Number> dof_values(dofs_per_cell);
3425 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3427 dof_values.end()),
3428 this->finite_element_output.shape_gradients,
3429 gradients);
3430}
3431
3432
3433
3434template <int dim, int spacedim>
3435template <class InputVector>
3436void
3438 const InputVector & fe_function,
3441 const
3442{
3443 using Number = typename InputVector::value_type;
3444 Assert(this->update_flags & update_gradients,
3445 ExcAccessToUninitializedField("update_gradients"));
3446 AssertDimension(fe->n_components(), 1);
3447 AssertDimension(indices.size(), dofs_per_cell);
3448
3449 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3450 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3451 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3453 dof_values.end()),
3454 this->finite_element_output.shape_gradients,
3455 gradients);
3456}
3457
3458
3459
3460template <int dim, int spacedim>
3461template <class InputVector>
3462void
3464 const InputVector &fe_function,
3465 std::vector<
3467 &gradients) const
3468{
3469 using Number = typename InputVector::value_type;
3470 Assert(this->update_flags & update_gradients,
3471 ExcAccessToUninitializedField("update_gradients"));
3472 Assert(present_cell.is_initialized(), ExcNotReinited());
3473 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3474
3475 // get function values of dofs on this cell
3476 Vector<Number> dof_values(dofs_per_cell);
3477 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3479 make_array_view(dof_values.begin(), dof_values.end()),
3480 this->finite_element_output.shape_gradients,
3481 *fe,
3482 this->finite_element_output.shape_function_to_row_table,
3483 make_array_view(gradients.begin(), gradients.end()));
3484}
3485
3486
3487
3488template <int dim, int spacedim>
3489template <class InputVector>
3490void
3492 const InputVector & fe_function,
3495 gradients,
3496 const bool quadrature_points_fastest) const
3497{
3498 using Number = typename InputVector::value_type;
3499 // Size of indices must be a multiple of dofs_per_cell such that an integer
3500 // number of function values is generated in each point.
3501 Assert(indices.size() % dofs_per_cell == 0,
3502 ExcNotMultiple(indices.size(), dofs_per_cell));
3503 Assert(this->update_flags & update_gradients,
3504 ExcAccessToUninitializedField("update_gradients"));
3505
3506 boost::container::small_vector<Number, 200> dof_values(indices.size());
3507 for (unsigned int i = 0; i < indices.size(); ++i)
3508 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3510 make_array_view(dof_values.begin(), dof_values.end()),
3511 this->finite_element_output.shape_gradients,
3512 *fe,
3513 this->finite_element_output.shape_function_to_row_table,
3514 make_array_view(gradients.begin(), gradients.end()),
3515 quadrature_points_fastest,
3516 indices.size() / dofs_per_cell);
3517}
3518
3519
3520
3521template <int dim, int spacedim>
3522template <class InputVector>
3523void
3525 const InputVector &fe_function,
3527 const
3528{
3529 using Number = typename InputVector::value_type;
3530 AssertDimension(fe->n_components(), 1);
3531 Assert(this->update_flags & update_hessians,
3532 ExcAccessToUninitializedField("update_hessians"));
3533 Assert(present_cell.is_initialized(), ExcNotReinited());
3534 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3535
3536 // get function values of dofs on this cell
3537 Vector<Number> dof_values(dofs_per_cell);
3538 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3540 dof_values.end()),
3541 this->finite_element_output.shape_hessians,
3542 hessians);
3543}
3544
3545
3546
3547template <int dim, int spacedim>
3548template <class InputVector>
3549void
3551 const InputVector & fe_function,
3554 const
3555{
3556 using Number = typename InputVector::value_type;
3557 Assert(this->update_flags & update_hessians,
3558 ExcAccessToUninitializedField("update_hessians"));
3559 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3560 AssertDimension(indices.size(), dofs_per_cell);
3561
3562 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3563 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3564 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3566 dof_values.end()),
3567 this->finite_element_output.shape_hessians,
3568 hessians);
3569}
3570
3571
3572
3573template <int dim, int spacedim>
3574template <class InputVector>
3575void
3577 const InputVector &fe_function,
3578 std::vector<
3580 & hessians,
3581 const bool quadrature_points_fastest) const
3582{
3583 using Number = typename InputVector::value_type;
3584 Assert(this->update_flags & update_hessians,
3585 ExcAccessToUninitializedField("update_hessians"));
3586 Assert(present_cell.is_initialized(), ExcNotReinited());
3587 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3588
3589 // get function values of dofs on this cell
3590 Vector<Number> dof_values(dofs_per_cell);
3591 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3593 make_array_view(dof_values.begin(), dof_values.end()),
3594 this->finite_element_output.shape_hessians,
3595 *fe,
3596 this->finite_element_output.shape_function_to_row_table,
3597 make_array_view(hessians.begin(), hessians.end()),
3598 quadrature_points_fastest);
3599}
3600
3601
3602
3603template <int dim, int spacedim>
3604template <class InputVector>
3605void
3607 const InputVector & fe_function,
3610 hessians,
3611 const bool quadrature_points_fastest) const
3612{
3613 using Number = typename InputVector::value_type;
3614 Assert(this->update_flags & update_hessians,
3615 ExcAccessToUninitializedField("update_hessians"));
3616 Assert(indices.size() % dofs_per_cell == 0,
3617 ExcNotMultiple(indices.size(), dofs_per_cell));
3618
3619 boost::container::small_vector<Number, 200> dof_values(indices.size());
3620 for (unsigned int i = 0; i < indices.size(); ++i)
3621 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3623 make_array_view(dof_values.begin(), dof_values.end()),
3624 this->finite_element_output.shape_hessians,
3625 *fe,
3626 this->finite_element_output.shape_function_to_row_table,
3627 make_array_view(hessians.begin(), hessians.end()),
3628 quadrature_points_fastest,
3629 indices.size() / dofs_per_cell);
3630}
3631
3632
3633
3634template <int dim, int spacedim>
3635template <class InputVector>
3636void
3638 const InputVector & fe_function,
3639 std::vector<typename InputVector::value_type> &laplacians) const
3640{
3641 using Number = typename InputVector::value_type;
3642 Assert(this->update_flags & update_hessians,
3643 ExcAccessToUninitializedField("update_hessians"));
3644 AssertDimension(fe->n_components(), 1);
3645 Assert(present_cell.is_initialized(), ExcNotReinited());
3646 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3647
3648 // get function values of dofs on this cell
3649 Vector<Number> dof_values(dofs_per_cell);
3650 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3652 dof_values.end()),
3653 this->finite_element_output.shape_hessians,
3654 laplacians);
3655}
3656
3657
3658
3659template <int dim, int spacedim>
3660template <class InputVector>
3661void
3663 const InputVector & fe_function,
3665 std::vector<typename InputVector::value_type> & laplacians) const
3666{
3667 using Number = typename InputVector::value_type;
3668 Assert(this->update_flags & update_hessians,
3669 ExcAccessToUninitializedField("update_hessians"));
3670 AssertDimension(fe->n_components(), 1);
3671 AssertDimension(indices.size(), dofs_per_cell);
3672
3673 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3674 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3675 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3677 dof_values.end()),
3678 this->finite_element_output.shape_hessians,
3679 laplacians);
3680}
3681
3682
3683
3684template <int dim, int spacedim>
3685template <class InputVector>
3686void
3688 const InputVector & fe_function,
3689 std::vector<Vector<typename InputVector::value_type>> &laplacians) const
3690{
3691 using Number = typename InputVector::value_type;
3692 Assert(present_cell.is_initialized(), ExcNotReinited());
3693 Assert(this->update_flags & update_hessians,
3694 ExcAccessToUninitializedField("update_hessians"));
3695 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3696
3697 // get function values of dofs on this cell
3698 Vector<Number> dof_values(dofs_per_cell);
3699 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3701 make_array_view(dof_values.begin(), dof_values.end()),
3702 this->finite_element_output.shape_hessians,
3703 *fe,
3704 this->finite_element_output.shape_function_to_row_table,
3705 laplacians);
3706}
3707
3708
3709
3710template <int dim, int spacedim>
3711template <class InputVector>
3712void
3714 const InputVector & fe_function,
3716 std::vector<Vector<typename InputVector::value_type>> &laplacians) const
3717{
3718 using Number = typename InputVector::value_type;
3719 // Size of indices must be a multiple of dofs_per_cell such that an integer
3720 // number of function values is generated in each point.
3721 Assert(indices.size() % dofs_per_cell == 0,
3722 ExcNotMultiple(indices.size(), dofs_per_cell));
3723 Assert(this->update_flags & update_hessians,
3724 ExcAccessToUninitializedField("update_hessians"));
3725
3726 boost::container::small_vector<Number, 200> dof_values(indices.size());
3727 for (unsigned int i = 0; i < indices.size(); ++i)
3728 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3730 make_array_view(dof_values.begin(), dof_values.end()),
3731 this->finite_element_output.shape_hessians,
3732 *fe,
3733 this->finite_element_output.shape_function_to_row_table,
3734 laplacians,
3735 false,
3736 indices.size() / dofs_per_cell);
3737}
3738
3739
3740
3741template <int dim, int spacedim>
3742template <class InputVector>
3743void
3745 const InputVector & fe_function,
3747 std::vector<std::vector<typename InputVector::value_type>> &laplacians,
3748 const bool quadrature_points_fastest) const
3749{
3750 using Number = typename InputVector::value_type;
3751 Assert(indices.size() % dofs_per_cell == 0,
3752 ExcNotMultiple(indices.size(), dofs_per_cell));
3753 Assert(this->update_flags & update_hessians,
3754 ExcAccessToUninitializedField("update_hessians"));
3755
3756 boost::container::small_vector<Number, 200> dof_values(indices.size());
3757 for (unsigned int i = 0; i < indices.size(); ++i)
3758 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3760 make_array_view(dof_values.begin(), dof_values.end()),
3761 this->finite_element_output.shape_hessians,
3762 *fe,
3763 this->finite_element_output.shape_function_to_row_table,
3764 laplacians,
3765 quadrature_points_fastest,
3766 indices.size() / dofs_per_cell);
3767}
3768
3769
3770
3771template <int dim, int spacedim>
3772template <class InputVector>
3773void
3775 const InputVector &fe_function,
3777 &third_derivatives) const
3778{
3779 using Number = typename InputVector::value_type;
3780 AssertDimension(fe->n_components(), 1);
3781 Assert(this->update_flags & update_3rd_derivatives,
3782 ExcAccessToUninitializedField("update_3rd_derivatives"));
3783 Assert(present_cell.is_initialized(), ExcNotReinited());
3784 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3785
3786 // get function values of dofs on this cell
3787 Vector<Number> dof_values(dofs_per_cell);
3788 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3790 make_array_view(dof_values.begin(), dof_values.end()),
3791 this->finite_element_output.shape_3rd_derivatives,
3792 third_derivatives);
3793}
3794
3795
3796
3797template <int dim, int spacedim>
3798template <class InputVector>
3799void
3801 const InputVector & fe_function,
3804 &third_derivatives) const
3805{
3806 using Number = typename InputVector::value_type;
3807 Assert(this->update_flags & update_3rd_derivatives,
3808 ExcAccessToUninitializedField("update_3rd_derivatives"));
3809 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3810 AssertDimension(indices.size(), dofs_per_cell);
3811
3812 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3813 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3814 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3816 make_array_view(dof_values.begin(), dof_values.end()),
3817 this->finite_element_output.shape_3rd_derivatives,
3818 third_derivatives);
3819}
3820
3821
3822
3823template <int dim, int spacedim>
3824template <class InputVector>
3825void
3827 const InputVector &fe_function,
3828 std::vector<
3830 & third_derivatives,
3831 const bool quadrature_points_fastest) const
3832{
3833 using Number = typename InputVector::value_type;
3834 Assert(this->update_flags & update_3rd_derivatives,
3835 ExcAccessToUninitializedField("update_3rd_derivatives"));
3836 Assert(present_cell.is_initialized(), ExcNotReinited());
3837 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3838
3839 // get function values of dofs on this cell
3840 Vector<Number> dof_values(dofs_per_cell);
3841 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3843 make_array_view(dof_values.begin(), dof_values.end()),
3844 this->finite_element_output.shape_3rd_derivatives,
3845 *fe,
3846 this->finite_element_output.shape_function_to_row_table,
3847 make_array_view(third_derivatives.begin(), third_derivatives.end()),
3848 quadrature_points_fastest);
3849}
3850
3851
3852
3853template <int dim, int spacedim>
3854template <class InputVector>
3855void
3857 const InputVector & fe_function,
3860 third_derivatives,
3861 const bool quadrature_points_fastest) const
3862{
3863 using Number = typename InputVector::value_type;
3864 Assert(this->update_flags & update_3rd_derivatives,
3865 ExcAccessToUninitializedField("update_3rd_derivatives"));
3866 Assert(indices.size() % dofs_per_cell == 0,
3867 ExcNotMultiple(indices.size(), dofs_per_cell));
3868
3869 boost::container::small_vector<Number, 200> dof_values(indices.size());
3870 for (unsigned int i = 0; i < indices.size(); ++i)
3871 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3873 make_array_view(dof_values.begin(), dof_values.end()),
3874 this->finite_element_output.shape_3rd_derivatives,
3875 *fe,
3876 this->finite_element_output.shape_function_to_row_table,
3877 make_array_view(third_derivatives.begin(), third_derivatives.end()),
3878 quadrature_points_fastest,
3879 indices.size() / dofs_per_cell);
3880}
3881
3882
3883
3884template <int dim, int spacedim>
3888 return present_cell;
3889}
3890
3891
3892
3893template <int dim, int spacedim>
3894const std::vector<Tensor<1, spacedim>> &
3897 Assert(this->update_flags & update_normal_vectors,
3899 "update_normal_vectors")));
3900
3901 return this->mapping_output.normal_vectors;
3902}
3903
3904
3905
3906template <int dim, int spacedim>
3907std::size_t
3909{
3910 return (sizeof(this->update_flags) +
3911 MemoryConsumption::memory_consumption(n_quadrature_points) +
3912 MemoryConsumption::memory_consumption(max_n_quadrature_points) +
3913 sizeof(cell_similarity) +
3922 MemoryConsumption::memory_consumption(finite_element_output));
3923}
3924
3925
3926
3927template <int dim, int spacedim>
3930 const UpdateFlags update_flags) const
3931{
3932 // first find out which objects need to be recomputed on each
3933 // cell we visit. this we have to ask the finite element and mapping.
3934 // elements are first since they might require update in mapping
3935 //
3936 // there is no need to iterate since mappings will never require
3937 // the finite element to compute something for them
3938 UpdateFlags flags = update_flags | fe->requires_update_flags(update_flags);
3939 flags |= mapping->requires_update_flags(flags);
3940
3941 return flags;
3942}
3943
3944
3945
3946template <int dim, int spacedim>
3947void
3950 // if there is no present cell, then we shouldn't be
3951 // connected via a signal to a triangulation
3952 Assert(present_cell.is_initialized(), ExcInternalError());
3953
3954 // so delete the present cell and
3955 // disconnect from the signal we have with
3956 // it
3957 tria_listener_refinement.disconnect();
3958 tria_listener_mesh_transform.disconnect();
3959 present_cell = {};
3960}
3961
3962
3963
3964template <int dim, int spacedim>
3965void
3968{
3969 if (present_cell.is_initialized())
3970 {
3971 if (&cell->get_triangulation() !=
3972 &present_cell
3973 .
3975 ->get_triangulation())
3976 {
3977 // the triangulations for the previous cell and the current cell
3978 // do not match. disconnect from the previous triangulation and
3979 // connect to the current one; also invalidate the previous
3980 // cell because we shouldn't be comparing cells from different
3981 // triangulations
3982 invalidate_present_cell();
3983 tria_listener_refinement =
3984 cell->get_triangulation().signals.any_change.connect(
3985 [this]() { this->invalidate_present_cell(); });
3986 tria_listener_mesh_transform =
3987 cell->get_triangulation().signals.mesh_movement.connect(
3988 [this]() { this->invalidate_present_cell(); });
3989 }
3990 }
3991 else
3992 {
3993 // if this FEValues has never been set to any cell at all, then
3994 // at least subscribe to the triangulation to get notified of
3995 // changes
3996 tria_listener_refinement =
3997 cell->get_triangulation().signals.post_refinement.connect(
3998 [this]() { this->invalidate_present_cell(); });
3999 tria_listener_mesh_transform =
4000 cell->get_triangulation().signals.mesh_movement.connect(
4001 [this]() { this->invalidate_present_cell(); });
4002 }
4003}
4004
4005
4006
4007template <int dim, int spacedim>
4008inline void
4011{
4012 // Unfortunately, the detection of simple geometries with CellSimilarity is
4013 // sensitive to the first cell detected. When doing this with multiple
4014 // threads, each thread will get its own scratch data object with an
4015 // FEValues object in the implementation framework from late 2013, which is
4016 // initialized to the first cell the thread sees. As this number might
4017 // different between different runs (after all, the tasks are scheduled
4018 // dynamically onto threads), this slight deviation leads to difference in
4019 // roundoff errors that propagate through the program. Therefore, we need to
4020 // disable CellSimilarity in case there is more than one thread in the
4021 // problem. This will likely not affect many MPI test cases as there
4022 // multithreading is disabled on default, but in many other situations
4023 // because we rarely explicitly set the number of threads.
4024 //
4025 // TODO: Is it reasonable to introduce a flag "unsafe" in the constructor of
4026 // FEValues to re-enable this feature?
4028 {
4029 cell_similarity = CellSimilarity::none;
4030 return;
4031 }
4032
4033 // case that there has not been any cell before
4034 if (this->present_cell.is_initialized() == false)
4035 cell_similarity = CellSimilarity::none;
4036 else
4037 // in MappingQ, data can have been modified during the last call. Then, we
4038 // can't use that data on the new cell.
4039 if (cell_similarity == CellSimilarity::invalid_next_cell)
4040 cell_similarity = CellSimilarity::none;
4041 else
4042 cell_similarity =
4043 (cell->is_translation_of(
4044 static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
4045 &>(this->present_cell)) ?
4048
4049 if ((dim < spacedim) && (cell_similarity == CellSimilarity::translation))
4050 {
4051 if (static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
4052 &>(this->present_cell)
4053 ->direction_flag() != cell->direction_flag())
4054 cell_similarity = CellSimilarity::inverted_translation;
4055 }
4056 // TODO: here, one could implement other checks for similarity, e.g. for
4057 // children of a parallelogram.
4058}
4059
4060
4061
4062template <int dim, int spacedim>
4065{
4066 return cell_similarity;
4067}
4068
4069
4070
4071template <int dim, int spacedim>
4073
4074
4075
4076template <int dim, int spacedim>
4078
4079/*------------------------------- FEValues -------------------------------*/
4080
4081template <int dim, int spacedim>
4083
4084
4085
4086template <int dim, int spacedim>
4089 const Quadrature<dim> & q,
4090 const UpdateFlags update_flags)
4091 : FEValuesBase<dim, spacedim>(q.size(),
4092 fe.n_dofs_per_cell(),
4094 mapping,
4095 fe)
4096 , quadrature(q)
4097{
4098 initialize(update_flags);
4099}
4100
4101
4102
4103template <int dim, int spacedim>
4106 const hp::QCollection<dim> & q,
4107 const UpdateFlags update_flags)
4108 : FEValues(mapping, fe, q[0], update_flags)
4109{
4110 AssertDimension(q.size(), 1);
4111}
4112
4113
4114
4115template <int dim, int spacedim>
4117 const Quadrature<dim> & q,
4118 const UpdateFlags update_flags)
4119 : FEValuesBase<dim, spacedim>(
4120 q.size(),
4121 fe.n_dofs_per_cell(),
4123 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4124 fe)
4125 , quadrature(q)
4126{
4127 initialize(update_flags);
4128}
4129
4130
4131
4132template <int dim, int spacedim>
4134 const hp::QCollection<dim> & q,
4135 const UpdateFlags update_flags)
4136 : FEValues(fe, q[0], update_flags)
4137{
4138 AssertDimension(q.size(), 1);
4139}
4140
4141
4142
4143template <int dim, int spacedim>
4144void
4146{
4147 // You can compute normal vectors to the cells only in the
4148 // codimension one case.
4149 if (dim != spacedim - 1)
4150 Assert((update_flags & update_normal_vectors) == false,
4151 ExcMessage("You can only pass the 'update_normal_vectors' "
4152 "flag to FEFaceValues or FESubfaceValues objects, "
4153 "but not to an FEValues object unless the "
4154 "triangulation it refers to is embedded in a higher "
4155 "dimensional space."));
4156
4157 const UpdateFlags flags = this->compute_update_flags(update_flags);
4158
4159 // initialize the base classes
4160 if (flags & update_mapping)
4161 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4162 this->finite_element_output.initialize(this->max_n_quadrature_points,
4163 *this->fe,
4164 flags);
4165
4166 // then get objects into which the FE and the Mapping can store
4167 // intermediate data used across calls to reinit. we can do this in parallel
4169 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4170 fe_get_data = Threads::new_task([&]() {
4171 return this->fe->get_data(flags,
4172 *this->mapping,
4173 quadrature,
4174 this->finite_element_output);
4175 });
4176
4178 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4179 mapping_get_data;
4180 if (flags & update_mapping)
4181 mapping_get_data = Threads::new_task(
4182 [&]() { return this->mapping->get_data(flags, quadrature); });
4183
4184 this->update_flags = flags;
4185
4186 // then collect answers from the two task above
4187 this->fe_data = std::move(fe_get_data.return_value());
4188 if (flags & update_mapping)
4189 this->mapping_data = std::move(mapping_get_data.return_value());
4190 else
4191 this->mapping_data =
4192 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4193}
4194
4195
4196
4197template <int dim, int spacedim>
4198void
4201{
4202 // Check that mapping and reference cell type are compatible:
4203 Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
4204 ExcMessage(
4205 "You are trying to call FEValues::reinit() with a cell of type " +
4206 cell->reference_cell().to_string() +
4207 " with a Mapping that is not compatible with it."));
4208
4209 // no FE in this cell, so no assertion
4210 // necessary here
4211 this->maybe_invalidate_previous_present_cell(cell);
4212 this->check_cell_similarity(cell);
4213
4214 this->present_cell = {cell};
4215
4216 // this was the part of the work that is dependent on the actual
4217 // data type of the iterator. now pass on to the function doing
4218 // the real work.
4219 do_reinit();
4220}
4221
4222
4223
4224template <int dim, int spacedim>
4225template <bool lda>
4226void
4229{
4230 // assert that the finite elements passed to the constructor and
4231 // used by the DoFHandler used by this cell, are the same
4232 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4233 static_cast<const FiniteElementData<dim> &>(cell->get_fe()),
4235
4236 // Check that mapping and reference cell type are compatible:
4237 Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
4238 ExcMessage(
4239 "You are trying to call FEValues::reinit() with a cell of type " +
4240 cell->reference_cell().to_string() +
4241 " with a Mapping that is not compatible with it."));
4242
4243 this->maybe_invalidate_previous_present_cell(cell);
4244 this->check_cell_similarity(cell);
4245
4246 this->present_cell = {cell};
4247
4248 // this was the part of the work that is dependent on the actual
4249 // data type of the iterator. now pass on to the function doing
4250 // the real work.
4251 do_reinit();
4252}
4253
4254
4255
4256template <int dim, int spacedim>
4257void
4259{
4260 // first call the mapping and let it generate the data
4261 // specific to the mapping. also let it inspect the
4262 // cell similarity flag and, if necessary, update
4263 // it
4264 if (this->update_flags & update_mapping)
4265 {
4266 this->cell_similarity =
4267 this->get_mapping().fill_fe_values(this->present_cell,
4268 this->cell_similarity,
4269 quadrature,
4270 *this->mapping_data,
4271 this->mapping_output);
4272 }
4273
4274 // then call the finite element and, with the data
4275 // already filled by the mapping, let it compute the
4276 // data for the mapped shape function values, gradients,
4277 // etc.
4278 this->get_fe().fill_fe_values(this->present_cell,
4279 this->cell_similarity,
4280 this->quadrature,
4281 this->get_mapping(),
4282 *this->mapping_data,
4283 this->mapping_output,
4284 *this->fe_data,
4285 this->finite_element_output);
4286}
4287
4288
4289
4290template <int dim, int spacedim>
4291std::size_t
4293{
4296}
4297
4298
4299/*------------------------------- FEFaceValuesBase --------------------------*/
4300
4301
4302template <int dim, int spacedim>
4304 const unsigned int dofs_per_cell,
4305 const UpdateFlags flags,
4306 const Mapping<dim, spacedim> & mapping,
4308 const Quadrature<dim - 1> & quadrature)
4309 : FEFaceValuesBase<dim, spacedim>(dofs_per_cell,
4310 flags,
4311 mapping,
4312 fe,
4313 hp::QCollection<dim - 1>(quadrature))
4314{}
4315
4316
4317
4318template <int dim, int spacedim>
4320 const unsigned int dofs_per_cell,
4321 const UpdateFlags,
4322 const Mapping<dim, spacedim> & mapping,
4324 const hp::QCollection<dim - 1> & quadrature)
4325 : FEValuesBase<dim, spacedim>(quadrature.max_n_quadrature_points(),
4326 dofs_per_cell,
4328 mapping,
4329 fe)
4330 , present_face_index(numbers::invalid_unsigned_int)
4331 , quadrature(quadrature)
4332{
4333 Assert(quadrature.size() == 1 ||
4334 quadrature.size() == fe.reference_cell().n_faces(),
4336}
4337
4338
4339
4340template <int dim, int spacedim>
4341const std::vector<Tensor<1, spacedim>> &
4343{
4344 Assert(this->update_flags & update_boundary_forms,
4346 "update_boundary_forms")));
4347 return this->mapping_output.boundary_forms;
4348}
4349
4350
4351
4352template <int dim, int spacedim>
4353std::size_t
4355{
4358}
4359
4360
4361/*------------------------------- FEFaceValues -------------------------------*/
4362
4363template <int dim, int spacedim>
4365
4366
4367
4368template <int dim, int spacedim>
4370
4371
4372
4373template <int dim, int spacedim>
4375 const Mapping<dim, spacedim> & mapping,
4377 const Quadrature<dim - 1> & quadrature,
4378 const UpdateFlags update_flags)
4379 : FEFaceValues<dim, spacedim>(mapping,
4380 fe,
4381 hp::QCollection<dim - 1>(quadrature),
4382 update_flags)
4383{}
4384
4385
4386
4387template <int dim, int spacedim>
4389 const Mapping<dim, spacedim> & mapping,
4391 const hp::QCollection<dim - 1> & quadrature,
4392 const UpdateFlags update_flags)
4393 : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
4394 update_flags,
4395 mapping,
4396 fe,
4397 quadrature)
4398{
4399 initialize(update_flags);
4400}
4401
4402
4403
4404template <int dim, int spacedim>
4407 const Quadrature<dim - 1> & quadrature,
4408 const UpdateFlags update_flags)
4409 : FEFaceValues<dim, spacedim>(fe,
4410 hp::QCollection<dim - 1>(quadrature),
4411 update_flags)
4412{}
4413
4414
4415
4416template <int dim, int spacedim>
4419 const hp::QCollection<dim - 1> & quadrature,
4420 const UpdateFlags update_flags)
4421 : FEFaceValuesBase<dim, spacedim>(
4422 fe.n_dofs_per_cell(),
4423 update_flags,
4424 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4425 fe,
4426 quadrature)
4427{
4428 initialize(update_flags);
4429}
4430
4431
4432
4433template <int dim, int spacedim>
4434void
4436{
4437 const UpdateFlags flags = this->compute_update_flags(update_flags);
4438
4439 // initialize the base classes
4440 if (flags & update_mapping)
4441 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4442 this->finite_element_output.initialize(this->max_n_quadrature_points,
4443 *this->fe,
4444 flags);
4445
4446 // then get objects into which the FE and the Mapping can store
4447 // intermediate data used across calls to reinit. this can be done in parallel
4448
4449 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> (
4450 FiniteElement<dim, spacedim>::*finite_element_get_face_data)(
4451 const UpdateFlags,
4452 const Mapping<dim, spacedim> &,
4455 spacedim>
4457
4458 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> (
4459 Mapping<dim, spacedim>::*mapping_get_face_data)(
4460 const UpdateFlags, const hp::QCollection<dim - 1> &) const =
4462
4463
4465 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4466 fe_get_data = Threads::new_task(finite_element_get_face_data,
4467 *this->fe,
4468 flags,
4469 *this->mapping,
4470 this->quadrature,
4471 this->finite_element_output);
4473 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4474 mapping_get_data;
4475 if (flags & update_mapping)
4476 mapping_get_data = Threads::new_task(mapping_get_face_data,
4477 *this->mapping,
4478 flags,
4479 this->quadrature);
4480
4481 this->update_flags = flags;
4482
4483 // then collect answers from the two task above
4484 this->fe_data = std::move(fe_get_data.return_value());
4485 if (flags & update_mapping)
4486 this->mapping_data = std::move(mapping_get_data.return_value());
4487 else
4488 this->mapping_data =
4489 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4490}
4491
4492
4493
4494template <int dim, int spacedim>
4495template <bool lda>
4496void
4499 const unsigned int face_no)
4500{
4501 // assert that the finite elements passed to the constructor and
4502 // used by the DoFHandler used by this cell, are the same
4503 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4504 static_cast<const FiniteElementData<dim> &>(
4505 cell->get_dof_handler().get_fe(cell->active_fe_index())),
4507
4509
4510 this->maybe_invalidate_previous_present_cell(cell);
4511 this->present_cell = {cell};
4512
4513 // this was the part of the work that is dependent on the actual
4514 // data type of the iterator. now pass on to the function doing
4515 // the real work.
4516 do_reinit(face_no);
4517}
4518
4519
4520
4521template <int dim, int spacedim>
4522template <bool lda>
4523void
4527{
4528 const auto face_n = cell->face_iterator_to_index(face);
4529 reinit(cell, face_n);
4530}
4531
4532
4533
4534template <int dim, int spacedim>
4535void
4538 const unsigned int face_no)
4539{
4541
4542 this->maybe_invalidate_previous_present_cell(cell);
4543 this->present_cell = {cell};
4544
4545 // this was the part of the work that is dependent on the actual
4546 // data type of the iterator. now pass on to the function doing
4547 // the real work.
4548 do_reinit(face_no);
4549}
4550
4551
4552
4553template <int dim, int spacedim>
4554void
4558{
4559 const auto face_n = cell->face_iterator_to_index(face);
4560 reinit(cell, face_n);
4561}
4562
4563
4564
4565template <int dim, int spacedim>
4566void
4567FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no)
4568{
4569 this->present_face_no = face_no;
4570
4571 // first of all, set the present_face_index (if available)
4573 this->present_cell;
4574 this->present_face_index = cell->face_index(face_no);
4575
4576 if (this->update_flags & update_mapping)
4577 {
4578 this->get_mapping().fill_fe_face_values(this->present_cell,
4579 face_no,
4580 this->quadrature,
4581 *this->mapping_data,
4582 this->mapping_output);
4583 }
4584
4585 this->get_fe().fill_fe_face_values(this->present_cell,
4586 face_no,
4587 this->quadrature,
4588 this->get_mapping(),
4589 *this->mapping_data,
4590 this->mapping_output,
4591 *this->fe_data,
4592 this->finite_element_output);
4593
4594 const_cast<unsigned int &>(this->n_quadrature_points) =
4595 this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size();
4596}
4597
4598
4599/* ---------------------------- FESubFaceValues ---------------------------- */
4600
4601
4602template <int dim, int spacedim>
4604
4605
4606
4607template <int dim, int spacedim>
4609
4610
4611
4612template <int dim, int spacedim>
4614 const Mapping<dim, spacedim> & mapping,
4616 const Quadrature<dim - 1> & quadrature,
4617 const UpdateFlags update_flags)
4618 : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
4619 update_flags,
4620 mapping,
4621 fe,
4622 quadrature)
4623{
4624 initialize(update_flags);
4625}
4626
4627
4628
4629template <int dim, int spacedim>
4631 const Mapping<dim, spacedim> & mapping,
4633 const hp::QCollection<dim - 1> & quadrature,
4634 const UpdateFlags update_flags)
4635 : FESubfaceValues(mapping, fe, quadrature[0], update_flags)
4636{
4637 AssertDimension(quadrature.size(), 1);
4638}
4639
4640
4641
4642template <int dim, int spacedim>
4645 const Quadrature<dim - 1> & quadrature,
4646 const UpdateFlags update_flags)
4647 : FEFaceValuesBase<dim, spacedim>(
4648 fe.n_dofs_per_cell(),
4649 update_flags,
4650 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4651 fe,
4652 quadrature)
4653{
4654 initialize(update_flags);
4655}
4656
4657
4658
4659template <int dim, int spacedim>
4662 const hp::QCollection<dim - 1> & quadrature,
4663 const UpdateFlags update_flags)
4664 : FESubfaceValues(fe, quadrature[0], update_flags)
4665{
4666 AssertDimension(quadrature.size(), 1);
4667}
4668
4669
4670
4671template <int dim, int spacedim>
4672void
4674{
4675 const UpdateFlags flags = this->compute_update_flags(update_flags);
4676
4677 // initialize the base classes
4678 if (flags & update_mapping)
4679 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4680 this->finite_element_output.initialize(this->max_n_quadrature_points,
4681 *this->fe,
4682 flags);
4683
4684 // then get objects into which the FE and the Mapping can store
4685 // intermediate data used across calls to reinit. this can be done
4686 // in parallel
4688 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4689 fe_get_data =
4691 *this->fe,
4692 flags,
4693 *this->mapping,
4694 this->quadrature[0],
4695 this->finite_element_output);
4697 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4698 mapping_get_data;
4699 if (flags & update_mapping)
4700 mapping_get_data =
4702 *this->mapping,
4703 flags,
4704 this->quadrature[0]);
4705
4706 this->update_flags = flags;
4707
4708 // then collect answers from the two task above
4709 this->fe_data = std::move(fe_get_data.return_value());
4710 if (flags & update_mapping)
4711 this->mapping_data = std::move(mapping_get_data.return_value());
4712 else
4713 this->mapping_data =
4714 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4715}
4716
4717
4718
4719template <int dim, int spacedim>
4720template <bool lda>
4721void
4724 const unsigned int face_no,
4725 const unsigned int subface_no)
4726{
4727 // assert that the finite elements passed to the constructor and
4728 // used by the DoFHandler used by this cell, are the same
4729 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4730 static_cast<const FiniteElementData<dim> &>(
4731 cell->get_dof_handler().get_fe(cell->active_fe_index())),
4734 // We would like to check for subface_no < cell->face(face_no)->n_children(),
4735 // but unfortunately the current function is also called for
4736 // faces without children (see tests/fe/mapping.cc). Therefore,
4737 // we must use following workaround of two separate assertions
4738 Assert(cell->face(face_no)->has_children() ||
4740 ExcIndexRange(subface_no,
4741 0,
4743 Assert(!cell->face(face_no)->has_children() ||
4744 subface_no < cell->face(face_no)->n_active_descendants(),
4745 ExcIndexRange(subface_no,
4746 0,
4747 cell->face(face_no)->n_active_descendants()));
4748 Assert(cell->has_children() == false,
4749 ExcMessage("You can't use subface data for cells that are "
4750 "already refined. Iterate over their children "
4751 "instead in these cases."));
4752
4753 this->maybe_invalidate_previous_present_cell(cell);
4754 this->present_cell = {cell};
4755
4756 // this was the part of the work that is dependent on the actual
4757 // data type of the iterator. now pass on to the function doing
4758 // the real work.
4759 do_reinit(face_no, subface_no);
4760}
4761
4762
4763
4764template <int dim, int spacedim>
4765template <bool lda>
4766void
4770 const typename Triangulation<dim, spacedim>::face_iterator &subface)
4771{
4772 reinit(cell,
4773 cell->face_iterator_to_index(face),
4774 face->child_iterator_to_index(subface));
4775}
4776
4777
4778
4779template <int dim, int spacedim>
4780void
4783 const unsigned int face_no,
4784 const unsigned int subface_no)
4785{
4787 // We would like to check for subface_no < cell->face(face_no)->n_children(),
4788 // but unfortunately the current function is also called for
4789 // faces without children for periodic faces, which have hanging nodes on
4790 // the other side (see include/deal.II/matrix_free/mapping_info.templates.h).
4791 AssertIndexRange(subface_no,
4792 (cell->has_periodic_neighbor(face_no) ?
4793 cell->periodic_neighbor(face_no)
4794 ->face(cell->periodic_neighbor_face_no(face_no))
4795 ->n_children() :
4796 cell->face(face_no)->n_children()));
4797
4798 this->maybe_invalidate_previous_present_cell(cell);
4799 this->present_cell = {cell};
4800
4801 // this was the part of the work that is dependent on the actual
4802 // data type of the iterator. now pass on to the function doing
4803 // the real work.
4804 do_reinit(face_no, subface_no);
4805}
4806
4807
4808
4809template <int dim, int spacedim>
4810void
4814 const typename Triangulation<dim, spacedim>::face_iterator &subface)
4815{
4816 reinit(cell,
4817 cell->face_iterator_to_index(face),
4818 face->child_iterator_to_index(subface));
4819}
4820
4821
4822
4823template <int dim, int spacedim>
4824void
4825FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no,
4826 const unsigned int subface_no)
4827{
4828 this->present_face_no = face_no;
4829
4830 // first of all, set the present_face_index (if available)
4832 this->present_cell;
4833
4834 if (!cell->face(face_no)->has_children())
4835 // no subfaces at all, so set present_face_index to this face rather
4836 // than any subface
4837 this->present_face_index = cell->face_index(face_no);
4838 else if (dim != 3)
4839 this->present_face_index = cell->face(face_no)->child_index(subface_no);
4840 else
4841 {
4842 // this is the same logic we use in cell->neighbor_child_on_subface(). See
4843 // there for an explanation of the different cases
4844 unsigned int subface_index = numbers::invalid_unsigned_int;
4845 switch (cell->subface_case(face_no))
4846 {
4850 subface_index = cell->face(face_no)->child_index(subface_no);
4851 break;
4854 subface_index = cell->face(face_no)
4855 ->child(subface_no / 2)
4856 ->child_index(subface_no % 2);
4857 break;
4860 switch (subface_no)
4861 {
4862 case 0:
4863 case 1:
4864 subface_index =
4865 cell->face(face_no)->child(0)->child_index(subface_no);
4866 break;
4867 case 2:
4868 subface_index = cell->face(face_no)->child_index(1);
4869 break;
4870 default:
4871 Assert(false, ExcInternalError());
4872 }
4873 break;
4876 switch (subface_no)
4877 {
4878 case 0:
4879 subface_index = cell->face(face_no)->child_index(0);
4880 break;
4881 case 1:
4882 case 2:
4883 subface_index =
4884 cell->face(face_no)->child(1)->child_index(subface_no - 1);
4885 break;
4886 default:
4887 Assert(false, ExcInternalError());
4888 }
4889 break;
4890 default:
4891 Assert(false, ExcInternalError());
4892 break;
4893 }
4894 Assert(subface_index != numbers::invalid_unsigned_int,
4896 this->present_face_index = subface_index;
4897 }
4898
4899 // now ask the mapping and the finite element to do the actual work
4900 if (this->update_flags & update_mapping)
4901 {
4902 this->get_mapping().fill_fe_subface_values(this->present_cell,
4903 face_no,
4904 subface_no,
4905 this->quadrature[0],
4906 *this->mapping_data,
4907 this->mapping_output);
4908 }
4909
4910 this->get_fe().fill_fe_subface_values(this->present_cell,
4911 face_no,
4912 subface_no,
4913 this->quadrature[0],
4914 this->get_mapping(),
4915 *this->mapping_data,
4916 this->mapping_output,
4917 *this->fe_data,
4918 this->finite_element_output);
4919}
4920
4921
4922/*------------------------------- Explicit Instantiations -------------*/
4923#define SPLIT_INSTANTIATIONS_COUNT 6
4924#ifndef SPLIT_INSTANTIATIONS_INDEX
4925# define SPLIT_INSTANTIATIONS_INDEX 0
4926#endif
4927#include "fe_values.inst"
4928
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:704
iterator begin() const
Definition array_view.h:594
iterator end() const
Definition array_view.h:603
std::size_t size() const
Definition array_view.h:576
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
void get_interpolated_dof_values(const InputVector &values, Vector< number > &interpolated_values, const types::fe_index fe_index=numbers::invalid_fe_index) const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
FEFaceValuesBase(const unsigned int dofs_per_cell, const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature)
std::size_t memory_consumption() const
const std::vector< Tensor< 1, spacedim > > & get_boundary_forms() const
void initialize(const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no)
void do_reinit(const unsigned int face_no)
FEFaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature, const UpdateFlags update_flags)
FESubfaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &face_quadrature, const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no, const unsigned int subface_no)
void initialize(const UpdateFlags update_flags)
void do_reinit(const unsigned int face_no, const unsigned int subface_no)
void get_interpolated_dof_values(const VectorType &in, Vector< typename VectorType::value_type > &out) const
types::global_dof_index n_dofs_for_dof_handler() const
CellSimilarity::Similarity cell_similarity
Definition fe_values.h:4019
CellIteratorContainer present_cell
Definition fe_values.h:3911
::internal::FEValuesViews::Cache< dim, spacedim > fe_values_views_cache
Definition fe_values.h:4034
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &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)
virtual ~FEValuesBase() override
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition fe_values.h:3956
const unsigned int dofs_per_cell
Definition fe_values.h:2451
void check_cell_similarity(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
UpdateFlags update_flags
Definition fe_values.h:4001
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
Definition fe_values.h:3979
const unsigned int n_quadrature_points
Definition fe_values.h:2433
CellSimilarity::Similarity get_cell_similarity() const
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
std::size_t memory_consumption() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
const Triangulation< dim, spacedim >::cell_iterator get_cell() const
void get_function_laplacians(const InputVector &fe_function, std::vector< typename InputVector::value_type > &laplacians) const
UpdateFlags compute_update_flags(const UpdateFlags update_flags) const
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
void invalidate_present_cell()
::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > finite_element_output
Definition fe_values.h:3995
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type > > &hessians) const
const FiniteElement< dim, spacedim > & get_fe() const
void maybe_invalidate_previous_present_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
const unsigned int max_n_quadrature_points
Definition fe_values.h:2444
void get_function_third_derivatives(const InputVector &fe_function, std::vector< Tensor< 3, spacedim, typename InputVector::value_type > > &third_derivatives) const
void get_function_third_derivatives_from_local_dof_values(const InputVector &dof_values, std::vector< solution_third_derivative_type< typename InputVector::value_type > > &third_derivatives) const
typename ProductType< Number, hessian_type >::type solution_hessian_type
Definition fe_values.h:215
const unsigned int component
Definition fe_values.h:635
void get_function_hessians_from_local_dof_values(const InputVector &dof_values, std::vector< solution_hessian_type< typename InputVector::value_type > > &hessians) const
void get_function_values_from_local_dof_values(const InputVector &dof_values, std::vector< solution_value_type< typename InputVector::value_type > > &values) const
std::vector< ShapeFunctionData > shape_function_data
Definition fe_values.h:640
void get_function_laplacians(const InputVector &fe_function, std::vector< solution_laplacian_type< typename InputVector::value_type > > &laplacians) const
typename ProductType< Number, value_type >::type solution_laplacian_type
Definition fe_values.h:205
void get_function_third_derivatives(const InputVector &fe_function, std::vector< solution_third_derivative_type< typename InputVector::value_type > > &third_derivatives) const
void get_function_hessians(const InputVector &fe_function, std::vector< solution_hessian_type< typename InputVector::value_type > > &hessians) const
typename ProductType< Number, third_derivative_type >::type solution_third_derivative_type
Definition fe_values.h:225
typename ProductType< Number, value_type >::type solution_value_type
Definition fe_values.h:185
typename ProductType< Number, gradient_type >::type solution_gradient_type
Definition fe_values.h:195
void get_function_laplacians_from_local_dof_values(const InputVector &dof_values, std::vector< solution_laplacian_type< typename InputVector::value_type > > &laplacians) const
void get_function_gradients(const InputVector &fe_function, std::vector< solution_gradient_type< typename InputVector::value_type > > &gradients) const
void get_function_gradients_from_local_dof_values(const InputVector &dof_values, std::vector< solution_gradient_type< typename InputVector::value_type > > &gradients) const
void get_function_values(const InputVector &fe_function, std::vector< solution_value_type< typename InputVector::value_type > > &values) const
typename ProductType< Number, value_type >::type solution_value_type
Definition fe_values.h:1506
typename ProductType< Number, divergence_type >::type solution_divergence_type
Definition fe_values.h:1516
typename ProductType< Number, value_type >::type solution_value_type
Definition fe_values.h:1842
typename ProductType< Number, divergence_type >::type solution_divergence_type
Definition fe_values.h:1852
typename ProductType< Number, gradient_type >::type solution_gradient_type
Definition fe_values.h:1862
void get_function_hessians(const InputVector &fe_function, std::vector< solution_hessian_type< typename InputVector::value_type > > &hessians) const
typename ProductType< Number, third_derivative_type >::type solution_third_derivative_type
Definition fe_values.h:812
void get_function_symmetric_gradients_from_local_dof_values(const InputVector &dof_values, std::vector< solution_symmetric_gradient_type< typename InputVector::value_type > > &symmetric_gradients) const
void get_function_gradients_from_local_dof_values(const InputVector &dof_values, std::vector< solution_gradient_type< typename InputVector::value_type > > &gradients) const
typename ProductType< Number, divergence_type >::type solution_divergence_type
Definition fe_values.h:773
typename ProductType< Number, hessian_type >::type solution_hessian_type
Definition fe_values.h:802
typename ProductType< Number, symmetric_gradient_type >::type solution_symmetric_gradient_type
Definition fe_values.h:763
void get_function_laplacians_from_local_dof_values(const InputVector &dof_values, std::vector< solution_laplacian_type< typename InputVector::value_type > > &laplacians) const
typename ProductType< Number, gradient_type >::type solution_gradient_type
Definition fe_values.h:753
void get_function_third_derivatives(const InputVector &fe_function, std::vector< solution_third_derivative_type< typename InputVector::value_type > > &third_derivatives) const
typename ProductType< Number, value_type >::type solution_value_type
Definition fe_values.h:743
void get_function_curls_from_local_dof_values(const InputVector &dof_values, std::vector< solution_curl_type< typename InputVector::value_type > > &curls) const
void get_function_values(const InputVector &fe_function, std::vector< solution_value_type< typename InputVector::value_type > > &values) const
const unsigned int first_vector_component
Definition fe_values.h:1442
typename ProductType< Number, curl_type >::type solution_curl_type
Definition fe_values.h:792
std::vector< ShapeFunctionData > shape_function_data
Definition fe_values.h:1447
void get_function_gradients(const InputVector &fe_function, std::vector< solution_gradient_type< typename InputVector::value_type > > &gradients) const
void get_function_symmetric_gradients(const InputVector &fe_function, std::vector< solution_symmetric_gradient_type< typename InputVector::value_type > > &symmetric_gradients) const
void get_function_hessians_from_local_dof_values(const InputVector &dof_values, std::vector< solution_hessian_type< typename InputVector::value_type > > &hessians) const
void get_function_laplacians(const InputVector &fe_function, std::vector< solution_laplacian_type< typename InputVector::value_type > > &laplacians) const
typename ProductType< Number, value_type >::type solution_laplacian_type
Definition fe_values.h:783
void get_function_values_from_local_dof_values(const InputVector &dof_values, std::vector< solution_value_type< typename InputVector::value_type > > &values) const
void get_function_curls(const InputVector &fe_function, std::vector< solution_curl_type< typename InputVector::value_type > > &curls) const
void get_function_divergences_from_local_dof_values(const InputVector &dof_values, std::vector< solution_divergence_type< typename InputVector::value_type > > &divergences) const
void get_function_third_derivatives_from_local_dof_values(const InputVector &dof_values, std::vector< solution_third_derivative_type< typename InputVector::value_type > > &third_derivatives) const
void get_function_divergences(const InputVector &fe_function, std::vector< solution_divergence_type< typename InputVector::value_type > > &divergences) const
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
void do_reinit()
void initialize(const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
std::size_t memory_consumption() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
virtual std::unique_ptr< InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
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
bool is_element(const size_type index) const
Definition index_set.h:1776
signed int value_type
Definition index_set.h:96
Abstract base class for mapping classes.
Definition mapping.h:317
static unsigned int n_threads()
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)
DEAL_II_HOST constexpr SymmetricTensor()=default
friend class Tensor
Definition tensor.h:907
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
iterator end()
iterator begin()
unsigned int size() const
Definition collection.h:265
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotMultiple(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotReinited()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition tria.h:1370
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_mapping
@ update_gradients
Shape function gradients.
@ update_default
No update.
@ update_boundary_forms
Outer normal vector, not normalized.
Task< RT > new_task(const std::function< RT()> &function)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:285
void do_function_laplacians(const ArrayView< 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.cc:522
void do_function_values(const ArrayView< 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.cc:440
void do_function_gradients(const ArrayView< 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_divergences(const ArrayView< 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.cc:761
void do_function_curls(const ArrayView< 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.cc:824
void do_function_symmetric_gradients(const ArrayView< 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.cc:691
void do_function_derivatives(const ArrayView< 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.cc:480
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
Definition hp.h:118
void do_function_laplacians(const ArrayView< Number2 > &dof_values, const ::Table< 2, Tensor< 2, spacedim > > &shape_hessians, std::vector< Number > &laplacians)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
VectorType::value_type get_vector_element(const VectorType &vector, const types::global_dof_index cell_number)
Definition fe_values.cc:60
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:80
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:213
T signaling_nan()
STL namespace.
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
Cache(const FEValuesBase< dim, spacedim > &fe_values)
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 > &)