Reference documentation for deal.II version 9.4.1
Searching...
No Matches
polynomials_raviart_thomas.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2021 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
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
37template <int dim>
39 : TensorPolynomialsBase<dim>(k, n_polynomials(k))
40 , polynomial_space(create_polynomials(k))
41{}
42
43
44
45template <int dim>
47 const PolynomialsRaviartThomas &other)
48 : TensorPolynomialsBase<dim>(other)
49 , polynomial_space(other.polynomial_space)
50// no need to copy the scratch data, or the mutex
51{}
52
53
54
55template <int dim>
56std::vector<std::vector<Polynomials::Polynomial<double>>>
58{
59 // Create a vector of polynomial spaces where the first element
60 // has degree k+1 and the rest has degree k. This corresponds to
61 // the space of single-variable polynomials from which we will create the
62 // space for the first component of the RT element by way of tensor
63 // product.
64 //
65 // The other components of the RT space can be created by rotating
66 // this vector of single-variable polynomials.
67 //
68 std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
69 if (k == 0)
70 {
71 // We need to treat the case k=0 differently because there,
72 // the polynomial in x has degree 1 and so has node points
73 // equal to the end points of the interval [0,1] (i.e., it
74 // is a "Lagrange" polynomial), whereas the y- and z-polynomials
75 // only have the interval midpoint as a node (i.e., they are
76 // a "Legendre" polynomial).
78 for (unsigned int d = 1; d < dim; ++d)
80 }
81 else
82 {
83 pols[0] =
85 for (unsigned int d = 1; d < dim; ++d)
87 }
88
89 return pols;
90}
91
92
93
94template <int dim>
95void
97 const Point<dim> & unit_point,
98 std::vector<Tensor<1, dim>> &values,
101 std::vector<Tensor<4, dim>> &third_derivatives,
102 std::vector<Tensor<5, dim>> &fourth_derivatives) const
103{
104 Assert(values.size() == this->n() || values.size() == 0,
105 ExcDimensionMismatch(values.size(), this->n()));
110 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
111 ExcDimensionMismatch(third_derivatives.size(), this->n()));
112 Assert(fourth_derivatives.size() == this->n() ||
113 fourth_derivatives.size() == 0,
114 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
115
116 // In the following, we will access the scratch arrays declared as 'mutable'
117 // in the class. In order to do so from this function, we have to make sure
118 // that we guard this access by a mutex so that the invocation of this 'const'
120 std::lock_guard<std::mutex> lock(scratch_arrays_mutex);
121
122 const unsigned int n_sub = polynomial_space.n();
123 scratch_values.resize((values.size() == 0) ? 0 : n_sub);
126 scratch_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
127 scratch_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 :
128 n_sub);
129
130 for (unsigned int d = 0; d < dim; ++d)
131 {
132 // First we copy the point. The
133 // polynomial space for
134 // component d consists of
135 // polynomials of degree k+1 in
136 // x_d and degree k in the
137 // other variables. in order to
138 // simplify this, we use the
139 // same AnisotropicPolynomial
140 // space and simply rotate the
141 // coordinates through all
142 // directions.
143 Point<dim> p;
144 for (unsigned int c = 0; c < dim; ++c)
145 p(c) = unit_point((c + d) % dim);
146
147 polynomial_space.evaluate(p,
148 scratch_values,
151 scratch_third_derivatives,
152 scratch_fourth_derivatives);
153
154 for (unsigned int i = 0; i < scratch_values.size(); ++i)
155 values[i + d * n_sub][d] = scratch_values[i];
156
157 for (unsigned int i = 0; i < scratch_grads.size(); ++i)
158 for (unsigned int d1 = 0; d1 < dim; ++d1)
159 grads[i + d * n_sub][d][(d1 + d) % dim] = scratch_grads[i][d1];
160
162 for (unsigned int d1 = 0; d1 < dim; ++d1)
163 for (unsigned int d2 = 0; d2 < dim; ++d2)
164 grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
166
167 for (unsigned int i = 0; i < scratch_third_derivatives.size(); ++i)
168 for (unsigned int d1 = 0; d1 < dim; ++d1)
169 for (unsigned int d2 = 0; d2 < dim; ++d2)
170 for (unsigned int d3 = 0; d3 < dim; ++d3)
171 third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
172 [(d2 + d) % dim][(d3 + d) % dim] =
173 scratch_third_derivatives[i][d1][d2][d3];
174
175 for (unsigned int i = 0; i < scratch_fourth_derivatives.size(); ++i)
176 for (unsigned int d1 = 0; d1 < dim; ++d1)
177 for (unsigned int d2 = 0; d2 < dim; ++d2)
178 for (unsigned int d3 = 0; d3 < dim; ++d3)
179 for (unsigned int d4 = 0; d4 < dim; ++d4)
180 fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
181 [(d2 + d) % dim][(d3 + d) % dim]
182 [(d4 + d) % dim] =
183 scratch_fourth_derivatives[i][d1][d2][d3]
184 [d4];
185 }
186}
187
188
189template <int dim>
190unsigned int
192{
193 if (dim == 1)
194 return k + 1;
195 if (dim == 2)
196 return 2 * (k + 1) * (k + 2);
197 if (dim == 3)
198 return 3 * (k + 1) * (k + 1) * (k + 2);
199
200 Assert(false, ExcNotImplemented());
201 return 0;
202}
203
204
205template <int dim>
206std::unique_ptr<TensorPolynomialsBase<dim>>
208{
209 return std::make_unique<PolynomialsRaviartThomas<dim>>(*this);
210}
211
212
213template class PolynomialsRaviartThomas<1>;
214template class PolynomialsRaviartThomas<2>;
215template class PolynomialsRaviartThomas<3>;
216
217
Definition: point.h:111
static unsigned int n_polynomials(const unsigned int degree)
PolynomialsRaviartThomas(const unsigned int k)
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials(const unsigned int k)
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
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:679
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:745
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)