Reference documentation for deal.II version GIT relicensing-462-g03c82a14bd 2024-04-21 16:20: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_trace.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2014 - 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#include <deal.II/base/config.h>
16
19
22#include <deal.II/fe/fe_q.h>
23#include <deal.II/fe/fe_tools.h>
24#include <deal.II/fe/fe_trace.h>
25
26#include <fstream>
27#include <iostream>
28#include <memory>
29#include <sstream>
30
32
33
34
35template <int dim, int spacedim>
36FE_TraceQ<dim, spacedim>::FE_TraceQ(const unsigned int degree)
37 : FE_PolyFace<TensorProductPolynomials<dim - 1>, 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 std::vector<bool>(1, true))
46 , fe_q(degree)
47{
48 Assert(degree > 0,
49 ExcMessage("FE_Trace can only be used for polynomial degrees "
50 "greater than zero"));
51 this->poly_space.set_numbering(
52 FETools::hierarchic_to_lexicographic_numbering<dim - 1>(degree));
53
54 // Initialize face support points
55 AssertDimension(this->n_unique_faces(), fe_q.n_unique_faces());
56 for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
57 this->unit_face_support_points[face_no] =
58 fe_q.get_unit_face_support_points(face_no);
59
60 // initialize unit support points (this makes it possible to assign initial
61 // values to FE_TraceQ). Note that we simply take the points of fe_q but
62 // skip the last ones which are associated with the interior of FE_Q.
63 this->unit_support_points.resize(this->n_dofs_per_cell());
64 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
65 this->unit_support_points[i] = fe_q.get_unit_support_points()[i];
66
67 // Initialize constraint matrices
68 this->interface_constraints = fe_q.constraints();
69}
70
71
72
73template <int dim, int spacedim>
74std::unique_ptr<FiniteElement<dim, spacedim>>
76{
77 return std::make_unique<FE_TraceQ<dim, spacedim>>(this->degree);
78}
79
80
81
82template <int dim, int spacedim>
83std::string
85{
86 // note that the FETools::get_fe_by_name function depends on the
87 // particular format of the string this function returns, so they have to be
88 // kept in synch
89
90 std::ostringstream namebuf;
91 namebuf << "FE_TraceQ<" << Utilities::dim_string(dim, spacedim) << ">("
92 << this->degree << ")";
93
94 return namebuf.str();
95}
96
97
98
99template <int dim, int spacedim>
100bool
102 const unsigned int shape_index,
103 const unsigned int face_index) const
104{
105 AssertIndexRange(shape_index, this->n_dofs_per_cell());
107
108 // FE_TraceQ shares the numbering of elemental degrees of freedom with FE_Q
109 // except for the missing interior ones (quad dofs in 2d and hex dofs in
110 // 3d). Therefore, it is safe to ask fe_q for the corresponding
111 // information. The assertion 'shape_index < this->n_dofs_per_cell()' will
112 // make sure that we only access the trace dofs.
113 return fe_q.has_support_on_face(shape_index, face_index);
114}
115
116
117
118template <int dim, int spacedim>
119std::pair<Table<2, bool>, std::vector<unsigned int>>
121{
122 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
123 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
124 constant_modes(0, i) = true;
125 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
126 constant_modes, std::vector<unsigned int>(1, 0));
127}
128
129template <int dim, int spacedim>
130void
133 const std::vector<Vector<double>> &support_point_values,
134 std::vector<double> &nodal_values) const
135{
136 AssertDimension(support_point_values.size(),
137 this->get_unit_support_points().size());
138 AssertDimension(support_point_values.size(), nodal_values.size());
139 AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
140
141 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
142 {
143 AssertDimension(support_point_values[i].size(), 1);
144
145 nodal_values[i] = support_point_values[i](0);
146 }
147}
148
149
150template <int dim, int spacedim>
151std::vector<unsigned int>
153{
154 // This constructs FE_TraceQ in exactly the same way as FE_Q except for the
155 // interior degrees of freedom that are not present here (line in 1d, quad
156 // in 2d, hex in 3d).
157 AssertThrow(deg > 0, ExcMessage("FE_TraceQ needs to be of degree > 0."));
158 std::vector<unsigned int> dpo(dim + 1, 1U);
159 dpo[dim] = 0;
160 dpo[0] = 1;
161 for (unsigned int i = 1; i < dim; ++i)
162 dpo[i] = dpo[i - 1] * (deg - 1);
163 return dpo;
164}
165
166
167
168template <int dim, int spacedim>
169bool
171{
172 return fe_q.hp_constraints_are_implemented();
173}
174
175
176template <int dim, int spacedim>
179 const FiniteElement<dim, spacedim> &fe_other,
180 const unsigned int codim) const
181{
182 Assert(codim <= dim, ExcImpossibleInDim(dim));
183 (void)codim;
184
185 // vertex/line/face/cell domination
186 // --------------------------------
187 if (const FE_TraceQ<dim, spacedim> *fe_traceq_other =
188 dynamic_cast<const FE_TraceQ<dim, spacedim> *>(&fe_other))
189 {
190 if (this->degree < fe_traceq_other->degree)
192 else if (this->degree == fe_traceq_other->degree)
194 else
196 }
197 else if (const FE_Nothing<dim> *fe_nothing =
198 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
199 {
200 if (fe_nothing->is_dominating())
202 else
203 // the FE_Nothing has no degrees of freedom and it is typically used
204 // in a context where we don't require any continuity along the
205 // interface
207 }
208
211}
212
213
214
215template <int dim, int spacedim>
216void
218 const FiniteElement<dim, spacedim> &source_fe,
219 FullMatrix<double> &interpolation_matrix,
220 const unsigned int face_no) const
221{
222 get_subface_interpolation_matrix(source_fe,
224 interpolation_matrix,
225 face_no);
226}
227
228
229
230template <int dim, int spacedim>
231void
233 const FiniteElement<dim, spacedim> &x_source_fe,
234 const unsigned int subface,
235 FullMatrix<double> &interpolation_matrix,
236 const unsigned int face_no) const
237{
238 // this is the code from FE_FaceQ
239 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
240 ExcDimensionMismatch(interpolation_matrix.n(),
241 this->n_dofs_per_face(face_no)));
242 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
243 ExcDimensionMismatch(interpolation_matrix.m(),
244 x_source_fe.n_dofs_per_face(face_no)));
245
246 // see if source is a FaceQ element
247 if (const FE_TraceQ<dim, spacedim> *source_fe =
248 dynamic_cast<const FE_TraceQ<dim, spacedim> *>(&x_source_fe))
249 {
250 fe_q.get_subface_interpolation_matrix(source_fe->fe_q,
251 subface,
252 interpolation_matrix,
253 face_no);
254 }
255 else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
256 {
257 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
258 }
259 else
261 false,
262 (typename FiniteElement<dim,
263 spacedim>::ExcInterpolationNotImplemented()));
264}
265
266
267
268template <int spacedim>
269FE_TraceQ<1, spacedim>::FE_TraceQ(const unsigned int degree)
270 : FE_FaceQ<1, spacedim>(degree)
271{}
272
273
274
275template <int spacedim>
276std::string
278{
279 // note that the FETools::get_fe_by_name function depends on the
280 // particular format of the string this function returns, so they have to be
281 // kept in synch
282 std::ostringstream namebuf;
283 namebuf << "FE_TraceQ<" << Utilities::dim_string(1, spacedim) << ">("
284 << this->degree << ")";
285
286 return namebuf.str();
287}
288
289
290
291// explicit instantiations
292#include "fe_trace.inst"
293
294
PolynomialType poly_space
FE_Q< dim, spacedim > fe_q
Definition fe_trace.h:140
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_trace.cc:75
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition fe_trace.cc:101
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition fe_trace.cc:120
virtual bool hp_constraints_are_implemented() const override
Definition fe_trace.cc:170
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition fe_trace.cc:178
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_trace.cc:132
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_trace.cc:232
virtual std::string get_name() const override
Definition fe_trace.cc:84
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition fe_trace.cc:152
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_trace.cc:217
FE_TraceQ(unsigned int p)
Definition fe_trace.cc:36
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition fe.h:2458
std::vector< Point< dim > > unit_support_points
Definition fe.h:2451
FullMatrix< double > interface_constraints
Definition fe.h:2439
size_type n() const
size_type m() 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)
#define AssertIndexRange(index, range)
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
static const unsigned int invalid_unsigned_int
Definition types.h:220
STL namespace.