Reference documentation for deal.II version 9.2.0
\(\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 - 2019 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 #include <deal.II/base/std_cxx14/memory.h>
22 
23 #include <iomanip>
24 #include <iostream>
25 
26 DEAL_II_NAMESPACE_OPEN
27 
28 
29 template <int dim>
30 PolynomialsBDM<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 
54 template <int dim>
55 void
56 PolynomialsBDM<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
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);
93  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
94  p_grad_grads.resize((grad_grads.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,
103  p_grads,
104  p_grad_grads,
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 
113  std::fill(grads.begin(), grads.end(), Tensor<2, dim>());
114  for (unsigned int i = 0; i < p_grads.size(); ++i)
115  for (unsigned int j = 0; j < dim; ++j)
116  grads[i + j * n_sub][j] = p_grads[i];
117 
118  std::fill(grad_grads.begin(), grad_grads.end(), Tensor<3, dim>());
119  for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
120  for (unsigned int j = 0; j < dim; ++j)
121  grad_grads[i + j * n_sub][j] = p_grad_grads[i];
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  }
145  if (grads.size() != 0)
146  {
147  grads[start][0][0] = monovali[0][1];
148  grads[start][0][1] = 0.;
149  grads[start][1][0] = -unit_point(1) * monovali[0][2];
150  grads[start][1][1] = -monovali[0][1];
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  }
156  if (grad_grads.size() != 0)
157  {
158  grad_grads[start][0][0][0] = monovali[0][2];
159  grad_grads[start][0][0][1] = 0.;
160  grad_grads[start][0][1][0] = 0.;
161  grad_grads[start][0][1][1] = 0.;
162  grad_grads[start][1][0][0] = -unit_point(1) * monovali[0][3];
163  grad_grads[start][1][0][1] = -monovali[0][2];
164  grad_grads[start][1][1][0] = -monovali[0][2];
165  grad_grads[start][1][1][1] = 0.;
166  grad_grads[start + 1][0][0][0] = 0;
167  grad_grads[start + 1][0][0][1] = monovali[1][2];
168  grad_grads[start + 1][0][1][0] = monovali[1][2];
169  grad_grads[start + 1][0][1][1] = unit_point(0) * monovali[1][3];
170  grad_grads[start + 1][1][0][0] = 0.;
171  grad_grads[start + 1][1][0][1] = 0.;
172  grad_grads[start + 1][1][1][0] = 0.;
173  grad_grads[start + 1][1][1][1] = -monovali[1][2];
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  }
222  if (grads.size() != 0)
223  {
224  grads[start][0][0] = monovali[1][1] * monovalk[2][0];
225  grads[start][0][1] =
226  unit_point(0) * monovali[1][2] * monovalk[2][0];
227  grads[start][0][2] =
228  unit_point(0) * monovali[1][1] * monovalk[2][1];
229  grads[start][1][0] = 0.;
230  grads[start][1][1] = -monovali[1][1] * monovalk[2][0];
231  grads[start][1][2] = -monovali[1][0] * monovalk[2][1];
232  grads[start][2][0] = 0.;
233  grads[start][2][1] = 0.;
234  grads[start][2][2] = 0.;
235 
236  grads[start + 1][1][1] = monovali[2][1] * monovalk[0][0];
237  grads[start + 1][1][2] =
238  unit_point(1) * monovali[2][2] * monovalk[0][0];
239  grads[start + 1][1][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];
249  grads[start + 2][2][0] =
250  unit_point(2) * monovali[0][2] * monovalk[1][0];
251  grads[start + 2][2][1] =
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  }
260  if (grad_grads.size() != 0)
261  {
262  grad_grads[start][0][0][0] = 0.;
263  grad_grads[start][0][0][1] = monovali[1][2] * monovalk[2][0];
264  grad_grads[start][0][0][2] = monovali[1][1] * monovalk[2][1];
265  grad_grads[start][0][1][0] = monovali[1][2] * monovalk[2][0];
266  grad_grads[start][0][1][1] =
267  unit_point(0) * monovali[1][3] * monovalk[2][0];
268  grad_grads[start][0][1][2] =
269  unit_point(0) * monovali[1][2] * monovalk[2][1];
270  grad_grads[start][0][2][0] = monovali[1][1] * monovalk[2][1];
271  grad_grads[start][0][2][1] =
272  unit_point(0) * monovali[1][2] * monovalk[2][1];
273  grad_grads[start][0][2][2] =
274  unit_point(0) * monovali[1][1] * monovalk[2][2];
275  grad_grads[start][1][0][0] = 0.;
276  grad_grads[start][1][0][1] = 0.;
277  grad_grads[start][1][0][2] = 0.;
278  grad_grads[start][1][1][0] = 0.;
279  grad_grads[start][1][1][1] = -monovali[1][2] * monovalk[2][0];
280  grad_grads[start][1][1][2] = -monovali[1][1] * monovalk[2][1];
281  grad_grads[start][1][2][0] = 0.;
282  grad_grads[start][1][2][1] = -monovali[1][1] * monovalk[2][1];
283  grad_grads[start][1][2][2] = -monovali[1][0] * monovalk[2][2];
284  grad_grads[start][2][0][0] = 0.;
285  grad_grads[start][2][0][1] = 0.;
286  grad_grads[start][2][0][2] = 0.;
287  grad_grads[start][2][1][0] = 0.;
288  grad_grads[start][2][1][1] = 0.;
289  grad_grads[start][2][1][2] = 0.;
290  grad_grads[start][2][2][0] = 0.;
291  grad_grads[start][2][2][1] = 0.;
292  grad_grads[start][2][2][2] = 0.;
293 
294  grad_grads[start + 1][0][0][0] = 0.;
295  grad_grads[start + 1][0][0][1] = 0.;
296  grad_grads[start + 1][0][0][2] = 0.;
297  grad_grads[start + 1][0][1][0] = 0.;
298  grad_grads[start + 1][0][1][1] = 0.;
299  grad_grads[start + 1][0][1][2] = 0.;
300  grad_grads[start + 1][0][2][0] = 0.;
301  grad_grads[start + 1][0][2][1] = 0.;
302  grad_grads[start + 1][0][2][2] = 0.;
303  grad_grads[start + 1][1][0][0] =
304  unit_point(1) * monovali[2][1] * monovalk[0][2];
305  grad_grads[start + 1][1][0][1] = monovali[2][1] * monovalk[0][1];
306  grad_grads[start + 1][1][0][2] =
307  unit_point(1) * monovali[2][2] * monovalk[0][1];
308  grad_grads[start + 1][1][1][0] = monovalk[0][1] * monovali[2][1];
309  grad_grads[start + 1][1][1][1] = 0.;
310  grad_grads[start + 1][1][1][2] = monovalk[0][0] * monovali[2][2];
311  grad_grads[start + 1][1][2][0] =
312  unit_point(1) * monovalk[0][1] * monovali[2][2];
313  grad_grads[start + 1][1][2][1] = monovalk[0][0] * monovali[2][2];
314  grad_grads[start + 1][1][2][2] =
315  unit_point(1) * monovalk[0][0] * monovali[2][3];
316  grad_grads[start + 1][2][0][0] = -monovalk[0][2] * monovali[2][0];
317  grad_grads[start + 1][2][0][1] = 0.;
318  grad_grads[start + 1][2][0][2] = -monovalk[0][1] * monovali[2][1];
319  grad_grads[start + 1][2][1][0] = 0.;
320  grad_grads[start + 1][2][1][1] = 0.;
321  grad_grads[start + 1][2][1][2] = 0.;
322  grad_grads[start + 1][2][2][0] = -monovalk[0][1] * monovali[2][1];
323  grad_grads[start + 1][2][2][1] = 0.;
324  grad_grads[start + 1][2][2][2] = -monovalk[0][0] * monovali[2][2];
325 
326  grad_grads[start + 2][0][0][0] = -monovali[0][2] * monovalk[1][0];
327  grad_grads[start + 2][0][0][1] = -monovali[0][1] * monovalk[1][1];
328  grad_grads[start + 2][0][0][2] = 0.;
329  grad_grads[start + 2][0][1][0] = -monovali[0][1] * monovalk[1][1];
330  grad_grads[start + 2][0][1][1] = -monovali[0][0] * monovalk[1][2];
331  grad_grads[start + 2][0][1][2] = 0.;
332  grad_grads[start + 2][0][2][0] = 0.;
333  grad_grads[start + 2][0][2][1] = 0.;
334  grad_grads[start + 2][0][2][2] = 0.;
335  grad_grads[start + 2][1][0][0] = 0.;
336  grad_grads[start + 2][1][0][1] = 0.;
337  grad_grads[start + 2][1][0][2] = 0.;
338  grad_grads[start + 2][1][1][0] = 0.;
339  grad_grads[start + 2][1][1][1] = 0.;
340  grad_grads[start + 2][1][1][2] = 0.;
341  grad_grads[start + 2][1][2][0] = 0.;
342  grad_grads[start + 2][1][2][1] = 0.;
343  grad_grads[start + 2][1][2][2] = 0.;
344  grad_grads[start + 2][2][0][0] =
345  unit_point(2) * monovali[0][3] * monovalk[1][0];
346  grad_grads[start + 2][2][0][1] =
347  unit_point(2) * monovali[0][2] * monovalk[1][1];
348  grad_grads[start + 2][2][0][2] = monovali[0][2] * monovalk[1][0];
349  grad_grads[start + 2][2][1][0] =
350  unit_point(2) * monovali[0][2] * monovalk[1][1];
351  grad_grads[start + 2][2][1][1] =
352  unit_point(2) * monovali[0][1] * monovalk[1][2];
353  grad_grads[start + 2][2][1][2] = monovali[0][1] * monovalk[1][1];
354  grad_grads[start + 2][2][2][0] = monovali[0][2] * monovalk[1][0];
355  grad_grads[start + 2][2][2][1] = monovali[0][1] * monovalk[1][1];
356  grad_grads[start + 2][2][2][2] = 0.;
357  }
358  }
359  Assert(start == this->n(), ExcInternalError());
360  }
361 }
362 
363 
364 /*
365 template <int dim>
366 void
367 PolynomialsBDM<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());
376  std::vector<Tensor<2,dim> > grads;
377  std::vector<Tensor<3,dim> > grad_grads;
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 
408  compute (p, values, grads, grad_grads);
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
415 square 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 
431 template <int dim>
432 unsigned int
433 PolynomialsBDM<dim>::n_polynomials(const unsigned int k)
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 
446 template <int dim>
447 std::unique_ptr<TensorPolynomialsBase<dim>>
448 PolynomialsBDM<dim>::clone() const
449 {
450  return std_cxx14::make_unique<PolynomialsBDM<dim>>(*this);
451 }
452 
453 
454 template class PolynomialsBDM<1>;
455 template class PolynomialsBDM<2>;
456 template class PolynomialsBDM<3>;
457 
458 
459 DEAL_II_NAMESPACE_CLOSE
Polynomials
Definition: polynomial.h:41
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
PolynomialsBDM::n_polynomials
static unsigned int n_polynomials(const unsigned int degree)
Definition: polynomials_bdm.cc:433
PolynomialsBDM::PolynomialsBDM
PolynomialsBDM(const unsigned int k)
Definition: polynomials_bdm.cc:30
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
quadrature_lib.h
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
polynomials_bdm.h
geometry_info.h
Tensor< 1, dim >
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
polynomial_space.h
TensorPolynomialsBase
Definition: tensor_polynomials_base.h:62
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
PolynomialsBDM
Definition: polynomials_bdm.h:101
Polynomials::Monomial
Definition: polynomial.h:299
PolynomialsBDM::monomials
std::vector< Polynomials::Polynomial< double > > monomials
Definition: polynomials_bdm.h:165
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
memory.h
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
PolynomialsBDM::clone
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
Definition: polynomials_bdm.cc:448
Generated by   doxygen 1.8.17