Reference documentation for deal.II version Git 932f7faded 2020-11-28 20:02:43 +0100
\(\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_adini.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 
18 
19 #include <memory>
20 
21 #define ENTER_COEFFICIENTS( \
22  koefs, z, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11) \
23  koefs(0, z) = a0; \
24  koefs(1, z) = a1; \
25  koefs(2, z) = a2; \
26  koefs(3, z) = a3; \
27  koefs(4, z) = a4; \
28  koefs(5, z) = a5; \
29  koefs(6, z) = a6; \
30  koefs(7, z) = a7; \
31  koefs(8, z) = a8; \
32  koefs(9, z) = a9; \
33  koefs(10, z) = a10; \
34  koefs(11, z) = a11;
35 
36 
38 
39 
40 
41 template <int dim>
43  : ScalarPolynomialsBase<dim>(3, 12)
44  , coef(12, 12)
45  , dx(12, 12)
46  , dy(12, 12)
47  , dxx(12, 12)
48  , dyy(12, 12)
49  , dxy(12, 12)
50 {
51  Assert(dim == 2, ExcNotImplemented());
52 
53  // 1 x y xx yy xy 3x 3y xyy xxy 3xy x3y
54  // 0 1 2 3 4 5 6 7 8 9 10 11
55  ENTER_COEFFICIENTS(coef, 0, 1, 0, 0, -3, -3, -1, 2, 2, 3, 3, -2, -2);
56  ENTER_COEFFICIENTS(coef, 1, 0, 1, 0, -2, 0, -1, 1, 0, 0, 2, -1, 0);
57  ENTER_COEFFICIENTS(coef, 2, 0, 0, 1, 0, -2, -1, 0, 1, 2, 0, 0, -1);
58  ENTER_COEFFICIENTS(coef, 3, 0, 0, 0, 3, 0, 1, -2, 0, -3, -3, 2, 2);
59  ENTER_COEFFICIENTS(coef, 4, 0, 0, 0, -1, 0, 0, 1, 0, 0, 1, -1, 0);
60  ENTER_COEFFICIENTS(coef, 5, 0, 0, 0, 0, 0, 1, 0, 0, -2, 0, 0, 1);
61  ENTER_COEFFICIENTS(coef, 6, 0, 0, 0, 0, 3, 1, 0, -2, -3, -3, 2, 2);
62  ENTER_COEFFICIENTS(coef, 7, 0, 0, 0, 0, 0, 1, 0, 0, 0, -2, 1, 0);
63  ENTER_COEFFICIENTS(coef, 8, 0, 0, 0, 0, -1, 0, 0, 1, 1, 0, 0, -1);
64  ENTER_COEFFICIENTS(coef, 9, 0, 0, 0, 0, 0, -1, 0, 0, 3, 3, -2, -2);
65  ENTER_COEFFICIENTS(coef, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0);
66  ENTER_COEFFICIENTS(coef, 11, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1);
67 
68  ENTER_COEFFICIENTS(dx, 0, 0, -6, -1, 6, 3, 6, 0, -2, 0, -6, 0, 0);
69  ENTER_COEFFICIENTS(dx, 1, 1, -4, -1, 3, 0, 4, 0, 0, 0, -3, 0, 0);
70  ENTER_COEFFICIENTS(dx, 2, 0, 0, -1, 0, 2, 0, 0, -1, 0, 0, 0, 0);
71  ENTER_COEFFICIENTS(dx, 3, 0, 6, 1, -6, -3, -6, 0, 2, 0, 6, 0, 0);
72  ENTER_COEFFICIENTS(dx, 4, 0, -2, 0, 3, 0, 2, 0, 0, 0, -3, 0, 0);
73  ENTER_COEFFICIENTS(dx, 5, 0, 0, 1, 0, -2, 0, 0, 1, 0, 0, 0, 0);
74  ENTER_COEFFICIENTS(dx, 6, 0, 0, 1, 0, -3, -6, 0, 2, 0, 6, 0, 0);
75  ENTER_COEFFICIENTS(dx, 7, 0, 0, 1, 0, 0, -4, 0, 0, 0, 3, 0, 0);
76  ENTER_COEFFICIENTS(dx, 8, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0);
77  ENTER_COEFFICIENTS(dx, 9, 0, 0, -1, 0, 3, 6, 0, -2, 0, -6, 0, 0);
78  ENTER_COEFFICIENTS(dx, 10, 0, 0, 0, 0, 0, -2, 0, 0, 0, 3, 0, 0);
79  ENTER_COEFFICIENTS(dx, 11, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0);
80 
81  ENTER_COEFFICIENTS(dy, 0, 0, -1, -6, 3, 6, 6, -2, 0, -6, 0, 0, 0);
82  ENTER_COEFFICIENTS(dy, 1, 0, -1, 0, 2, 0, 0, -1, 0, 0, 0, 0, 0);
83  ENTER_COEFFICIENTS(dy, 2, 1, -1, -4, 0, 3, 4, 0, 0, -3, 0, 0, 0);
84  ENTER_COEFFICIENTS(dy, 3, 0, 1, 0, -3, 0, -6, 2, 0, 6, 0, 0, 0);
85  ENTER_COEFFICIENTS(dy, 4, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0);
86  ENTER_COEFFICIENTS(dy, 5, 0, 1, 0, 0, 0, -4, 0, 0, 3, 0, 0, 0);
87  ENTER_COEFFICIENTS(dy, 6, 0, 1, 6, -3, -6, -6, 2, 0, 6, 0, 0, 0);
88  ENTER_COEFFICIENTS(dy, 7, 0, 1, 0, -2, 0, 0, 1, 0, 0, 0, 0, 0);
89  ENTER_COEFFICIENTS(dy, 8, 0, 0, -2, 0, 3, 2, 0, 0, -3, 0, 0, 0);
90  ENTER_COEFFICIENTS(dy, 9, 0, -1, 0, 3, 0, 6, -2, 0, -6, 0, 0, 0);
91  ENTER_COEFFICIENTS(dy, 10, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0);
92  ENTER_COEFFICIENTS(dy, 11, 0, 0, 0, 0, 0, -2, 0, 0, 3, 0, 0, 0);
93 
94  ENTER_COEFFICIENTS(dxx, 0, -6, 12, 6, 0, 0, -12, 0, 0, 0, 0, 0, 0);
95  ENTER_COEFFICIENTS(dxx, 1, -4, 6, 4, 0, 0, -6, 0, 0, 0, 0, 0, 0);
96  ENTER_COEFFICIENTS(dxx, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
97  ENTER_COEFFICIENTS(dxx, 3, 6, -12, -6, 0, 0, 12, 0, 0, 0, 0, 0, 0);
98  ENTER_COEFFICIENTS(dxx, 4, -2, 6, 2, 0, 0, -6, 0, 0, 0, 0, 0, 0);
99  ENTER_COEFFICIENTS(dxx, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
100  ENTER_COEFFICIENTS(dxx, 6, 0, 0, -6, 0, 0, 12, 0, 0, 0, 0, 0, 0);
101  ENTER_COEFFICIENTS(dxx, 7, 0, 0, -4, 0, 0, 6, 0, 0, 0, 0, 0, 0);
102  ENTER_COEFFICIENTS(dxx, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
103  ENTER_COEFFICIENTS(dxx, 9, 0, 0, 6, 0, 0, -12, 0, 0, 0, 0, 0, 0);
104  ENTER_COEFFICIENTS(dxx, 10, 0, 0, -2, 0, 0, 6, 0, 0, 0, 0, 0, 0);
105  ENTER_COEFFICIENTS(dxx, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
106 
107  ENTER_COEFFICIENTS(dyy, 0, -6, 6, 12, 0, 0, -12, 0, 0, 0, 0, 0, 0);
108  ENTER_COEFFICIENTS(dyy, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
109  ENTER_COEFFICIENTS(dyy, 2, -4, 4, 6, 0, 0, -6, 0, 0, 0, 0, 0, 0);
110  ENTER_COEFFICIENTS(dyy, 3, 0, -6, 0, 0, 0, 12, 0, 0, 0, 0, 0, 0);
111  ENTER_COEFFICIENTS(dyy, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
112  ENTER_COEFFICIENTS(dyy, 5, 0, -4, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0);
113  ENTER_COEFFICIENTS(dyy, 6, 6, -6, -12, 0, 0, 12, 0, 0, 0, 0, 0, 0);
114  ENTER_COEFFICIENTS(dyy, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
115  ENTER_COEFFICIENTS(dyy, 8, -2, 2, 6, 0, 0, -6, 0, -0, 0, 0, 0, 0);
116  ENTER_COEFFICIENTS(dyy, 9, 0, 6, 0, 0, 0, -12, 0, 0, 0, 0, 0, 0);
117  ENTER_COEFFICIENTS(dyy, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
118  ENTER_COEFFICIENTS(dyy, 11, 0, -2, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0);
119 
120  ENTER_COEFFICIENTS(dxy, 0, -1, 6, 6, -6, -6, 0, 0, 0, 0, 0, 0, 0);
121  ENTER_COEFFICIENTS(dxy, 1, -1, 4, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0);
122  ENTER_COEFFICIENTS(dxy, 2, -1, 0, 4, 0, -3, 0, 0, 0, 0, 0, 0, 0);
123  ENTER_COEFFICIENTS(dxy, 3, 1, -6, -6, 6, 6, 0, 0, 0, 0, 0, 0, 0);
124  ENTER_COEFFICIENTS(dxy, 4, 0, 2, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0);
125  ENTER_COEFFICIENTS(dxy, 5, 1, 0, -4, 0, 3, 0, 0, 0, 0, 0, 0, 0);
126  ENTER_COEFFICIENTS(dxy, 6, 1, -6, -6, 6, 6, 0, 0, 0, 0, 0, 0, 0);
127  ENTER_COEFFICIENTS(dxy, 7, 1, -4, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0);
128  ENTER_COEFFICIENTS(dxy, 8, 0, 0, 2, 0, -3, 0, 0, 0, 0, 0, 0, 0);
129  ENTER_COEFFICIENTS(dxy, 9, -1, 6, 6, -6, -6, 0, 0, 0, 0, 0, 0, 0);
130  ENTER_COEFFICIENTS(dxy, 10, 0, -2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0);
131  ENTER_COEFFICIENTS(dxy, 11, 0, 0, -2, 0, 3, 0, 0, 0, 0, 0, 0, 0);
132 }
133 
134 
135 
136 template <int dim>
137 void
139  const Point<dim> & unit_point,
140  std::vector<double> & values,
141  std::vector<Tensor<1, dim>> &grads,
142  std::vector<Tensor<2, dim>> &grad_grads,
143  std::vector<Tensor<3, dim>> &third_derivatives,
144  std::vector<Tensor<4, dim>> &fourth_derivatives) const
145 {
146  const unsigned int n_pols = this->n();
147  (void)n_pols;
148 
149  Assert(values.size() == n_pols || values.size() == 0,
150  ExcDimensionMismatch(values.size(), n_pols));
151  Assert(grads.size() == n_pols || grads.size() == 0,
152  ExcDimensionMismatch(grads.size(), n_pols));
153  Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
154  ExcDimensionMismatch(grad_grads.size(), n_pols));
155  (void)third_derivatives;
156  Assert(third_derivatives.size() == n_pols || third_derivatives.size() == 0,
157  ExcDimensionMismatch(third_derivatives.size(), n_pols));
158  (void)fourth_derivatives;
159  Assert(fourth_derivatives.size() == n_pols || fourth_derivatives.size() == 0,
160  ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
161 
162  if (values.empty() == false) // do not bother if empty
163  {
164  for (unsigned int i = 0; i < values.size(); ++i)
165  {
166  values[i] = compute_value(i, unit_point);
167  }
168  }
169 
170  if (grads.empty() == false) // do not bother if empty
171  {
172  for (unsigned int i = 0; i < grads.size(); ++i)
173  {
174  grads[i] = compute_grad(i, unit_point);
175  }
176  }
177 
178  if (grad_grads.empty() == false) // do not bother if empty
179  {
180  for (unsigned int i = 0; i < grad_grads.size(); ++i)
181  {
182  grad_grads[i] = compute_grad_grad(i, unit_point);
183  }
184  }
185 
186  return;
187 }
188 
189 
190 
191 template <int dim>
192 double
194  const Point<dim> & p) const
195 {
196  const double x = p(0);
197  const double y = p(1);
198  return coef(0, i) + coef(1, i) * x + coef(2, i) * y + coef(3, i) * x * x +
199  coef(4, i) * y * y + coef(5, i) * x * y + coef(6, i) * x * x * x +
200  coef(7, i) * y * y * y + coef(8, i) * x * y * y +
201  coef(9, i) * x * x * y + coef(10, i) * x * x * x * y +
202  coef(11, i) * x * y * y * y;
203 }
204 
205 
206 
207 template <int dim>
210  const Point<dim> & p) const
211 {
212  const double x = p(0);
213  const double y = p(1);
214  Tensor<1, dim> tensor;
215  tensor[0] = dx(0, i) + dx(1, i) * x + dx(2, i) * y + dx(3, i) * x * x +
216  dx(4, i) * y * y + dx(5, i) * x * y + dx(6, i) * x * x * x +
217  dx(7, i) * y * y * y + dx(8, i) * x * y * y +
218  dx(9, i) * x * x * y + dx(10, i) * x * x * x * y +
219  dx(11, i) * x * y * y * y;
220 
221  tensor[1] = dy(0, i) + dy(1, i) * x + dy(2, i) * y + dy(3, i) * x * x +
222  dy(4, i) * y * y + dy(5, i) * x * y + dy(6, i) * x * x * x +
223  dy(7, i) * y * y * y + dy(8, i) * x * y * y +
224  dy(9, i) * x * x * y + dy(10, i) * x * x * x * y +
225  dy(11, i) * x * y * y * y;
226  return tensor;
227 }
228 
229 
230 
231 template <int dim>
234  const Point<dim> & p) const
235 {
236  const double x = p(0);
237  const double y = p(1);
238  Tensor<2, dim> tensor;
239  tensor[0][0] = dxx(0, i) + dxx(1, i) * x + dxx(2, i) * y + dxx(3, i) * x * x +
240  dxx(4, i) * y * y + dxx(5, i) * x * y + dxx(6, i) * x * x * x +
241  dxx(7, i) * y * y * y + dxx(8, i) * x * y * y +
242  dxx(9, i) * x * x * y + dxx(10, i) * x * x * x * y +
243  dxx(11, i) * x * y * y * y;
244  tensor[0][1] = dxy(0, i) + dxy(1, i) * x + dxy(2, i) * y + dxy(3, i) * x * x +
245  dxy(4, i) * y * y + dxy(5, i) * x * y + dxy(6, i) * x * x * x +
246  dxy(7, i) * y * y * y + dxy(8, i) * x * y * y +
247  dxy(9, i) * x * x * y + dxy(10, i) * x * x * x * y +
248  dxy(11, i) * x * y * y * y;
249  tensor[1][0] = tensor[0][1];
250  tensor[1][1] = dyy(0, i) + dyy(1, i) * x + dyy(2, i) * y + dyy(3, i) * x * x +
251  dyy(4, i) * y * y + dyy(5, i) * x * y + dyy(6, i) * x * x * x +
252  dyy(7, i) * y * y * y + dyy(8, i) * x * y * y +
253  dyy(9, i) * x * x * y + dyy(10, i) * x * x * x * y +
254  dyy(11, i) * x * y * y * y;
255  return tensor;
256 }
257 
258 
259 
260 template <int dim>
261 std::unique_ptr<ScalarPolynomialsBase<dim>>
263 {
264  return std::make_unique<PolynomialsAdini<dim>>(*this);
265 }
266 
267 
268 
269 template class PolynomialsAdini<0>;
270 template class PolynomialsAdini<1>;
271 template class PolynomialsAdini<2>;
272 template class PolynomialsAdini<3>;
273 
Table< 2, double > dyy
Table< 2, double > dx
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
#define ENTER_COEFFICIENTS( koefs, z, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11)
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Table< 2, double > coef
double compute_value(const unsigned int i, const Point< dim > &p) const override
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Table< 2, double > dxx
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
Table< 2, double > dy
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
Table< 2, double > dxy
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
static ::ExceptionBase & ExcNotImplemented()
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim >> &grads, std::vector< Tensor< 2, dim >> &grad_grads, std::vector< Tensor< 3, dim >> &third_derivatives, std::vector< Tensor< 4, dim >> &fourth_derivatives) const override