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