Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00: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\}}\)
fe_q.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2023 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 
18 
19 #include <deal.II/fe/fe_dgq.h>
20 #include <deal.II/fe/fe_hermite.h>
21 #include <deal.II/fe/fe_nothing.h>
23 #include <deal.II/fe/fe_q.h>
26 #include <deal.II/fe/fe_wedge_p.h>
27 
28 #include <deal.II/lac/vector.h>
29 
30 #include <memory>
31 #include <sstream>
32 #include <vector>
33 
35 
36 
37 namespace internal
38 {
39  namespace FE_Q
40  {
41  namespace
42  {
43  std::vector<Point<1>>
44  get_QGaussLobatto_points(const unsigned int degree)
45  {
46  if (degree > 0)
47  return QGaussLobatto<1>(degree + 1).get_points();
48  else
49  {
50  using FEQ = ::FE_Q_Base<1, 1>;
51  AssertThrow(false, FEQ::ExcFEQCannotHaveDegree0());
52  }
53  return std::vector<Point<1>>();
54  }
55  } // namespace
56  } // namespace FE_Q
57 } // namespace internal
58 
59 
60 
61 template <int dim, int spacedim>
62 FE_Q<dim, spacedim>::FE_Q(const unsigned int degree)
63  : FE_Q_Base<dim, spacedim>(
66  internal::FE_Q::get_QGaussLobatto_points(degree))),
67  FiniteElementData<dim>(this->get_dpo_vector(degree),
68  1,
69  degree,
70  FiniteElementData<dim>::H1),
71  std::vector<bool>(1, false))
72 {
73  this->initialize(internal::FE_Q::get_QGaussLobatto_points(degree));
74 }
75 
76 
77 
78 template <int dim, int spacedim>
80  : FE_Q_Base<dim, spacedim>(
82  Polynomials::generate_complete_Lagrange_basis(points.get_points())),
83  FiniteElementData<dim>(this->get_dpo_vector(points.size() - 1),
84  1,
85  points.size() - 1,
86  FiniteElementData<dim>::H1),
87  std::vector<bool>(1, false))
88 {
89  this->initialize(points.get_points());
90 }
91 
92 
93 
94 template <int dim, int spacedim>
95 std::string
97 {
98  // note that the FETools::get_fe_by_name function depends on the
99  // particular format of the string this function returns, so they have to be
100  // kept in synch
101 
102  std::ostringstream namebuf;
103  bool equidistant = true;
104  std::vector<double> points(this->degree + 1);
105 
106  // Decode the support points in one coordinate direction.
107  TensorProductPolynomials<dim> *poly_space_derived_ptr =
108  dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
109  std::vector<unsigned int> lexicographic =
110  poly_space_derived_ptr->get_numbering_inverse();
111  for (unsigned int j = 0; j <= this->degree; ++j)
112  points[j] = this->unit_support_points[lexicographic[j]][0];
113 
114  // Check whether the support points are equidistant.
115  for (unsigned int j = 0; j <= this->degree; ++j)
116  if (std::fabs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
117  {
118  equidistant = false;
119  break;
120  }
121 
122  if (equidistant == true)
123  {
124  if (this->degree > 2)
125  namebuf << "FE_Q<" << Utilities::dim_string(dim, spacedim)
126  << ">(QIterated(QTrapezoid()," << this->degree << "))";
127  else
128  namebuf << "FE_Q<" << Utilities::dim_string(dim, spacedim) << ">("
129  << this->degree << ")";
130  }
131  else
132  {
133  // Check whether the support points come from QGaussLobatto.
134  const QGaussLobatto<1> points_gl(this->degree + 1);
135  bool gauss_lobatto = true;
136  for (unsigned int j = 0; j <= this->degree; ++j)
137  if (points[j] != points_gl.point(j)(0))
138  {
139  gauss_lobatto = false;
140  break;
141  }
142  if (gauss_lobatto == true)
143  namebuf << "FE_Q<" << Utilities::dim_string(dim, spacedim) << ">("
144  << this->degree << ")";
145  else
146  namebuf << "FE_Q<" << Utilities::dim_string(dim, spacedim)
147  << ">(QUnknownNodes(" << this->degree << "))";
148  }
149  return namebuf.str();
150 }
151 
152 
153 
154 template <int dim, int spacedim>
155 void
157  const std::vector<Vector<double>> &support_point_values,
158  std::vector<double> &nodal_values) const
159 {
160  AssertDimension(support_point_values.size(),
161  this->get_unit_support_points().size());
162  AssertDimension(support_point_values.size(), nodal_values.size());
163  AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
164 
165  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
166  {
167  AssertDimension(support_point_values[i].size(), 1);
168 
169  nodal_values[i] = support_point_values[i](0);
170  }
171 }
172 
173 
174 
175 template <int dim, int spacedim>
176 std::unique_ptr<FiniteElement<dim, spacedim>>
178 {
179  return std::make_unique<FE_Q<dim, spacedim>>(*this);
180 }
181 
182 
183 
184 template <int dim, int spacedim>
187  const FiniteElement<dim, spacedim> &fe_other,
188  const unsigned int codim) const
189 {
190  Assert(codim <= dim, ExcImpossibleInDim(dim));
191 
192  // vertex/line/face domination
193  // (if fe_other is derived from FE_DGQ)
194  // ------------------------------------
195  if (codim > 0)
196  if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
197  // there are no requirements between continuous and discontinuous elements
199 
200  // vertex/line/face domination
201  // (if fe_other is not derived from FE_DGQ)
202  // & cell domination
203  // ----------------------------------------
204  if (const FE_Q<dim, spacedim> *fe_q_other =
205  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
206  {
207  if (this->degree < fe_q_other->degree)
209  else if (this->degree == fe_q_other->degree)
211  else
213  }
214  else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
215  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
216  {
217  if (this->degree < fe_p_other->degree)
219  else if (this->degree == fe_p_other->degree)
221  else
223  }
224  else if (const FE_WedgeP<dim, spacedim> *fe_wp_other =
225  dynamic_cast<const FE_WedgeP<dim, spacedim> *>(&fe_other))
226  {
227  if (this->degree < fe_wp_other->degree)
229  else if (this->degree == fe_wp_other->degree)
231  else
233  }
234  else if (const FE_PyramidP<dim, spacedim> *fe_pp_other =
235  dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other))
236  {
237  if (this->degree < fe_pp_other->degree)
239  else if (this->degree == fe_pp_other->degree)
241  else
243  }
244  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
245  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
246  {
247  if (fe_nothing->is_dominating())
249  else
250  // the FE_Nothing has no degrees of freedom and it is typically used
251  // in a context where we don't require any continuity along the
252  // interface
254  }
255  else if (const FE_Hermite<dim, spacedim> *fe_hermite_other =
256  dynamic_cast<const FE_Hermite<dim, spacedim> *>(&fe_other))
257  {
258  if (this->degree == 1)
259  {
260  if (fe_hermite_other->degree > 1)
262  else
264  }
265  else if (this->degree >= fe_hermite_other->degree)
267  else
269  }
270 
271  Assert(false, ExcNotImplemented());
273 }
274 
275 
276 // explicit instantiations
277 #include "fe_q.inst"
278 
Definition: fe_dgq.h:113
void initialize(const std::vector< Point< 1 >> &support_points_1d)
Definition: fe_q_base.cc:437
Definition: fe_q.h:551
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.cc:156
FE_Q(const unsigned int p)
Definition: fe_q.cc:62
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_q.cc:177
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_q.cc:186
virtual std::string get_name() const override
Definition: fe_q.cc:96
const unsigned int degree
Definition: fe_data.h:453
const std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) const
const std::vector< unsigned int > & get_numbering_inverse() const
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Expression fabs(const Expression &x)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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:556
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)