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