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