Reference documentation for deal.II version GIT d1e6af55b6 2023-02-07 07:50:02+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\}}\)
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
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 
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())
35  , p_grads(polynomial_space.n())
36  , p_grad_grads(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,
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 
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)