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