Reference documentation for deal.II version Git 932f7faded 2020-11-28 20:02:43 +0100
\(\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_abf.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 
19 
20 #include <iomanip>
21 #include <iostream>
22 #include <memory>
23 
24 
26 
27 
28 
29 namespace
30 {
31  template <int dim>
32  std::vector<std::vector<Polynomials::Polynomial<double>>>
33  get_abf_polynomials(const unsigned int k)
34  {
35  std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
37 
38  if (k == 0)
39  for (unsigned int d = 1; d < dim; ++d)
41  else
42  for (unsigned int d = 1; d < dim; ++d)
44 
45  return pols;
46  }
47 } // namespace
48 
49 template <int dim>
51  : TensorPolynomialsBase<dim>(k, n_polynomials(k))
52  , polynomial_space(get_abf_polynomials<dim>(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() == this->n() || values.size() == 0,
73  ExcDimensionMismatch(values.size(), this->n()));
74  Assert(grads.size() == this->n() || grads.size() == 0,
75  ExcDimensionMismatch(grads.size(), this->n()));
76  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
77  ExcDimensionMismatch(grad_grads.size(), this->n()));
78  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
79  ExcDimensionMismatch(third_derivatives.size(), this->n()));
80  Assert(fourth_derivatives.size() == this->n() ||
81  fourth_derivatives.size() == 0,
82  ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
83 
84  const unsigned int n_sub = polynomial_space.n();
85  // guard access to the scratch
86  // arrays in the following block
87  // using a mutex to make sure they
88  // are not used by multiple threads
89  // at once
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  p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
96  p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
97 
98  for (unsigned int d = 0; d < dim; ++d)
99  {
100  // First we copy the point. The
101  // polynomial space for
102  // component d consists of
103  // polynomials of degree k+1 in
104  // x_d and degree k in the
105  // other variables. in order to
106  // simplify this, we use the
107  // same AnisotropicPolynomial
108  // space and simply rotate the
109  // coordinates through all
110  // directions.
111  Point<dim> p;
112  for (unsigned int c = 0; c < dim; ++c)
113  p(c) = unit_point((c + d) % dim);
114 
115  polynomial_space.evaluate(p,
116  p_values,
117  p_grads,
118  p_grad_grads,
121 
122  for (unsigned int i = 0; i < p_values.size(); ++i)
123  values[i + d * n_sub][d] = p_values[i];
124 
125  for (unsigned int i = 0; i < p_grads.size(); ++i)
126  for (unsigned int d1 = 0; d1 < dim; ++d1)
127  grads[i + d * n_sub][d][(d1 + d) % dim] = p_grads[i][d1];
128 
129  for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
130  for (unsigned int d1 = 0; d1 < dim; ++d1)
131  for (unsigned int d2 = 0; d2 < dim; ++d2)
132  grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
133  p_grad_grads[i][d1][d2];
134 
135  for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
136  for (unsigned int d1 = 0; d1 < dim; ++d1)
137  for (unsigned int d2 = 0; d2 < dim; ++d2)
138  for (unsigned int d3 = 0; d3 < dim; ++d3)
139  third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
140  [(d2 + d) % dim][(d3 + d) % dim] =
141  p_third_derivatives[i][d1][d2][d3];
142 
143  for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
144  for (unsigned int d1 = 0; d1 < dim; ++d1)
145  for (unsigned int d2 = 0; d2 < dim; ++d2)
146  for (unsigned int d3 = 0; d3 < dim; ++d3)
147  for (unsigned int d4 = 0; d4 < dim; ++d4)
148  fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
149  [(d2 + d) % dim][(d3 + d) % dim]
150  [(d4 + d) % dim] =
151  p_fourth_derivatives[i][d1][d2][d3][d4];
152  }
153 }
154 
155 
156 template <int dim>
157 unsigned int
159 {
160  switch (dim)
161  {
162  case 1:
163  // in 1d, we simply have Q_{k+2}, which has dimension k+3
164  return k + 3;
165 
166  case 2:
167  // the polynomial space is Q_{k+2,k} \times Q_{k,k+2}, which has
168  // 2(k+3)(k+1) DoFs
169  return 2 * (k + 3) * (k + 1);
170 
171  case 3:
172  // the polynomial space is Q_{k+2,k,k} \times Q_{k,k+2,k} \times
173  // Q_{k,k,k+2}, which has 3(k+3)(k+1)(k+1) DoFs
174  return 3 * (k + 3) * (k + 1) * (k + 1);
175 
176  default:
177  Assert(false, ExcNotImplemented());
178  }
179 
180  return 0;
181 }
182 
183 
184 template <int dim>
185 std::unique_ptr<TensorPolynomialsBase<dim>>
187 {
188  return std::make_unique<PolynomialsABF<dim>>(*this);
189 }
190 
191 
192 template class PolynomialsABF<1>;
193 template class PolynomialsABF<2>;
194 template class PolynomialsABF<3>;
195 
196 
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:744
std::vector< double > p_values
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:678
std::vector< Tensor< 3, dim > > p_third_derivatives
std::vector< Tensor< 1, dim > > p_grads
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
std::vector< Tensor< 4, dim > > p_fourth_derivatives
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
static unsigned int n_polynomials(const unsigned int degree)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
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
static ::ExceptionBase & ExcNotImplemented()
const AnisotropicPolynomials< dim > polynomial_space
std::vector< Tensor< 2, dim > > p_grad_grads
PolynomialsABF(const unsigned int k)
Threads::Mutex mutex
static ::ExceptionBase & ExcInternalError()