Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2012 - 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
20#include <memory>
21
23
24
25
26/* ------------------- TensorProductPolynomialsBubbles -------------- */
27
28
29
30template <int dim>
31void
33{
34 std::array<unsigned int, dim> ix;
35 for (unsigned int i = 0; i < tensor_polys.n(); ++i)
36 {
37 tensor_polys.compute_index(i, ix);
38 out << i << "\t";
39 for (unsigned int d = 0; d < dim; ++d)
40 out << ix[d] << " ";
41 out << std::endl;
42 }
43}
44
45
46
47template <int dim>
48void
50 const std::vector<unsigned int> &renumber)
51{
52 Assert(renumber.size() == index_map.size(),
53 ExcDimensionMismatch(renumber.size(), index_map.size()));
54
55 index_map = renumber;
56 for (unsigned int i = 0; i < index_map.size(); ++i)
57 index_map_inverse[index_map[i]] = i;
58
59 std::vector<unsigned int> renumber_base;
60 for (unsigned int i = 0; i < tensor_polys.n(); ++i)
61 renumber_base.push_back(renumber[i]);
62
63 tensor_polys.set_numbering(renumber_base);
64}
65
66
67template <int dim>
68double
70 const Point<dim> & p) const
71{
72 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
73 const unsigned int max_q_indices = tensor_polys.n();
74 Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
76
77 // treat the regular basis functions
78 if (i < max_q_indices)
79 return tensor_polys.compute_value(i, p);
80
81 const unsigned int comp = i - tensor_polys.n();
82
83 // Compute \prod_{i=1}^d 4*x_i*(1-x_i)
84 double value = 1.;
85 for (unsigned int j = 0; j < dim; ++j)
86 value *= 4 * p(j) * (1 - p(j));
87
88 // Then multiply with (2x_i-1)^{r-1}. Since q_degree is generally a
89 // small integer, using a loop is likely faster than using std::pow.
90 for (unsigned int i = 0; i < q_degree - 1; ++i)
91 value *= (2 * p(comp) - 1);
92 return value;
93}
94
95
96
97template <>
98double
100 const Point<0> &) const
101{
102 Assert(false, ExcNotImplemented());
103 return 0.;
104}
105
106
107
108template <int dim>
111 const Point<dim> & p) const
112{
113 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
114 const unsigned int max_q_indices = tensor_polys.n();
115 Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
117
118 // treat the regular basis functions
119 if (i < max_q_indices)
120 return tensor_polys.compute_grad(i, p);
121
122 const unsigned int comp = i - tensor_polys.n();
123 Tensor<1, dim> grad;
124
125 for (unsigned int d = 0; d < dim; ++d)
126 {
127 grad[d] = 1.;
128 // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
129 for (unsigned j = 0; j < dim; ++j)
130 grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
131 // and multiply with (2*x_i-1)^{r-1}
132 for (unsigned int i = 0; i < q_degree - 1; ++i)
133 grad[d] *= 2 * p(comp) - 1;
134 }
135
136 if (q_degree >= 2)
137 {
138 // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
139 double value = 1.;
140 for (unsigned int j = 0; j < dim; ++j)
141 value *= 4 * p(j) * (1 - p(j));
142 // and multiply with grad(2*x_i-1)^{r-1}
143 double tmp = value * 2 * (q_degree - 1);
144 for (unsigned int i = 0; i < q_degree - 2; ++i)
145 tmp *= 2 * p(comp) - 1;
146 grad[comp] += tmp;
147 }
148
149 return grad;
150}
151
152
153
154template <int dim>
157 const unsigned int i,
158 const Point<dim> & p) const
159{
160 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
161 const unsigned int max_q_indices = tensor_polys.n();
162 Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
164
165 // treat the regular basis functions
166 if (i < max_q_indices)
167 return tensor_polys.compute_grad_grad(i, p);
168
169 const unsigned int comp = i - tensor_polys.n();
170
171 double v[dim + 1][3];
172 {
173 for (unsigned int c = 0; c < dim; ++c)
174 {
175 v[c][0] = 4 * p(c) * (1 - p(c));
176 v[c][1] = 4 * (1 - 2 * p(c));
177 v[c][2] = -8;
178 }
179
180 double tmp = 1.;
181 for (unsigned int i = 0; i < q_degree - 1; ++i)
182 tmp *= 2 * p(comp) - 1;
183 v[dim][0] = tmp;
184
185 if (q_degree >= 2)
186 {
187 double tmp = 2 * (q_degree - 1);
188 for (unsigned int i = 0; i < q_degree - 2; ++i)
189 tmp *= 2 * p(comp) - 1;
190 v[dim][1] = tmp;
191 }
192 else
193 v[dim][1] = 0.;
194
195 if (q_degree >= 3)
196 {
197 double tmp = 4 * (q_degree - 2) * (q_degree - 1);
198 for (unsigned int i = 0; i < q_degree - 3; ++i)
199 tmp *= 2 * p(comp) - 1;
200 v[dim][2] = tmp;
201 }
202 else
203 v[dim][2] = 0.;
204 }
205
206 // calculate (\partial_j \partial_k \psi) * monomial
207 Tensor<2, dim> grad_grad_1;
208 for (unsigned int d1 = 0; d1 < dim; ++d1)
209 for (unsigned int d2 = 0; d2 < dim; ++d2)
210 {
211 grad_grad_1[d1][d2] = v[dim][0];
212 for (unsigned int x = 0; x < dim; ++x)
213 {
214 unsigned int derivative = 0;
215 if (d1 == x || d2 == x)
216 {
217 if (d1 == d2)
218 derivative = 2;
219 else
220 derivative = 1;
221 }
222 grad_grad_1[d1][d2] *= v[x][derivative];
223 }
224 }
225
226 // calculate (\partial_j \psi) *(\partial_k monomial)
227 // and (\partial_k \psi) *(\partial_j monomial)
228 Tensor<2, dim> grad_grad_2;
229 Tensor<2, dim> grad_grad_3;
230 for (unsigned int d = 0; d < dim; ++d)
231 {
232 grad_grad_2[d][comp] = v[dim][1];
233 grad_grad_3[comp][d] = v[dim][1];
234 for (unsigned int x = 0; x < dim; ++x)
235 {
236 grad_grad_2[d][comp] *= v[x][d == x];
237 grad_grad_3[comp][d] *= v[x][d == x];
238 }
239 }
240
241 // calculate \psi *(\partial j \partial_k monomial) and sum
242 Tensor<2, dim> grad_grad;
243 double psi_value = 1.;
244 for (unsigned int x = 0; x < dim; ++x)
245 psi_value *= v[x][0];
246
247 for (unsigned int d1 = 0; d1 < dim; ++d1)
248 for (unsigned int d2 = 0; d2 < dim; ++d2)
249 grad_grad[d1][d2] =
250 grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
251 grad_grad[comp][comp] += psi_value * v[dim][2];
252
253 return grad_grad;
254}
255
256
257
258template <int dim>
259void
261 const Point<dim> & p,
262 std::vector<double> & values,
263 std::vector<Tensor<1, dim>> &grads,
264 std::vector<Tensor<2, dim>> &grad_grads,
265 std::vector<Tensor<3, dim>> &third_derivatives,
266 std::vector<Tensor<4, dim>> &fourth_derivatives) const
267{
268 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
269 const unsigned int max_q_indices = tensor_polys.n();
270 (void)max_q_indices;
271 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
272 Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
273 ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
274 Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
275 ExcDimensionMismatch2(grads.size(), max_q_indices + n_bubbles, 0));
276 Assert(
277 grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
278 ExcDimensionMismatch2(grad_grads.size(), max_q_indices + n_bubbles, 0));
279 Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
280 third_derivatives.size() == 0,
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.size() == 0,
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:112
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#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)