Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_trace.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 #include <deal.II/base/config.h>
17 
18 #include <deal.II/base/qprojector.h>
19 #include <deal.II/base/quadrature_lib.h>
20 #include <deal.II/base/std_cxx14/memory.h>
21 #include <deal.II/base/tensor_product_polynomials.h>
22 
23 #include <deal.II/fe/fe_nothing.h>
24 #include <deal.II/fe/fe_poly_face.h>
25 #include <deal.II/fe/fe_q.h>
26 #include <deal.II/fe/fe_tools.h>
27 #include <deal.II/fe/fe_trace.h>
28 
29 #include <fstream>
30 #include <iostream>
31 #include <sstream>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 
36 
37 template <int dim, int spacedim>
38 FE_TraceQ<dim, spacedim>::FE_TraceQ(const unsigned int degree)
39  : FE_PolyFace<TensorProductPolynomials<dim - 1>, dim, spacedim>(
40  TensorProductPolynomials<dim - 1>(
41  Polynomials::generate_complete_Lagrange_basis(
42  QGaussLobatto<1>(degree + 1).get_points())),
43  FiniteElementData<dim>(get_dpo_vector(degree),
44  1,
45  degree,
46  FiniteElementData<dim>::L2),
47  std::vector<bool>(1, true))
48  , fe_q(degree)
49 {
50  Assert(degree > 0,
51  ExcMessage("FE_Trace can only be used for polynomial degrees "
52  "greater than zero"));
53  std::vector<unsigned int> renumber(this->dofs_per_face);
55  this->poly_space.set_numbering(renumber);
56 
57  // Initialize face support points
58  this->unit_face_support_points = fe_q.get_unit_face_support_points();
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->dofs_per_cell);
64  for (unsigned int i = 0; i < this->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 
73 template <int dim, int spacedim>
74 std::unique_ptr<FiniteElement<dim, spacedim>>
76 {
77  return std_cxx14::make_unique<FE_TraceQ<dim, spacedim>>(this->degree);
78 }
79 
80 
81 
82 template <int dim, int spacedim>
83 std::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 
99 template <int dim, int spacedim>
100 bool
102  const unsigned int shape_index,
103  const unsigned int face_index) const
104 {
105  Assert(shape_index < this->dofs_per_cell,
106  ExcIndexRange(shape_index, 0, this->dofs_per_cell));
109 
110  // FE_TraceQ shares the numbering of elemental degrees of freedom with FE_Q
111  // except for the missing interior ones (quad dofs in 2D and hex dofs in
112  // 3D). Therefore, it is safe to ask fe_q for the corresponding
113  // information. The assertion 'shape_index < this->dofs_per_cell' will make
114  // sure that we only access the trace dofs.
115  return fe_q.has_support_on_face(shape_index, face_index);
116 }
117 
118 
119 
120 template <int dim, int spacedim>
121 std::pair<Table<2, bool>, std::vector<unsigned int>>
123 {
124  Table<2, bool> constant_modes(1, this->dofs_per_cell);
125  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
126  constant_modes(0, i) = true;
127  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
128  constant_modes, std::vector<unsigned int>(1, 0));
129 }
130 
131 template <int dim, int spacedim>
132 void
135  const std::vector<Vector<double>> &support_point_values,
136  std::vector<double> & nodal_values) const
137 {
138  AssertDimension(support_point_values.size(),
139  this->get_unit_support_points().size());
140  AssertDimension(support_point_values.size(), nodal_values.size());
141  AssertDimension(this->dofs_per_cell, nodal_values.size());
142 
143  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
144  {
145  AssertDimension(support_point_values[i].size(), 1);
146 
147  nodal_values[i] = support_point_values[i](0);
148  }
149 }
150 
151 
152 template <int dim, int spacedim>
153 std::vector<unsigned int>
155 {
156  // This constructs FE_TraceQ in exactly the same way as FE_Q except for the
157  // interior degrees of freedom that are not present here (line in 1D, quad
158  // in 2D, hex in 3D).
159  AssertThrow(deg > 0, ExcMessage("FE_TraceQ needs to be of degree > 0."));
160  std::vector<unsigned int> dpo(dim + 1, 1U);
161  dpo[dim] = 0;
162  dpo[0] = 1;
163  for (unsigned int i = 1; i < dim; ++i)
164  dpo[i] = dpo[i - 1] * (deg - 1);
165  return dpo;
166 }
167 
168 
169 
170 template <int dim, int spacedim>
171 bool
173 {
174  return fe_q.hp_constraints_are_implemented();
175 }
176 
177 
178 template <int dim, int spacedim>
181  const FiniteElement<dim, spacedim> &fe_other,
182  const unsigned int codim) const
183 {
184  Assert(codim <= dim, ExcImpossibleInDim(dim));
185  (void)codim;
186 
187  // vertex/line/face/cell domination
188  // --------------------------------
189  if (const FE_TraceQ<dim, spacedim> *fe_traceq_other =
190  dynamic_cast<const FE_TraceQ<dim, spacedim> *>(&fe_other))
191  {
192  if (this->degree < fe_traceq_other->degree)
194  else if (this->degree == fe_traceq_other->degree)
196  else
198  }
199  else if (const FE_Nothing<dim> *fe_nothing =
200  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
201  {
202  if (fe_nothing->is_dominating())
204  else
205  // the FE_Nothing has no degrees of freedom and it is typically used
206  // in a context where we don't require any continuity along the
207  // interface
209  }
210 
211  Assert(false, ExcNotImplemented());
213 }
214 
215 
216 
217 template <int dim, int spacedim>
218 void
220  const FiniteElement<dim, spacedim> &source_fe,
221  FullMatrix<double> & interpolation_matrix) const
222 {
223  get_subface_interpolation_matrix(source_fe,
225  interpolation_matrix);
226 }
227 
228 
229 
230 template <int dim, int spacedim>
231 void
233  const FiniteElement<dim, spacedim> &x_source_fe,
234  const unsigned int subface,
235  FullMatrix<double> & interpolation_matrix) const
236 {
237  // this is the code from FE_FaceQ
238  Assert(interpolation_matrix.n() == this->dofs_per_face,
239  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
240  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
241  ExcDimensionMismatch(interpolation_matrix.m(),
242  x_source_fe.dofs_per_face));
243 
244  // see if source is a FaceQ element
245  if (const FE_TraceQ<dim, spacedim> *source_fe =
246  dynamic_cast<const FE_TraceQ<dim, spacedim> *>(&x_source_fe))
247  {
248  fe_q.get_subface_interpolation_matrix(source_fe->fe_q,
249  subface,
250  interpolation_matrix);
251  }
252  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
253  {
254  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
255  }
256  else
257  AssertThrow(
258  false,
259  (typename FiniteElement<dim,
260  spacedim>::ExcInterpolationNotImplemented()));
261 }
262 
263 
264 
265 template <int spacedim>
266 FE_TraceQ<1, spacedim>::FE_TraceQ(const unsigned int degree)
267  : FE_FaceQ<1, spacedim>(degree)
268 {}
269 
270 
271 
272 template <int spacedim>
273 std::string
275 {
276  // note that the FETools::get_fe_by_name function depends on the
277  // particular format of the string this function returns, so they have to be
278  // kept in synch
279  std::ostringstream namebuf;
280  namebuf << "FE_TraceQ<" << Utilities::dim_string(1, spacedim) << ">("
281  << this->degree << ")";
282 
283  return namebuf.str();
284 }
285 
286 
287 
288 // explicit instantiations
289 #include "fe_trace.inst"
290 
291 
292 DEAL_II_NAMESPACE_CLOSE
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_trace.cc:122
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
FullMatrix< double > interface_constraints
Definition: fe.h:2468
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_trace.cc:101
FE_TraceQ(unsigned int p)
Definition: fe_trace.cc:38
void hierarchic_to_lexicographic_numbering(unsigned int degree, std::vector< unsigned int > &h2l)
const unsigned int degree
Definition: fe_base.h:298
void set_numbering(const std::vector< unsigned int > &renumber)
std::vector< Point< dim - 1 > > unit_face_support_points
Definition: fe.h:2487
STL namespace.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2480
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_trace.cc:219
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:180
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
FE_Q< dim, spacedim > fe_q
Definition: fe_trace.h:138
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:134
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:457
const unsigned int dofs_per_cell
Definition: fe_base.h:282
virtual bool hp_constraints_are_implemented() const override
Definition: fe_trace.cc:172
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_trace.cc:75
const unsigned int dofs_per_face
Definition: fe_base.h:275
static ::ExceptionBase & ExcNotImplemented()
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition: fe_trace.cc:154
Definition: table.h:37
virtual std::string get_name() const override
Definition: fe_trace.cc:84
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_trace.cc:232