Reference documentation for deal.II version Git 2618e0f 2017-11-23 17:25:26 +0100
fe_values.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2017 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 at
12 // the top level of the deal.II distribution.
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/differentiation/ad.h>
23 #include <deal.II/lac/vector.h>
24 #include <deal.II/lac/block_vector.h>
25 #include <deal.II/lac/la_vector.h>
26 #include <deal.II/lac/la_parallel_vector.h>
27 #include <deal.II/lac/la_parallel_block_vector.h>
28 #include <deal.II/lac/vector_element_access.h>
29 #include <deal.II/lac/petsc_parallel_vector.h>
30 #include <deal.II/lac/petsc_parallel_block_vector.h>
31 #include <deal.II/lac/trilinos_vector.h>
32 #include <deal.II/lac/trilinos_parallel_block_vector.h>
33 #include <deal.II/grid/tria_iterator.h>
34 #include <deal.II/grid/tria_accessor.h>
35 #include <deal.II/grid/tria_boundary.h>
36 #include <deal.II/dofs/dof_accessor.h>
37 #include <deal.II/fe/mapping_q1.h>
38 #include <deal.II/fe/fe_values.h>
39 #include <deal.II/fe/fe.h>
40 
41 #include <iomanip>
42 #include <memory>
43 #include <type_traits>
44 
45 
46 #include <boost/container/small_vector.hpp>
47 
48 DEAL_II_NAMESPACE_OPEN
49 
50 
51 namespace
52 {
53  template <class VectorType>
54  typename VectorType::value_type
55  get_vector_element (const VectorType &vector,
56  const types::global_dof_index cell_number)
57  {
58  return internal::ElementAccess<VectorType>::get(vector,cell_number);
59  }
60 
61 
63  get_vector_element (const IndexSet &is,
64  const types::global_dof_index cell_number)
65  {
66  return (is.is_element(cell_number) ? 1 : 0);
67  }
68 }
69 
70 
71 namespace
72 {
73  template <int dim, int spacedim>
74  inline
75  std::vector<unsigned int>
76  make_shape_function_to_row_table (const FiniteElement<dim,spacedim> &fe)
77  {
78  std::vector<unsigned int> shape_function_to_row_table (fe.dofs_per_cell * fe.n_components(),
80  unsigned int row = 0;
81  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
82  {
83  // loop over all components that are nonzero for this particular
84  // shape function. if a component is zero then we leave the
85  // value in the table unchanged (at the invalid value)
86  // otherwise it is mapped to the next free entry
87  unsigned int nth_nonzero_component = 0;
88  for (unsigned int c=0; c<fe.n_components(); ++c)
89  if (fe.get_nonzero_components(i)[c] == true)
90  {
91  shape_function_to_row_table[i*fe.n_components()+c] = row + nth_nonzero_component;
92  ++nth_nonzero_component;
93  }
94  row += fe.n_nonzero_components (i);
95  }
96 
97  return shape_function_to_row_table;
98  }
99 }
100 
101 
102 
103 namespace FEValuesViews
104 {
105  template <int dim, int spacedim>
107  const unsigned int component)
108  :
109  fe_values (&fe_values),
110  component (component),
111  shape_function_data (this->fe_values->fe->dofs_per_cell)
112  {
113  const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
114  Assert (component < fe.n_components(),
115  ExcIndexRange(component, 0, fe.n_components()));
116 
117 //TODO: we'd like to use the fields with the same name as these
118 // variables from FEValuesBase, but they aren't initialized yet
119 // at the time we get here, so re-create it all
120  const std::vector<unsigned int> shape_function_to_row_table
121  = make_shape_function_to_row_table (fe);
122 
123  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
124  {
125  const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
126 
127  if (is_primitive == true)
128  shape_function_data[i].is_nonzero_shape_function_component
129  = (component ==
130  fe.system_to_component_index(i).first);
131  else
132  shape_function_data[i].is_nonzero_shape_function_component
134  == true);
135 
136  if (shape_function_data[i].is_nonzero_shape_function_component == true)
137  shape_function_data[i].row_index
138  = shape_function_to_row_table[i*fe.n_components()+component];
139  else
141  }
142  }
143 
144 
145 
146  template <int dim, int spacedim>
148  :
149  fe_values (nullptr),
150  component (numbers::invalid_unsigned_int)
151  {}
152 
153 
154  template <int dim, int spacedim>
157  {
158  // we shouldn't be copying these objects
159  Assert (false, ExcInternalError());
160  return *this;
161  }
162 
163 
164 
165  template <int dim, int spacedim>
167  const unsigned int first_vector_component)
168  :
169  fe_values (&fe_values),
170  first_vector_component (first_vector_component),
171  shape_function_data (this->fe_values->fe->dofs_per_cell)
172  {
173  const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
174  Assert (first_vector_component+spacedim-1 < fe.n_components(),
175  ExcIndexRange(first_vector_component+spacedim-1, 0,
176  fe.n_components()));
177 
178 //TODO: we'd like to use the fields with the same name as these
179 // variables from FEValuesBase, but they aren't initialized yet
180 // at the time we get here, so re-create it all
181  const std::vector<unsigned int> shape_function_to_row_table
182  = make_shape_function_to_row_table (fe);
183 
184  for (unsigned int d=0; d<spacedim; ++d)
185  {
186  const unsigned int component = first_vector_component + d;
187 
188  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
189  {
190  const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
191 
192  if (is_primitive == true)
193  shape_function_data[i].is_nonzero_shape_function_component[d]
194  = (component ==
195  fe.system_to_component_index(i).first);
196  else
197  shape_function_data[i].is_nonzero_shape_function_component[d]
198  = (fe.get_nonzero_components(i)[component]
199  == true);
200 
201  if (shape_function_data[i].is_nonzero_shape_function_component[d]
202  == true)
203  shape_function_data[i].row_index[d]
204  = shape_function_to_row_table[i*fe.n_components()+component];
205  else
206  shape_function_data[i].row_index[d]
208  }
209  }
210 
211  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
212  {
213  unsigned int n_nonzero_components = 0;
214  for (unsigned int d=0; d<spacedim; ++d)
215  if (shape_function_data[i].is_nonzero_shape_function_component[d]
216  == true)
217  ++n_nonzero_components;
218 
219  if (n_nonzero_components == 0)
220  shape_function_data[i].single_nonzero_component = -2;
221  else if (n_nonzero_components > 1)
222  shape_function_data[i].single_nonzero_component = -1;
223  else
224  {
225  for (unsigned int d=0; d<spacedim; ++d)
226  if (shape_function_data[i].is_nonzero_shape_function_component[d]
227  == true)
228  {
229  shape_function_data[i].single_nonzero_component
230  = shape_function_data[i].row_index[d];
231  shape_function_data[i].single_nonzero_component_index
232  = d;
233  break;
234  }
235  }
236  }
237  }
238 
239 
240  template <int dim, int spacedim>
242  :
243  fe_values (nullptr),
244  first_vector_component (numbers::invalid_unsigned_int)
245  {}
246 
247 
248 
249  template <int dim, int spacedim>
252  {
253  // we shouldn't be copying these objects
254  Assert (false, ExcInternalError());
255  return *this;
256  }
257 
258 
259  template <int dim, int spacedim>
262  const unsigned int first_tensor_component)
263  :
264  fe_values(&fe_values),
265  first_tensor_component(first_tensor_component),
266  shape_function_data(this->fe_values->fe->dofs_per_cell)
267  {
268  const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
269  Assert(first_tensor_component + (dim*dim+dim)/2 - 1
270  <
271  fe.n_components(),
272  ExcIndexRange(first_tensor_component +
274  0,
275  fe.n_components()));
276 //TODO: we'd like to use the fields with the same name as these
277 // variables from FEValuesBase, but they aren't initialized yet
278 // at the time we get here, so re-create it all
279  const std::vector<unsigned int> shape_function_to_row_table
280  = make_shape_function_to_row_table (fe);
281 
282  for (unsigned int d = 0; d < ::SymmetricTensor<2,dim>::n_independent_components; ++d)
283  {
284  const unsigned int component = first_tensor_component + d;
285 
286  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
287  {
288  const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
289 
290  if (is_primitive == true)
291  shape_function_data[i].is_nonzero_shape_function_component[d]
292  = (component ==
293  fe.system_to_component_index(i).first);
294  else
295  shape_function_data[i].is_nonzero_shape_function_component[d]
296  = (fe.get_nonzero_components(i)[component]
297  == true);
298 
299  if (shape_function_data[i].is_nonzero_shape_function_component[d]
300  == true)
301  shape_function_data[i].row_index[d]
302  = shape_function_to_row_table[i*fe.n_components()+component];
303  else
304  shape_function_data[i].row_index[d]
306  }
307  }
308 
309  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
310  {
311  unsigned int n_nonzero_components = 0;
312  for (unsigned int d = 0; d < ::SymmetricTensor<2,dim>::n_independent_components; ++d)
313  if (shape_function_data[i].is_nonzero_shape_function_component[d]
314  == true)
315  ++n_nonzero_components;
316 
317  if (n_nonzero_components == 0)
318  shape_function_data[i].single_nonzero_component = -2;
319  else if (n_nonzero_components > 1)
320  shape_function_data[i].single_nonzero_component = -1;
321  else
322  {
323  for (unsigned int d = 0; d < ::SymmetricTensor<2,dim>::n_independent_components; ++d)
324  if (shape_function_data[i].is_nonzero_shape_function_component[d]
325  == true)
326  {
327  shape_function_data[i].single_nonzero_component
328  = shape_function_data[i].row_index[d];
329  shape_function_data[i].single_nonzero_component_index
330  = d;
331  break;
332  }
333  }
334  }
335  }
336 
337 
338 
339  template <int dim, int spacedim>
341  :
342  fe_values(nullptr),
343  first_tensor_component(numbers::invalid_unsigned_int)
344  {}
345 
346 
347 
348  template <int dim, int spacedim>
351  {
352  // we shouldn't be copying these objects
353  Assert(false, ExcInternalError());
354  return *this;
355  }
356 
357 
358  template <int dim, int spacedim>
361  const unsigned int first_tensor_component)
362  :
363  fe_values(&fe_values),
364  first_tensor_component(first_tensor_component),
365  shape_function_data(this->fe_values->fe->dofs_per_cell)
366  {
367  const FiniteElement<dim, spacedim> &fe = *this->fe_values->fe;
368  Assert(first_tensor_component + dim*dim - 1
369  <
370  fe.n_components(),
371  ExcIndexRange(first_tensor_component +
372  dim*dim - 1,
373  0,
374  fe.n_components()));
375 //TODO: we'd like to use the fields with the same name as these
376 // variables from FEValuesBase, but they aren't initialized yet
377 // at the time we get here, so re-create it all
378  const std::vector<unsigned int> shape_function_to_row_table
379  = make_shape_function_to_row_table (fe);
380 
381  for (unsigned int d = 0; d < dim*dim; ++d)
382  {
383  const unsigned int component = first_tensor_component + d;
384 
385  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
386  {
387  const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
388 
389  if (is_primitive == true)
390  shape_function_data[i].is_nonzero_shape_function_component[d]
391  = (component ==
392  fe.system_to_component_index(i).first);
393  else
394  shape_function_data[i].is_nonzero_shape_function_component[d]
395  = (fe.get_nonzero_components(i)[component]
396  == true);
397 
398  if (shape_function_data[i].is_nonzero_shape_function_component[d]
399  == true)
400  shape_function_data[i].row_index[d]
401  = shape_function_to_row_table[i*fe.n_components()+component];
402  else
403  shape_function_data[i].row_index[d]
405  }
406  }
407 
408  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
409  {
410  unsigned int n_nonzero_components = 0;
411  for (unsigned int d = 0; d < dim*dim; ++d)
412  if (shape_function_data[i].is_nonzero_shape_function_component[d]
413  == true)
414  ++n_nonzero_components;
415 
416  if (n_nonzero_components == 0)
417  shape_function_data[i].single_nonzero_component = -2;
418  else if (n_nonzero_components > 1)
419  shape_function_data[i].single_nonzero_component = -1;
420  else
421  {
422  for (unsigned int d = 0; d < dim*dim; ++d)
423  if (shape_function_data[i].is_nonzero_shape_function_component[d]
424  == true)
425  {
426  shape_function_data[i].single_nonzero_component
427  = shape_function_data[i].row_index[d];
428  shape_function_data[i].single_nonzero_component_index
429  = d;
430  break;
431  }
432  }
433  }
434  }
435 
436 
437 
438  template <int dim, int spacedim>
440  :
441  fe_values(nullptr),
442  first_tensor_component(numbers::invalid_unsigned_int)
443  {}
444 
445 
446 
447  template <int dim, int spacedim>
450  {
451  // we shouldn't be copying these objects
452  Assert(false, ExcInternalError());
453  return *this;
454  }
455 
456 
457  namespace internal
458  {
459  // Given values of degrees of freedom, evaluate the
460  // values/gradients/... at quadrature points
461 
462  // ------------------------- scalar functions --------------------------
463  template <int dim, int spacedim, typename Number>
464  void
465  do_function_values (const ArrayView<Number> &dof_values,
466  const Table<2,double> &shape_values,
467  const std::vector<typename Scalar<dim,spacedim>::ShapeFunctionData> &shape_function_data,
468  std::vector<typename ProductType<Number,double>::type> &values)
469  {
470  const unsigned int dofs_per_cell = dof_values.size();
471  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
472  shape_values.n_cols() : values.size();
473  AssertDimension (values.size(), n_quadrature_points);
474 
475  std::fill (values.begin(), values.end(),
477 
478  for (unsigned int shape_function=0;
479  shape_function<dofs_per_cell; ++shape_function)
480  if (shape_function_data[shape_function].is_nonzero_shape_function_component)
481  {
482  const Number &value = dof_values[shape_function];
483  // For auto-differentiable numbers, the fact that a DoF value is zero
484  // does not imply that its derivatives are zero as well. So we
485  // can't filter by value for these number types.
487  if (value == ::internal::NumberType<Number>::value(0.0) )
488  continue;
489 
490  const double *shape_value_ptr =
491  &shape_values(shape_function_data[shape_function].row_index, 0);
492  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
493  values[q_point] += value * (*shape_value_ptr++);
494  }
495  }
496 
497 
498 
499  // same code for gradient and Hessian, template argument 'order' to give
500  // the order of the derivative (= rank of gradient/Hessian tensor)
501  template <int order, int dim, int spacedim, typename Number>
502  void
503  do_function_derivatives (const ArrayView<Number> &dof_values,
504  const Table<2,::Tensor<order,spacedim> > &shape_derivatives,
505  const std::vector<typename Scalar<dim,spacedim>::ShapeFunctionData> &shape_function_data,
506  std::vector<typename ProductType<Number,::Tensor<order,spacedim> >::type> &derivatives)
507  {
508  const unsigned int dofs_per_cell = dof_values.size();
509  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
510  shape_derivatives[0].size() : derivatives.size();
511  AssertDimension (derivatives.size(), n_quadrature_points);
512 
513  std::fill (derivatives.begin(), derivatives.end(),
515 
516  for (unsigned int shape_function=0;
517  shape_function<dofs_per_cell; ++shape_function)
518  if (shape_function_data[shape_function].is_nonzero_shape_function_component)
519  {
520  const Number &value = dof_values[shape_function];
521  // For auto-differentiable numbers, the fact that a DoF value is zero
522  // does not imply that its derivatives are zero as well. So we
523  // can't filter by value for these number types.
525  if (value == ::internal::NumberType<Number>::value(0.0) )
526  continue;
527 
528  const ::Tensor<order,spacedim> *shape_derivative_ptr =
529  &shape_derivatives[shape_function_data[shape_function].row_index][0];
530  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
531  derivatives[q_point] += value *
532  ::Tensor<order,spacedim>(*shape_derivative_ptr++);
533  }
534  }
535 
536 
537 
538  template <int dim, int spacedim, typename Number>
539  void
540  do_function_laplacians (const ArrayView<Number> &dof_values,
541  const Table<2,::Tensor<2,spacedim> > &shape_hessians,
542  const std::vector<typename Scalar<dim,spacedim>::ShapeFunctionData> &shape_function_data,
543  std::vector<typename Scalar<dim,spacedim>::template OutputType<Number>::laplacian_type> &laplacians)
544  {
545  const unsigned int dofs_per_cell = dof_values.size();
546  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
547  shape_hessians[0].size() : laplacians.size();
548  AssertDimension (laplacians.size(), n_quadrature_points);
549 
550  std::fill (laplacians.begin(), laplacians.end(),
551  typename Scalar<dim,spacedim>::template OutputType<Number>::laplacian_type());
552 
553  for (unsigned int shape_function=0;
554  shape_function<dofs_per_cell; ++shape_function)
555  if (shape_function_data[shape_function].is_nonzero_shape_function_component)
556  {
557  const Number &value = dof_values[shape_function];
558  // For auto-differentiable numbers, the fact that a DoF value is zero
559  // does not imply that its derivatives are zero as well. So we
560  // can't filter by value for these number types.
562  if (value == ::internal::NumberType<Number>::value(0.0))
563  continue;
564 
565  const ::Tensor<2,spacedim> *shape_hessian_ptr =
566  &shape_hessians[shape_function_data[shape_function].row_index][0];
567  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
568  laplacians[q_point] += value * trace(*shape_hessian_ptr++);
569  }
570  }
571 
572 
573 
574  // ----------------------------- vector part ---------------------------
575 
576  template <int dim, int spacedim, typename Number>
577  void do_function_values (const ArrayView<Number> &dof_values,
578  const Table<2,double> &shape_values,
579  const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
580  std::vector<typename ProductType<Number,::Tensor<1,spacedim> >::type> &values)
581  {
582  const unsigned int dofs_per_cell = dof_values.size();
583  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
584  shape_values.n_cols() : values.size();
585  AssertDimension (values.size(), n_quadrature_points);
586 
587  std::fill (values.begin(), values.end(),
588  typename ProductType<Number,::Tensor<1,spacedim> >::type());
589 
590  for (unsigned int shape_function=0;
591  shape_function<dofs_per_cell; ++shape_function)
592  {
593  const int snc = shape_function_data[shape_function].single_nonzero_component;
594 
595  if (snc == -2)
596  // shape function is zero for the selected components
597  continue;
598 
599  const Number &value = dof_values[shape_function];
600  // For auto-differentiable numbers, the fact that a DoF value is zero
601  // does not imply that its derivatives are zero as well. So we
602  // can't filter by value for these number types.
604  if (value == ::internal::NumberType<Number>::value(0.0))
605  continue;
606 
607  if (snc != -1)
608  {
609  const unsigned int comp =
610  shape_function_data[shape_function].single_nonzero_component_index;
611  const double *shape_value_ptr = &shape_values(snc,0);
612  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
613  values[q_point][comp] += value * (*shape_value_ptr++);
614  }
615  else
616  for (unsigned int d=0; d<spacedim; ++d)
617  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
618  {
619  const double *shape_value_ptr =
620  &shape_values(shape_function_data[shape_function].row_index[d],0);
621  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
622  values[q_point][d] += value * (*shape_value_ptr++);
623  }
624  }
625  }
626 
627 
628 
629  template <int order, int dim, int spacedim, typename Number>
630  void
631  do_function_derivatives (const ArrayView<Number> &dof_values,
632  const Table<2,::Tensor<order,spacedim> > &shape_derivatives,
633  const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
634  std::vector<typename ProductType<Number,::Tensor<order+1,spacedim> >::type> &derivatives)
635  {
636  const unsigned int dofs_per_cell = dof_values.size();
637  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
638  shape_derivatives[0].size() : derivatives.size();
639  AssertDimension (derivatives.size(), n_quadrature_points);
640 
641  std::fill (derivatives.begin(), derivatives.end(),
643 
644  for (unsigned int shape_function=0;
645  shape_function<dofs_per_cell; ++shape_function)
646  {
647  const int snc = shape_function_data[shape_function].single_nonzero_component;
648 
649  if (snc == -2)
650  // shape function is zero for the selected components
651  continue;
652 
653  const Number &value = dof_values[shape_function];
654  // For auto-differentiable numbers, the fact that a DoF value is zero
655  // does not imply that its derivatives are zero as well. So we
656  // can't filter by value for these number types.
658  if (value == ::internal::NumberType<Number>::value(0.0))
659  continue;
660 
661  if (snc != -1)
662  {
663  const unsigned int comp =
664  shape_function_data[shape_function].single_nonzero_component_index;
665  const ::Tensor<order,spacedim> *shape_derivative_ptr =
666  &shape_derivatives[snc][0];
667  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
668  derivatives[q_point][comp] += value *
669  ::Tensor<order,spacedim>(*shape_derivative_ptr++);
670  }
671  else
672  for (unsigned int d=0; d<spacedim; ++d)
673  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
674  {
675  const ::Tensor<order,spacedim> *shape_derivative_ptr =
676  &shape_derivatives[shape_function_data[shape_function].
677  row_index[d]][0];
678  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
679  derivatives[q_point][d] += value *
680  ::Tensor<order,spacedim>(*shape_derivative_ptr++);
681  }
682  }
683  }
684 
685 
686 
687  template <int dim, int spacedim, typename Number>
688  void
689  do_function_symmetric_gradients (const ArrayView<Number> &dof_values,
690  const Table<2,::Tensor<1,spacedim> > &shape_gradients,
691  const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
692  std::vector<typename ProductType<Number,::SymmetricTensor<2,spacedim> >::type> &symmetric_gradients)
693  {
694  const unsigned int dofs_per_cell = dof_values.size();
695  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
696  shape_gradients[0].size() : symmetric_gradients.size();
697  AssertDimension (symmetric_gradients.size(), n_quadrature_points);
698 
699  std::fill (symmetric_gradients.begin(), symmetric_gradients.end(),
701 
702  for (unsigned int shape_function=0;
703  shape_function<dofs_per_cell; ++shape_function)
704  {
705  const int snc = shape_function_data[shape_function].single_nonzero_component;
706 
707  if (snc == -2)
708  // shape function is zero for the selected components
709  continue;
710 
711  const Number &value = dof_values[shape_function];
712  // For auto-differentiable numbers, the fact that a DoF value is zero
713  // does not imply that its derivatives are zero as well. So we
714  // can't filter by value for these number types.
716  if (value == ::internal::NumberType<Number>::value(0.0))
717  continue;
718 
719  if (snc != -1)
720  {
721  const unsigned int comp =
722  shape_function_data[shape_function].single_nonzero_component_index;
723  const ::Tensor<1,spacedim> *shape_gradient_ptr =
724  &shape_gradients[snc][0];
725  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
726  symmetric_gradients[q_point] += value *
727  ::SymmetricTensor<2,spacedim> (symmetrize_single_row(comp, *shape_gradient_ptr++));
728  }
729  else
730  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
731  {
732  typename ProductType<Number,::Tensor<2,spacedim> >::type grad;
733  for (unsigned int d=0; d<spacedim; ++d)
734  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
735  grad[d] = value *
736  shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
737  symmetric_gradients[q_point] += symmetrize(grad);
738  }
739  }
740  }
741 
742 
743 
744  template <int dim, int spacedim, typename Number>
745  void
746  do_function_divergences (const ArrayView<Number> &dof_values,
747  const Table<2,::Tensor<1,spacedim> > &shape_gradients,
748  const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
749  std::vector<typename Vector<dim,spacedim>::template OutputType<Number>::divergence_type> &divergences)
750  {
751  const unsigned int dofs_per_cell = dof_values.size();
752  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
753  shape_gradients[0].size() : divergences.size();
754  AssertDimension (divergences.size(), n_quadrature_points);
755 
756  std::fill (divergences.begin(), divergences.end(),
757  typename Vector<dim,spacedim>::template OutputType<Number>::divergence_type());
758 
759  for (unsigned int shape_function=0;
760  shape_function<dofs_per_cell; ++shape_function)
761  {
762  const int snc = shape_function_data[shape_function].single_nonzero_component;
763 
764  if (snc == -2)
765  // shape function is zero for the selected components
766  continue;
767 
768  const Number &value = dof_values[shape_function];
769  // For auto-differentiable numbers, the fact that a DoF value is zero
770  // does not imply that its derivatives are zero as well. So we
771  // can't filter by value for these number types.
773  if (value == ::internal::NumberType<Number>::value(0.0))
774  continue;
775 
776  if (snc != -1)
777  {
778  const unsigned int comp =
779  shape_function_data[shape_function].single_nonzero_component_index;
780  const ::Tensor<1,spacedim> *shape_gradient_ptr = &shape_gradients[snc][0];
781  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
782  divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
783  }
784  else
785  for (unsigned int d=0; d<spacedim; ++d)
786  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
787  {
788  const ::Tensor<1,spacedim> *shape_gradient_ptr =
789  &shape_gradients[shape_function_data[shape_function].
790  row_index[d]][0];
791  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
792  divergences[q_point] += value * (*shape_gradient_ptr++)[d];
793  }
794  }
795  }
796 
797 
798 
799  template <int dim, int spacedim, typename Number>
800  void
801  do_function_curls (const ArrayView<Number> &dof_values,
802  const Table<2,::Tensor<1,spacedim> > &shape_gradients,
803  const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
804  std::vector<typename ProductType<Number,typename ::internal::CurlType<spacedim>::type>::type> &curls)
805  {
806  const unsigned int dofs_per_cell = dof_values.size();
807  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
808  shape_gradients[0].size() : curls.size();
809  AssertDimension (curls.size(), n_quadrature_points);
810 
811  std::fill (curls.begin(), curls.end(),
813 
814  switch (spacedim)
815  {
816  case 1:
817  {
818  Assert (false, ExcMessage("Computing the curl in 1d is not a useful operation"));
819  break;
820  }
821 
822  case 2:
823  {
824  for (unsigned int shape_function = 0;
825  shape_function < dofs_per_cell; ++shape_function)
826  {
827  const int snc = shape_function_data[shape_function].single_nonzero_component;
828 
829  if (snc == -2)
830  // shape function is zero for the selected components
831  continue;
832 
833  const Number &value = dof_values[shape_function];
834  // For auto-differentiable numbers, the fact that a DoF value is zero
835  // does not imply that its derivatives are zero as well. So we
836  // can't filter by value for these number types.
838  if (value == ::internal::NumberType<Number>::value(0.0))
839  continue;
840 
841  if (snc != -1)
842  {
843  const ::Tensor<1, spacedim> *shape_gradient_ptr =
844  &shape_gradients[snc][0];
845 
846  Assert (shape_function_data[shape_function].single_nonzero_component >= 0,
847  ExcInternalError());
848  // we're in 2d, so the formula for the curl is simple:
849  if (shape_function_data[shape_function].single_nonzero_component_index == 0)
850  for (unsigned int q_point = 0;
851  q_point < n_quadrature_points; ++q_point)
852  curls[q_point][0] -= value * (*shape_gradient_ptr++)[1];
853  else
854  for (unsigned int q_point = 0;
855  q_point < n_quadrature_points; ++q_point)
856  curls[q_point][0] += value * (*shape_gradient_ptr++)[0];
857  }
858  else
859  // we have multiple non-zero components in the shape functions. not
860  // all of them must necessarily be within the 2-component window
861  // this FEValuesViews::Vector object considers, however.
862  {
863  if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
864  {
865  const ::Tensor<1,spacedim> *shape_gradient_ptr =
866  &shape_gradients[shape_function_data[shape_function].row_index[0]][0];
867 
868  for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
869  curls[q_point][0] -= value * (*shape_gradient_ptr++)[1];
870  }
871 
872  if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
873  {
874  const ::Tensor<1,spacedim> *shape_gradient_ptr =
875  &shape_gradients[shape_function_data[shape_function].row_index[1]][0];
876 
877  for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
878  curls[q_point][0] += value * (*shape_gradient_ptr++)[0];
879  }
880  }
881  }
882  break;
883  }
884 
885  case 3:
886  {
887  for (unsigned int shape_function = 0;
888  shape_function < dofs_per_cell; ++shape_function)
889  {
890  const int snc = shape_function_data[shape_function].single_nonzero_component;
891 
892  if (snc == -2)
893  // shape function is zero for the selected components
894  continue;
895 
896  const Number &value = dof_values[shape_function];
897  // For auto-differentiable numbers, the fact that a DoF value is zero
898  // does not imply that its derivatives are zero as well. So we
899  // can't filter by value for these number types.
901  if (value == ::internal::NumberType<Number>::value(0.0))
902  continue;
903 
904  if (snc != -1)
905  {
906  const ::Tensor<1, spacedim> *shape_gradient_ptr = &shape_gradients[snc][0];
907 
908  switch (shape_function_data[shape_function].single_nonzero_component_index)
909  {
910  case 0:
911  {
912  for (unsigned int q_point = 0;
913  q_point < n_quadrature_points; ++q_point)
914  {
915  curls[q_point][1] += value * (*shape_gradient_ptr)[2];
916  curls[q_point][2] -= value * (*shape_gradient_ptr++)[1];
917  }
918 
919  break;
920  }
921 
922  case 1:
923  {
924  for (unsigned int q_point = 0;
925  q_point < n_quadrature_points; ++q_point)
926  {
927  curls[q_point][0] -= value * (*shape_gradient_ptr)[2];
928  curls[q_point][2] += value * (*shape_gradient_ptr++)[0];
929  }
930 
931  break;
932  }
933 
934  case 2:
935  {
936  for (unsigned int q_point = 0;
937  q_point < n_quadrature_points; ++q_point)
938  {
939  curls[q_point][0] += value * (*shape_gradient_ptr)[1];
940  curls[q_point][1] -= value * (*shape_gradient_ptr++)[0];
941  }
942  break;
943  }
944 
945  default:
946  Assert (false, ExcInternalError());
947  }
948  }
949 
950  else
951  // we have multiple non-zero components in the shape functions. not
952  // all of them must necessarily be within the 3-component window
953  // this FEValuesViews::Vector object considers, however.
954  {
955  if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
956  {
957  const ::Tensor<1,spacedim> *shape_gradient_ptr =
958  &shape_gradients[shape_function_data[shape_function].row_index[0]][0];
959 
960  for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
961  {
962  curls[q_point][1] += value * (*shape_gradient_ptr)[2];
963  curls[q_point][2] -= value * (*shape_gradient_ptr++)[1];
964  }
965  }
966 
967  if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
968  {
969  const ::Tensor<1,spacedim> *shape_gradient_ptr =
970  &shape_gradients[shape_function_data[shape_function].row_index[1]][0];
971 
972  for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
973  {
974  curls[q_point][0] -= value * (*shape_gradient_ptr)[2];
975  curls[q_point][2] += value * (*shape_gradient_ptr++)[0];
976  }
977  }
978 
979  if (shape_function_data[shape_function].is_nonzero_shape_function_component[2])
980  {
981  const ::Tensor<1,spacedim> *shape_gradient_ptr =
982  &shape_gradients[shape_function_data[shape_function].row_index[2]][0];
983 
984  for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
985  {
986  curls[q_point][0] += value * (*shape_gradient_ptr)[1];
987  curls[q_point][1] -= value * (*shape_gradient_ptr++)[0];
988  }
989  }
990  }
991  }
992  }
993  }
994  }
995 
996 
997 
998  template <int dim, int spacedim, typename Number>
999  void
1000  do_function_laplacians (const ArrayView<Number> &dof_values,
1001  const Table<2,::Tensor<2,spacedim> > &shape_hessians,
1002  const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
1003  std::vector<typename Vector<dim,spacedim>::template OutputType<Number>::laplacian_type> &laplacians)
1004  {
1005  const unsigned int dofs_per_cell = dof_values.size();
1006  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
1007  shape_hessians[0].size() : laplacians.size();
1008  AssertDimension (laplacians.size(), n_quadrature_points);
1009 
1010  std::fill (laplacians.begin(), laplacians.end(),
1011  typename Vector<dim,spacedim>::template OutputType<Number>::laplacian_type());
1012 
1013  for (unsigned int shape_function=0;
1014  shape_function<dofs_per_cell; ++shape_function)
1015  {
1016  const int snc = shape_function_data[shape_function].single_nonzero_component;
1017 
1018  if (snc == -2)
1019  // shape function is zero for the selected components
1020  continue;
1021 
1022  const Number &value = dof_values[shape_function];
1023  // For auto-differentiable numbers, the fact that a DoF value is zero
1024  // does not imply that its derivatives are zero as well. So we
1025  // can't filter by value for these number types.
1027  if (value == ::internal::NumberType<Number>::value(0.0))
1028  continue;
1029 
1030  if (snc != -1)
1031  {
1032  const unsigned int comp =
1033  shape_function_data[shape_function].single_nonzero_component_index;
1034  const ::Tensor<2,spacedim> *shape_hessian_ptr =
1035  &shape_hessians[snc][0];
1036  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
1037  laplacians[q_point][comp] += value * trace(*shape_hessian_ptr++);
1038  }
1039  else
1040  for (unsigned int d=0; d<spacedim; ++d)
1041  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
1042  {
1043  const ::Tensor<2,spacedim> *shape_hessian_ptr =
1044  &shape_hessians[shape_function_data[shape_function].
1045  row_index[d]][0];
1046  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
1047  laplacians[q_point][d] += value * trace(*shape_hessian_ptr++);
1048  }
1049  }
1050  }
1051 
1052 
1053 
1054  // ---------------------- symmetric tensor part ------------------------
1055 
1056  template <int dim, int spacedim, typename Number>
1057  void
1058  do_function_values (const ArrayView<Number> &dof_values,
1059  const ::Table<2,double> &shape_values,
1060  const std::vector<typename SymmetricTensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
1061  std::vector<typename ProductType<Number,::SymmetricTensor<2,spacedim> >::type> &values)
1062  {
1063  const unsigned int dofs_per_cell = dof_values.size();
1064  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
1065  shape_values.n_cols() : values.size();
1066  AssertDimension (values.size(), n_quadrature_points);
1067 
1068  std::fill (values.begin(), values.end(),
1070 
1071  for (unsigned int shape_function=0;
1072  shape_function<dofs_per_cell; ++shape_function)
1073  {
1074  const int snc = shape_function_data[shape_function].single_nonzero_component;
1075 
1076  if (snc == -2)
1077  // shape function is zero for the selected components
1078  continue;
1079 
1080  const Number &value = dof_values[shape_function];
1081  // For auto-differentiable numbers, the fact that a DoF value is zero
1082  // does not imply that its derivatives are zero as well. So we
1083  // can't filter by value for these number types.
1085  if (value == ::internal::NumberType<Number>::value(0.0))
1086  continue;
1087 
1088  if (snc != -1)
1089  {
1090  const TableIndices<2> comp =
1092  (shape_function_data[shape_function].single_nonzero_component_index);
1093  const double *shape_value_ptr = &shape_values(snc,0);
1094  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
1095  values[q_point][comp] += value * (*shape_value_ptr++);
1096  }
1097  else
1098  for (unsigned int d=0;
1099  d<::SymmetricTensor<2,spacedim>::n_independent_components; ++d)
1100  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
1101  {
1102  const TableIndices<2> comp =
1104  const double *shape_value_ptr =
1105  &shape_values(shape_function_data[shape_function].row_index[d],0);
1106  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
1107  values[q_point][comp] += value * (*shape_value_ptr++);
1108  }
1109  }
1110  }
1111 
1112 
1113 
1114  template <int dim, int spacedim, typename Number>
1115  void
1116  do_function_divergences (const ArrayView<Number> &dof_values,
1117  const Table<2,::Tensor<1,spacedim> > &shape_gradients,
1118  const std::vector<typename SymmetricTensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
1119  std::vector<typename SymmetricTensor<2,dim,spacedim>::template OutputType<Number>::divergence_type> &divergences)
1120  {
1121  const unsigned int dofs_per_cell = dof_values.size();
1122  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
1123  shape_gradients[0].size() : divergences.size();
1124  AssertDimension (divergences.size(), n_quadrature_points);
1125 
1126  std::fill (divergences.begin(), divergences.end(),
1127  typename SymmetricTensor<2,dim,spacedim>::template OutputType<Number>::divergence_type());
1128 
1129  for (unsigned int shape_function=0;
1130  shape_function<dofs_per_cell; ++shape_function)
1131  {
1132  const int snc = shape_function_data[shape_function].single_nonzero_component;
1133 
1134  if (snc == -2)
1135  // shape function is zero for the selected components
1136  continue;
1137 
1138  const Number &value = dof_values[shape_function];
1139  // For auto-differentiable numbers, the fact that a DoF value is zero
1140  // does not imply that its derivatives are zero as well. So we
1141  // can't filter by value for these number types.
1143  if (value == ::internal::NumberType<Number>::value(0.0))
1144  continue;
1145 
1146  if (snc != -1)
1147  {
1148  const unsigned int comp =
1149  shape_function_data[shape_function].single_nonzero_component_index;
1150 
1151  const ::Tensor < 1, spacedim> *shape_gradient_ptr =
1152  &shape_gradients[snc][0];
1153 
1154  const unsigned int ii = ::SymmetricTensor<2,spacedim>::
1156  const unsigned int jj = ::SymmetricTensor<2,spacedim>::
1158 
1159  for (unsigned int q_point = 0; q_point < n_quadrature_points;
1160  ++q_point, ++shape_gradient_ptr)
1161  {
1162  divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
1163 
1164  if (ii != jj)
1165  divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
1166  }
1167  }
1168  else
1169  {
1170  for (unsigned int d = 0;
1171  d < ::SymmetricTensor<2,spacedim>::n_independent_components; ++d)
1172  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
1173  {
1174  Assert (false, ExcNotImplemented());
1175 
1176  // the following implementation needs to be looked over -- I
1177  // think it can't be right, because we are in a case where
1178  // there is no single nonzero component
1179  //
1180  // the following is not implemented! we need to consider the
1181  // interplay between multiple non-zero entries in shape
1182  // function and the representation as a symmetric
1183  // second-order tensor
1184  const unsigned int comp =
1185  shape_function_data[shape_function].single_nonzero_component_index;
1186 
1187  const ::Tensor < 1, spacedim> *shape_gradient_ptr =
1188  &shape_gradients[shape_function_data[shape_function].
1189  row_index[d]][0];
1190  for (unsigned int q_point = 0; q_point < n_quadrature_points;
1191  ++q_point, ++shape_gradient_ptr)
1192  {
1193  for (unsigned int j = 0; j < spacedim; ++j)
1194  {
1195  const unsigned int vector_component = ::SymmetricTensor<2,spacedim>::component_to_unrolled_index (TableIndices<2>(comp,j));
1196  divergences[q_point][vector_component] += value * (*shape_gradient_ptr++)[j];
1197  }
1198  }
1199  }
1200  }
1201  }
1202  }
1203 
1204  // ---------------------- non-symmetric tensor part ------------------------
1205 
1206  template <int dim, int spacedim, typename Number>
1207  void
1208  do_function_values (const ArrayView<Number> &dof_values,
1209  const ::Table<2,double> &shape_values,
1210  const std::vector<typename Tensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
1211  std::vector<typename ProductType<Number,::Tensor<2,spacedim> >::type> &values)
1212  {
1213  const unsigned int dofs_per_cell = dof_values.size();
1214  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
1215  shape_values.n_cols() : values.size();
1216  AssertDimension (values.size(), n_quadrature_points);
1217 
1218  std::fill (values.begin(), values.end(),
1219  typename ProductType<Number,::Tensor<2,spacedim> >::type());
1220 
1221  for (unsigned int shape_function=0;
1222  shape_function<dofs_per_cell; ++shape_function)
1223  {
1224  const int snc = shape_function_data[shape_function].single_nonzero_component;
1225 
1226  if (snc == -2)
1227  // shape function is zero for the selected components
1228  continue;
1229 
1230  const Number &value = dof_values[shape_function];
1231  // For auto-differentiable numbers, the fact that a DoF value is zero
1232  // does not imply that its derivatives are zero as well. So we
1233  // can't filter by value for these number types.
1235  if (value == ::internal::NumberType<Number>::value(0.0))
1236  continue;
1237 
1238  if (snc != -1)
1239  {
1240  const unsigned int comp =
1241  shape_function_data[shape_function].single_nonzero_component_index;
1242 
1244 
1245  const double *shape_value_ptr = &shape_values(snc,0);
1246  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
1247  values[q_point][indices] += value * (*shape_value_ptr++);
1248  }
1249  else
1250  for (unsigned int d=0;
1251  d<dim*dim; ++d)
1252  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
1253  {
1255 
1256  const double *shape_value_ptr =
1257  &shape_values(shape_function_data[shape_function].row_index[d],0);
1258  for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
1259  values[q_point][indices] += value * (*shape_value_ptr++);
1260  }
1261  }
1262  }
1263 
1264 
1265 
1266  template <int dim, int spacedim, typename Number>
1267  void
1268  do_function_divergences (const ArrayView<Number> &dof_values,
1269  const Table<2,::Tensor<1,spacedim> > &shape_gradients,
1270  const std::vector<typename Tensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
1271  std::vector<typename Tensor<2,dim,spacedim>::template OutputType<Number>::divergence_type> &divergences)
1272  {
1273  const unsigned int dofs_per_cell = dof_values.size();
1274  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
1275  shape_gradients[0].size() : divergences.size();
1276  AssertDimension (divergences.size(), n_quadrature_points);
1277 
1278  std::fill (divergences.begin(), divergences.end(),
1279  typename Tensor<2,dim,spacedim>::template OutputType<Number>::divergence_type());
1280 
1281  for (unsigned int shape_function=0;
1282  shape_function<dofs_per_cell; ++shape_function)
1283  {
1284  const int snc = shape_function_data[shape_function].single_nonzero_component;
1285 
1286  if (snc == -2)
1287  // shape function is zero for the selected components
1288  continue;
1289 
1290  const Number &value = dof_values[shape_function];
1291  // For auto-differentiable numbers, the fact that a DoF value is zero
1292  // does not imply that its derivatives are zero as well. So we
1293  // can't filter by value for these number types.
1295  if (value == ::internal::NumberType<Number>::value(0.0))
1296  continue;
1297 
1298  if (snc != -1)
1299  {
1300  const unsigned int comp =
1301  shape_function_data[shape_function].single_nonzero_component_index;
1302 
1303  const ::Tensor < 1, spacedim> *shape_gradient_ptr =
1304  &shape_gradients[snc][0];
1305 
1307  const unsigned int ii = indices[0];
1308  const unsigned int jj = indices[1];
1309 
1310  for (unsigned int q_point = 0; q_point < n_quadrature_points;
1311  ++q_point, ++shape_gradient_ptr)
1312  {
1313  divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
1314  }
1315  }
1316  else
1317  {
1318  for (unsigned int d = 0;
1319  d < dim*dim; ++d)
1320  if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
1321  {
1322  Assert (false, ExcNotImplemented());
1323  }
1324  }
1325  }
1326  }
1327 
1328  } // end of namespace internal
1329 
1330 
1331 
1332  template <int dim, int spacedim>
1333  template <class InputVector>
1334  void
1336  get_function_values (const InputVector &fe_function,
1337  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const
1338  {
1339  Assert (fe_values->update_flags & update_values,
1340  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
1341  Assert (fe_values->present_cell.get() != nullptr,
1342  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1343  AssertDimension (fe_function.size(),
1344  fe_values->present_cell->n_dofs_for_dof_handler());
1345 
1346  // get function values of dofs on this cell and call internal worker function
1347  ::Vector<typename InputVector::value_type> dof_values(fe_values->dofs_per_cell);
1348  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1349  internal::do_function_values<dim,spacedim>
1350  (make_array_view(dof_values.begin(), dof_values.end()),
1351  fe_values->finite_element_output.shape_values, shape_function_data, values);
1352  }
1353 
1354  template <int dim, int spacedim>
1355  template <class InputVector>
1356  void
1358  get_function_values_from_local_dof_values (const InputVector &dof_values,
1359  std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const
1360  {
1361  Assert (fe_values->update_flags & update_values,
1362  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
1363  Assert (fe_values->present_cell.get() != nullptr,
1364  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1365  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1366 
1367  internal::do_function_values<dim,spacedim>
1368  (make_array_view(dof_values.begin(), dof_values.end()),
1369  fe_values->finite_element_output.shape_values, shape_function_data, values);
1370  }
1371 
1372 
1373 
1374  template <int dim, int spacedim>
1375  template <class InputVector>
1376  void
1378  get_function_gradients (const InputVector &fe_function,
1379  std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const
1380  {
1381  Assert (fe_values->update_flags & update_gradients,
1382  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1383  Assert (fe_values->present_cell.get() != nullptr,
1384  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1385  AssertDimension (fe_function.size(),
1386  fe_values->present_cell->n_dofs_for_dof_handler());
1387 
1388  // get function values of dofs on this cell
1389  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1390  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1391  internal::do_function_derivatives<1,dim,spacedim>
1392  (make_array_view(dof_values.begin(), dof_values.end()),
1393  fe_values->finite_element_output.shape_gradients, shape_function_data, gradients);
1394  }
1395 
1396 
1397 
1398  template <int dim, int spacedim>
1399  template <class InputVector>
1400  void
1403  std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const
1404  {
1405  Assert (fe_values->update_flags & update_gradients,
1406  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1407  Assert (fe_values->present_cell.get() != nullptr,
1408  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1409  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1410 
1411  internal::do_function_derivatives<1,dim,spacedim>
1412  (make_array_view(dof_values.begin(), dof_values.end()),
1413  fe_values->finite_element_output.shape_gradients, shape_function_data, gradients);
1414  }
1415 
1416 
1417 
1418  template <int dim, int spacedim>
1419  template <class InputVector>
1420  void
1422  get_function_hessians (const InputVector &fe_function,
1423  std::vector<typename ProductType<hessian_type,typename InputVector::value_type>::type> &hessians) const
1424  {
1425  Assert (fe_values->update_flags & update_hessians,
1426  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1427  Assert (fe_values->present_cell.get() != nullptr,
1428  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1429  AssertDimension (fe_function.size(),
1430  fe_values->present_cell->n_dofs_for_dof_handler());
1431 
1432  // get function values of dofs on this cell
1433  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1434  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1435  internal::do_function_derivatives<2,dim,spacedim>
1436  (make_array_view(dof_values.begin(), dof_values.end()),
1437  fe_values->finite_element_output.shape_hessians, shape_function_data, hessians);
1438  }
1439 
1440 
1441 
1442  template <int dim, int spacedim>
1443  template <class InputVector>
1444  void
1446  get_function_hessians_from_local_dof_values(const InputVector &dof_values,
1447  std::vector<typename OutputType<typename InputVector::value_type>::hessian_type> &hessians) const
1448  {
1449  Assert (fe_values->update_flags & update_hessians,
1450  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1451  Assert (fe_values->present_cell.get() != nullptr,
1452  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1453  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1454 
1455  internal::do_function_derivatives<2,dim,spacedim>
1456  (make_array_view(dof_values.begin(), dof_values.end()),
1457  fe_values->finite_element_output.shape_hessians, shape_function_data, hessians);
1458  }
1459 
1460 
1461 
1462  template <int dim, int spacedim>
1463  template <class InputVector>
1464  void
1466  get_function_laplacians (const InputVector &fe_function,
1467  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &laplacians) const
1468  {
1469  Assert (fe_values->update_flags & update_hessians,
1470  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1471  Assert (fe_values->present_cell.get() != nullptr,
1472  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1473  AssertDimension (fe_function.size(),
1474  fe_values->present_cell->n_dofs_for_dof_handler());
1475 
1476  // get function values of dofs on this cell
1477  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1478  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1479  internal::do_function_laplacians<dim,spacedim>
1480  (make_array_view(dof_values.begin(), dof_values.end()),
1481  fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians);
1482  }
1483 
1484 
1485 
1486  template <int dim, int spacedim>
1487  template <class InputVector>
1488  void
1491  std::vector<typename OutputType<typename InputVector::value_type>::laplacian_type> &laplacians) const
1492  {
1493  Assert (fe_values->update_flags & update_hessians,
1494  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1495  Assert (fe_values->present_cell.get() != nullptr,
1496  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1497  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1498 
1499  internal::do_function_laplacians<dim,spacedim>
1500  (make_array_view(dof_values.begin(), dof_values.end()),
1501  fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians);
1502  }
1503 
1504 
1505 
1506  template <int dim, int spacedim>
1507  template <class InputVector>
1508  void
1510  get_function_third_derivatives (const InputVector &fe_function,
1511  std::vector<typename ProductType<third_derivative_type,typename InputVector::value_type>::type> &third_derivatives) const
1512  {
1513  Assert (fe_values->update_flags & update_3rd_derivatives,
1514  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_3rd_derivatives")));
1515  Assert (fe_values->present_cell.get() != nullptr,
1516  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1517  AssertDimension (fe_function.size(),
1518  fe_values->present_cell->n_dofs_for_dof_handler());
1519 
1520  // get function values of dofs on this cell
1521  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1522  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1523  internal::do_function_derivatives<3,dim,spacedim>
1524  (make_array_view(dof_values.begin(), dof_values.end()),
1525  fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives);
1526  }
1527 
1528 
1529 
1530  template <int dim, int spacedim>
1531  template <class InputVector>
1532  void
1535  std::vector<typename OutputType<typename InputVector::value_type>::third_derivative_type> &third_derivatives) const
1536  {
1537  Assert (fe_values->update_flags & update_3rd_derivatives,
1538  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_3rd_derivatives")));
1539  Assert (fe_values->present_cell.get() != nullptr,
1540  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1541  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1542 
1543  internal::do_function_derivatives<3,dim,spacedim>
1544  (make_array_view(dof_values.begin(), dof_values.end()),
1545  fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives);
1546  }
1547 
1548 
1549 
1550  template <int dim, int spacedim>
1551  template <class InputVector>
1552  void
1554  get_function_values (const InputVector &fe_function,
1555  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const
1556  {
1557  Assert (fe_values->update_flags & update_values,
1558  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
1559  Assert (fe_values->present_cell.get() != nullptr,
1560  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1561  AssertDimension (fe_function.size(),
1562  fe_values->present_cell->n_dofs_for_dof_handler());
1563 
1564  // get function values of dofs on this cell
1565  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1566  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1567  internal::do_function_values<dim,spacedim>
1568  (make_array_view(dof_values.begin(), dof_values.end()),
1569  fe_values->finite_element_output.shape_values, shape_function_data, values);
1570  }
1571 
1572 
1573 
1574  template <int dim, int spacedim>
1575  template <class InputVector>
1576  void
1578  get_function_values_from_local_dof_values (const InputVector &dof_values,
1579  std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const
1580  {
1581  Assert (fe_values->update_flags & update_values,
1582  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
1583  Assert (fe_values->present_cell.get() != nullptr,
1584  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1585  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1586 
1587  internal::do_function_values<dim,spacedim>
1588  (make_array_view(dof_values.begin(), dof_values.end()),
1589  fe_values->finite_element_output.shape_values, shape_function_data, values);
1590  }
1591 
1592 
1593 
1594  template <int dim, int spacedim>
1595  template <class InputVector>
1596  void
1598  get_function_gradients (const InputVector &fe_function,
1599  std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const
1600  {
1601  Assert (fe_values->update_flags & update_gradients,
1602  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1603  Assert (fe_values->present_cell.get() != nullptr,
1604  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1605  AssertDimension (fe_function.size(),
1606  fe_values->present_cell->n_dofs_for_dof_handler());
1607 
1608  // get function values of dofs on this cell
1609  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1610  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1611  internal::do_function_derivatives<1,dim,spacedim>
1612  (make_array_view(dof_values.begin(), dof_values.end()),
1613  fe_values->finite_element_output.shape_gradients, shape_function_data, gradients);
1614  }
1615 
1616 
1617 
1618  template <int dim, int spacedim>
1619  template <class InputVector>
1620  void
1622  get_function_gradients_from_local_dof_values (const InputVector &dof_values,
1623  std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const
1624  {
1625  Assert (fe_values->update_flags & update_gradients,
1626  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1627  Assert (fe_values->present_cell.get() != nullptr,
1628  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1629  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1630 
1631  internal::do_function_derivatives<1,dim,spacedim>
1632  (make_array_view(dof_values.begin(), dof_values.end()),
1633  fe_values->finite_element_output.shape_gradients, shape_function_data, gradients);
1634  }
1635 
1636 
1637 
1638  template <int dim, int spacedim>
1639  template <class InputVector>
1640  void
1642  get_function_symmetric_gradients (const InputVector &fe_function,
1643  std::vector<typename ProductType<symmetric_gradient_type,typename InputVector::value_type>::type> &symmetric_gradients) const
1644  {
1645  Assert (fe_values->update_flags & update_gradients,
1646  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1647  Assert (fe_values->present_cell.get() != nullptr,
1648  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1649  AssertDimension (fe_function.size(),
1650  fe_values->present_cell->n_dofs_for_dof_handler());
1651 
1652  // get function values of dofs on this cell
1653  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1654  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1655  internal::do_function_symmetric_gradients<dim,spacedim>
1656  (make_array_view(dof_values.begin(), dof_values.end()),
1657  fe_values->finite_element_output.shape_gradients, shape_function_data,
1658  symmetric_gradients);
1659  }
1660 
1661 
1662 
1663  template <int dim, int spacedim>
1664  template <class InputVector>
1665  void
1668  std::vector<typename OutputType<typename InputVector::value_type>::symmetric_gradient_type> &symmetric_gradients) const
1669  {
1670  Assert (fe_values->update_flags & update_gradients,
1671  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1672  Assert (fe_values->present_cell.get() != nullptr,
1673  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1674  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1675 
1676  internal::do_function_symmetric_gradients<dim,spacedim>
1677  (make_array_view(dof_values.begin(), dof_values.end()),
1678  fe_values->finite_element_output.shape_gradients, shape_function_data,
1679  symmetric_gradients);
1680  }
1681 
1682 
1683 
1684  template <int dim, int spacedim>
1685  template <class InputVector>
1686  void
1688  get_function_divergences (const InputVector &fe_function,
1689  std::vector<typename ProductType<divergence_type,typename InputVector::value_type>::type> &divergences) const
1690  {
1691  Assert (fe_values->update_flags & update_gradients,
1692  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1693  Assert (fe_values->present_cell.get() != nullptr,
1694  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1695  AssertDimension (fe_function.size(),
1696  fe_values->present_cell->n_dofs_for_dof_handler());
1697 
1698  // get function values of dofs
1699  // on this cell
1700  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1701  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1702  internal::do_function_divergences<dim,spacedim>
1703  (make_array_view(dof_values.begin(), dof_values.end()),
1704  fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
1705  }
1706 
1707 
1708 
1709  template <int dim, int spacedim>
1710  template <class InputVector>
1711  void
1714  std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const
1715  {
1716  Assert (fe_values->update_flags & update_gradients,
1717  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1718  Assert (fe_values->present_cell.get() != nullptr,
1719  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1720  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1721 
1722  internal::do_function_divergences<dim,spacedim>
1723  (make_array_view(dof_values.begin(), dof_values.end()),
1724  fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
1725  }
1726 
1727 
1728 
1729  template <int dim, int spacedim>
1730  template <class InputVector>
1731  void
1733  get_function_curls (const InputVector &fe_function,
1734  std::vector<typename ProductType<curl_type,typename InputVector::value_type>::type> &curls) const
1735  {
1736  Assert (fe_values->update_flags & update_gradients,
1737  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1738  Assert (fe_values->present_cell.get () != nullptr,
1739  ExcMessage ("FEValues object is not reinited to any cell"));
1740  AssertDimension (fe_function.size (),
1741  fe_values->present_cell->n_dofs_for_dof_handler ());
1742 
1743  // get function values of dofs on this cell
1744  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1745  fe_values->present_cell->get_interpolated_dof_values (fe_function, dof_values);
1746  internal::do_function_curls<dim,spacedim>
1747  (make_array_view(dof_values.begin(), dof_values.end()),
1748  fe_values->finite_element_output.shape_gradients, shape_function_data, curls);
1749  }
1750 
1751 
1752 
1753  template <int dim, int spacedim>
1754  template <class InputVector>
1755  void
1757  get_function_curls_from_local_dof_values(const InputVector &dof_values,
1758  std::vector<typename OutputType<typename InputVector::value_type>::curl_type> &curls) const
1759  {
1760  Assert (fe_values->update_flags & update_gradients,
1761  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1762  Assert (fe_values->present_cell.get () != nullptr,
1763  ExcMessage ("FEValues object is not reinited to any cell"));
1764  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1765 
1766  internal::do_function_curls<dim,spacedim>
1767  (make_array_view(dof_values.begin(), dof_values.end()),
1768  fe_values->finite_element_output.shape_gradients, shape_function_data, curls);
1769  }
1770 
1771 
1772 
1773  template <int dim, int spacedim>
1774  template <class InputVector>
1775  void
1777  get_function_hessians (const InputVector &fe_function,
1778  std::vector<typename ProductType<hessian_type,typename InputVector::value_type>::type> &hessians) const
1779  {
1780  Assert (fe_values->update_flags & update_hessians,
1781  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1782  Assert (fe_values->present_cell.get() != nullptr,
1783  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1784  AssertDimension (fe_function.size(),
1785  fe_values->present_cell->n_dofs_for_dof_handler());
1786 
1787  // get function values of dofs on this cell
1788  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1789  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1790  internal::do_function_derivatives<2,dim,spacedim>
1791  (make_array_view(dof_values.begin(), dof_values.end()),
1792  fe_values->finite_element_output.shape_hessians, shape_function_data, hessians);
1793  }
1794 
1795 
1796 
1797  template <int dim, int spacedim>
1798  template <class InputVector>
1799  void
1801  get_function_hessians_from_local_dof_values (const InputVector &dof_values,
1802  std::vector<typename OutputType<typename InputVector::value_type>::hessian_type> &hessians) const
1803  {
1804  Assert (fe_values->update_flags & update_hessians,
1805  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1806  Assert (fe_values->present_cell.get() != nullptr,
1807  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1808  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1809 
1810  internal::do_function_derivatives<2,dim,spacedim>
1811  (make_array_view(dof_values.begin(), dof_values.end()),
1812  fe_values->finite_element_output.shape_hessians, shape_function_data, hessians);
1813  }
1814 
1815 
1816 
1817  template <int dim, int spacedim>
1818  template <class InputVector>
1819  void
1821  get_function_laplacians (const InputVector &fe_function,
1822  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &laplacians) const
1823  {
1824  Assert (fe_values->update_flags & update_hessians,
1825  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1826  Assert (laplacians.size() == fe_values->n_quadrature_points,
1827  ExcDimensionMismatch(laplacians.size(), fe_values->n_quadrature_points));
1828  Assert (fe_values->present_cell.get() != nullptr,
1829  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1830  Assert (fe_function.size() == fe_values->present_cell->n_dofs_for_dof_handler(),
1831  ExcDimensionMismatch(fe_function.size(),
1832  fe_values->present_cell->n_dofs_for_dof_handler()));
1833 
1834  // get function values of dofs on this cell
1835  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1836  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1837  internal::do_function_laplacians<dim,spacedim>
1838  (make_array_view(dof_values.begin(), dof_values.end()),
1839  fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians);
1840  }
1841 
1842 
1843 
1844 
1845  template <int dim, int spacedim>
1846  template <class InputVector>
1847  void
1850  std::vector<typename OutputType<typename InputVector::value_type>::laplacian_type> &laplacians) const
1851  {
1852  Assert (fe_values->update_flags & update_hessians,
1853  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_hessians")));
1854  Assert (laplacians.size() == fe_values->n_quadrature_points,
1855  ExcDimensionMismatch(laplacians.size(), fe_values->n_quadrature_points));
1856  Assert (fe_values->present_cell.get() != nullptr,
1857  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1858  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1859 
1860  internal::do_function_laplacians<dim,spacedim>
1861  (make_array_view(dof_values.begin(), dof_values.end()),
1862  fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians);
1863  }
1864 
1865 
1866  template <int dim, int spacedim>
1867  template <class InputVector>
1868  void
1870  get_function_third_derivatives (const InputVector &fe_function,
1871  std::vector<typename ProductType<third_derivative_type,typename InputVector::value_type>::type> &third_derivatives) const
1872  {
1873  Assert (fe_values->update_flags & update_3rd_derivatives,
1874  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_3rd_derivatives")));
1875  Assert (fe_values->present_cell.get() != nullptr,
1876  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1877  AssertDimension (fe_function.size(),
1878  fe_values->present_cell->n_dofs_for_dof_handler());
1879 
1880  // get function values of dofs on this cell
1881  ::Vector<typename InputVector::value_type> dof_values (fe_values->dofs_per_cell);
1882  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1883  internal::do_function_derivatives<3,dim,spacedim>
1884  (make_array_view(dof_values.begin(), dof_values.end()),
1885  fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives);
1886  }
1887 
1888 
1889 
1890  template <int dim, int spacedim>
1891  template <class InputVector>
1892  void
1895  std::vector<typename OutputType<typename InputVector::value_type>::third_derivative_type> &third_derivatives) const
1896  {
1897  Assert (fe_values->update_flags & update_3rd_derivatives,
1898  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_3rd_derivatives")));
1899  Assert (fe_values->present_cell.get() != nullptr,
1900  ExcMessage ("FEValues object is not reinit'ed to any cell"));
1901  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1902 
1903  internal::do_function_derivatives<3,dim,spacedim>
1904  (make_array_view(dof_values.begin(), dof_values.end()),
1905  fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives);
1906  }
1907 
1908 
1909 
1910  template <int dim, int spacedim>
1911  template <class InputVector>
1912  void
1914  get_function_values(const InputVector &fe_function,
1915  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const
1916  {
1917  Assert(fe_values->update_flags & update_values,
1918  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
1919  Assert(fe_values->present_cell.get() != nullptr,
1920  ExcMessage("FEValues object is not reinit'ed to any cell"));
1921  AssertDimension(fe_function.size(),
1922  fe_values->present_cell->n_dofs_for_dof_handler());
1923 
1924  // get function values of dofs on this cell
1925  ::Vector<typename InputVector::value_type> dof_values(fe_values->dofs_per_cell);
1926  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1927  internal::do_function_values<dim,spacedim>
1928  (make_array_view(dof_values.begin(), dof_values.end()),
1929  fe_values->finite_element_output.shape_values, shape_function_data, values);
1930  }
1931 
1932 
1933 
1934  template <int dim, int spacedim>
1935  template <class InputVector>
1936  void
1938  get_function_values_from_local_dof_values(const InputVector &dof_values,
1939  std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const
1940  {
1941  Assert(fe_values->update_flags & update_values,
1942  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
1943  Assert(fe_values->present_cell.get() != nullptr,
1944  ExcMessage("FEValues object is not reinit'ed to any cell"));
1945  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1946 
1947  internal::do_function_values<dim,spacedim>
1948  (make_array_view(dof_values.begin(), dof_values.end()),
1949  fe_values->finite_element_output.shape_values, shape_function_data, values);
1950  }
1951 
1952 
1953 
1954  template <int dim, int spacedim>
1955  template <class InputVector>
1956  void
1958  get_function_divergences(const InputVector &fe_function,
1959  std::vector<typename ProductType<divergence_type,typename InputVector::value_type>::type> &divergences) const
1960  {
1961  Assert(fe_values->update_flags & update_gradients,
1962  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1963  Assert(fe_values->present_cell.get() != nullptr,
1964  ExcMessage("FEValues object is not reinit'ed to any cell"));
1965  AssertDimension(fe_function.size(),
1966  fe_values->present_cell->n_dofs_for_dof_handler());
1967 
1968  // get function values of dofs
1969  // on this cell
1970  ::Vector<typename InputVector::value_type> dof_values(fe_values->dofs_per_cell);
1971  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
1972  internal::do_function_divergences<dim,spacedim>
1973  (make_array_view(dof_values.begin(), dof_values.end()),
1974  fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
1975  }
1976 
1977 
1978 
1979  template <int dim, int spacedim>
1980  template <class InputVector>
1981  void
1984  std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const
1985  {
1986  Assert(fe_values->update_flags & update_gradients,
1987  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
1988  Assert(fe_values->present_cell.get() != nullptr,
1989  ExcMessage("FEValues object is not reinit'ed to any cell"));
1990  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
1991 
1992  internal::do_function_divergences<dim,spacedim>
1993  (make_array_view(dof_values.begin(), dof_values.end()),
1994  fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
1995  }
1996 
1997 
1998 
1999  template <int dim, int spacedim>
2000  template <class InputVector>
2001  void
2003  get_function_values(const InputVector &fe_function,
2004  std::vector<typename ProductType<value_type,typename InputVector::value_type>::type> &values) const
2005  {
2006  Assert(fe_values->update_flags & update_values,
2007  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
2008  Assert(fe_values->present_cell.get() != nullptr,
2009  ExcMessage("FEValues object is not reinit'ed to any cell"));
2010  AssertDimension(fe_function.size(),
2011  fe_values->present_cell->n_dofs_for_dof_handler());
2012 
2013  // get function values of dofs on this cell
2014  ::Vector<typename InputVector::value_type> dof_values(fe_values->dofs_per_cell);
2015  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
2016  internal::do_function_values<dim,spacedim>
2017  (make_array_view(dof_values.begin(), dof_values.end()),
2018  fe_values->finite_element_output.shape_values, shape_function_data, values);
2019  }
2020 
2021 
2022 
2023  template <int dim, int spacedim>
2024  template <class InputVector>
2025  void
2027  get_function_values_from_local_dof_values (const InputVector &dof_values,
2028  std::vector<typename OutputType<typename InputVector::value_type>::value_type> &values) const
2029  {
2030  Assert(fe_values->update_flags & update_values,
2031  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_values")));
2032  Assert(fe_values->present_cell.get() != nullptr,
2033  ExcMessage("FEValues object is not reinit'ed to any cell"));
2034  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
2035 
2036  internal::do_function_values<dim,spacedim>
2037  (make_array_view(dof_values.begin(), dof_values.end()),
2038  fe_values->finite_element_output.shape_values, shape_function_data, values);
2039  }
2040 
2041 
2042 
2043  template <int dim, int spacedim>
2044  template <class InputVector>
2045  void
2047  get_function_divergences(const InputVector &fe_function,
2048  std::vector<typename ProductType<divergence_type,typename InputVector::value_type>::type> &divergences) const
2049  {
2050  Assert(fe_values->update_flags & update_gradients,
2051  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
2052  Assert(fe_values->present_cell.get() != nullptr,
2053  ExcMessage("FEValues object is not reinit'ed to any cell"));
2054  AssertDimension(fe_function.size(),
2055  fe_values->present_cell->n_dofs_for_dof_handler());
2056 
2057  // get function values of dofs
2058  // on this cell
2059  ::Vector<typename InputVector::value_type> dof_values(fe_values->dofs_per_cell);
2060  fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
2061  internal::do_function_divergences<dim,spacedim>
2062  (make_array_view(dof_values.begin(), dof_values.end()),
2063  fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
2064  }
2065 
2066 
2067 
2068  template <int dim, int spacedim>
2069  template <class InputVector>
2070  void
2073  std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const
2074  {
2075  Assert(fe_values->update_flags & update_gradients,
2076  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_gradients")));
2077  Assert(fe_values->present_cell.get() != nullptr,
2078  ExcMessage("FEValues object is not reinit'ed to any cell"));
2079  AssertDimension (dof_values.size(), fe_values->dofs_per_cell);
2080 
2081  internal::do_function_divergences<dim,spacedim>
2082  (make_array_view(dof_values.begin(), dof_values.end()),
2083  fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
2084  }
2085 }
2086 
2087 
2088 namespace internal
2089 {
2090  namespace FEValuesViews
2091  {
2092  template <int dim, int spacedim>
2094  {
2095  const FiniteElement<dim,spacedim> &fe = fe_values.get_fe();
2096 
2097  // create the views objects: Allocate a bunch of default-constructed ones
2098  // then destroy them again and do in-place construction of those we
2099  // actually want to use.
2100  const unsigned int n_scalars = fe.n_components();
2101  scalars.resize (n_scalars);
2102  for (unsigned int component=0; component<n_scalars; ++component)
2103  {
2104  // Use a typedef here to work around an issue with gcc-4.1:
2105  typedef ::FEValuesViews::Scalar<dim,spacedim> ScalarView;
2106  scalars[component].ScalarView::~ScalarView ();
2107 
2108  new (&scalars[component])
2110  component);
2111  }
2112 
2113  // compute number of vectors
2114  // that we can fit into
2115  // this finite element. note
2116  // that this is based on the
2117  // dimensionality 'dim' of the
2118  // manifold, not 'spacedim' of
2119  // the output vector
2120  const unsigned int n_vectors = (fe.n_components() >= spacedim ?
2121  fe.n_components()-spacedim+1 :
2122  0);
2123  vectors.resize (n_vectors);
2124  for (unsigned int component=0; component<n_vectors; ++component)
2125  {
2126  // Use a typedef here to work around an issue with gcc-4.1:
2127  typedef ::FEValuesViews::Vector<dim,spacedim> VectorView;
2128  vectors[component].VectorView::~VectorView ();
2129 
2130  new (&vectors[component])
2132  component);
2133  }
2134 
2135  // compute number of symmetric
2136  // tensors in the same way as above
2137  const unsigned int n_symmetric_second_order_tensors
2138  = (fe.n_components() >= (dim*dim + dim)/2 ?
2139  fe.n_components() - (dim*dim + dim)/2 + 1 :
2140  0);
2141  symmetric_second_order_tensors.resize(n_symmetric_second_order_tensors);
2142  for (unsigned int component = 0; component < n_symmetric_second_order_tensors; ++component)
2143  {
2144  // Use a typedef here to work around an issue with gcc-4.1:
2145  typedef ::FEValuesViews::SymmetricTensor<2, dim, spacedim> SymmetricTensorView;
2146  symmetric_second_order_tensors[component].SymmetricTensorView::~SymmetricTensorView();
2147 
2148  new (&symmetric_second_order_tensors[component])
2150  component);
2151  }
2152 
2153 
2154  // compute number of symmetric
2155  // tensors in the same way as above
2156  const unsigned int n_second_order_tensors
2157  = (fe.n_components() >= dim*dim ?
2158  fe.n_components() - dim*dim + 1 :
2159  0);
2160  second_order_tensors.resize(n_second_order_tensors);
2161  for (unsigned int component = 0; component < n_second_order_tensors; ++component)
2162  {
2163  // Use a typedef here to work around an issue with gcc-4.1:
2164  typedef ::FEValuesViews::Tensor<2, dim, spacedim> TensorView;
2165  second_order_tensors[component].TensorView::~TensorView();
2166 
2167  new (&second_order_tensors[component])
2169  component);
2170  }
2171  }
2172  }
2173 }
2174 
2175 
2176 /* ---------------- FEValuesBase<dim,spacedim>::CellIteratorBase --------- */
2177 
2178 template <int dim, int spacedim>
2179 class FEValuesBase<dim,spacedim>::CellIteratorBase
2180 {
2181 public:
2188  virtual ~CellIteratorBase () = default;
2189 
2203  virtual
2204  operator typename Triangulation<dim,spacedim>::cell_iterator () const = 0;
2205 
2213  virtual
2215  n_dofs_for_dof_handler () const = 0;
2216 
2217 #include "fe_values.decl.1.inst"
2218 
2223  virtual
2224  void
2225  get_interpolated_dof_values (const IndexSet &in,
2226  Vector<IndexSet::value_type> &out) const = 0;
2227 };
2228 
2229 /* ---------------- classes derived from FEValuesBase<dim,spacedim>::CellIteratorBase --------- */
2230 
2231 
2242 template <int dim, int spacedim>
2243 template <typename CI>
2244 class FEValuesBase<dim,spacedim>::CellIterator : public FEValuesBase<dim,spacedim>::CellIteratorBase
2245 {
2246 public:
2252  CellIterator (const CI &cell);
2253 
2267  virtual
2268  operator typename Triangulation<dim,spacedim>::cell_iterator () const;
2269 
2277  virtual
2279  n_dofs_for_dof_handler () const;
2280 
2281 #include "fe_values.decl.2.inst"
2282 
2287  virtual
2288  void
2289  get_interpolated_dof_values (const IndexSet &in,
2290  Vector<IndexSet::value_type> &out) const;
2291 
2292 private:
2297  const CI cell;
2298 };
2299 
2300 
2344 template <int dim, int spacedim>
2345 class FEValuesBase<dim,spacedim>::TriaCellIterator : public FEValuesBase<dim,spacedim>::CellIteratorBase
2346 {
2347 public:
2354 
2370  virtual
2371  operator typename Triangulation<dim,spacedim>::cell_iterator () const;
2372 
2380  virtual
2382  n_dofs_for_dof_handler () const;
2383 
2384 #include "fe_values.decl.2.inst"
2385 
2390  virtual
2391  void
2392  get_interpolated_dof_values (const IndexSet &in,
2393  Vector<IndexSet::value_type> &out) const;
2394 
2395 private:
2401 
2411  static const char *const message_string;
2412 };
2413 
2414 
2415 
2416 
2417 /* ---------------- FEValuesBase<dim,spacedim>::CellIterator<CI> --------- */
2418 
2419 
2420 template <int dim, int spacedim>
2421 template <typename CI>
2423  :
2424  cell(cell)
2425 {}
2426 
2427 
2428 
2429 template <int dim, int spacedim>
2430 template <typename CI>
2433 {
2434  return cell;
2435 }
2436 
2437 
2438 
2439 template <int dim, int spacedim>
2440 template <typename CI>
2443 {
2444  return cell->get_dof_handler().n_dofs();
2445 }
2446 
2447 
2448 
2449 #include "fe_values.impl.1.inst"
2450 
2451 
2452 template <int dim, int spacedim>
2453 template <typename CI>
2454 void
2457  Vector<IndexSet::value_type> &out) const
2458 {
2459  Assert (cell->has_children() == false, ExcNotImplemented());
2460 
2461  std::vector<types::global_dof_index> dof_indices (cell->get_fe().dofs_per_cell);
2462  cell->get_dof_indices (dof_indices);
2463 
2464  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
2465  out[i] = (in.is_element (dof_indices[i]) ? 1 : 0);
2466 }
2467 
2468 
2469 /* ---------------- FEValuesBase<dim,spacedim>::TriaCellIterator --------- */
2470 
2471 template <int dim, int spacedim>
2472 const char *const
2474  = ("You have previously called the FEValues::reinit function with a\n"
2475  "cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,\n"
2476  "when you do this, you cannot call some functions in the FEValues\n"
2477  "class, such as the get_function_values/gradients/hessians/third_derivatives\n"
2478  "functions. If you need these functions, then you need to call\n"
2479  "FEValues::reinit with an iterator type that allows to extract\n"
2480  "degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.");
2481 
2482 
2483 template <int dim, int spacedim>
2486  :
2487  cell(cell)
2488 {}
2489 
2490 
2491 
2492 template <int dim, int spacedim>
2495 {
2496  return cell;
2497 }
2498 
2499 
2500 
2501 template <int dim, int spacedim>
2504 {
2505  Assert (false, ExcMessage (message_string));
2506  return 0;
2507 }
2508 
2509 
2510 #include "fe_values.impl.2.inst"
2511 
2512 
2513 template <int dim, int spacedim>
2514 void
2518 {
2519  Assert (false, ExcMessage (message_string));
2520 }
2521 
2522 
2523 
2524 namespace internal
2525 {
2526  namespace FEValues
2527  {
2528  template <int dim, int spacedim>
2529  void
2530  MappingRelatedData<dim,spacedim>::initialize (const unsigned int n_quadrature_points,
2531  const UpdateFlags flags)
2532  {
2533  if (flags & update_quadrature_points)
2534  this->quadrature_points.resize(n_quadrature_points,
2536 
2537  if (flags & update_JxW_values)
2538  this->JxW_values.resize(n_quadrature_points,
2539  numbers::signaling_nan<double>());
2540 
2541  if (flags & update_jacobians)
2542  this->jacobians.resize(n_quadrature_points,
2544 
2545  if (flags & update_jacobian_grads)
2546  this->jacobian_grads.resize(n_quadrature_points,
2548 
2550  this->jacobian_pushed_forward_grads.resize(n_quadrature_points,
2552 
2553  if (flags & update_jacobian_2nd_derivatives)
2554  this->jacobian_2nd_derivatives.resize(n_quadrature_points,
2556 
2558  this->jacobian_pushed_forward_2nd_derivatives.resize(n_quadrature_points,
2560 
2561  if (flags & update_jacobian_3rd_derivatives)
2562  this->jacobian_3rd_derivatives.resize(n_quadrature_points);
2563 
2565  this->jacobian_pushed_forward_3rd_derivatives.resize(n_quadrature_points,
2567 
2568  if (flags & update_inverse_jacobians)
2569  this->inverse_jacobians.resize(n_quadrature_points,
2571 
2572  if (flags & update_boundary_forms)
2573  this->boundary_forms.resize(n_quadrature_points,
2575 
2576  if (flags & update_normal_vectors)
2577  this->normal_vectors.resize(n_quadrature_points,
2579  }
2580 
2581 
2582 
2583  template <int dim, int spacedim>
2584  std::size_t
2586  {
2587  return (MemoryConsumption::memory_consumption (JxW_values) +
2589  MemoryConsumption::memory_consumption (jacobian_grads) +
2590  MemoryConsumption::memory_consumption (jacobian_pushed_forward_grads) +
2591  MemoryConsumption::memory_consumption (jacobian_2nd_derivatives) +
2592  MemoryConsumption::memory_consumption (jacobian_pushed_forward_2nd_derivatives) +
2593  MemoryConsumption::memory_consumption (jacobian_3rd_derivatives) +
2594  MemoryConsumption::memory_consumption (jacobian_pushed_forward_3rd_derivatives) +
2595  MemoryConsumption::memory_consumption (inverse_jacobians) +
2596  MemoryConsumption::memory_consumption (quadrature_points) +
2597  MemoryConsumption::memory_consumption (normal_vectors) +
2598  MemoryConsumption::memory_consumption (boundary_forms));
2599  }
2600 
2601 
2602 
2603 
2604  template <int dim, int spacedim>
2605  void
2606  FiniteElementRelatedData<dim,spacedim>::initialize (const unsigned int n_quadrature_points,
2607  const FiniteElement<dim,spacedim> &fe,
2608  const UpdateFlags flags)
2609  {
2610  // initialize the table mapping from shape function number to
2611  // the rows in the tables storing the data by shape function and
2612  // nonzero component
2613  this->shape_function_to_row_table
2614  = make_shape_function_to_row_table (fe);
2615 
2616  // count the total number of non-zero components accumulated
2617  // over all shape functions
2618  unsigned int n_nonzero_shape_components = 0;
2619  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
2620  n_nonzero_shape_components += fe.n_nonzero_components (i);
2621  Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
2622  ExcInternalError());
2623 
2624  // with the number of rows now
2625  // known, initialize those fields
2626  // that we will need to their
2627  // correct size
2628  if (flags & update_values)
2629  {
2630  this->shape_values.reinit(n_nonzero_shape_components,
2631  n_quadrature_points);
2632  this->shape_values.fill(numbers::signaling_nan<double>());
2633  }
2634 
2635  if (flags & update_gradients)
2636  {
2637  this->shape_gradients.reinit(n_nonzero_shape_components,
2638  n_quadrature_points);
2639  this->shape_gradients.fill (numbers::signaling_nan<Tensor<1,spacedim> >());
2640  }
2641 
2642  if (flags & update_hessians)
2643  {
2644  this->shape_hessians.reinit(n_nonzero_shape_components,
2645  n_quadrature_points);
2646  this->shape_hessians.fill (numbers::signaling_nan<Tensor<2,spacedim> >());
2647  }
2648 
2649  if (flags & update_3rd_derivatives)
2650  {
2651  this->shape_3rd_derivatives.reinit(n_nonzero_shape_components,
2652  n_quadrature_points);
2653  this->shape_3rd_derivatives.fill (numbers::signaling_nan<Tensor<3,spacedim> >());
2654  }
2655  }
2656 
2657 
2658 
2659 
2660  template <int dim, int spacedim>
2661  std::size_t
2663  {
2664  return (MemoryConsumption::memory_consumption (shape_values) +
2665  MemoryConsumption::memory_consumption (shape_gradients) +
2666  MemoryConsumption::memory_consumption (shape_hessians) +
2667  MemoryConsumption::memory_consumption (shape_3rd_derivatives) +
2668  MemoryConsumption::memory_consumption (shape_function_to_row_table));
2669  }
2670  }
2671 }
2672 
2673 
2674 
2675 /*------------------------------- FEValuesBase ---------------------------*/
2676 
2677 
2678 template <int dim, int spacedim>
2679 FEValuesBase<dim,spacedim>::FEValuesBase (const unsigned int n_q_points,
2680  const unsigned int dofs_per_cell,
2681  const UpdateFlags flags,
2683  const FiniteElement<dim,spacedim> &fe)
2684  :
2685  n_quadrature_points (n_q_points),
2686  dofs_per_cell (dofs_per_cell),
2687  mapping(&mapping, typeid(*this).name()),
2688  fe(&fe, typeid(*this).name()),
2689  fe_values_views_cache (*this)
2690 {
2691  Assert (n_q_points > 0,
2692  ExcMessage ("There is nothing useful you can do with an FEValues "
2693  "object when using a quadrature formula with zero "
2694  "quadrature points!"));
2695  this->update_flags = flags;
2696 }
2697 
2698 
2699 
2700 template <int dim, int spacedim>
2702 {
2703  tria_listener_refinement.disconnect ();
2704  tria_listener_mesh_transform.disconnect ();
2705 }
2706 
2707 
2708 
2709 namespace internal
2710 {
2711  // put shape function part of get_function_xxx methods into separate
2712  // internal functions. this allows us to reuse the same code for several
2713  // functions (e.g. both the versions with and without indices) as well as
2714  // the same code for gradients and Hessians. Moreover, this speeds up
2715  // compilation and reduces the size of the final file since all the
2716  // different global vectors get channeled through the same code.
2717 
2718  template <typename Number, typename Number2>
2719  void
2720  do_function_values (const Number2 *dof_values_ptr,
2721  const ::Table<2,double> &shape_values,
2722  std::vector<Number> &values)
2723  {
2724  // scalar finite elements, so shape_values.size() == dofs_per_cell
2725  const unsigned int dofs_per_cell = shape_values.n_rows();
2726  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
2727  shape_values.n_cols() : values.size();
2728  AssertDimension(values.size(), n_quadrature_points);
2729 
2730  // initialize with zero
2731  std::fill_n (values.begin(), n_quadrature_points,
2733 
2734  // add up contributions of trial functions. note that here we deal with
2735  // scalar finite elements, so no need to check for non-primitivity of
2736  // shape functions. in order to increase the speed of this function, we
2737  // directly access the data in the shape_values array, and increment
2738  // pointers for accessing the data. this saves some lookup time and
2739  // indexing. moreover, the order of the loops is such that we can access
2740  // the shape_values data stored contiguously
2741  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
2742  {
2743  const Number2 value = dof_values_ptr[shape_func];
2744  // For auto-differentiable numbers, the fact that a DoF value is zero
2745  // does not imply that its derivatives are zero as well. So we
2746  // can't filter by value for these number types.
2748  if (value == ::internal::NumberType<Number2>::value(0.0))
2749  continue;
2750 
2751  const double *shape_value_ptr = &shape_values(shape_func, 0);
2752  for (unsigned int point=0; point<n_quadrature_points; ++point)
2753  values[point] += value * (*shape_value_ptr++);
2754  }
2755  }
2756 
2757  template <int dim, int spacedim, typename VectorType>
2758  void
2759  do_function_values (const typename VectorType::value_type *dof_values_ptr,
2760  const ::Table<2,double> &shape_values,
2761  const FiniteElement<dim,spacedim> &fe,
2762  const std::vector<unsigned int> &shape_function_to_row_table,
2763  ArrayView<VectorType> values,
2764  const bool quadrature_points_fastest = false,
2765  const unsigned int component_multiple = 1)
2766  {
2767  typedef typename VectorType::value_type Number;
2768  // initialize with zero
2769  for (unsigned int i=0; i<values.size(); ++i)
2770  std::fill_n (values[i].begin(), values[i].size(),
2771  typename VectorType::value_type());
2772 
2773  // see if there the current cell has DoFs at all, and if not
2774  // then there is nothing else to do.
2775  const unsigned int dofs_per_cell = fe.dofs_per_cell;
2776  if (dofs_per_cell == 0)
2777  return;
2778 
2779  const unsigned int n_quadrature_points = shape_values.n_cols();
2780  const unsigned int n_components = fe.n_components();
2781 
2782  // Assert that we can write all components into the result vectors
2783  const unsigned result_components = n_components * component_multiple;
2784  (void)result_components;
2785  if (quadrature_points_fastest)
2786  {
2787  AssertDimension(values.size(), result_components);
2788  for (unsigned int i=0; i<values.size(); ++i)
2789  AssertDimension (values[i].size(), n_quadrature_points);
2790  }
2791  else
2792  {
2793  AssertDimension(values.size(), n_quadrature_points);
2794  for (unsigned int i=0; i<values.size(); ++i)
2795  AssertDimension (values[i].size(), result_components);
2796  }
2797 
2798  // add up contributions of trial functions. now check whether the shape
2799  // function is primitive or not. if it is, then set its only non-zero
2800  // component, otherwise loop over components
2801  for (unsigned int mc = 0; mc < component_multiple; ++mc)
2802  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
2803  {
2804  const Number &value = dof_values_ptr[shape_func+mc*dofs_per_cell];
2805  // For auto-differentiable numbers, the fact that a DoF value is zero
2806  // does not imply that its derivatives are zero as well. So we
2807  // can't filter by value for these number types.
2809  if (value == ::internal::NumberType<Number>::value(0.0))
2810  continue;
2811 
2812  if (fe.is_primitive(shape_func))
2813  {
2814  const unsigned int comp =
2815  fe.system_to_component_index(shape_func).first
2816  + mc * n_components;
2817  const unsigned int
2818  row = shape_function_to_row_table[shape_func*n_components+comp];
2819 
2820  const double *shape_value_ptr = &shape_values(row, 0);
2821 
2822  if (quadrature_points_fastest)
2823  {
2824  VectorType &values_comp = values[comp];
2825  for (unsigned int point=0; point<n_quadrature_points; ++point)
2826  values_comp[point] += value * (*shape_value_ptr++);
2827  }
2828  else
2829  for (unsigned int point=0; point<n_quadrature_points; ++point)
2830  values[point][comp] += value * (*shape_value_ptr++);
2831  }
2832  else
2833  for (unsigned int c=0; c<n_components; ++c)
2834  {
2835  if (fe.get_nonzero_components(shape_func)[c] == false)
2836  continue;
2837 
2838  const unsigned int
2839  row = shape_function_to_row_table[shape_func*n_components+c];
2840 
2841  const double *shape_value_ptr = &shape_values(row, 0);
2842  const unsigned int comp = c + mc * n_components;
2843 
2844  if (quadrature_points_fastest)
2845  {
2846  VectorType &values_comp = values[comp];
2847  for (unsigned int point=0; point<n_quadrature_points;
2848  ++point)
2849  values_comp[point] += value * (*shape_value_ptr++);
2850  }
2851  else
2852  for (unsigned int point=0; point<n_quadrature_points; ++point)
2853  values[point][comp] += value * (*shape_value_ptr++);
2854  }
2855  }
2856  }
2857 
2858  // use the same implementation for gradients and Hessians, distinguish them
2859  // by the rank of the tensors
2860  template <int order, int spacedim, typename Number>
2861  void
2862  do_function_derivatives (const Number *dof_values_ptr,
2863  const ::Table<2,Tensor<order,spacedim> > &shape_derivatives,
2864  std::vector<Tensor<order,spacedim,Number> > &derivatives)
2865  {
2866  const unsigned int dofs_per_cell = shape_derivatives.size()[0];
2867  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
2868  shape_derivatives[0].size() : derivatives.size();
2869  AssertDimension(derivatives.size(), n_quadrature_points);
2870 
2871  // initialize with zero
2872  std::fill_n (derivatives.begin(), n_quadrature_points, Tensor<order,spacedim,Number>());
2873 
2874  // add up contributions of trial functions. note that here we deal with
2875  // scalar finite elements, so no need to check for non-primitivity of
2876  // shape functions. in order to increase the speed of this function, we
2877  // directly access the data in the shape_gradients/hessians array, and
2878  // increment pointers for accessing the data. this saves some lookup time
2879  // and indexing. moreover, the order of the loops is such that we can
2880  // access the shape_gradients/hessians data stored contiguously
2881  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
2882  {
2883  const Number &value = dof_values_ptr[shape_func];
2884  // For auto-differentiable numbers, the fact that a DoF value is zero
2885  // does not imply that its derivatives are zero as well. So we
2886  // can't filter by value for these number types.
2888  if (value == ::internal::NumberType<Number>::value(0.0))
2889  continue;
2890 
2891  const Tensor<order,spacedim> *shape_derivative_ptr
2892  = &shape_derivatives[shape_func][0];
2893  for (unsigned int point=0; point<n_quadrature_points; ++point)
2894  derivatives[point] += value *
2895  ::Tensor<order,spacedim,Number>(*shape_derivative_ptr++);
2896  }
2897  }
2898 
2899  template <int order, int dim, int spacedim, typename Number>
2900  void
2901  do_function_derivatives (const Number *dof_values_ptr,
2902  const ::Table<2,Tensor<order,spacedim> > &shape_derivatives,
2903  const FiniteElement<dim,spacedim> &fe,
2904  const std::vector<unsigned int> &shape_function_to_row_table,
2905  ArrayView<std::vector<Tensor<order,spacedim,Number> > > derivatives,
2906  const bool quadrature_points_fastest = false,
2907  const unsigned int component_multiple = 1)
2908  {
2909  // initialize with zero
2910  for (unsigned int i=0; i<derivatives.size(); ++i)
2911  std::fill_n (derivatives[i].begin(), derivatives[i].size(),
2913 
2914  // see if there the current cell has DoFs at all, and if not
2915  // then there is nothing else to do.
2916  const unsigned int dofs_per_cell = fe.dofs_per_cell;
2917  if (dofs_per_cell == 0)
2918  return;
2919 
2920 
2921  const unsigned int n_quadrature_points = shape_derivatives[0].size();
2922  const unsigned int n_components = fe.n_components();
2923 
2924  // Assert that we can write all components into the result vectors
2925  const unsigned result_components = n_components * component_multiple;
2926  (void)result_components;
2927  if (quadrature_points_fastest)
2928  {
2929  AssertDimension(derivatives.size(), result_components);
2930  for (unsigned int i=0; i<derivatives.size(); ++i)
2931  AssertDimension (derivatives[i].size(), n_quadrature_points);
2932  }
2933  else
2934  {
2935  AssertDimension(derivatives.size(), n_quadrature_points);
2936  for (unsigned int i=0; i<derivatives.size(); ++i)
2937  AssertDimension (derivatives[i].size(), result_components);
2938  }
2939 
2940  // add up contributions of trial functions. now check whether the shape
2941  // function is primitive or not. if it is, then set its only non-zero
2942  // component, otherwise loop over components
2943  for (unsigned int mc = 0; mc < component_multiple; ++mc)
2944  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
2945  {
2946  const Number &value = dof_values_ptr[shape_func+mc*dofs_per_cell];
2947  // For auto-differentiable numbers, the fact that a DoF value is zero
2948  // does not imply that its derivatives are zero as well. So we
2949  // can't filter by value for these number types.
2951  if (value == ::internal::NumberType<Number>::value(0.0))
2952  continue;
2953 
2954  if (fe.is_primitive(shape_func))
2955  {
2956  const unsigned int comp =
2957  fe.system_to_component_index(shape_func).first
2958  + mc * n_components;
2959  const unsigned int
2960  row = shape_function_to_row_table[shape_func*n_components+comp];
2961 
2962  const Tensor<order,spacedim> *shape_derivative_ptr =
2963  &shape_derivatives[row][0];
2964 
2965  if (quadrature_points_fastest)
2966  for (unsigned int point=0; point<n_quadrature_points; ++point)
2967  derivatives[comp][point] += value *
2968  ::Tensor<order,spacedim,Number>(*shape_derivative_ptr++);
2969  else
2970  for (unsigned int point=0; point<n_quadrature_points; ++point)
2971  derivatives[point][comp] += value *
2972  ::Tensor<order,spacedim,Number>(*shape_derivative_ptr++);
2973  }
2974  else
2975  for (unsigned int c=0; c<n_components; ++c)
2976  {
2977  if (fe.get_nonzero_components(shape_func)[c] == false)
2978  continue;
2979 
2980  const unsigned int
2981  row = shape_function_to_row_table[shape_func*n_components+c];
2982 
2983  const Tensor<order,spacedim> *shape_derivative_ptr =
2984  &shape_derivatives[row][0];
2985  const unsigned int comp = c + mc * n_components;
2986 
2987  if (quadrature_points_fastest)
2988  for (unsigned int point=0; point<n_quadrature_points; ++point)
2989  derivatives[comp][point] += value *
2990  ::Tensor<order,spacedim,Number>(*shape_derivative_ptr++);
2991  else
2992  for (unsigned int point=0; point<n_quadrature_points; ++point)
2993  derivatives[point][comp] += value *
2994  ::Tensor<order,spacedim,Number>(*shape_derivative_ptr++);
2995  }
2996  }
2997  }
2998 
2999  template <int spacedim, typename Number, typename Number2>
3000  void
3001  do_function_laplacians (const Number2 *dof_values_ptr,
3002  const ::Table<2,Tensor<2,spacedim> > &shape_hessians,
3003  std::vector<Number> &laplacians)
3004  {
3005  const unsigned int dofs_per_cell = shape_hessians.size()[0];
3006  const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
3007  shape_hessians[0].size() : laplacians.size();
3008  AssertDimension(laplacians.size(), n_quadrature_points);
3009 
3010  // initialize with zero
3011  std::fill_n (laplacians.begin(), n_quadrature_points,
3013 
3014  // add up contributions of trial functions. note that here we deal with
3015  // scalar finite elements and also note that the Laplacian is
3016  // the trace of the Hessian.
3017  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
3018  {
3019  const Number2 value = dof_values_ptr[shape_func];
3020  // For auto-differentiable numbers, the fact that a DoF value is zero
3021  // does not imply that its derivatives are zero as well. So we
3022  // can't filter by value for these number types.
3024  if (value == ::internal::NumberType<Number2>::value(0.0))
3025  continue;
3026 
3027  const Tensor<2,spacedim> *shape_hessian_ptr
3028  = &shape_hessians[shape_func][0];
3029  for (unsigned int point=0; point<n_quadrature_points; ++point)
3030  laplacians[point] += value * trace(*shape_hessian_ptr++);
3031  }
3032  }
3033 
3034  template <int dim, int spacedim, typename VectorType, typename Number>
3035  void
3036  do_function_laplacians (const Number *dof_values_ptr,
3037  const ::Table<2,Tensor<2,spacedim> > &shape_hessians,
3038  const FiniteElement<dim,spacedim> &fe,
3039  const std::vector<unsigned int> &shape_function_to_row_table,
3040  std::vector<VectorType> &laplacians,
3041  const bool quadrature_points_fastest = false,
3042  const unsigned int component_multiple = 1)
3043  {
3044  // initialize with zero
3045  for (unsigned int i=0; i<laplacians.size(); ++i)
3046  std::fill_n (laplacians[i].begin(), laplacians[i].size(),
3047  typename VectorType::value_type());
3048 
3049  // see if there the current cell has DoFs at all, and if not
3050  // then there is nothing else to do.
3051  const unsigned int dofs_per_cell = fe.dofs_per_cell;
3052  if (dofs_per_cell == 0)
3053  return;
3054 
3055 
3056  const unsigned int n_quadrature_points = shape_hessians[0].size();
3057  const unsigned int n_components = fe.n_components();
3058 
3059  // Assert that we can write all components into the result vectors
3060  const unsigned result_components = n_components * component_multiple;
3061  (void)result_components;
3062  if (quadrature_points_fastest)
3063  {
3064  AssertDimension(laplacians.size(), result_components);
3065  for (unsigned int i=0; i<laplacians.size(); ++i)
3066  AssertDimension (laplacians[i].size(), n_quadrature_points);
3067  }
3068  else
3069  {
3070  AssertDimension(laplacians.size(), n_quadrature_points);
3071  for (unsigned int i=0; i<laplacians.size(); ++i)
3072  AssertDimension (laplacians[i].size(), result_components);
3073  }
3074 
3075  // add up contributions of trial functions. now check whether the shape
3076  // function is primitive or not. if it is, then set its only non-zero
3077  // component, otherwise loop over components
3078  for (unsigned int mc = 0; mc < component_multiple; ++mc)
3079  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
3080  {
3081  const Number &value = dof_values_ptr[shape_func+mc*dofs_per_cell];
3082  // For auto-differentiable numbers, the fact that a DoF value is zero
3083  // does not imply that its derivatives are zero as well. So we
3084  // can't filter by value for these number types.
3086  if (value == ::internal::NumberType<Number>::value(0.0))
3087  continue;
3088 
3089  if (fe.is_primitive(shape_func))
3090  {
3091  const unsigned int comp =
3092  fe.system_to_component_index(shape_func).first
3093  + mc * n_components;
3094  const unsigned int
3095  row = shape_function_to_row_table[shape_func*n_components+comp];
3096 
3097  const Tensor<2,spacedim> *shape_hessian_ptr =
3098  &shape_hessians[row][0];
3099  if (quadrature_points_fastest)
3100  {
3101  VectorType &laplacians_comp = laplacians[comp];
3102  for (unsigned int point=0; point<n_quadrature_points; ++point)
3103  laplacians_comp[point] += value * trace(*shape_hessian_ptr++);
3104  }
3105  else
3106  for (unsigned int point=0; point<n_quadrature_points; ++point)
3107  laplacians[point][comp] += value * trace(*shape_hessian_ptr++);
3108  }
3109  else
3110  for (unsigned int c=0; c<n_components; ++c)
3111  {
3112  if (fe.get_nonzero_components(shape_func)[c] == false)
3113  continue;
3114 
3115  const unsigned int
3116  row = shape_function_to_row_table[shape_func*n_components+c];
3117 
3118  const Tensor<2,spacedim> *shape_hessian_ptr =
3119  &shape_hessians[row][0];
3120  const unsigned int comp = c + mc * n_components;
3121 
3122  if (quadrature_points_fastest)
3123  {
3124  VectorType &laplacians_comp = laplacians[comp];
3125  for (unsigned int point=0; point<n_quadrature_points;
3126  ++point)
3127  laplacians_comp[point] += value * trace(*shape_hessian_ptr++);
3128  }
3129  else
3130  for (unsigned int point=0; point<n_quadrature_points; ++point)
3131  laplacians[point][comp] += value * trace(*shape_hessian_ptr++);
3132  }
3133  }
3134  }
3135 }
3136 
3137 
3138 
3139 template <int dim, int spacedim>
3140 template <class InputVector>
3142  const InputVector &fe_function,
3143  std::vector<typename InputVector::value_type> &values) const
3144 {
3145  typedef typename InputVector::value_type Number;
3146  Assert (this->update_flags & update_values,
3147  ExcAccessToUninitializedField("update_values"));
3148  AssertDimension (fe->n_components(), 1);
3149  Assert (present_cell.get() != nullptr,
3150  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3151  AssertDimension (fe_function.size(),
3152  present_cell->n_dofs_for_dof_handler());
3153 
3154  // get function values of dofs on this cell
3155  Vector<Number> dof_values (dofs_per_cell);
3156  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3157  internal::do_function_values (dof_values.begin(), this->finite_element_output.shape_values,
3158  values);
3159 }
3160 
3161 
3162 
3163 template <int dim, int spacedim>
3164 template <class InputVector>
3166  const InputVector &fe_function,
3167  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3168  std::vector<typename InputVector::value_type> &values) const
3169 {
3170  typedef typename InputVector::value_type Number;
3171  Assert (this->update_flags & update_values,
3172  ExcAccessToUninitializedField("update_values"));
3173  AssertDimension (fe->n_components(), 1);
3174  AssertDimension (indices.size(), dofs_per_cell);
3175 
3176  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3177  for (unsigned int i=0; i<dofs_per_cell; ++i)
3178  dof_values[i] = get_vector_element (fe_function, indices[i]);
3179  internal::do_function_values(dof_values.data(), this->finite_element_output.shape_values, values);
3180 }
3181 
3182 
3183 
3184 template <int dim, int spacedim>
3185 template <class InputVector>
3187  const InputVector &fe_function,
3188  std::vector<Vector<typename InputVector::value_type> > &values) const
3189 {
3190  typedef typename InputVector::value_type Number;
3191  Assert (present_cell.get() != nullptr,
3192  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3193 
3194  Assert (this->update_flags & update_values,
3195  ExcAccessToUninitializedField("update_values"));
3196  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3197 
3198  // get function values of dofs on this cell
3199  Vector<Number> dof_values (dofs_per_cell);
3200  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3201  internal::do_function_values(dof_values.begin(),
3202  this->finite_element_output.shape_values,
3203  *fe,
3204  this->finite_element_output.shape_function_to_row_table,
3205  make_array_view(values.begin(), values.end()));
3206 }
3207 
3208 
3209 
3210 template <int dim, int spacedim>
3211 template <class InputVector>
3213  const InputVector &fe_function,
3214  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3215  std::vector<Vector<typename InputVector::value_type> > &values) const
3216 {
3217  typedef typename InputVector::value_type Number;
3218  // Size of indices must be a multiple of dofs_per_cell such that an integer
3219  // number of function values is generated in each point.
3220  Assert (indices.size() % dofs_per_cell == 0,
3221  ExcNotMultiple(indices.size(), dofs_per_cell));
3222  Assert (this->update_flags & update_values,
3223  ExcAccessToUninitializedField("update_values"));
3224 
3225  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3226  for (unsigned int i=0; i<dofs_per_cell; ++i)
3227  dof_values[i] = get_vector_element (fe_function, indices[i]);
3228  internal::do_function_values(dof_values.data(),
3229  this->finite_element_output.shape_values,
3230  *fe,
3231  this->finite_element_output.shape_function_to_row_table,
3232  make_array_view(values.begin(), values.end()),
3233  false,
3234  indices.size()/dofs_per_cell);
3235 }
3236 
3237 
3238 
3239 template <int dim, int spacedim>
3240 template <class InputVector>
3242  const InputVector &fe_function,
3243  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3244  VectorSlice<std::vector<std::vector<typename InputVector::value_type> > > values,
3245  bool quadrature_points_fastest) const
3246 {
3247  typedef typename InputVector::value_type Number;
3248  Assert (this->update_flags & update_values,
3249  ExcAccessToUninitializedField("update_values"));
3250 
3251  // Size of indices must be a multiple of dofs_per_cell such that an integer
3252  // number of function values is generated in each point.
3253  Assert (indices.size() % dofs_per_cell == 0,
3254  ExcNotMultiple(indices.size(), dofs_per_cell));
3255 
3256  boost::container::small_vector<Number, 200> dof_values(indices.size());
3257  for (unsigned int i=0; i<indices.size(); ++i)
3258  dof_values[i] = get_vector_element (fe_function, indices[i]);
3259  internal::do_function_values(dof_values.data(),
3260  this->finite_element_output.shape_values,
3261  *fe,
3262  this->finite_element_output.shape_function_to_row_table,
3263  make_array_view(values.begin(), values.end()),
3264  quadrature_points_fastest,
3265  indices.size()/dofs_per_cell);
3266 }
3267 
3268 
3269 
3270 template <int dim, int spacedim>
3271 template <class InputVector>
3272 void
3274  const InputVector &fe_function,
3275  std::vector<Tensor<1,spacedim,typename InputVector::value_type> > &gradients) const
3276 {
3277  typedef typename InputVector::value_type Number;
3278  Assert (this->update_flags & update_gradients,
3279  ExcAccessToUninitializedField("update_gradients"));
3280  AssertDimension (fe->n_components(), 1);
3281  Assert (present_cell.get() != nullptr,
3282  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3283  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3284 
3285  // get function values of dofs on this cell
3286  Vector<Number> dof_values (dofs_per_cell);
3287  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3288  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_gradients,
3289  gradients);
3290 }
3291 
3292 
3293 
3294 template <int dim, int spacedim>
3295 template <class InputVector>
3297  const InputVector &fe_function,
3298  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3299  std::vector<Tensor<1,spacedim,typename InputVector::value_type> > &gradients) const
3300 {
3301  typedef typename InputVector::value_type Number;
3302  Assert (this->update_flags & update_gradients,
3303  ExcAccessToUninitializedField("update_gradients"));
3304  AssertDimension (fe->n_components(), 1);
3305  AssertDimension (indices.size(), dofs_per_cell);
3306 
3307  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3308  for (unsigned int i=0; i<dofs_per_cell; ++i)
3309  dof_values[i] = get_vector_element (fe_function, indices[i]);
3310  internal::do_function_derivatives(dof_values.data(), this->finite_element_output.shape_gradients,
3311  gradients);
3312 }
3313 
3314 
3315 
3316 
3317 template <int dim, int spacedim>
3318 template <class InputVector>
3319 void
3321  const InputVector &fe_function,
3322  std::vector<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > > &gradients) const
3323 {
3324  typedef typename InputVector::value_type Number;
3325  Assert (this->update_flags & update_gradients,
3326  ExcAccessToUninitializedField("update_gradients"));
3327  Assert (present_cell.get() != nullptr,
3328  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3329  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3330 
3331  // get function values of dofs on this cell
3332  Vector<Number> dof_values (dofs_per_cell);
3333  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3334  internal::do_function_derivatives(dof_values.begin(),
3335  this->finite_element_output.shape_gradients,
3336  *fe,
3337  this->finite_element_output.shape_function_to_row_table,
3338  make_array_view(gradients.begin(), gradients.end()));
3339 }
3340 
3341 
3342 
3343 template <int dim, int spacedim>
3344 template <class InputVector>
3346  const InputVector &fe_function,
3347  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3348  VectorSlice<std::vector<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > > > gradients,
3349  bool quadrature_points_fastest) const
3350 {
3351  typedef typename InputVector::value_type Number;
3352  // Size of indices must be a multiple of dofs_per_cell such that an integer
3353  // number of function values is generated in each point.
3354  Assert (indices.size() % dofs_per_cell == 0,
3355  ExcNotMultiple(indices.size(), dofs_per_cell));
3356  Assert (this->update_flags & update_gradients,
3357  ExcAccessToUninitializedField("update_gradients"));
3358 
3359  boost::container::small_vector<Number, 200> dof_values(indices.size());
3360  for (unsigned int i=0; i<indices.size(); ++i)
3361  dof_values[i] = get_vector_element (fe_function, indices[i]);
3362  internal::do_function_derivatives(dof_values.data(),
3363  this->finite_element_output.shape_gradients,
3364  *fe,
3365  this->finite_element_output.shape_function_to_row_table,
3366  make_array_view(gradients.begin(), gradients.end()),
3367  quadrature_points_fastest,
3368  indices.size()/dofs_per_cell);
3369 }
3370 
3371 
3372 
3373 template <int dim, int spacedim>
3374 template <class InputVector>
3375 void
3377 get_function_hessians (const InputVector &fe_function,
3378  std::vector<Tensor<2,spacedim,typename InputVector::value_type> > &hessians) const
3379 {
3380  typedef typename InputVector::value_type Number;
3381  AssertDimension (fe->n_components(), 1);
3382  Assert (this->update_flags & update_hessians,
3383  ExcAccessToUninitializedField("update_hessians"));
3384  Assert (present_cell.get() != nullptr,
3385  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3386  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3387 
3388  // get function values of dofs on this cell
3389  Vector<Number> dof_values (dofs_per_cell);
3390  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3391  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_hessians,
3392  hessians);
3393 }
3394 
3395 
3396 
3397 template <int dim, int spacedim>
3398 template <class InputVector>
3400  const InputVector &fe_function,
3401  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3402  std::vector<Tensor<2,spacedim,typename InputVector::value_type> > &hessians) const
3403 {
3404  typedef typename InputVector::value_type Number;
3405  Assert (this->update_flags & update_hessians,
3406  ExcAccessToUninitializedField("update_hessians"));
3407  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3408  AssertDimension (indices.size(), dofs_per_cell);
3409 
3410  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3411  for (unsigned int i=0; i<dofs_per_cell; ++i)
3412  dof_values[i] = get_vector_element (fe_function, indices[i]);
3413  internal::do_function_derivatives(dof_values.data(), this->finite_element_output.shape_hessians,
3414  hessians);
3415 }
3416 
3417 
3418 
3419 
3420 template <int dim, int spacedim>
3421 template <class InputVector>
3422 void
3424 get_function_hessians (const InputVector &fe_function,
3425  std::vector<std::vector<Tensor<2,spacedim,typename InputVector::value_type> > > &hessians,
3426  bool quadrature_points_fastest) const
3427 {
3428  typedef typename InputVector::value_type Number;
3429  Assert (this->update_flags & update_hessians,
3430  ExcAccessToUninitializedField("update_hessians"));
3431  Assert (present_cell.get() != nullptr,
3432  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3433  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3434 
3435  // get function values of dofs on this cell
3436  Vector<Number> dof_values (dofs_per_cell);
3437  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3438  internal::do_function_derivatives(dof_values.begin(),
3439  this->finite_element_output.shape_hessians,
3440  *fe,
3441  this->finite_element_output.shape_function_to_row_table,
3442  make_array_view(hessians.begin(), hessians.end()),
3443  quadrature_points_fastest);
3444 }
3445 
3446 
3447 
3448 template <int dim, int spacedim>
3449 template <class InputVector>
3451  const InputVector &fe_function,
3452  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3453  VectorSlice<std::vector<std::vector<Tensor<2,spacedim,typename InputVector::value_type> > > > hessians,
3454  bool quadrature_points_fastest) const
3455 {
3456  typedef typename InputVector::value_type Number;
3457  Assert (this->update_flags & update_hessians,
3458  ExcAccessToUninitializedField("update_hessians"));
3459  Assert (indices.size() % dofs_per_cell == 0,
3460  ExcNotMultiple(indices.size(), dofs_per_cell));
3461 
3462  boost::container::small_vector<Number, 200> dof_values(indices.size());
3463  for (unsigned int i=0; i<indices.size(); ++i)
3464  dof_values[i] = get_vector_element (fe_function, indices[i]);
3465  internal::do_function_derivatives(dof_values.data(),
3466  this->finite_element_output.shape_hessians,
3467  *fe,
3468  this->finite_element_output.shape_function_to_row_table,
3469  make_array_view(hessians.begin(), hessians.end()),
3470  quadrature_points_fastest,
3471  indices.size()/dofs_per_cell);
3472 }
3473 
3474 
3475 
3476 template <int dim, int spacedim>
3477 template <class InputVector>
3479  const InputVector &fe_function,
3480  std::vector<typename InputVector::value_type> &laplacians) const
3481 {
3482  typedef typename InputVector::value_type Number;
3483  Assert (this->update_flags & update_hessians,
3484  ExcAccessToUninitializedField("update_hessians"));
3485  AssertDimension (fe->n_components(), 1);
3486  Assert (present_cell.get() != nullptr,
3487  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3488  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3489 
3490  // get function values of dofs on this cell
3491  Vector<Number> dof_values (dofs_per_cell);
3492  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3493  internal::do_function_laplacians(dof_values.begin(), this->finite_element_output.shape_hessians,
3494  laplacians);
3495 }
3496 
3497 
3498 
3499 template <int dim, int spacedim>
3500 template <class InputVector>
3502  const InputVector &fe_function,
3503  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3504  std::vector<typename InputVector::value_type> &laplacians) const
3505 {
3506  typedef typename InputVector::value_type Number;
3507  Assert (this->update_flags & update_hessians,
3508  ExcAccessToUninitializedField("update_hessians"));
3509  AssertDimension (fe->n_components(), 1);
3510  AssertDimension (indices.size(), dofs_per_cell);
3511 
3512  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3513  for (unsigned int i=0; i<dofs_per_cell; ++i)
3514  dof_values[i] = get_vector_element (fe_function, indices[i]);
3515  internal::do_function_laplacians(dof_values.data(), this->finite_element_output.shape_hessians,
3516  laplacians);
3517 }
3518 
3519 
3520 
3521 template <int dim, int spacedim>
3522 template <class InputVector>
3524  const InputVector &fe_function,
3525  std::vector<Vector<typename InputVector::value_type> > &laplacians) const
3526 {
3527  typedef typename InputVector::value_type Number;
3528  Assert (present_cell.get() != nullptr,
3529  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3530  Assert (this->update_flags & update_hessians,
3531  ExcAccessToUninitializedField("update_hessians"));
3532  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3533 
3534  // get function values of dofs on this cell
3535  Vector<Number> dof_values (dofs_per_cell);
3536  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3537  internal::do_function_laplacians(dof_values.begin(), this->finite_element_output.shape_hessians,
3538  *fe, this->finite_element_output.shape_function_to_row_table,
3539  laplacians);
3540 }
3541 
3542 
3543 
3544 template <int dim, int spacedim>
3545 template <class InputVector>
3547  const InputVector &fe_function,
3548  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3549  std::vector<Vector<typename InputVector::value_type> > &laplacians) const
3550 {
3551  typedef typename InputVector::value_type Number;
3552  // Size of indices must be a multiple of dofs_per_cell such that an integer
3553  // number of function values is generated in each point.
3554  Assert (indices.size() % dofs_per_cell == 0,
3555  ExcNotMultiple(indices.size(), dofs_per_cell));
3556  Assert (this->update_flags & update_hessians,
3557  ExcAccessToUninitializedField("update_hessians"));
3558 
3559  boost::container::small_vector<Number, 200> dof_values(indices.size());
3560  for (unsigned int i=0; i<indices.size(); ++i)
3561  dof_values[i] = get_vector_element (fe_function, indices[i]);
3562  internal::do_function_laplacians(dof_values.data(), this->finite_element_output.shape_hessians,
3563  *fe, this->finite_element_output.shape_function_to_row_table,
3564  laplacians, false,
3565  indices.size()/dofs_per_cell);
3566 }
3567 
3568 
3569 
3570 template <int dim, int spacedim>
3571 template <class InputVector>
3573  const InputVector &fe_function,
3574  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3575  std::vector<std::vector<typename InputVector::value_type> > &laplacians,
3576  bool quadrature_points_fastest) const
3577 {
3578  typedef typename InputVector::value_type Number;
3579  Assert (indices.size() % dofs_per_cell == 0,
3580  ExcNotMultiple(indices.size(), dofs_per_cell));
3581  Assert (this->update_flags & update_hessians,
3582  ExcAccessToUninitializedField("update_hessians"));
3583 
3584  boost::container::small_vector<Number, 200> dof_values(indices.size());
3585  for (unsigned int i=0; i<indices.size(); ++i)
3586  dof_values[i] = get_vector_element (fe_function, indices[i]);
3587  internal::do_function_laplacians(dof_values.data(), this->finite_element_output.shape_hessians,
3588  *fe, this->finite_element_output.shape_function_to_row_table,
3589  laplacians, quadrature_points_fastest,
3590  indices.size()/dofs_per_cell);
3591 }
3592 
3593 
3594 
3595 template <int dim, int spacedim>
3596 template <class InputVector>
3597 void
3599 get_function_third_derivatives (const InputVector &fe_function,
3600  std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const
3601 {
3602  typedef typename InputVector::value_type Number;
3603  AssertDimension (fe->n_components(), 1);
3604  Assert (this->update_flags & update_3rd_derivatives,
3605  ExcAccessToUninitializedField("update_3rd_derivatives"));
3606  Assert (present_cell.get() != nullptr,
3607  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3608  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3609 
3610  // get function values of dofs on this cell
3611  Vector<Number> dof_values (dofs_per_cell);
3612  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3613  internal::do_function_derivatives(dof_values.begin(), this->finite_element_output.shape_3rd_derivatives,
3614  third_derivatives);
3615 }
3616 
3617 
3618 
3619 template <int dim, int spacedim>
3620 template <class InputVector>
3622  const InputVector &fe_function,
3623  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3624  std::vector<Tensor<3,spacedim,typename InputVector::value_type> > &third_derivatives) const
3625 {
3626  typedef typename InputVector::value_type Number;
3627  Assert (this->update_flags & update_3rd_derivatives,
3628  ExcAccessToUninitializedField("update_3rd_derivatives"));
3629  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3630  AssertDimension (indices.size(), dofs_per_cell);
3631 
3632  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
3633  for (unsigned int i=0; i<dofs_per_cell; ++i)
3634  dof_values[i] = get_vector_element (fe_function, indices[i]);
3635  internal::do_function_derivatives(dof_values.data(), this->finite_element_output.shape_3rd_derivatives,
3636  third_derivatives);
3637 }
3638 
3639 
3640 
3641 
3642 template <int dim, int spacedim>
3643 template <class InputVector>
3644 void
3646 get_function_third_derivatives (const InputVector &fe_function,
3647  std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > &third_derivatives,
3648  bool quadrature_points_fastest) const
3649 {
3650  typedef typename InputVector::value_type Number;
3651  Assert (this->update_flags & update_3rd_derivatives,
3652  ExcAccessToUninitializedField("update_3rd_derivatives"));
3653  Assert (present_cell.get() != nullptr,
3654  ExcMessage ("FEValues object is not reinit'ed to any cell"));
3655  AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
3656 
3657  // get function values of dofs on this cell
3658  Vector<Number> dof_values (dofs_per_cell);
3659  present_cell->get_interpolated_dof_values(fe_function, dof_values);
3660  internal::do_function_derivatives(dof_values.begin(),
3661  this->finite_element_output.shape_3rd_derivatives,
3662  *fe,
3663  this->finite_element_output.shape_function_to_row_table,
3664  make_array_view(third_derivatives.begin(), third_derivatives.end()),
3665  quadrature_points_fastest);
3666 }
3667 
3668 
3669 
3670 template <int dim, int spacedim>
3671 template <class InputVector>
3673  const InputVector &fe_function,
3674  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
3675  VectorSlice<std::vector<std::vector<Tensor<3,spacedim,typename InputVector::value_type> > > > third_derivatives,
3676  bool quadrature_points_fastest) const
3677 {
3678  typedef typename InputVector::value_type Number;
3679  Assert (this->update_flags & update_3rd_derivatives,
3680  ExcAccessToUninitializedField("update_3rd_derivatives"));
3681  Assert (indices.size() % dofs_per_cell == 0,
3682  ExcNotMultiple(indices.size(), dofs_per_cell));
3683 
3684  boost::container::small_vector<Number, 200> dof_values(indices.size());
3685  for (unsigned int i=0; i<indices.size(); ++i)
3686  dof_values[i] = get_vector_element (fe_function, indices[i]);
3687  internal::do_function_derivatives(dof_values.data(),
3688  this->finite_element_output.shape_3rd_derivatives,
3689  *fe,
3690  this->finite_element_output.shape_function_to_row_table,
3691  make_array_view(third_derivatives.begin(), third_derivatives.end()),
3692  quadrature_points_fastest,
3693  indices.size()/dofs_per_cell);
3694 }
3695 
3696 
3697 
3698 template <int dim, int spacedim>
3701 {
3702  return *present_cell;
3703 }
3704 
3705 
3706 
3707 template <int dim, int spacedim>
3708 const std::vector<Tensor<1,spacedim> > &
3710 {
3711  Assert (this->update_flags & update_normal_vectors,
3712  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_normal_vectors")));
3713  return get_normal_vectors();
3714 }
3715 
3716 
3717 
3718 template <int dim, int spacedim>
3719 const std::vector<Tensor<1,spacedim> > &
3721 {
3722  Assert (this->update_flags & update_normal_vectors,
3723  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_normal_vectors")));
3724 
3725  return this->mapping_output.normal_vectors;
3726 }
3727 
3728 
3729 
3730 template <int dim, int spacedim>
3731 std::size_t
3733 {
3734  return (sizeof(this->update_flags) +
3735  MemoryConsumption::memory_consumption (n_quadrature_points) +
3736  sizeof (cell_similarity) +
3737  MemoryConsumption::memory_consumption (dofs_per_cell) +
3739  MemoryConsumption::memory_consumption (mapping_data) +
3740  MemoryConsumption::memory_consumption (*mapping_data) +
3741  MemoryConsumption::memory_consumption (mapping_output) +
3745  MemoryConsumption::memory_consumption (finite_element_output));
3746 }
3747 
3748 
3749 
3750 template <int dim, int spacedim>
3753 {
3754  // first find out which objects need to be recomputed on each
3755  // cell we visit. this we have to ask the finite element and mapping.
3756  // elements are first since they might require update in mapping
3757  //
3758  // there is no need to iterate since mappings will never require
3759  // the finite element to compute something for them
3760  UpdateFlags flags = update_flags
3761  | fe->requires_update_flags (update_flags);
3762  flags |= mapping->requires_update_flags (flags);
3763 
3764  return flags;
3765 }
3766 
3767 
3768 template <int dim, int spacedim>
3769 void
3771 {
3772  // if there is no present cell, then we shouldn't be
3773  // connected via a signal to a triangulation
3774  Assert (present_cell.get() != nullptr, ExcInternalError());
3775 
3776  // so delete the present cell and
3777  // disconnect from the signal we have with
3778  // it
3779  tria_listener_refinement.disconnect ();
3780  tria_listener_mesh_transform.disconnect ();
3781  present_cell.reset ();
3782 }
3783 
3784 
3785 template <int dim, int spacedim>
3786 void
3789 {
3790  if (present_cell.get() != nullptr)
3791  {
3792  if (&cell->get_triangulation() !=
3793  &present_cell->operator typename Triangulation<dim,spacedim>::cell_iterator()
3794  ->get_triangulation())
3795  {
3796  // the triangulations for the previous cell and the current cell
3797  // do not match. disconnect from the previous triangulation and
3798  // connect to the current one; also invalidate the previous
3799  // cell because we shouldn't be comparing cells from different
3800  // triangulations
3801  invalidate_present_cell();
3802  tria_listener_refinement =
3803  cell->get_triangulation().signals.any_change.connect
3805  std::ref(static_cast<FEValuesBase<dim,spacedim>&>(*this))));
3806  tria_listener_mesh_transform =
3807  cell->get_triangulation().signals.mesh_movement.connect
3809  std::ref(static_cast<FEValuesBase<dim,spacedim>&>(*this))));
3810  }
3811  }
3812  else
3813  {
3814  // if this FEValues has never been set to any cell at all, then
3815  // at least subscribe to the triangulation to get notified of
3816  // changes
3817  tria_listener_refinement =
3818  cell->get_triangulation().signals.post_refinement.connect
3820  std::ref(static_cast<FEValuesBase<dim,spacedim>&>(*this))));
3821  tria_listener_mesh_transform =
3822  cell->get_triangulation().signals.mesh_movement.connect
3824  std::ref(static_cast<FEValuesBase<dim,spacedim>&>(*this))));
3825  }
3826 }
3827 
3828 
3829 template <int dim, int spacedim>
3830 inline
3831 void
3834 {
3835  // Unfortunately, the detection of simple geometries with CellSimilarity is
3836  // sensitive to the first cell detected. When doing this with multiple
3837  // threads, each thread will get its own scratch data object with an
3838  // FEValues object in the implementation framework from late 2013, which is
3839  // initialized to the first cell the thread sees. As this number might
3840  // different between different runs (after all, the tasks are scheduled
3841  // dynamically onto threads), this slight deviation leads to difference in
3842  // roundoff errors that propagate through the program. Therefore, we need to
3843  // disable CellSimilarity in case there is more than one thread in the
3844  // problem. This will likely not affect many MPI test cases as there
3845  // multithreading is disabled on default, but in many other situations
3846  // because we rarely explicitly set the number of threads.
3847  //
3848  // TODO: Is it reasonable to introduce a flag "unsafe" in the constructor of
3849  // FEValues to re-enable this feature?
3850  if (MultithreadInfo::n_threads() > 1)
3851  {
3852  cell_similarity = CellSimilarity::none;
3853  return;
3854  }
3855 
3856  // case that there has not been any cell before
3857  if (this->present_cell.get() == nullptr)
3858  cell_similarity = CellSimilarity::none;
3859  else
3860  // in MappingQ, data can have been modified during the last call. Then, we
3861  // can't use that data on the new cell.
3862  if (cell_similarity == CellSimilarity::invalid_next_cell)
3863  cell_similarity = CellSimilarity::none;
3864  else
3865  cell_similarity = (cell->is_translation_of
3866  (static_cast<const typename Triangulation<dim,spacedim>::cell_iterator &>(*this->present_cell))
3867  ?
3869  :
3871 
3872  if ( (dim<spacedim) && (cell_similarity == CellSimilarity::translation) )
3873  {
3874  if (static_cast<const typename Triangulation<dim,spacedim>::cell_iterator &>
3875  (*this->present_cell)->direction_flag()
3876  != cell->direction_flag() )
3877  cell_similarity = CellSimilarity::inverted_translation;
3878  }
3879  // TODO: here, one could implement other checks for similarity, e.g. for
3880  // children of a parallelogram.
3881 }
3882 
3883 
3884 
3885 template <int dim, int spacedim>
3888 {
3889  return cell_similarity;
3890 }
3891 
3892 
3893 template <int dim, int spacedim>
3894 const unsigned int FEValuesBase<dim,spacedim>::dimension;
3895 
3896 
3897 template <int dim, int spacedim>
3899 
3900 /*------------------------------- FEValues -------------------------------*/
3901 
3902 template <int dim, int spacedim>
3904 
3905 
3906 
3907 
3908 template <int dim, int spacedim>
3910  const FiniteElement<dim,spacedim> &fe,
3911  const Quadrature<dim> &q,
3912  const UpdateFlags update_flags)
3913  :
3914  FEValuesBase<dim,spacedim> (q.size(),
3915  fe.dofs_per_cell,
3917  mapping,
3918  fe),
3919  quadrature (q)
3920 {
3921  initialize (update_flags);
3922 }
3923 
3924 
3925 
3926 template <int dim, int spacedim>
3928  const Quadrature<dim> &q,
3929  const UpdateFlags update_flags)
3930  :
3931  FEValuesBase<dim,spacedim> (q.size(),
3932  fe.dofs_per_cell,
3934  StaticMappingQ1<dim,spacedim>::mapping,
3935  fe),
3936  quadrature (q)
3937 {
3938  initialize (update_flags);
3939 }
3940 
3941 
3942 
3943 template <int dim, int spacedim>
3944 void
3946 {
3947  // You can compute normal vectors
3948  // to the cells only in the
3949  // codimension one case.
3950  if (dim != spacedim-1)
3951  Assert ((update_flags & update_normal_vectors) == false,
3952  ExcMessage ("You can only pass the 'update_normal_vectors' "
3953  "flag to FEFaceValues or FESubfaceValues objects, "
3954  "but not to an FEValues object unless the "
3955  "triangulation it refers to is embedded in a higher "
3956  "dimensional space."));
3957 
3958  const UpdateFlags flags = this->compute_update_flags (update_flags);
3959 
3960  // initialize the base classes
3961  if (flags & update_mapping)
3962  this->mapping_output.initialize(this->n_quadrature_points, flags);
3963  this->finite_element_output.initialize(this->n_quadrature_points, *this->fe, flags);
3964 
3965  // then get objects into which the FE and the Mapping can store
3966  // intermediate data used across calls to reinit. we can do this in parallel
3969  *this->fe,
3970  flags,
3971  *this->mapping,
3972  quadrature,
3973  this->finite_element_output);
3975  mapping_get_data;
3976  if (flags & update_mapping)
3978  *this->mapping,
3979  flags,
3980  quadrature);
3981 
3982  this->update_flags = flags;
3983 
3984  // then collect answers from the two task above
3985  this->fe_data.reset (fe_get_data.return_value());
3986  if (flags & update_mapping)
3987  this->mapping_data.reset (mapping_get_data.return_value());
3988  else
3989  this->mapping_data.reset (new typename Mapping<dim,spacedim>::InternalDataBase());
3990 }
3991 
3992 
3993 namespace
3994 {
3995  // Reset a unique_ptr. If we can, do not de-allocate the previously
3996  // held memory but re-use it for the next item to avoid the repeated
3997  // memory allocation. We do this because FEValues objects are heavily
3998  // used in multithreaded contexts where memory allocations are evil.
3999  template <typename Type, typename Pointer, typename Iterator>
4000  void
4001  reset_pointer_in_place_if_possible
4002  (std::unique_ptr<Pointer> &present_cell,
4003  const Iterator &new_cell)
4004  {
4005  // see if the existing pointer is non-null and if the type of
4006  // the old object pointed to matches that of the one we'd
4007  // like to create
4008  if (present_cell.get()
4009  &&
4010  (typeid(*present_cell.get()) == typeid(Type)))
4011  {
4012  // call destructor of the old object
4013  static_cast<const Type *>(present_cell.get())->~Type();
4014 
4015  // then construct a new object in-place
4016  new (const_cast<void *>(static_cast<const void *>(present_cell.get()))) Type(new_cell);
4017  }
4018  else
4019  // if the types don't match, there is nothing we can do here
4020  present_cell.reset (new Type(new_cell));
4021  }
4022 }
4023 
4024 
4025 template <int dim, int spacedim>
4027 {
4028  // no FE in this cell, so no assertion
4029  // necessary here
4030  this->maybe_invalidate_previous_present_cell (cell);
4031  this->check_cell_similarity(cell);
4032 
4033  reset_pointer_in_place_if_possible<typename FEValuesBase<dim,spacedim>::TriaCellIterator>
4034  (this->present_cell, cell);
4035 
4036  // this was the part of the work
4037  // that is dependent on the actual
4038  // data type of the iterator. now
4039  // pass on to the function doing
4040  // the real work.
4041  do_reinit ();
4042 }
4043 
4044 
4045 
4046 template <int dim, int spacedim>
4047 template <template <int, int> class DoFHandlerType, bool lda>
4048 void
4050 (const TriaIterator<DoFCellAccessor<DoFHandlerType<dim,spacedim>, lda> > &cell)
4051 {
4052  // assert that the finite elements
4053  // passed to the constructor and
4054  // used by the DoFHandler used by
4055  // this cell, are the same
4056  Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
4057  static_cast<const FiniteElementData<dim>&>(cell->get_fe()),
4059 
4060  this->maybe_invalidate_previous_present_cell (cell);
4061  this->check_cell_similarity(cell);
4062 
4063  reset_pointer_in_place_if_possible<typename FEValuesBase<dim,spacedim>::template
4065  lda> > > >
4066  (this->present_cell, cell);
4067 
4068  // this was the part of the work
4069  // that is dependent on the actual
4070  // data type of the iterator. now
4071  // pass on to the function doing
4072  // the real work.
4073  do_reinit ();
4074 }
4075 
4076 
4077 
4078 template <int dim, int spacedim>
4080 {
4081  // first call the mapping and let it generate the data
4082  // specific to the mapping. also let it inspect the
4083  // cell similarity flag and, if necessary, update
4084  // it
4085  if (this->update_flags & update_mapping)
4086  {
4087  this->cell_similarity
4088  = this->get_mapping().fill_fe_values(*this->present_cell,
4089  this->cell_similarity,
4090  quadrature,
4091  *this->mapping_data,
4092  this->mapping_output);
4093  }
4094 
4095  // then call the finite element and, with the data
4096  // already filled by the mapping, let it compute the
4097  // data for the mapped shape function values, gradients,
4098  // etc.
4099  this->get_fe().fill_fe_values(*this->present_cell,
4100  this->cell_similarity,
4101  this->quadrature,
4102  this->get_mapping(),
4103  *this->mapping_data,
4104  this->mapping_output,
4105  *this->fe_data,
4106  this->finite_element_output);
4107 }
4108 
4109 
4110 
4111 template <int dim, int spacedim>
4112 std::size_t
4114 {
4117 }
4118 
4119 
4120 /*------------------------------- FEFaceValuesBase --------------------------*/
4121 
4122 
4123 template <int dim, int spacedim>
4125  const unsigned int dofs_per_cell,
4126  const UpdateFlags,
4127  const Mapping<dim,spacedim> &mapping,
4128  const FiniteElement<dim,spacedim> &fe,
4129  const Quadrature<dim-1>& quadrature)
4130  :
4131  FEValuesBase<dim,spacedim> (n_q_points,
4132  dofs_per_cell,
4134  mapping,
4135  fe),
4136  present_face_index (numbers::invalid_unsigned_int),
4137  quadrature(quadrature)
4138 {}
4139 
4140 
4141 
4142 template <int dim, int spacedim>
4143 const std::vector<Tensor<1,spacedim> > &
4145 {
4146  Assert (this->update_flags & update_boundary_forms,
4147  (typename FEValuesBase<dim,spacedim>::ExcAccessToUninitializedField("update_boundary_forms")));
4148  return this->mapping_output.boundary_forms;
4149 }
4150 
4151 
4152 
4153 template <int dim, int spacedim>
4154 std::size_t
4156 {
4159 }
4160 
4161 
4162 /*------------------------------- FEFaceValues -------------------------------*/
4163 
4164 template <int dim, int spacedim>
4165 const unsigned int FEFaceValues<dim,spacedim>::dimension;
4166 
4167 template <int dim, int spacedim>
4169 
4170 
4171 template <int dim, int spacedim>
4173  const FiniteElement<dim,spacedim> &fe,
4174  const Quadrature<dim-1> &quadrature,
4175  const UpdateFlags update_flags)
4176  :
4177  FEFaceValuesBase<dim,spacedim> (quadrature.size(),
4178  fe.dofs_per_cell,
4179  update_flags,
4180  mapping,
4181  fe, quadrature)
4182 {
4183  initialize (update_flags);
4184 }
4185 
4186 
4187 
4188 template <int dim, int spacedim>
4190  const Quadrature<dim-1> &quadrature,
4191  const UpdateFlags update_flags)
4192  :
4193  FEFaceValuesBase<dim,spacedim> (quadrature.size(),
4194  fe.dofs_per_cell,
4195  update_flags,
4196  StaticMappingQ1<dim,spacedim>::mapping,
4197  fe, quadrature)
4198 {
4199  initialize (update_flags);
4200 }
4201 
4202 
4203 
4204 template <int dim, int spacedim>
4205 void
4207 {
4208  const UpdateFlags flags = this->compute_update_flags (update_flags);
4209 
4210  // initialize the base classes
4211  if (flags & update_mapping)
4212  this->mapping_output.initialize(this->n_quadrature_points, flags);
4213  this->finite_element_output.initialize(this->n_quadrature_points, *this->fe, flags);
4214 
4215  // then get objects into which the FE and the Mapping can store
4216  // intermediate data used across calls to reinit. this can be done in parallel
4219  *this->fe,
4220  flags,
4221  *this->mapping,
4222  this->quadrature,
4223  this->finite_element_output);
4225  mapping_get_data;
4226  if (flags & update_mapping)
4228  *this->mapping,
4229  flags,
4230  this->quadrature);
4231 
4232  this->update_flags = flags;
4233 
4234  // then collect answers from the two task above
4235  this->fe_data.reset (fe_get_data.return_value());
4236  if (flags & update_mapping)
4237  this->mapping_data.reset (mapping_get_data.return_value());
4238  else
4239  this->mapping_data.reset (new typename Mapping<dim,spacedim>::InternalDataBase());
4240 }
4241 
4242 
4243 
4244 template <int dim, int spacedim>
4245 template <template <int, int> class DoFHandlerType, bool lda>
4246 void
4248 (const TriaIterator<DoFCellAccessor<DoFHandlerType<dim,spacedim>, lda> > &cell,
4249  const unsigned int face_no)
4250 {
4251  // assert that the finite elements
4252  // passed to the constructor and
4253  // used by the DoFHandler used by
4254  // this cell, are the same
4255  Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
4256  static_cast<const FiniteElementData<dim>&>(
4257  cell->get_dof_handler().get_fe(cell->active_fe_index ())),
4259 
4262 
4263  this->maybe_invalidate_previous_present_cell (cell);
4264  reset_pointer_in_place_if_possible<typename FEValuesBase<dim,spacedim>::template
4266  lda> > > >
4267  (this->present_cell, cell);
4268 
4269  // this was the part of the work
4270  // that is dependent on the actual
4271  // data type of the iterator. now
4272  // pass on to the function doing
4273  // the real work.
4274  do_reinit (face_no);
4275 }
4276 
4277 
4278 
4279 template <int dim, int spacedim>
4281  const unsigned int face_no)
4282 {
4285 
4286  this->maybe_invalidate_previous_present_cell (cell);
4287  reset_pointer_in_place_if_possible<typename FEValuesBase<dim,spacedim>::TriaCellIterator>
4288  (this->present_cell, cell);
4289 
4290  // this was the part of the work
4291  // that is dependent on the actual
4292  // data type of the iterator. now
4293  // pass on to the function doing
4294  // the real work.
4295  do_reinit (face_no);
4296 }
4297 
4298 
4299 
4300 template <int dim, int spacedim>
4301 void FEFaceValues<dim,spacedim>::do_reinit (const unsigned int face_no)
4302 {
4303  // first of all, set the present_face_index (if available)
4304  const typename Triangulation<dim,spacedim>::cell_iterator cell=*this->present_cell;
4305  this->present_face_index=cell->face_index(face_no);
4306 
4307  if (this->update_flags & update_mapping)
4308  {
4309  this->get_mapping().fill_fe_face_values(*this->present_cell,
4310  face_no,
4311  this->quadrature,
4312  *this->mapping_data,
4313  this->mapping_output);
4314  }
4315 
4316  this->get_fe().fill_fe_face_values(*this->present_cell,
4317  face_no,
4318  this->quadrature,
4319  this->get_mapping(),
4320  *this->mapping_data,
4321  this->mapping_output,
4322  *this->fe_data,
4323  this->finite_element_output);
4324 }
4325 
4326 
4327 /*------------------------------- FESubFaceValues -------------------------------*/
4328 
4329 
4330 template <int dim, int spacedim>
4331 const unsigned int FESubfaceValues<dim,spacedim>::dimension;
4332 
4333 template <int dim, int spacedim>
4335 
4336 
4337 
4338 template <int dim, int spacedim>
4340  const FiniteElement<dim,spacedim> &fe,
4341  const Quadrature<dim-1> &quadrature,
4342  const UpdateFlags update_flags)
4343  :
4344  FEFaceValuesBase<dim,spacedim> (quadrature.size(),
4345  fe.dofs_per_cell,
4346  update_flags,
4347  mapping,
4348  fe, quadrature)
4349 {
4350  initialize (update_flags);
4351 }
4352 
4353 
4354 
4355 template <int dim, int spacedim>
4357  const Quadrature<dim-1> &quadrature,
4358  const UpdateFlags update_flags)
4359  :
4360  FEFaceValuesBase<dim,spacedim> (quadrature.size(),
4361  fe.dofs_per_cell,
4362  update_flags,
4363  StaticMappingQ1<dim,spacedim>::mapping,
4364  fe, quadrature)
4365 {
4366  initialize (update_flags);
4367 }
4368 
4369 
4370 
4371 template <int dim, int spacedim>
4372 void
4374 {
4375  const UpdateFlags flags = this->compute_update_flags (update_flags);
4376 
4377  // initialize the base classes
4378  if (flags & update_mapping)
4379  this->mapping_output.initialize(this->n_quadrature_points, flags);
4380  this->finite_element_output.initialize(this->n_quadrature_points, *this->fe, flags);
4381 
4382  // then get objects into which the FE and the Mapping can store
4383  // intermediate data used across calls to reinit. this can be done
4384  // in parallel
4387  *this->fe,
4388  flags,
4389  *this->mapping,
4390  this->quadrature,
4391  this->finite_element_output);
4393  mapping_get_data;
4394  if (flags & update_mapping)
4396  *this->mapping,
4397  flags,
4398  this->quadrature);
4399 
4400  this->update_flags = flags;
4401 
4402  // then collect answers from the two task above
4403  this->fe_data.reset (fe_get_data.return_value());
4404  if (flags & update_mapping)
4405  this->mapping_data.reset (mapping_get_data.return_value());
4406  else
4407  this->mapping_data.reset (new typename Mapping<dim,spacedim>::InternalDataBase());
4408 }
4409 
4410 
4411 template <int dim, int spacedim>
4412 template <template <int, int> class DoFHandlerType, bool lda>
4414 (const TriaIterator<DoFCellAccessor<DoFHandlerType<dim,spacedim>, lda> > &cell,
4415  const unsigned int face_no,
4416  const unsigned int subface_no)
4417 {
4418  // assert that the finite elements
4419  // passed to the constructor and
4420  // used by the hp::DoFHandler used by
4421  // this cell, are the same
4422  Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
4423  static_cast<const FiniteElementData<dim>&>(
4424  cell->get_dof_handler().get_fe(cell->active_fe_index ())),
4428  // We would like to check for
4429  // subface_no < cell->face(face_no)->n_children(),
4430  // but unfortunately the current
4431  // function is also called for
4432  // faces without children (see
4433  // tests/fe/mapping.cc). Therefore,
4434  // we must use following workaround
4435  // of two separate assertions
4436  Assert (cell->face(face_no)->has_children() ||
4437  subface_no < GeometryInfo<dim>::max_children_per_face,
4439  Assert (!cell->face(face_no)->has_children() ||
4440  subface_no < cell->face(face_no)->number_of_children(),
4441  ExcIndexRange (subface_no, 0, cell->face(face_no)->number_of_children()));
4442  Assert (cell->has_children() == false,
4443  ExcMessage ("You can't use subface data for cells that are "
4444  "already refined. Iterate over their children "
4445  "instead in these cases."));
4446 
4447  this->maybe_invalidate_previous_present_cell (cell);
4448  reset_pointer_in_place_if_possible<typename FEValuesBase<dim,spacedim>::template
4450  lda> > > >
4451  (this->present_cell, cell);
4452 
4453  // this was the part of the work
4454  // that is dependent on the actual
4455  // data type of the iterator. now
4456  // pass on to the function doing
4457  // the real work.
4458  do_reinit (face_no, subface_no);
4459 }
4460 
4461 
4462 template <int dim, int spacedim>
4464  const unsigned int face_no,
4465  const unsigned int subface_no)
4466 {
4469  Assert (subface_no < cell->face(face_no)->n_children(),
4470  ExcIndexRange (subface_no, 0, cell->face(face_no)->n_children()));
4471 
4472  this->maybe_invalidate_previous_present_cell (cell);
4473  reset_pointer_in_place_if_possible<typename FEValuesBase<dim,spacedim>::TriaCellIterator>
4474  (this->present_cell, cell);
4475 
4476  // this was the part of the work
4477  // that is dependent on the actual
4478  // data type of the iterator. now
4479  // pass on to the function doing
4480  // the real work.
4481  do_reinit (face_no, subface_no);
4482 }
4483 
4484 
4485 
4486 template <int dim, int spacedim>
4487 void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
4488  const unsigned int subface_no)
4489 {
4490  // first of all, set the present_face_index
4491  // (if available)
4492  const typename Triangulation<dim,spacedim>::cell_iterator cell=*this->present_cell;
4493 
4494  if (!cell->face(face_no)->has_children())
4495  // no subfaces at all, so set
4496  // present_face_index to this face rather
4497  // than any subface
4498  this->present_face_index=cell->face_index(face_no);
4499  else if (dim!=3)
4500  this->present_face_index=cell->face(face_no)->child_index(subface_no);
4501  else
4502  {
4503  // this is the same logic we use in
4504  // cell->neighbor_child_on_subface(). See
4505  // there for an explanation of the
4506  // different cases
4507  unsigned int subface_index=numbers::invalid_unsigned_int;
4508  switch (cell->subface_case(face_no))
4509  {
4513  subface_index=cell->face(face_no)->child_index(subface_no);
4514  break;
4517  subface_index=cell->face(face_no)->child(subface_no/2)->child_index(subface_no%2);
4518  break;
4521  switch (subface_no)
4522  {
4523  case 0:
4524  case 1:
4525  subface_index=cell->face(face_no)->child(0)->child_index(subface_no);
4526  break;
4527  case 2:
4528  subface_index=cell->face(face_no)->child_index(1);
4529  break;
4530  default:
4531  Assert(false, ExcInternalError());
4532  }
4533  break;
4536  switch (subface_no)
4537  {
4538  case 0:
4539  subface_index=cell->face(face_no)->child_index(0);
4540  break;
4541  case 1:
4542  case 2:
4543  subface_index=cell->face(face_no)->child(1)->child_index(subface_no-1);
4544  break;
4545  default:
4546  Assert(false, ExcInternalError());
4547  }
4548  break;
4549  default:
4550  Assert(false, ExcInternalError());
4551  break;
4552  }
4553  Assert(subface_index!=numbers::invalid_unsigned_int,
4554  ExcInternalError());
4555  this->present_face_index=subface_index;
4556  }
4557 
4558  // now ask the mapping and the finite element to do the actual work
4559  if (this->update_flags & update_mapping)
4560  {
4561  this->get_mapping().fill_fe_subface_values(*this->present_cell,
4562  face_no,
4563  subface_no,
4564  this->quadrature,
4565  *this->mapping_data,
4566  this->mapping_output);
4567  }
4568 
4569  this->get_fe().fill_fe_subface_values(*this->present_cell,
4570  face_no,
4571  subface_no,
4572  this->quadrature,
4573  this->get_mapping(),
4574  *this->mapping_data,
4575  this->mapping_output,
4576  *this->fe_data,
4577  this->finite_element_output);
4578 }
4579 
4580 
4581 /*------------------------------- Explicit Instantiations -------------*/
4582 #define SPLIT_INSTANTIATIONS_COUNT 2
4583 #ifndef SPLIT_INSTANTIATIONS_INDEX
4584 #define SPLIT_INSTANTIATIONS_INDEX 0
4585 #endif
4586 #include "fe_values.inst"
4587 
4588 DEAL_II_NAMESPACE_CLOSE
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:1801
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3141
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:1733
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:2047
ProductType< Number, typename Scalar< dim, spacedim >::third_derivative_type >::type third_derivative_type
Definition: fe_values.h:205
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:1466
DEAL_II_CUDA_HOST_DEV Tensor()
Definition: tensor.h:825
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:1336
static const unsigned int invalid_unsigned_int
Definition: types.h:173
ProductType< Number, typename Vector< dim, spacedim >::third_derivative_type >::type third_derivative_type
Definition: fe_values.h:640
std::size_t size() const
Definition: array_view.h:382
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:2003
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:3599
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:1938
Cache(const FEValuesBase< dim, spacedim > &fe_values)
Definition: fe_values.cc:2093
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1183
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:1958
std::size_t memory_consumption() const
Definition: fe_values.cc:4155
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
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:1422
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:1534
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:1667
Task< RT > new_task(const std::function< RT()> &function)
iterator end() const
Definition: array_view.h:400
static const unsigned int n_independent_components
const unsigned int component
Definition: fe_values.h:486
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:2679
signed int value_type
Definition: index_set.h:98
Volume element.
Outer normal vector, not normalized.
static::ExceptionBase & ExcFEDontMatch()
Scalar & operator=(const Scalar< dim, spacedim > &)
Definition: fe_values.cc:156
TriaCellIterator(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Definition: fe_values.cc:2485
iterator begin() const
Definition: array_view.h:390
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:1757
void do_reinit(const unsigned int face_no)
Definition: fe_values.cc:4301
const Triangulation< dim, spacedim >::cell_iterator get_cell() const
Definition: fe_values.cc:3700
ProductType< Number, typename SymmetricTensor< 2, dim, spacedim >::divergence_type >::type divergence_type
Definition: fe_values.h:1166
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:1578
bool is_primitive() const
Definition: fe.h:3199
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type > > &hessians) const
Definition: fe_values.cc:3377
void check_cell_similarity(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Definition: fe_values.cc:3833
::internal::FEValuesViews::Cache< dim, spacedim > fe_values_views_cache
Definition: fe_values.h:3057
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:1490
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:1510
Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
ProductType< Number, typename Vector< dim, spacedim >::hessian_type >::type hessian_type
Definition: fe_values.h:634
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:1642
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:1598
void do_reinit(const unsigned int face_no, const unsigned int subface_no)
Definition: fe_values.cc:4487
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
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:4414
ProductType< Number, typename Vector< dim, spacedim >::gradient_type >::type gradient_type
Definition: fe_values.h:604
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:1849
UpdateFlags compute_update_flags(const UpdateFlags update_flags) const
Definition: fe_values.cc:3752
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:1622
ProductType< Number, typename Vector< dim, spacedim >::value_type >::type value_type
Definition: fe_values.h:598
ProductType< Number, typename Scalar< dim, spacedim >::gradient_type >::type gradient_type
Definition: fe_values.h:187
void get_function_laplacians(const InputVector &fe_function, std::vector< typename InputVector::value_type > &laplacians) const
Definition: fe_values.cc:3478
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
Definition: fe_values.h:2995
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:1777
const std::vector< Tensor< 1, spacedim > > & get_boundary_forms() const
Definition: fe_values.cc:4144
ProductType< Number, typename Scalar< dim, spacedim >::hessian_type >::type hessian_type
Definition: fe_values.h:199
virtual types::global_dof_index n_dofs_for_dof_handler() const
Definition: fe_values.cc:2503
No update.
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition: tria.h:1487
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:1983
unsigned int global_dof_index
Definition: types.h:88
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:1378
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:1446
#define Assert(cond, exc)
Definition: exceptions.h:337
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:1870
unsigned int n_nonzero_components(const unsigned int i) const
Definition: fe.h:3188
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access > > &cell, const unsigned int face_no)
Definition: fe_values.cc:4248
Abstract base class for mapping classes.
Definition: dof_tools.h:46
std::vector< ShapeFunctionData > shape_function_data
Definition: fe_values.h:1095
std::size_t memory_consumption() const
Definition: fe_values.cc:3732
const Triangulation< dim, spacedim >::cell_iterator cell
Definition: fe_values.cc:2400
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:1402
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:1914
void invalidate_present_cell()
Definition: fe_values.cc:3770
ProductType< Number, typename Scalar< dim, spacedim >::value_type >::type value_type
Definition: fe_values.h:181
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition: fe_values.h:2975
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:3177
Vector & operator=(const Vector< dim, spacedim > &)
Definition: fe_values.cc:251
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Definition: fe_values.cc:3720
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:1713
static const char *const message_string
Definition: fe_values.cc:2411
unsigned int n_components() const
virtual void get_interpolated_dof_values(const IndexSet &in, Vector< IndexSet::value_type > &out) const
Definition: fe_values.cc:2516
ProductType< Number, typename Vector< dim, spacedim >::divergence_type >::type divergence_type
Definition: fe_values.h:616
Second derivatives of shape functions.
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:1139
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:1554
ProductType< Number, typename Vector< dim, spacedim >::symmetric_gradient_type >::type symmetric_gradient_type
Definition: fe_values.h:610
const unsigned int dofs_per_cell
Definition: fe_base.h:295
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:2027
ProductType< Number, typename Vector< dim, spacedim >::value_type >::type laplacian_type
Definition: fe_values.h:622
void get_function_divergences_from_local_dof_values(const InputVector &dof_values, std::vector< typename OutputType< typename InputVector::value_type >::divergence_type > &values) const
Definition: fe_values.cc:2072
ArrayView< ElementType > make_array_view(std::vector< ElementType > &vector)
Definition: array_view.h:455
Definition: mpi.h:53
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:4124
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3003
ProductType< Number, typename SymmetricTensor< 2, dim, spacedim >::value_type >::type value_type
Definition: fe_values.h:1160
void initialize(const UpdateFlags update_flags)
Definition: fe_values.cc:3945
Shape function gradients.
Normal vectors.
void maybe_invalidate_previous_present_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell)
Definition: fe_values.cc:3788
T signaling_nan()
value_type * data() const noexcept
Definition: array_view.h:358
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:1688
Definition: fe.h:33
void initialize(const UpdateFlags update_flags)
Definition: fe_values.cc:4373
const std::vector< Tensor< 1, spacedim > > & get_all_normal_vectors() const
Definition: fe_values.cc:3709
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:1821
static::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1623
ProductType< Number, typename Vector< dim, spacedim >::curl_type >::type curl_type
Definition: fe_values.h:628
const FiniteElement< dim, spacedim > & get_fe() const
static unsigned int n_threads()
unsigned int size(const unsigned int i) const
void initialize(const UpdateFlags update_flags)
Definition: fe_values.cc:4206
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:3909
void do_reinit()
Definition: fe_values.cc:4079
Point< 3 > point(const gp_Pnt &p)
Definition: utilities.cc:174
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access > > &cell)
Definition: fe_values.cc:4050
ProductType< Number, typename Tensor< 2, dim, spacedim >::divergence_type >::type divergence_type
Definition: fe_values.h:1422
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:4339
ProductType< Number, typename Scalar< dim, spacedim >::value_type >::type laplacian_type
Definition: fe_values.h:193
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:1358
FEFaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim-1 > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:4172
CellSimilarity::Similarity get_cell_similarity() const
Definition: fe_values.cc:3887
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:1894
std::vector< ShapeFunctionData > shape_function_data
Definition: fe_values.h:491
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
UpdateFlags update_flags
Definition: fe_values.h:3014
std::size_t memory_consumption() const
Definition: fe_values.cc:4113
static::ExceptionBase & ExcInternalError()
ProductType< Number, typename Tensor< 2, dim, spacedim >::value_type >::type value_type
Definition: fe_values.h:1416
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
Definition: fe_values.cc:3273