Reference documentation for deal.II version 9.3.3
\(\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\}}\)
fe_simplex_p_bubbles.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 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
29namespace
30{
35 template <int dim>
36 std::vector<Point<dim>>
37 unit_support_points_fe_poly_bubbles(const unsigned int degree)
38 {
39 std::vector<Point<dim>> unit_points;
40
41 // Piecewise constants are a special case: use a support point at the
42 // centroid and only the centroid
43 if (degree == 0)
44 {
45 Point<dim> centroid;
46 std::fill(centroid.begin_raw(),
47 centroid.end_raw(),
48 1.0 / double(dim + 1));
49 unit_points.emplace_back(centroid);
50 return unit_points;
51 }
52
53 if (dim == 1)
54 {
55 // We don't really have dim = 1 support for simplex elements yet, but
56 // its convenient for populating the face array
57 Assert(degree <= 2, ExcNotImplemented());
58 if (degree >= 1)
59 {
60 unit_points.emplace_back(0.0);
61 unit_points.emplace_back(1.0);
62
63 if (degree == 2)
64 unit_points.emplace_back(0.5);
65 }
66 }
67 else if (dim == 2)
68 {
69 Assert(degree <= 2, ExcNotImplemented());
70 if (degree >= 1)
71 {
72 unit_points.emplace_back(0.0, 0.0);
73 unit_points.emplace_back(1.0, 0.0);
74 unit_points.emplace_back(0.0, 1.0);
75
76 if (degree == 2)
77 {
78 unit_points.emplace_back(0.5, 0.0);
79 unit_points.emplace_back(0.5, 0.5);
80 unit_points.emplace_back(0.0, 0.5);
81 }
82 }
83 }
84 else if (dim == 3)
85 {
86 Assert(degree <= 2, ExcNotImplemented());
87 if (degree >= 1)
88 {
89 unit_points.emplace_back(0.0, 0.0, 0.0);
90 unit_points.emplace_back(1.0, 0.0, 0.0);
91 unit_points.emplace_back(0.0, 1.0, 0.0);
92 unit_points.emplace_back(0.0, 0.0, 1.0);
93
94 if (degree == 2)
95 {
96 unit_points.emplace_back(0.5, 0.0, 0.0);
97 unit_points.emplace_back(0.5, 0.5, 0.0);
98 unit_points.emplace_back(0.0, 0.5, 0.0);
99 unit_points.emplace_back(0.0, 0.0, 0.5);
100 unit_points.emplace_back(0.5, 0.0, 0.5);
101 unit_points.emplace_back(0.0, 0.5, 0.5);
102 }
103 }
104 }
105 else
106 {
107 Assert(false, ExcNotImplemented());
108 }
109
110 return unit_points;
111 }
112} // namespace
113
115{
116 template <int dim>
117 std::vector<unsigned int>
118 get_dpo_vector(const unsigned int degree)
119 {
120 std::vector<unsigned int> dpo(dim + 1);
121 if (degree == 0)
122 {
123 dpo[dim] = 1; // single interior dof
124 }
125 else
126 {
127 Assert(degree == 1 || degree == 2, ExcNotImplemented());
128 dpo[0] = 1; // vertex dofs
129
130 if (degree == 2)
131 {
132 dpo[1] = 1; // line dofs
133
134 if (dim > 1)
135 dpo[dim] = 1; // the internal bubble function
136 if (dim == 3)
137 dpo[dim - 1] = 1; // face bubble functions
138 }
139 }
140
141 return dpo;
142 }
143
144
145
146 template <int dim>
147 std::vector<Point<dim>>
148 unit_support_points(const unsigned int degree)
149 {
150 Assert(degree < 3, ExcNotImplemented());
151 std::vector<Point<dim>> points =
152 unit_support_points_fe_poly_bubbles<dim>(degree);
153
154 Point<dim> centroid;
155 std::fill(centroid.begin_raw(), centroid.end_raw(), 1.0 / double(dim + 1));
156
157 switch (dim)
158 {
159 case 1:
160 // nothing more to do
161 return points;
162 case 2:
163 {
164 if (degree == 2)
165 points.push_back(centroid);
166 return points;
167 }
168 case 3:
169 {
170 if (degree == 2)
171 {
172 const double q13 = 1.0 / 3.0;
173 points.emplace_back(q13, q13, 0.0);
174 points.emplace_back(q13, 0.0, q13);
175 points.emplace_back(0.0, q13, q13);
176 points.emplace_back(q13, q13, q13);
177 points.push_back(centroid);
178 }
179 return points;
180 }
181 default:
182 Assert(false, ExcNotImplemented());
183 }
184 return points;
185 }
186
187
188
189 template <int dim>
191 get_basis(const unsigned int degree)
192 {
193 Point<dim> centroid;
194 std::fill(centroid.begin_raw(), centroid.end_raw(), 1.0 / double(dim + 1));
195
196 auto M = [](const unsigned int d) {
198 };
199
200 switch (degree)
201 {
202 // we don't need to add bubbles to P0 or P1
203 case 0:
204 case 1:
206 case 2:
207 {
208 const auto fe_p =
210 // no further work is needed in 1D
211 if (dim == 1)
212 return fe_p;
213
214 // in 2D and 3D we add a centroid bubble function
215 auto c_bubble = BarycentricPolynomial<dim>() + 1;
216 for (unsigned int d = 0; d < dim + 1; ++d)
217 c_bubble = c_bubble * M(d);
218 c_bubble = c_bubble / c_bubble.value(centroid);
219
220 std::vector<BarycentricPolynomial<dim>> bubble_functions;
221 if (dim == 2)
222 {
223 bubble_functions.push_back(c_bubble);
224 }
225 else if (dim == 3)
226 {
227 // need 'face bubble' functions in addition to the centroid.
228 // Furthermore we need to subtract them off from the other
229 // functions so that we end up with an interpolatory basis
230 auto b0 = 27 * M(0) * M(1) * M(2);
231 bubble_functions.push_back(b0 - b0.value(centroid) * c_bubble);
232 auto b1 = 27 * M(0) * M(1) * M(3);
233 bubble_functions.push_back(b1 - b1.value(centroid) * c_bubble);
234 auto b2 = 27 * M(0) * M(2) * M(3);
235 bubble_functions.push_back(b2 - b2.value(centroid) * c_bubble);
236 auto b3 = 27 * M(1) * M(2) * M(3);
237 bubble_functions.push_back(b3 - b3.value(centroid) * c_bubble);
238
239 bubble_functions.push_back(c_bubble);
240 }
241
242 // Extract out the support points for the extra bubble (both
243 // volume and face) functions:
244 const std::vector<Point<dim>> support_points =
245 unit_support_points<dim>(degree);
246 const std::vector<Point<dim>> bubble_support_points(
247 support_points.begin() + fe_p.n(), support_points.end());
248 Assert(bubble_support_points.size() == bubble_functions.size(),
250 const unsigned int n_bubbles = bubble_support_points.size();
251
252 // Assemble the final basis:
253 std::vector<BarycentricPolynomial<dim>> lump_polys;
254 for (unsigned int i = 0; i < fe_p.n(); ++i)
255 {
256 BarycentricPolynomial<dim> p = fe_p[i];
257
258 for (unsigned int j = 0; j < n_bubbles; ++j)
259 {
260 p = p -
261 p.value(bubble_support_points[j]) * bubble_functions[j];
262 }
263
264 lump_polys.push_back(p);
265 }
266
267 for (auto &p : bubble_functions)
268 lump_polys.push_back(std::move(p));
269
270 // Sanity check:
271#ifdef DEBUG
273 for (const auto &p : lump_polys)
274 unity = unity + p;
275
276 Point<dim> test;
277 for (unsigned int d = 0; d < dim; ++d)
278 test[d] = 2.0;
279 Assert(std::abs(unity.value(test) - 1.0) < 1e-10,
281#endif
282
283 return BarycentricPolynomials<dim>(lump_polys);
284 }
285 default:
286 Assert(degree < 3, ExcNotImplemented());
287 }
288
289 Assert(degree < 3, ExcNotImplemented());
290 // bogus return to placate compilers
292 }
293
294
295
296 template <int dim>
298 get_fe_data(const unsigned int degree)
299 {
300 // It's not efficient, but delegate computation of the degree of the
301 // finite element (which is different from the input argument) to the
302 // basis.
303 const auto polys = get_basis<dim>(degree);
304 return FiniteElementData<dim>(get_dpo_vector<dim>(degree),
305 ReferenceCells::get_simplex<dim>(),
306 1, // n_components
307 polys.degree(),
309 }
310} // namespace FE_P_BubblesImplementation
311
312
313
314template <int dim, int spacedim>
316 const unsigned int degree)
317 : ::FE_Poly<dim, spacedim>(
320 std::vector<bool>(
321 FE_P_BubblesImplementation::get_fe_data<dim>(degree).dofs_per_cell,
322 true),
323 std::vector<ComponentMask>(
324 FE_P_BubblesImplementation::get_fe_data<dim>(degree).dofs_per_cell,
325 std::vector<bool>(1, true)))
326 , approximation_degree(degree)
327{
328 this->unit_support_points =
329 FE_P_BubblesImplementation::unit_support_points<dim>(degree);
330
331 // TODO
332 // this->unit_face_support_points =
333 // unit_face_support_points_fe_poly<dim>(degree);
334}
335
336
337
338template <int dim, int spacedim>
339std::string
341{
342 return "FE_SimplexP_Bubbles<" + Utilities::dim_string(dim, spacedim) + ">" +
343 "(" + std::to_string(approximation_degree) + ")";
344}
345
346
347
348template <int dim, int spacedim>
349void
352 const std::vector<Vector<double>> &support_point_values,
353 std::vector<double> & nodal_values) const
354{
355 AssertDimension(support_point_values.size(),
356 this->get_unit_support_points().size());
357 AssertDimension(support_point_values.size(), nodal_values.size());
358 AssertDimension(this->dofs_per_cell, nodal_values.size());
359
360 for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
361 {
362 AssertDimension(support_point_values[i].size(), 1);
363
364 nodal_values[i] = support_point_values[i](0);
365 }
366}
367
368
369
370template <int dim, int spacedim>
371std::unique_ptr<FiniteElement<dim, spacedim>>
373{
374 return std::make_unique<FE_SimplexP_Bubbles<dim, spacedim>>(*this);
375}
376
377// explicit instantiations
378#include "fe_simplex_p_bubbles.inst"
379
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 void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
virtual std::string get_name() const override
const unsigned int degree
Definition: fe_base.h:435
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2442
Definition: point.h:111
Number * begin_raw()
Number * end_raw()
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
std::string to_string(const T &t)
Definition: patterns.h:2329
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
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)
FiniteElementData< dim > get_fe_data(const unsigned int degree)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:558
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)