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