Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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
37template <int dim>
39 : TensorPolynomialsBase<dim>(k, n_polynomials(k))
40 , polynomial_space(create_polynomials(k))
41{}
42
43
44
45template <int dim>
46std::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
62template <int dim>
63void
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
171template <int dim>
172unsigned 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
187template <int dim>
188std::unique_ptr<TensorPolynomialsBase<dim>>
190{
191 return std::make_unique<PolynomialsRaviartThomas<dim>>(*this);
192}
193
194
195template class PolynomialsRaviartThomas<1>;
196template class PolynomialsRaviartThomas<2>;
197template class PolynomialsRaviartThomas<3>;
198
199
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:678
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:744
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
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)