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