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