Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
fe_poly.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: fe_poly.templates.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2006 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #include <deal.II/base/qprojector.h>
19 #include <deal.II/base/polynomial_space.h>
20 #include <deal.II/base/tensor_product_polynomials.h>
21 #include <deal.II/fe/fe_values.h>
22 #include <deal.II/fe/fe_poly.h>
23 
24 
26 
27 template <class POLY, int dim, int spacedim>
28 FE_Poly<POLY,dim,spacedim>::FE_Poly (const POLY &poly_space,
29  const FiniteElementData<dim> &fe_data,
30  const std::vector<bool> &restriction_is_additive_flags,
31  const std::vector<ComponentMask> &nonzero_components):
32  FiniteElement<dim,spacedim> (fe_data,
33  restriction_is_additive_flags,
34  nonzero_components),
35  poly_space(poly_space)
36 {
37  AssertDimension(dim, POLY::dimension);
38 }
39 
40 
41 template <class POLY, int dim, int spacedim>
42 unsigned int
44 {
45  return this->degree;
46 }
47 
48 
49 template <class POLY, int dim, int spacedim>
50 double
52  const Point<dim> &p) const
53 {
54  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
55  return poly_space.compute_value(i, p);
56 }
57 
58 
59 template <class POLY, int dim, int spacedim>
60 double
62  const Point<dim> &p,
63  const unsigned int component) const
64 {
65  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
66  Assert (component == 0, ExcIndexRange (component, 0, 1));
67  return poly_space.compute_value(i, p);
68 }
69 
70 
71 
72 template <class POLY, int dim, int spacedim>
75  const Point<dim> &p) const
76 {
77  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
78  return poly_space.compute_grad(i, p);
79 }
80 
81 
82 
83 template <class POLY, int dim, int spacedim>
86  const Point<dim> &p,
87  const unsigned int component) const
88 {
89  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
90  Assert (component == 0, ExcIndexRange (component, 0, 1));
91  return poly_space.compute_grad(i, p);
92 }
93 
94 
95 
96 template <class POLY, int dim, int spacedim>
99  const Point<dim> &p) const
100 {
101  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
102  return poly_space.compute_grad_grad(i, p);
103 }
104 
105 
106 
107 template <class POLY, int dim, int spacedim>
110  const Point<dim> &p,
111  const unsigned int component) const
112 {
113  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
114  Assert (component == 0, ExcIndexRange (component, 0, 1));
115  return poly_space.compute_grad_grad(i, p);
116 }
117 
118 
119 
120 //---------------------------------------------------------------------------
121 // Auxiliary functions
122 //---------------------------------------------------------------------------
123 
124 
125 
126 
127 template <class POLY, int dim, int spacedim>
130 {
131  // for this kind of elements, only
132  // the values can be precomputed
133  // once and for all. set this flag
134  // if the values are requested at
135  // all
136  return (update_default | (flags & update_values));
137 }
138 
139 
140 
141 template <class POLY, int dim, int spacedim>
144 {
146 
147  if (flags & update_gradients)
148  out |= update_gradients | update_covariant_transformation;
149  if (flags & update_hessians)
150  out |= update_hessians | update_covariant_transformation;
151  if (flags & update_cell_normal_vectors)
152  out |= update_cell_normal_vectors | update_JxW_values;
153 
154  return out;
155 }
156 
157 
158 
159 //---------------------------------------------------------------------------
160 // Data field initialization
161 //---------------------------------------------------------------------------
162 
163 template <class POLY, int dim, int spacedim>
166  const Mapping<dim,spacedim> &mapping,
167  const Quadrature<dim> &quadrature) const
168 {
169  // generate a new data object and
170  // initialize some fields
171  InternalData *data = new InternalData;
172 
173  // check what needs to be
174  // initialized only once and what
175  // on every cell/face/subface we
176  // visit
177  data->update_once = update_once(update_flags);
178  data->update_each = update_each(update_flags);
179  data->update_flags = data->update_once | data->update_each;
180 
181  const UpdateFlags flags(data->update_flags);
182  const unsigned int n_q_points = quadrature.size();
183 
184  // some scratch arrays
185  std::vector<double> values(0);
186  std::vector<Tensor<1,dim> > grads(0);
187  std::vector<Tensor<2,dim> > grad_grads(0);
188 
189  // initialize fields only if really
190  // necessary. otherwise, don't
191  // allocate memory
192  if (flags & update_values)
193  {
194  values.resize (this->dofs_per_cell);
195  data->shape_values.resize (this->dofs_per_cell,
196  std::vector<double> (n_q_points));
197  }
198 
199  if (flags & update_gradients)
200  {
201  grads.resize (this->dofs_per_cell);
202  data->shape_gradients.resize (this->dofs_per_cell,
203  std::vector<Tensor<1,dim> > (n_q_points));
204  }
205 
206  // if second derivatives through
207  // finite differencing is required,
208  // then initialize some objects for
209  // that
210  if (flags & update_hessians)
211  data->initialize_2nd (this, mapping, quadrature);
212 
213  // next already fill those fields
214  // of which we have information by
215  // now. note that the shape
216  // gradients are only those on the
217  // unit cell, and need to be
218  // transformed when visiting an
219  // actual cell
220  if (flags & (update_values | update_gradients))
221  for (unsigned int i=0; i<n_q_points; ++i)
222  {
223  poly_space.compute(quadrature.point(i),
224  values, grads, grad_grads);
225 
226  if (flags & update_values)
227  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
228  data->shape_values[k][i] = values[k];
229 
230  if (flags & update_gradients)
231  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
232  data->shape_gradients[k][i] = grads[k];
233  }
234  return data;
235 }
236 
237 
238 
239 
240 //---------------------------------------------------------------------------
241 // Fill data of FEValues
242 //---------------------------------------------------------------------------
243 
244 
245 template <class POLY, int dim, int spacedim>
246 void
248 (const Mapping<dim,spacedim> &mapping,
249  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
250  const Quadrature<dim> &quadrature,
251  typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
254  CellSimilarity::Similarity &cell_similarity) const
255 {
256  // convert data object to internal
257  // data for this class. fails with
258  // an exception if that is not
259  // possible
260  Assert (dynamic_cast<InternalData *> (&fedata) != 0, ExcInternalError());
261  InternalData &fe_data = static_cast<InternalData &> (fedata);
262 
263  const UpdateFlags flags(fe_data.current_update_flags());
264 
265  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
266  {
267  if (flags & update_values)
268  for (unsigned int i=0; i<quadrature.size(); ++i)
269  data.shape_values(k,i) = fe_data.shape_values[k][i];
270 
271  if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
272  mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
273  mapping_data, mapping_covariant);
274  }
275 
276  if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
277  this->compute_2nd (mapping, cell, QProjector<dim>::DataSetDescriptor::cell(),
278  mapping_data, fe_data, data);
279 }
280 
281 
282 
283 template <class POLY, int dim, int spacedim>
284 void
287  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
288  const unsigned int face,
289  const Quadrature<dim-1> &quadrature,
290  typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
292  FEValuesData<dim,spacedim> &data) const
293 {
294  // convert data object to internal
295  // data for this class. fails with
296  // an exception if that is not
297  // possible
298  Assert (dynamic_cast<InternalData *> (&fedata) != 0, ExcInternalError());
299  InternalData &fe_data = static_cast<InternalData &> (fedata);
300 
301  // offset determines which data set
302  // to take (all data sets for all
303  // faces are stored contiguously)
304 
305  const typename QProjector<dim>::DataSetDescriptor offset
307  cell->face_orientation(face),
308  cell->face_flip(face),
309  cell->face_rotation(face),
310  quadrature.size());
311 
312  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
313 
314  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
315  {
316  if (flags & update_values)
317  for (unsigned int i=0; i<quadrature.size(); ++i)
318  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
319 
320  if (flags & update_gradients)
321  mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
322  data.shape_gradients[k],
323  mapping_data, mapping_covariant);
324  }
325 
326  if (flags & update_hessians)
327  this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
328 }
329 
330 
331 //codimension 1
332 // template <>
333 // inline
334 // void
335 // FE_Poly<TensorProductPolynomials<1>,1,2>::fill_fe_subface_values (const Mapping<1,2> &,
336 // const Triangulation<1,2>::cell_iterator &,
337 // const unsigned int,
338 // const unsigned int,
339 // const Quadrature<0> &,
340 // Mapping<1,2>::InternalDataBase &,
341 // Mapping<1,2>::InternalDataBase &,
342 // FEValuesData<1,2> &) const
343 // {
344 // AssertThrow(false, ExcNotImplemented());
345 // }
346 
347 
348 // template <>
349 // inline
350 // void
351 // FE_Poly<TensorProductPolynomials<2>,2,3>::fill_fe_subface_values (const Mapping<2,3> &,
352 // const Triangulation<2,3>::cell_iterator &,
353 // const unsigned int,
354 // const unsigned int,
355 // const Quadrature<1> &,
356 // Mapping<2,3>::InternalDataBase &,
357 // Mapping<2,3>::InternalDataBase &,
358 // FEValuesData<2,3> &) const
359 // {
360 // AssertThrow(false, ExcNotImplemented());
361 // }
362 
363 
364 // template <>
365 // inline
366 // void
367 // FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_subface_values (const Mapping<1,2> &,
368 // const Triangulation<1,2>::cell_iterator &,
369 // const unsigned int,
370 // const unsigned int,
371 // const Quadrature<0> &,
372 // Mapping<1,2>::InternalDataBase &,
373 // Mapping<1,2>::InternalDataBase &,
374 // FEValuesData<1,2> &) const
375 // {
376 // AssertThrow(false, ExcNotImplemented());
377 // }
378 
379 
380 // template <>
381 // inline
382 // void
383 // FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_subface_values (const Mapping<2,3> &,
384 // const Triangulation<2,3>::cell_iterator &,
385 // const unsigned int,
386 // const unsigned int,
387 // const Quadrature<1> &,
388 // Mapping<2,3>::InternalDataBase &,
389 // Mapping<2,3>::InternalDataBase &,
390 // FEValuesData<2,3> &) const
391 // {
392 // AssertThrow(false, ExcNotImplemented());
393 // }
394 
395 
396 
397 
398 template <class POLY, int dim, int spacedim>
399 void
401  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
402  const unsigned int face,
403  const unsigned int subface,
404  const Quadrature<dim-1> &quadrature,
405  typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
407  FEValuesData<dim,spacedim> &data) const
408 {
409  // convert data object to internal
410  // data for this class. fails with
411  // an exception if that is not
412  // possible
413  Assert (dynamic_cast<InternalData *> (&fedata) != 0, ExcInternalError());
414  InternalData &fe_data = static_cast<InternalData &> (fedata);
415 
416  // offset determines which data set
417  // to take (all data sets for all
418  // sub-faces are stored contiguously)
419 
420  const typename QProjector<dim>::DataSetDescriptor offset
422  cell->face_orientation(face),
423  cell->face_flip(face),
424  cell->face_rotation(face),
425  quadrature.size(),
426  cell->subface_case(face));
427 
428  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
429 
430  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
431  {
432  if (flags & update_values)
433  for (unsigned int i=0; i<quadrature.size(); ++i)
434  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
435 
436  if (flags & update_gradients)
437  mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
438  data.shape_gradients[k],
439  mapping_data, mapping_covariant);
440  }
441 
442  if (flags & update_hessians)
443  this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
444 }
445 
446 
447 
448 namespace internal
449 {
450  template <class POLY>
451  inline
452  std::vector<unsigned int>
453  get_poly_space_numbering (const POLY &)
454  {
455  Assert (false, ExcNotImplemented());
456  return std::vector<unsigned int>();
457  }
458 
459  template <class POLY>
460  inline
461  std::vector<unsigned int>
462  get_poly_space_numbering_inverse (const POLY &)
463  {
464  Assert (false, ExcNotImplemented());
465  return std::vector<unsigned int>();
466  }
467 
468  template <int dim, typename POLY>
469  inline
470  std::vector<unsigned int>
471  get_poly_space_numbering (const TensorProductPolynomials<dim,POLY> &poly)
472  {
473  return poly.get_numbering();
474  }
475 
476  template <int dim, typename POLY>
477  inline
478  std::vector<unsigned int>
479  get_poly_space_numbering_inverse (const TensorProductPolynomials<dim,POLY> &poly)
480  {
481  return poly.get_numbering_inverse();
482  }
483 }
484 
485 
486 
487 template <class POLY, int dim, int spacedim>
488 std::vector<unsigned int>
490 {
491  return internal::get_poly_space_numbering (poly_space);
492 }
493 
494 
495 
496 
497 template <class POLY, int dim, int spacedim>
498 std::vector<unsigned int>
500 {
501  return internal::get_poly_space_numbering_inverse (poly_space);
502 }
503 
504 
505 
506 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
Shape function values.
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void initialize_2nd(const FiniteElement< dim, spacedim > *element, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature)
virtual void fill_fe_values(const Mapping< dim, spacedim > &mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Quadrature< dim > &quadrature, typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, typename Mapping< dim, spacedim >::InternalDataBase &fe_internal, FEValuesData< dim, spacedim > &data, CellSimilarity::Similarity &cell_similarity) const
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
virtual Mapping< dim, spacedim >::InternalDataBase * get_data(const UpdateFlags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature) const
virtual void fill_fe_subface_values(const Mapping< dim, spacedim > &mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim-1 > &quadrature, typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, typename Mapping< dim, spacedim >::InternalDataBase &fe_internal, FEValuesData< dim, spacedim > &data) const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
std::vector< unsigned int > get_poly_space_numbering() const
std::vector< std::vector< Tensor< 1, dim > > > shape_gradients
Definition: fe_poly.h:278
unsigned int get_degree() const
virtual void transform(const VectorSlice< const std::vector< Tensor< 1, dim > > > input, VectorSlice< std::vector< Tensor< 1, spacedim > > > output, const InternalDataBase &internal, const MappingType type) const =0
No update.
#define Assert(cond, exc)
Definition: exceptions.h:299
UpdateFlags
UpdateFlags current_update_flags() const
virtual UpdateFlags update_once(const UpdateFlags flags) const
UpdateFlags update_flags
Definition: mapping.h:353
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const
ShapeVector shape_values
Definition: fe_values.h:1171
Second derivatives of shape functions.
std::vector< std::vector< double > > shape_values
Definition: fe_poly.h:267
GradientVector shape_gradients
Definition: fe_values.h:1178
unsigned int size() const
const std::vector< unsigned int > & get_numbering_inverse() const
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
const std::vector< unsigned int > & get_numbering() const
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
virtual void fill_fe_face_values(const Mapping< dim, spacedim > &mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim-1 > &quadrature, typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, typename Mapping< dim, spacedim >::InternalDataBase &fe_internal, FEValuesData< dim, spacedim > &data) const
const VectorSlice< const VECTOR > make_slice(VECTOR &v)
Definition: vector_slice.h:143
std::vector< unsigned int > get_poly_space_numbering_inverse() const
Shape function gradients.
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Covariant mapping (see Mapping::transform() for details)
Definition: mapping.h:63
FE_Poly(const POLY &poly_space, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
::ExceptionBase & ExcNotImplemented()
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
UpdateFlags update_once
Definition: mapping.h:359
virtual UpdateFlags update_each(const UpdateFlags flags) const
::ExceptionBase & ExcInternalError()
const Point< dim > & point(const unsigned int i) const
Covariant transformation.
UpdateFlags update_each
Definition: mapping.h:365