Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10: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\}}\)
Loading...
Searching...
No Matches
fe_simplex_p_bubbles.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 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#include <deal.II/base/config.h>
16
19
20#include <deal.II/fe/fe_dgq.h>
22#include <deal.II/fe/fe_q.h>
24#include <deal.II/fe/fe_tools.h>
25
27
29{
30 template <int dim>
31 std::vector<unsigned int>
32 get_dpo_vector(const unsigned int degree)
33 {
34 std::vector<unsigned int> dpo(dim + 1);
35 if (degree == 0)
36 {
37 dpo[dim] = 1; // single interior dof
38 }
39 else
40 {
41 Assert(degree == 1 || degree == 2, ExcNotImplemented());
42 dpo[0] = 1; // vertex dofs
43
44 if (degree == 2)
45 {
46 dpo[1] = 1; // line dofs
47
48 if (dim > 1)
49 dpo[dim] = 1; // the internal bubble function
50 if (dim == 3)
51 dpo[dim - 1] = 1; // face bubble functions
52 }
53 }
54
55 return dpo;
56 }
57
58
59
60 template <int dim>
61 std::vector<Point<dim>>
62 unit_support_points(const unsigned int degree)
63 {
64 Assert(degree < 3, ExcNotImplemented());
65 // Start with the points used by FE_SimplexP, and then add bubbles.
66 FE_SimplexP<dim> fe_p(degree);
67 std::vector<Point<dim>> points = fe_p.get_unit_support_points();
68
69 const auto reference_cell = fe_p.reference_cell();
70 const Point<dim> centroid = reference_cell.template barycenter<dim>();
71
72 switch (dim)
73 {
74 case 1:
75 // nothing more to do
76 return points;
77 case 2:
78 {
79 if (degree == 2)
80 points.push_back(centroid);
81 return points;
82 }
83 case 3:
84 {
85 if (degree == 2)
86 {
87 for (const auto &face_no : reference_cell.face_indices())
88 {
89 Point<dim> midpoint;
90 for (const auto face_vertex_no :
91 reference_cell.face_reference_cell(0).vertex_indices())
92 {
93 const auto vertex_no =
94 reference_cell.face_to_cell_vertices(
95 face_no,
96 face_vertex_no,
98
99 midpoint +=
100 reference_cell.template vertex<dim>(vertex_no);
101 }
102
103 midpoint /=
104 reference_cell.face_reference_cell(0).n_vertices();
105 points.push_back(midpoint);
106 }
107
108 points.push_back(centroid);
109 }
110 return points;
111 }
112 default:
114 }
115 return points;
116 }
117
118
119
120 template <>
121 std::vector<Point<0>>
122 unit_support_points<0>(const unsigned int /*degree*/)
123 {
124 std::vector<Point<0>> points;
125 points.emplace_back();
126 return points;
127 }
128
129
130
131 template <int dim>
133 get_basis(const unsigned int degree)
134 {
135 const auto reference_cell = ReferenceCells::get_simplex<dim>();
136 const Point<dim> centroid = reference_cell.template barycenter<dim>();
137
138 auto M = [](const unsigned int d) {
140 };
141
142 switch (degree)
143 {
144 // we don't need to add bubbles to P0 or P1
145 case 0:
146 case 1:
148 case 2:
149 {
150 const auto fe_p =
152 // no further work is needed in 1d
153 if (dim == 1)
154 return fe_p;
155
156 // in 2d and 3d we add a centroid bubble function
157 auto c_bubble = BarycentricPolynomial<dim>() + 1;
158 for (const auto &vertex : reference_cell.vertex_indices())
159 c_bubble = c_bubble * M(vertex);
160 c_bubble = c_bubble / c_bubble.value(centroid);
161
162 std::vector<BarycentricPolynomial<dim>> bubble_functions;
163 if (dim == 2)
164 {
165 bubble_functions.push_back(c_bubble);
166 }
167 else if (dim == 3)
168 {
169 // need 'face bubble' functions in addition to the centroid.
170 // Furthermore we need to subtract them off from the other
171 // functions so that we end up with an interpolatory basis
172 for (const auto &face_no : reference_cell.face_indices())
173 {
174 std::vector<unsigned int> vertices;
175 for (const auto face_vertex_no :
176 reference_cell.face_reference_cell(0).vertex_indices())
177 vertices.push_back(reference_cell.face_to_cell_vertices(
178 face_no,
179 face_vertex_no,
181
182 Assert(vertices.size() == 3, ExcInternalError());
183 auto b =
184 27.0 * M(vertices[0]) * M(vertices[1]) * M(vertices[2]);
185 bubble_functions.push_back(b -
186 b.value(centroid) * c_bubble);
187 }
188
189 bubble_functions.push_back(c_bubble);
190 }
191
192 // Extract out the support points for the extra bubble (both
193 // volume and face) functions:
194 const std::vector<Point<dim>> support_points =
195 unit_support_points<dim>(degree);
196 const std::vector<Point<dim>> bubble_support_points(
197 support_points.begin() + fe_p.n(), support_points.end());
198 Assert(bubble_support_points.size() == bubble_functions.size(),
200 const unsigned int n_bubbles = bubble_support_points.size();
201
202 // Assemble the final basis:
203 std::vector<BarycentricPolynomial<dim>> lump_polys;
204 for (unsigned int i = 0; i < fe_p.n(); ++i)
205 {
206 BarycentricPolynomial<dim> p = fe_p[i];
207
208 for (unsigned int j = 0; j < n_bubbles; ++j)
209 {
210 p = p -
211 p.value(bubble_support_points[j]) * bubble_functions[j];
212 }
213
214 lump_polys.push_back(p);
215 }
216
217 for (auto &p : bubble_functions)
218 lump_polys.push_back(std::move(p));
219
220 // Sanity check:
221#ifdef DEBUG
223 for (const auto &p : lump_polys)
224 unity = unity + p;
225
226 Point<dim> test;
227 for (unsigned int d = 0; d < dim; ++d)
228 test[d] = 2.0;
229 Assert(std::abs(unity.value(test) - 1.0) < 1e-10,
231#endif
232
233 return BarycentricPolynomials<dim>(lump_polys);
234 }
235 default:
236 Assert(degree < 3, ExcNotImplemented());
237 }
238
239 Assert(degree < 3, ExcNotImplemented());
240 // bogus return to placate compilers
242 }
243
244
245
246 template <int dim>
248 get_fe_data(const unsigned int degree)
249 {
250 // It's not efficient, but delegate computation of the degree of the
251 // finite element (which is different from the input argument) to the
252 // basis.
253 const auto polys = get_basis<dim>(degree);
254 return FiniteElementData<dim>(get_dpo_vector<dim>(degree),
255 ReferenceCells::get_simplex<dim>(),
256 1, // n_components
257 polys.degree(),
259 }
260} // namespace FE_P_BubblesImplementation
261
262
263
264template <int dim, int spacedim>
266 const unsigned int degree)
267 : FE_SimplexPoly<dim, spacedim>(
268 FE_P_BubblesImplementation::get_basis<dim>(degree),
269 FE_P_BubblesImplementation::get_fe_data<dim>(degree),
270 FE_P_BubblesImplementation::unit_support_points<dim>(degree),
272 // Interface constraints are not yet implemented
274 , approximation_degree(degree)
275{}
276
277
278
279template <int dim, int spacedim>
280std::string
282{
283 return "FE_SimplexP_Bubbles<" + Utilities::dim_string(dim, spacedim) + ">" +
284 "(" + std::to_string(approximation_degree) + ")";
285}
286
287
288
289template <int dim, int spacedim>
290std::unique_ptr<FiniteElement<dim, spacedim>>
292{
293 return std::make_unique<FE_SimplexP_Bubbles<dim, spacedim>>(*this);
294}
295
296// explicit instantiations
297#include "fe_simplex_p_bubbles.inst"
298
Number value(const Point< dim > &point) const
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
FE_SimplexP_Bubbles(const unsigned int degree)
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::string get_name() const override
const unsigned int degree
Definition fe_data.h:452
ReferenceCell reference_cell() const
const std::vector< Point< dim > > & get_unit_support_points() const
Definition point.h:111
static constexpr unsigned char default_combined_face_orientation()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define DEAL_II_NOT_IMPLEMENTED()
BarycentricPolynomials< dim > get_basis(const unsigned int degree)
std::vector< Point< dim > > unit_support_points(const unsigned int degree)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
std::vector< Point< 0 > > unit_support_points< 0 >(const unsigned int)
FiniteElementData< dim > get_fe_data(const unsigned int degree)
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:555
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)