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