Reference documentation for deal.II version 8.1.0
shape_info.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: shape_info.templates.h 31495 2013-10-31 12:58:57Z kronbichler @f$
3 //
4 // Copyright (C) 2011 - 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/utilities.h>
19 #include <deal.II/base/memory_consumption.h>
20 #include <deal.II/base/polynomial.h>
21 #include <deal.II/base/tensor_product_polynomials.h>
22 #include <deal.II/base/polynomials_piecewise.h>
23 #include <deal.II/fe/fe_poly.h>
24 
25 #include <deal.II/matrix_free/shape_info.h>
26 
27 
29 
30 
31 namespace internal
32 {
33  namespace MatrixFreeFunctions
34  {
35 
36  // ----------------- actual ShapeInfo functions --------------------
37 
38  template <typename Number>
40  :
41  n_q_points (0),
42  dofs_per_cell (0)
43  {}
44 
45 
46 
47  template <typename Number>
48  template <int dim>
49  void
51  const FiniteElement<dim> &fe)
52  {
53  Assert (fe.n_components() == 1,
54  ExcMessage("FEEvaluation only works for scalar finite elements."));
55 
56  const unsigned int n_dofs_1d = fe.degree+1,
57  n_q_points_1d = quad.size();
58  AssertDimension(fe.dofs_per_cell, Utilities::fixed_power<dim>(n_dofs_1d));
59  std::vector<unsigned int> lexicographic (fe.dofs_per_cell);
60 
61  // renumber (this is necessary for FE_Q, for example, since there the
62  // vertex DoFs come first, which is incompatible with the lexicographic
63  // ordering necessary to apply tensor products efficiently)
64  {
65  const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
66  dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(&fe);
68  PiecewisePolynomial<double> >,dim,dim> *fe_poly_piece =
69  dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
71  Assert (fe_poly != 0 || fe_poly_piece, ExcNotImplemented());
72  lexicographic = fe_poly != 0 ?
74  fe_poly_piece->get_poly_space_numbering_inverse();
75 
76  // to evaluate 1D polynomials, evaluate along the line where y=z=0,
77  // assuming that shape_value(0,Point<dim>()) == 1. otherwise, need
78  // other entry point (e.g. generating a 1D element by reading the
79  // name, as done before r29356)
80  Assert(std::fabs(fe.shape_value(lexicographic[0], Point<dim>())-1) < 1e-13,
82  }
83 
84  n_q_points = Utilities::fixed_power<dim>(n_q_points_1d);
85  dofs_per_cell = Utilities::fixed_power<dim>(n_dofs_1d);
86  n_q_points_face = dim>1?Utilities::fixed_power<dim-1>(n_q_points_1d):1;
87  dofs_per_face = dim>1?Utilities::fixed_power<dim-1>(n_dofs_1d):1;
88 
89  const unsigned int array_size = n_dofs_1d*n_q_points_1d;
90  this->shape_gradients.resize_fast (array_size);
91  this->shape_values.resize_fast (array_size);
92  this->shape_hessians.resize_fast (array_size);
93 
94  this->face_value[0].resize(n_dofs_1d);
95  this->face_gradient[0].resize(n_dofs_1d);
96  this->subface_value[0].resize(array_size);
97  this->face_value[1].resize(n_dofs_1d);
98  this->face_gradient[1].resize(n_dofs_1d);
99  this->subface_value[1].resize(array_size);
100  this->shape_values_number.resize (array_size);
101  this->shape_gradient_number.resize (array_size);
102 
103  for (unsigned int i=0; i<n_dofs_1d; ++i)
104  {
105  // need to reorder from hierarchical to lexicographic to get the
106  // DoFs correct
107  const unsigned int my_i = lexicographic[i];
108  for (unsigned int q=0; q<n_q_points_1d; ++q)
109  {
110  // fill both vectors with
111  // VectorizedArray<Number>::n_array_elements
112  // copies for the shape information and
113  // non-vectorized fields
114  Point<dim> q_point;
115  q_point[0] = quad.get_points()[q][0];
116  shape_values_number[i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
117  shape_gradient_number[i*n_q_points_1d+q] = fe.shape_grad (my_i,q_point)[0];
118  shape_values [i*n_q_points_1d+q] =
119  shape_values_number [i*n_q_points_1d+q];
120  shape_gradients[i*n_q_points_1d+q] =
121  shape_gradient_number[i*n_q_points_1d+q];
122  shape_hessians[i*n_q_points_1d+q] =
123  fe.shape_grad_grad(my_i,q_point)[0][0];
124  q_point[0] *= 0.5;
125  subface_value[0][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
126  q_point[0] += 0.5;
127  subface_value[1][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
128  }
129  Point<dim> q_point;
130  this->face_value[0][i] = fe.shape_value(my_i,q_point);
131  this->face_gradient[0][i] = fe.shape_grad(my_i,q_point)[0];
132  q_point[0] = 1;
133  this->face_value[1][i] = fe.shape_value(my_i,q_point);
134  this->face_gradient[1][i] = fe.shape_grad(my_i,q_point)[0];
135  }
136 
137  // face information
138  unsigned int n_faces = GeometryInfo<dim>::faces_per_cell;
139  this->face_indices.reinit(n_faces, this->dofs_per_face);
140  switch (dim)
141  {
142  case 3:
143  {
144  for (unsigned int i=0; i<this->dofs_per_face; i++)
145  {
146  const unsigned int jump_term =
147  this->dofs_per_face*((i*n_dofs_1d)/this->dofs_per_face);
148  this->face_indices(0,i) = i*n_dofs_1d;
149  this->face_indices(1,i) = i*n_dofs_1d + n_dofs_1d-1;
150  this->face_indices(2,i) = i%n_dofs_1d + jump_term;
151  this->face_indices(3,i) = (i%n_dofs_1d + jump_term +
152  (n_dofs_1d-1)*n_dofs_1d);
153  this->face_indices(4,i) = i;
154  this->face_indices(5,i) = (n_dofs_1d-1)*this->dofs_per_face+i;
155  }
156  break;
157  }
158  case 2:
159  {
160  for (unsigned int i=0; i<n_dofs_1d; i++)
161  {
162  this->face_indices(0,i) = n_dofs_1d*i;
163  this->face_indices(1,i) = n_dofs_1d*i + n_dofs_1d-1;
164  this->face_indices(2,i) = i;
165  this->face_indices(3,i) = (n_dofs_1d-1)*n_dofs_1d+i;
166  }
167  break;
168  }
169  case 1:
170  {
171  this->face_indices(0,0) = 0;
172  this->face_indices(1,0) = n_dofs_1d-1;
173  break;
174  }
175  default:
176  Assert (false, ExcNotImplemented());
177  }
178  }
179 
180 
181 
182  template <typename Number>
183  std::size_t
185  {
186  std::size_t memory = sizeof(*this);
187  memory += MemoryConsumption::memory_consumption(shape_values);
188  memory += MemoryConsumption::memory_consumption(shape_gradients);
189  memory += MemoryConsumption::memory_consumption(shape_hessians);
190  memory += face_indices.memory_consumption();
191  for (unsigned int i=0; i<2; ++i)
192  {
193  memory += MemoryConsumption::memory_consumption(face_value[i]);
194  memory += MemoryConsumption::memory_consumption(face_gradient[i]);
195  }
196  memory += MemoryConsumption::memory_consumption(shape_values_number);
197  memory += MemoryConsumption::memory_consumption(shape_gradient_number);
198  return memory;
199  }
200 
201  // end of functions for ShapeInfo
202 
203  } // end of namespace MatrixFreeFunctions
204 } // end of namespace internal
205 
206 
207 DEAL_II_NAMESPACE_CLOSE
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
::ExceptionBase & ExcMessage(std::string arg1)
const unsigned int degree
Definition: fe_base.h:287
T fixed_power(const T t)
Definition: utilities.h:705
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
unsigned int n_components() const
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const
void reinit(const Quadrature< 1 > &quad, const FiniteElement< dim > &fe_dim)
const unsigned int dofs_per_cell
Definition: fe_base.h:271
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const
std::vector< unsigned int > get_poly_space_numbering_inverse() const
::ExceptionBase & ExcNotImplemented()
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
::ExceptionBase & ExcInternalError()