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