Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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_q_dg0.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2022 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 
20 
22 
23 #include <deal.II/fe/fe_dgq.h>
24 #include <deal.II/fe/fe_nothing.h>
25 #include <deal.II/fe/fe_q_dg0.h>
26 
27 #include <memory>
28 #include <sstream>
29 #include <vector>
30 
32 
33 
34 template <int dim, int spacedim>
35 FE_Q_DG0<dim, spacedim>::FE_Q_DG0(const unsigned int degree)
36  : FE_Q_Base<dim, spacedim>(TensorProductPolynomialsConst<dim>(
38  QGaussLobatto<1>(degree + 1).get_points())),
39  FiniteElementData<dim>(get_dpo_vector(degree),
40  1,
41  degree,
42  FiniteElementData<dim>::L2),
43  get_riaf_vector(degree))
44 {
45  Assert(degree > 0,
46  ExcMessage("This element can only be used for polynomial degrees "
47  "greater than zero"));
48 
49  this->initialize(QGaussLobatto<1>(degree + 1).get_points());
50 
51  // adjust unit support point for discontinuous node
53  for (unsigned int d = 0; d < dim; ++d)
54  point[d] = 0.5;
55  this->unit_support_points.push_back(point);
57 }
58 
59 
60 
61 template <int dim, int spacedim>
63  : FE_Q_Base<dim, spacedim>(
65  Polynomials::generate_complete_Lagrange_basis(points.get_points())),
66  FiniteElementData<dim>(get_dpo_vector(points.size() - 1),
67  1,
68  points.size() - 1,
69  FiniteElementData<dim>::L2),
70  get_riaf_vector(points.size() - 1))
71 {
72  const int degree = points.size() - 1;
73  (void)degree;
74 
75  Assert(degree > 0,
76  ExcMessage("This element can only be used for polynomial degrees "
77  "at least zero"));
78 
79  this->initialize(points.get_points());
80 
81  // adjust unit support point for discontinuous node
83  for (unsigned int d = 0; d < dim; ++d)
84  point[d] = 0.5;
85  this->unit_support_points.push_back(point);
87 }
88 
89 
90 
91 template <int dim, int spacedim>
92 std::string
94 {
95  // note that the FETools::get_fe_by_name function depends on the
96  // particular format of the string this function returns, so they have to be
97  // kept in synch
98 
99  std::ostringstream namebuf;
100  bool type = true;
101  const unsigned int n_points = this->degree + 1;
102  std::vector<double> points(n_points);
103  const unsigned int dofs_per_cell = this->n_dofs_per_cell();
104  const std::vector<Point<dim>> &unit_support_points =
105  this->unit_support_points;
106  unsigned int index = 0;
107 
108  // Decode the support points in one coordinate direction.
109  for (unsigned int j = 0; j < dofs_per_cell; ++j)
110  {
111  if ((dim > 1) ? (unit_support_points[j](1) == 0 &&
112  ((dim > 2) ? unit_support_points[j](2) == 0 : true)) :
113  true)
114  {
115  if (index == 0)
116  points[index] = unit_support_points[j](0);
117  else if (index == 1)
118  points[n_points - 1] = unit_support_points[j](0);
119  else
120  points[index - 1] = unit_support_points[j](0);
121 
122  index++;
123  }
124  }
125  // Do not consider the discontinuous node for dimension 1
126  Assert(index == n_points || (dim == 1 && index == n_points + 1),
127  ExcMessage(
128  "Could not decode support points in one coordinate direction."));
129 
130  // Check whether the support points are equidistant.
131  for (unsigned int j = 0; j < n_points; ++j)
132  if (std::fabs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
133  {
134  type = false;
135  break;
136  }
137 
138  if (type == true)
139  {
140  if (this->degree > 2)
141  namebuf << "FE_Q_DG0<" << Utilities::dim_string(dim, spacedim)
142  << ">(QIterated(QTrapezoid()," << this->degree << "))";
143  else
144  namebuf << "FE_Q_DG0<" << Utilities::dim_string(dim, spacedim) << ">("
145  << this->degree << ")";
146  }
147  else
148  {
149  // Check whether the support points come from QGaussLobatto.
150  const QGaussLobatto<1> points_gl(n_points);
151  type = true;
152  for (unsigned int j = 0; j < n_points; ++j)
153  if (points[j] != points_gl.point(j)(0))
154  {
155  type = false;
156  break;
157  }
158  if (type == true)
159  namebuf << "FE_Q_DG0<" << Utilities::dim_string(dim, spacedim) << ">("
160  << this->degree << ")";
161  else
162  namebuf << "FE_Q_DG0<" << Utilities::dim_string(dim, spacedim)
163  << ">(QUnknownNodes(" << this->degree << "))";
164  }
165  return namebuf.str();
166 }
167 
168 
169 
170 template <int dim, int spacedim>
171 std::unique_ptr<FiniteElement<dim, spacedim>>
173 {
174  return std::make_unique<FE_Q_DG0<dim, spacedim>>(*this);
175 }
176 
177 
178 
179 template <int dim, int spacedim>
180 void
182  const std::vector<Vector<double>> &support_point_values,
183  std::vector<double> & nodal_dofs) const
184 {
185  Assert(support_point_values.size() == this->unit_support_points.size(),
186  ExcDimensionMismatch(support_point_values.size(),
187  this->unit_support_points.size()));
188  Assert(nodal_dofs.size() == this->n_dofs_per_cell(),
189  ExcDimensionMismatch(nodal_dofs.size(), this->n_dofs_per_cell()));
190  Assert(support_point_values[0].size() == this->n_components(),
191  ExcDimensionMismatch(support_point_values[0].size(),
192  this->n_components()));
193 
194  for (unsigned int i = 0; i < this->n_dofs_per_cell() - 1; ++i)
195  {
196  const std::pair<unsigned int, unsigned int> index =
197  this->system_to_component_index(i);
198  nodal_dofs[i] = support_point_values[i](index.first);
199  }
200 
201  // We don't need the discontinuous function for local interpolation
202  nodal_dofs.back() = 0.;
203 }
204 
205 
206 
207 template <int dim, int spacedim>
208 void
210  const FiniteElement<dim, spacedim> &x_source_fe,
211  FullMatrix<double> & interpolation_matrix) const
212 {
213  // this is only implemented, if the source FE is also a Q_DG0 element
214  using FEQDG0 = FE_Q_DG0<dim, spacedim>;
215 
216  AssertThrow(
217  (x_source_fe.get_name().find("FE_Q_DG0<") == 0) ||
218  (dynamic_cast<const FEQDG0 *>(&x_source_fe) != nullptr),
220 
221  Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
222  ExcDimensionMismatch(interpolation_matrix.m(),
223  this->n_dofs_per_cell()));
224  Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
225  ExcDimensionMismatch(interpolation_matrix.m(),
226  x_source_fe.n_dofs_per_cell()));
227 
229  x_source_fe, interpolation_matrix);
230 }
231 
232 
233 
234 template <int dim, int spacedim>
235 std::vector<bool>
237 {
238  std::vector<bool> riaf(Utilities::fixed_power<dim>(deg + 1) + 1, false);
239  riaf.back() = true;
240  return riaf;
241 }
242 
243 
244 
245 template <int dim, int spacedim>
246 std::vector<unsigned int>
248 {
249  std::vector<unsigned int> dpo(dim + 1, 1U);
250  for (unsigned int i = 1; i < dpo.size(); ++i)
251  dpo[i] = dpo[i - 1] * (deg - 1);
252 
253  dpo[dim]++; // we need an additional DG0-node for a dim-dimensional object
254  return dpo;
255 }
256 
257 
258 
259 template <int dim, int spacedim>
260 bool
262  const unsigned int shape_index,
263  const unsigned int face_index) const
264 {
265  // discontinuous function has support on all faces
266  if (shape_index == this->n_dofs_per_cell() - 1)
267  return true;
268  else
270  face_index);
271 }
272 
273 
274 
275 template <int dim, int spacedim>
276 std::pair<Table<2, bool>, std::vector<unsigned int>>
278 {
279  Table<2, bool> constant_modes(2, this->n_dofs_per_cell());
280 
281  // 1 represented by FE_Q part
282  for (unsigned int i = 0; i < this->n_dofs_per_cell() - 1; ++i)
283  constant_modes(0, i) = true;
284 
285  // 1 represented by DG0 part
286  constant_modes(1, this->n_dofs_per_cell() - 1) = true;
287 
288  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
289  constant_modes, std::vector<unsigned int>(2, 0));
290 }
291 
292 
293 
294 template <int dim, int spacedim>
297  const FiniteElement<dim, spacedim> &fe_other,
298  const unsigned int codim) const
299 {
300  Assert(codim <= dim, ExcImpossibleInDim(dim));
301 
302  // vertex/line/face domination
303  // (if fe_other is derived from FE_DGQ)
304  // ------------------------------------
305  if (codim > 0)
306  if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
307  // there are no requirements between continuous and discontinuous elements
309 
310  // vertex/line/face domination
311  // (if fe_other is not derived from FE_DGQ)
312  // & cell domination
313  // ----------------------------------------
314  if (const FE_Q_DG0<dim, spacedim> *fe_dg0_other =
315  dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe_other))
316  {
317  if (this->degree < fe_dg0_other->degree)
319  else if (this->degree == fe_dg0_other->degree)
321  else
323  }
324  else if (const FE_Nothing<dim> *fe_nothing =
325  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
326  {
327  if (fe_nothing->is_dominating())
329  else
330  // the FE_Nothing has no degrees of freedom and it is typically used
331  // in a context where we don't require any continuity along the
332  // interface
334  }
335 
336  Assert(false, ExcNotImplemented());
338 }
339 
340 
341 // explicit instantiations
342 #include "fe_q_dg0.inst"
343 
Definition: fe_dgq.h:113
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_q_base.cc:512
void initialize(const std::vector< Point< 1 >> &support_points_1d)
Definition: fe_q_base.cc:436
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_q_base.cc:1592
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_q_dg0.cc:209
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_q_dg0.cc:172
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_q_dg0.cc:277
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
Definition: fe_q_dg0.cc:181
FE_Q_DG0(const unsigned int p)
Definition: fe_q_dg0.cc:35
virtual std::string get_name() const override
Definition: fe_q_dg0.cc:93
static std::vector< bool > get_riaf_vector(const unsigned int degree)
Definition: fe_q_dg0.cc:236
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_q_dg0.cc:261
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_q_dg0.cc:296
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_q_dg0.cc:247
const unsigned int degree
Definition: fe_data.h:449
unsigned int n_dofs_per_cell() const
virtual std::string get_name() const =0
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2440
size_type n() const
size_type m() const
Definition: point.h:110
const std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
Expression fabs(const Expression &x)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:189
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)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:702
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:556
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)