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