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