Reference documentation for deal.II version 9.5.0
Searching...
No Matches
polynomials_bernardi_raugel.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2023 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
18
19#include <memory>
20
22
23
24template <int dim>
26 : TensorPolynomialsBase<dim>(k + 1, n_polynomials(k))
27 , polynomial_space_Q(create_polynomials_Q())
28 , polynomial_space_bubble(create_polynomials_bubble())
29{}
30
31
32template <int dim>
33std::vector<std::vector<Polynomials::Polynomial<double>>>
35{
36 std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
37 std::vector<Polynomials::Polynomial<double>> bubble_shapes;
38 bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
39 bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
40 bubble_shapes.push_back(Polynomials::LagrangeEquidistant(2, 1));
41
42 for (unsigned int d = 0; d < dim; ++d)
43 pols.push_back(bubble_shapes);
44 // In 2d, the only q_ij polynomials we will use are 31,32,13,23
45 // where ij corresponds to index (i-1)+3*(j-1) (2,5,6,7)
46
47 // In 3d, the only q_ijk polynomials we will use are 331,332,313,323,133,233
48 // where ijk corresponds to index (i-1)+3*(j-1)+9*(k-1) (8,17,20,23,24,25)
49 return pols;
50}
51
52
53
54template <int dim>
55std::vector<std::vector<Polynomials::Polynomial<double>>>
57{
58 std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
59 std::vector<Polynomials::Polynomial<double>> Q_shapes;
60 Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
61 Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
62 for (unsigned int d = 0; d < dim; ++d)
63 pols.push_back(Q_shapes);
64
65 return pols;
66}
67
68
69template <int dim>
70void
72 const Point<dim> & unit_point,
73 std::vector<Tensor<1, dim>> &values,
76 std::vector<Tensor<4, dim>> &third_derivatives,
77 std::vector<Tensor<5, dim>> &fourth_derivatives) const
78{
79 Assert(values.size() == this->n() || values.size() == 0,
80 ExcDimensionMismatch(values.size(), this->n()));
85 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
86 ExcDimensionMismatch(third_derivatives.size(), this->n()));
87 Assert(fourth_derivatives.size() == this->n() ||
88 fourth_derivatives.size() == 0,
89 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
90
91 std::vector<double> Q_values;
94 std::vector<Tensor<3, dim>> Q_third_derivatives;
95 std::vector<Tensor<4, dim>> Q_fourth_derivatives;
96 std::vector<double> bubble_values;
99 std::vector<Tensor<3, dim>> bubble_third_derivatives;
100 std::vector<Tensor<4, dim>> bubble_fourth_derivatives;
101
102 constexpr int n_bubbles =
103 Utilities::pow(3, dim); // size for create_polynomials_bubble
104 constexpr int n_q = 1 << dim; // size for create_polynomials_q
105
106 // don't resize if the provided vector has 0 length
107 Q_values.resize((values.size() == 0) ? 0 : n_q);
110 Q_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_q);
111 Q_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_q);
112 bubble_values.resize((values.size() == 0) ? 0 : n_bubbles);
115 bubble_third_derivatives.resize((third_derivatives.size() == 0) ? 0 :
116 n_bubbles);
117 bubble_fourth_derivatives.resize(
118 (fourth_derivatives.size() == 0) ? 0 : n_bubbles);
119
120 // 1 normal vector per face, ordering consistent with GeometryInfo
121 // Normal vectors point in the +x, +y, and +z directions for
122 // consistent orientation across edges
123 std::vector<Tensor<1, dim>> normals;
124 for (const unsigned int i : GeometryInfo<dim>::face_indices())
125 {
126 Tensor<1, dim> normal;
127 normal[i / 2] = 1;
128 normals.push_back(normal);
129 }
130
131 // dim standard basis vectors for R^dim, usual ordering
132 std::vector<Tensor<1, dim>> units;
133 for (unsigned int i = 0; i < dim; ++i)
134 {
135 Tensor<1, dim> unit;
136 unit[i] = 1;
137 units.push_back(unit);
138 }
139
140 // set indices for the anisotropic polynomials to find
141 // them after polynomial_space_bubble.evaluate is called
142 std::vector<int> aniso_indices;
143 if (dim == 2)
144 {
145 aniso_indices.push_back(6);
146 aniso_indices.push_back(7);
147 aniso_indices.push_back(2);
148 aniso_indices.push_back(5);
149 }
150 else if (dim == 3)
151 {
152 aniso_indices.push_back(24);
153 aniso_indices.push_back(25);
154 aniso_indices.push_back(20);
155 aniso_indices.push_back(23);
156 aniso_indices.push_back(8);
157 aniso_indices.push_back(17);
158 }
159
160 polynomial_space_bubble.evaluate(unit_point,
161 bubble_values,
164 bubble_third_derivatives,
165 bubble_fourth_derivatives);
166 polynomial_space_Q.evaluate(unit_point,
167 Q_values,
170 Q_third_derivatives,
171 Q_fourth_derivatives);
172
173 // first dim*vertices_per_cell functions are Q_1^2 functions
174 for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
175 {
176 if (values.size() != 0)
177 {
178 values[i] = units[i % dim] * Q_values[i / dim];
179 }
181 {
183 }
185 {
187 }
188 if (third_derivatives.size() != 0)
189 {
190 third_derivatives[i] =
191 outer_product(units[i % dim], Q_third_derivatives[i / dim]);
192 }
193 if (fourth_derivatives.size() != 0)
194 {
195 fourth_derivatives[i] =
196 outer_product(units[i % dim], Q_fourth_derivatives[i / dim]);
197 }
198 }
199
200 // last faces_per_cell functions are bubble functions
201 for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
202 i < dim * GeometryInfo<dim>::vertices_per_cell +
204 ++i)
205 {
206 unsigned int j =
207 i -
208 dim *
209 GeometryInfo<dim>::vertices_per_cell; // ranges 0 to faces_per_cell-1
210 if (values.size() != 0)
211 {
212 values[i] = normals[j] * bubble_values[aniso_indices[j]];
213 }
215 {
217 }
219 {
222 }
223 if (third_derivatives.size() != 0)
224 {
225 third_derivatives[i] =
226 outer_product(normals[j],
227 bubble_third_derivatives[aniso_indices[j]]);
228 }
229 if (fourth_derivatives.size() != 0)
230 {
231 fourth_derivatives[i] =
232 outer_product(normals[j],
233 bubble_fourth_derivatives[aniso_indices[j]]);
234 }
235 }
236}
237
238template <int dim>
239unsigned int
241{
242 (void)k;
243 Assert(k == 1, ExcNotImplemented());
244 if (dim == 2 || dim == 3)
247 // 2*4+4=12 polynomials in 2d and 3*8+6=30 polynomials in 3d
248
249 Assert(false, ExcNotImplemented());
250 return 0;
251}
252
253
254template <int dim>
255std::unique_ptr<TensorPolynomialsBase<dim>>
257{
258 return std::make_unique<PolynomialsBernardiRaugel<dim>>(*this);
259}
260
261template class PolynomialsBernardiRaugel<1>; // to prevent errors
262template class PolynomialsBernardiRaugel<2>;
263template class PolynomialsBernardiRaugel<3>;
264
265
Definition point.h:112
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials_bubble()
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
PolynomialsBernardiRaugel(const unsigned int k)
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials_Q()
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
static unsigned int n_polynomials(const unsigned int k)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)