Reference documentation for deal.II version 9.4.0
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
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>
25#include <deal.II/fe/fe_q_dg0.h>
26
27#include <memory>
28#include <sstream>
29#include <vector>
30
32
33
34template <int dim, int spacedim>
35FE_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())),
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
61template <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
91template <int dim, int spacedim>
92std::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 =
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),
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
170template <int dim, int spacedim>
171std::unique_ptr<FiniteElement<dim, spacedim>>
173{
174 return std::make_unique<FE_Q_DG0<dim, spacedim>>(*this);
175}
176
177
178
179template <int dim, int spacedim>
180void
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
207template <int dim, int spacedim>
208void
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
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
234template <int dim, int spacedim>
235std::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
245template <int dim, int spacedim>
246std::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
259template <int dim, int spacedim>
260bool
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
275template <int dim, int spacedim>
276std::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
294template <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:111
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
void initialize(const std::vector< Point< 1 > > &support_points_1d)
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
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
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 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
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:2451
size_type n() const
size_type m() const
Definition: point.h:111
const Point< dim > & point(const unsigned int i) const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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:190
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:558
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)