Reference documentation for deal.II version GIT d1e6af55b6 2023-02-07 07:50:02+00:00
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
29 template <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
54 template <int dim>
55 void
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
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);
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,
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
108  for (unsigned int i = 0; i < p_grads.size(); ++i)
109  for (unsigned int j = 0; j < dim; ++j)
111
114  for (unsigned int j = 0; j < dim; ++j)
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  }
145  {
148  grads[start][1][0] = -unit_point(1) * monovali[0][2];
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  }
156  {
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
219  {
220  grads[start][0][0] = monovali[1][1] * monovalk[2][0];
222  unit_point(0) * monovali[1][2] * monovalk[2][0];
224  unit_point(0) * monovali[1][1] * monovalk[2][1];
226  grads[start][1][1] = -monovali[1][1] * monovalk[2][0];
227  grads[start][1][2] = -monovali[1][0] * monovalk[2][1];
231
232  grads[start + 1][1][1] = monovali[2][1] * monovalk[0][0];
234  unit_point(1) * monovali[2][2] * monovalk[0][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];
246  unit_point(2) * monovali[0][2] * monovalk[1][0];
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
258  {
264  unit_point(0) * monovali[1][3] * monovalk[2][0];
266  unit_point(0) * monovali[1][2] * monovalk[2][1];
269  unit_point(0) * monovali[1][2] * monovalk[2][1];
271  unit_point(0) * monovali[1][1] * monovalk[2][2];
290
301  unit_point(1) * monovali[2][1] * monovalk[0][2];
304  unit_point(1) * monovali[2][2] * monovalk[0][1];
309  unit_point(1) * monovalk[0][1] * monovali[2][2];
312  unit_point(1) * monovalk[0][0] * monovali[2][3];
322
342  unit_point(2) * monovali[0][3] * monovalk[1][0];
344  unit_point(2) * monovali[0][2] * monovalk[1][1];
347  unit_point(2) * monovali[0][2] * monovalk[1][1];
349  unit_point(2) * monovali[0][1] * monovalk[1][2];
354  }
355  }
356  Assert(start == this->n(), ExcInternalError());
357  }
358 }
359
360
361 template <int dim>
362 unsigned int
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
376 template <int dim>
377 std::unique_ptr<TensorPolynomialsBase<dim>>
379 {
380  return std::make_unique<PolynomialsBDM<dim>>(*this);
381 }
382
383
384 template class PolynomialsBDM<1>;
385 template class PolynomialsBDM<2>;
386 template class PolynomialsBDM<3>;
387
388
Definition: point.h:111
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
std::vector< Polynomials::Polynomial< double > > monomials
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:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)