Reference documentation for deal.II version 9.4.0
polynomials_pyramid.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 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
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
19
21
22namespace
23{
24 unsigned int
25 compute_n_polynomials_pyramid(const unsigned int dim,
26 const unsigned int degree)
27 {
28 if (dim == 3)
29 {
30 if (degree == 1)
31 return 5;
32 }
33
34 Assert(false, ExcNotImplemented());
35
36 return 0;
37 }
38} // namespace
39
40
41
42template <int dim>
44 const unsigned int degree)
45 : ScalarPolynomialsBase<dim>(degree,
46 compute_n_polynomials_pyramid(dim, degree))
47{}
48
49
50template <int dim>
51double
53 const Point<dim> & p) const
54{
55 AssertDimension(dim, 3);
56 AssertIndexRange(this->degree(), 2);
57
58 const double Q14 = 0.25;
59 double ration;
60
61 const double r = p[0];
62 const double s = p[1];
63 const double t = p[2];
64
65 if (fabs(t - 1.0) > 1.0e-14)
66 {
67 ration = (r * s * t) / (1.0 - t);
68 }
69 else
70 {
71 ration = 0.0;
72 }
73
74 if (i == 0)
75 return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
76 if (i == 1)
77 return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
78 if (i == 2)
79 return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
80 if (i == 3)
81 return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
82 else
83 return t;
84}
85
86
87
88template <int dim>
91 const Point<dim> & p) const
92{
93 AssertDimension(dim, 3);
94 AssertIndexRange(this->degree(), 4);
95
97
98 if (this->degree() == 1)
99 {
100 const double Q14 = 0.25;
101
102 const double r = p[0];
103 const double s = p[1];
104 const double t = p[2];
105
106 double rationdr;
107 double rationds;
108 double rationdt;
109
110 if (fabs(t - 1.0) > 1.0e-14)
111 {
112 rationdr = s * t / (1.0 - t);
113 rationds = r * t / (1.0 - t);
114 rationdt = r * s / ((1.0 - t) * (1.0 - t));
115 }
116 else
117 {
118 rationdr = 1.0;
119 rationds = 1.0;
120 rationdt = 1.0;
121 }
122
123
124 if (i == 0)
125 {
126 grad[0] = Q14 * (-1.0 * (1.0 - s) + rationdr);
127 grad[1] = Q14 * (-1.0 * (1.0 - r) + rationds);
128 grad[2] = Q14 * (rationdt - 1.0);
129 }
130 else if (i == 1)
131 {
132 grad[0] = Q14 * (1.0 * (1.0 - s) - rationdr);
133 grad[1] = Q14 * (-1.0 * (1.0 + r) - rationds);
134 grad[2] = Q14 * (-1.0 * rationdt - 1.0);
135 }
136 else if (i == 2)
137 {
138 grad[0] = Q14 * (-1.0 * (1.0 + s) - rationdr);
139 grad[1] = Q14 * (1.0 * (1.0 - r) - rationds);
140 grad[2] = Q14 * (-1.0 * rationdt - 1.0);
141 }
142 else if (i == 3)
143 {
144 grad[0] = Q14 * (1.0 * (1.0 + s) + rationdr);
145 grad[1] = Q14 * (1.0 * (1.0 + r) + rationds);
146 grad[2] = Q14 * (rationdt - 1.0);
147 }
148 else if (i == 4)
149 {
153 }
154 else
155 {
156 Assert(false, ExcNotImplemented());
157 }
158 }
159
161}
162
163
164
165template <int dim>
168 const unsigned int i,
169 const Point<dim> & p) const
170{
171 (void)i;
172 (void)p;
173
174 Assert(false, ExcNotImplemented());
175 return Tensor<2, dim>();
176}
177
178
179
180template <int dim>
181void
183 const Point<dim> & unit_point,
184 std::vector<double> & values,
187 std::vector<Tensor<3, dim>> &third_derivatives,
188 std::vector<Tensor<4, dim>> &fourth_derivatives) const
189{
192 (void)third_derivatives;
193 (void)fourth_derivatives;
194
195 if (values.size() == this->n())
196 for (unsigned int i = 0; i < this->n(); ++i)
197 values[i] = compute_value(i, unit_point);
198
200 for (unsigned int i = 0; i < this->n(); ++i)
202}
203
204
205
206template <int dim>
209 const unsigned int i,
210 const Point<dim> & p) const
211{
213}
214
215
216
217template <int dim>
220 const unsigned int i,
221 const Point<dim> & p) const
222{
223 (void)i;
224 (void)p;
225
226 Assert(false, ExcNotImplemented());
227
228 return {};
229}
230
231
232
233template <int dim>
236 const unsigned int i,
237 const Point<dim> & p) const
238{
239 (void)i;
240 (void)p;
241
242 Assert(false, ExcNotImplemented());
243
244 return {};
245}
246
247
248
249template <int dim>
252 const unsigned int i,
253 const Point<dim> & p) const
254{
255 (void)i;
256 (void)p;
257
258 Assert(false, ExcNotImplemented());
259
260 return {};
261}
262
263
264
265template <int dim>
266std::string
268{
269 return "ScalarLagrangePolynomialPyramid";
270}
271
272
273
274template <int dim>
275std::unique_ptr<ScalarPolynomialsBase<dim>>
277{
278 return std::make_unique<ScalarLagrangePolynomialPyramid<dim>>(*this);
279}
280
281
282
286
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
Expression fabs(const Expression &x)