deal.II version GIT relicensing-2167-g9622207b8f 2024-11-21 12: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_bdm.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 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
16#include <deal.II/base/geometry_info.h>
17#include <deal.II/base/polynomial_space.h>
18#include <deal.II/base/polynomials_bdm.h>
19#include <deal.II/base/quadrature_lib.h>
20
21#include <iomanip>
22#include <iostream>
23#include <memory>
24
25DEAL_II_NAMESPACE_OPEN
26
27
28template <int dim>
29PolynomialsBDM<dim>::PolynomialsBDM(const unsigned int k)
30 : TensorPolynomialsBase<dim>(k + 1, n_polynomials(k))
31 , polynomial_space(Polynomials::Legendre::generate_complete_basis(k))
32 , monomials((dim == 2) ? (1) : (k + 2))
33 , p_values(polynomial_space.n())
34 , p_grads(polynomial_space.n())
35 , p_grad_grads(polynomial_space.n())
36{
37 switch (dim)
38 {
39 case 2:
40 monomials[0] = Polynomials::Monomial<double>(k + 1);
41 break;
42 case 3:
43 for (unsigned int i = 0; i < monomials.size(); ++i)
44 monomials[i] = Polynomials::Monomial<double>(i);
45 break;
46 default:
47 DEAL_II_NOT_IMPLEMENTED();
48 }
49}
50
51
52
53template <int dim>
54void
55PolynomialsBDM<dim>::evaluate(
56 const Point<dim> &unit_point,
57 std::vector<Tensor<1, dim>> &values,
58 std::vector<Tensor<2, dim>> &grads,
59 std::vector<Tensor<3, dim>> &grad_grads,
60 std::vector<Tensor<4, dim>> &third_derivatives,
61 std::vector<Tensor<5, dim>> &fourth_derivatives) const
62{
63 Assert(values.size() == this->n() || values.empty(),
64 ExcDimensionMismatch(values.size(), this->n()));
65 Assert(grads.size() == this->n() || grads.empty(),
66 ExcDimensionMismatch(grads.size(), this->n()));
67 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
68 ExcDimensionMismatch(grad_grads.size(), this->n()));
69 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
70 ExcDimensionMismatch(third_derivatives.size(), this->n()));
71 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
72 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
73
74 // third and fourth derivatives not implemented
75 (void)third_derivatives;
76 Assert(third_derivatives.empty(), ExcNotImplemented());
77 (void)fourth_derivatives;
78 Assert(fourth_derivatives.empty(), ExcNotImplemented());
79
80 const unsigned int n_sub = polynomial_space.n();
81
82 // guard access to the scratch arrays in the following block using a
83 // mutex to make sure they are not used by multiple threads at once
84 {
85 std::lock_guard<std::mutex> lock(mutex);
86
87 p_values.resize((values.empty()) ? 0 : n_sub);
88 p_grads.resize((grads.empty()) ? 0 : n_sub);
89 p_grad_grads.resize((grad_grads.empty()) ? 0 : n_sub);
90
91 // Compute values of complete space and insert into tensors. Result
92 // will have first all polynomials in the x-component, then y and z.
93 polynomial_space.evaluate(unit_point,
94 p_values,
95 p_grads,
96 p_grad_grads,
97 p_third_derivatives,
98 p_fourth_derivatives);
99
100 std::fill(values.begin(), values.end(), Tensor<1, dim>());
101 for (unsigned int i = 0; i < p_values.size(); ++i)
102 for (unsigned int j = 0; j < dim; ++j)
103 values[i + j * n_sub][j] = p_values[i];
104
105 std::fill(grads.begin(), grads.end(), Tensor<2, dim>());
106 for (unsigned int i = 0; i < p_grads.size(); ++i)
107 for (unsigned int j = 0; j < dim; ++j)
108 grads[i + j * n_sub][j] = p_grads[i];
109
110 std::fill(grad_grads.begin(), grad_grads.end(), Tensor<3, dim>());
111 for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
112 for (unsigned int j = 0; j < dim; ++j)
113 grad_grads[i + j * n_sub][j] = p_grad_grads[i];
114 }
115
116 // This is the first polynomial not covered by the P_k subspace
117 unsigned int start = dim * n_sub;
118
119 // Store values of auxiliary polynomials and their three derivatives
120 std::vector<std::vector<double>> monovali(dim, std::vector<double>(4));
121 std::vector<std::vector<double>> monovalk(dim, std::vector<double>(4));
122
123 if (dim == 1)
124 {
125 // Despite the fact that we are instantiating this class for 1, 2 and
126 // 3 space dimensions we only support dimension 2 and 3.
127 Assert(false,
128 ::ExcMessage("PolynomialsBDF::evaluate is only "
129 "available for dim == 2, or dim == 3"));
130 }
131 else if (dim == 2)
132 {
133 for (unsigned int d = 0; d < dim; ++d)
134 monomials[0].value(unit_point[d], monovali[d]);
135 if (values.size() != 0)
136 {
137 values[start][0] = monovali[0][0];
138 values[start][1] = -unit_point[1] * monovali[0][1];
139 values[start + 1][0] = unit_point[0] * monovali[1][1];
140 values[start + 1][1] = -monovali[1][0];
141 }
142 if (grads.size() != 0)
143 {
144 grads[start][0][0] = monovali[0][1];
145 grads[start][0][1] = 0.;
146 grads[start][1][0] = -unit_point[1] * monovali[0][2];
147 grads[start][1][1] = -monovali[0][1];
148 grads[start + 1][0][0] = monovali[1][1];
149 grads[start + 1][0][1] = unit_point[0] * monovali[1][2];
150 grads[start + 1][1][0] = 0.;
151 grads[start + 1][1][1] = -monovali[1][1];
152 }
153 if (grad_grads.size() != 0)
154 {
155 grad_grads[start][0][0][0] = monovali[0][2];
156 grad_grads[start][0][0][1] = 0.;
157 grad_grads[start][0][1][0] = 0.;
158 grad_grads[start][0][1][1] = 0.;
159 grad_grads[start][1][0][0] = -unit_point[1] * monovali[0][3];
160 grad_grads[start][1][0][1] = -monovali[0][2];
161 grad_grads[start][1][1][0] = -monovali[0][2];
162 grad_grads[start][1][1][1] = 0.;
163 grad_grads[start + 1][0][0][0] = 0;
164 grad_grads[start + 1][0][0][1] = monovali[1][2];
165 grad_grads[start + 1][0][1][0] = monovali[1][2];
166 grad_grads[start + 1][0][1][1] = unit_point[0] * monovali[1][3];
167 grad_grads[start + 1][1][0][0] = 0.;
168 grad_grads[start + 1][1][0][1] = 0.;
169 grad_grads[start + 1][1][1][0] = 0.;
170 grad_grads[start + 1][1][1][1] = -monovali[1][2];
171 }
172 }
173 else if (dim == 3)
174 {
175 // The number of curls in each component. Note that the table in
176 // BrezziFortin91 has a typo, but the text has the right basis
177
178 // Note that the next basis function is always obtained from the
179 // previous by cyclic rotation of the coordinates
180 const unsigned int n_curls = monomials.size() - 1;
181 for (unsigned int i = 0; i < n_curls; ++i, start += dim)
182 {
183 for (unsigned int d = 0; d < dim; ++d)
184 {
185 // p(t) = t^(i+1)
186 monomials[i + 1].value(unit_point[d], monovali[d]);
187 // q(t) = t^(k-i)
188 monomials[this->degree() - 1 - i].value(unit_point[d],
189 monovalk[d]);
190 }
191
192 if (values.size() != 0)
193 {
194 // x p'(y) q(z)
195 values[start][0] =
196 unit_point[0] * monovali[1][1] * monovalk[2][0];
197 // - p(y) q(z)
198 values[start][1] = -monovali[1][0] * monovalk[2][0];
199 values[start][2] = 0.;
200
201 // y p'(z) q(x)
202 values[start + 1][1] =
203 unit_point[1] * monovali[2][1] * monovalk[0][0];
204 // - p(z) q(x)
205 values[start + 1][2] = -monovali[2][0] * monovalk[0][0];
206 values[start + 1][0] = 0.;
207
208 // z p'(x) q(y)
209 values[start + 2][2] =
210 unit_point[2] * monovali[0][1] * monovalk[1][0];
211 // -p(x) q(y)
212 values[start + 2][0] = -monovali[0][0] * monovalk[1][0];
213 values[start + 2][1] = 0.;
214 }
215
216 if (grads.size() != 0)
217 {
218 grads[start][0][0] = monovali[1][1] * monovalk[2][0];
219 grads[start][0][1] =
220 unit_point[0] * monovali[1][2] * monovalk[2][0];
221 grads[start][0][2] =
222 unit_point[0] * monovali[1][1] * monovalk[2][1];
223 grads[start][1][0] = 0.;
224 grads[start][1][1] = -monovali[1][1] * monovalk[2][0];
225 grads[start][1][2] = -monovali[1][0] * monovalk[2][1];
226 grads[start][2][0] = 0.;
227 grads[start][2][1] = 0.;
228 grads[start][2][2] = 0.;
229
230 grads[start + 1][1][1] = monovali[2][1] * monovalk[0][0];
231 grads[start + 1][1][2] =
232 unit_point[1] * monovali[2][2] * monovalk[0][0];
233 grads[start + 1][1][0] =
234 unit_point[1] * monovali[2][1] * monovalk[0][1];
235 grads[start + 1][2][1] = 0.;
236 grads[start + 1][2][2] = -monovali[2][1] * monovalk[0][0];
237 grads[start + 1][2][0] = -monovali[2][0] * monovalk[0][1];
238 grads[start + 1][0][1] = 0.;
239 grads[start + 1][0][2] = 0.;
240 grads[start + 1][0][0] = 0.;
241
242 grads[start + 2][2][2] = monovali[0][1] * monovalk[1][0];
243 grads[start + 2][2][0] =
244 unit_point[2] * monovali[0][2] * monovalk[1][0];
245 grads[start + 2][2][1] =
246 unit_point[2] * monovali[0][1] * monovalk[1][1];
247 grads[start + 2][0][2] = 0.;
248 grads[start + 2][0][0] = -monovali[0][1] * monovalk[1][0];
249 grads[start + 2][0][1] = -monovali[0][0] * monovalk[1][1];
250 grads[start + 2][1][2] = 0.;
251 grads[start + 2][1][0] = 0.;
252 grads[start + 2][1][1] = 0.;
253 }
254
255 if (grad_grads.size() != 0)
256 {
257 grad_grads[start][0][0][0] = 0.;
258 grad_grads[start][0][0][1] = monovali[1][2] * monovalk[2][0];
259 grad_grads[start][0][0][2] = monovali[1][1] * monovalk[2][1];
260 grad_grads[start][0][1][0] = monovali[1][2] * monovalk[2][0];
261 grad_grads[start][0][1][1] =
262 unit_point[0] * monovali[1][3] * monovalk[2][0];
263 grad_grads[start][0][1][2] =
264 unit_point[0] * monovali[1][2] * monovalk[2][1];
265 grad_grads[start][0][2][0] = monovali[1][1] * monovalk[2][1];
266 grad_grads[start][0][2][1] =
267 unit_point[0] * monovali[1][2] * monovalk[2][1];
268 grad_grads[start][0][2][2] =
269 unit_point[0] * monovali[1][1] * monovalk[2][2];
270 grad_grads[start][1][0][0] = 0.;
271 grad_grads[start][1][0][1] = 0.;
272 grad_grads[start][1][0][2] = 0.;
273 grad_grads[start][1][1][0] = 0.;
274 grad_grads[start][1][1][1] = -monovali[1][2] * monovalk[2][0];
275 grad_grads[start][1][1][2] = -monovali[1][1] * monovalk[2][1];
276 grad_grads[start][1][2][0] = 0.;
277 grad_grads[start][1][2][1] = -monovali[1][1] * monovalk[2][1];
278 grad_grads[start][1][2][2] = -monovali[1][0] * monovalk[2][2];
279 grad_grads[start][2][0][0] = 0.;
280 grad_grads[start][2][0][1] = 0.;
281 grad_grads[start][2][0][2] = 0.;
282 grad_grads[start][2][1][0] = 0.;
283 grad_grads[start][2][1][1] = 0.;
284 grad_grads[start][2][1][2] = 0.;
285 grad_grads[start][2][2][0] = 0.;
286 grad_grads[start][2][2][1] = 0.;
287 grad_grads[start][2][2][2] = 0.;
288
289 grad_grads[start + 1][0][0][0] = 0.;
290 grad_grads[start + 1][0][0][1] = 0.;
291 grad_grads[start + 1][0][0][2] = 0.;
292 grad_grads[start + 1][0][1][0] = 0.;
293 grad_grads[start + 1][0][1][1] = 0.;
294 grad_grads[start + 1][0][1][2] = 0.;
295 grad_grads[start + 1][0][2][0] = 0.;
296 grad_grads[start + 1][0][2][1] = 0.;
297 grad_grads[start + 1][0][2][2] = 0.;
298 grad_grads[start + 1][1][0][0] =
299 unit_point[1] * monovali[2][1] * monovalk[0][2];
300 grad_grads[start + 1][1][0][1] = monovali[2][1] * monovalk[0][1];
301 grad_grads[start + 1][1][0][2] =
302 unit_point[1] * monovali[2][2] * monovalk[0][1];
303 grad_grads[start + 1][1][1][0] = monovalk[0][1] * monovali[2][1];
304 grad_grads[start + 1][1][1][1] = 0.;
305 grad_grads[start + 1][1][1][2] = monovalk[0][0] * monovali[2][2];
306 grad_grads[start + 1][1][2][0] =
307 unit_point[1] * monovalk[0][1] * monovali[2][2];
308 grad_grads[start + 1][1][2][1] = monovalk[0][0] * monovali[2][2];
309 grad_grads[start + 1][1][2][2] =
310 unit_point[1] * monovalk[0][0] * monovali[2][3];
311 grad_grads[start + 1][2][0][0] = -monovalk[0][2] * monovali[2][0];
312 grad_grads[start + 1][2][0][1] = 0.;
313 grad_grads[start + 1][2][0][2] = -monovalk[0][1] * monovali[2][1];
314 grad_grads[start + 1][2][1][0] = 0.;
315 grad_grads[start + 1][2][1][1] = 0.;
316 grad_grads[start + 1][2][1][2] = 0.;
317 grad_grads[start + 1][2][2][0] = -monovalk[0][1] * monovali[2][1];
318 grad_grads[start + 1][2][2][1] = 0.;
319 grad_grads[start + 1][2][2][2] = -monovalk[0][0] * monovali[2][2];
320
321 grad_grads[start + 2][0][0][0] = -monovali[0][2] * monovalk[1][0];
322 grad_grads[start + 2][0][0][1] = -monovali[0][1] * monovalk[1][1];
323 grad_grads[start + 2][0][0][2] = 0.;
324 grad_grads[start + 2][0][1][0] = -monovali[0][1] * monovalk[1][1];
325 grad_grads[start + 2][0][1][1] = -monovali[0][0] * monovalk[1][2];
326 grad_grads[start + 2][0][1][2] = 0.;
327 grad_grads[start + 2][0][2][0] = 0.;
328 grad_grads[start + 2][0][2][1] = 0.;
329 grad_grads[start + 2][0][2][2] = 0.;
330 grad_grads[start + 2][1][0][0] = 0.;
331 grad_grads[start + 2][1][0][1] = 0.;
332 grad_grads[start + 2][1][0][2] = 0.;
333 grad_grads[start + 2][1][1][0] = 0.;
334 grad_grads[start + 2][1][1][1] = 0.;
335 grad_grads[start + 2][1][1][2] = 0.;
336 grad_grads[start + 2][1][2][0] = 0.;
337 grad_grads[start + 2][1][2][1] = 0.;
338 grad_grads[start + 2][1][2][2] = 0.;
339 grad_grads[start + 2][2][0][0] =
340 unit_point[2] * monovali[0][3] * monovalk[1][0];
341 grad_grads[start + 2][2][0][1] =
342 unit_point[2] * monovali[0][2] * monovalk[1][1];
343 grad_grads[start + 2][2][0][2] = monovali[0][2] * monovalk[1][0];
344 grad_grads[start + 2][2][1][0] =
345 unit_point[2] * monovali[0][2] * monovalk[1][1];
346 grad_grads[start + 2][2][1][1] =
347 unit_point[2] * monovali[0][1] * monovalk[1][2];
348 grad_grads[start + 2][2][1][2] = monovali[0][1] * monovalk[1][1];
349 grad_grads[start + 2][2][2][0] = monovali[0][2] * monovalk[1][0];
350 grad_grads[start + 2][2][2][1] = monovali[0][1] * monovalk[1][1];
351 grad_grads[start + 2][2][2][2] = 0.;
352 }
353 }
354 Assert(start == this->n(), ExcInternalError());
355 }
356}
357
358
359template <int dim>
360unsigned int
361PolynomialsBDM<dim>::n_polynomials(const unsigned int k)
362{
363 if (dim == 1)
364 return k + 1;
365 if (dim == 2)
366 return (k + 1) * (k + 2) + 2;
367 if (dim == 3)
368 return ((k + 1) * (k + 2) * (k + 3)) / 2 + 3 * (k + 1);
369 DEAL_II_NOT_IMPLEMENTED();
370 return 0;
371}
372
373
374template <int dim>
375std::unique_ptr<TensorPolynomialsBase<dim>>
376PolynomialsBDM<dim>::clone() const
377{
378 return std::make_unique<PolynomialsBDM<dim>>(*this);
379}
380
381
382template class PolynomialsBDM<1>;
383template class PolynomialsBDM<2>;
384template class PolynomialsBDM<3>;
385
386
387DEAL_II_NAMESPACE_CLOSE
Point
Definition point.h:111
PolynomialsBDM::monomials
std::vector< Polynomials::Polynomial< double > > monomials
Definition polynomials_bdm.h:157
PolynomialsBDM::evaluate
void evaluate(const Point< dim > &unit_point, std::vector< Tensor< 1, dim > > &values, std::vector< Tensor< 2, dim > > &grads, std::vector< Tensor< 3, dim > > &grad_grads, std::vector< Tensor< 4, dim > > &third_derivatives, std::vector< Tensor< 5, dim > > &fourth_derivatives) const override
Definition polynomials_bdm.cc:55
PolynomialsBDM::clone
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
Definition polynomials_bdm.cc:376
PolynomialsBDM::PolynomialsBDM
PolynomialsBDM(const unsigned int k)
Definition polynomials_bdm.cc:29
PolynomialsBDM::n_polynomials
static unsigned int n_polynomials(const unsigned int degree)
Definition polynomials_bdm.cc:361
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
geometry_info.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
polynomial_space.h
polynomials_bdm.h
quadrature_lib.h
Generated by doxygen 1.9.8