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