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