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