deal.II version GIT relicensing-1927-g3de9220933 2024-10-03 08:40: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_pyramid.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 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
18
20
21namespace
22{
23 unsigned int
24 compute_n_polynomials_pyramid(const unsigned int dim,
25 const unsigned int degree)
26 {
27 if (dim == 3)
28 {
29 if (degree == 1)
30 return 5;
31 }
32
34
35 return 0;
36 }
37} // namespace
38
39
40
41template <int dim>
43 const unsigned int degree)
44 : ScalarPolynomialsBase<dim>(degree,
45 compute_n_polynomials_pyramid(dim, degree))
46{}
47
48
49template <int dim>
50double
52 const Point<dim> &p) const
53{
54 AssertDimension(dim, 3);
55 AssertIndexRange(this->degree(), 2);
56
57 const double Q14 = 0.25;
58 double ration;
59
60 const double r = p[0];
61 const double s = p[1];
62 const double t = p[2];
63
64 if (fabs(t - 1.0) > 1.0e-14)
65 {
66 ration = (r * s * t) / (1.0 - t);
67 }
68 else
69 {
70 ration = 0.0;
71 }
72
73 if (i == 0)
74 return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
75 if (i == 1)
76 return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
77 if (i == 2)
78 return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
79 if (i == 3)
80 return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
81 else
82 return t;
83}
84
85
86
87template <int dim>
90 const Point<dim> &p) const
91{
92 AssertDimension(dim, 3);
93 AssertIndexRange(this->degree(), 4);
94
95 Tensor<1, dim> grad;
96
97 if (this->degree() == 1)
98 {
99 const double Q14 = 0.25;
100
101 const double r = p[0];
102 const double s = p[1];
103 const double t = p[2];
104
105 double rationdr;
106 double rationds;
107 double rationdt;
108
109 if (fabs(t - 1.0) > 1.0e-14)
110 {
111 rationdr = s * t / (1.0 - t);
112 rationds = r * t / (1.0 - t);
113 rationdt = r * s / ((1.0 - t) * (1.0 - t));
114 }
115 else
116 {
117 rationdr = 1.0;
118 rationds = 1.0;
119 rationdt = 1.0;
120 }
121
122
123 if (i == 0)
124 {
125 grad[0] = Q14 * (-1.0 * (1.0 - s) + rationdr);
126 grad[1] = Q14 * (-1.0 * (1.0 - r) + rationds);
127 grad[2] = Q14 * (rationdt - 1.0);
128 }
129 else if (i == 1)
130 {
131 grad[0] = Q14 * (1.0 * (1.0 - s) - rationdr);
132 grad[1] = Q14 * (-1.0 * (1.0 + r) - rationds);
133 grad[2] = Q14 * (-1.0 * rationdt - 1.0);
134 }
135 else if (i == 2)
136 {
137 grad[0] = Q14 * (-1.0 * (1.0 + s) - rationdr);
138 grad[1] = Q14 * (1.0 * (1.0 - r) - rationds);
139 grad[2] = Q14 * (-1.0 * rationdt - 1.0);
140 }
141 else if (i == 3)
142 {
143 grad[0] = Q14 * (1.0 * (1.0 + s) + rationdr);
144 grad[1] = Q14 * (1.0 * (1.0 + r) + rationds);
145 grad[2] = Q14 * (rationdt - 1.0);
146 }
147 else if (i == 4)
148 {
149 grad[0] = 0.0;
150 grad[1] = 0.0;
151 grad[2] = 1.0;
152 }
153 else
154 {
156 }
157 }
158
159 return grad;
160}
161
162
163
164template <int dim>
167 const unsigned int i,
168 const Point<dim> &p) const
169{
170 (void)i;
171 (void)p;
172
174 return Tensor<2, dim>();
175}
176
177
178
179template <int dim>
180void
182 const Point<dim> &unit_point,
183 std::vector<double> &values,
184 std::vector<Tensor<1, dim>> &grads,
185 std::vector<Tensor<2, dim>> &grad_grads,
186 std::vector<Tensor<3, dim>> &third_derivatives,
187 std::vector<Tensor<4, dim>> &fourth_derivatives) const
188{
189 (void)grads;
190 (void)grad_grads;
191 (void)third_derivatives;
192 (void)fourth_derivatives;
193
194 if (values.size() == this->n())
195 for (unsigned int i = 0; i < this->n(); ++i)
196 values[i] = compute_value(i, unit_point);
197
198 if (grads.size() == this->n())
199 for (unsigned int i = 0; i < this->n(); ++i)
200 grads[i] = compute_grad(i, unit_point);
201}
202
203
204
205template <int dim>
208 const unsigned int i,
209 const Point<dim> &p) const
210{
211 return compute_grad(i, p);
212}
213
214
215
216template <int dim>
219 const unsigned int i,
220 const Point<dim> &p) const
221{
222 (void)i;
223 (void)p;
224
226
227 return {};
228}
229
230
231
232template <int dim>
235 const unsigned int i,
236 const Point<dim> &p) const
237{
238 (void)i;
239 (void)p;
240
242
243 return {};
244}
245
246
247
248template <int dim>
251 const unsigned int i,
252 const Point<dim> &p) const
253{
254 (void)i;
255 (void)p;
256
258
259 return {};
260}
261
262
263
264template <int dim>
265std::string
267{
268 return "ScalarLagrangePolynomialPyramid";
269}
270
271
272
273template <int dim>
274std::unique_ptr<ScalarPolynomialsBase<dim>>
276{
277 return std::make_unique<ScalarLagrangePolynomialPyramid<dim>>(*this);
278}
279
280
281
285
Definition point.h:111
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
ScalarLagrangePolynomialPyramid(const unsigned int degree)
Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
double compute_value(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
std::string name() const override
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DEAL_II_NOT_IMPLEMENTED()