Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3124-g7c413cc06e 2025-04-23 12:20:01+00:00
\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 Concepts
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 pols.reserve(dim);
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 pols.reserve(dim);
63 for (unsigned int d = 0; d < dim; ++d)
64 pols.push_back(Q_shapes);
65
66 return pols;
67}
68
69
70template <int dim>
71void
73 const Point<dim> &unit_point,
74 std::vector<Tensor<1, dim>> &values,
75 std::vector<Tensor<2, dim>> &grads,
76 std::vector<Tensor<3, dim>> &grad_grads,
77 std::vector<Tensor<4, dim>> &third_derivatives,
78 std::vector<Tensor<5, dim>> &fourth_derivatives) const
79{
80 Assert(values.size() == this->n() || values.empty(),
81 ExcDimensionMismatch(values.size(), this->n()));
82 Assert(grads.size() == this->n() || grads.empty(),
83 ExcDimensionMismatch(grads.size(), this->n()));
84 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
85 ExcDimensionMismatch(grad_grads.size(), this->n()));
86 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
87 ExcDimensionMismatch(third_derivatives.size(), this->n()));
88 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
89 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
90
91 std::vector<double> Q_values;
92 std::vector<Tensor<1, dim>> Q_grads;
93 std::vector<Tensor<2, dim>> Q_grad_grads;
94 std::vector<Tensor<3, dim>> Q_third_derivatives;
95 std::vector<Tensor<4, dim>> Q_fourth_derivatives;
96 std::vector<double> bubble_values;
97 std::vector<Tensor<1, dim>> bubble_grads;
98 std::vector<Tensor<2, dim>> bubble_grad_grads;
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.empty()) ? 0 : n_q);
108 Q_grads.resize((grads.empty()) ? 0 : n_q);
109 Q_grad_grads.resize((grad_grads.empty()) ? 0 : n_q);
110 Q_third_derivatives.resize((third_derivatives.empty()) ? 0 : n_q);
111 Q_fourth_derivatives.resize((fourth_derivatives.empty()) ? 0 : n_q);
112 bubble_values.resize((values.empty()) ? 0 : n_bubbles);
113 bubble_grads.resize((grads.empty()) ? 0 : n_bubbles);
114 bubble_grad_grads.resize((grad_grads.empty()) ? 0 : n_bubbles);
115 bubble_third_derivatives.resize((third_derivatives.empty()) ? 0 : n_bubbles);
116 bubble_fourth_derivatives.resize((fourth_derivatives.empty()) ? 0 :
117 n_bubbles);
118
119 // 1 normal vector per face, ordering consistent with GeometryInfo
120 // Normal vectors point in the +x, +y, and +z directions for
121 // consistent orientation across edges
122 std::vector<Tensor<1, dim>> normals;
123 for (const unsigned int i : GeometryInfo<dim>::face_indices())
124 {
125 Tensor<1, dim> normal;
126 normal[i / 2] = 1;
127 normals.push_back(normal);
128 }
129
130 // dim standard basis vectors for R^dim, usual ordering
131 std::vector<Tensor<1, dim>> units;
132 for (unsigned int i = 0; i < dim; ++i)
133 {
134 Tensor<1, dim> unit;
135 unit[i] = 1;
136 units.push_back(unit);
137 }
138
139 // set indices for the anisotropic polynomials to find
140 // them after polynomial_space_bubble.evaluate is called
141 std::vector<int> aniso_indices;
142 if (dim == 2)
143 {
144 aniso_indices.push_back(6);
145 aniso_indices.push_back(7);
146 aniso_indices.push_back(2);
147 aniso_indices.push_back(5);
148 }
149 else if (dim == 3)
150 {
151 aniso_indices.push_back(24);
152 aniso_indices.push_back(25);
153 aniso_indices.push_back(20);
154 aniso_indices.push_back(23);
155 aniso_indices.push_back(8);
156 aniso_indices.push_back(17);
157 }
158
159 polynomial_space_bubble.evaluate(unit_point,
160 bubble_values,
161 bubble_grads,
162 bubble_grad_grads,
163 bubble_third_derivatives,
164 bubble_fourth_derivatives);
165 polynomial_space_Q.evaluate(unit_point,
166 Q_values,
167 Q_grads,
168 Q_grad_grads,
169 Q_third_derivatives,
170 Q_fourth_derivatives);
171
172 // first dim*vertices_per_cell functions are Q_1^2 functions
173 for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
174 {
175 if (values.size() != 0)
176 {
177 values[i] = units[i % dim] * Q_values[i / dim];
178 }
179 if (grads.size() != 0)
180 {
181 grads[i] = outer_product(units[i % dim], Q_grads[i / dim]);
182 }
183 if (grad_grads.size() != 0)
184 {
185 grad_grads[i] = outer_product(units[i % dim], Q_grad_grads[i / dim]);
186 }
187 if (third_derivatives.size() != 0)
188 {
189 third_derivatives[i] =
190 outer_product(units[i % dim], Q_third_derivatives[i / dim]);
191 }
192 if (fourth_derivatives.size() != 0)
193 {
194 fourth_derivatives[i] =
195 outer_product(units[i % dim], Q_fourth_derivatives[i / dim]);
196 }
197 }
198
199 // last faces_per_cell functions are bubble functions
200 for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
201 i < dim * GeometryInfo<dim>::vertices_per_cell +
203 ++i)
204 {
205 unsigned int j =
206 i -
207 dim *
208 GeometryInfo<dim>::vertices_per_cell; // ranges 0 to faces_per_cell-1
209 if (values.size() != 0)
210 {
211 values[i] = normals[j] * bubble_values[aniso_indices[j]];
212 }
213 if (grads.size() != 0)
214 {
215 grads[i] = outer_product(normals[j], bubble_grads[aniso_indices[j]]);
216 }
217 if (grad_grads.size() != 0)
218 {
219 grad_grads[i] =
220 outer_product(normals[j], bubble_grad_grads[aniso_indices[j]]);
221 }
222 if (third_derivatives.size() != 0)
223 {
224 third_derivatives[i] =
225 outer_product(normals[j],
226 bubble_third_derivatives[aniso_indices[j]]);
227 }
228 if (fourth_derivatives.size() != 0)
229 {
230 fourth_derivatives[i] =
231 outer_product(normals[j],
232 bubble_fourth_derivatives[aniso_indices[j]]);
233 }
234 }
235}
236
237template <int dim>
238unsigned int
240{
241 (void)k;
242 Assert(k == 1, ExcNotImplemented());
243 if (dim == 2 || dim == 3)
246 // 2*4+4=12 polynomials in 2d and 3*8+6=30 polynomials in 3d
247
249 return 0;
250}
251
252
253template <int dim>
254std::unique_ptr<TensorPolynomialsBase<dim>>
256{
257 return std::make_unique<PolynomialsBernardiRaugel<dim>>(*this);
258}
259
260template class PolynomialsBernardiRaugel<1>; // to prevent errors
261template class PolynomialsBernardiRaugel<2>;
262template class PolynomialsBernardiRaugel<3>;
263
264
Definition point.h:113
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:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
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:967
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)