Reference documentation for deal.II version Git f102a78320 2021-10-16 14:49:02 -0400
\(\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>
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  Point<dim> centroid;
71  std::fill(centroid.begin_raw(), centroid.end_raw(), 1.0 / double(dim + 1));
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  const double q13 = 1.0 / 3.0;
89  points.emplace_back(q13, q13, 0.0);
90  points.emplace_back(q13, 0.0, q13);
91  points.emplace_back(0.0, q13, q13);
92  points.emplace_back(q13, q13, q13);
93  points.push_back(centroid);
94  }
95  return points;
96  }
97  default:
98  Assert(false, ExcNotImplemented());
99  }
100  return points;
101  }
102 
103 
104 
105  template <>
106  std::vector<Point<0>>
107  unit_support_points<0>(const unsigned int /*degree*/)
108  {
109  std::vector<Point<0>> points;
110  points.emplace_back();
111  return points;
112  }
113 
114 
115 
116  template <int dim>
118  get_basis(const unsigned int degree)
119  {
120  Point<dim> centroid;
121  std::fill(centroid.begin_raw(), centroid.end_raw(), 1.0 / double(dim + 1));
122 
123  auto M = [](const unsigned int d) {
125  };
126 
127  switch (degree)
128  {
129  // we don't need to add bubbles to P0 or P1
130  case 0:
131  case 1:
133  case 2:
134  {
135  const auto fe_p =
137  // no further work is needed in 1D
138  if (dim == 1)
139  return fe_p;
140 
141  // in 2D and 3D we add a centroid bubble function
142  auto c_bubble = BarycentricPolynomial<dim>() + 1;
143  for (unsigned int d = 0; d < dim + 1; ++d)
144  c_bubble = c_bubble * M(d);
145  c_bubble = c_bubble / c_bubble.value(centroid);
146 
147  std::vector<BarycentricPolynomial<dim>> bubble_functions;
148  if (dim == 2)
149  {
150  bubble_functions.push_back(c_bubble);
151  }
152  else if (dim == 3)
153  {
154  // need 'face bubble' functions in addition to the centroid.
155  // Furthermore we need to subtract them off from the other
156  // functions so that we end up with an interpolatory basis
157  auto b0 = 27 * M(0) * M(1) * M(2);
158  bubble_functions.push_back(b0 - b0.value(centroid) * c_bubble);
159  auto b1 = 27 * M(0) * M(1) * M(3);
160  bubble_functions.push_back(b1 - b1.value(centroid) * c_bubble);
161  auto b2 = 27 * M(0) * M(2) * M(3);
162  bubble_functions.push_back(b2 - b2.value(centroid) * c_bubble);
163  auto b3 = 27 * M(1) * M(2) * M(3);
164  bubble_functions.push_back(b3 - b3.value(centroid) * c_bubble);
165 
166  bubble_functions.push_back(c_bubble);
167  }
168 
169  // Extract out the support points for the extra bubble (both
170  // volume and face) functions:
171  const std::vector<Point<dim>> support_points =
172  unit_support_points<dim>(degree);
173  const std::vector<Point<dim>> bubble_support_points(
174  support_points.begin() + fe_p.n(), support_points.end());
175  Assert(bubble_support_points.size() == bubble_functions.size(),
176  ExcInternalError());
177  const unsigned int n_bubbles = bubble_support_points.size();
178 
179  // Assemble the final basis:
180  std::vector<BarycentricPolynomial<dim>> lump_polys;
181  for (unsigned int i = 0; i < fe_p.n(); ++i)
182  {
183  BarycentricPolynomial<dim> p = fe_p[i];
184 
185  for (unsigned int j = 0; j < n_bubbles; ++j)
186  {
187  p = p -
188  p.value(bubble_support_points[j]) * bubble_functions[j];
189  }
190 
191  lump_polys.push_back(p);
192  }
193 
194  for (auto &p : bubble_functions)
195  lump_polys.push_back(std::move(p));
196 
197  // Sanity check:
198 #ifdef DEBUG
200  for (const auto &p : lump_polys)
201  unity = unity + p;
202 
203  Point<dim> test;
204  for (unsigned int d = 0; d < dim; ++d)
205  test[d] = 2.0;
206  Assert(std::abs(unity.value(test) - 1.0) < 1e-10,
207  ExcInternalError());
208 #endif
209 
210  return BarycentricPolynomials<dim>(lump_polys);
211  }
212  default:
213  Assert(degree < 3, ExcNotImplemented());
214  }
215 
216  Assert(degree < 3, ExcNotImplemented());
217  // bogus return to placate compilers
219  }
220 
221 
222 
223  template <int dim>
225  get_fe_data(const unsigned int degree)
226  {
227  // It's not efficient, but delegate computation of the degree of the
228  // finite element (which is different from the input argument) to the
229  // basis.
230  const auto polys = get_basis<dim>(degree);
231  return FiniteElementData<dim>(get_dpo_vector<dim>(degree),
232  ReferenceCells::get_simplex<dim>(),
233  1, // n_components
234  polys.degree(),
236  }
237 } // namespace FE_P_BubblesImplementation
238 
239 
240 
241 template <int dim, int spacedim>
243  const unsigned int degree)
244  : FE_SimplexPoly<dim, spacedim>(
249  // Interface contraints are not yet implemented
252 {}
253 
254 
255 
256 template <int dim, int spacedim>
257 std::string
259 {
260  return "FE_SimplexP_Bubbles<" + Utilities::dim_string(dim, spacedim) + ">" +
262 }
263 
264 
265 
266 template <int dim, int spacedim>
267 std::unique_ptr<FiniteElement<dim, spacedim>>
269 {
270  return std::make_unique<FE_SimplexP_Bubbles<dim, spacedim>>(*this);
271 }
272 
273 // explicit instantiations
274 #include "fe_simplex_p_bubbles.inst"
275 
Number * begin_raw()
std::vector< Point< 0 > > unit_support_points< 0 >(const unsigned int)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
const unsigned int degree
Definition: fe_base.h:436
FiniteElementData< dim > get_fe_data(const unsigned int degree)
Number * end_raw()
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
#define Assert(cond, exc)
Definition: exceptions.h:1461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
std::string to_string(const T &t)
Definition: patterns.h:2329
unsigned int approximation_degree
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
BarycentricPolynomials< dim > get_basis(const unsigned int degree)
Number value(const Point< dim > &point) const
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
const std::vector< Point< dim > > & get_unit_support_points() const
Definition: fe.cc:1049
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
FE_SimplexP_Bubbles(const unsigned int degree)
static ::ExceptionBase & ExcNotImplemented()
std::vector< Point< dim > > unit_support_points(const unsigned int degree)
virtual std::string get_name() const override
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
static ::ExceptionBase & ExcInternalError()