Reference documentation for deal.II version 9.3.3
polynomials_bdm.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2020 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
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
21
22#include <iomanip>
23#include <iostream>
24#include <memory>
25
27
28
29template <int dim>
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())
37{
38 switch (dim)
39 {
40 case 2:
42 break;
43 case 3:
44 for (unsigned int i = 0; i < monomials.size(); ++i)
46 break;
47 default:
48 Assert(false, ExcNotImplemented());
49 }
50}
51
52
53
54template <int dim>
55void
57 const Point<dim> & unit_point,
58 std::vector<Tensor<1, dim>> &values,
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()));
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
85 // arrays in the following block
86 // using a mutex to make sure they
87 // are not used by multiple threads
88 // at once
89 {
90 std::lock_guard<std::mutex> lock(mutex);
91
92 p_values.resize((values.size() == 0) ? 0 : n_sub);
95
96 // Compute values of complete space
97 // and insert into tensors. Result
98 // will have first all polynomials
99 // in the x-component, then y and
100 // z.
101 polynomial_space.evaluate(unit_point,
102 p_values,
105 p_third_derivatives,
106 p_fourth_derivatives);
107
108 std::fill(values.begin(), values.end(), Tensor<1, dim>());
109 for (unsigned int i = 0; i < p_values.size(); ++i)
110 for (unsigned int j = 0; j < dim; ++j)
111 values[i + j * n_sub][j] = p_values[i];
112
114 for (unsigned int i = 0; i < p_grads.size(); ++i)
115 for (unsigned int j = 0; j < dim; ++j)
117
120 for (unsigned int j = 0; j < dim; ++j)
122 }
123
124 // This is the first polynomial not
125 // covered by the P_k subspace
126 unsigned int start = dim * n_sub;
127
128 // Store values of auxiliary
129 // polynomials and their three
130 // derivatives
131 std::vector<std::vector<double>> monovali(dim, std::vector<double>(4));
132 std::vector<std::vector<double>> monovalk(dim, std::vector<double>(4));
133
134 if (dim == 2)
135 {
136 for (unsigned int d = 0; d < dim; ++d)
137 monomials[0].value(unit_point(d), monovali[d]);
138 if (values.size() != 0)
139 {
140 values[start][0] = monovali[0][0];
141 values[start][1] = -unit_point(1) * monovali[0][1];
142 values[start + 1][0] = unit_point(0) * monovali[1][1];
143 values[start + 1][1] = -monovali[1][0];
144 }
146 {
149 grads[start][1][0] = -unit_point(1) * monovali[0][2];
151 grads[start + 1][0][0] = monovali[1][1];
152 grads[start + 1][0][1] = unit_point(0) * monovali[1][2];
153 grads[start + 1][1][0] = 0.;
154 grads[start + 1][1][1] = -monovali[1][1];
155 }
157 {
174 }
175 }
176 else // dim == 3
177 {
178 // The number of curls in each
179 // component. Note that the
180 // table in BrezziFortin91 has
181 // a typo, but the text has the
182 // right basis
183
184 // Note that the next basis
185 // function is always obtained
186 // from the previous by cyclic
187 // rotation of the coordinates
188 const unsigned int n_curls = monomials.size() - 1;
189 for (unsigned int i = 0; i < n_curls; ++i, start += dim)
190 {
191 for (unsigned int d = 0; d < dim; ++d)
192 {
193 // p(t) = t^(i+1)
194 monomials[i + 1].value(unit_point(d), monovali[d]);
195 // q(t) = t^(k-i)
196 monomials[this->degree() - 1 - i].value(unit_point(d),
197 monovalk[d]);
198 }
199 if (values.size() != 0)
200 {
201 // x p'(y) q(z)
202 values[start][0] =
203 unit_point(0) * monovali[1][1] * monovalk[2][0];
204 // - p(y) q(z)
205 values[start][1] = -monovali[1][0] * monovalk[2][0];
206 values[start][2] = 0.;
207
208 // y p'(z) q(x)
209 values[start + 1][1] =
210 unit_point(1) * monovali[2][1] * monovalk[0][0];
211 // - p(z) q(x)
212 values[start + 1][2] = -monovali[2][0] * monovalk[0][0];
213 values[start + 1][0] = 0.;
214
215 // z p'(x) q(y)
216 values[start + 2][2] =
217 unit_point(2) * monovali[0][1] * monovalk[1][0];
218 // -p(x) q(y)
219 values[start + 2][0] = -monovali[0][0] * monovalk[1][0];
220 values[start + 2][1] = 0.;
221 }
223 {
224 grads[start][0][0] = monovali[1][1] * monovalk[2][0];
226 unit_point(0) * monovali[1][2] * monovalk[2][0];
228 unit_point(0) * monovali[1][1] * monovalk[2][1];
230 grads[start][1][1] = -monovali[1][1] * monovalk[2][0];
231 grads[start][1][2] = -monovali[1][0] * monovalk[2][1];
235
236 grads[start + 1][1][1] = monovali[2][1] * monovalk[0][0];
238 unit_point(1) * monovali[2][2] * monovalk[0][0];
240 unit_point(1) * monovali[2][1] * monovalk[0][1];
241 grads[start + 1][2][1] = 0.;
242 grads[start + 1][2][2] = -monovali[2][1] * monovalk[0][0];
243 grads[start + 1][2][0] = -monovali[2][0] * monovalk[0][1];
244 grads[start + 1][0][1] = 0.;
245 grads[start + 1][0][2] = 0.;
246 grads[start + 1][0][0] = 0.;
247
248 grads[start + 2][2][2] = monovali[0][1] * monovalk[1][0];
250 unit_point(2) * monovali[0][2] * monovalk[1][0];
252 unit_point(2) * monovali[0][1] * monovalk[1][1];
253 grads[start + 2][0][2] = 0.;
254 grads[start + 2][0][0] = -monovali[0][1] * monovalk[1][0];
255 grads[start + 2][0][1] = -monovali[0][0] * monovalk[1][1];
256 grads[start + 2][1][2] = 0.;
257 grads[start + 2][1][0] = 0.;
258 grads[start + 2][1][1] = 0.;
259 }
261 {
267 unit_point(0) * monovali[1][3] * monovalk[2][0];
269 unit_point(0) * monovali[1][2] * monovalk[2][1];
272 unit_point(0) * monovali[1][2] * monovalk[2][1];
274 unit_point(0) * monovali[1][1] * monovalk[2][2];
293
304 unit_point(1) * monovali[2][1] * monovalk[0][2];
307 unit_point(1) * monovali[2][2] * monovalk[0][1];
312 unit_point(1) * monovalk[0][1] * monovali[2][2];
315 unit_point(1) * monovalk[0][0] * monovali[2][3];
325
345 unit_point(2) * monovali[0][3] * monovalk[1][0];
347 unit_point(2) * monovali[0][2] * monovalk[1][1];
350 unit_point(2) * monovali[0][2] * monovalk[1][1];
352 unit_point(2) * monovali[0][1] * monovalk[1][2];
357 }
358 }
359 Assert(start == this->n(), ExcInternalError());
360 }
361}
362
363
364/*
365template <int dim>
366void
367PolynomialsBDM<dim>::compute_node_matrix (Table<2,double>& A) const
368{
369 std::vector<Polynomial<double> > moment_weight(2);
370 for (unsigned int i=0;i<moment_weight.size();++i)
371 moment_weight[i] = Monomial<double>(i);
372
373 QGauss<dim-1> qface(polynomial_space.degree()+1);
374
375 std::vector<Tensor<1,dim> > values(n());
378 values.resize(n());
379
380 for (unsigned int face=0;face<2*dim;++face)
381 {
382 double orientation = 1.;
383 if ((face==0) || (face==3))
384 orientation = -1.;
385
386 for (unsigned int k=0;k<qface.size();++k)
387 {
388 const double w = qface.weight(k) * orientation;
389 const double x = qface.point(k)(0);
390 Point<dim> p;
391 switch (face)
392 {
393 case 2:
394 p(1) = 1.;
395 case 0:
396 p(0) = x;
397 break;
398 case 1:
399 p(0) = 1.;
400 case 3:
401 p(1) = x;
402 break;
403 }
404// std::cerr << p
405// << '\t' << moment_weight[0].value(x)
406// << '\t' << moment_weight[1].value(x);
407
409
410 for (unsigned int i=0;i<n();++i)
411 {
412// std::cerr << '\t' << std::setw(6) << values[i][1-face%2];
413 // Integrate normal component.
414 // This is easy on the unit
415square for (unsigned int j=0;j<moment_weight.size();++j)
416 A(moment_weight.size()*face+j,i)
417 += w * values[i][1-face%2] * moment_weight[j].value(x);
418 }
419// std::cerr << std::endl;
420 }
421 }
422
423 // Volume integrals are missing
424 //
425 // This degree is one larger
426 Assert (polynomial_space.degree() <= 2,
427 ExcNotImplemented());
428}
429*/
430
431template <int dim>
432unsigned int
434{
435 if (dim == 1)
436 return k + 1;
437 if (dim == 2)
438 return (k + 1) * (k + 2) + 2;
439 if (dim == 3)
440 return ((k + 1) * (k + 2) * (k + 3)) / 2 + 3 * (k + 1);
441 Assert(false, ExcNotImplemented());
442 return 0;
443}
444
445
446template <int dim>
447std::unique_ptr<TensorPolynomialsBase<dim>>
449{
450 return std::make_unique<PolynomialsBDM<dim>>(*this);
451}
452
453
454template class PolynomialsBDM<1>;
455template class PolynomialsBDM<2>;
456template class PolynomialsBDM<3>;
457
458
Definition: point.h:111
std::vector< Polynomials::Polynomial< double > > monomials
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
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
PolynomialsBDM(const unsigned int k)
static unsigned int n_polynomials(const unsigned int degree)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)