Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
polynomials_abf.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2018 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/polynomials_abf.h>
18 #include <deal.II/base/quadrature_lib.h>
19 
20 #include <iomanip>
21 #include <iostream>
22 
23 
24 DEAL_II_NAMESPACE_OPEN
25 
26 
27 
28 namespace
29 {
30  template <int dim>
31  std::vector<std::vector<Polynomials::Polynomial<double>>>
32  get_abf_polynomials(const unsigned int k)
33  {
34  std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
36 
37  if (k == 0)
38  for (unsigned int d = 1; d < dim; ++d)
40  else
41  for (unsigned int d = 1; d < dim; ++d)
43 
44  return pols;
45  }
46 } // namespace
47 
48 template <int dim>
50  : my_degree(k)
51  , polynomial_space(get_abf_polynomials<dim>(k))
52  , n_pols(compute_n_pols(k))
53 {
54  // check that the dimensions match. we only store one of the 'dim'
55  // anisotropic polynomials that make up the vector-valued space, so
56  // multiply by 'dim'
58 }
59 
60 
61 
62 template <int dim>
63 void
65  const Point<dim> & unit_point,
66  std::vector<Tensor<1, dim>> &values,
67  std::vector<Tensor<2, dim>> &grads,
68  std::vector<Tensor<3, dim>> &grad_grads,
69  std::vector<Tensor<4, dim>> &third_derivatives,
70  std::vector<Tensor<5, dim>> &fourth_derivatives) const
71 {
72  Assert(values.size() == n_pols || values.size() == 0,
73  ExcDimensionMismatch(values.size(), n_pols));
74  Assert(grads.size() == n_pols || grads.size() == 0,
75  ExcDimensionMismatch(grads.size(), n_pols));
76  Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
77  ExcDimensionMismatch(grad_grads.size(), n_pols));
78  Assert(third_derivatives.size() == n_pols || third_derivatives.size() == 0,
79  ExcDimensionMismatch(third_derivatives.size(), n_pols));
80  Assert(fourth_derivatives.size() == n_pols || fourth_derivatives.size() == 0,
81  ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
82 
83  const unsigned int n_sub = polynomial_space.n();
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  std::lock_guard<std::mutex> lock(mutex);
90 
91  p_values.resize((values.size() == 0) ? 0 : n_sub);
92  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
93  p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
94  p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
95  p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
96 
97  for (unsigned int d = 0; d < dim; ++d)
98  {
99  // First we copy the point. The
100  // polynomial space for
101  // component d consists of
102  // polynomials of degree k+1 in
103  // x_d and degree k in the
104  // other variables. in order to
105  // simplify this, we use the
106  // same AnisotropicPolynomial
107  // space and simply rotate the
108  // coordinates through all
109  // directions.
110  Point<dim> p;
111  for (unsigned int c = 0; c < dim; ++c)
112  p(c) = unit_point((c + d) % dim);
113 
114  polynomial_space.compute(p,
115  p_values,
116  p_grads,
117  p_grad_grads,
118  p_third_derivatives,
119  p_fourth_derivatives);
120 
121  for (unsigned int i = 0; i < p_values.size(); ++i)
122  values[i + d * n_sub][d] = p_values[i];
123 
124  for (unsigned int i = 0; i < p_grads.size(); ++i)
125  for (unsigned int d1 = 0; d1 < dim; ++d1)
126  grads[i + d * n_sub][d][(d1 + d) % dim] = p_grads[i][d1];
127 
128  for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
129  for (unsigned int d1 = 0; d1 < dim; ++d1)
130  for (unsigned int d2 = 0; d2 < dim; ++d2)
131  grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
132  p_grad_grads[i][d1][d2];
133 
134  for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
135  for (unsigned int d1 = 0; d1 < dim; ++d1)
136  for (unsigned int d2 = 0; d2 < dim; ++d2)
137  for (unsigned int d3 = 0; d3 < dim; ++d3)
138  third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
139  [(d2 + d) % dim][(d3 + d) % dim] =
140  p_third_derivatives[i][d1][d2][d3];
141 
142  for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
143  for (unsigned int d1 = 0; d1 < dim; ++d1)
144  for (unsigned int d2 = 0; d2 < dim; ++d2)
145  for (unsigned int d3 = 0; d3 < dim; ++d3)
146  for (unsigned int d4 = 0; d4 < dim; ++d4)
147  fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
148  [(d2 + d) % dim][(d3 + d) % dim]
149  [(d4 + d) % dim] =
150  p_fourth_derivatives[i][d1][d2][d3][d4];
151  }
152 }
153 
154 
155 template <int dim>
156 unsigned int
158 {
159  switch (dim)
160  {
161  case 1:
162  // in 1d, we simply have Q_{k+2}, which has dimension k+3
163  return k + 3;
164 
165  case 2:
166  // the polynomial space is Q_{k+2,k} \times Q_{k,k+2}, which has
167  // 2(k+3)(k+1) DoFs
168  return 2 * (k + 3) * (k + 1);
169 
170  case 3:
171  // the polynomial space is Q_{k+2,k,k} \times Q_{k,k+2,k} \times
172  // Q_{k,k,k+2}, which has 3(k+3)(k+1)(k+1) DoFs
173  return 3 * (k + 3) * (k + 1) * (k + 1);
174 
175  default:
176  Assert(false, ExcNotImplemented());
177  }
178 
179  return 0;
180 }
181 
182 
183 template class PolynomialsABF<1>;
184 template class PolynomialsABF<2>;
185 template class PolynomialsABF<3>;
186 
187 
188 DEAL_II_NAMESPACE_CLOSE
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:866
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:800
#define Assert(cond, exc)
Definition: exceptions.h:1407
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)
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
static unsigned int compute_n_pols(unsigned int degree)
static ::ExceptionBase & ExcNotImplemented()
const AnisotropicPolynomials< dim > polynomial_space
PolynomialsABF(const unsigned int k)
static ::ExceptionBase & ExcInternalError()