Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
polynomials_bernardi_raugel.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2018 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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/polynomials_bernardi_raugel.h>
18 #include <deal.II/base/std_cxx14/memory.h>
19 
20 DEAL_II_NAMESPACE_OPEN
21 
22 
23 template <int dim>
25  : TensorPolynomialsBase<dim>(k + 1, n_polynomials(k))
26  , polynomial_space_Q(create_polynomials_Q())
27  , polynomial_space_bubble(create_polynomials_bubble())
28 {}
29 
30 
31 template <int dim>
32 std::vector<std::vector<Polynomials::Polynomial<double>>>
34 {
35  std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
36  std::vector<Polynomials::Polynomial<double>> bubble_shapes;
37  bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
38  bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
39  bubble_shapes.push_back(Polynomials::LagrangeEquidistant(2, 1));
40 
41  for (unsigned int d = 0; d < dim; ++d)
42  pols.push_back(bubble_shapes);
43  // In 2D, the only q_ij polynomials we will use are 31,32,13,23
44  // where ij corresponds to index (i-1)+3*(j-1) (2,5,6,7)
45 
46  // In 3D, the only q_ijk polynomials we will use are 331,332,313,323,133,233
47  // where ijk corresponds to index (i-1)+3*(j-1)+9*(k-1) (8,17,20,23,24,25)
48  return pols;
49 }
50 
51 
52 
53 template <int dim>
54 std::vector<std::vector<Polynomials::Polynomial<double>>>
56 {
57  std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
58  std::vector<Polynomials::Polynomial<double>> Q_shapes;
59  Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
60  Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
61  for (unsigned int d = 0; d < dim; ++d)
62  pols.push_back(Q_shapes);
63 
64  return pols;
65 }
66 
67 
68 template <int dim>
69 void
71  const Point<dim> & unit_point,
72  std::vector<Tensor<1, dim>> &values,
73  std::vector<Tensor<2, dim>> &grads,
74  std::vector<Tensor<3, dim>> &grad_grads,
75  std::vector<Tensor<4, dim>> &third_derivatives,
76  std::vector<Tensor<5, dim>> &fourth_derivatives) const
77 {
78  Assert(values.size() == this->n() || values.size() == 0,
79  ExcDimensionMismatch(values.size(), this->n()));
80  Assert(grads.size() == this->n() || grads.size() == 0,
81  ExcDimensionMismatch(grads.size(), this->n()));
82  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
83  ExcDimensionMismatch(grad_grads.size(), this->n()));
84  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
85  ExcDimensionMismatch(third_derivatives.size(), this->n()));
86  Assert(fourth_derivatives.size() == this->n() ||
87  fourth_derivatives.size() == 0,
88  ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
89 
90  std::vector<double> Q_values;
91  std::vector<Tensor<1, dim>> Q_grads;
92  std::vector<Tensor<2, dim>> Q_grad_grads;
93  std::vector<Tensor<3, dim>> Q_third_derivatives;
94  std::vector<Tensor<4, dim>> Q_fourth_derivatives;
95  std::vector<double> bubble_values;
96  std::vector<Tensor<1, dim>> bubble_grads;
97  std::vector<Tensor<2, dim>> bubble_grad_grads;
98  std::vector<Tensor<3, dim>> bubble_third_derivatives;
99  std::vector<Tensor<4, dim>> bubble_fourth_derivatives;
100 
101  constexpr int n_bubbles =
102  Utilities::pow(3, dim); // size for create_polynomials_bubble
103  constexpr int n_q = 1 << dim; // size for create_polynomials_q
104 
105  // don't resize if the provided vector has 0 length
106  Q_values.resize((values.size() == 0) ? 0 : n_q);
107  Q_grads.resize((grads.size() == 0) ? 0 : n_q);
108  Q_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_q);
109  Q_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_q);
110  Q_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_q);
111  bubble_values.resize((values.size() == 0) ? 0 : n_bubbles);
112  bubble_grads.resize((grads.size() == 0) ? 0 : n_bubbles);
113  bubble_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_bubbles);
114  bubble_third_derivatives.resize((third_derivatives.size() == 0) ? 0 :
115  n_bubbles);
116  bubble_fourth_derivatives.resize(
117  (fourth_derivatives.size() == 0) ? 0 : n_bubbles);
118 
119  // 1 normal vector per face, ordering consistent with GeometryInfo
120  // Normal vectors point in the +x, +y, and +z directions for
121  // consistent orientation across edges
122  std::vector<Tensor<1, dim>> normals;
123  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
124  {
125  Tensor<1, dim> normal;
126  normal[i / 2] = 1;
127  normals.push_back(normal);
128  }
129 
130  // dim standard basis vectors for R^dim, usual ordering
131  std::vector<Tensor<1, dim>> units;
132  for (unsigned int i = 0; i < dim; ++i)
133  {
134  Tensor<1, dim> unit;
135  unit[i] = 1;
136  units.push_back(unit);
137  }
138 
139  // set indices for the anisotropic polynomials to find
140  // them after polynomial_space_bubble.evaluate is called
141  std::vector<int> aniso_indices;
142  if (dim == 2)
143  {
144  aniso_indices.push_back(6);
145  aniso_indices.push_back(7);
146  aniso_indices.push_back(2);
147  aniso_indices.push_back(5);
148  }
149  else if (dim == 3)
150  {
151  aniso_indices.push_back(24);
152  aniso_indices.push_back(25);
153  aniso_indices.push_back(20);
154  aniso_indices.push_back(23);
155  aniso_indices.push_back(8);
156  aniso_indices.push_back(17);
157  }
158 
159  polynomial_space_bubble.evaluate(unit_point,
160  bubble_values,
161  bubble_grads,
162  bubble_grad_grads,
163  bubble_third_derivatives,
164  bubble_fourth_derivatives);
165  polynomial_space_Q.evaluate(unit_point,
166  Q_values,
167  Q_grads,
168  Q_grad_grads,
169  Q_third_derivatives,
170  Q_fourth_derivatives);
171 
172  // first dim*vertices_per_cell functions are Q_1^2 functions
173  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
174  {
175  if (values.size() != 0)
176  {
177  values[i] = units[i % dim] * Q_values[i / dim];
178  }
179  if (grads.size() != 0)
180  {
181  grads[i] = outer_product(units[i % dim], Q_grads[i / dim]);
182  }
183  if (grad_grads.size() != 0)
184  {
185  grad_grads[i] = outer_product(units[i % dim], Q_grad_grads[i / dim]);
186  }
187  if (third_derivatives.size() != 0)
188  {
189  third_derivatives[i] =
190  outer_product(units[i % dim], Q_third_derivatives[i / dim]);
191  }
192  if (fourth_derivatives.size() != 0)
193  {
194  fourth_derivatives[i] =
195  outer_product(units[i % dim], Q_fourth_derivatives[i / dim]);
196  }
197  }
198 
199  // last faces_per_cell functions are bubble functions
200  for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
201  i < dim * GeometryInfo<dim>::vertices_per_cell +
203  ++i)
204  {
205  unsigned int j =
206  i -
207  dim *
208  GeometryInfo<dim>::vertices_per_cell; // ranges 0 to faces_per_cell-1
209  if (values.size() != 0)
210  {
211  values[i] = normals[j] * bubble_values[aniso_indices[j]];
212  }
213  if (grads.size() != 0)
214  {
215  grads[i] = outer_product(normals[j], bubble_grads[aniso_indices[j]]);
216  }
217  if (grad_grads.size() != 0)
218  {
219  grad_grads[i] =
220  outer_product(normals[j], bubble_grad_grads[aniso_indices[j]]);
221  }
222  if (third_derivatives.size() != 0)
223  {
224  third_derivatives[i] =
225  outer_product(normals[j],
226  bubble_third_derivatives[aniso_indices[j]]);
227  }
228  if (fourth_derivatives.size() != 0)
229  {
230  fourth_derivatives[i] =
231  outer_product(normals[j],
232  bubble_fourth_derivatives[aniso_indices[j]]);
233  }
234  }
235 }
236 
237 template <int dim>
238 unsigned int
240 {
241  (void)k;
242  Assert(k == 1, ExcNotImplemented());
243  if (dim == 2 || dim == 3)
246  // 2*4+4=12 polynomials in 2D and 3*8+6=30 polynomials in 3D
247 
248  Assert(false, ExcNotImplemented());
249  return 0;
250 }
251 
252 
253 template <int dim>
254 std::unique_ptr<TensorPolynomialsBase<dim>>
256 {
257  return std_cxx14::make_unique<PolynomialsBernardiRaugel<dim>>(*this);
258 }
259 
260 template class PolynomialsBernardiRaugel<1>; // to prevent errors
261 template class PolynomialsBernardiRaugel<2>;
262 template class PolynomialsBernardiRaugel<3>;
263 
264 
265 DEAL_II_NAMESPACE_CLOSE
void evaluate(const Point< dim > &unit_point, std::vector< Tensor< 1, dim >> &values, std::vector< Tensor< 2, dim >> &grads, std::vector< Tensor< 3, dim >> &grad_grads, std::vector< Tensor< 4, dim >> &third_derivatives, std::vector< Tensor< 5, dim >> &fourth_derivatives) const override
const AnisotropicPolynomials< dim > polynomial_space_Q
static unsigned int n_polynomials(const unsigned int k)
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials_Q()
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:430
DEAL_II_CONSTEXPR SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials_bubble()
PolynomialsBernardiRaugel(const unsigned int k)
static ::ExceptionBase & ExcNotImplemented()
const AnisotropicPolynomials< dim > polynomial_space_bubble