Reference documentation for deal.II version 9.3.3
\(\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\}}\)
fe_values.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2021 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
22
24
26
27#include <deal.II/fe/fe.h>
30
33
44#include <deal.II/lac/vector.h>
46
48#include <boost/container/small_vector.hpp>
50
51#include <iomanip>
52#include <memory>
53#include <type_traits>
54
56
57
58namespace internal
59{
60 template <class VectorType>
61 typename VectorType::value_type inline get_vector_element(
62 const VectorType & vector,
63 const types::global_dof_index cell_number)
64 {
65 return internal::ElementAccess<VectorType>::get(vector, cell_number);
66 }
67
68
69
71 const IndexSet & is,
72 const types::global_dof_index cell_number)
73 {
74 return (is.is_element(cell_number) ? 1 : 0);
75 }
76
77
78
79 template <int dim, int spacedim>
80 inline std::vector<unsigned int>
82 {
83 std::vector<unsigned int> shape_function_to_row_table(
85 unsigned int row = 0;
86 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
87 {
88 // loop over all components that are nonzero for this particular
89 // shape function. if a component is zero then we leave the
90 // value in the table unchanged (at the invalid value)
91 // otherwise it is mapped to the next free entry
92 unsigned int nth_nonzero_component = 0;
93 for (unsigned int c = 0; c < fe.n_components(); ++c)
94 if (fe.get_nonzero_components(i)[c] == true)
95 {
96 shape_function_to_row_table[i * fe.n_components() + c] =
97 row + nth_nonzero_component;
98 ++nth_nonzero_component;
99 }
100 row += fe.n_nonzero_components(i);
101 }
102
103 return shape_function_to_row_table;
104 }
105
106 namespace
107 {
108 // Check to see if a DoF value is zero, implying that subsequent operations
109 // with the value have no effect.
110 template <typename Number, typename T = void>
111 struct CheckForZero
112 {
113 static bool
114 value(const Number &value)
115 {
116 return value == ::internal::NumberType<Number>::value(0.0);
117 }
118 };
119
120 // For auto-differentiable numbers, the fact that a DoF value is zero
121 // does not imply that its derivatives are zero as well. So we
122 // can't filter by value for these number types.
123 // Note that we also want to avoid actually checking the value itself,
124 // since some AD numbers are not contextually convertible to booleans.
125 template <typename Number>
126 struct CheckForZero<
127 Number,
128 typename std::enable_if<
129 Differentiation::AD::is_ad_number<Number>::value>::type>
130 {
131 static bool
132 value(const Number & /*value*/)
133 {
134 return false;
135 }
136 };
137 } // namespace
138} // namespace internal
139
140
141
142namespace FEValuesViews
143{
144 template <int dim, int spacedim>
146 const unsigned int component)
147 : fe_values(&fe_values)
148 , component(component)
149 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
150 {
151 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
153
154 // TODO: we'd like to use the fields with the same name as these
155 // variables from FEValuesBase, but they aren't initialized yet
156 // at the time we get here, so re-create it all
157 const std::vector<unsigned int> shape_function_to_row_table =
159
160 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
161 {
162 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
163
164 if (is_primitive == true)
165 shape_function_data[i].is_nonzero_shape_function_component =
166 (component == fe.system_to_component_index(i).first);
167 else
168 shape_function_data[i].is_nonzero_shape_function_component =
169 (fe.get_nonzero_components(i)[component] == true);
170
171 if (shape_function_data[i].is_nonzero_shape_function_component == true)
172 shape_function_data[i].row_index =
173 shape_function_to_row_table[i * fe.n_components() + component];
174 else
176 }
177 }
178
179
180
181 template <int dim, int spacedim>
183 : fe_values(nullptr)
184 , component(numbers::invalid_unsigned_int)
185 {}
186
187
188
189 template <int dim, int spacedim>
191 const unsigned int first_vector_component)
192 : fe_values(&fe_values)
193 , first_vector_component(first_vector_component)
194 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
195 {
196 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
198
199 // TODO: we'd like to use the fields with the same name as these
200 // variables from FEValuesBase, but they aren't initialized yet
201 // at the time we get here, so re-create it all
202 const std::vector<unsigned int> shape_function_to_row_table =
204
205 for (unsigned int d = 0; d < spacedim; ++d)
206 {
207 const unsigned int component = first_vector_component + d;
208
209 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
210 {
211 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
212
213 if (is_primitive == true)
214 shape_function_data[i].is_nonzero_shape_function_component[d] =
215 (component == fe.system_to_component_index(i).first);
216 else
217 shape_function_data[i].is_nonzero_shape_function_component[d] =
218 (fe.get_nonzero_components(i)[component] == true);
219
220 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
221 true)
222 shape_function_data[i].row_index[d] =
223 shape_function_to_row_table[i * fe.n_components() + component];
224 else
225 shape_function_data[i].row_index[d] =
227 }
228 }
229
230 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
231 {
232 unsigned int n_nonzero_components = 0;
233 for (unsigned int d = 0; d < spacedim; ++d)
234 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
235 true)
236 ++n_nonzero_components;
237
238 if (n_nonzero_components == 0)
239 shape_function_data[i].single_nonzero_component = -2;
240 else if (n_nonzero_components > 1)
241 shape_function_data[i].single_nonzero_component = -1;
242 else
243 {
244 for (unsigned int d = 0; d < spacedim; ++d)
246 .is_nonzero_shape_function_component[d] == true)
247 {
248 shape_function_data[i].single_nonzero_component =
249 shape_function_data[i].row_index[d];
250 shape_function_data[i].single_nonzero_component_index = d;
251 break;
252 }
253 }
254 }
255 }
256
257
258
259 template <int dim, int spacedim>
261 : fe_values(nullptr)
262 , first_vector_component(numbers::invalid_unsigned_int)
263 {}
264
265
266
267 template <int dim, int spacedim>
269 const FEValuesBase<dim, spacedim> &fe_values,
270 const unsigned int first_tensor_component)
271 : fe_values(&fe_values)
272 , first_tensor_component(first_tensor_component)
273 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
274 {
275 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
276 Assert(first_tensor_component + (dim * dim + dim) / 2 - 1 <
277 fe.n_components(),
279 first_tensor_component +
281 0,
282 fe.n_components()));
283 // TODO: we'd like to use the fields with the same name as these
284 // variables from FEValuesBase, but they aren't initialized yet
285 // at the time we get here, so re-create it all
286 const std::vector<unsigned int> shape_function_to_row_table =
288
289 for (unsigned int d = 0;
290 d < ::SymmetricTensor<2, dim>::n_independent_components;
291 ++d)
292 {
293 const unsigned int component = first_tensor_component + d;
294
295 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
296 {
297 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
298
299 if (is_primitive == true)
300 shape_function_data[i].is_nonzero_shape_function_component[d] =
301 (component == fe.system_to_component_index(i).first);
302 else
303 shape_function_data[i].is_nonzero_shape_function_component[d] =
304 (fe.get_nonzero_components(i)[component] == true);
305
306 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
307 true)
308 shape_function_data[i].row_index[d] =
309 shape_function_to_row_table[i * fe.n_components() + component];
310 else
311 shape_function_data[i].row_index[d] =
314 }
315
316 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
317 {
318 unsigned int n_nonzero_components = 0;
319 for (unsigned int d = 0;
320 d < ::SymmetricTensor<2, dim>::n_independent_components;
321 ++d)
322 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
323 true)
324 ++n_nonzero_components;
325
326 if (n_nonzero_components == 0)
327 shape_function_data[i].single_nonzero_component = -2;
328 else if (n_nonzero_components > 1)
329 shape_function_data[i].single_nonzero_component = -1;
330 else
331 {
332 for (unsigned int d = 0;
333 d < ::SymmetricTensor<2, dim>::n_independent_components;
334 ++d)
335 if (shape_function_data[i]
336 .is_nonzero_shape_function_component[d] == true)
337 {
338 shape_function_data[i].single_nonzero_component =
339 shape_function_data[i].row_index[d];
340 shape_function_data[i].single_nonzero_component_index = d;
341 break;
342 }
343 }
344 }
345 }
346
347
348
349 template <int dim, int spacedim>
351 : fe_values(nullptr)
352 , first_tensor_component(numbers::invalid_unsigned_int)
353 {}
354
355
356
357 template <int dim, int spacedim>
359 const unsigned int first_tensor_component)
360 : fe_values(&fe_values)
361 , first_tensor_component(first_tensor_component)
362 , shape_function_data(this->fe_values->fe->n_dofs_per_cell())
363 {
364 const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
365 AssertIndexRange(first_tensor_component + dim * dim - 1, fe.n_components());
366 // TODO: we'd like to use the fields with the same name as these
367 // variables from FEValuesBase, but they aren't initialized yet
368 // at the time we get here, so re-create it all
369 const std::vector<unsigned int> shape_function_to_row_table =
371
372 for (unsigned int d = 0; d < dim * dim; ++d)
373 {
374 const unsigned int component = first_tensor_component + d;
375
376 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
377 {
378 const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
379
380 if (is_primitive == true)
381 shape_function_data[i].is_nonzero_shape_function_component[d] =
382 (component == fe.system_to_component_index(i).first);
383 else
384 shape_function_data[i].is_nonzero_shape_function_component[d] =
385 (fe.get_nonzero_components(i)[component] == true);
386
387 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
388 true)
389 shape_function_data[i].row_index[d] =
390 shape_function_to_row_table[i * fe.n_components() + component];
391 else
392 shape_function_data[i].row_index[d] =
394 }
395 }
396
397 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
398 {
399 unsigned int n_nonzero_components = 0;
400 for (unsigned int d = 0; d < dim * dim; ++d)
401 if (shape_function_data[i].is_nonzero_shape_function_component[d] ==
402 true)
403 ++n_nonzero_components;
404
405 if (n_nonzero_components == 0)
406 shape_function_data[i].single_nonzero_component = -2;
407 else if (n_nonzero_components > 1)
408 shape_function_data[i].single_nonzero_component = -1;
409 else
410 {
411 for (unsigned int d = 0; d < dim * dim; ++d)
412 if (shape_function_data[i]
413 .is_nonzero_shape_function_component[d] == true)
414 {
415 shape_function_data[i].single_nonzero_component =
416 shape_function_data[i].row_index[d];
417 shape_function_data[i].single_nonzero_component_index = d;
418 break;
419 }
420 }
421 }
422 }
424
425
426 template <int dim, int spacedim>
428 : fe_values(nullptr)
429 , first_tensor_component(numbers::invalid_unsigned_int)
430 {}
431
432
433
434 namespace internal
435 {
436 // Given values of degrees of freedom, evaluate the
437 // values/gradients/... at quadrature points
438
439 // ------------------------- scalar functions --------------------------
440 template <int dim, int spacedim, typename Number>
441 void
443 const ArrayView<Number> &dof_values,
444 const Table<2, double> & shape_values,
445 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
446 &shape_function_data,
447 std::vector<typename ProductType<Number, double>::type> &values)
448 {
449 const unsigned int dofs_per_cell = dof_values.size();
450 const unsigned int n_quadrature_points = values.size();
451
452 std::fill(values.begin(),
453 values.end(),
455
456 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
457 ++shape_function)
458 if (shape_function_data[shape_function]
459 .is_nonzero_shape_function_component)
460 {
461 const Number &value = dof_values[shape_function];
462 // For auto-differentiable numbers, the fact that a DoF value is
463 // zero does not imply that its derivatives are zero as well. So we
464 // can't filter by value for these number types.
465 if (::internal::CheckForZero<Number>::value(value) == true)
466 continue;
467
468 const double *shape_value_ptr =
469 &shape_values(shape_function_data[shape_function].row_index, 0);
470 for (unsigned int q_point = 0; q_point < n_quadrature_points;
471 ++q_point)
472 values[q_point] += value * (*shape_value_ptr++);
473 }
474 }
475
476
477
478 // same code for gradient and Hessian, template argument 'order' to give
479 // the order of the derivative (= rank of gradient/Hessian tensor)
480 template <int order, int dim, int spacedim, typename Number>
481 void
483 const ArrayView<Number> & dof_values,
484 const Table<2, ::Tensor<order, spacedim>> &shape_derivatives,
485 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
486 &shape_function_data,
487 std::vector<
488 typename ProductType<Number, ::Tensor<order, spacedim>>::type>
489 &derivatives)
490 {
491 const unsigned int dofs_per_cell = dof_values.size();
492 const unsigned int n_quadrature_points = derivatives.size();
493
494 std::fill(
495 derivatives.begin(),
496 derivatives.end(),
497 typename ProductType<Number, ::Tensor<order, spacedim>>::type());
498
499 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
500 ++shape_function)
501 if (shape_function_data[shape_function]
502 .is_nonzero_shape_function_component)
503 {
504 const Number &value = dof_values[shape_function];
505 // For auto-differentiable numbers, the fact that a DoF value is
506 // zero does not imply that its derivatives are zero as well. So we
507 // can't filter by value for these number types.
508 if (::internal::CheckForZero<Number>::value(value) == true)
509 continue;
510
511 const ::Tensor<order, spacedim> *shape_derivative_ptr =
512 &shape_derivatives[shape_function_data[shape_function].row_index]
513 [0];
514 for (unsigned int q_point = 0; q_point < n_quadrature_points;
515 ++q_point)
516 derivatives[q_point] += value * (*shape_derivative_ptr++);
517 }
518 }
519
520
521
522 template <int dim, int spacedim, typename Number>
523 void
525 const ArrayView<Number> & dof_values,
526 const Table<2, ::Tensor<2, spacedim>> &shape_hessians,
527 const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
528 &shape_function_data,
529 std::vector<typename Scalar<dim, spacedim>::
530 template solution_laplacian_type<Number>> &laplacians)
531 {
532 const unsigned int dofs_per_cell = dof_values.size();
533 const unsigned int n_quadrature_points = laplacians.size();
534
535 std::fill(
536 laplacians.begin(),
537 laplacians.end(),
538 typename Scalar<dim,
539 spacedim>::template solution_laplacian_type<Number>());
540
541 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
542 ++shape_function)
543 if (shape_function_data[shape_function]
544 .is_nonzero_shape_function_component)
545 {
546 const Number &value = dof_values[shape_function];
547 // For auto-differentiable numbers, the fact that a DoF value is
548 // zero does not imply that its derivatives are zero as well. So we
549 // can't filter by value for these number types.
550 if (::internal::CheckForZero<Number>::value(value) == true)
551 continue;
552
553 const ::Tensor<2, spacedim> *shape_hessian_ptr =
554 &shape_hessians[shape_function_data[shape_function].row_index][0];
555 for (unsigned int q_point = 0; q_point < n_quadrature_points;
556 ++q_point)
557 laplacians[q_point] += value * trace(*shape_hessian_ptr++);
558 }
559 }
560
561
562
563 // ----------------------------- vector part ---------------------------
565 template <int dim, int spacedim, typename Number>
566 void
568 const ArrayView<Number> &dof_values,
569 const Table<2, double> & shape_values,
570 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
571 &shape_function_data,
572 std::vector<
573 typename ProductType<Number, ::Tensor<1, spacedim>>::type>
574 &values)
575 {
576 const unsigned int dofs_per_cell = dof_values.size();
577 const unsigned int n_quadrature_points = values.size();
578
579 std::fill(
580 values.begin(),
581 values.end(),
582 typename ProductType<Number, ::Tensor<1, spacedim>>::type());
583
584 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
585 ++shape_function)
586 {
587 const int snc =
588 shape_function_data[shape_function].single_nonzero_component;
589
590 if (snc == -2)
591 // shape function is zero for the selected components
592 continue;
593
594 const Number &value = dof_values[shape_function];
595 // For auto-differentiable numbers, the fact that a DoF value is zero
596 // does not imply that its derivatives are zero as well. So we
597 // can't filter by value for these number types.
598 if (::internal::CheckForZero<Number>::value(value) == true)
599 continue;
600
601 if (snc != -1)
602 {
603 const unsigned int comp = shape_function_data[shape_function]
604 .single_nonzero_component_index;
605 const double *shape_value_ptr = &shape_values(snc, 0);
606 for (unsigned int q_point = 0; q_point < n_quadrature_points;
607 ++q_point)
608 values[q_point][comp] += value * (*shape_value_ptr++);
609 }
610 else
611 for (unsigned int d = 0; d < spacedim; ++d)
612 if (shape_function_data[shape_function]
613 .is_nonzero_shape_function_component[d])
614 {
615 const double *shape_value_ptr = &shape_values(
616 shape_function_data[shape_function].row_index[d], 0);
617 for (unsigned int q_point = 0; q_point < n_quadrature_points;
618 ++q_point)
619 values[q_point][d] += value * (*shape_value_ptr++);
620 }
621 }
622 }
623
624
625
626 template <int order, int dim, int spacedim, typename Number>
627 void
629 const ArrayView<Number> & dof_values,
630 const Table<2, ::Tensor<order, spacedim>> &shape_derivatives,
631 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
632 &shape_function_data,
633 std::vector<
634 typename ProductType<Number, ::Tensor<order + 1, spacedim>>::type>
635 &derivatives)
636 {
637 const unsigned int dofs_per_cell = dof_values.size();
638 const unsigned int n_quadrature_points = derivatives.size();
639
640 std::fill(
641 derivatives.begin(),
642 derivatives.end(),
643 typename ProductType<Number,
645
646 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
647 ++shape_function)
648 {
649 const int snc =
650 shape_function_data[shape_function].single_nonzero_component;
651
652 if (snc == -2)
653 // shape function is zero for the selected components
654 continue;
655
656 const Number &value = dof_values[shape_function];
657 // For auto-differentiable numbers, the fact that a DoF value is zero
658 // does not imply that its derivatives are zero as well. So we
659 // can't filter by value for these number types.
660 if (::internal::CheckForZero<Number>::value(value) == true)
661 continue;
662
663 if (snc != -1)
664 {
665 const unsigned int comp = shape_function_data[shape_function]
666 .single_nonzero_component_index;
667 const ::Tensor<order, spacedim> *shape_derivative_ptr =
668 &shape_derivatives[snc][0];
669 for (unsigned int q_point = 0; q_point < n_quadrature_points;
670 ++q_point)
671 derivatives[q_point][comp] += value * (*shape_derivative_ptr++);
672 }
673 else
674 for (unsigned int d = 0; d < spacedim; ++d)
675 if (shape_function_data[shape_function]
676 .is_nonzero_shape_function_component[d])
677 {
678 const ::Tensor<order, spacedim> *shape_derivative_ptr =
679 &shape_derivatives[shape_function_data[shape_function]
680 .row_index[d]][0];
681 for (unsigned int q_point = 0; q_point < n_quadrature_points;
682 ++q_point)
683 derivatives[q_point][d] +=
684 value * (*shape_derivative_ptr++);
685 }
686 }
687 }
688
689
690
691 template <int dim, int spacedim, typename Number>
692 void
694 const ArrayView<Number> & dof_values,
695 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
696 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
697 &shape_function_data,
698 std::vector<
699 typename ProductType<Number,
701 &symmetric_gradients)
702 {
703 const unsigned int dofs_per_cell = dof_values.size();
704 const unsigned int n_quadrature_points = symmetric_gradients.size();
705
706 std::fill(
707 symmetric_gradients.begin(),
708 symmetric_gradients.end(),
709 typename ProductType<Number,
711
712 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
713 ++shape_function)
714 {
715 const int snc =
716 shape_function_data[shape_function].single_nonzero_component;
717
718 if (snc == -2)
719 // shape function is zero for the selected components
720 continue;
721
722 const Number &value = dof_values[shape_function];
723 // For auto-differentiable numbers, the fact that a DoF value is zero
724 // does not imply that its derivatives are zero as well. So we
725 // can't filter by value for these number types.
726 if (::internal::CheckForZero<Number>::value(value) == true)
727 continue;
728
729 if (snc != -1)
730 {
731 const unsigned int comp = shape_function_data[shape_function]
732 .single_nonzero_component_index;
733 const ::Tensor<1, spacedim> *shape_gradient_ptr =
734 &shape_gradients[snc][0];
735 for (unsigned int q_point = 0; q_point < n_quadrature_points;
736 ++q_point)
737 symmetric_gradients[q_point] +=
739 symmetrize_single_row(comp, *shape_gradient_ptr++));
740 }
741 else
742 for (unsigned int q_point = 0; q_point < n_quadrature_points;
743 ++q_point)
744 {
746 grad;
747 for (unsigned int d = 0; d < spacedim; ++d)
748 if (shape_function_data[shape_function]
749 .is_nonzero_shape_function_component[d])
750 grad[d] =
751 value *
752 shape_gradients[shape_function_data[shape_function]
753 .row_index[d]][q_point];
754 symmetric_gradients[q_point] += symmetrize(grad);
755 }
756 }
757 }
758
759
760
761 template <int dim, int spacedim, typename Number>
762 void
764 const ArrayView<Number> & dof_values,
765 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
766 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
767 &shape_function_data,
768 std::vector<typename Vector<dim, spacedim>::
769 template solution_divergence_type<Number>> &divergences)
770 {
771 const unsigned int dofs_per_cell = dof_values.size();
772 const unsigned int n_quadrature_points = divergences.size();
773
774 std::fill(
775 divergences.begin(),
776 divergences.end(),
777 typename Vector<dim,
778 spacedim>::template solution_divergence_type<Number>());
779
780 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
781 ++shape_function)
782 {
783 const int snc =
784 shape_function_data[shape_function].single_nonzero_component;
785
786 if (snc == -2)
787 // shape function is zero for the selected components
788 continue;
789
790 const Number &value = dof_values[shape_function];
791 // For auto-differentiable numbers, the fact that a DoF value is zero
792 // does not imply that its derivatives are zero as well. So we
793 // can't filter by value for these number types.
794 if (::internal::CheckForZero<Number>::value(value) == true)
795 continue;
796
797 if (snc != -1)
798 {
799 const unsigned int comp = shape_function_data[shape_function]
800 .single_nonzero_component_index;
801 const ::Tensor<1, spacedim> *shape_gradient_ptr =
802 &shape_gradients[snc][0];
803 for (unsigned int q_point = 0; q_point < n_quadrature_points;
804 ++q_point)
805 divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
806 }
807 else
808 for (unsigned int d = 0; d < spacedim; ++d)
809 if (shape_function_data[shape_function]
810 .is_nonzero_shape_function_component[d])
811 {
812 const ::Tensor<1, spacedim> *shape_gradient_ptr =
813 &shape_gradients[shape_function_data[shape_function]
814 .row_index[d]][0];
815 for (unsigned int q_point = 0; q_point < n_quadrature_points;
816 ++q_point)
817 divergences[q_point] += value * (*shape_gradient_ptr++)[d];
818 }
819 }
820 }
821
822
823
824 template <int dim, int spacedim, typename Number>
825 void
827 const ArrayView<Number> & dof_values,
828 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
829 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
830 &shape_function_data,
831 std::vector<typename ProductType<
832 Number,
833 typename ::internal::CurlType<spacedim>::type>::type> &curls)
834 {
835 const unsigned int dofs_per_cell = dof_values.size();
836 const unsigned int n_quadrature_points = curls.size();
837
838 std::fill(curls.begin(),
839 curls.end(),
840 typename ProductType<
841 Number,
842 typename ::internal::CurlType<spacedim>::type>::type());
843
844 switch (spacedim)
845 {
846 case 1:
847 {
848 Assert(false,
850 "Computing the curl in 1d is not a useful operation"));
851 break;
852 }
853
854 case 2:
855 {
856 for (unsigned int shape_function = 0;
857 shape_function < dofs_per_cell;
858 ++shape_function)
859 {
860 const int snc = shape_function_data[shape_function]
861 .single_nonzero_component;
862
863 if (snc == -2)
864 // shape function is zero for the selected components
865 continue;
866
867 const Number &value = dof_values[shape_function];
868 // For auto-differentiable numbers, the fact that a DoF value
869 // is zero does not imply that its derivatives are zero as
870 // well. So we can't filter by value for these number types.
871 if (::internal::CheckForZero<Number>::value(value) ==
872 true)
873 continue;
874
875 if (snc != -1)
876 {
877 const ::Tensor<1, spacedim> *shape_gradient_ptr =
878 &shape_gradients[snc][0];
879
880 Assert(shape_function_data[shape_function]
881 .single_nonzero_component >= 0,
883 // we're in 2d, so the formula for the curl is simple:
884 if (shape_function_data[shape_function]
885 .single_nonzero_component_index == 0)
886 for (unsigned int q_point = 0;
887 q_point < n_quadrature_points;
888 ++q_point)
889 curls[q_point][0] -=
890 value * (*shape_gradient_ptr++)[1];
891 else
892 for (unsigned int q_point = 0;
893 q_point < n_quadrature_points;
894 ++q_point)
895 curls[q_point][0] +=
896 value * (*shape_gradient_ptr++)[0];
897 }
898 else
899 // we have multiple non-zero components in the shape
900 // functions. not all of them must necessarily be within the
901 // 2-component window this FEValuesViews::Vector object
902 // considers, however.
903 {
904 if (shape_function_data[shape_function]
905 .is_nonzero_shape_function_component[0])
906 {
907 const ::Tensor<1,
908 spacedim> *shape_gradient_ptr =
909 &shape_gradients[shape_function_data[shape_function]
910 .row_index[0]][0];
911
912 for (unsigned int q_point = 0;
913 q_point < n_quadrature_points;
914 ++q_point)
915 curls[q_point][0] -=
916 value * (*shape_gradient_ptr++)[1];
917 }
918
919 if (shape_function_data[shape_function]
920 .is_nonzero_shape_function_component[1])
921 {
922 const ::Tensor<1,
923 spacedim> *shape_gradient_ptr =
924 &shape_gradients[shape_function_data[shape_function]
925 .row_index[1]][0];
926
927 for (unsigned int q_point = 0;
928 q_point < n_quadrature_points;
929 ++q_point)
930 curls[q_point][0] +=
931 value * (*shape_gradient_ptr++)[0];
932 }
933 }
934 }
935 break;
936 }
937
938 case 3:
940 for (unsigned int shape_function = 0;
941 shape_function < dofs_per_cell;
942 ++shape_function)
943 {
944 const int snc = shape_function_data[shape_function]
945 .single_nonzero_component;
946
947 if (snc == -2)
948 // shape function is zero for the selected components
949 continue;
950
951 const Number &value = dof_values[shape_function];
952 // For auto-differentiable numbers, the fact that a DoF value
953 // is zero does not imply that its derivatives are zero as
954 // well. So we can't filter by value for these number types.
955 if (::internal::CheckForZero<Number>::value(value) ==
956 true)
957 continue;
958
959 if (snc != -1)
960 {
961 const ::Tensor<1, spacedim> *shape_gradient_ptr =
962 &shape_gradients[snc][0];
963
964 switch (shape_function_data[shape_function]
965 .single_nonzero_component_index)
966 {
967 case 0:
968 {
969 for (unsigned int q_point = 0;
970 q_point < n_quadrature_points;
971 ++q_point)
972 {
973 curls[q_point][1] +=
974 value * (*shape_gradient_ptr)[2];
975 curls[q_point][2] -=
976 value * (*shape_gradient_ptr++)[1];
977 }
978
979 break;
980 }
981
982 case 1:
983 {
984 for (unsigned int q_point = 0;
985 q_point < n_quadrature_points;
986 ++q_point)
987 {
988 curls[q_point][0] -=
989 value * (*shape_gradient_ptr)[2];
990 curls[q_point][2] +=
991 value * (*shape_gradient_ptr++)[0];
992 }
993
994 break;
995 }
996
997 case 2:
998 {
999 for (unsigned int q_point = 0;
1000 q_point < n_quadrature_points;
1001 ++q_point)
1002 {
1003 curls[q_point][0] +=
1004 value * (*shape_gradient_ptr)[1];
1005 curls[q_point][1] -=
1006 value * (*shape_gradient_ptr++)[0];
1007 }
1008 break;
1009 }
1010
1011 default:
1012 Assert(false, ExcInternalError());
1013 }
1014 }
1015
1016 else
1017 // we have multiple non-zero components in the shape
1018 // functions. not all of them must necessarily be within the
1019 // 3-component window this FEValuesViews::Vector object
1020 // considers, however.
1021 {
1022 if (shape_function_data[shape_function]
1023 .is_nonzero_shape_function_component[0])
1024 {
1025 const ::Tensor<1,
1026 spacedim> *shape_gradient_ptr =
1027 &shape_gradients[shape_function_data[shape_function]
1028 .row_index[0]][0];
1029
1030 for (unsigned int q_point = 0;
1031 q_point < n_quadrature_points;
1032 ++q_point)
1033 {
1034 curls[q_point][1] +=
1035 value * (*shape_gradient_ptr)[2];
1036 curls[q_point][2] -=
1037 value * (*shape_gradient_ptr++)[1];
1038 }
1039 }
1040
1041 if (shape_function_data[shape_function]
1042 .is_nonzero_shape_function_component[1])
1043 {
1044 const ::Tensor<1,
1045 spacedim> *shape_gradient_ptr =
1046 &shape_gradients[shape_function_data[shape_function]
1047 .row_index[1]][0];
1048
1049 for (unsigned int q_point = 0;
1050 q_point < n_quadrature_points;
1051 ++q_point)
1052 {
1053 curls[q_point][0] -=
1054 value * (*shape_gradient_ptr)[2];
1055 curls[q_point][2] +=
1056 value * (*shape_gradient_ptr++)[0];
1057 }
1058 }
1059
1060 if (shape_function_data[shape_function]
1061 .is_nonzero_shape_function_component[2])
1062 {
1063 const ::Tensor<1,
1064 spacedim> *shape_gradient_ptr =
1065 &shape_gradients[shape_function_data[shape_function]
1066 .row_index[2]][0];
1067
1068 for (unsigned int q_point = 0;
1069 q_point < n_quadrature_points;
1070 ++q_point)
1071 {
1072 curls[q_point][0] +=
1073 value * (*shape_gradient_ptr)[1];
1074 curls[q_point][1] -=
1075 value * (*shape_gradient_ptr++)[0];
1076 }
1077 }
1078 }
1079 }
1080 }
1081 }
1082 }
1083
1084
1085
1086 template <int dim, int spacedim, typename Number>
1087 void
1089 const ArrayView<Number> & dof_values,
1090 const Table<2, ::Tensor<2, spacedim>> &shape_hessians,
1091 const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
1092 &shape_function_data,
1093 std::vector<typename Vector<dim, spacedim>::
1094 template solution_laplacian_type<Number>> &laplacians)
1095 {
1096 const unsigned int dofs_per_cell = dof_values.size();
1097 const unsigned int n_quadrature_points = laplacians.size();
1098
1099 std::fill(
1100 laplacians.begin(),
1101 laplacians.end(),
1102 typename Vector<dim,
1103 spacedim>::template solution_laplacian_type<Number>());
1104
1105 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1106 ++shape_function)
1107 {
1108 const int snc =
1109 shape_function_data[shape_function].single_nonzero_component;
1110
1111 if (snc == -2)
1112 // shape function is zero for the selected components
1113 continue;
1114
1115 const Number &value = dof_values[shape_function];
1116 // For auto-differentiable numbers, the fact that a DoF value is zero
1117 // does not imply that its derivatives are zero as well. So we
1118 // can't filter by value for these number types.
1119 if (::internal::CheckForZero<Number>::value(value) == true)
1120 continue;
1121
1122 if (snc != -1)
1123 {
1124 const unsigned int comp = shape_function_data[shape_function]
1125 .single_nonzero_component_index;
1126 const ::Tensor<2, spacedim> *shape_hessian_ptr =
1127 &shape_hessians[snc][0];
1128 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1129 ++q_point)
1130 laplacians[q_point][comp] +=
1131 value * trace(*shape_hessian_ptr++);
1132 }
1133 else
1134 for (unsigned int d = 0; d < spacedim; ++d)
1135 if (shape_function_data[shape_function]
1136 .is_nonzero_shape_function_component[d])
1137 {
1138 const ::Tensor<2, spacedim> *shape_hessian_ptr =
1139 &shape_hessians[shape_function_data[shape_function]
1140 .row_index[d]][0];
1141 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1142 ++q_point)
1143 laplacians[q_point][d] +=
1144 value * trace(*shape_hessian_ptr++);
1145 }
1146 }
1147 }
1148
1149
1150
1151 // ---------------------- symmetric tensor part ------------------------
1152
1153 template <int dim, int spacedim, typename Number>
1154 void
1156 const ArrayView<Number> & dof_values,
1157 const ::Table<2, double> &shape_values,
1158 const std::vector<
1160 &shape_function_data,
1161 std::vector<
1162 typename ProductType<Number,
1164 &values)
1165 {
1166 const unsigned int dofs_per_cell = dof_values.size();
1167 const unsigned int n_quadrature_points = values.size();
1168
1169 std::fill(
1170 values.begin(),
1171 values.end(),
1172 typename ProductType<Number,
1174
1175 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1176 ++shape_function)
1178 const int snc =
1179 shape_function_data[shape_function].single_nonzero_component;
1180
1181 if (snc == -2)
1182 // shape function is zero for the selected components
1183 continue;
1184
1185 const Number &value = dof_values[shape_function];
1186 // For auto-differentiable numbers, the fact that a DoF value is zero
1187 // does not imply that its derivatives are zero as well. So we
1188 // can't filter by value for these number types.
1189 if (::internal::CheckForZero<Number>::value(value) == true)
1190 continue;
1191
1192 if (snc != -1)
1193 {
1194 const TableIndices<2> comp = ::
1196 shape_function_data[shape_function]
1197 .single_nonzero_component_index);
1198 const double *shape_value_ptr = &shape_values(snc, 0);
1199 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1200 ++q_point)
1201 values[q_point][comp] += value * (*shape_value_ptr++);
1202 }
1203 else
1204 for (unsigned int d = 0;
1205 d <
1207 ++d)
1208 if (shape_function_data[shape_function]
1209 .is_nonzero_shape_function_component[d])
1210 {
1211 const TableIndices<2> comp =
1214 const double *shape_value_ptr = &shape_values(
1215 shape_function_data[shape_function].row_index[d], 0);
1216 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1217 ++q_point)
1218 values[q_point][comp] += value * (*shape_value_ptr++);
1219 }
1221 }
1222
1223
1224
1225 template <int dim, int spacedim, typename Number>
1226 void
1228 const ArrayView<Number> & dof_values,
1229 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1230 const std::vector<
1232 &shape_function_data,
1233 std::vector<typename SymmetricTensor<2, dim, spacedim>::
1234 template solution_divergence_type<Number>> &divergences)
1235 {
1236 const unsigned int dofs_per_cell = dof_values.size();
1237 const unsigned int n_quadrature_points = divergences.size();
1238
1239 std::fill(divergences.begin(),
1240 divergences.end(),
1243
1244 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1245 ++shape_function)
1246 {
1247 const int snc =
1248 shape_function_data[shape_function].single_nonzero_component;
1249
1250 if (snc == -2)
1251 // shape function is zero for the selected components
1252 continue;
1253
1254 const Number &value = dof_values[shape_function];
1255 // For auto-differentiable numbers, the fact that a DoF value is zero
1256 // does not imply that its derivatives are zero as well. So we
1257 // can't filter by value for these number types.
1258 if (::internal::CheckForZero<Number>::value(value) == true)
1259 continue;
1261 if (snc != -1)
1262 {
1263 const unsigned int comp = shape_function_data[shape_function]
1264 .single_nonzero_component_index;
1265
1266 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1267 &shape_gradients[snc][0];
1268
1269 const unsigned int ii = ::SymmetricTensor<2, spacedim>::
1271 const unsigned int jj = ::SymmetricTensor<2, spacedim>::
1274 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1275 ++q_point, ++shape_gradient_ptr)
1276 {
1277 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1278
1279 if (ii != jj)
1280 divergences[q_point][jj] +=
1281 value * (*shape_gradient_ptr)[ii];
1282 }
1283 }
1284 else
1285 {
1286 for (unsigned int d = 0;
1287 d <
1289 spacedim>::n_independent_components;
1290 ++d)
1291 if (shape_function_data[shape_function]
1292 .is_nonzero_shape_function_component[d])
1293 {
1294 Assert(false, ExcNotImplemented());
1295
1296 // the following implementation needs to be looked over -- I
1297 // think it can't be right, because we are in a case where
1298 // there is no single nonzero component
1299 //
1300 // the following is not implemented! we need to consider the
1301 // interplay between multiple non-zero entries in shape
1302 // function and the representation as a symmetric
1303 // second-order tensor
1304 const unsigned int comp =
1305 shape_function_data[shape_function]
1306 .single_nonzero_component_index;
1307
1308 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1309 &shape_gradients[shape_function_data[shape_function]
1310 .row_index[d]][0];
1311 for (unsigned int q_point = 0;
1312 q_point < n_quadrature_points;
1313 ++q_point, ++shape_gradient_ptr)
1314 {
1315 for (unsigned int j = 0; j < spacedim; ++j)
1316 {
1317 const unsigned int vector_component =
1320 TableIndices<2>(comp, j));
1321 divergences[q_point][vector_component] +=
1322 value * (*shape_gradient_ptr++)[j];
1323 }
1324 }
1325 }
1326 }
1327 }
1328 }
1329
1330 // ---------------------- non-symmetric tensor part ------------------------
1331
1332 template <int dim, int spacedim, typename Number>
1333 void
1335 const ArrayView<Number> & dof_values,
1336 const ::Table<2, double> &shape_values,
1337 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1338 &shape_function_data,
1339 std::vector<
1340 typename ProductType<Number, ::Tensor<2, spacedim>>::type>
1341 &values)
1342 {
1343 const unsigned int dofs_per_cell = dof_values.size();
1344 const unsigned int n_quadrature_points = values.size();
1345
1346 std::fill(
1347 values.begin(),
1348 values.end(),
1349 typename ProductType<Number, ::Tensor<2, spacedim>>::type());
1350
1351 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1352 ++shape_function)
1353 {
1354 const int snc =
1355 shape_function_data[shape_function].single_nonzero_component;
1356
1357 if (snc == -2)
1358 // shape function is zero for the selected components
1359 continue;
1360
1361 const Number &value = dof_values[shape_function];
1362 // For auto-differentiable numbers, the fact that a DoF value is zero
1363 // does not imply that its derivatives are zero as well. So we
1364 // can't filter by value for these number types.
1365 if (::internal::CheckForZero<Number>::value(value) == true)
1366 continue;
1367
1368 if (snc != -1)
1369 {
1370 const unsigned int comp = shape_function_data[shape_function]
1371 .single_nonzero_component_index;
1372
1373 const TableIndices<2> indices =
1375 comp);
1376
1377 const double *shape_value_ptr = &shape_values(snc, 0);
1378 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1379 ++q_point)
1380 values[q_point][indices] += value * (*shape_value_ptr++);
1381 }
1382 else
1383 for (unsigned int d = 0; d < dim * dim; ++d)
1384 if (shape_function_data[shape_function]
1385 .is_nonzero_shape_function_component[d])
1387 const TableIndices<2> indices =
1389 d);
1390
1391 const double *shape_value_ptr = &shape_values(
1392 shape_function_data[shape_function].row_index[d], 0);
1393 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1394 ++q_point)
1395 values[q_point][indices] += value * (*shape_value_ptr++);
1396 }
1397 }
1398 }
1399
1400
1401
1402 template <int dim, int spacedim, typename Number>
1403 void
1405 const ArrayView<Number> & dof_values,
1406 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1407 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1408 &shape_function_data,
1409 std::vector<typename Tensor<2, dim, spacedim>::
1410 template solution_divergence_type<Number>> &divergences)
1412 const unsigned int dofs_per_cell = dof_values.size();
1413 const unsigned int n_quadrature_points = divergences.size();
1414
1415 std::fill(
1416 divergences.begin(),
1417 divergences.end(),
1419 Number>());
1420
1421 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1422 ++shape_function)
1423 {
1424 const int snc =
1425 shape_function_data[shape_function].single_nonzero_component;
1426
1427 if (snc == -2)
1428 // shape function is zero for the selected components
1429 continue;
1430
1431 const Number &value = dof_values[shape_function];
1432 // For auto-differentiable numbers, the fact that a DoF value is zero
1433 // does not imply that its derivatives are zero as well. So we
1434 // can't filter by value for these number types.
1435 if (::internal::CheckForZero<Number>::value(value) == true)
1436 continue;
1437
1438 if (snc != -1)
1439 {
1440 const unsigned int comp = shape_function_data[shape_function]
1441 .single_nonzero_component_index;
1442
1443 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1444 &shape_gradients[snc][0];
1445
1446 const TableIndices<2> indices =
1448 comp);
1449 const unsigned int ii = indices[0];
1450 const unsigned int jj = indices[1];
1451
1452 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1453 ++q_point, ++shape_gradient_ptr)
1454 {
1455 divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1456 }
1457 }
1458 else
1459 {
1460 for (unsigned int d = 0; d < dim * dim; ++d)
1461 if (shape_function_data[shape_function]
1462 .is_nonzero_shape_function_component[d])
1463 {
1464 Assert(false, ExcNotImplemented());
1465 }
1466 }
1467 }
1468 }
1469
1470
1471
1472 template <int dim, int spacedim, typename Number>
1473 void
1475 const ArrayView<Number> & dof_values,
1476 const Table<2, ::Tensor<1, spacedim>> &shape_gradients,
1477 const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
1478 &shape_function_data,
1479 std::vector<typename Tensor<2, dim, spacedim>::
1480 template solution_gradient_type<Number>> &gradients)
1481 {
1482 const unsigned int dofs_per_cell = dof_values.size();
1483 const unsigned int n_quadrature_points = gradients.size();
1484
1485 std::fill(
1486 gradients.begin(),
1487 gradients.end(),
1488 typename Tensor<2, dim, spacedim>::template solution_gradient_type<
1489 Number>());
1490
1491 for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
1492 ++shape_function)
1493 {
1494 const int snc =
1495 shape_function_data[shape_function].single_nonzero_component;
1496
1497 if (snc == -2)
1498 // shape function is zero for the selected components
1499 continue;
1500
1501 const Number &value = dof_values[shape_function];
1502 // For auto-differentiable numbers, the fact that a DoF value is zero
1503 // does not imply that its derivatives are zero as well. So we
1504 // can't filter by value for these number types.
1505 if (::internal::CheckForZero<Number>::value(value) == true)
1506 continue;
1507
1508 if (snc != -1)
1509 {
1510 const unsigned int comp = shape_function_data[shape_function]
1511 .single_nonzero_component_index;
1512
1513 const ::Tensor<1, spacedim> *shape_gradient_ptr =
1514 &shape_gradients[snc][0];
1515
1516 const TableIndices<2> indices =
1518 comp);
1519 const unsigned int ii = indices[0];
1520 const unsigned int jj = indices[1];
1521
1522 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1523 ++q_point, ++shape_gradient_ptr)
1524 {
1525 gradients[q_point][ii][jj] += value * (*shape_gradient_ptr);
1526 }
1527 }
1528 else
1529 {
1530 for (unsigned int d = 0; d < dim * dim; ++d)
1531 if (shape_function_data[shape_function]
1532 .is_nonzero_shape_function_component[d])
1533 {
1534 Assert(false, ExcNotImplemented());
1535 }
1536 }
1537 }
1538 }
1539
1540 } // end of namespace internal
1541
1542
1543
1544 template <int dim, int spacedim>
1545 template <class InputVector>
1546 void
1548 const InputVector &fe_function,
1550 const
1551 {
1552 Assert(fe_values->update_flags & update_values,
1554 "update_values")));
1555 Assert(fe_values->present_cell.get() != nullptr,
1556 ExcMessage("FEValues object is not reinit'ed to any cell"));
1557 AssertDimension(fe_function.size(),
1558 fe_values->present_cell->n_dofs_for_dof_handler());
1559
1560 // get function values of dofs on this cell and call internal worker
1561 // function
1563 fe_values->dofs_per_cell);
1564 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1565 dof_values);
1566 internal::do_function_values<dim, spacedim>(
1567 make_array_view(dof_values.begin(), dof_values.end()),
1569 shape_function_data,
1570 values);
1571 }
1572
1573
1574
1575 template <int dim, int spacedim>
1576 template <class InputVector>
1577 void
1579 const InputVector &dof_values,
1581 const
1582 {
1583 Assert(fe_values->update_flags & update_values,
1585 "update_values")));
1586 Assert(fe_values->present_cell.get() != nullptr,
1587 ExcMessage("FEValues object is not reinit'ed to any cell"));
1588 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1589
1590 internal::do_function_values<dim, spacedim>(
1591 make_array_view(dof_values.begin(), dof_values.end()),
1593 shape_function_data,
1594 values);
1595 }
1596
1597
1598
1599 template <int dim, int spacedim>
1600 template <class InputVector>
1601 void
1603 const InputVector &fe_function,
1605 &gradients) const
1606 {
1609 "update_gradients")));
1610 Assert(fe_values->present_cell.get() != nullptr,
1611 ExcMessage("FEValues object is not reinit'ed to any cell"));
1612 AssertDimension(fe_function.size(),
1613 fe_values->present_cell->n_dofs_for_dof_handler());
1614
1615 // get function values of dofs on this cell
1617 fe_values->dofs_per_cell);
1618 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1619 dof_values);
1620 internal::do_function_derivatives<1, dim, spacedim>(
1621 make_array_view(dof_values.begin(), dof_values.end()),
1623 shape_function_data,
1624 gradients);
1625 }
1626
1627
1628
1629 template <int dim, int spacedim>
1630 template <class InputVector>
1631 void
1633 const InputVector &dof_values,
1635 &gradients) const
1636 {
1639 "update_gradients")));
1640 Assert(fe_values->present_cell.get() != nullptr,
1641 ExcMessage("FEValues object is not reinit'ed to any cell"));
1642 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1643
1644 internal::do_function_derivatives<1, dim, spacedim>(
1645 make_array_view(dof_values.begin(), dof_values.end()),
1647 shape_function_data,
1648 gradients);
1649 }
1650
1651
1652
1653 template <int dim, int spacedim>
1654 template <class InputVector>
1655 void
1657 const InputVector &fe_function,
1659 &hessians) const
1660 {
1661 Assert(fe_values->update_flags & update_hessians,
1663 "update_hessians")));
1664 Assert(fe_values->present_cell.get() != nullptr,
1665 ExcMessage("FEValues object is not reinit'ed to any cell"));
1666 AssertDimension(fe_function.size(),
1667 fe_values->present_cell->n_dofs_for_dof_handler());
1668
1669 // get function values of dofs on this cell
1671 fe_values->dofs_per_cell);
1672 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1673 dof_values);
1674 internal::do_function_derivatives<2, dim, spacedim>(
1675 make_array_view(dof_values.begin(), dof_values.end()),
1677 shape_function_data,
1678 hessians);
1679 }
1680
1681
1682
1683 template <int dim, int spacedim>
1684 template <class InputVector>
1685 void
1687 const InputVector &dof_values,
1689 &hessians) const
1690 {
1691 Assert(fe_values->update_flags & update_hessians,
1693 "update_hessians")));
1694 Assert(fe_values->present_cell.get() != nullptr,
1695 ExcMessage("FEValues object is not reinit'ed to any cell"));
1696 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1697
1698 internal::do_function_derivatives<2, dim, spacedim>(
1699 make_array_view(dof_values.begin(), dof_values.end()),
1701 shape_function_data,
1702 hessians);
1703 }
1704
1705
1706
1707 template <int dim, int spacedim>
1708 template <class InputVector>
1709 void
1711 const InputVector &fe_function,
1713 &laplacians) const
1714 {
1715 Assert(fe_values->update_flags & update_hessians,
1717 "update_hessians")));
1718 Assert(fe_values->present_cell.get() != nullptr,
1719 ExcMessage("FEValues object is not reinit'ed to any cell"));
1720 AssertDimension(fe_function.size(),
1721 fe_values->present_cell->n_dofs_for_dof_handler());
1722
1723 // get function values of dofs on this cell
1725 fe_values->dofs_per_cell);
1726 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1727 dof_values);
1728 internal::do_function_laplacians<dim, spacedim>(
1729 make_array_view(dof_values.begin(), dof_values.end()),
1731 shape_function_data,
1732 laplacians);
1733 }
1734
1735
1736
1737 template <int dim, int spacedim>
1738 template <class InputVector>
1739 void
1741 const InputVector &dof_values,
1743 &laplacians) const
1744 {
1745 Assert(fe_values->update_flags & update_hessians,
1747 "update_hessians")));
1748 Assert(fe_values->present_cell.get() != nullptr,
1749 ExcMessage("FEValues object is not reinit'ed to any cell"));
1750 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1751
1752 internal::do_function_laplacians<dim, spacedim>(
1753 make_array_view(dof_values.begin(), dof_values.end()),
1755 shape_function_data,
1756 laplacians);
1757 }
1758
1759
1760
1761 template <int dim, int spacedim>
1762 template <class InputVector>
1763 void
1765 const InputVector &fe_function,
1766 std::vector<
1768 &third_derivatives) const
1769 {
1772 "update_3rd_derivatives")));
1773 Assert(fe_values->present_cell.get() != nullptr,
1774 ExcMessage("FEValues object is not reinit'ed to any cell"));
1775 AssertDimension(fe_function.size(),
1776 fe_values->present_cell->n_dofs_for_dof_handler());
1777
1778 // get function values of dofs on this cell
1780 fe_values->dofs_per_cell);
1781 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1782 dof_values);
1783 internal::do_function_derivatives<3, dim, spacedim>(
1784 make_array_view(dof_values.begin(), dof_values.end()),
1786 shape_function_data,
1787 third_derivatives);
1788 }
1789
1790
1791
1792 template <int dim, int spacedim>
1793 template <class InputVector>
1794 void
1796 const InputVector &dof_values,
1797 std::vector<
1799 &third_derivatives) const
1800 {
1803 "update_3rd_derivatives")));
1804 Assert(fe_values->present_cell.get() != nullptr,
1805 ExcMessage("FEValues object is not reinit'ed to any cell"));
1806 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1807
1808 internal::do_function_derivatives<3, dim, spacedim>(
1809 make_array_view(dof_values.begin(), dof_values.end()),
1811 shape_function_data,
1812 third_derivatives);
1813 }
1814
1815
1816
1817 template <int dim, int spacedim>
1818 template <class InputVector>
1819 void
1821 const InputVector &fe_function,
1823 const
1824 {
1825 Assert(fe_values->update_flags & update_values,
1827 "update_values")));
1828 Assert(fe_values->present_cell.get() != nullptr,
1829 ExcMessage("FEValues object is not reinit'ed to any cell"));
1830 AssertDimension(fe_function.size(),
1831 fe_values->present_cell->n_dofs_for_dof_handler());
1832
1833 // get function values of dofs on this cell
1835 fe_values->dofs_per_cell);
1836 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1837 dof_values);
1838 internal::do_function_values<dim, spacedim>(
1839 make_array_view(dof_values.begin(), dof_values.end()),
1841 shape_function_data,
1842 values);
1843 }
1844
1845
1846
1847 template <int dim, int spacedim>
1848 template <class InputVector>
1849 void
1851 const InputVector &dof_values,
1853 const
1854 {
1855 Assert(fe_values->update_flags & update_values,
1857 "update_values")));
1858 Assert(fe_values->present_cell.get() != nullptr,
1859 ExcMessage("FEValues object is not reinit'ed to any cell"));
1860 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1861
1862 internal::do_function_values<dim, spacedim>(
1863 make_array_view(dof_values.begin(), dof_values.end()),
1865 shape_function_data,
1866 values);
1867 }
1868
1869
1870
1871 template <int dim, int spacedim>
1872 template <class InputVector>
1873 void
1875 const InputVector &fe_function,
1877 &gradients) const
1878 {
1881 "update_gradients")));
1882 Assert(fe_values->present_cell.get() != nullptr,
1883 ExcMessage("FEValues object is not reinit'ed to any cell"));
1884 AssertDimension(fe_function.size(),
1885 fe_values->present_cell->n_dofs_for_dof_handler());
1886
1887 // get function values of dofs on this cell
1889 fe_values->dofs_per_cell);
1890 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1891 dof_values);
1892 internal::do_function_derivatives<1, dim, spacedim>(
1893 make_array_view(dof_values.begin(), dof_values.end()),
1895 shape_function_data,
1896 gradients);
1897 }
1898
1899
1900
1901 template <int dim, int spacedim>
1902 template <class InputVector>
1903 void
1905 const InputVector &dof_values,
1907 &gradients) const
1908 {
1911 "update_gradients")));
1912 Assert(fe_values->present_cell.get() != nullptr,
1913 ExcMessage("FEValues object is not reinit'ed to any cell"));
1914 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1915
1916 internal::do_function_derivatives<1, dim, spacedim>(
1917 make_array_view(dof_values.begin(), dof_values.end()),
1919 shape_function_data,
1920 gradients);
1921 }
1922
1923
1924
1925 template <int dim, int spacedim>
1926 template <class InputVector>
1927 void
1929 const InputVector &fe_function,
1930 std::vector<
1932 &symmetric_gradients) const
1933 {
1936 "update_gradients")));
1937 Assert(fe_values->present_cell.get() != nullptr,
1938 ExcMessage("FEValues object is not reinit'ed to any cell"));
1939 AssertDimension(fe_function.size(),
1940 fe_values->present_cell->n_dofs_for_dof_handler());
1941
1942 // get function values of dofs on this cell
1944 fe_values->dofs_per_cell);
1945 fe_values->present_cell->get_interpolated_dof_values(fe_function,
1946 dof_values);
1947 internal::do_function_symmetric_gradients<dim, spacedim>(
1948 make_array_view(dof_values.begin(), dof_values.end()),
1950 shape_function_data,
1951 symmetric_gradients);
1952 }
1953
1954
1955
1956 template <int dim, int spacedim>
1957 template <class InputVector>
1958 void
1960 const InputVector &dof_values,
1961 std::vector<
1963 &symmetric_gradients) const
1964 {
1967 "update_gradients")));
1968 Assert(fe_values->present_cell.get() != nullptr,
1969 ExcMessage("FEValues object is not reinit'ed to any cell"));
1970 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
1971
1972 internal::do_function_symmetric_gradients<dim, spacedim>(
1973 make_array_view(dof_values.begin(), dof_values.end()),
1975 shape_function_data,
1976 symmetric_gradients);
1977 }
1978
1979
1980
1981 template <int dim, int spacedim>
1982 template <class InputVector>
1983 void
1985 const InputVector &fe_function,
1987 &divergences) const
1988 {
1991 "update_gradients")));
1992 Assert(fe_values->present_cell.get() != nullptr,
1993 ExcMessage("FEValues object is not reinit'ed to any cell"));
1994 AssertDimension(fe_function.size(),
1995 fe_values->present_cell->n_dofs_for_dof_handler());
1996
1997 // get function values of dofs
1998 // on this cell
2000 fe_values->dofs_per_cell);
2001 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2002 dof_values);
2003 internal::do_function_divergences<dim, spacedim>(
2004 make_array_view(dof_values.begin(), dof_values.end()),
2006 shape_function_data,
2007 divergences);
2008 }
2009
2010
2011
2012 template <int dim, int spacedim>
2013 template <class InputVector>
2014 void
2016 const InputVector &dof_values,
2018 &divergences) const
2019 {
2022 "update_gradients")));
2023 Assert(fe_values->present_cell.get() != nullptr,
2024 ExcMessage("FEValues object is not reinit'ed to any cell"));
2025 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2026
2027 internal::do_function_divergences<dim, spacedim>(
2028 make_array_view(dof_values.begin(), dof_values.end()),
2030 shape_function_data,
2031 divergences);
2032 }
2033
2034
2035
2036 template <int dim, int spacedim>
2037 template <class InputVector>
2038 void
2040 const InputVector &fe_function,
2042 const
2043 {
2046 "update_gradients")));
2047 Assert(fe_values->present_cell.get() != nullptr,
2048 ExcMessage("FEValues object is not reinited to any cell"));
2049 AssertDimension(fe_function.size(),
2050 fe_values->present_cell->n_dofs_for_dof_handler());
2051
2052 // get function values of dofs on this cell
2054 fe_values->dofs_per_cell);
2055 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2056 dof_values);
2057 internal::do_function_curls<dim, spacedim>(
2058 make_array_view(dof_values.begin(), dof_values.end()),
2060 shape_function_data,
2061 curls);
2062 }
2063
2064
2065
2066 template <int dim, int spacedim>
2067 template <class InputVector>
2068 void
2070 const InputVector &dof_values,
2072 const
2073 {
2076 "update_gradients")));
2077 Assert(fe_values->present_cell.get() != nullptr,
2078 ExcMessage("FEValues object is not reinited to any cell"));
2079 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2080
2081 internal::do_function_curls<dim, spacedim>(
2082 make_array_view(dof_values.begin(), dof_values.end()),
2084 shape_function_data,
2085 curls);
2086 }
2087
2088
2089
2090 template <int dim, int spacedim>
2091 template <class InputVector>
2092 void
2094 const InputVector &fe_function,
2096 &hessians) const
2097 {
2098 Assert(fe_values->update_flags & update_hessians,
2100 "update_hessians")));
2101 Assert(fe_values->present_cell.get() != nullptr,
2102 ExcMessage("FEValues object is not reinit'ed to any cell"));
2103 AssertDimension(fe_function.size(),
2104 fe_values->present_cell->n_dofs_for_dof_handler());
2105
2106 // get function values of dofs on this cell
2108 fe_values->dofs_per_cell);
2109 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2110 dof_values);
2111 internal::do_function_derivatives<2, dim, spacedim>(
2112 make_array_view(dof_values.begin(), dof_values.end()),
2114 shape_function_data,
2115 hessians);
2116 }
2117
2118
2119
2120 template <int dim, int spacedim>
2121 template <class InputVector>
2122 void
2124 const InputVector &dof_values,
2126 &hessians) const
2127 {
2128 Assert(fe_values->update_flags & update_hessians,
2130 "update_hessians")));
2131 Assert(fe_values->present_cell.get() != nullptr,
2132 ExcMessage("FEValues object is not reinit'ed to any cell"));
2133 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2134
2135 internal::do_function_derivatives<2, dim, spacedim>(
2136 make_array_view(dof_values.begin(), dof_values.end()),
2138 shape_function_data,
2139 hessians);
2140 }
2141
2142
2143
2144 template <int dim, int spacedim>
2145 template <class InputVector>
2146 void
2148 const InputVector &fe_function,
2150 &laplacians) const
2151 {
2152 Assert(fe_values->update_flags & update_hessians,
2154 "update_hessians")));
2155 Assert(laplacians.size() == fe_values->n_quadrature_points,
2156 ExcDimensionMismatch(laplacians.size(),
2157 fe_values->n_quadrature_points));
2158 Assert(fe_values->present_cell.get() != nullptr,
2159 ExcMessage("FEValues object is not reinit'ed to any cell"));
2160 Assert(
2161 fe_function.size() == fe_values->present_cell->n_dofs_for_dof_handler(),
2162 ExcDimensionMismatch(fe_function.size(),
2163 fe_values->present_cell->n_dofs_for_dof_handler()));
2164
2165 // get function values of dofs on this cell
2167 fe_values->dofs_per_cell);
2168 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2169 dof_values);
2170 internal::do_function_laplacians<dim, spacedim>(
2171 make_array_view(dof_values.begin(), dof_values.end()),
2173 shape_function_data,
2174 laplacians);
2175 }
2176
2177
2178
2179 template <int dim, int spacedim>
2180 template <class InputVector>
2181 void
2183 const InputVector &dof_values,
2185 &laplacians) const
2186 {
2187 Assert(fe_values->update_flags & update_hessians,
2189 "update_hessians")));
2190 Assert(laplacians.size() == fe_values->n_quadrature_points,
2191 ExcDimensionMismatch(laplacians.size(),
2192 fe_values->n_quadrature_points));
2193 Assert(fe_values->present_cell.get() != nullptr,
2194 ExcMessage("FEValues object is not reinit'ed to any cell"));
2195 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2196
2197 internal::do_function_laplacians<dim, spacedim>(
2198 make_array_view(dof_values.begin(), dof_values.end()),
2200 shape_function_data,
2201 laplacians);
2202 }
2203
2204
2205
2206 template <int dim, int spacedim>
2207 template <class InputVector>
2208 void
2210 const InputVector &fe_function,
2211 std::vector<
2213 &third_derivatives) const
2214 {
2217 "update_3rd_derivatives")));
2218 Assert(fe_values->present_cell.get() != nullptr,
2219 ExcMessage("FEValues object is not reinit'ed to any cell"));
2220 AssertDimension(fe_function.size(),
2221 fe_values->present_cell->n_dofs_for_dof_handler());
2222
2223 // get function values of dofs on this cell
2225 fe_values->dofs_per_cell);
2226 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2227 dof_values);
2228 internal::do_function_derivatives<3, dim, spacedim>(
2229 make_array_view(dof_values.begin(), dof_values.end()),
2231 shape_function_data,
2232 third_derivatives);
2233 }
2234
2235
2236
2237 template <int dim, int spacedim>
2238 template <class InputVector>
2239 void
2241 const InputVector &dof_values,
2242 std::vector<
2244 &third_derivatives) const
2245 {
2248 "update_3rd_derivatives")));
2249 Assert(fe_values->present_cell.get() != nullptr,
2250 ExcMessage("FEValues object is not reinit'ed to any cell"));
2251 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2252
2253 internal::do_function_derivatives<3, dim, spacedim>(
2254 make_array_view(dof_values.begin(), dof_values.end()),
2256 shape_function_data,
2257 third_derivatives);
2258 }
2259
2260
2261
2262 template <int dim, int spacedim>
2263 template <class InputVector>
2264 void
2266 const InputVector &fe_function,
2268 const
2269 {
2270 Assert(fe_values->update_flags & update_values,
2272 "update_values")));
2273 Assert(fe_values->present_cell.get() != nullptr,
2274 ExcMessage("FEValues object is not reinit'ed to any cell"));
2275 AssertDimension(fe_function.size(),
2276 fe_values->present_cell->n_dofs_for_dof_handler());
2277
2278 // get function values of dofs on this cell
2280 fe_values->dofs_per_cell);
2281 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2282 dof_values);
2283 internal::do_function_values<dim, spacedim>(
2284 make_array_view(dof_values.begin(), dof_values.end()),
2286 shape_function_data,
2287 values);
2288 }
2289
2290
2291
2292 template <int dim, int spacedim>
2293 template <class InputVector>
2294 void
2296 const InputVector &dof_values,
2298 const
2299 {
2300 Assert(fe_values->update_flags & update_values,
2302 "update_values")));
2303 Assert(fe_values->present_cell.get() != nullptr,
2304 ExcMessage("FEValues object is not reinit'ed to any cell"));
2305 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2306
2307 internal::do_function_values<dim, spacedim>(
2308 make_array_view(dof_values.begin(), dof_values.end()),
2310 shape_function_data,
2311 values);
2312 }
2313
2314
2315
2316 template <int dim, int spacedim>
2317 template <class InputVector>
2318 void
2320 const InputVector &fe_function,
2322 &divergences) const
2323 {
2326 "update_gradients")));
2327 Assert(fe_values->present_cell.get() != nullptr,
2328 ExcMessage("FEValues object is not reinit'ed to any cell"));
2329 AssertDimension(fe_function.size(),
2330 fe_values->present_cell->n_dofs_for_dof_handler());
2331
2332 // get function values of dofs
2333 // on this cell
2335 fe_values->dofs_per_cell);
2336 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2337 dof_values);
2338 internal::do_function_divergences<dim, spacedim>(
2339 make_array_view(dof_values.begin(), dof_values.end()),
2341 shape_function_data,
2342 divergences);
2343 }
2344
2345
2346
2347 template <int dim, int spacedim>
2348 template <class InputVector>
2349 void
2352 const InputVector &dof_values,
2354 &divergences) const
2355 {
2358 "update_gradients")));
2359 Assert(fe_values->present_cell.get() != nullptr,
2360 ExcMessage("FEValues object is not reinit'ed to any cell"));
2361 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2362
2363 internal::do_function_divergences<dim, spacedim>(
2364 make_array_view(dof_values.begin(), dof_values.end()),
2366 shape_function_data,
2367 divergences);
2368 }
2369
2370
2371
2372 template <int dim, int spacedim>
2373 template <class InputVector>
2374 void
2376 const InputVector &fe_function,
2378 const
2379 {
2380 Assert(fe_values->update_flags & update_values,
2382 "update_values")));
2383 Assert(fe_values->present_cell.get() != nullptr,
2384 ExcMessage("FEValues object is not reinit'ed to any cell"));
2385 AssertDimension(fe_function.size(),
2386 fe_values->present_cell->n_dofs_for_dof_handler());
2387
2388 // get function values of dofs on this cell
2390 fe_values->dofs_per_cell);
2391 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2392 dof_values);
2393 internal::do_function_values<dim, spacedim>(
2394 make_array_view(dof_values.begin(), dof_values.end()),
2396 shape_function_data,
2397 values);
2398 }
2399
2400
2401
2402 template <int dim, int spacedim>
2403 template <class InputVector>
2404 void
2406 const InputVector &dof_values,
2408 const
2409 {
2410 Assert(fe_values->update_flags & update_values,
2412 "update_values")));
2413 Assert(fe_values->present_cell.get() != nullptr,
2414 ExcMessage("FEValues object is not reinit'ed to any cell"));
2415 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2416
2417 internal::do_function_values<dim, spacedim>(
2418 make_array_view(dof_values.begin(), dof_values.end()),
2420 shape_function_data,
2421 values);
2422 }
2423
2424
2425
2426 template <int dim, int spacedim>
2427 template <class InputVector>
2428 void
2430 const InputVector &fe_function,
2432 &divergences) const
2433 {
2436 "update_gradients")));
2437 Assert(fe_values->present_cell.get() != nullptr,
2438 ExcMessage("FEValues object is not reinit'ed to any cell"));
2439 AssertDimension(fe_function.size(),
2440 fe_values->present_cell->n_dofs_for_dof_handler());
2441
2442 // get function values of dofs
2443 // on this cell
2445 fe_values->dofs_per_cell);
2446 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2447 dof_values);
2448 internal::do_function_divergences<dim, spacedim>(
2449 make_array_view(dof_values.begin(), dof_values.end()),
2451 shape_function_data,
2452 divergences);
2453 }
2454
2455
2456
2457 template <int dim, int spacedim>
2458 template <class InputVector>
2459 void
2461 const InputVector &dof_values,
2462 std::vector<solution_divergence_type<typename InputVector::value_type>>
2463 &divergences) const
2464 {
2467 "update_gradients")));
2468 Assert(fe_values->present_cell.get() != nullptr,
2469 ExcMessage("FEValues object is not reinit'ed to any cell"));
2470 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2471
2472 internal::do_function_divergences<dim, spacedim>(
2473 make_array_view(dof_values.begin(), dof_values.end()),
2475 shape_function_data,
2476 divergences);
2477 }
2478
2479
2480
2481 template <int dim, int spacedim>
2482 template <class InputVector>
2483 void
2485 const InputVector &fe_function,
2487 &gradients) const
2488 {
2491 "update_gradients")));
2492 Assert(fe_values->present_cell.get() != nullptr,
2493 ExcMessage("FEValues object is not reinit'ed to any cell"));
2494 AssertDimension(fe_function.size(),
2495 fe_values->present_cell->n_dofs_for_dof_handler());
2496
2497 // get function values of dofs
2498 // on this cell
2500 fe_values->dofs_per_cell);
2501 fe_values->present_cell->get_interpolated_dof_values(fe_function,
2502 dof_values);
2503 internal::do_function_gradients<dim, spacedim>(
2504 make_array_view(dof_values.begin(), dof_values.end()),
2506 shape_function_data,
2507 gradients);
2508 }
2509
2510
2511
2512 template <int dim, int spacedim>
2513 template <class InputVector>
2514 void
2516 const InputVector &dof_values,
2518 &gradients) const
2519 {
2522 "update_gradients")));
2523 Assert(fe_values->present_cell.get() != nullptr,
2524 ExcMessage("FEValues object is not reinit'ed to any cell"));
2525 AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
2526
2527 internal::do_function_gradients<dim, spacedim>(
2528 make_array_view(dof_values.begin(), dof_values.end()),
2530 shape_function_data,
2531 gradients);
2532 }
2533
2534} // namespace FEValuesViews
2535
2536
2537namespace internal
2538{
2539 namespace FEValuesViews
2540 {
2541 template <int dim, int spacedim>
2543 {
2544 const FiniteElement<dim, spacedim> &fe = fe_values.get_fe();
2545
2546 const unsigned int n_scalars = fe.n_components();
2547 scalars.reserve(n_scalars);
2548 for (unsigned int component = 0; component < n_scalars; ++component)
2549 scalars.emplace_back(fe_values, component);
2550
2551 // compute number of vectors that we can fit into this finite element.
2552 // note that this is based on the dimensionality 'dim' of the manifold,
2553 // not 'spacedim' of the output vector
2554 const unsigned int n_vectors =
2557 1 :
2558 0);
2559 vectors.reserve(n_vectors);
2560 for (unsigned int component = 0; component < n_vectors; ++component)
2561 vectors.emplace_back(fe_values, component);
2562
2563 // compute number of symmetric tensors in the same way as above
2564 const unsigned int n_symmetric_second_order_tensors =
2565 (fe.n_components() >=
2567 fe.n_components() -
2569 0);
2570 symmetric_second_order_tensors.reserve(n_symmetric_second_order_tensors);
2571 for (unsigned int component = 0;
2572 component < n_symmetric_second_order_tensors;
2573 ++component)
2574 symmetric_second_order_tensors.emplace_back(fe_values, component);
2575
2576
2577 // compute number of symmetric tensors in the same way as above
2578 const unsigned int n_second_order_tensors =
2581 1 :
2582 0);
2583 second_order_tensors.reserve(n_second_order_tensors);
2584 for (unsigned int component = 0; component < n_second_order_tensors;
2585 ++component)
2586 second_order_tensors.emplace_back(fe_values, component);
2587 }
2588 } // namespace FEValuesViews
2589} // namespace internal
2590
2591
2592/* ---------------- FEValuesBase<dim,spacedim>::CellIteratorBase --------- */
2593
2594template <int dim, int spacedim>
2595class FEValuesBase<dim, spacedim>::CellIteratorBase
2596{
2597public:
2602 virtual ~CellIteratorBase() = default;
2603
2610 virtual
2611 operator typename Triangulation<dim, spacedim>::cell_iterator() const = 0;
2612
2619
2620#include "fe_values.decl.1.inst"
2621
2626 virtual void
2628 Vector<IndexSet::value_type> &out) const = 0;
2629};
2630
2631/* --- classes derived from FEValuesBase<dim,spacedim>::CellIteratorBase --- */
2632
2633
2638template <int dim, int spacedim>
2639template <typename CI>
2640class FEValuesBase<dim, spacedim>::CellIterator
2641 : public FEValuesBase<dim, spacedim>::CellIteratorBase
2642{
2643public:
2647 CellIterator(const CI &cell);
2648
2656 const override;
2657
2663 n_dofs_for_dof_handler() const override;
2664
2665#include "fe_values.decl.2.inst"
2666
2671 virtual void
2672 get_interpolated_dof_values(const IndexSet & in,
2673 Vector<IndexSet::value_type> &out) const override;
2674
2675private:
2679 const CI cell;
2680};
2681
2682
2701template <int dim, int spacedim>
2702class FEValuesBase<dim, spacedim>::TriaCellIterator
2703 : public FEValuesBase<dim, spacedim>::CellIteratorBase
2704{
2705public:
2709 TriaCellIterator(
2710 const typename Triangulation<dim, spacedim>::cell_iterator &cell);
2711
2720 const override;
2721
2727 n_dofs_for_dof_handler() const override;
2728
2729#include "fe_values.decl.2.inst"
2730
2735 virtual void
2736 get_interpolated_dof_values(const IndexSet & in,
2737 Vector<IndexSet::value_type> &out) const override;
2739private:
2744
2750 static const char *const message_string;
2751};
2752
2753
2754
2755/* ---------------- FEValuesBase<dim,spacedim>::CellIterator<CI> --------- */
2756
2757
2758template <int dim, int spacedim>
2759template <typename CI>
2761 : cell(cell)
2762{}
2763
2764
2765
2766template <int dim, int spacedim>
2767template <typename CI>
2769operator typename Triangulation<dim, spacedim>::cell_iterator() const
2770{
2771 return cell;
2772}
2773
2774
2775
2776template <int dim, int spacedim>
2777template <typename CI>
2780{
2781 return cell->get_dof_handler().n_dofs();
2782}
2783
2784
2785
2786#include "fe_values.impl.1.inst"
2787
2788
2789
2790template <int dim, int spacedim>
2791template <typename CI>
2792void
2794 const IndexSet & in,
2796{
2797 Assert(cell->is_active(), ExcNotImplemented());
2798
2799 std::vector<types::global_dof_index> dof_indices(
2800 cell->get_fe().n_dofs_per_cell());
2801 cell->get_dof_indices(dof_indices);
2802
2803 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
2804 out[i] = (in.is_element(dof_indices[i]) ? 1 : 0);
2805}
2806
2807
2808/* ---------------- FEValuesBase<dim,spacedim>::TriaCellIterator --------- */
2809
2810template <int dim, int spacedim>
2811const char *const FEValuesBase<dim,
2812 spacedim>::TriaCellIterator::message_string =
2813 ("You have previously called the FEValues::reinit function with a\n"
2814 "cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,\n"
2815 "when you do this, you cannot call some functions in the FEValues\n"
2816 "class, such as the get_function_values/gradients/hessians/third_derivatives\n"
2817 "functions. If you need these functions, then you need to call\n"
2818 "FEValues::reinit with an iterator type that allows to extract\n"
2819 "degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.");
2820
2821
2822
2823template <int dim, int spacedim>
2826 : cell(cell)
2827{}
2828
2829
2830
2831template <int dim, int spacedim>
2833operator typename Triangulation<dim, spacedim>::cell_iterator() const
2834{
2835 return cell;
2836}
2837
2838
2839
2840template <int dim, int spacedim>
2843{
2844 Assert(false, ExcMessage(message_string));
2845 return 0;
2846}
2847
2848
2849
2850#include "fe_values.impl.2.inst"
2851
2852
2853
2854template <int dim, int spacedim>
2855void
2857 const IndexSet &,
2859{
2860 Assert(false, ExcMessage(message_string));
2861}
2862
2863
2864
2865namespace internal
2866{
2867 namespace FEValuesImplementation
2868 {
2869 template <int dim, int spacedim>
2870 void
2872 const unsigned int n_quadrature_points,
2873 const UpdateFlags flags)
2874 {
2875 if (flags & update_quadrature_points)
2876 this->quadrature_points.resize(
2877 n_quadrature_points,
2879
2880 if (flags & update_JxW_values)
2881 this->JxW_values.resize(n_quadrature_points,
2882 numbers::signaling_nan<double>());
2883
2884 if (flags & update_jacobians)
2885 this->jacobians.resize(
2888
2889 if (flags & update_jacobian_grads)
2890 this->jacobian_grads.resize(
2893
2895 this->jacobian_pushed_forward_grads.resize(
2897
2899 this->jacobian_2nd_derivatives.resize(
2902
2904 this->jacobian_pushed_forward_2nd_derivatives.resize(
2906
2908 this->jacobian_3rd_derivatives.resize(n_quadrature_points);
2909
2911 this->jacobian_pushed_forward_3rd_derivatives.resize(
2913
2914 if (flags & update_inverse_jacobians)
2915 this->inverse_jacobians.resize(
2918
2919 if (flags & update_boundary_forms)
2920 this->boundary_forms.resize(
2922
2923 if (flags & update_normal_vectors)
2924 this->normal_vectors.resize(
2926 }
2927
2928
2929
2930 template <int dim, int spacedim>
2931 std::size_t
2933 {
2934 return (
2938 MemoryConsumption::memory_consumption(jacobian_pushed_forward_grads) +
2939 MemoryConsumption::memory_consumption(jacobian_2nd_derivatives) +
2941 jacobian_pushed_forward_2nd_derivatives) +
2942 MemoryConsumption::memory_consumption(jacobian_3rd_derivatives) +
2944 jacobian_pushed_forward_3rd_derivatives) +
2945 MemoryConsumption::memory_consumption(inverse_jacobians) +
2949 }
2950
2951
2952
2953 template <int dim, int spacedim>
2954 void
2956 const unsigned int n_quadrature_points,
2958 const UpdateFlags flags)
2959 {
2960 // initialize the table mapping from shape function number to
2961 // the rows in the tables storing the data by shape function and
2962 // nonzero component
2963 this->shape_function_to_row_table =
2965
2966 // count the total number of non-zero components accumulated
2967 // over all shape functions
2968 unsigned int n_nonzero_shape_components = 0;
2969 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
2970 n_nonzero_shape_components += fe.n_nonzero_components(i);
2971 Assert(n_nonzero_shape_components >= fe.n_dofs_per_cell(),
2973
2974 // with the number of rows now known, initialize those fields
2975 // that we will need to their correct size
2976 if (flags & update_values)
2977 {
2978 this->shape_values.reinit(n_nonzero_shape_components,
2980 this->shape_values.fill(numbers::signaling_nan<double>());
2981 }
2982
2983 if (flags & update_gradients)
2984 {
2985 this->shape_gradients.reinit(n_nonzero_shape_components,
2987 this->shape_gradients.fill(
2989 }
2990
2991 if (flags & update_hessians)
2992 {
2993 this->shape_hessians.reinit(n_nonzero_shape_components,
2995 this->shape_hessians.fill(
2997 }
2998
2999 if (flags & update_3rd_derivatives)
3001 this->shape_3rd_derivatives.reinit(n_nonzero_shape_components,
3003 this->shape_3rd_derivatives.fill(
3005 }
3006 }
3007
3008
3009
3010 template <int dim, int spacedim>
3011 std::size_t
3013 {
3014 return (
3016 MemoryConsumption::memory_consumption(shape_gradients) +
3018 MemoryConsumption::memory_consumption(shape_3rd_derivatives) +
3019 MemoryConsumption::memory_consumption(shape_function_to_row_table));
3020 }
3021 } // namespace FEValuesImplementation
3022} // namespace internal
3023
3025
3026/*------------------------------- FEValuesBase ---------------------------*/
3027
3028
3029template <int dim, int spacedim>
3031 const unsigned int n_q_points,
3032 const unsigned int dofs_per_cell,
3033 const UpdateFlags flags,
3036 : n_quadrature_points(n_q_points)
3037 , max_n_quadrature_points(n_q_points)
3039 , mapping(&mapping, typeid(*this).name())
3040 , fe(&fe, typeid(*this).name())
3042 , fe_values_views_cache(*this)
3043{
3044 Assert(n_q_points > 0,
3045 ExcMessage("There is nothing useful you can do with an FEValues "
3046 "object when using a quadrature formula with zero "
3047 "quadrature points!"));
3048 this->update_flags = flags;
3049}
3050
3051
3052
3053template <int dim, int spacedim>
3055{
3056 tria_listener_refinement.disconnect();
3058}
3059
3060
3061
3062namespace internal
3063{
3064 // put shape function part of get_function_xxx methods into separate
3065 // internal functions. this allows us to reuse the same code for several
3066 // functions (e.g. both the versions with and without indices) as well as
3067 // the same code for gradients and Hessians. Moreover, this speeds up
3068 // compilation and reduces the size of the final file since all the
3069 // different global vectors get channeled through the same code.
3070
3071 template <typename Number, typename Number2>
3072 void
3073 do_function_values(const Number2 * dof_values_ptr,
3074 const ::Table<2, double> &shape_values,
3075 std::vector<Number> & values)
3076 {
3077 // scalar finite elements, so shape_values.size() == dofs_per_cell
3078 const unsigned int dofs_per_cell = shape_values.n_rows();
3079 const unsigned int n_quadrature_points = values.size();
3080
3081 // initialize with zero
3082 std::fill_n(values.begin(),
3085
3086 // add up contributions of trial functions. note that here we deal with
3087 // scalar finite elements, so no need to check for non-primitivity of
3088 // shape functions. in order to increase the speed of this function, we
3089 // directly access the data in the shape_values array, and increment
3090 // pointers for accessing the data. this saves some lookup time and
3091 // indexing. moreover, the order of the loops is such that we can access
3092 // the shape_values data stored contiguously
3093 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
3094 {
3095 const Number2 value = dof_values_ptr[shape_func];
3096 // For auto-differentiable numbers, the fact that a DoF value is zero
3097 // does not imply that its derivatives are zero as well. So we
3098 // can't filter by value for these number types.
3101 continue;
3102
3103 const double *shape_value_ptr = &shape_values(shape_func, 0);
3104 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3105 values[point] += value * (*shape_value_ptr++);
3106 }
3108
3109
3110
3111 template <int dim, int spacedim, typename VectorType>
3112 void
3114 const typename VectorType::value_type *dof_values_ptr,
3115 const ::Table<2, double> & shape_values,
3117 const std::vector<unsigned int> & shape_function_to_row_table,
3119 const bool quadrature_points_fastest = false,
3120 const unsigned int component_multiple = 1)
3121 {
3122 using Number = typename VectorType::value_type;
3123 // initialize with zero
3124 for (unsigned int i = 0; i < values.size(); ++i)
3125 std::fill_n(values[i].begin(),
3126 values[i].size(),
3127 typename VectorType::value_type());
3128
3129 // see if there the current cell has DoFs at all, and if not
3130 // then there is nothing else to do.
3131 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3133 return;
3134
3135 const unsigned int n_quadrature_points =
3136 quadrature_points_fastest ? values[0].size() : values.size();
3137 const unsigned int n_components = fe.n_components();
3138
3139 // Assert that we can write all components into the result vectors
3140 const unsigned result_components = n_components * component_multiple;
3141 (void)result_components;
3142 if (quadrature_points_fastest)
3143 {
3144 AssertDimension(values.size(), result_components);
3145 for (unsigned int i = 0; i < values.size(); ++i)
3147 }
3148 else
3149 {
3151 for (unsigned int i = 0; i < values.size(); ++i)
3152 AssertDimension(values[i].size(), result_components);
3153 }
3154
3155 // add up contributions of trial functions. now check whether the shape
3156 // function is primitive or not. if it is, then set its only non-zero
3157 // component, otherwise loop over components
3158 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3159 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3160 ++shape_func)
3162 const Number &value = dof_values_ptr[shape_func + mc * dofs_per_cell];
3163 // For auto-differentiable numbers, the fact that a DoF value is zero
3164 // does not imply that its derivatives are zero as well. So we
3165 // can't filter by value for these number types.
3166 if (::internal::CheckForZero<Number>::value(value) == true)
3167 continue;
3168
3169 if (fe.is_primitive(shape_func))
3170 {
3171 const unsigned int comp =
3172 fe.system_to_component_index(shape_func).first +
3173 mc * n_components;
3174 const unsigned int row =
3175 shape_function_to_row_table[shape_func * n_components + comp];
3177 const double *shape_value_ptr = &shape_values(row, 0);
3178
3179 if (quadrature_points_fastest)
3180 {
3181 VectorType &values_comp = values[comp];
3182 for (unsigned int point = 0; point < n_quadrature_points;
3183 ++point)
3184 values_comp[point] += value * (*shape_value_ptr++);
3185 }
3186 else
3187 for (unsigned int point = 0; point < n_quadrature_points;
3188 ++point)
3189 values[point][comp] += value * (*shape_value_ptr++);
3190 }
3191 else
3192 for (unsigned int c = 0; c < n_components; ++c)
3193 {
3194 if (fe.get_nonzero_components(shape_func)[c] == false)
3195 continue;
3196
3197 const unsigned int row =
3198 shape_function_to_row_table[shape_func * n_components + c];
3199
3200 const double * shape_value_ptr = &shape_values(row, 0);
3201 const unsigned int comp = c + mc * n_components;
3202
3203 if (quadrature_points_fastest)
3204 {
3205 VectorType &values_comp = values[comp];
3206 for (unsigned int point = 0; point < n_quadrature_points;
3207 ++point)
3208 values_comp[point] += value * (*shape_value_ptr++);
3209 }
3210 else
3211 for (unsigned int point = 0; point < n_quadrature_points;
3212 ++point)
3213 values[point][comp] += value * (*shape_value_ptr++);
3214 }
3215 }
3216 }
3217
3218
3219
3220 // use the same implementation for gradients and Hessians, distinguish them
3221 // by the rank of the tensors
3222 template <int order, int spacedim, typename Number>
3223 void
3225 const Number * dof_values_ptr,
3226 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
3227 std::vector<Tensor<order, spacedim, Number>> & derivatives)
3228 {
3229 const unsigned int dofs_per_cell = shape_derivatives.size()[0];
3230 const unsigned int n_quadrature_points = derivatives.size();
3231
3232 // initialize with zero
3233 std::fill_n(derivatives.begin(),
3236
3237 // add up contributions of trial functions. note that here we deal with
3238 // scalar finite elements, so no need to check for non-primitivity of
3239 // shape functions. in order to increase the speed of this function, we
3240 // directly access the data in the shape_gradients/hessians array, and
3241 // increment pointers for accessing the data. this saves some lookup time
3242 // and indexing. moreover, the order of the loops is such that we can
3243 // access the shape_gradients/hessians data stored contiguously
3244 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
3245 {
3246 const Number &value = dof_values_ptr[shape_func];
3247 // For auto-differentiable numbers, the fact that a DoF value is zero
3248 // does not imply that its derivatives are zero as well. So we
3249 // can't filter by value for these number types.
3250 if (::internal::CheckForZero<Number>::value(value) == true)
3251 continue;
3252
3253 const Tensor<order, spacedim> *shape_derivative_ptr =
3254 &shape_derivatives[shape_func][0];
3255 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3256 derivatives[point] += value * (*shape_derivative_ptr++);
3257 }
3258 }
3259
3260
3261
3262 template <int order, int dim, int spacedim, typename Number>
3263 void
3265 const Number * dof_values_ptr,
3266 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
3268 const std::vector<unsigned int> &shape_function_to_row_table,
3269 ArrayView<std::vector<Tensor<order, spacedim, Number>>> derivatives,
3270 const bool quadrature_points_fastest = false,
3271 const unsigned int component_multiple = 1)
3272 {
3273 // initialize with zero
3274 for (unsigned int i = 0; i < derivatives.size(); ++i)
3275 std::fill_n(derivatives[i].begin(),
3276 derivatives[i].size(),
3278
3279 // see if there the current cell has DoFs at all, and if not
3280 // then there is nothing else to do.
3281 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3282 if (dofs_per_cell == 0)
3283 return;
3285
3286 const unsigned int n_quadrature_points =
3287 quadrature_points_fastest ? derivatives[0].size() : derivatives.size();
3288 const unsigned int n_components = fe.n_components();
3289
3290 // Assert that we can write all components into the result vectors
3291 const unsigned result_components = n_components * component_multiple;
3292 (void)result_components;
3293 if (quadrature_points_fastest)
3294 {
3295 AssertDimension(derivatives.size(), result_components);
3296 for (unsigned int i = 0; i < derivatives.size(); ++i)
3297 AssertDimension(derivatives[i].size(), n_quadrature_points);
3298 }
3299 else
3300 {
3301 AssertDimension(derivatives.size(), n_quadrature_points);
3302 for (unsigned int i = 0; i < derivatives.size(); ++i)
3303 AssertDimension(derivatives[i].size(), result_components);
3304 }
3305
3306 // add up contributions of trial functions. now check whether the shape
3307 // function is primitive or not. if it is, then set its only non-zero
3308 // component, otherwise loop over components
3309 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3310 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3311 ++shape_func)
3312 {
3313 const Number &value = dof_values_ptr[shape_func + mc * dofs_per_cell];
3314 // For auto-differentiable numbers, the fact that a DoF value is zero
3315 // does not imply that its derivatives are zero as well. So we
3316 // can't filter by value for these number types.
3317 if (::internal::CheckForZero<Number>::value(value) == true)
3318 continue;
3319
3320 if (fe.is_primitive(shape_func))
3321 {
3322 const unsigned int comp =
3323 fe.system_to_component_index(shape_func).first +
3324 mc * n_components;
3325 const unsigned int row =
3326 shape_function_to_row_table[shape_func * n_components + comp];
3327
3328 const Tensor<order, spacedim> *shape_derivative_ptr =
3329 &shape_derivatives[row][0];
3330
3331 if (quadrature_points_fastest)
3332 for (unsigned int point = 0; point < n_quadrature_points;
3333 ++point)
3334 derivatives[comp][point] += value * (*shape_derivative_ptr++);
3335 else
3336 for (unsigned int point = 0; point < n_quadrature_points;
3337 ++point)
3338 derivatives[point][comp] += value * (*shape_derivative_ptr++);
3339 }
3340 else
3341 for (unsigned int c = 0; c < n_components; ++c)
3342 {
3343 if (fe.get_nonzero_components(shape_func)[c] == false)
3344 continue;
3345
3346 const unsigned int row =
3347 shape_function_to_row_table[shape_func * n_components + c];
3348
3349 const Tensor<order, spacedim> *shape_derivative_ptr =
3350 &shape_derivatives[row][0];
3351 const unsigned int comp = c + mc * n_components;
3352
3353 if (quadrature_points_fastest)
3354 for (unsigned int point = 0; point < n_quadrature_points;
3355 ++point)
3356 derivatives[comp][point] +=
3357 value * (*shape_derivative_ptr++);
3358 else
3359 for (unsigned int point = 0; point < n_quadrature_points;
3360 ++point)
3361 derivatives[point][comp] +=
3362 value * (*shape_derivative_ptr++);
3363 }
3364 }
3365 }
3366
3367
3368
3369 template <int spacedim, typename Number, typename Number2>
3370 void
3372 const Number2 * dof_values_ptr,
3373 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
3374 std::vector<Number> & laplacians)
3375 {
3376 const unsigned int dofs_per_cell = shape_hessians.size()[0];
3377 const unsigned int n_quadrature_points = laplacians.size();
3378
3379 // initialize with zero
3380 std::fill_n(laplacians.begin(),
3383
3384 // add up contributions of trial functions. note that here we deal with
3385 // scalar finite elements and also note that the Laplacian is
3386 // the trace of the Hessian.
3387 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
3388 {
3389 const Number2 value = dof_values_ptr[shape_func];
3390 // For auto-differentiable numbers, the fact that a DoF value is zero
3391 // does not imply that its derivatives are zero as well. So we
3392 // can't filter by value for these number types.
3395 continue;
3396
3397 const Tensor<2, spacedim> *shape_hessian_ptr =
3398 &shape_hessians[shape_func][0];
3399 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3400 laplacians[point] += value * trace(*shape_hessian_ptr++);
3401 }
3402 }
3403
3404
3405
3406 template <int dim, int spacedim, typename VectorType, typename Number>
3407 void
3409 const Number * dof_values_ptr,
3410 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
3412 const std::vector<unsigned int> & shape_function_to_row_table,
3413 std::vector<VectorType> & laplacians,
3414 const bool quadrature_points_fastest = false,
3415 const unsigned int component_multiple = 1)
3416 {
3417 // initialize with zero
3418 for (unsigned int i = 0; i < laplacians.size(); ++i)
3419 std::fill_n(laplacians[i].begin(),
3420 laplacians[i].size(),
3421 typename VectorType::value_type());
3422
3423 // see if there the current cell has DoFs at all, and if not
3424 // then there is nothing else to do.
3425 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3426 if (dofs_per_cell == 0)
3427 return;
3428
3429
3430 const unsigned int n_quadrature_points = laplacians.size();
3431 const unsigned int n_components = fe.n_components();
3432
3433 // Assert that we can write all components into the result vectors
3434 const unsigned result_components = n_components * component_multiple;
3435 (void)result_components;
3436 if (quadrature_points_fastest)
3437 {
3438 AssertDimension(laplacians.size(), result_components);
3439 for (unsigned int i = 0; i < laplacians.size(); ++i)
3440 AssertDimension(laplacians[i].size(), n_quadrature_points);
3441 }
3442 else
3443 {
3444 AssertDimension(laplacians.size(), n_quadrature_points);
3445 for (unsigned int i = 0; i < laplacians.size(); ++i)
3446 AssertDimension(laplacians[i].size(), result_components);
3447 }
3448
3449 // add up contributions of trial functions. now check whether the shape
3450 // function is primitive or not. if it is, then set its only non-zero
3451 // component, otherwise loop over components
3452 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3453 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3454 ++shape_func)
3455 {
3456 const Number &value = dof_values_ptr[shape_func + mc * dofs_per_cell];
3457 // For auto-differentiable numbers, the fact that a DoF value is zero
3458 // does not imply that its derivatives are zero as well. So we
3459 // can't filter by value for these number types.
3460 if (::internal::CheckForZero<Number>::value(value) == true)
3461 continue;
3462
3463 if (fe.is_primitive(shape_func))
3464 {
3465 const unsigned int comp =
3466 fe.system_to_component_index(shape_func).first +
3467 mc * n_components;
3468 const unsigned int row =
3469 shape_function_to_row_table[shape_func * n_components + comp];
3470
3471 const Tensor<2, spacedim> *shape_hessian_ptr =
3472 &shape_hessians[row][0];
3473 if (quadrature_points_fastest)
3474 {
3475 VectorType &laplacians_comp = laplacians[comp];
3476 for (unsigned int point = 0; point < n_quadrature_points;
3477 ++point)
3478 laplacians_comp[point] +=
3479 value * trace(*shape_hessian_ptr++);
3480 }
3481 else
3482 for (unsigned int point = 0; point < n_quadrature_points;
3483 ++point)
3484 laplacians[point][comp] +=
3485 value * trace(*shape_hessian_ptr++);
3486 }
3487 else
3488 for (unsigned int c = 0; c < n_components; ++c)
3489 {
3490 if (fe.get_nonzero_components(shape_func)[c] == false)
3491 continue;
3492
3493 const unsigned int row =
3494 shape_function_to_row_table[shape_func * n_components + c];
3495
3496 const Tensor<2, spacedim> *shape_hessian_ptr =
3497 &shape_hessians[row][0];
3498 const unsigned int comp = c + mc * n_components;
3499
3500 if (quadrature_points_fastest)
3501 {
3502 VectorType &laplacians_comp = laplacians[comp];
3503 for (unsigned int point = 0; point < n_quadrature_points;
3504 ++point)
3505 laplacians_comp[point] +=
3506 value * trace(*shape_hessian_ptr++);
3507 }
3508 else
3509 for (unsigned int point = 0; point < n_quadrature_points;
3510 ++point)
3511 laplacians[point][comp] +=
3512 value * trace(*shape_hessian_ptr++);
3513 }
3514 }
3515 }
3516} // namespace internal
3517
3518
3519
3520template <int dim, int spacedim>
3521template <class InputVector>
3522void
3524 const InputVector & fe_function,
3525 std::vector<typename InputVector::value_type> &values) const
3526{
3527 using Number = typename InputVector::value_type;
3529 ExcAccessToUninitializedField("update_values"));
3530 AssertDimension(fe->n_components(), 1);
3531 Assert(present_cell.get() != nullptr,
3532 ExcMessage("FEValues object is not reinit'ed to any cell"));
3533 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3534
3535 // get function values of dofs on this cell
3536 Vector<Number> dof_values(dofs_per_cell);
3537 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3539 this->finite_element_output.shape_values,
3540 values);
3541}
3542
3543
3544
3545template <int dim, int spacedim>
3546template <class InputVector>
3547void
3549 const InputVector & fe_function,
3551 std::vector<typename InputVector::value_type> & values) const
3552{
3553 using Number = typename InputVector::value_type;
3555 ExcAccessToUninitializedField("update_values"));
3556 AssertDimension(fe->n_components(), 1);
3558
3559 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3560 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3561 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3562 internal::do_function_values(dof_values.data(),
3563 this->finite_element_output.shape_values,
3564 values);
3565}
3566
3567
3568
3569template <int dim, int spacedim>
3570template <class InputVector>
3571void
3573 const InputVector & fe_function,
3575{
3576 using Number = typename InputVector::value_type;
3577 Assert(present_cell.get() != nullptr,
3578 ExcMessage("FEValues object is not reinit'ed to any cell"));
3579
3581 ExcAccessToUninitializedField("update_values"));
3582 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3583
3584 // get function values of dofs on this cell
3585 Vector<Number> dof_values(dofs_per_cell);
3586 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3588 dof_values.begin(),
3589 this->finite_element_output.shape_values,
3590 *fe,
3591 this->finite_element_output.shape_function_to_row_table,
3592 make_array_view(values.begin(), values.end()));
3593}
3594
3595
3596
3597template <int dim, int spacedim>
3598template <class InputVector>
3599void
3601 const InputVector & fe_function,
3604{
3605 using Number = typename InputVector::value_type;
3606 // Size of indices must be a multiple of dofs_per_cell such that an integer
3607 // number of function values is generated in each point.
3608 Assert(indices.size() % dofs_per_cell == 0,
3609 ExcNotMultiple(indices.size(), dofs_per_cell));
3611 ExcAccessToUninitializedField("update_values"));
3612
3613 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3614 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3615 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3617 dof_values.data(),
3618 this->finite_element_output.shape_values,
3619 *fe,
3620 this->finite_element_output.shape_function_to_row_table,
3621 make_array_view(values.begin(), values.end()),
3622 false,
3623 indices.size() / dofs_per_cell);
3624}
3625
3626
3627
3628template <int dim, int spacedim>
3629template <class InputVector>
3630void
3632 const InputVector & fe_function,
3634 ArrayView<std::vector<typename InputVector::value_type>> values,
3635 const bool quadrature_points_fastest) const
3636{
3637 using Number = typename InputVector::value_type;
3639 ExcAccessToUninitializedField("update_values"));
3640
3641 // Size of indices must be a multiple of dofs_per_cell such that an integer
3642 // number of function values is generated in each point.
3643 Assert(indices.size() % dofs_per_cell == 0,
3644 ExcNotMultiple(indices.size(), dofs_per_cell));
3645
3646 boost::container::small_vector<Number, 200> dof_values(indices.size());
3647 for (unsigned int i = 0; i < indices.size(); ++i)
3648 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3650 dof_values.data(),
3651 this->finite_element_output.shape_values,
3652 *fe,
3653 this->finite_element_output.shape_function_to_row_table,
3654 make_array_view(values.begin(), values.end()),
3655 quadrature_points_fastest,
3656 indices.size() / dofs_per_cell);
3657}
3658
3659
3660
3661template <int dim, int spacedim>
3662template <class InputVector>
3663void
3665 const InputVector &fe_function,
3667 const
3668{
3669 using Number = typename InputVector::value_type;
3671 ExcAccessToUninitializedField("update_gradients"));
3672 AssertDimension(fe->n_components(), 1);
3673 Assert(present_cell.get() != nullptr,
3674 ExcMessage("FEValues object is not reinit'ed to any cell"));
3675 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3676
3677 // get function values of dofs on this cell
3678 Vector<Number> dof_values(dofs_per_cell);
3679 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3681 this->finite_element_output.shape_gradients,
3682 gradients);
3683}
3684
3685
3686
3687template <int dim, int spacedim>
3688template <class InputVector>
3689void
3691 const InputVector & fe_function,
3694 const
3695{
3696 using Number = typename InputVector::value_type;
3698 ExcAccessToUninitializedField("update_gradients"));
3699 AssertDimension(fe->n_components(), 1);
3701
3702 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3703 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3704 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3705 internal::do_function_derivatives(dof_values.data(),
3706 this->finite_element_output.shape_gradients,
3707 gradients);
3708}
3709
3710
3711
3712template <int dim, int spacedim>
3713template <class InputVector>
3714void
3716 const InputVector &fe_function,
3717 std::vector<
3719 &gradients) const
3720{
3721 using Number = typename InputVector::value_type;
3723 ExcAccessToUninitializedField("update_gradients"));
3724 Assert(present_cell.get() != nullptr,
3725 ExcMessage("FEValues object is not reinit'ed to any cell"));
3726 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3727
3728 // get function values of dofs on this cell
3729 Vector<Number> dof_values(dofs_per_cell);
3730 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3732 dof_values.begin(),
3733 this->finite_element_output.shape_gradients,
3735 this->finite_element_output.shape_function_to_row_table,
3736 make_array_view(gradients.begin(), gradients.end()));
3737}
3738
3739
3740
3741template <int dim, int spacedim>
3742template <class InputVector>
3743void
3745 const InputVector & fe_function,
3748 gradients,
3749 const bool quadrature_points_fastest) const
3750{
3751 using Number = typename InputVector::value_type;
3752 // Size of indices must be a multiple of dofs_per_cell such that an integer
3753 // number of function values is generated in each point.
3754 Assert(indices.size() % dofs_per_cell == 0,
3755 ExcNotMultiple(indices.size(), dofs_per_cell));
3756 Assert(this->update_flags & update_gradients,
3757 ExcAccessToUninitializedField("update_gradients"));
3758
3759 boost::container::small_vector<Number, 200> dof_values(indices.size());
3760 for (unsigned int i = 0; i < indices.size(); ++i)
3761 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3763 dof_values.data(),
3764 this->finite_element_output.shape_gradients,
3765 *fe,
3766 this->finite_element_output.shape_function_to_row_table,
3767 make_array_view(gradients.begin(), gradients.end()),
3768 quadrature_points_fastest,
3769 indices.size() / dofs_per_cell);
3770}
3771
3772
3773
3774template <int dim, int spacedim>
3775template <class InputVector>
3776void
3778 const InputVector &fe_function,
3780 const
3781{
3782 using Number = typename InputVector::value_type;
3785 ExcAccessToUninitializedField("update_hessians"));
3786 Assert(present_cell.get() != nullptr,
3787 ExcMessage("FEValues object is not reinit'ed to any cell"));
3788 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3789
3790 // get function values of dofs on this cell
3791 Vector<Number> dof_values(dofs_per_cell);
3792 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3794 this->finite_element_output.shape_hessians,
3795 hessians);
3796}
3797
3798
3799
3800template <int dim, int spacedim>
3801template <class InputVector>
3802void
3804 const InputVector & fe_function,
3807 const
3808{
3809 using Number = typename InputVector::value_type;
3811 ExcAccessToUninitializedField("update_hessians"));
3812 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3814
3815 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3816 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3817 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3818 internal::do_function_derivatives(dof_values.data(),
3819 this->finite_element_output.shape_hessians,
3820 hessians);
3821}
3822
3823
3824
3825template <int dim, int spacedim>
3826template <class InputVector>
3827void
3829 const InputVector &fe_function,
3830 std::vector<
3832 & hessians,
3833 const bool quadrature_points_fastest) const
3834{
3835 using Number = typename InputVector::value_type;
3837 ExcAccessToUninitializedField("update_hessians"));
3838 Assert(present_cell.get() != nullptr,
3839 ExcMessage("FEValues object is not reinit'ed to any cell"));
3840 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3841
3842 // get function values of dofs on this cell
3843 Vector<Number> dof_values(dofs_per_cell);
3844 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3846 dof_values.begin(),
3847 this->finite_element_output.shape_hessians,
3848 *fe,
3849 this->finite_element_output.shape_function_to_row_table,
3850 make_array_view(hessians.begin(), hessians.end()),
3851 quadrature_points_fastest);
3852}
3853
3854
3855
3856template <int dim, int spacedim>
3857template <class InputVector>
3858void
3860 const InputVector & fe_function,
3863 hessians,
3864 const bool quadrature_points_fastest) const
3865{
3866 using Number = typename InputVector::value_type;
3868 ExcAccessToUninitializedField("update_hessians"));
3869 Assert(indices.size() % dofs_per_cell == 0,
3870 ExcNotMultiple(indices.size(), dofs_per_cell));
3871
3872 boost::container::small_vector<Number, 200> dof_values(indices.size());
3873 for (unsigned int i = 0; i < indices.size(); ++i)
3874 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3876 dof_values.data(),
3877 this->finite_element_output.shape_hessians,
3878 *fe,
3879 this->finite_element_output.shape_function_to_row_table,
3880 make_array_view(hessians.begin(), hessians.end()),
3881 quadrature_points_fastest,
3882 indices.size() / dofs_per_cell);
3883}
3884
3885
3886
3887template <int dim, int spacedim>
3888template <class InputVector>
3889void
3891 const InputVector & fe_function,
3892 std::vector<typename InputVector::value_type> &laplacians) const
3893{
3894 using Number = typename InputVector::value_type;
3895 Assert(this->update_flags & update_hessians,
3896 ExcAccessToUninitializedField("update_hessians"));
3898 Assert(present_cell.get() != nullptr,
3899 ExcMessage("FEValues object is not reinit'ed to any cell"));
3900 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3901
3902 // get function values of dofs on this cell
3903 Vector<Number> dof_values(dofs_per_cell);
3904 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3906 this->finite_element_output.shape_hessians,
3907 laplacians);
3908}
3909
3910
3911
3912template <int dim, int spacedim>
3913template <class InputVector>
3914void
3916 const InputVector & fe_function,
3918 std::vector<typename InputVector::value_type> & laplacians) const
3919{
3920 using Number = typename InputVector::value_type;
3921 Assert(this->update_flags & update_hessians,
3922 ExcAccessToUninitializedField("update_hessians"));
3925
3926 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3927 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3928 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3929 internal::do_function_laplacians(dof_values.data(),
3930 this->finite_element_output.shape_hessians,
3931 laplacians);
3932}
3933
3934
3935
3936template <int dim, int spacedim>
3937template <class InputVector>
3938void
3940 const InputVector & fe_function,
3941 std::vector<Vector<typename InputVector::value_type>> &laplacians) const
3942{
3943 using Number = typename InputVector::value_type;
3944 Assert(present_cell.get() != nullptr,
3945 ExcMessage("FEValues object is not reinit'ed to any cell"));
3946 Assert(this->update_flags & update_hessians,
3947 ExcAccessToUninitializedField("update_hessians"));
3948 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
3949
3950 // get function values of dofs on this cell
3951 Vector<Number> dof_values(dofs_per_cell);
3952 present_cell->get_interpolated_dof_values(fe_function, dof_values);
3954 dof_values.begin(),
3955 this->finite_element_output.shape_hessians,
3956 *fe,
3957 this->finite_element_output.shape_function_to_row_table,
3958 laplacians);
3959}
3961
3962
3963template <int dim, int spacedim>
3964template <class InputVector>
3965void
3967 const InputVector & fe_function,
3969 std::vector<Vector<typename InputVector::value_type>> &laplacians) const
3970{
3971 using Number = typename InputVector::value_type;
3972 // Size of indices must be a multiple of dofs_per_cell such that an integer
3973 // number of function values is generated in each point.
3974 Assert(indices.size() % dofs_per_cell == 0,
3975 ExcNotMultiple(indices.size(), dofs_per_cell));
3976 Assert(this->update_flags & update_hessians,
3977 ExcAccessToUninitializedField("update_hessians"));
3978
3979 boost::container::small_vector<Number, 200> dof_values(indices.size());
3980 for (unsigned int i = 0; i < indices.size(); ++i)
3981 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3983 dof_values.data(),
3984 this->finite_element_output.shape_hessians,
3985 *fe,
3986 this->finite_element_output.shape_function_to_row_table,
3987 laplacians,
3988 false,
3989 indices.size() / dofs_per_cell);
3990}
3991
3992
3993
3994template <int dim, int spacedim>
3995template <class InputVector>
3996void
3998 const InputVector & fe_function,
4000 std::vector<std::vector<typename InputVector::value_type>> &laplacians,
4001 const bool quadrature_points_fastest) const
4002{
4003 using Number = typename InputVector::value_type;
4004 Assert(indices.size() % dofs_per_cell == 0,
4005 ExcNotMultiple(indices.size(), dofs_per_cell));
4006 Assert(this->update_flags & update_hessians,
4007 ExcAccessToUninitializedField("update_hessians"));
4008
4009 boost::container::small_vector<Number, 200> dof_values(indices.size());
4010 for (unsigned int i = 0; i < indices.size(); ++i)
4011 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
4013 dof_values.data(),
4014 this->finite_element_output.shape_hessians,
4015 *fe,
4016 this->finite_element_output.shape_function_to_row_table,
4017 laplacians,
4018 quadrature_points_fastest,
4019 indices.size() / dofs_per_cell);
4020}
4021
4022
4023
4024template <int dim, int spacedim>
4025template <class InputVector>
4026void
4028 const InputVector &fe_function,
4030 &third_derivatives) const
4031{
4032 using Number = typename InputVector::value_type;
4035 ExcAccessToUninitializedField("update_3rd_derivatives"));
4036 Assert(present_cell.get() != nullptr,
4037 ExcMessage("FEValues object is not reinit'ed to any cell"));
4038 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
4039
4040 // get function values of dofs on this cell
4041 Vector<Number> dof_values(dofs_per_cell);
4042 present_cell->get_interpolated_dof_values(fe_function, dof_values);
4044 dof_values.begin(),
4045 this->finite_element_output.shape_3rd_derivatives,
4046 third_derivatives);
4047}
4048
4049
4050
4051template <int dim, int spacedim>
4052template <class InputVector>
4053void
4055 const InputVector & fe_function,
4058 &third_derivatives) const
4059{
4060 using Number = typename InputVector::value_type;
4062 ExcAccessToUninitializedField("update_3rd_derivatives"));
4063 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
4065
4066 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
4067 for (unsigned int i = 0; i < dofs_per_cell; ++i)
4068 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
4070 dof_values.data(),
4071 this->finite_element_output.shape_3rd_derivatives,
4072 third_derivatives);
4073}
4074
4075
4076
4077template <int dim, int spacedim>
4078template <class InputVector>
4079void
4081 const InputVector &fe_function,
4082 std::vector<
4084 & third_derivatives,
4085 const bool quadrature_points_fastest) const
4086{
4087 using Number = typename InputVector::value_type;
4089 ExcAccessToUninitializedField("update_3rd_derivatives"));
4090 Assert(present_cell.get() != nullptr,
4091 ExcMessage("FEValues object is not reinit'ed to any cell"));
4092 AssertDimension(fe_function.size(), present_cell->n_dofs_for_dof_handler());
4093
4094 // get function values of dofs on this cell
4095 Vector<Number> dof_values(dofs_per_cell);
4096 present_cell->get_interpolated_dof_values(fe_function, dof_values);
4098 dof_values.begin(),
4099 this->finite_element_output.shape_3rd_derivatives,
4100 *fe,
4101 this->finite_element_output.shape_function_to_row_table,
4102 make_array_view(third_derivatives.begin(), third_derivatives.end()),
4103 quadrature_points_fastest);
4104}
4105
4106
4107
4108template <int dim, int spacedim>
4109template <class InputVector>
4110void
4112 const InputVector & fe_function,
4115 third_derivatives,
4116 const bool quadrature_points_fastest) const
4117{
4118 using Number = typename InputVector::value_type;
4120 ExcAccessToUninitializedField("update_3rd_derivatives"));
4121 Assert(indices.size() % dofs_per_cell == 0,
4122 ExcNotMultiple(indices.size(), dofs_per_cell));
4123
4124 boost::container::small_vector<Number, 200> dof_values(indices.size());
4125 for (unsigned int i = 0; i < indices.size(); ++i)
4126 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
4128 dof_values.data(),
4129 this->finite_element_output.shape_3rd_derivatives,
4130 *fe,
4131 this->finite_element_output.shape_function_to_row_table,
4132 make_array_view(third_derivatives.begin(), third_derivatives.end()),
4133 quadrature_points_fastest,
4134 indices.size() / dofs_per_cell);
4135}
4136
4137
4138
4139template <int dim, int spacedim>
4142{
4143 return *present_cell;
4144}
4145
4146
4147
4148template <int dim, int spacedim>
4149const std::vector<Tensor<1, spacedim>> &
4151{
4152 Assert(this->update_flags & update_normal_vectors,
4154 "update_normal_vectors")));
4155
4156 return this->mapping_output.normal_vectors;
4157}
4158
4159
4160
4161template <int dim, int spacedim>
4162std::size_t
4164{
4165 return (sizeof(this->update_flags) +
4166 MemoryConsumption::memory_consumption(n_quadrature_points) +
4168 sizeof(cell_similarity) +
4178}
4179
4180
4181
4182template <int dim, int spacedim>
4185 const UpdateFlags update_flags) const
4186{
4187 // first find out which objects need to be recomputed on each
4188 // cell we visit. this we have to ask the finite element and mapping.
4189 // elements are first since they might require update in mapping
4190 //
4191 // there is no need to iterate since mappings will never require
4192 // the finite element to compute something for them
4194 flags |= mapping->requires_update_flags(flags);
4195
4196 return flags;
4197}
4198
4199
4200
4201template <int dim, int spacedim>
4202void
4204{
4205 // if there is no present cell, then we shouldn't be
4206 // connected via a signal to a triangulation
4207 Assert(present_cell.get() != nullptr, ExcInternalError());
4208
4209 // so delete the present cell and
4210 // disconnect from the signal we have with
4211 // it
4212 tria_listener_refinement.disconnect();
4213 tria_listener_mesh_transform.disconnect();
4214 present_cell.reset();
4215}
4216
4217
4218
4219template <int dim, int spacedim>
4220void
4223{
4224 if (present_cell.get() != nullptr)
4225 {
4226 if (&cell->get_triangulation() !=
4228 ->
4231 {
4232 // the triangulations for the previous cell and the current cell
4233 // do not match. disconnect from the previous triangulation and
4234 // connect to the current one; also invalidate the previous
4235 // cell because we shouldn't be comparing cells from different
4236 // triangulations
4239 cell->get_triangulation().signals.any_change.connect(
4240 [this]() { this->invalidate_present_cell(); });
4242 cell->get_triangulation().signals.mesh_movement.connect(
4243 [this]() { this->invalidate_present_cell(); });
4244 }
4245 }
4246 else
4247 {
4248 // if this FEValues has never been set to any cell at all, then
4249 // at least subscribe to the triangulation to get notified of
4250 // changes
4252 cell->get_triangulation().signals.post_refinement.connect(
4253 [this]() { this->invalidate_present_cell(); });
4255 cell->get_triangulation().signals.mesh_movement.connect(
4256 [this]() { this->invalidate_present_cell(); });
4257 }
4258}
4259
4260
4261
4262template <int dim, int spacedim>
4263inline void
4266{
4267 // Unfortunately, the detection of simple geometries with CellSimilarity is
4268 // sensitive to the first cell detected. When doing this with multiple
4269 // threads, each thread will get its own scratch data object with an
4270 // FEValues object in the implementation framework from late 2013, which is
4271 // initialized to the first cell the thread sees. As this number might
4272 // different between different runs (after all, the tasks are scheduled
4273 // dynamically onto threads), this slight deviation leads to difference in
4274 // roundoff errors that propagate through the program. Therefore, we need to
4275 // disable CellSimilarity in case there is more than one thread in the
4276 // problem. This will likely not affect many MPI test cases as there
4277 // multithreading is disabled on default, but in many other situations
4278 // because we rarely explicitly set the number of threads.
4279 //
4280 // TODO: Is it reasonable to introduce a flag "unsafe" in the constructor of
4281 // FEValues to re-enable this feature?
4283 {
4285 return;
4286 }
4287
4288 // case that there has not been any cell before
4289 if (this->present_cell.get() == nullptr)
4291 else
4292 // in MappingQ, data can have been modified during the last call. Then, we
4293 // can't use that data on the new cell.
4296 else
4298 (cell->is_translation_of(
4299 static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
4300 &>(*this->present_cell)) ?
4303
4304 if ((dim < spacedim) && (cell_similarity == CellSimilarity::translation))
4305 {
4306 if (static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
4307 &>(*this->present_cell)
4308 ->direction_flag() != cell->direction_flag())
4310 }
4311 // TODO: here, one could implement other checks for similarity, e.g. for
4312 // children of a parallelogram.
4313}
4314
4315
4316
4317template <int dim, int spacedim>
4320{
4321 return cell_similarity;
4322}
4323
4324
4325
4326template <int dim, int spacedim>
4328
4329
4330
4331template <int dim, int spacedim>
4333
4334/*------------------------------- FEValues -------------------------------*/
4335
4336template <int dim, int spacedim>
4338
4339
4340
4341template <int dim, int spacedim>
4344 const Quadrature<dim> & q,
4346 : FEValuesBase<dim, spacedim>(q.size(),
4347 fe.n_dofs_per_cell(),
4349 mapping,
4350 fe)
4351 , quadrature(q)
4352{
4353 initialize(update_flags);
4354}
4355
4356
4357
4358template <int dim, int spacedim>
4361 const hp::QCollection<dim> & q,
4362 const UpdateFlags update_flags)
4363 : FEValues(mapping, fe, q[0], update_flags)
4364{
4365 AssertDimension(q.size(), 1);
4366}
4367
4368
4369
4370template <int dim, int spacedim>
4372 const Quadrature<dim> & q,
4373 const UpdateFlags update_flags)
4374 : FEValuesBase<dim, spacedim>(
4375 q.size(),
4376 fe.n_dofs_per_cell(),
4378 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4379 fe)
4380 , quadrature(q)
4381{
4382 initialize(update_flags);
4383}
4384
4385
4386
4387template <int dim, int spacedim>
4389 const hp::QCollection<dim> & q,
4390 const UpdateFlags update_flags)
4391 : FEValues(fe, q[0], update_flags)
4392{
4393 AssertDimension(q.size(), 1);
4394}
4395
4396
4397
4398template <int dim, int spacedim>
4399void
4401{
4402 // You can compute normal vectors to the cells only in the
4403 // codimension one case.
4404 if (dim != spacedim - 1)
4405 Assert((update_flags & update_normal_vectors) == false,
4406 ExcMessage("You can only pass the 'update_normal_vectors' "
4407 "flag to FEFaceValues or FESubfaceValues objects, "
4408 "but not to an FEValues object unless the "
4409 "triangulation it refers to is embedded in a higher "
4410 "dimensional space."));
4411
4412 const UpdateFlags flags = this->compute_update_flags(update_flags);
4413
4414 // initialize the base classes
4415 if (flags & update_mapping)
4416 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4417 this->finite_element_output.initialize(this->max_n_quadrature_points,
4418 *this->fe,
4419 flags);
4420
4421 // then get objects into which the FE and the Mapping can store
4422 // intermediate data used across calls to reinit. we can do this in parallel
4424 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4425 fe_get_data = Threads::new_task([&]() {
4426 return this->fe->get_data(flags,
4427 *this->mapping,
4428 quadrature,
4429 this->finite_element_output);
4430 });
4431
4433 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4434 mapping_get_data;
4435 if (flags & update_mapping)
4436 mapping_get_data = Threads::new_task(
4437 [&]() { return this->mapping->get_data(flags, quadrature); });
4438
4439 this->update_flags = flags;
4440
4441 // then collect answers from the two task above
4442 this->fe_data = std::move(fe_get_data.return_value());
4443 if (flags & update_mapping)
4444 this->mapping_data = std::move(mapping_get_data.return_value());
4445 else
4446 this->mapping_data =
4447 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4448}
4449
4450
4451
4452namespace
4453{
4454 // Reset a unique_ptr. If we can, do not de-allocate the previously
4455 // held memory but re-use it for the next item to avoid the repeated
4456 // memory allocation. We do this because FEValues objects are heavily
4457 // used in multithreaded contexts where memory allocations are evil.
4458 template <typename Type, typename Pointer, typename Iterator>
4459 void
4460 reset_pointer_in_place_if_possible(std::unique_ptr<Pointer> &present_cell,
4461 const Iterator & new_cell)
4462 {
4463 // see if the existing pointer is non-null and if the type of
4464 // the old object pointed to matches that of the one we'd
4465 // like to create
4466 if (present_cell.get() && (typeid(*present_cell.get()) == typeid(Type)))
4467 {
4468 // call destructor of the old object
4469 static_cast<const Type *>(present_cell.get())->~Type();
4470
4471 // then construct a new object in-place
4472 new (const_cast<void *>(static_cast<const void *>(present_cell.get())))
4473 Type(new_cell);
4474 }
4475 else
4476 // if the types don't match, there is nothing we can do here
4477 present_cell = std::make_unique<Type>(new_cell);
4478 }
4479} // namespace
4480
4481
4482
4483template <int dim, int spacedim>
4484void
4487{
4488 // Check that mapping and reference cell type are compatible:
4489 Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
4490 ExcMessage(
4491 "You are trying to call FEValues::reinit() with a cell of type " +
4492 cell->reference_cell().to_string() +
4493 " with a Mapping that is not compatible with it."));
4494
4495 // no FE in this cell, so no assertion
4496 // necessary here
4497 this->maybe_invalidate_previous_present_cell(cell);
4498 this->check_cell_similarity(cell);
4499
4500 reset_pointer_in_place_if_possible<
4501 typename FEValuesBase<dim, spacedim>::TriaCellIterator>(this->present_cell,
4502 cell);
4503
4504 // this was the part of the work that is dependent on the actual
4505 // data type of the iterator. now pass on to the function doing
4506 // the real work.
4507 do_reinit();
4508}
4509
4510
4511
4512template <int dim, int spacedim>
4513template <bool lda>
4514void
4517{
4518 // assert that the finite elements passed to the constructor and
4519 // used by the DoFHandler used by this cell, are the same
4520 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4521 static_cast<const FiniteElementData<dim> &>(cell->get_fe()),
4523
4524 // Check that mapping and reference cell type are compatible:
4525 Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
4526 ExcMessage(
4527 "You are trying to call FEValues::reinit() with a cell of type " +
4528 cell->reference_cell().to_string() +
4529 " with a Mapping that is not compatible with it."));
4530
4531 this->maybe_invalidate_previous_present_cell(cell);
4532 this->check_cell_similarity(cell);
4533
4534 reset_pointer_in_place_if_possible<
4535 typename FEValuesBase<dim, spacedim>::template CellIterator<
4537 cell);
4538
4539 // this was the part of the work that is dependent on the actual
4540 // data type of the iterator. now pass on to the function doing
4541 // the real work.
4542 do_reinit();
4543}
4544
4545
4546
4547template <int dim, int spacedim>
4548void
4550{
4551 // first call the mapping and let it generate the data
4552 // specific to the mapping. also let it inspect the
4553 // cell similarity flag and, if necessary, update
4554 // it
4555 if (this->update_flags & update_mapping)
4556 {
4557 this->cell_similarity =
4558 this->get_mapping().fill_fe_values(*this->present_cell,
4559 this->cell_similarity,
4560 quadrature,
4561 *this->mapping_data,
4562 this->mapping_output);
4563 }
4564
4565 // then call the finite element and, with the data
4566 // already filled by the mapping, let it compute the
4567 // data for the mapped shape function values, gradients,
4568 // etc.
4569 this->get_fe().fill_fe_values(*this->present_cell,
4570 this->cell_similarity,
4571 this->quadrature,
4572 this->get_mapping(),
4573 *this->mapping_data,
4574 this->mapping_output,
4575 *this->fe_data,
4576 this->finite_element_output);
4577}
4578
4579
4580
4581template <int dim, int spacedim>
4582std::size_t
4584{
4587}
4588
4589
4590/*------------------------------- FEFaceValuesBase --------------------------*/
4591
4592
4593template <int dim, int spacedim>
4595 const unsigned int dofs_per_cell,
4596 const UpdateFlags flags,
4597 const Mapping<dim, spacedim> & mapping,
4599 const Quadrature<dim - 1> & quadrature)
4600 : FEFaceValuesBase<dim, spacedim>(dofs_per_cell,
4601 flags,
4602 mapping,
4603 fe,
4604 hp::QCollection<dim - 1>(quadrature))
4605{}
4606
4607
4608
4609template <int dim, int spacedim>
4611 const unsigned int dofs_per_cell,
4612 const UpdateFlags,
4613 const Mapping<dim, spacedim> & mapping,
4615 const hp::QCollection<dim - 1> & quadrature)
4616 : FEValuesBase<dim, spacedim>(quadrature.max_n_quadrature_points(),
4617 dofs_per_cell,
4619 mapping,
4620 fe)
4621 , present_face_index(numbers::invalid_unsigned_int)
4622 , quadrature(quadrature)
4623{
4624 Assert(quadrature.size() == 1 ||
4625 quadrature.size() == fe.reference_cell().n_faces(),
4627}
4628
4629
4630
4631template <int dim, int spacedim>
4632const std::vector<Tensor<1, spacedim>> &
4634{
4635 Assert(this->update_flags & update_boundary_forms,
4637 "update_boundary_forms")));
4638 return this->mapping_output.boundary_forms;
4639}
4640
4641
4642
4643template <int dim, int spacedim>
4644std::size_t
4646{
4649}
4650
4651
4652/*------------------------------- FEFaceValues -------------------------------*/
4653
4654template <int dim, int spacedim>
4656
4657
4658
4659template <int dim, int spacedim>
4661
4662
4663
4664template <int dim, int spacedim>
4666 const Mapping<dim, spacedim> & mapping,
4668 const Quadrature<dim - 1> & quadrature,
4669 const UpdateFlags update_flags)
4670 : FEFaceValues<dim, spacedim>(mapping,
4671 fe,
4672 hp::QCollection<dim - 1>(quadrature),
4673 update_flags)
4674{}
4675
4676
4677
4678template <int dim, int spacedim>
4680 const Mapping<dim, spacedim> & mapping,
4682 const hp::QCollection<dim - 1> & quadrature,
4683 const UpdateFlags update_flags)
4684 : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
4685 update_flags,
4686 mapping,
4687 fe,
4688 quadrature)
4689{
4690 initialize(update_flags);
4691}
4692
4693
4694
4695template <int dim, int spacedim>
4698 const Quadrature<dim - 1> & quadrature,
4699 const UpdateFlags update_flags)
4700 : FEFaceValues<dim, spacedim>(fe,
4701 hp::QCollection<dim - 1>(quadrature),
4702 update_flags)
4703{}
4704
4705
4706
4707template <int dim, int spacedim>
4710 const hp::QCollection<dim - 1> & quadrature,
4711 const UpdateFlags update_flags)
4712 : FEFaceValuesBase<dim, spacedim>(
4713 fe.n_dofs_per_cell(),
4714 update_flags,
4715 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4716 fe,
4717 quadrature)
4718{
4719 initialize(update_flags);
4720}
4721
4722
4723
4724template <int dim, int spacedim>
4725void
4727{
4728 const UpdateFlags flags = this->compute_update_flags(update_flags);
4729
4730 // initialize the base classes
4731 if (flags & update_mapping)
4732 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4733 this->finite_element_output.initialize(this->max_n_quadrature_points,
4734 *this->fe,
4735 flags);
4736
4737 // then get objects into which the FE and the Mapping can store
4738 // intermediate data used across calls to reinit. this can be done in parallel
4739
4740 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> (
4741 FiniteElement<dim, spacedim>::*finite_element_get_face_data)(
4742 const UpdateFlags,
4743 const Mapping<dim, spacedim> &,
4746 spacedim>
4748
4749 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> (
4750 Mapping<dim, spacedim>::*mapping_get_face_data)(
4751 const UpdateFlags, const hp::QCollection<dim - 1> &) const =
4753
4754
4756 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4757 fe_get_data = Threads::new_task(finite_element_get_face_data,
4758 *this->fe,
4759 flags,
4760 *this->mapping,
4761 this->quadrature,
4762 this->finite_element_output);
4764 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4765 mapping_get_data;
4766 if (flags & update_mapping)
4767 mapping_get_data = Threads::new_task(mapping_get_face_data,
4768 *this->mapping,
4769 flags,
4770 this->quadrature);
4771
4772 this->update_flags = flags;
4773
4774 // then collect answers from the two task above
4775 this->fe_data = std::move(fe_get_data.return_value());
4776 if (flags & update_mapping)
4777 this->mapping_data = std::move(mapping_get_data.return_value());
4778 else
4779 this->mapping_data =
4780 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4781}
4782
4783
4784
4785template <int dim, int spacedim>
4786template <bool lda>
4787void
4790 const unsigned int face_no)
4791{
4792 // assert that the finite elements passed to the constructor and
4793 // used by the DoFHandler used by this cell, are the same
4794 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4795 static_cast<const FiniteElementData<dim> &>(
4796 cell->get_dof_handler().get_fe(cell->active_fe_index())),
4798
4800
4801 this->maybe_invalidate_previous_present_cell(cell);
4802 reset_pointer_in_place_if_possible<
4803 typename FEValuesBase<dim, spacedim>::template CellIterator<
4805 cell);
4806
4807 // this was the part of the work that is dependent on the actual
4808 // data type of the iterator. now pass on to the function doing
4809 // the real work.
4810 do_reinit(face_no);
4811}
4812
4813
4814
4815template <int dim, int spacedim>
4816template <bool lda>
4817void
4821{
4822 const auto face_n = cell->face_iterator_to_index(face);
4823 reinit(cell, face_n);
4824}
4825
4826
4827
4828template <int dim, int spacedim>
4829void
4832 const unsigned int face_no)
4833{
4835
4836 this->maybe_invalidate_previous_present_cell(cell);
4837 reset_pointer_in_place_if_possible<
4838 typename FEValuesBase<dim, spacedim>::TriaCellIterator>(this->present_cell,
4839 cell);
4840
4841 // this was the part of the work that is dependent on the actual
4842 // data type of the iterator. now pass on to the function doing
4843 // the real work.
4844 do_reinit(face_no);
4845}
4846
4847
4848
4849template <int dim, int spacedim>
4850void
4854{
4855 const auto face_n = cell->face_iterator_to_index(face);
4856 reinit(cell, face_n);
4857}
4858
4859
4860
4861template <int dim, int spacedim>
4862void
4863FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no)
4864{
4865 this->present_face_no = face_no;
4866
4867 // first of all, set the present_face_index (if available)
4869 *this->present_cell;
4870 this->present_face_index = cell->face_index(face_no);
4871
4872 if (this->update_flags & update_mapping)
4873 {
4874 this->get_mapping().fill_fe_face_values(*this->present_cell,
4875 face_no,
4876 this->quadrature,
4877 *this->mapping_data,
4878 this->mapping_output);
4879 }
4880
4881 this->get_fe().fill_fe_face_values(*this->present_cell,
4882 face_no,
4883 this->quadrature,
4884 this->get_mapping(),
4885 *this->mapping_data,
4886 this->mapping_output,
4887 *this->fe_data,
4888 this->finite_element_output);
4889
4890 const_cast<unsigned int &>(this->n_quadrature_points) =
4891 this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size();
4892}
4893
4894
4895/* ---------------------------- FESubFaceValues ---------------------------- */
4896
4897
4898template <int dim, int spacedim>
4900
4901
4902
4903template <int dim, int spacedim>
4905
4906
4907
4908template <int dim, int spacedim>
4910 const Mapping<dim, spacedim> & mapping,
4912 const Quadrature<dim - 1> & quadrature,
4913 const UpdateFlags update_flags)
4914 : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
4915 update_flags,
4916 mapping,
4917 fe,
4918 quadrature)
4919{
4920 initialize(update_flags);
4921}
4922
4923
4924
4925template <int dim, int spacedim>
4927 const Mapping<dim, spacedim> & mapping,
4929 const hp::QCollection<dim - 1> & quadrature,
4930 const UpdateFlags update_flags)
4931 : FESubfaceValues(mapping, fe, quadrature[0], update_flags)
4932{
4933 AssertDimension(quadrature.size(), 1);
4934}
4935
4936
4937
4938template <int dim, int spacedim>
4941 const Quadrature<dim - 1> & quadrature,
4942 const UpdateFlags update_flags)
4943 : FEFaceValuesBase<dim, spacedim>(
4944 fe.n_dofs_per_cell(),
4945 update_flags,
4946 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4947 fe,
4948 quadrature)
4949{
4950 initialize(update_flags);
4951}
4952
4953
4954
4955template <int dim, int spacedim>
4958 const hp::QCollection<dim - 1> & quadrature,
4959 const UpdateFlags update_flags)
4960 : FESubfaceValues(fe, quadrature[0], update_flags)
4961{
4962 AssertDimension(quadrature.size(), 1);
4963}
4964
4965
4966
4967template <int dim, int spacedim>
4968void
4970{
4971 const UpdateFlags flags = this->compute_update_flags(update_flags);
4972
4973 // initialize the base classes
4974 if (flags & update_mapping)
4977 *this->fe,
4978 flags);
4979
4980 // then get objects into which the FE and the Mapping can store
4981 // intermediate data used across calls to reinit. this can be done
4982 // in parallel
4984 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4985 fe_get_data =
4987 *this->fe,
4988 flags,
4989 *this->mapping,
4990 this->quadrature[0],
4991 this->finite_element_output);
4993 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4994 mapping_get_data;
4995 if (flags & update_mapping)
4996 mapping_get_data =
4998 *this->mapping,
4999 flags,
5000 this->quadrature[0]);
5001
5002 this->update_flags = flags;
5003
5004 // then collect answers from the two task above
5005 this->fe_data = std::move(fe_get_data.return_value());
5006 if (flags & update_mapping)
5007 this->mapping_data = std::move(mapping_get_data.return_value());
5008 else
5009 this->mapping_data =
5010 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
5011}
5012
5013
5014
5015template <int dim, int spacedim>
5016template <bool lda>
5017void
5020 const unsigned int face_no,
5021 const unsigned int subface_no)
5022{
5023 // assert that the finite elements passed to the constructor and
5024 // used by the DoFHandler used by this cell, are the same
5025 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
5026 static_cast<const FiniteElementData<dim> &>(
5027 cell->get_dof_handler().get_fe(cell->active_fe_index())),
5030 // We would like to check for subface_no < cell->face(face_no)->n_children(),
5031 // but unfortunately the current function is also called for
5032 // faces without children (see tests/fe/mapping.cc). Therefore,
5033 // we must use following workaround of two separate assertions
5034 Assert(cell->face(face_no)->has_children() ||
5036 ExcIndexRange(subface_no,
5037 0,
5039 Assert(!cell->face(face_no)->has_children() ||
5040 subface_no < cell->face(face_no)->n_active_descendants(),
5041 ExcIndexRange(subface_no,
5042 0,
5043 cell->face(face_no)->n_active_descendants()));
5044 Assert(cell->has_children() == false,
5045 ExcMessage("You can't use subface data for cells that are "
5046 "already refined. Iterate over their children "
5047 "instead in these cases."));
5048
5050 reset_pointer_in_place_if_possible<
5051 typename FEValuesBase<dim, spacedim>::template CellIterator<
5053 cell);
5054
5055 // this was the part of the work that is dependent on the actual
5056 // data type of the iterator. now pass on to the function doing
5057 // the real work.
5058 do_reinit(face_no, subface_no);
5059}
5060
5061
5062
5063template <int dim, int spacedim>
5064template <bool lda>
5065void
5069 const typename Triangulation<dim, spacedim>::face_iterator &subface)
5070{
5071 reinit(cell,
5072 cell->face_iterator_to_index(face),
5073 face->child_iterator_to_index(subface));
5074}
5075
5076
5077
5078template <int dim, int spacedim>
5079void
5082 const unsigned int face_no,
5083 const unsigned int subface_no)
5084{
5086 // We would like to check for subface_no < cell->face(face_no)->n_children(),
5087 // but unfortunately the current function is also called for
5088 // faces without children for periodic faces, which have hanging nodes on
5089 // the other side (see include/deal.II/matrix_free/mapping_info.templates.h).
5090 AssertIndexRange(subface_no,
5091 (cell->has_periodic_neighbor(face_no) ?
5092 cell->periodic_neighbor(face_no)
5093 ->face(cell->periodic_neighbor_face_no(face_no))
5094 ->n_children() :
5095 cell->face(face_no)->n_children()));
5096
5098 reset_pointer_in_place_if_possible<
5100 cell);
5101
5102 // this was the part of the work that is dependent on the actual
5103 // data type of the iterator. now pass on to the function doing
5104 // the real work.
5105 do_reinit(face_no, subface_no);
5106}
5107
5108
5109
5110template <int dim, int spacedim>
5111void
5115 const typename Triangulation<dim, spacedim>::face_iterator &subface)
5116{
5117 reinit(cell,
5118 cell->face_iterator_to_index(face),
5119 face->child_iterator_to_index(subface));
5120}
5121
5122
5123
5124template <int dim, int spacedim>
5125void
5126FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no,
5127 const unsigned int subface_no)
5128{
5129 this->present_face_no = face_no;
5130
5131 // first of all, set the present_face_index (if available)
5133 *this->present_cell;
5134
5135 if (!cell->face(face_no)->has_children())
5136 // no subfaces at all, so set present_face_index to this face rather
5137 // than any subface
5138 this->present_face_index = cell->face_index(face_no);
5139 else if (dim != 3)
5140 this->present_face_index = cell->face(face_no)->child_index(subface_no);
5141 else
5142 {
5143 // this is the same logic we use in cell->neighbor_child_on_subface(). See
5144 // there for an explanation of the different cases
5145 unsigned int subface_index = numbers::invalid_unsigned_int;
5146 switch (cell->subface_case(face_no))
5147 {
5151 subface_index = cell->face(face_no)->child_index(subface_no);
5152 break;
5155 subface_index = cell->face(face_no)
5156 ->child(subface_no / 2)
5157 ->child_index(subface_no % 2);
5158 break;
5161 switch (subface_no)
5162 {
5163 case 0:
5164 case 1:
5165 subface_index =
5166 cell->face(face_no)->child(0)->child_index(subface_no);
5167 break;
5168 case 2:
5169 subface_index = cell->face(face_no)->child_index(1);
5170 break;
5171 default:
5172 Assert(false, ExcInternalError());
5173 }
5174 break;
5177 switch (subface_no)
5178 {
5179 case 0:
5180 subface_index = cell->face(face_no)->child_index(0);
5181 break;
5182 case 1:
5183 case 2:
5184 subface_index =
5185 cell->face(face_no)->child(1)->child_index(subface_no - 1);
5186 break;
5187 default:
5188 Assert(false, ExcInternalError());
5189 }
5190 break;
5191 default:
5192 Assert(false, ExcInternalError());
5193 break;
5194 }
5195 Assert(subface_index != numbers::invalid_unsigned_int,
5197 this->present_face_index = subface_index;
5198 }
5199
5200 // now ask the mapping and the finite element to do the actual work
5201 if (this->update_flags & update_mapping)
5202 {
5203 this->get_mapping().fill_fe_subface_values(*this->present_cell,
5204 face_no,
5205 subface_no,
5206 this->quadrature[0],
5207 *this->mapping_data,
5208 this->mapping_output);
5209 }
5210
5211 this->get_fe().fill_fe_subface_values(*this->present_cell,
5212 face_no,
5213 subface_no,
5214 this->quadrature[0],
5215 this->get_mapping(),
5216 *this->mapping_data,
5217 this->mapping_output,
5218 *this->fe_data,
5219 this->finite_element_output);
5220}
5221
5222
5223/*------------------------------- Explicit Instantiations -------------*/
5224#define SPLIT_INSTANTIATIONS_COUNT 6
5225#ifndef SPLIT_INSTANTIATIONS_INDEX
5226# define SPLIT_INSTANTIATIONS_INDEX 0
5227#endif
5228#include "fe_values.inst"
5229
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:697
iterator begin() const
Definition: array_view.h:583
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:697
value_type * data() const noexcept
Definition: array_view.h:551
std::size_t size() const
Definition: array_view.h:574
unsigned int present_face_index
Definition: fe_values.h:4221
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
const hp::QCollection< dim - 1 > quadrature
Definition: fe_values.h:4226
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)
virtual void get_interpolated_dof_values(const IndexSet &in, Vector< IndexSet::value_type > &out) const =0
virtual types::global_dof_index n_dofs_for_dof_handler() const =0
virtual ~CellIteratorBase()=default
TriaCellIterator(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Definition: fe_values.cc:2824
virtual types::global_dof_index n_dofs_for_dof_handler() const override
Definition: fe_values.cc:2842
static const char *const message_string
Definition: fe_values.cc:2750
virtual void get_interpolated_dof_values(const IndexSet &in, Vector< IndexSet::value_type > &out) const override
Definition: fe_values.cc:2856
const Triangulation< dim, spacedim >::cell_iterator cell
Definition: fe_values.cc:2743
CellSimilarity::Similarity cell_similarity
Definition: fe_values.h:3952
::internal::FEValuesViews::Cache< dim, spacedim > fe_values_views_cache
Definition: fe_values.h:3967
boost::signals2::connection tria_listener_mesh_transform
Definition: fe_values.h:3861
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3523
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)
Definition: fe_values.cc:3030
virtual ~FEValuesBase() override
Definition: fe_values.cc:3054
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition: fe_values.h:3888
const unsigned int dofs_per_cell
Definition: fe_values.h:2450
void check_cell_similarity(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Definition: fe_values.cc:4264
UpdateFlags update_flags
Definition: fe_values.h:3934
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
Definition: fe_values.h:3912
const Mapping< dim, spacedim > & get_mapping() const
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
CellSimilarity::Similarity get_cell_similarity() const
Definition: fe_values.cc:4319
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Definition: fe_values.cc:4150
std::size_t memory_consumption() const
Definition: fe_values.cc:4163
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
const Triangulation< dim, spacedim >::cell_iterator get_cell() const
Definition: fe_values.cc:4141
void get_function_laplacians(const InputVector &fe_function, std::vector< typename InputVector::value_type > &laplacians) const
Definition: fe_values.cc:3890
::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > mapping_output
Definition: fe_values.h:3903
UpdateFlags compute_update_flags(const UpdateFlags update_flags) const
Definition: fe_values.cc:4184
std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > mapping_data
Definition: fe_values.h:3896
boost::signals2::connection tria_listener_refinement
Definition: fe_values.h:3852
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
Definition: fe_values.cc:3664
void invalidate_present_cell()
Definition: fe_values.cc:4203
std::unique_ptr< const CellIteratorBase > present_cell
Definition: fe_values.h:3843
::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > finite_element_output
Definition: fe_values.h:3928
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type > > &hessians) const
Definition: fe_values.cc:3777
std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > fe_data
Definition: fe_values.h:3920
const FiniteElement< dim, spacedim > & get_fe() const
void maybe_invalidate_previous_present_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Definition: fe_values.cc:4221
const unsigned int max_n_quadrature_points
Definition: fe_values.h:2443
void get_function_third_derivatives(const InputVector &fe_function, std::vector< Tensor< 3, spacedim, typename InputVector::value_type > > &third_derivatives) const
Definition: fe_values.cc:4027
typename ProductType< Number, hessian_type >::type solution_hessian_type
Definition: fe_values.h:214
const unsigned int component
Definition: fe_values.h:634
std::vector< ShapeFunctionData > shape_function_data
Definition: fe_values.h:639
typename ProductType< Number, value_type >::type solution_laplacian_type
Definition: fe_values.h:204
typename ProductType< Number, third_derivative_type >::type solution_third_derivative_type
Definition: fe_values.h:224
typename ProductType< Number, value_type >::type solution_value_type
Definition: fe_values.h:184
typename ProductType< Number, gradient_type >::type solution_gradient_type
Definition: fe_values.h:194
<