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