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_raviart_thomas.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 
20 
21 #include <iomanip>
22 #include <iostream>
23 #include <memory>
24 
25 // TODO[WB]: This class is not thread-safe: it uses mutable member variables
26 // that contain temporary state. this is not what one would want when one uses a
27 // finite element object in a number of different contexts on different threads:
28 // finite element objects should be stateless
29 // TODO:[GK] This can be achieved by writing a function in Polynomial space
30 // which does the rotated fill performed below and writes the data into the
31 // right data structures. The same function would be used by Nedelec
32 // polynomials.
33 
35 
36 
37 template <int dim>
39  : TensorPolynomialsBase<dim>(k, n_polynomials(k))
40  , polynomial_space(create_polynomials(k))
41 {}
42 
43 
44 
45 template <int dim>
46 std::vector<std::vector<Polynomials::Polynomial<double>>>
48 {
49  std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
51  if (k == 0)
52  for (unsigned int d = 1; d < dim; ++d)
54  else
55  for (unsigned int d = 1; d < dim; ++d)
57 
58  return pols;
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  // have a few scratch
85  // arrays. because we don't want to
86  // re-allocate them every time this
87  // function is called, we make them
88  // static. however, in return we
89  // have to ensure that the calls to
90  // the use of these variables is
91  // locked with a mutex. if the
92  // mutex is removed, several tests
93  // (notably
94  // deal.II/create_mass_matrix_05)
95  // will start to produce random
96  // results in multithread mode
97  static std::mutex mutex;
98  std::lock_guard<std::mutex> lock(mutex);
99 
100  static std::vector<double> p_values;
101  static std::vector<Tensor<1, dim>> p_grads;
102  static std::vector<Tensor<2, dim>> p_grad_grads;
103  static std::vector<Tensor<3, dim>> p_third_derivatives;
104  static std::vector<Tensor<4, dim>> p_fourth_derivatives;
105 
106  const unsigned int n_sub = polynomial_space.n();
107  p_values.resize((values.size() == 0) ? 0 : n_sub);
108  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
109  p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
110  p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
111  p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
112 
113  for (unsigned int d = 0; d < dim; ++d)
114  {
115  // First we copy the point. The
116  // polynomial space for
117  // component d consists of
118  // polynomials of degree k+1 in
119  // x_d and degree k in the
120  // other variables. in order to
121  // simplify this, we use the
122  // same AnisotropicPolynomial
123  // space and simply rotate the
124  // coordinates through all
125  // directions.
126  Point<dim> p;
127  for (unsigned int c = 0; c < dim; ++c)
128  p(c) = unit_point((c + d) % dim);
129 
130  polynomial_space.evaluate(p,
131  p_values,
132  p_grads,
133  p_grad_grads,
134  p_third_derivatives,
135  p_fourth_derivatives);
136 
137  for (unsigned int i = 0; i < p_values.size(); ++i)
138  values[i + d * n_sub][d] = p_values[i];
139 
140  for (unsigned int i = 0; i < p_grads.size(); ++i)
141  for (unsigned int d1 = 0; d1 < dim; ++d1)
142  grads[i + d * n_sub][d][(d1 + d) % dim] = p_grads[i][d1];
143 
144  for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
145  for (unsigned int d1 = 0; d1 < dim; ++d1)
146  for (unsigned int d2 = 0; d2 < dim; ++d2)
147  grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
148  p_grad_grads[i][d1][d2];
149 
150  for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
151  for (unsigned int d1 = 0; d1 < dim; ++d1)
152  for (unsigned int d2 = 0; d2 < dim; ++d2)
153  for (unsigned int d3 = 0; d3 < dim; ++d3)
154  third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
155  [(d2 + d) % dim][(d3 + d) % dim] =
156  p_third_derivatives[i][d1][d2][d3];
157 
158  for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
159  for (unsigned int d1 = 0; d1 < dim; ++d1)
160  for (unsigned int d2 = 0; d2 < dim; ++d2)
161  for (unsigned int d3 = 0; d3 < dim; ++d3)
162  for (unsigned int d4 = 0; d4 < dim; ++d4)
163  fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
164  [(d2 + d) % dim][(d3 + d) % dim]
165  [(d4 + d) % dim] =
166  p_fourth_derivatives[i][d1][d2][d3][d4];
167  }
168 }
169 
170 
171 template <int dim>
172 unsigned int
174 {
175  if (dim == 1)
176  return k + 1;
177  if (dim == 2)
178  return 2 * (k + 1) * (k + 2);
179  if (dim == 3)
180  return 3 * (k + 1) * (k + 1) * (k + 2);
181 
182  Assert(false, ExcNotImplemented());
183  return 0;
184 }
185 
186 
187 template <int dim>
188 std::unique_ptr<TensorPolynomialsBase<dim>>
190 {
191  return std::make_unique<PolynomialsRaviartThomas<dim>>(*this);
192 }
193 
194 
195 template class PolynomialsRaviartThomas<1>;
196 template class PolynomialsRaviartThomas<2>;
197 template class PolynomialsRaviartThomas<3>;
198 
199 
static unsigned int n_polynomials(const unsigned int degree)
PolynomialsRaviartThomas(const unsigned int k)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:744
const AnisotropicPolynomials< dim > polynomial_space
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:678
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
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
#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
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials(const unsigned int k)
static ::ExceptionBase & ExcNotImplemented()