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