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