Reference documentation for deal.II version 9.1.0-pre
polynomials_bdm.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2015 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 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 template <int dim>
29 PolynomialsBDM<dim>::PolynomialsBDM(const unsigned int k)
30  : polynomial_space(Polynomials::Legendre::generate_complete_basis(k))
31  , monomials((dim == 2) ? (1) : (k + 2))
32  , n_pols(compute_n_pols(k))
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  Assert(false, ExcNotImplemented());
48  }
49 }
50 
51 
52 
53 template <int dim>
54 void
55 PolynomialsBDM<dim>::compute(
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() == n_pols || values.size() == 0,
64  ExcDimensionMismatch(values.size(), n_pols));
65  Assert(grads.size() == n_pols || grads.size() == 0,
66  ExcDimensionMismatch(grads.size(), n_pols));
67  Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
68  ExcDimensionMismatch(grad_grads.size(), n_pols));
69  Assert(third_derivatives.size() == n_pols || third_derivatives.size() == 0,
70  ExcDimensionMismatch(third_derivatives.size(), n_pols));
71  Assert(fourth_derivatives.size() == n_pols || fourth_derivatives.size() == 0,
72  ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
73 
74  // third and fourth derivatives not implemented
75  (void)third_derivatives;
76  Assert(third_derivatives.size() == 0, ExcNotImplemented());
77  (void)fourth_derivatives;
78  Assert(fourth_derivatives.size() == 0, ExcNotImplemented());
79 
80  const unsigned int n_sub = polynomial_space.n();
81 
82  // guard access to the scratch
83  // arrays in the following block
84  // using a mutex to make sure they
85  // are not used by multiple threads
86  // at once
87  {
88  Threads::Mutex::ScopedLock lock(mutex);
89 
90  p_values.resize((values.size() == 0) ? 0 : n_sub);
91  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
92  p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
93 
94  // Compute values of complete space
95  // and insert into tensors. Result
96  // will have first all polynomials
97  // in the x-component, then y and
98  // z.
99  polynomial_space.compute(unit_point,
100  p_values,
101  p_grads,
102  p_grad_grads,
103  p_third_derivatives,
104  p_fourth_derivatives);
105 
106  std::fill(values.begin(), values.end(), Tensor<1, dim>());
107  for (unsigned int i = 0; i < p_values.size(); ++i)
108  for (unsigned int j = 0; j < dim; ++j)
109  values[i + j * n_sub][j] = p_values[i];
110 
111  std::fill(grads.begin(), grads.end(), Tensor<2, dim>());
112  for (unsigned int i = 0; i < p_grads.size(); ++i)
113  for (unsigned int j = 0; j < dim; ++j)
114  grads[i + j * n_sub][j] = p_grads[i];
115 
116  std::fill(grad_grads.begin(), grad_grads.end(), Tensor<3, dim>());
117  for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
118  for (unsigned int j = 0; j < dim; ++j)
119  grad_grads[i + j * n_sub][j] = p_grad_grads[i];
120  }
121 
122  // This is the first polynomial not
123  // covered by the P_k subspace
124  unsigned int start = dim * n_sub;
125 
126  // Store values of auxiliary
127  // polynomials and their three
128  // derivatives
129  std::vector<std::vector<double>> monovali(dim, std::vector<double>(4));
130  std::vector<std::vector<double>> monovalk(dim, std::vector<double>(4));
131 
132  if (dim == 2)
133  {
134  for (unsigned int d = 0; d < dim; ++d)
135  monomials[0].value(unit_point(d), monovali[d]);
136  if (values.size() != 0)
137  {
138  values[start][0] = monovali[0][0];
139  values[start][1] = -unit_point(1) * monovali[0][1];
140  values[start + 1][0] = unit_point(0) * monovali[1][1];
141  values[start + 1][1] = -monovali[1][0];
142  }
143  if (grads.size() != 0)
144  {
145  grads[start][0][0] = monovali[0][1];
146  grads[start][0][1] = 0.;
147  grads[start][1][0] = -unit_point(1) * monovali[0][2];
148  grads[start][1][1] = -monovali[0][1];
149  grads[start + 1][0][0] = monovali[1][1];
150  grads[start + 1][0][1] = unit_point(0) * monovali[1][2];
151  grads[start + 1][1][0] = 0.;
152  grads[start + 1][1][1] = -monovali[1][1];
153  }
154  if (grad_grads.size() != 0)
155  {
156  grad_grads[start][0][0][0] = monovali[0][2];
157  grad_grads[start][0][0][1] = 0.;
158  grad_grads[start][0][1][0] = 0.;
159  grad_grads[start][0][1][1] = 0.;
160  grad_grads[start][1][0][0] = -unit_point(1) * monovali[0][3];
161  grad_grads[start][1][0][1] = -monovali[0][2];
162  grad_grads[start][1][1][0] = -monovali[0][2];
163  grad_grads[start][1][1][1] = 0.;
164  grad_grads[start + 1][0][0][0] = 0;
165  grad_grads[start + 1][0][0][1] = monovali[1][2];
166  grad_grads[start + 1][0][1][0] = monovali[1][2];
167  grad_grads[start + 1][0][1][1] = unit_point(0) * monovali[1][3];
168  grad_grads[start + 1][1][0][0] = 0.;
169  grad_grads[start + 1][1][0][1] = 0.;
170  grad_grads[start + 1][1][1][0] = 0.;
171  grad_grads[start + 1][1][1][1] = -monovali[1][2];
172  }
173  }
174  else // dim == 3
175  {
176  // The number of curls in each
177  // component. Note that the
178  // table in BrezziFortin91 has
179  // a typo, but the text has the
180  // right basis
181 
182  // Note that the next basis
183  // function is always obtained
184  // from the previous by cyclic
185  // rotation of the coordinates
186  const unsigned int n_curls = monomials.size() - 1;
187  for (unsigned int i = 0; i < n_curls; ++i, start += dim)
188  {
189  for (unsigned int d = 0; d < dim; ++d)
190  {
191  // p(t) = t^(i+1)
192  monomials[i + 1].value(unit_point(d), monovali[d]);
193  // q(t) = t^(k-i)
194  monomials[degree() - i].value(unit_point(d), monovalk[d]);
195  }
196  if (values.size() != 0)
197  {
198  // x p'(y) q(z)
199  values[start][0] =
200  unit_point(0) * monovali[1][1] * monovalk[2][0];
201  // - p(y) q(z)
202  values[start][1] = -monovali[1][0] * monovalk[2][0];
203  values[start][2] = 0.;
204 
205  // y p'(z) q(x)
206  values[start + 1][1] =
207  unit_point(1) * monovali[2][1] * monovalk[0][0];
208  // - p(z) q(x)
209  values[start + 1][2] = -monovali[2][0] * monovalk[0][0];
210  values[start + 1][0] = 0.;
211 
212  // z p'(x) q(y)
213  values[start + 2][2] =
214  unit_point(2) * monovali[0][1] * monovalk[1][0];
215  // -p(x) q(y)
216  values[start + 2][0] = -monovali[0][0] * monovalk[1][0];
217  values[start + 2][1] = 0.;
218  }
219  if (grads.size() != 0)
220  {
221  grads[start][0][0] = monovali[1][1] * monovalk[2][0];
222  grads[start][0][1] =
223  unit_point(0) * monovali[1][2] * monovalk[2][0];
224  grads[start][0][2] =
225  unit_point(0) * monovali[1][1] * monovalk[2][1];
226  grads[start][1][0] = 0.;
227  grads[start][1][1] = -monovali[1][1] * monovalk[2][0];
228  grads[start][1][2] = -monovali[1][0] * monovalk[2][1];
229  grads[start][2][0] = 0.;
230  grads[start][2][1] = 0.;
231  grads[start][2][2] = 0.;
232 
233  grads[start + 1][1][1] = monovali[2][1] * monovalk[0][0];
234  grads[start + 1][1][2] =
235  unit_point(1) * monovali[2][2] * monovalk[0][0];
236  grads[start + 1][1][0] =
237  unit_point(1) * monovali[2][1] * monovalk[0][1];
238  grads[start + 1][2][1] = 0.;
239  grads[start + 1][2][2] = -monovali[2][1] * monovalk[0][0];
240  grads[start + 1][2][0] = -monovali[2][0] * monovalk[0][1];
241  grads[start + 1][0][1] = 0.;
242  grads[start + 1][0][2] = 0.;
243  grads[start + 1][0][0] = 0.;
244 
245  grads[start + 2][2][2] = monovali[0][1] * monovalk[1][0];
246  grads[start + 2][2][0] =
247  unit_point(2) * monovali[0][2] * monovalk[1][0];
248  grads[start + 2][2][1] =
249  unit_point(2) * monovali[0][1] * monovalk[1][1];
250  grads[start + 2][0][2] = 0.;
251  grads[start + 2][0][0] = -monovali[0][1] * monovalk[1][0];
252  grads[start + 2][0][1] = -monovali[0][0] * monovalk[1][1];
253  grads[start + 2][1][2] = 0.;
254  grads[start + 2][1][0] = 0.;
255  grads[start + 2][1][1] = 0.;
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 == n_pols, ExcInternalError());
357  }
358 }
359 
360 
361 /*
362 template <int dim>
363 void
364 PolynomialsBDM<dim>::compute_node_matrix (Table<2,double>& A) const
365 {
366  std::vector<Polynomial<double> > moment_weight(2);
367  for (unsigned int i=0;i<moment_weight.size();++i)
368  moment_weight[i] = Monomial<double>(i);
369 
370  QGauss<dim-1> qface(polynomial_space.degree()+1);
371 
372  std::vector<Tensor<1,dim> > values(n());
373  std::vector<Tensor<2,dim> > grads;
374  std::vector<Tensor<3,dim> > grad_grads;
375  values.resize(n());
376 
377  for (unsigned int face=0;face<2*dim;++face)
378  {
379  double orientation = 1.;
380  if ((face==0) || (face==3))
381  orientation = -1.;
382 
383  for (unsigned int k=0;k<qface.size();++k)
384  {
385  const double w = qface.weight(k) * orientation;
386  const double x = qface.point(k)(0);
387  Point<dim> p;
388  switch (face)
389  {
390  case 2:
391  p(1) = 1.;
392  case 0:
393  p(0) = x;
394  break;
395  case 1:
396  p(0) = 1.;
397  case 3:
398  p(1) = x;
399  break;
400  }
401 // std::cerr << p
402 // << '\t' << moment_weight[0].value(x)
403 // << '\t' << moment_weight[1].value(x);
404 
405  compute (p, values, grads, grad_grads);
406 
407  for (unsigned int i=0;i<n();++i)
408  {
409 // std::cerr << '\t' << std::setw(6) << values[i][1-face%2];
410  // Integrate normal component.
411  // This is easy on the unit
412 square for (unsigned int j=0;j<moment_weight.size();++j)
413  A(moment_weight.size()*face+j,i)
414  += w * values[i][1-face%2] * moment_weight[j].value(x);
415  }
416 // std::cerr << std::endl;
417  }
418  }
419 
420  // Volume integrals are missing
421  //
422  // This degree is one larger
423  Assert (polynomial_space.degree() <= 2,
424  ExcNotImplemented());
425 }
426 */
427 
428 template <int dim>
429 unsigned int
430 PolynomialsBDM<dim>::compute_n_pols(unsigned int k)
431 {
432  if (dim == 1)
433  return k + 1;
434  if (dim == 2)
435  return (k + 1) * (k + 2) + 2;
436  if (dim == 3)
437  return ((k + 1) * (k + 2) * (k + 3)) / 2 + 3 * (k + 1);
438  Assert(false, ExcNotImplemented());
439  return 0;
440 }
441 
442 
443 template class PolynomialsBDM<1>;
444 template class PolynomialsBDM<2>;
445 template class PolynomialsBDM<3>;
446 
447 
448 DEAL_II_NAMESPACE_CLOSE
PolynomialsBDM::PolynomialsBDM
PolynomialsBDM(const unsigned int k)
Definition: polynomials_bdm.cc:29
PolynomialsBDM::monomials
std::vector< Polynomials::Polynomial< double > > monomials
Definition: polynomials_bdm.h:172
PolynomialsBDM::degree
unsigned int degree() const
Definition: polynomials_bdm.h:221
PolynomialsBDM::p_third_derivatives
std::vector< Tensor< 3, dim > > p_third_derivatives
Definition: polynomials_bdm.h:202
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1227
StandardExceptions::ExcDimensionMismatch
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
PolynomialsBDM::p_values
std::vector< double > p_values
Definition: polynomials_bdm.h:187
PolynomialsBDM::compute
void compute(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
Definition: polynomials_bdm.cc:55
PolynomialsBDM::p_grads
std::vector< Tensor< 1, dim > > p_grads
Definition: polynomials_bdm.h:192
PolynomialsBDM::compute_n_pols
static unsigned int compute_n_pols(unsigned int degree)
Definition: polynomials_bdm.cc:430
PolynomialsBDM::p_fourth_derivatives
std::vector< Tensor< 4, dim > > p_fourth_derivatives
Definition: polynomials_bdm.h:207
StandardExceptions::ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
PolynomialsBDM::p_grad_grads
std::vector< Tensor< 2, dim > > p_grad_grads
Definition: polynomials_bdm.h:197
PolynomialsBDM::n_pols
unsigned int n_pols
Definition: polynomials_bdm.h:177
PolynomialsBDM::mutex
Threads::Mutex mutex
Definition: polynomials_bdm.h:182
PolynomialsBDM::polynomial_space
const PolynomialSpace< dim > polynomial_space
Definition: polynomials_bdm.h:166
StandardExceptions::ExcInternalError
static::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.11