Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
polynomials_bernardi_raugel.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17
18#include <memory>
19
21
22
23template <int dim>
25 : TensorPolynomialsBase<dim>(k + 1, n_polynomials(k))
26 , polynomial_space_Q(create_polynomials_Q())
27 , polynomial_space_bubble(create_polynomials_bubble())
28{}
29
30
31template <int dim>
32std::vector<std::vector<Polynomials::Polynomial<double>>>
34{
35 std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
36 std::vector<Polynomials::Polynomial<double>> bubble_shapes;
37 bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
38 bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
39 bubble_shapes.push_back(Polynomials::LagrangeEquidistant(2, 1));
40
41 for (unsigned int d = 0; d < dim; ++d)
42 pols.push_back(bubble_shapes);
43 // In 2d, the only q_ij polynomials we will use are 31,32,13,23
44 // where ij corresponds to index (i-1)+3*(j-1) (2,5,6,7)
45
46 // In 3d, the only q_ijk polynomials we will use are 331,332,313,323,133,233
47 // where ijk corresponds to index (i-1)+3*(j-1)+9*(k-1) (8,17,20,23,24,25)
48 return pols;
49}
50
51
52
53template <int dim>
54std::vector<std::vector<Polynomials::Polynomial<double>>>
56{
57 std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
58 std::vector<Polynomials::Polynomial<double>> Q_shapes;
59 Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
60 Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
61 for (unsigned int d = 0; d < dim; ++d)
62 pols.push_back(Q_shapes);
63
64 return pols;
65}
66
67
68template <int dim>
69void
71 const Point<dim> &unit_point,
72 std::vector<Tensor<1, dim>> &values,
73 std::vector<Tensor<2, dim>> &grads,
74 std::vector<Tensor<3, dim>> &grad_grads,
75 std::vector<Tensor<4, dim>> &third_derivatives,
76 std::vector<Tensor<5, dim>> &fourth_derivatives) const
77{
78 Assert(values.size() == this->n() || values.empty(),
79 ExcDimensionMismatch(values.size(), this->n()));
80 Assert(grads.size() == this->n() || grads.empty(),
81 ExcDimensionMismatch(grads.size(), this->n()));
82 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
83 ExcDimensionMismatch(grad_grads.size(), this->n()));
84 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
85 ExcDimensionMismatch(third_derivatives.size(), this->n()));
86 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
87 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
88
89 std::vector<double> Q_values;
90 std::vector<Tensor<1, dim>> Q_grads;
91 std::vector<Tensor<2, dim>> Q_grad_grads;
92 std::vector<Tensor<3, dim>> Q_third_derivatives;
93 std::vector<Tensor<4, dim>> Q_fourth_derivatives;
94 std::vector<double> bubble_values;
95 std::vector<Tensor<1, dim>> bubble_grads;
96 std::vector<Tensor<2, dim>> bubble_grad_grads;
97 std::vector<Tensor<3, dim>> bubble_third_derivatives;
98 std::vector<Tensor<4, dim>> bubble_fourth_derivatives;
99
100 constexpr int n_bubbles =
101 Utilities::pow(3, dim); // size for create_polynomials_bubble
102 constexpr int n_q = 1 << dim; // size for create_polynomials_q
103
104 // don't resize if the provided vector has 0 length
105 Q_values.resize((values.empty()) ? 0 : n_q);
106 Q_grads.resize((grads.empty()) ? 0 : n_q);
107 Q_grad_grads.resize((grad_grads.empty()) ? 0 : n_q);
108 Q_third_derivatives.resize((third_derivatives.empty()) ? 0 : n_q);
109 Q_fourth_derivatives.resize((fourth_derivatives.empty()) ? 0 : n_q);
110 bubble_values.resize((values.empty()) ? 0 : n_bubbles);
111 bubble_grads.resize((grads.empty()) ? 0 : n_bubbles);
112 bubble_grad_grads.resize((grad_grads.empty()) ? 0 : n_bubbles);
113 bubble_third_derivatives.resize((third_derivatives.empty()) ? 0 : n_bubbles);
114 bubble_fourth_derivatives.resize((fourth_derivatives.empty()) ? 0 :
115 n_bubbles);
116
117 // 1 normal vector per face, ordering consistent with GeometryInfo
118 // Normal vectors point in the +x, +y, and +z directions for
119 // consistent orientation across edges
120 std::vector<Tensor<1, dim>> normals;
121 for (const unsigned int i : GeometryInfo<dim>::face_indices())
122 {
123 Tensor<1, dim> normal;
124 normal[i / 2] = 1;
125 normals.push_back(normal);
126 }
127
128 // dim standard basis vectors for R^dim, usual ordering
129 std::vector<Tensor<1, dim>> units;
130 for (unsigned int i = 0; i < dim; ++i)
131 {
132 Tensor<1, dim> unit;
133 unit[i] = 1;
134 units.push_back(unit);
135 }
136
137 // set indices for the anisotropic polynomials to find
138 // them after polynomial_space_bubble.evaluate is called
139 std::vector<int> aniso_indices;
140 if (dim == 2)
141 {
142 aniso_indices.push_back(6);
143 aniso_indices.push_back(7);
144 aniso_indices.push_back(2);
145 aniso_indices.push_back(5);
146 }
147 else if (dim == 3)
148 {
149 aniso_indices.push_back(24);
150 aniso_indices.push_back(25);
151 aniso_indices.push_back(20);
152 aniso_indices.push_back(23);
153 aniso_indices.push_back(8);
154 aniso_indices.push_back(17);
155 }
156
157 polynomial_space_bubble.evaluate(unit_point,
158 bubble_values,
159 bubble_grads,
160 bubble_grad_grads,
161 bubble_third_derivatives,
162 bubble_fourth_derivatives);
163 polynomial_space_Q.evaluate(unit_point,
164 Q_values,
165 Q_grads,
166 Q_grad_grads,
167 Q_third_derivatives,
168 Q_fourth_derivatives);
169
170 // first dim*vertices_per_cell functions are Q_1^2 functions
171 for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
172 {
173 if (values.size() != 0)
174 {
175 values[i] = units[i % dim] * Q_values[i / dim];
176 }
177 if (grads.size() != 0)
178 {
179 grads[i] = outer_product(units[i % dim], Q_grads[i / dim]);
180 }
181 if (grad_grads.size() != 0)
182 {
183 grad_grads[i] = outer_product(units[i % dim], Q_grad_grads[i / dim]);
184 }
185 if (third_derivatives.size() != 0)
186 {
187 third_derivatives[i] =
188 outer_product(units[i % dim], Q_third_derivatives[i / dim]);
189 }
190 if (fourth_derivatives.size() != 0)
191 {
192 fourth_derivatives[i] =
193 outer_product(units[i % dim], Q_fourth_derivatives[i / dim]);
194 }
195 }
196
197 // last faces_per_cell functions are bubble functions
198 for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
199 i < dim * GeometryInfo<dim>::vertices_per_cell +
201 ++i)
202 {
203 unsigned int j =
204 i -
205 dim *
206 GeometryInfo<dim>::vertices_per_cell; // ranges 0 to faces_per_cell-1
207 if (values.size() != 0)
208 {
209 values[i] = normals[j] * bubble_values[aniso_indices[j]];
210 }
211 if (grads.size() != 0)
212 {
213 grads[i] = outer_product(normals[j], bubble_grads[aniso_indices[j]]);
214 }
215 if (grad_grads.size() != 0)
216 {
217 grad_grads[i] =
218 outer_product(normals[j], bubble_grad_grads[aniso_indices[j]]);
219 }
220 if (third_derivatives.size() != 0)
221 {
222 third_derivatives[i] =
223 outer_product(normals[j],
224 bubble_third_derivatives[aniso_indices[j]]);
225 }
226 if (fourth_derivatives.size() != 0)
227 {
228 fourth_derivatives[i] =
229 outer_product(normals[j],
230 bubble_fourth_derivatives[aniso_indices[j]]);
231 }
232 }
233}
234
235template <int dim>
236unsigned int
238{
239 (void)k;
240 Assert(k == 1, ExcNotImplemented());
241 if (dim == 2 || dim == 3)
244 // 2*4+4=12 polynomials in 2d and 3*8+6=30 polynomials in 3d
245
247 return 0;
248}
249
250
251template <int dim>
252std::unique_ptr<TensorPolynomialsBase<dim>>
254{
255 return std::make_unique<PolynomialsBernardiRaugel<dim>>(*this);
256}
257
258template class PolynomialsBernardiRaugel<1>; // to prevent errors
259template class PolynomialsBernardiRaugel<2>;
260template class PolynomialsBernardiRaugel<3>;
261
262
Definition point.h:111
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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NOT_IMPLEMENTED()
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
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)