Reference documentation for deal.II version GIT f6a5d312c9 2023-10-04 08:50:02+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\}}\)
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>
22 #include <deal.II/fe/fe_nothing.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(),
200  ExcInternalError());
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,
231  ExcInternalError());
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 
265 template <int dim, int spacedim>
267  const unsigned int degree)
268  : FE_SimplexPoly<dim, spacedim>(
273  // Interface constraints are not yet implemented
275  , approximation_degree(degree)
276 {}
277 
278 
279 
280 template <int dim, int spacedim>
281 std::string
283 {
284  return "FE_SimplexP_Bubbles<" + Utilities::dim_string(dim, spacedim) + ">" +
285  "(" + std::to_string(approximation_degree) + ")";
286 }
287 
288 
289 
290 template <int dim, int spacedim>
291 std::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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
std::vector< Point< 0 > > unit_support_points< 0 >(const unsigned int)
BarycentricPolynomials< dim > get_basis(const unsigned int degree)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
FiniteElementData< dim > get_fe_data(const unsigned int degree)
std::vector< Point< dim > > unit_support_points(const unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::string to_string(const T &t)
Definition: patterns.h:2391
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:556