Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00: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\}}\)
Loading...
Searching...
No Matches
tensor_product_polynomials_bubbles.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 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
19#include <memory>
20
22
23
24
25/* ------------------- TensorProductPolynomialsBubbles -------------- */
26
27
28
29template <int dim>
30void
32{
33 std::array<unsigned int, dim> ix;
34 for (unsigned int i = 0; i < tensor_polys.n(); ++i)
35 {
36 tensor_polys.compute_index(i, ix);
37 out << i << "\t";
38 for (unsigned int d = 0; d < dim; ++d)
39 out << ix[d] << " ";
40 out << std::endl;
41 }
42}
43
44
45
46template <int dim>
47void
49 const std::vector<unsigned int> &renumber)
50{
51 Assert(renumber.size() == index_map.size(),
52 ExcDimensionMismatch(renumber.size(), index_map.size()));
53
54 index_map = renumber;
55 for (unsigned int i = 0; i < index_map.size(); ++i)
56 index_map_inverse[index_map[i]] = i;
57
58 std::vector<unsigned int> renumber_base;
59 for (unsigned int i = 0; i < tensor_polys.n(); ++i)
60 renumber_base.push_back(renumber[i]);
61
62 tensor_polys.set_numbering(renumber_base);
63}
64
65
66template <int dim>
67double
69 const Point<dim> &p) const
70{
71 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
72 const unsigned int max_q_indices = tensor_polys.n();
73 Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
75
76 // treat the regular basis functions
77 if (i < max_q_indices)
78 return tensor_polys.compute_value(i, p);
79
80 const unsigned int comp = i - tensor_polys.n();
81
82 // Compute \prod_{i=1}^d 4*x_i*(1-x_i)
83 double value = 1.;
84 for (unsigned int j = 0; j < dim; ++j)
85 value *= 4 * p[j] * (1 - p[j]);
86
87 // Then multiply with (2x_i-1)^{r-1}. Since q_degree is generally a
88 // small integer, using a loop is likely faster than using std::pow.
89 for (unsigned int i = 0; i < q_degree - 1; ++i)
90 value *= (2 * p[comp] - 1);
91 return value;
92}
93
94
95
96template <>
97double
99 const Point<0> &) const
100{
102 return 0.;
103}
104
105
106
107template <int dim>
110 const Point<dim> &p) const
111{
112 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
113 const unsigned int max_q_indices = tensor_polys.n();
114 Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
116
117 // treat the regular basis functions
118 if (i < max_q_indices)
119 return tensor_polys.compute_grad(i, p);
120
121 const unsigned int comp = i - tensor_polys.n();
122 Tensor<1, dim> grad;
123
124 for (unsigned int d = 0; d < dim; ++d)
125 {
126 grad[d] = 1.;
127 // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
128 for (unsigned j = 0; j < dim; ++j)
129 grad[d] *= (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j]));
130 // and multiply with (2*x_i-1)^{r-1}
131 for (unsigned int i = 0; i < q_degree - 1; ++i)
132 grad[d] *= 2 * p[comp] - 1;
133 }
134
135 if (q_degree >= 2)
136 {
137 // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
138 double value = 1.;
139 for (unsigned int j = 0; j < dim; ++j)
140 value *= 4 * p[j] * (1 - p[j]);
141 // and multiply with grad(2*x_i-1)^{r-1}
142 double tmp = value * 2 * (q_degree - 1);
143 for (unsigned int i = 0; i < q_degree - 2; ++i)
144 tmp *= 2 * p[comp] - 1;
145 grad[comp] += tmp;
146 }
147
148 return grad;
149}
150
151
152
153template <int dim>
156 const unsigned int i,
157 const Point<dim> &p) const
158{
159 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
160 const unsigned int max_q_indices = tensor_polys.n();
161 Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
163
164 // treat the regular basis functions
165 if (i < max_q_indices)
166 return tensor_polys.compute_grad_grad(i, p);
167
168 const unsigned int comp = i - tensor_polys.n();
169
170 double v[dim + 1][3];
171 {
172 for (unsigned int c = 0; c < dim; ++c)
173 {
174 v[c][0] = 4 * p[c] * (1 - p[c]);
175 v[c][1] = 4 * (1 - 2 * p[c]);
176 v[c][2] = -8;
177 }
178
179 double tmp = 1.;
180 for (unsigned int i = 0; i < q_degree - 1; ++i)
181 tmp *= 2 * p[comp] - 1;
182 v[dim][0] = tmp;
183
184 if (q_degree >= 2)
185 {
186 double tmp = 2 * (q_degree - 1);
187 for (unsigned int i = 0; i < q_degree - 2; ++i)
188 tmp *= 2 * p[comp] - 1;
189 v[dim][1] = tmp;
190 }
191 else
192 v[dim][1] = 0.;
193
194 if (q_degree >= 3)
195 {
196 double tmp = 4 * (q_degree - 2) * (q_degree - 1);
197 for (unsigned int i = 0; i < q_degree - 3; ++i)
198 tmp *= 2 * p[comp] - 1;
199 v[dim][2] = tmp;
200 }
201 else
202 v[dim][2] = 0.;
203 }
204
205 // calculate (\partial_j \partial_k \psi) * monomial
206 Tensor<2, dim> grad_grad_1;
207 for (unsigned int d1 = 0; d1 < dim; ++d1)
208 for (unsigned int d2 = 0; d2 < dim; ++d2)
209 {
210 grad_grad_1[d1][d2] = v[dim][0];
211 for (unsigned int x = 0; x < dim; ++x)
212 {
213 unsigned int derivative = 0;
214 if (d1 == x || d2 == x)
215 {
216 if (d1 == d2)
217 derivative = 2;
218 else
219 derivative = 1;
220 }
221 grad_grad_1[d1][d2] *= v[x][derivative];
222 }
223 }
224
225 // calculate (\partial_j \psi) *(\partial_k monomial)
226 // and (\partial_k \psi) *(\partial_j monomial)
227 Tensor<2, dim> grad_grad_2;
228 Tensor<2, dim> grad_grad_3;
229 for (unsigned int d = 0; d < dim; ++d)
230 {
231 grad_grad_2[d][comp] = v[dim][1];
232 grad_grad_3[comp][d] = v[dim][1];
233 for (unsigned int x = 0; x < dim; ++x)
234 {
235 grad_grad_2[d][comp] *= v[x][d == x];
236 grad_grad_3[comp][d] *= v[x][d == x];
237 }
238 }
239
240 // calculate \psi *(\partial j \partial_k monomial) and sum
241 Tensor<2, dim> grad_grad;
242 double psi_value = 1.;
243 for (unsigned int x = 0; x < dim; ++x)
244 psi_value *= v[x][0];
245
246 for (unsigned int d1 = 0; d1 < dim; ++d1)
247 for (unsigned int d2 = 0; d2 < dim; ++d2)
248 grad_grad[d1][d2] =
249 grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
250 grad_grad[comp][comp] += psi_value * v[dim][2];
251
252 return grad_grad;
253}
254
255
256
257template <int dim>
258void
260 const Point<dim> &p,
261 std::vector<double> &values,
262 std::vector<Tensor<1, dim>> &grads,
263 std::vector<Tensor<2, dim>> &grad_grads,
264 std::vector<Tensor<3, dim>> &third_derivatives,
265 std::vector<Tensor<4, dim>> &fourth_derivatives) const
266{
267 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
268 const unsigned int max_q_indices = tensor_polys.n();
269 (void)max_q_indices;
270 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
271 Assert(values.size() == max_q_indices + n_bubbles || values.empty(),
272 ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
273 Assert(grads.size() == max_q_indices + n_bubbles || grads.empty(),
274 ExcDimensionMismatch2(grads.size(), max_q_indices + n_bubbles, 0));
275 Assert(grad_grads.size() == max_q_indices + n_bubbles || grad_grads.empty(),
276 ExcDimensionMismatch2(grad_grads.size(),
277 max_q_indices + n_bubbles,
278 0));
279 Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
280 third_derivatives.empty(),
281 ExcDimensionMismatch2(third_derivatives.size(),
282 max_q_indices + n_bubbles,
283 0));
284 Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
285 fourth_derivatives.empty(),
286 ExcDimensionMismatch2(fourth_derivatives.size(),
287 max_q_indices + n_bubbles,
288 0));
289
290 bool do_values = false, do_grads = false, do_grad_grads = false;
291 bool do_3rd_derivatives = false, do_4th_derivatives = false;
292 if (values.empty() == false)
293 {
294 values.resize(tensor_polys.n());
295 do_values = true;
296 }
297 if (grads.empty() == false)
298 {
299 grads.resize(tensor_polys.n());
300 do_grads = true;
301 }
302 if (grad_grads.empty() == false)
303 {
304 grad_grads.resize(tensor_polys.n());
305 do_grad_grads = true;
306 }
307 if (third_derivatives.empty() == false)
308 {
309 third_derivatives.resize(tensor_polys.n());
310 do_3rd_derivatives = true;
311 }
312 if (fourth_derivatives.empty() == false)
313 {
314 fourth_derivatives.resize(tensor_polys.n());
315 do_4th_derivatives = true;
316 }
317
318 tensor_polys.evaluate(
319 p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
320
321 for (unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i)
322 {
323 if (do_values)
324 values.push_back(compute_value(i, p));
325 if (do_grads)
326 grads.push_back(compute_grad(i, p));
327 if (do_grad_grads)
328 grad_grads.push_back(compute_grad_grad(i, p));
329 if (do_3rd_derivatives)
330 third_derivatives.push_back(compute_derivative<3>(i, p));
331 if (do_4th_derivatives)
332 fourth_derivatives.push_back(compute_derivative<4>(i, p));
333 }
334}
335
336
337
338template <int dim>
339std::unique_ptr<ScalarPolynomialsBase<dim>>
341{
342 return std::make_unique<TensorProductPolynomialsBubbles<dim>>(*this);
343}
344
345
346/* ------------------- explicit instantiations -------------- */
350
Definition point.h:111
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< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
double compute_value(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NOT_IMPLEMENTED()