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