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