Reference documentation for deal.II version 9.2.0
\(\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 - 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 
22 
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 
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>(
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
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
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 
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 
FE_DGQ
Definition: fe_dgq.h:112
FE_Nothing< dim >
FE_Q_DG0::FE_Q_DG0
FE_Q_DG0(const unsigned int p)
Definition: fe_q_dg0.cc:36
Polynomials
Definition: polynomial.h:41
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
FiniteElementData
Definition: fe_base.h:147
FE_Q_DG0
Definition: fe_q_dg0.h:239
FullMatrix::m
size_type m() const
FE_Q_Base
Definition: fe_q_base.h:44
FE_Q_DG0::get_riaf_vector
static std::vector< bool > get_riaf_vector(const unsigned int degree)
Definition: fe_q_dg0.cc:237
LAPACKSupport::U
static const char U
Definition: lapack_support.h:167
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
FullMatrix::n
size_type n() const
FE_Q_DG0::clone
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_q_dg0.cc:174
FE_Q_DG0::has_support_on_face
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_q_dg0.cc:262
FETools::get_interpolation_matrix
void get_interpolation_matrix(const FiniteElement< dim, spacedim > &fe1, const FiniteElement< dim, spacedim > &fe2, FullMatrix< number > &interpolation_matrix)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
quadrature_lib.h
fe_nothing.h
Polynomials::generate_complete_Lagrange_basis
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:834
FiniteElementDomination::either_element_can_dominate
@ either_element_can_dominate
Definition: fe_base.h:100
Table
Definition: table.h:699
FE_Q_DG0::get_name
virtual std::string get_name() const override
Definition: fe_q_dg0.cc:95
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition: symengine_math.cc:273
FiniteElementData::degree
const unsigned int degree
Definition: fe_base.h:298
FiniteElementDomination::no_requirements
@ no_requirements
Definition: fe_base.h:104
Quadrature::point
const Point< dim > & point(const unsigned int i) const
FiniteElementDomination::neither_element_dominates
@ neither_element_dominates
Definition: fe_base.h:96
FiniteElementDomination::other_element_dominates
@ other_element_dominates
Definition: fe_base.h:92
FE_Q_DG0::get_interpolation_matrix
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_q_dg0.cc:211
FE_Q_DG0::get_dpo_vector
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_q_dg0.cc:248
FiniteElementDomination::Domination
Domination
Definition: fe_base.h:83
FE_Q_DG0::convert_generalized_support_point_values_to_dof_values
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
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_base.h:282
Utilities::dim_string
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:559
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
fe_dgq.h
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
QGaussLobatto
Definition: quadrature_lib.h:76
FE_Q_Base< TensorProductPolynomialsConst< dim >, dim, dim >::initialize
void initialize(const std::vector< Point< 1 >> &support_points_1d)
Definition: fe_q_base.cc:425
FiniteElement
Definition: fe.h:648
LocalIntegrators::L2::L2
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:171
FE_Q_DG0::get_constant_modes
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_q_dg0.cc:278
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
dof_accessor.h
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
FiniteElement::get_name
virtual std::string get_name() const =0
quadrature.h
template_constraints.h
FiniteElementDomination::this_element_dominates
@ this_element_dominates
Definition: fe_base.h:88
Quadrature::size
unsigned int size() const
FullMatrix< double >
memory.h
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Quadrature< 1 >
Vector< double >
TensorProductPolynomialsConst
Definition: tensor_product_polynomials.h:39
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
qprojector.h
FE_Q_DG0::compare_for_domination
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
fe_q_dg0.h
FiniteElement::unit_support_points
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2452