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