Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_values.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2022 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
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.is_initialized(),
1557 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1612 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1666 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1720 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1775 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1830 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1884 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1939 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
1994 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2048 ExcMessage("FEValues object is not reinited to any cell"));
2049 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2103 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
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.is_initialized(),
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.is_initialized(),
2220 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2275 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2329 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2385 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2439 AssertDimension(fe_function.size(),
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.is_initialized(),
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.is_initialized(),
2494 AssertDimension(fe_function.size(),
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.is_initialized(),
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>::CellIteratorContainer ---------
2593 */
2594
2595template <int dim, int spacedim>
2597 : initialized(false)
2598 , cell(typename Triangulation<dim, spacedim>::cell_iterator(nullptr, -1, -1))
2599 , dof_handler(nullptr)
2600 , level_dof_access(false)
2601{}
2602
2603
2604
2605template <int dim, int spacedim>
2608 : initialized(true)
2609 , cell(cell)
2610 , dof_handler(nullptr)
2611 , level_dof_access(false)
2612{}
2613
2614
2615
2616template <int dim, int spacedim>
2617bool
2619{
2620 return initialized;
2621}
2622
2623
2624
2625template <int dim, int spacedim>
2627operator typename Triangulation<dim, spacedim>::cell_iterator() const
2628{
2629 Assert(is_initialized(), ExcNotReinited());
2630
2631 return cell;
2632}
2633
2634
2635
2636template <int dim, int spacedim>
2639 const
2640{
2641 Assert(is_initialized(), ExcNotReinited());
2642 Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
2643
2644 return dof_handler->n_dofs();
2645}
2646
2647
2648
2649template <int dim, int spacedim>
2650template <typename VectorType>
2651void
2653 const VectorType & in,
2655{
2656 Assert(is_initialized(), ExcNotReinited());
2657 Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
2658
2659 if (level_dof_access)
2660 DoFCellAccessor<dim, spacedim, true>(&cell->get_triangulation(),
2661 cell->level(),
2662 cell->index(),
2663 dof_handler)
2665 else
2666 DoFCellAccessor<dim, spacedim, false>(&cell->get_triangulation(),
2667 cell->level(),
2668 cell->index(),
2669 dof_handler)
2671}
2672
2673
2674
2675template <int dim, int spacedim>
2676void
2678 const IndexSet & in,
2680{
2681 Assert(is_initialized(), ExcNotReinited());
2682 Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
2683 Assert(level_dof_access == false, ExcNotImplemented());
2684
2686 &cell->get_triangulation(), cell->level(), cell->index(), dof_handler);
2687
2688 std::vector<types::global_dof_index> dof_indices(
2689 cell_dofs.get_fe().n_dofs_per_cell());
2690 cell_dofs.get_dof_indices(dof_indices);
2691
2692 for (unsigned int i = 0; i < cell_dofs.get_fe().n_dofs_per_cell(); ++i)
2693 out[i] = (in.is_element(dof_indices[i]) ? 1 : 0);
2694}
2695
2696
2697
2698namespace internal
2699{
2700 namespace FEValuesImplementation
2701 {
2702 template <int dim, int spacedim>
2703 void
2705 const unsigned int n_quadrature_points,
2706 const UpdateFlags flags)
2707 {
2708 if (flags & update_quadrature_points)
2709 this->quadrature_points.resize(
2710 n_quadrature_points,
2712
2713 if (flags & update_JxW_values)
2714 this->JxW_values.resize(n_quadrature_points,
2715 numbers::signaling_nan<double>());
2716
2717 if (flags & update_jacobians)
2718 this->jacobians.resize(
2721
2722 if (flags & update_jacobian_grads)
2723 this->jacobian_grads.resize(
2726
2728 this->jacobian_pushed_forward_grads.resize(
2730
2732 this->jacobian_2nd_derivatives.resize(
2735
2737 this->jacobian_pushed_forward_2nd_derivatives.resize(
2739
2741 this->jacobian_3rd_derivatives.resize(n_quadrature_points);
2742
2744 this->jacobian_pushed_forward_3rd_derivatives.resize(
2746
2747 if (flags & update_inverse_jacobians)
2748 this->inverse_jacobians.resize(
2751
2752 if (flags & update_boundary_forms)
2753 this->boundary_forms.resize(
2755
2756 if (flags & update_normal_vectors)
2757 this->normal_vectors.resize(
2759 }
2760
2761
2762
2763 template <int dim, int spacedim>
2764 std::size_t
2766 {
2767 return (
2771 MemoryConsumption::memory_consumption(jacobian_pushed_forward_grads) +
2772 MemoryConsumption::memory_consumption(jacobian_2nd_derivatives) +
2774 jacobian_pushed_forward_2nd_derivatives) +
2775 MemoryConsumption::memory_consumption(jacobian_3rd_derivatives) +
2777 jacobian_pushed_forward_3rd_derivatives) +
2778 MemoryConsumption::memory_consumption(inverse_jacobians) +
2782 }
2783
2784
2785
2786 template <int dim, int spacedim>
2787 void
2789 const unsigned int n_quadrature_points,
2791 const UpdateFlags flags)
2792 {
2793 // initialize the table mapping from shape function number to
2794 // the rows in the tables storing the data by shape function and
2795 // nonzero component
2796 this->shape_function_to_row_table =
2798
2799 // count the total number of non-zero components accumulated
2800 // over all shape functions
2801 unsigned int n_nonzero_shape_components = 0;
2802 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
2803 n_nonzero_shape_components += fe.n_nonzero_components(i);
2804 Assert(n_nonzero_shape_components >= fe.n_dofs_per_cell(),
2806
2807 // with the number of rows now known, initialize those fields
2808 // that we will need to their correct size
2809 if (flags & update_values)
2810 {
2811 this->shape_values.reinit(n_nonzero_shape_components,
2813 this->shape_values.fill(numbers::signaling_nan<double>());
2814 }
2816 if (flags & update_gradients)
2817 {
2818 this->shape_gradients.reinit(n_nonzero_shape_components,
2820 this->shape_gradients.fill(
2822 }
2823
2824 if (flags & update_hessians)
2825 {
2826 this->shape_hessians.reinit(n_nonzero_shape_components,
2828 this->shape_hessians.fill(
2830 }
2831
2832 if (flags & update_3rd_derivatives)
2833 {
2834 this->shape_3rd_derivatives.reinit(n_nonzero_shape_components,
2836 this->shape_3rd_derivatives.fill(
2838 }
2839 }
2840
2841
2842
2843 template <int dim, int spacedim>
2844 std::size_t
2846 {
2847 return (
2849 MemoryConsumption::memory_consumption(shape_gradients) +
2851 MemoryConsumption::memory_consumption(shape_3rd_derivatives) +
2852 MemoryConsumption::memory_consumption(shape_function_to_row_table));
2853 }
2854 } // namespace FEValuesImplementation
2855} // namespace internal
2856
2857
2858
2859/*------------------------------- FEValuesBase ---------------------------*/
2860
2861
2862template <int dim, int spacedim>
2864 const unsigned int n_q_points,
2865 const unsigned int dofs_per_cell,
2866 const UpdateFlags flags,
2869 : n_quadrature_points(n_q_points)
2870 , max_n_quadrature_points(n_q_points)
2872 , mapping(&mapping, typeid(*this).name())
2873 , fe(&fe, typeid(*this).name())
2875 , fe_values_views_cache(*this)
2876{
2877 Assert(n_q_points > 0,
2878 ExcMessage("There is nothing useful you can do with an FEValues "
2879 "object when using a quadrature formula with zero "
2880 "quadrature points!"));
2881 this->update_flags = flags;
2882}
2883
2884
2885
2886template <int dim, int spacedim>
2888{
2889 tria_listener_refinement.disconnect();
2890 tria_listener_mesh_transform.disconnect();
2891}
2893
2894
2895namespace internal
2896{
2897 // put shape function part of get_function_xxx methods into separate
2898 // internal functions. this allows us to reuse the same code for several
2899 // functions (e.g. both the versions with and without indices) as well as
2900 // the same code for gradients and Hessians. Moreover, this speeds up
2901 // compilation and reduces the size of the final file since all the
2902 // different global vectors get channeled through the same code.
2903
2904 template <typename Number, typename Number2>
2905 void
2906 do_function_values(const Number2 * dof_values_ptr,
2907 const ::Table<2, double> &shape_values,
2908 std::vector<Number> & values)
2909 {
2910 // scalar finite elements, so shape_values.size() == dofs_per_cell
2911 const unsigned int dofs_per_cell = shape_values.n_rows();
2912 const unsigned int n_quadrature_points = values.size();
2913
2914 // initialize with zero
2915 std::fill_n(values.begin(),
2916 n_quadrature_points,
2918
2919 // add up contributions of trial functions. note that here we deal with
2920 // scalar finite elements, so no need to check for non-primitivity of
2921 // shape functions. in order to increase the speed of this function, we
2922 // directly access the data in the shape_values array, and increment
2923 // pointers for accessing the data. this saves some lookup time and
2924 // indexing. moreover, the order of the loops is such that we can access
2925 // the shape_values data stored contiguously
2926 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
2927 {
2928 const Number2 value = dof_values_ptr[shape_func];
2929 // For auto-differentiable numbers, the fact that a DoF value is zero
2930 // does not imply that its derivatives are zero as well. So we
2931 // can't filter by value for these number types.
2934 continue;
2935
2936 const double *shape_value_ptr = &shape_values(shape_func, 0);
2937 for (unsigned int point = 0; point < n_quadrature_points; ++point)
2938 values[point] += value * (*shape_value_ptr++);
2939 }
2940 }
2941
2942
2943
2944 template <int dim, int spacedim, typename VectorType>
2945 void
2947 const typename VectorType::value_type *dof_values_ptr,
2948 const ::Table<2, double> & shape_values,
2950 const std::vector<unsigned int> & shape_function_to_row_table,
2952 const bool quadrature_points_fastest = false,
2953 const unsigned int component_multiple = 1)
2954 {
2955 using Number = typename VectorType::value_type;
2956 // initialize with zero
2957 for (unsigned int i = 0; i < values.size(); ++i)
2958 std::fill_n(values[i].begin(),
2959 values[i].size(),
2960 typename VectorType::value_type());
2961
2962 // see if there the current cell has DoFs at all, and if not
2963 // then there is nothing else to do.
2964 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
2965 if (dofs_per_cell == 0)
2966 return;
2967
2968 const unsigned int n_quadrature_points =
2969 quadrature_points_fastest ? values[0].size() : values.size();
2970 const unsigned int n_components = fe.n_components();
2971
2972 // Assert that we can write all components into the result vectors
2973 const unsigned result_components = n_components * component_multiple;
2974 (void)result_components;
2975 if (quadrature_points_fastest)
2976 {
2977 AssertDimension(values.size(), result_components);
2978 for (unsigned int i = 0; i < values.size(); ++i)
2979 AssertDimension(values[i].size(), n_quadrature_points);
2980 }
2981 else
2982 {
2983 AssertDimension(values.size(), n_quadrature_points);
2984 for (unsigned int i = 0; i < values.size(); ++i)
2985 AssertDimension(values[i].size(), result_components);
2986 }
2987
2988 // add up contributions of trial functions. now check whether the shape
2989 // function is primitive or not. if it is, then set its only non-zero
2990 // component, otherwise loop over components
2991 for (unsigned int mc = 0; mc < component_multiple; ++mc)
2992 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
2993 ++shape_func)
2994 {
2995 const Number &value = dof_values_ptr[shape_func + mc * dofs_per_cell];
2996 // For auto-differentiable numbers, the fact that a DoF value is zero
2997 // does not imply that its derivatives are zero as well. So we
2998 // can't filter by value for these number types.
2999 if (::internal::CheckForZero<Number>::value(value) == true)
3000 continue;
3001
3002 if (fe.is_primitive(shape_func))
3003 {
3004 const unsigned int comp =
3005 fe.system_to_component_index(shape_func).first +
3006 mc * n_components;
3007 const unsigned int row =
3008 shape_function_to_row_table[shape_func * n_components + comp];
3009
3010 const double *shape_value_ptr = &shape_values(row, 0);
3011
3012 if (quadrature_points_fastest)
3013 {
3014 VectorType &values_comp = values[comp];
3015 for (unsigned int point = 0; point < n_quadrature_points;
3016 ++point)
3017 values_comp[point] += value * (*shape_value_ptr++);
3018 }
3019 else
3020 for (unsigned int point = 0; point < n_quadrature_points;
3021 ++point)
3022 values[point][comp] += value * (*shape_value_ptr++);
3023 }
3024 else
3025 for (unsigned int c = 0; c < n_components; ++c)
3026 {
3027 if (fe.get_nonzero_components(shape_func)[c] == false)
3028 continue;
3029
3030 const unsigned int row =
3031 shape_function_to_row_table[shape_func * n_components + c];
3032
3033 const double * shape_value_ptr = &shape_values(row, 0);
3034 const unsigned int comp = c + mc * n_components;
3035
3036 if (quadrature_points_fastest)
3037 {
3038 VectorType &values_comp = values[comp];
3039 for (unsigned int point = 0; point < n_quadrature_points;
3040 ++point)
3041 values_comp[point] += value * (*shape_value_ptr++);
3042 }
3043 else
3044 for (unsigned int point = 0; point < n_quadrature_points;
3045 ++point)
3046 values[point][comp] += value * (*shape_value_ptr++);
3047 }
3048 }
3049 }
3050
3051
3052
3053 // use the same implementation for gradients and Hessians, distinguish them
3054 // by the rank of the tensors
3055 template <int order, int spacedim, typename Number>
3056 void
3058 const Number * dof_values_ptr,
3059 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
3060 std::vector<Tensor<order, spacedim, Number>> & derivatives)
3061 {
3062 const unsigned int dofs_per_cell = shape_derivatives.size()[0];
3063 const unsigned int n_quadrature_points = derivatives.size();
3064
3065 // initialize with zero
3066 std::fill_n(derivatives.begin(),
3067 n_quadrature_points,
3069
3070 // add up contributions of trial functions. note that here we deal with
3071 // scalar finite elements, so no need to check for non-primitivity of
3072 // shape functions. in order to increase the speed of this function, we
3073 // directly access the data in the shape_gradients/hessians array, and
3074 // increment pointers for accessing the data. this saves some lookup time
3075 // and indexing. moreover, the order of the loops is such that we can
3076 // access the shape_gradients/hessians data stored contiguously
3077 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
3078 {
3079 const Number &value = dof_values_ptr[shape_func];
3080 // For auto-differentiable numbers, the fact that a DoF value is zero
3081 // does not imply that its derivatives are zero as well. So we
3082 // can't filter by value for these number types.
3083 if (::internal::CheckForZero<Number>::value(value) == true)
3084 continue;
3085
3086 const Tensor<order, spacedim> *shape_derivative_ptr =
3087 &shape_derivatives[shape_func][0];
3088 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3089 derivatives[point] += value * (*shape_derivative_ptr++);
3090 }
3091 }
3092
3093
3094
3095 template <int order, int dim, int spacedim, typename Number>
3096 void
3098 const Number * dof_values_ptr,
3099 const ::Table<2, Tensor<order, spacedim>> &shape_derivatives,
3101 const std::vector<unsigned int> &shape_function_to_row_table,
3102 ArrayView<std::vector<Tensor<order, spacedim, Number>>> derivatives,
3103 const bool quadrature_points_fastest = false,
3104 const unsigned int component_multiple = 1)
3105 {
3106 // initialize with zero
3107 for (unsigned int i = 0; i < derivatives.size(); ++i)
3108 std::fill_n(derivatives[i].begin(),
3109 derivatives[i].size(),
3111
3112 // see if there the current cell has DoFs at all, and if not
3113 // then there is nothing else to do.
3114 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3115 if (dofs_per_cell == 0)
3116 return;
3117
3118
3119 const unsigned int n_quadrature_points =
3120 quadrature_points_fastest ? derivatives[0].size() : derivatives.size();
3121 const unsigned int n_components = fe.n_components();
3122
3123 // Assert that we can write all components into the result vectors
3124 const unsigned result_components = n_components * component_multiple;
3125 (void)result_components;
3126 if (quadrature_points_fastest)
3127 {
3128 AssertDimension(derivatives.size(), result_components);
3129 for (unsigned int i = 0; i < derivatives.size(); ++i)
3130 AssertDimension(derivatives[i].size(), n_quadrature_points);
3131 }
3132 else
3133 {
3134 AssertDimension(derivatives.size(), n_quadrature_points);
3135 for (unsigned int i = 0; i < derivatives.size(); ++i)
3136 AssertDimension(derivatives[i].size(), result_components);
3137 }
3138
3139 // add up contributions of trial functions. now check whether the shape
3140 // function is primitive or not. if it is, then set its only non-zero
3141 // component, otherwise loop over components
3142 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3143 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3144 ++shape_func)
3145 {
3146 const Number &value = dof_values_ptr[shape_func + mc * dofs_per_cell];
3147 // For auto-differentiable numbers, the fact that a DoF value is zero
3148 // does not imply that its derivatives are zero as well. So we
3149 // can't filter by value for these number types.
3150 if (::internal::CheckForZero<Number>::value(value) == true)
3151 continue;
3152
3153 if (fe.is_primitive(shape_func))
3154 {
3155 const unsigned int comp =
3156 fe.system_to_component_index(shape_func).first +
3157 mc * n_components;
3158 const unsigned int row =
3159 shape_function_to_row_table[shape_func * n_components + comp];
3160
3161 const Tensor<order, spacedim> *shape_derivative_ptr =
3162 &shape_derivatives[row][0];
3163
3164 if (quadrature_points_fastest)
3165 for (unsigned int point = 0; point < n_quadrature_points;
3166 ++point)
3167 derivatives[comp][point] += value * (*shape_derivative_ptr++);
3168 else
3169 for (unsigned int point = 0; point < n_quadrature_points;
3170 ++point)
3171 derivatives[point][comp] += value * (*shape_derivative_ptr++);
3172 }
3173 else
3174 for (unsigned int c = 0; c < n_components; ++c)
3175 {
3176 if (fe.get_nonzero_components(shape_func)[c] == false)
3177 continue;
3178
3179 const unsigned int row =
3180 shape_function_to_row_table[shape_func * n_components + c];
3181
3182 const Tensor<order, spacedim> *shape_derivative_ptr =
3183 &shape_derivatives[row][0];
3184 const unsigned int comp = c + mc * n_components;
3185
3186 if (quadrature_points_fastest)
3187 for (unsigned int point = 0; point < n_quadrature_points;
3188 ++point)
3189 derivatives[comp][point] +=
3190 value * (*shape_derivative_ptr++);
3191 else
3192 for (unsigned int point = 0; point < n_quadrature_points;
3193 ++point)
3194 derivatives[point][comp] +=
3195 value * (*shape_derivative_ptr++);
3196 }
3197 }
3198 }
3199
3200
3201
3202 template <int spacedim, typename Number, typename Number2>
3203 void
3205 const Number2 * dof_values_ptr,
3206 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
3207 std::vector<Number> & laplacians)
3208 {
3209 const unsigned int dofs_per_cell = shape_hessians.size()[0];
3210 const unsigned int n_quadrature_points = laplacians.size();
3211
3212 // initialize with zero
3213 std::fill_n(laplacians.begin(),
3214 n_quadrature_points,
3216
3217 // add up contributions of trial functions. note that here we deal with
3218 // scalar finite elements and also note that the Laplacian is
3219 // the trace of the Hessian.
3220 for (unsigned int shape_func = 0; shape_func < dofs_per_cell; ++shape_func)
3221 {
3222 const Number2 value = dof_values_ptr[shape_func];
3223 // For auto-differentiable numbers, the fact that a DoF value is zero
3224 // does not imply that its derivatives are zero as well. So we
3225 // can't filter by value for these number types.
3228 continue;
3229
3230 const Tensor<2, spacedim> *shape_hessian_ptr =
3231 &shape_hessians[shape_func][0];
3232 for (unsigned int point = 0; point < n_quadrature_points; ++point)
3233 laplacians[point] += value * trace(*shape_hessian_ptr++);
3234 }
3235 }
3236
3237
3238
3239 template <int dim, int spacedim, typename VectorType, typename Number>
3240 void
3242 const Number * dof_values_ptr,
3243 const ::Table<2, Tensor<2, spacedim>> &shape_hessians,
3245 const std::vector<unsigned int> & shape_function_to_row_table,
3246 std::vector<VectorType> & laplacians,
3247 const bool quadrature_points_fastest = false,
3248 const unsigned int component_multiple = 1)
3249 {
3250 // initialize with zero
3251 for (unsigned int i = 0; i < laplacians.size(); ++i)
3252 std::fill_n(laplacians[i].begin(),
3253 laplacians[i].size(),
3254 typename VectorType::value_type());
3255
3256 // see if there the current cell has DoFs at all, and if not
3257 // then there is nothing else to do.
3258 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
3259 if (dofs_per_cell == 0)
3260 return;
3261
3262
3263 const unsigned int n_quadrature_points = laplacians.size();
3264 const unsigned int n_components = fe.n_components();
3265
3266 // Assert that we can write all components into the result vectors
3267 const unsigned result_components = n_components * component_multiple;
3268 (void)result_components;
3269 if (quadrature_points_fastest)
3270 {
3271 AssertDimension(laplacians.size(), result_components);
3272 for (unsigned int i = 0; i < laplacians.size(); ++i)
3273 AssertDimension(laplacians[i].size(), n_quadrature_points);
3274 }
3275 else
3276 {
3277 AssertDimension(laplacians.size(), n_quadrature_points);
3278 for (unsigned int i = 0; i < laplacians.size(); ++i)
3279 AssertDimension(laplacians[i].size(), result_components);
3280 }
3281
3282 // add up contributions of trial functions. now check whether the shape
3283 // function is primitive or not. if it is, then set its only non-zero
3284 // component, otherwise loop over components
3285 for (unsigned int mc = 0; mc < component_multiple; ++mc)
3286 for (unsigned int shape_func = 0; shape_func < dofs_per_cell;
3287 ++shape_func)
3288 {
3289 const Number &value = dof_values_ptr[shape_func + mc * dofs_per_cell];
3290 // For auto-differentiable numbers, the fact that a DoF value is zero
3291 // does not imply that its derivatives are zero as well. So we
3292 // can't filter by value for these number types.
3293 if (::internal::CheckForZero<Number>::value(value) == true)
3294 continue;
3295
3296 if (fe.is_primitive(shape_func))
3297 {
3298 const unsigned int comp =
3299 fe.system_to_component_index(shape_func).first +
3300 mc * n_components;
3301 const unsigned int row =
3302 shape_function_to_row_table[shape_func * n_components + comp];
3303
3304 const Tensor<2, spacedim> *shape_hessian_ptr =
3305 &shape_hessians[row][0];
3306 if (quadrature_points_fastest)
3307 {
3308 VectorType &laplacians_comp = laplacians[comp];
3309 for (unsigned int point = 0; point < n_quadrature_points;
3310 ++point)
3311 laplacians_comp[point] +=
3312 value * trace(*shape_hessian_ptr++);
3313 }
3314 else
3315 for (unsigned int point = 0; point < n_quadrature_points;
3316 ++point)
3317 laplacians[point][comp] +=
3318 value * trace(*shape_hessian_ptr++);
3319 }
3320 else
3321 for (unsigned int c = 0; c < n_components; ++c)
3322 {
3323 if (fe.get_nonzero_components(shape_func)[c] == false)
3324 continue;
3325
3326 const unsigned int row =
3327 shape_function_to_row_table[shape_func * n_components + c];
3328
3329 const Tensor<2, spacedim> *shape_hessian_ptr =
3330 &shape_hessians[row][0];
3331 const unsigned int comp = c + mc * n_components;
3332
3333 if (quadrature_points_fastest)
3334 {
3335 VectorType &laplacians_comp = laplacians[comp];
3336 for (unsigned int point = 0; point < n_quadrature_points;
3337 ++point)
3338 laplacians_comp[point] +=
3339 value * trace(*shape_hessian_ptr++);
3340 }
3341 else
3342 for (unsigned int point = 0; point < n_quadrature_points;
3343 ++point)
3344 laplacians[point][comp] +=
3345 value * trace(*shape_hessian_ptr++);
3346 }
3347 }
3348 }
3349} // namespace internal
3350
3351
3352
3353template <int dim, int spacedim>
3354template <class InputVector>
3355void
3357 const InputVector & fe_function,
3358 std::vector<typename InputVector::value_type> &values) const
3359{
3360 using Number = typename InputVector::value_type;
3361 Assert(this->update_flags & update_values,
3362 ExcAccessToUninitializedField("update_values"));
3363 AssertDimension(fe->n_components(), 1);
3364 Assert(present_cell.is_initialized(), ExcNotReinited());
3365 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3366
3367 // get function values of dofs on this cell
3368 Vector<Number> dof_values(dofs_per_cell);
3369 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3371 this->finite_element_output.shape_values,
3372 values);
3373}
3374
3375
3376
3377template <int dim, int spacedim>
3378template <class InputVector>
3379void
3381 const InputVector & fe_function,
3383 std::vector<typename InputVector::value_type> & values) const
3384{
3385 using Number = typename InputVector::value_type;
3386 Assert(this->update_flags & update_values,
3387 ExcAccessToUninitializedField("update_values"));
3388 AssertDimension(fe->n_components(), 1);
3389 AssertDimension(indices.size(), dofs_per_cell);
3390
3391 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3392 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3393 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3394 internal::do_function_values(dof_values.data(),
3395 this->finite_element_output.shape_values,
3396 values);
3397}
3398
3399
3400
3401template <int dim, int spacedim>
3402template <class InputVector>
3403void
3405 const InputVector & fe_function,
3407{
3408 using Number = typename InputVector::value_type;
3409 Assert(present_cell.is_initialized(), ExcNotReinited());
3410
3411 Assert(this->update_flags & update_values,
3412 ExcAccessToUninitializedField("update_values"));
3413 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3414
3415 // get function values of dofs on this cell
3416 Vector<Number> dof_values(dofs_per_cell);
3417 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3419 dof_values.begin(),
3420 this->finite_element_output.shape_values,
3421 *fe,
3422 this->finite_element_output.shape_function_to_row_table,
3423 make_array_view(values.begin(), values.end()));
3424}
3425
3426
3427
3428template <int dim, int spacedim>
3429template <class InputVector>
3430void
3432 const InputVector & fe_function,
3435{
3436 using Number = typename InputVector::value_type;
3437 // Size of indices must be a multiple of dofs_per_cell such that an integer
3438 // number of function values is generated in each point.
3439 Assert(indices.size() % dofs_per_cell == 0,
3440 ExcNotMultiple(indices.size(), dofs_per_cell));
3441 Assert(this->update_flags & update_values,
3442 ExcAccessToUninitializedField("update_values"));
3443
3444 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3445 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3446 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3448 dof_values.data(),
3449 this->finite_element_output.shape_values,
3450 *fe,
3451 this->finite_element_output.shape_function_to_row_table,
3452 make_array_view(values.begin(), values.end()),
3453 false,
3454 indices.size() / dofs_per_cell);
3455}
3456
3457
3458
3459template <int dim, int spacedim>
3460template <class InputVector>
3461void
3463 const InputVector & fe_function,
3465 ArrayView<std::vector<typename InputVector::value_type>> values,
3466 const bool quadrature_points_fastest) const
3467{
3468 using Number = typename InputVector::value_type;
3469 Assert(this->update_flags & update_values,
3470 ExcAccessToUninitializedField("update_values"));
3471
3472 // Size of indices must be a multiple of dofs_per_cell such that an integer
3473 // number of function values is generated in each point.
3474 Assert(indices.size() % dofs_per_cell == 0,
3475 ExcNotMultiple(indices.size(), dofs_per_cell));
3476
3477 boost::container::small_vector<Number, 200> dof_values(indices.size());
3478 for (unsigned int i = 0; i < indices.size(); ++i)
3479 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3481 dof_values.data(),
3482 this->finite_element_output.shape_values,
3483 *fe,
3484 this->finite_element_output.shape_function_to_row_table,
3485 make_array_view(values.begin(), values.end()),
3486 quadrature_points_fastest,
3487 indices.size() / dofs_per_cell);
3488}
3489
3490
3491
3492template <int dim, int spacedim>
3493template <class InputVector>
3494void
3496 const InputVector &fe_function,
3498 const
3499{
3500 using Number = typename InputVector::value_type;
3501 Assert(this->update_flags & update_gradients,
3502 ExcAccessToUninitializedField("update_gradients"));
3503 AssertDimension(fe->n_components(), 1);
3504 Assert(present_cell.is_initialized(), ExcNotReinited());
3505 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3506
3507 // get function values of dofs on this cell
3508 Vector<Number> dof_values(dofs_per_cell);
3509 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3511 this->finite_element_output.shape_gradients,
3512 gradients);
3513}
3514
3515
3516
3517template <int dim, int spacedim>
3518template <class InputVector>
3519void
3521 const InputVector & fe_function,
3524 const
3525{
3526 using Number = typename InputVector::value_type;
3527 Assert(this->update_flags & update_gradients,
3528 ExcAccessToUninitializedField("update_gradients"));
3529 AssertDimension(fe->n_components(), 1);
3530 AssertDimension(indices.size(), dofs_per_cell);
3531
3532 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3533 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3534 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3535 internal::do_function_derivatives(dof_values.data(),
3536 this->finite_element_output.shape_gradients,
3537 gradients);
3538}
3539
3540
3541
3542template <int dim, int spacedim>
3543template <class InputVector>
3544void
3546 const InputVector &fe_function,
3547 std::vector<
3549 &gradients) const
3550{
3551 using Number = typename InputVector::value_type;
3552 Assert(this->update_flags & update_gradients,
3553 ExcAccessToUninitializedField("update_gradients"));
3554 Assert(present_cell.is_initialized(), ExcNotReinited());
3555 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3556
3557 // get function values of dofs on this cell
3558 Vector<Number> dof_values(dofs_per_cell);
3559 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3561 dof_values.begin(),
3562 this->finite_element_output.shape_gradients,
3563 *fe,
3564 this->finite_element_output.shape_function_to_row_table,
3565 make_array_view(gradients.begin(), gradients.end()));
3566}
3567
3568
3569
3570template <int dim, int spacedim>
3571template <class InputVector>
3572void
3574 const InputVector & fe_function,
3577 gradients,
3578 const bool quadrature_points_fastest) const
3579{
3580 using Number = typename InputVector::value_type;
3581 // Size of indices must be a multiple of dofs_per_cell such that an integer
3582 // number of function values is generated in each point.
3583 Assert(indices.size() % dofs_per_cell == 0,
3584 ExcNotMultiple(indices.size(), dofs_per_cell));
3585 Assert(this->update_flags & update_gradients,
3586 ExcAccessToUninitializedField("update_gradients"));
3587
3588 boost::container::small_vector<Number, 200> dof_values(indices.size());
3589 for (unsigned int i = 0; i < indices.size(); ++i)
3590 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3592 dof_values.data(),
3593 this->finite_element_output.shape_gradients,
3594 *fe,
3595 this->finite_element_output.shape_function_to_row_table,
3596 make_array_view(gradients.begin(), gradients.end()),
3597 quadrature_points_fastest,
3598 indices.size() / dofs_per_cell);
3599}
3600
3601
3602
3603template <int dim, int spacedim>
3604template <class InputVector>
3605void
3607 const InputVector &fe_function,
3609 const
3610{
3611 using Number = typename InputVector::value_type;
3612 AssertDimension(fe->n_components(), 1);
3613 Assert(this->update_flags & update_hessians,
3614 ExcAccessToUninitializedField("update_hessians"));
3615 Assert(present_cell.is_initialized(), ExcNotReinited());
3616 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3617
3618 // get function values of dofs on this cell
3619 Vector<Number> dof_values(dofs_per_cell);
3620 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3622 this->finite_element_output.shape_hessians,
3623 hessians);
3624}
3625
3626
3627
3628template <int dim, int spacedim>
3629template <class InputVector>
3630void
3632 const InputVector & fe_function,
3635 const
3636{
3637 using Number = typename InputVector::value_type;
3638 Assert(this->update_flags & update_hessians,
3639 ExcAccessToUninitializedField("update_hessians"));
3640 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3641 AssertDimension(indices.size(), dofs_per_cell);
3642
3643 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3644 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3645 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3647 this->finite_element_output.shape_hessians,
3648 hessians);
3649}
3650
3651
3652
3653template <int dim, int spacedim>
3654template <class InputVector>
3655void
3657 const InputVector &fe_function,
3658 std::vector<
3660 & hessians,
3661 const bool quadrature_points_fastest) const
3662{
3663 using Number = typename InputVector::value_type;
3664 Assert(this->update_flags & update_hessians,
3665 ExcAccessToUninitializedField("update_hessians"));
3666 Assert(present_cell.is_initialized(), ExcNotReinited());
3667 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3668
3669 // get function values of dofs on this cell
3670 Vector<Number> dof_values(dofs_per_cell);
3671 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3673 dof_values.begin(),
3674 this->finite_element_output.shape_hessians,
3675 *fe,
3676 this->finite_element_output.shape_function_to_row_table,
3677 make_array_view(hessians.begin(), hessians.end()),
3678 quadrature_points_fastest);
3679}
3680
3681
3682
3683template <int dim, int spacedim>
3684template <class InputVector>
3685void
3687 const InputVector & fe_function,
3690 hessians,
3691 const bool quadrature_points_fastest) const
3692{
3693 using Number = typename InputVector::value_type;
3694 Assert(this->update_flags & update_hessians,
3695 ExcAccessToUninitializedField("update_hessians"));
3696 Assert(indices.size() % dofs_per_cell == 0,
3697 ExcNotMultiple(indices.size(), dofs_per_cell));
3698
3699 boost::container::small_vector<Number, 200> dof_values(indices.size());
3700 for (unsigned int i = 0; i < indices.size(); ++i)
3701 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3703 dof_values.data(),
3704 this->finite_element_output.shape_hessians,
3705 *fe,
3706 this->finite_element_output.shape_function_to_row_table,
3707 make_array_view(hessians.begin(), hessians.end()),
3708 quadrature_points_fastest,
3709 indices.size() / dofs_per_cell);
3710}
3711
3712
3713
3714template <int dim, int spacedim>
3715template <class InputVector>
3716void
3718 const InputVector & fe_function,
3719 std::vector<typename InputVector::value_type> &laplacians) const
3720{
3721 using Number = typename InputVector::value_type;
3722 Assert(this->update_flags & update_hessians,
3723 ExcAccessToUninitializedField("update_hessians"));
3724 AssertDimension(fe->n_components(), 1);
3725 Assert(present_cell.is_initialized(), ExcNotReinited());
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 this->finite_element_output.shape_hessians,
3733 laplacians);
3735
3736
3737
3738template <int dim, int spacedim>
3739template <class InputVector>
3740void
3742 const InputVector & fe_function,
3744 std::vector<typename InputVector::value_type> & laplacians) const
3745{
3746 using Number = typename InputVector::value_type;
3747 Assert(this->update_flags & update_hessians,
3748 ExcAccessToUninitializedField("update_hessians"));
3749 AssertDimension(fe->n_components(), 1);
3750 AssertDimension(indices.size(), dofs_per_cell);
3751
3752 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3753 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3754 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3755 internal::do_function_laplacians(dof_values.data(),
3756 this->finite_element_output.shape_hessians,
3757 laplacians);
3758}
3759
3760
3761
3762template <int dim, int spacedim>
3763template <class InputVector>
3764void
3766 const InputVector & fe_function,
3767 std::vector<Vector<typename InputVector::value_type>> &laplacians) const
3768{
3769 using Number = typename InputVector::value_type;
3770 Assert(present_cell.is_initialized(), ExcNotReinited());
3771 Assert(this->update_flags & update_hessians,
3772 ExcAccessToUninitializedField("update_hessians"));
3773 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3774
3775 // get function values of dofs on this cell
3776 Vector<Number> dof_values(dofs_per_cell);
3777 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3779 dof_values.begin(),
3780 this->finite_element_output.shape_hessians,
3781 *fe,
3782 this->finite_element_output.shape_function_to_row_table,
3783 laplacians);
3784}
3785
3786
3787
3788template <int dim, int spacedim>
3789template <class InputVector>
3790void
3792 const InputVector & fe_function,
3794 std::vector<Vector<typename InputVector::value_type>> &laplacians) const
3795{
3796 using Number = typename InputVector::value_type;
3797 // Size of indices must be a multiple of dofs_per_cell such that an integer
3798 // number of function values is generated in each point.
3799 Assert(indices.size() % dofs_per_cell == 0,
3800 ExcNotMultiple(indices.size(), dofs_per_cell));
3801 Assert(this->update_flags & update_hessians,
3802 ExcAccessToUninitializedField("update_hessians"));
3803
3804 boost::container::small_vector<Number, 200> dof_values(indices.size());
3805 for (unsigned int i = 0; i < indices.size(); ++i)
3806 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3808 dof_values.data(),
3809 this->finite_element_output.shape_hessians,
3810 *fe,
3811 this->finite_element_output.shape_function_to_row_table,
3812 laplacians,
3813 false,
3814 indices.size() / dofs_per_cell);
3815}
3816
3817
3818
3819template <int dim, int spacedim>
3820template <class InputVector>
3821void
3823 const InputVector & fe_function,
3825 std::vector<std::vector<typename InputVector::value_type>> &laplacians,
3826 const bool quadrature_points_fastest) const
3827{
3828 using Number = typename InputVector::value_type;
3829 Assert(indices.size() % dofs_per_cell == 0,
3830 ExcNotMultiple(indices.size(), dofs_per_cell));
3831 Assert(this->update_flags & update_hessians,
3832 ExcAccessToUninitializedField("update_hessians"));
3833
3834 boost::container::small_vector<Number, 200> dof_values(indices.size());
3835 for (unsigned int i = 0; i < indices.size(); ++i)
3836 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3838 dof_values.data(),
3839 this->finite_element_output.shape_hessians,
3840 *fe,
3841 this->finite_element_output.shape_function_to_row_table,
3842 laplacians,
3843 quadrature_points_fastest,
3844 indices.size() / dofs_per_cell);
3845}
3846
3848
3849template <int dim, int spacedim>
3850template <class InputVector>
3851void
3853 const InputVector &fe_function,
3855 &third_derivatives) const
3856{
3857 using Number = typename InputVector::value_type;
3858 AssertDimension(fe->n_components(), 1);
3859 Assert(this->update_flags & update_3rd_derivatives,
3860 ExcAccessToUninitializedField("update_3rd_derivatives"));
3861 Assert(present_cell.is_initialized(), ExcNotReinited());
3862 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3863
3864 // get function values of dofs on this cell
3865 Vector<Number> dof_values(dofs_per_cell);
3866 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3868 dof_values.begin(),
3869 this->finite_element_output.shape_3rd_derivatives,
3870 third_derivatives);
3871}
3872
3873
3874
3875template <int dim, int spacedim>
3876template <class InputVector>
3877void
3879 const InputVector & fe_function,
3882 &third_derivatives) const
3883{
3884 using Number = typename InputVector::value_type;
3885 Assert(this->update_flags & update_3rd_derivatives,
3886 ExcAccessToUninitializedField("update_3rd_derivatives"));
3887 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3888 AssertDimension(indices.size(), dofs_per_cell);
3889
3890 boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3891 for (unsigned int i = 0; i < dofs_per_cell; ++i)
3892 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3894 dof_values.data(),
3895 this->finite_element_output.shape_3rd_derivatives,
3896 third_derivatives);
3897}
3898
3899
3900
3901template <int dim, int spacedim>
3902template <class InputVector>
3903void
3905 const InputVector &fe_function,
3906 std::vector<
3908 & third_derivatives,
3909 const bool quadrature_points_fastest) const
3910{
3911 using Number = typename InputVector::value_type;
3912 Assert(this->update_flags & update_3rd_derivatives,
3913 ExcAccessToUninitializedField("update_3rd_derivatives"));
3914 Assert(present_cell.is_initialized(), ExcNotReinited());
3915 AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
3916
3917 // get function values of dofs on this cell
3918 Vector<Number> dof_values(dofs_per_cell);
3919 present_cell.get_interpolated_dof_values(fe_function, dof_values);
3921 dof_values.begin(),
3922 this->finite_element_output.shape_3rd_derivatives,
3923 *fe,
3924 this->finite_element_output.shape_function_to_row_table,
3925 make_array_view(third_derivatives.begin(), third_derivatives.end()),
3926 quadrature_points_fastest);
3927}
3928
3929
3930
3931template <int dim, int spacedim>
3932template <class InputVector>
3933void
3935 const InputVector & fe_function,
3938 third_derivatives,
3939 const bool quadrature_points_fastest) const
3940{
3941 using Number = typename InputVector::value_type;
3942 Assert(this->update_flags & update_3rd_derivatives,
3943 ExcAccessToUninitializedField("update_3rd_derivatives"));
3944 Assert(indices.size() % dofs_per_cell == 0,
3945 ExcNotMultiple(indices.size(), dofs_per_cell));
3946
3947 boost::container::small_vector<Number, 200> dof_values(indices.size());
3948 for (unsigned int i = 0; i < indices.size(); ++i)
3949 dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
3951 dof_values.data(),
3952 this->finite_element_output.shape_3rd_derivatives,
3953 *fe,
3954 this->finite_element_output.shape_function_to_row_table,
3955 make_array_view(third_derivatives.begin(), third_derivatives.end()),
3956 quadrature_points_fastest,
3957 indices.size() / dofs_per_cell);
3958}
3959
3960
3961
3962template <int dim, int spacedim>
3965{
3966 return present_cell;
3967}
3968
3969
3970
3971template <int dim, int spacedim>
3972const std::vector<Tensor<1, spacedim>> &
3974{
3975 Assert(this->update_flags & update_normal_vectors,
3977 "update_normal_vectors")));
3978
3979 return this->mapping_output.normal_vectors;
3980}
3981
3982
3983
3984template <int dim, int spacedim>
3985std::size_t
3987{
3988 return (sizeof(this->update_flags) +
3989 MemoryConsumption::memory_consumption(n_quadrature_points) +
3990 MemoryConsumption::memory_consumption(max_n_quadrature_points) +
3991 sizeof(cell_similarity) +
4000 MemoryConsumption::memory_consumption(finite_element_output));
4001}
4002
4003
4004
4005template <int dim, int spacedim>
4008 const UpdateFlags update_flags) const
4009{
4010 // first find out which objects need to be recomputed on each
4011 // cell we visit. this we have to ask the finite element and mapping.
4012 // elements are first since they might require update in mapping
4013 //
4014 // there is no need to iterate since mappings will never require
4015 // the finite element to compute something for them
4016 UpdateFlags flags = update_flags | fe->requires_update_flags(update_flags);
4017 flags |= mapping->requires_update_flags(flags);
4018
4019 return flags;
4020}
4021
4022
4023
4024template <int dim, int spacedim>
4025void
4027{
4028 // if there is no present cell, then we shouldn't be
4029 // connected via a signal to a triangulation
4030 Assert(present_cell.is_initialized(), ExcInternalError());
4031
4032 // so delete the present cell and
4033 // disconnect from the signal we have with
4034 // it
4035 tria_listener_refinement.disconnect();
4036 tria_listener_mesh_transform.disconnect();
4037 present_cell = {};
4038}
4039
4040
4041
4042template <int dim, int spacedim>
4043void
4046{
4047 if (present_cell.is_initialized())
4048 {
4049 if (&cell->get_triangulation() !=
4050 &present_cell
4051 .
4054 {
4055 // the triangulations for the previous cell and the current cell
4056 // do not match. disconnect from the previous triangulation and
4057 // connect to the current one; also invalidate the previous
4058 // cell because we shouldn't be comparing cells from different
4059 // triangulations
4060 invalidate_present_cell();
4061 tria_listener_refinement =
4062 cell->get_triangulation().signals.any_change.connect(
4063 [this]() { this->invalidate_present_cell(); });
4064 tria_listener_mesh_transform =
4065 cell->get_triangulation().signals.mesh_movement.connect(
4066 [this]() { this->invalidate_present_cell(); });
4067 }
4068 }
4069 else
4070 {
4071 // if this FEValues has never been set to any cell at all, then
4072 // at least subscribe to the triangulation to get notified of
4073 // changes
4074 tria_listener_refinement =
4075 cell->get_triangulation().signals.post_refinement.connect(
4076 [this]() { this->invalidate_present_cell(); });
4077 tria_listener_mesh_transform =
4078 cell->get_triangulation().signals.mesh_movement.connect(
4079 [this]() { this->invalidate_present_cell(); });
4080 }
4081}
4082
4083
4084
4085template <int dim, int spacedim>
4086inline void
4089{
4090 // Unfortunately, the detection of simple geometries with CellSimilarity is
4091 // sensitive to the first cell detected. When doing this with multiple
4092 // threads, each thread will get its own scratch data object with an
4093 // FEValues object in the implementation framework from late 2013, which is
4094 // initialized to the first cell the thread sees. As this number might
4095 // different between different runs (after all, the tasks are scheduled
4096 // dynamically onto threads), this slight deviation leads to difference in
4097 // roundoff errors that propagate through the program. Therefore, we need to
4098 // disable CellSimilarity in case there is more than one thread in the
4099 // problem. This will likely not affect many MPI test cases as there
4100 // multithreading is disabled on default, but in many other situations
4101 // because we rarely explicitly set the number of threads.
4102 //
4103 // TODO: Is it reasonable to introduce a flag "unsafe" in the constructor of
4104 // FEValues to re-enable this feature?
4106 {
4107 cell_similarity = CellSimilarity::none;
4108 return;
4109 }
4110
4111 // case that there has not been any cell before
4112 if (this->present_cell.is_initialized() == false)
4113 cell_similarity = CellSimilarity::none;
4114 else
4115 // in MappingQ, data can have been modified during the last call. Then, we
4116 // can't use that data on the new cell.
4117 if (cell_similarity == CellSimilarity::invalid_next_cell)
4118 cell_similarity = CellSimilarity::none;
4119 else
4120 cell_similarity =
4121 (cell->is_translation_of(
4122 static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
4123 &>(this->present_cell)) ?
4126
4127 if ((dim < spacedim) && (cell_similarity == CellSimilarity::translation))
4128 {
4129 if (static_cast<const typename Triangulation<dim, spacedim>::cell_iterator
4130 &>(this->present_cell)
4131 ->direction_flag() != cell->direction_flag())
4132 cell_similarity = CellSimilarity::inverted_translation;
4133 }
4134 // TODO: here, one could implement other checks for similarity, e.g. for
4135 // children of a parallelogram.
4136}
4137
4138
4139
4140template <int dim, int spacedim>
4143{
4144 return cell_similarity;
4145}
4146
4147
4148
4149template <int dim, int spacedim>
4151
4152
4153
4154template <int dim, int spacedim>
4156
4157/*------------------------------- FEValues -------------------------------*/
4158
4159template <int dim, int spacedim>
4161
4162
4163
4164template <int dim, int spacedim>
4167 const Quadrature<dim> & q,
4168 const UpdateFlags update_flags)
4169 : FEValuesBase<dim, spacedim>(q.size(),
4170 fe.n_dofs_per_cell(),
4172 mapping,
4173 fe)
4174 , quadrature(q)
4175{
4176 initialize(update_flags);
4177}
4178
4179
4180
4181template <int dim, int spacedim>
4184 const hp::QCollection<dim> & q,
4185 const UpdateFlags update_flags)
4186 : FEValues(mapping, fe, q[0], update_flags)
4187{
4188 AssertDimension(q.size(), 1);
4189}
4190
4191
4192
4193template <int dim, int spacedim>
4195 const Quadrature<dim> & q,
4196 const UpdateFlags update_flags)
4197 : FEValuesBase<dim, spacedim>(
4198 q.size(),
4199 fe.n_dofs_per_cell(),
4201 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4202 fe)
4203 , quadrature(q)
4204{
4205 initialize(update_flags);
4206}
4207
4208
4209
4210template <int dim, int spacedim>
4212 const hp::QCollection<dim> & q,
4213 const UpdateFlags update_flags)
4214 : FEValues(fe, q[0], update_flags)
4215{
4216 AssertDimension(q.size(), 1);
4217}
4218
4219
4220
4221template <int dim, int spacedim>
4222void
4224{
4225 // You can compute normal vectors to the cells only in the
4226 // codimension one case.
4227 if (dim != spacedim - 1)
4228 Assert((update_flags & update_normal_vectors) == false,
4229 ExcMessage("You can only pass the 'update_normal_vectors' "
4230 "flag to FEFaceValues or FESubfaceValues objects, "
4231 "but not to an FEValues object unless the "
4232 "triangulation it refers to is embedded in a higher "
4233 "dimensional space."));
4234
4235 const UpdateFlags flags = this->compute_update_flags(update_flags);
4236
4237 // initialize the base classes
4238 if (flags & update_mapping)
4239 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4240 this->finite_element_output.initialize(this->max_n_quadrature_points,
4241 *this->fe,
4242 flags);
4243
4244 // then get objects into which the FE and the Mapping can store
4245 // intermediate data used across calls to reinit. we can do this in parallel
4247 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4248 fe_get_data = Threads::new_task([&]() {
4249 return this->fe->get_data(flags,
4250 *this->mapping,
4251 quadrature,
4252 this->finite_element_output);
4253 });
4254
4256 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4257 mapping_get_data;
4258 if (flags & update_mapping)
4259 mapping_get_data = Threads::new_task(
4260 [&]() { return this->mapping->get_data(flags, quadrature); });
4261
4262 this->update_flags = flags;
4263
4264 // then collect answers from the two task above
4265 this->fe_data = std::move(fe_get_data.return_value());
4266 if (flags & update_mapping)
4267 this->mapping_data = std::move(mapping_get_data.return_value());
4268 else
4269 this->mapping_data =
4270 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4271}
4272
4273
4274
4275template <int dim, int spacedim>
4276void
4279{
4280 // Check that mapping and reference cell type are compatible:
4281 Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
4282 ExcMessage(
4283 "You are trying to call FEValues::reinit() with a cell of type " +
4284 cell->reference_cell().to_string() +
4285 " with a Mapping that is not compatible with it."));
4286
4287 // no FE in this cell, so no assertion
4288 // necessary here
4289 this->maybe_invalidate_previous_present_cell(cell);
4290 this->check_cell_similarity(cell);
4291
4292 this->present_cell = {cell};
4293
4294 // this was the part of the work that is dependent on the actual
4295 // data type of the iterator. now pass on to the function doing
4296 // the real work.
4297 do_reinit();
4298}
4299
4300
4301
4302template <int dim, int spacedim>
4303template <bool lda>
4304void
4307{
4308 // assert that the finite elements passed to the constructor and
4309 // used by the DoFHandler used by this cell, are the same
4310 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4311 static_cast<const FiniteElementData<dim> &>(cell->get_fe()),
4313
4314 // Check that mapping and reference cell type are compatible:
4315 Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
4316 ExcMessage(
4317 "You are trying to call FEValues::reinit() with a cell of type " +
4318 cell->reference_cell().to_string() +
4319 " with a Mapping that is not compatible with it."));
4320
4321 this->maybe_invalidate_previous_present_cell(cell);
4322 this->check_cell_similarity(cell);
4323
4324 this->present_cell = {cell};
4325
4326 // this was the part of the work that is dependent on the actual
4327 // data type of the iterator. now pass on to the function doing
4328 // the real work.
4329 do_reinit();
4330}
4331
4332
4333
4334template <int dim, int spacedim>
4335void
4337{
4338 // first call the mapping and let it generate the data
4339 // specific to the mapping. also let it inspect the
4340 // cell similarity flag and, if necessary, update
4341 // it
4342 if (this->update_flags & update_mapping)
4343 {
4344 this->cell_similarity =
4345 this->get_mapping().fill_fe_values(this->present_cell,
4346 this->cell_similarity,
4347 quadrature,
4348 *this->mapping_data,
4349 this->mapping_output);
4350 }
4351
4352 // then call the finite element and, with the data
4353 // already filled by the mapping, let it compute the
4354 // data for the mapped shape function values, gradients,
4355 // etc.
4356 this->get_fe().fill_fe_values(this->present_cell,
4357 this->cell_similarity,
4358 this->quadrature,
4359 this->get_mapping(),
4360 *this->mapping_data,
4361 this->mapping_output,
4362 *this->fe_data,
4363 this->finite_element_output);
4364}
4365
4366
4367
4368template <int dim, int spacedim>
4369std::size_t
4371{
4374}
4375
4376
4377/*------------------------------- FEFaceValuesBase --------------------------*/
4378
4379
4380template <int dim, int spacedim>
4382 const unsigned int dofs_per_cell,
4383 const UpdateFlags flags,
4384 const Mapping<dim, spacedim> & mapping,
4386 const Quadrature<dim - 1> & quadrature)
4387 : FEFaceValuesBase<dim, spacedim>(dofs_per_cell,
4388 flags,
4389 mapping,
4390 fe,
4391 hp::QCollection<dim - 1>(quadrature))
4392{}
4393
4394
4395
4396template <int dim, int spacedim>
4398 const unsigned int dofs_per_cell,
4399 const UpdateFlags,
4400 const Mapping<dim, spacedim> & mapping,
4402 const hp::QCollection<dim - 1> & quadrature)
4403 : FEValuesBase<dim, spacedim>(quadrature.max_n_quadrature_points(),
4404 dofs_per_cell,
4406 mapping,
4407 fe)
4408 , present_face_index(numbers::invalid_unsigned_int)
4409 , quadrature(quadrature)
4410{
4411 Assert(quadrature.size() == 1 ||
4412 quadrature.size() == fe.reference_cell().n_faces(),
4414}
4415
4416
4417
4418template <int dim, int spacedim>
4419const std::vector<Tensor<1, spacedim>> &
4421{
4422 Assert(this->update_flags & update_boundary_forms,
4424 "update_boundary_forms")));
4425 return this->mapping_output.boundary_forms;
4426}
4427
4428
4429
4430template <int dim, int spacedim>
4431std::size_t
4433{
4436}
4437
4438
4439/*------------------------------- FEFaceValues -------------------------------*/
4440
4441template <int dim, int spacedim>
4443
4444
4445
4446template <int dim, int spacedim>
4448
4449
4450
4451template <int dim, int spacedim>
4453 const Mapping<dim, spacedim> & mapping,
4455 const Quadrature<dim - 1> & quadrature,
4456 const UpdateFlags update_flags)
4457 : FEFaceValues<dim, spacedim>(mapping,
4458 fe,
4459 hp::QCollection<dim - 1>(quadrature),
4460 update_flags)
4461{}
4462
4463
4464
4465template <int dim, int spacedim>
4467 const Mapping<dim, spacedim> & mapping,
4469 const hp::QCollection<dim - 1> & quadrature,
4470 const UpdateFlags update_flags)
4471 : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
4472 update_flags,
4473 mapping,
4474 fe,
4475 quadrature)
4476{
4477 initialize(update_flags);
4478}
4479
4480
4481
4482template <int dim, int spacedim>
4485 const Quadrature<dim - 1> & quadrature,
4486 const UpdateFlags update_flags)
4487 : FEFaceValues<dim, spacedim>(fe,
4488 hp::QCollection<dim - 1>(quadrature),
4489 update_flags)
4490{}
4491
4492
4493
4494template <int dim, int spacedim>
4497 const hp::QCollection<dim - 1> & quadrature,
4498 const UpdateFlags update_flags)
4499 : FEFaceValuesBase<dim, spacedim>(
4500 fe.n_dofs_per_cell(),
4501 update_flags,
4502 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4503 fe,
4504 quadrature)
4505{
4506 initialize(update_flags);
4507}
4508
4509
4510
4511template <int dim, int spacedim>
4512void
4514{
4515 const UpdateFlags flags = this->compute_update_flags(update_flags);
4516
4517 // initialize the base classes
4518 if (flags & update_mapping)
4519 this->mapping_output.initialize(this->max_n_quadrature_points, flags);
4520 this->finite_element_output.initialize(this->max_n_quadrature_points,
4521 *this->fe,
4522 flags);
4523
4524 // then get objects into which the FE and the Mapping can store
4525 // intermediate data used across calls to reinit. this can be done in parallel
4526
4527 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> (
4528 FiniteElement<dim, spacedim>::*finite_element_get_face_data)(
4529 const UpdateFlags,
4530 const Mapping<dim, spacedim> &,
4533 spacedim>
4535
4536 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> (
4537 Mapping<dim, spacedim>::*mapping_get_face_data)(
4538 const UpdateFlags, const hp::QCollection<dim - 1> &) const =
4540
4541
4543 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4544 fe_get_data = Threads::new_task(finite_element_get_face_data,
4545 *this->fe,
4546 flags,
4547 *this->mapping,
4548 this->quadrature,
4549 this->finite_element_output);
4551 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4552 mapping_get_data;
4553 if (flags & update_mapping)
4554 mapping_get_data = Threads::new_task(mapping_get_face_data,
4555 *this->mapping,
4556 flags,
4557 this->quadrature);
4558
4559 this->update_flags = flags;
4560
4561 // then collect answers from the two task above
4562 this->fe_data = std::move(fe_get_data.return_value());
4563 if (flags & update_mapping)
4564 this->mapping_data = std::move(mapping_get_data.return_value());
4565 else
4566 this->mapping_data =
4567 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4568}
4569
4570
4571
4572template <int dim, int spacedim>
4573template <bool lda>
4574void
4577 const unsigned int face_no)
4578{
4579 // assert that the finite elements passed to the constructor and
4580 // used by the DoFHandler used by this cell, are the same
4581 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4582 static_cast<const FiniteElementData<dim> &>(
4583 cell->get_dof_handler().get_fe(cell->active_fe_index())),
4585
4587
4588 this->maybe_invalidate_previous_present_cell(cell);
4589 this->present_cell = {cell};
4590
4591 // this was the part of the work that is dependent on the actual
4592 // data type of the iterator. now pass on to the function doing
4593 // the real work.
4594 do_reinit(face_no);
4595}
4596
4597
4598
4599template <int dim, int spacedim>
4600template <bool lda>
4601void
4605{
4606 const auto face_n = cell->face_iterator_to_index(face);
4607 reinit(cell, face_n);
4608}
4609
4610
4611
4612template <int dim, int spacedim>
4613void
4616 const unsigned int face_no)
4617{
4619
4620 this->maybe_invalidate_previous_present_cell(cell);
4621 this->present_cell = {cell};
4622
4623 // this was the part of the work that is dependent on the actual
4624 // data type of the iterator. now pass on to the function doing
4625 // the real work.
4626 do_reinit(face_no);
4627}
4628
4629
4630
4631template <int dim, int spacedim>
4632void
4636{
4637 const auto face_n = cell->face_iterator_to_index(face);
4638 reinit(cell, face_n);
4639}
4640
4641
4642
4643template <int dim, int spacedim>
4644void
4645FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no)
4646{
4647 this->present_face_no = face_no;
4648
4649 // first of all, set the present_face_index (if available)
4651 this->present_cell;
4652 this->present_face_index = cell->face_index(face_no);
4653
4654 if (this->update_flags & update_mapping)
4655 {
4656 this->get_mapping().fill_fe_face_values(this->present_cell,
4657 face_no,
4658 this->quadrature,
4659 *this->mapping_data,
4660 this->mapping_output);
4661 }
4662
4663 this->get_fe().fill_fe_face_values(this->present_cell,
4664 face_no,
4665 this->quadrature,
4666 this->get_mapping(),
4667 *this->mapping_data,
4668 this->mapping_output,
4669 *this->fe_data,
4670 this->finite_element_output);
4671
4672 const_cast<unsigned int &>(this->n_quadrature_points) =
4673 this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size();
4674}
4675
4676
4677/* ---------------------------- FESubFaceValues ---------------------------- */
4678
4679
4680template <int dim, int spacedim>
4682
4683
4684
4685template <int dim, int spacedim>
4687
4688
4689
4690template <int dim, int spacedim>
4692 const Mapping<dim, spacedim> & mapping,
4694 const Quadrature<dim - 1> & quadrature,
4695 const UpdateFlags update_flags)
4696 : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
4697 update_flags,
4698 mapping,
4699 fe,
4700 quadrature)
4701{
4702 initialize(update_flags);
4703}
4704
4705
4706
4707template <int dim, int spacedim>
4709 const Mapping<dim, spacedim> & mapping,
4711 const hp::QCollection<dim - 1> & quadrature,
4712 const UpdateFlags update_flags)
4713 : FESubfaceValues(mapping, fe, quadrature[0], update_flags)
4714{
4715 AssertDimension(quadrature.size(), 1);
4716}
4717
4718
4719
4720template <int dim, int spacedim>
4723 const Quadrature<dim - 1> & quadrature,
4724 const UpdateFlags update_flags)
4725 : FEFaceValuesBase<dim, spacedim>(
4726 fe.n_dofs_per_cell(),
4727 update_flags,
4728 fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
4729 fe,
4730 quadrature)
4731{
4732 initialize(update_flags);
4733}
4734
4735
4736
4737template <int dim, int spacedim>
4740 const hp::QCollection<dim - 1> & quadrature,
4741 const UpdateFlags update_flags)
4742 : FESubfaceValues(fe, quadrature[0], update_flags)
4743{
4744 AssertDimension(quadrature.size(), 1);
4745}
4746
4747
4748
4749template <int dim, int spacedim>
4750void
4752{
4753 const UpdateFlags flags = this->compute_update_flags(update_flags);
4754
4755 // initialize the base classes
4756 if (flags & update_mapping)
4759 *this->fe,
4760 flags);
4761
4762 // then get objects into which the FE and the Mapping can store
4763 // intermediate data used across calls to reinit. this can be done
4764 // in parallel
4766 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
4767 fe_get_data =
4769 *this->fe,
4770 flags,
4771 *this->mapping,
4772 this->quadrature[0],
4773 this->finite_element_output);
4775 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
4776 mapping_get_data;
4777 if (flags & update_mapping)
4778 mapping_get_data =
4780 *this->mapping,
4781 flags,
4782 this->quadrature[0]);
4783
4784 this->update_flags = flags;
4785
4786 // then collect answers from the two task above
4787 this->fe_data = std::move(fe_get_data.return_value());
4788 if (flags & update_mapping)
4789 this->mapping_data = std::move(mapping_get_data.return_value());
4790 else
4791 this->mapping_data =
4792 std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
4793}
4794
4795
4796
4797template <int dim, int spacedim>
4798template <bool lda>
4799void
4802 const unsigned int face_no,
4803 const unsigned int subface_no)
4804{
4805 // assert that the finite elements passed to the constructor and
4806 // used by the DoFHandler used by this cell, are the same
4807 Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
4808 static_cast<const FiniteElementData<dim> &>(
4809 cell->get_dof_handler().get_fe(cell->active_fe_index())),
4812 // We would like to check for subface_no < cell->face(face_no)->n_children(),
4813 // but unfortunately the current function is also called for
4814 // faces without children (see tests/fe/mapping.cc). Therefore,
4815 // we must use following workaround of two separate assertions
4816 Assert(cell->face(face_no)->has_children() ||
4818 ExcIndexRange(subface_no,
4819 0,
4821 Assert(!cell->face(face_no)->has_children() ||
4822 subface_no < cell->face(face_no)->n_active_descendants(),
4823 ExcIndexRange(subface_no,
4824 0,
4825 cell->face(face_no)->n_active_descendants()));
4826 Assert(cell->has_children() == false,
4827 ExcMessage("You can't use subface data for cells that are "
4828 "already refined. Iterate over their children "
4829 "instead in these cases."));
4830
4832 this->present_cell = {cell};
4833
4834 // this was the part of the work that is dependent on the actual
4835 // data type of the iterator. now pass on to the function doing
4836 // the real work.
4837 do_reinit(face_no, subface_no);
4838}
4839
4840
4841
4842template <int dim, int spacedim>
4843template <bool lda>
4844void
4848 const typename Triangulation<dim, spacedim>::face_iterator &subface)
4849{
4850 reinit(cell,
4851 cell->face_iterator_to_index(face),
4852 face->child_iterator_to_index(subface));
4853}
4854
4855
4856
4857template <int dim, int spacedim>
4858void
4861 const unsigned int face_no,
4862 const unsigned int subface_no)
4863{
4865 // We would like to check for subface_no < cell->face(face_no)->n_children(),
4866 // but unfortunately the current function is also called for
4867 // faces without children for periodic faces, which have hanging nodes on
4868 // the other side (see include/deal.II/matrix_free/mapping_info.templates.h).
4869 AssertIndexRange(subface_no,
4870 (cell->has_periodic_neighbor(face_no) ?
4871 cell->periodic_neighbor(face_no)
4872 ->face(cell->periodic_neighbor_face_no(face_no))
4873 ->n_children() :
4874 cell->face(face_no)->n_children()));
4875
4877 this->present_cell = {cell};
4878
4879 // this was the part of the work that is dependent on the actual
4880 // data type of the iterator. now pass on to the function doing
4881 // the real work.
4882 do_reinit(face_no, subface_no);
4883}
4884
4885
4886
4887template <int dim, int spacedim>
4888void
4892 const typename Triangulation<dim, spacedim>::face_iterator &subface)
4893{
4894 reinit(cell,
4895 cell->face_iterator_to_index(face),
4896 face->child_iterator_to_index(subface));
4897}
4898
4899
4900
4901template <int dim, int spacedim>
4902void
4903FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no,
4904 const unsigned int subface_no)
4905{
4906 this->present_face_no = face_no;
4907
4908 // first of all, set the present_face_index (if available)
4910 this->present_cell;
4911
4912 if (!cell->face(face_no)->has_children())
4913 // no subfaces at all, so set present_face_index to this face rather
4914 // than any subface
4915 this->present_face_index = cell->face_index(face_no);
4916 else if (dim != 3)
4917 this->present_face_index = cell->face(face_no)->child_index(subface_no);
4918 else
4919 {
4920 // this is the same logic we use in cell->neighbor_child_on_subface(). See
4921 // there for an explanation of the different cases
4922 unsigned int subface_index = numbers::invalid_unsigned_int;
4923 switch (cell->subface_case(face_no))
4924 {
4928 subface_index = cell->face(face_no)->child_index(subface_no);
4929 break;
4932 subface_index = cell->face(face_no)
4933 ->child(subface_no / 2)
4934 ->child_index(subface_no % 2);
4935 break;
4938 switch (subface_no)
4939 {
4940 case 0:
4941 case 1:
4942 subface_index =
4943 cell->face(face_no)->child(0)->child_index(subface_no);
4944 break;
4945 case 2:
4946 subface_index = cell->face(face_no)->child_index(1);
4947 break;
4948 default:
4949 Assert(false, ExcInternalError());
4950 }
4951 break;
4954 switch (subface_no)
4955 {
4956 case 0:
4957 subface_index = cell->face(face_no)->child_index(0);
4958 break;
4959 case 1:
4960 case 2:
4961 subface_index =
4962 cell->face(face_no)->child(1)->child_index(subface_no - 1);
4963 break;
4964 default:
4965 Assert(false, ExcInternalError());
4966 }
4967 break;
4968 default:
4969 Assert(false, ExcInternalError());
4970 break;
4971 }
4972 Assert(subface_index != numbers::invalid_unsigned_int,
4974 this->present_face_index = subface_index;
4975 }
4976
4977 // now ask the mapping and the finite element to do the actual work
4978 if (this->update_flags & update_mapping)
4979 {
4981 face_no,
4982 subface_no,
4983 this->quadrature[0],
4984 *this->mapping_data,
4985 this->mapping_output);
4986 }
4987
4989 face_no,
4990 subface_no,
4991 this->quadrature[0],
4992 this->get_mapping(),
4993 *this->mapping_data,
4994 this->mapping_output,
4995 *this->fe_data,
4996 this->finite_element_output);
4997}
4998
4999
5000/*------------------------------- Explicit Instantiations -------------*/
5001#define SPLIT_INSTANTIATIONS_COUNT 6
5002#ifndef SPLIT_INSTANTIATIONS_INDEX
5003# define SPLIT_INSTANTIATIONS_INDEX 0
5004#endif
5005#include "fe_values.inst"
5006
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:699
iterator begin() const
Definition: array_view.h:585
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:699
value_type * data() const noexcept
Definition: array_view.h:553
std::size_t size() const
Definition: array_view.h:576
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
void get_interpolated_dof_values(const InputVector &values, Vector< number > &interpolated_values, const unsigned int fe_index=DoFHandler< dimension_, space_dimension_ >::invalid_fe_index) const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
unsigned int present_face_index
Definition: fe_values.h:4279
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:4284
void initialize(const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no)
void do_reinit(const unsigned int face_no)
FEFaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature, const UpdateFlags update_flags)
FESubfaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &face_quadrature, const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no, const unsigned int subface_no)
void initialize(const UpdateFlags update_flags)
void do_reinit(const unsigned int face_no, const unsigned int subface_no)
void get_interpolated_dof_values(const VectorType &in, Vector< typename VectorType::value_type > &out) const
Definition: fe_values.cc:2652
types::global_dof_index n_dofs_for_dof_handler() const
Definition: fe_values.cc:2638
CellSimilarity::Similarity cell_similarity
Definition: fe_values.h:4003
CellIteratorContainer present_cell
Definition: fe_values.h:3894
::internal::FEValuesViews::Cache< dim, spacedim > fe_values_views_cache
Definition: fe_values.h:4018
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3356
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:2863
virtual ~FEValuesBase() override
Definition: fe_values.cc:2887
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition: fe_values.h:3939
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:4087
UpdateFlags update_flags
Definition: fe_values.h:3985
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
Definition: fe_values.h:3963
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:4142
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Definition: fe_values.cc:3973
std::size_t memory_consumption() const
Definition: fe_values.cc:3986
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:3964
void get_function_laplacians(const InputVector &fe_function, std::vector< typename InputVector::value_type > &laplacians) const
Definition: fe_values.cc:3717
::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > mapping_output
Definition: fe_values.h:3954
UpdateFlags compute_update_flags(const UpdateFlags update_flags) const
Definition: fe_values.cc:4007
std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > mapping_data
Definition: fe_values.h:3947
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
Definition: fe_values.cc:3495
void invalidate_present_cell()
Definition: fe_values.cc:4026
::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > finite_element_output
Definition: fe_values.h:3979
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type > > &hessians) const
Definition: fe_values.cc:3606
std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > fe_data
Definition: fe_values.h:3971
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:4044
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:3852
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
void get_function_values(const InputVector &fe_function, std::vector< solution_value_type< typename InputVector::value_type > > &values) const
Definition: fe_values.cc:1547
typename ProductType< Number, value_type >::type solution_value_type
Definition: fe_values.h:1505
typename ProductType< Number, divergence_type >::type solution_divergence_type
Definition: fe_values.h:1515
typename ProductType< Number, value_type >::type solution_value_type
Definition: fe_values.h:1841
typename ProductType< Number, divergence_type >::type solution_divergence_type
Definition: fe_values.h:1851
typename ProductType< Number, gradient_type >::type solution_gradient_type
Definition: fe_values.h:1861
typename ProductType< Number, third_derivative_type >::type solution_third_derivative_type
Definition: fe_values.h:811
typename ProductType< Number, divergence_type >::type solution_divergence_type
Definition: fe_values.h:772
typename ProductType< Number, hessian_type >::type solution_hessian_type
Definition: fe_values.h:801
typename ProductType< Number, symmetric_gradient_type >::type solution_symmetric_gradient_type
Definition: fe_values.h:762
typename ProductType< Number, gradient_type >::type solution_gradient_type
Definition: fe_values.h:752
typename ProductType< Number, value_type >::type solution_value_type
Definition: fe_values.h:742
const unsigned int first_vector_component
Definition: fe_values.h:1441
typename ProductType< Number, curl_type >::type solution_curl_type
Definition: fe_values.h:791
std::vector< ShapeFunctionData > shape_function_data
Definition: fe_values.h:1446
typename ProductType< Number, value_type >::type solution_laplacian_type
Definition: fe_values.h:782
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
void do_reinit()
void initialize(const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
std::size_t memory_consumption() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
virtual std::unique_ptr< InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const =0
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
unsigned int n_nonzero_components(const unsigned int i) const
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
bool is_element(const size_type index) const
Definition: index_set.h:1767
signed int value_type
Definition: index_set.h:99
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const =0
static unsigned int n_threads()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
friend class SymmetricTensor
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:503
friend class Tensor
Definition: tensor.h:888
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Triangulation< dim, spacedim > & get_triangulation()
Signals signals
Definition: tria.h:2448
Definition: vector.h:109
unsigned int size() const
Definition: collection.h:264
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:2788
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
Definition: fe_values.cc:2704
std::vector< Tensor< 1, spacedim > > boundary_forms
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_grads
@ update_hessians
Second derivatives of shape functions.
@ update_jacobian_3rd_derivatives
@ update_values
Shape function values.
@ update_jacobian_grads
Gradient of volume element.
@ update_normal_vectors
Normal vectors.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
@ update_mapping
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ update_default
No update.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
static ::ExceptionBase & ExcNotMultiple(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotReinited()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition: tria.h:1355
iterator end()
iterator begin()
Task< RT > new_task(const std::function< RT()> &function)