Reference documentation for deal.II version GIT 921d917bf4 2023-02-06 18:40:02+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\}}\)
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
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 
16 
19 
21 
22 namespace
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 
42 template <int dim>
44  const unsigned int degree)
45  : ScalarPolynomialsBase<dim>(degree,
46  compute_n_polynomials_pyramid(dim, degree))
47 {}
48 
49 
50 template <int dim>
51 double
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 
88 template <int dim>
91  const Point<dim> & p) const
92 {
93  AssertDimension(dim, 3);
94  AssertIndexRange(this->degree(), 4);
95 
96  Tensor<1, dim> grad;
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  {
150  grad[0] = 0.0;
151  grad[1] = 0.0;
152  grad[2] = 1.0;
153  }
154  else
155  {
156  Assert(false, ExcNotImplemented());
157  }
158  }
159 
160  return grad;
161 }
162 
163 
164 
165 template <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 
180 template <int dim>
181 void
183  const Point<dim> & unit_point,
184  std::vector<double> & values,
185  std::vector<Tensor<1, dim>> &grads,
186  std::vector<Tensor<2, dim>> &grad_grads,
187  std::vector<Tensor<3, dim>> &third_derivatives,
188  std::vector<Tensor<4, dim>> &fourth_derivatives) const
189 {
190  (void)grads;
191  (void)grad_grads;
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 
199  if (grads.size() == this->n())
200  for (unsigned int i = 0; i < this->n(); ++i)
201  grads[i] = compute_grad(i, unit_point);
202 }
203 
204 
205 
206 template <int dim>
209  const unsigned int i,
210  const Point<dim> & p) const
211 {
212  return compute_grad(i, p);
213 }
214 
215 
216 
217 template <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 
233 template <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 
249 template <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 
265 template <int dim>
266 std::string
268 {
269  return "ScalarLagrangePolynomialPyramid";
270 }
271 
272 
273 
274 template <int dim>
275 std::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
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
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1756
#define AssertIndexRange(index, range)
Definition: exceptions.h:1821
Expression fabs(const Expression &x)