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