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