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