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