Reference documentation for deal.II version GIT c9976103bc 2022-12-09 17:30: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_wedge_p.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2022 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_dgq.h>
22 #include <deal.II/fe/fe_nothing.h>
24 #include <deal.II/fe/fe_q.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/fe_wedge_p.h>
29 
31 
32 namespace
33 {
38  get_dpo_vector_fe_wedge_p(const unsigned int degree)
39  {
41 
42  if (degree == 1)
43  {
44  dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}};
45  dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}};
46  dpo.object_index = {{}, {6}, {6}, {6}};
47  dpo.first_object_index_on_face = {{}, {3, 3, 4, 4, 4}, {3, 3, 4, 4, 4}};
48  }
49  else if (degree == 2)
50  {
51  dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}};
52  dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}};
53  dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}};
54  dpo.first_object_index_on_face = {{}, {3, 3, 4, 4, 4}, {6, 6, 8, 8, 8}};
55  }
56  else
57  {
58  Assert(false, ExcNotImplemented());
59  }
60 
61  return dpo;
62  }
63 
68  get_dpo_vector_fe_wedge_dgp(const unsigned int degree)
69  {
70  unsigned int n_dofs = 0;
71 
72  if (degree == 1)
73  n_dofs = 6;
74  else if (degree == 2)
75  n_dofs = 18;
76  else
77  Assert(false, ExcNotImplemented());
78 
79  return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Wedge);
80  }
81 } // namespace
82 
83 
84 
85 template <int dim, int spacedim>
87  const unsigned int degree,
88  const internal::GenericDoFsPerObject & dpos,
89  const typename FiniteElementData<dim>::Conformity conformity)
90  : ::FE_Poly<dim, spacedim>(
91  ScalarLagrangePolynomialWedge<dim>(degree),
92  FiniteElementData<dim>(dpos,
94  1,
95  degree,
96  conformity),
97  std::vector<bool>(
98  FiniteElementData<dim>(dpos, ReferenceCells::Wedge, 1, degree)
99  .dofs_per_cell,
100  true),
101  std::vector<ComponentMask>(
102  FiniteElementData<dim>(dpos, ReferenceCells::Wedge, 1, degree)
103  .dofs_per_cell,
104  std::vector<bool>(1, true)))
105 {
106  AssertDimension(dim, 3);
107 
108  Assert(1 <= degree && degree <= 2, ExcNotImplemented());
109 
110  FE_SimplexP<2> fe_triangle(degree);
111  FE_Q<1> fe_line(degree);
112  FE_Q<2> fe_quad(degree);
113 
114  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
115  {
116  const auto pair = this->degree == 1 ? internal::wedge_table_1[i] :
118 
119  this->unit_support_points.emplace_back(
120  fe_triangle.get_unit_support_points()[pair[0]][0],
121  fe_triangle.get_unit_support_points()[pair[0]][1],
122  fe_line.get_unit_support_points()[pair[1]][0]);
123  }
124 
125  this->unit_face_support_points.resize(this->reference_cell().n_faces());
126 
127  for (const auto f : this->reference_cell().face_indices())
128  if (this->reference_cell().face_reference_cell(f) ==
130  for (const auto &p : fe_triangle.get_unit_support_points())
131  this->unit_face_support_points[f].emplace_back(p[0], p[1]);
132  else if (this->reference_cell().face_reference_cell(f) ==
134  for (const auto &p : fe_quad.get_unit_support_points())
135  this->unit_face_support_points[f].emplace_back(p[0], p[1]);
136  else
137  Assert(false, ExcInternalError());
138 }
139 
140 
141 
142 template <int dim, int spacedim>
143 void
146  const std::vector<Vector<double>> &support_point_values,
147  std::vector<double> & nodal_values) const
148 {
149  AssertDimension(support_point_values.size(),
150  this->get_unit_support_points().size());
151  AssertDimension(support_point_values.size(), nodal_values.size());
152  AssertDimension(this->dofs_per_cell, nodal_values.size());
153 
154  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
155  {
156  AssertDimension(support_point_values[i].size(), 1);
157 
158  nodal_values[i] = support_point_values[i](0);
159  }
160 }
161 
162 
163 
164 template <int dim, int spacedim>
165 FE_WedgeP<dim, spacedim>::FE_WedgeP(const unsigned int degree)
166  : FE_WedgePoly<dim, spacedim>(degree,
167  get_dpo_vector_fe_wedge_p(degree),
168  FiniteElementData<dim>::H1)
169 {}
170 
171 
172 
173 template <int dim, int spacedim>
174 std::unique_ptr<FiniteElement<dim, spacedim>>
176 {
177  return std::make_unique<FE_WedgeP<dim, spacedim>>(*this);
178 }
179 
180 
181 
182 template <int dim, int spacedim>
183 std::string
185 {
186  std::ostringstream namebuf;
187  namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")";
188 
189  return namebuf.str();
190 }
191 
192 
193 
194 template <int dim, int spacedim>
197  const FiniteElement<dim, spacedim> &fe_other,
198  const unsigned int codim) const
199 {
200  Assert(codim <= dim, ExcImpossibleInDim(dim));
201 
202  // vertex/line/face domination
203  // (if fe_other is derived from FE_SimplexDGP)
204  // ------------------------------------
205  if (codim > 0)
206  if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
207  nullptr)
208  // there are no requirements between continuous and discontinuous
209  // elements
211 
212 
213  // vertex/line/face domination
214  // (if fe_other is not derived from FE_SimplexDGP)
215  // & cell domination
216  // ----------------------------------------
217  if (const FE_WedgeP<dim, spacedim> *fe_wp_other =
218  dynamic_cast<const FE_WedgeP<dim, spacedim> *>(&fe_other))
219  {
220  if (this->degree < fe_wp_other->degree)
222  else if (this->degree == fe_wp_other->degree)
224  else
226  }
227  else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
228  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
229  {
230  if (this->degree < fe_p_other->degree)
232  else if (this->degree == fe_p_other->degree)
234  else
236  }
237  else if (const FE_Q<dim, spacedim> *fe_q_other =
238  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
239  {
240  if (this->degree < fe_q_other->degree)
242  else if (this->degree == fe_q_other->degree)
244  else
246  }
247  else if (const FE_Nothing<dim> *fe_nothing =
248  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
249  {
250  if (fe_nothing->is_dominating())
252  else
253  // the FE_Nothing has no degrees of freedom and it is typically used
254  // in a context where we don't require any continuity along the
255  // interface
257  }
258 
259  Assert(false, ExcNotImplemented());
261 }
262 
263 
264 
265 template <int dim, int spacedim>
266 std::vector<std::pair<unsigned int, unsigned int>>
268  const FiniteElement<dim, spacedim> &fe_other) const
269 {
270  (void)fe_other;
271 
272  Assert((dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other)) ||
273  (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
275 
276  return {{0, 0}};
277 }
278 
279 
280 
281 template <int dim, int spacedim>
282 std::vector<std::pair<unsigned int, unsigned int>>
284  const FiniteElement<dim, spacedim> &fe_other) const
285 {
286  (void)fe_other;
287 
288  Assert((dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other)) ||
289  (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
291 
292  std::vector<std::pair<unsigned int, unsigned int>> result;
293 
294  for (unsigned int i = 0; i < this->degree - 1; ++i)
295  result.emplace_back(i, i);
296 
297  return result;
298 }
299 
300 
301 
302 template <int dim, int spacedim>
303 std::vector<std::pair<unsigned int, unsigned int>>
305  const FiniteElement<dim, spacedim> &fe_other,
306  const unsigned int face_no) const
307 {
308  (void)fe_other;
309 
310  AssertIndexRange(face_no, 5);
311 
312  if (face_no < 2)
313  {
314  Assert((dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other)),
316  }
317  else
318  {
319  Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
321  }
322 
323  std::vector<std::pair<unsigned int, unsigned int>> result;
324 
325  for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i)
326  result.emplace_back(i, i);
327 
328  return result;
329 }
330 
331 
332 
333 template <int dim, int spacedim>
334 FE_WedgeDGP<dim, spacedim>::FE_WedgeDGP(const unsigned int degree)
335  : FE_WedgePoly<dim, spacedim>(degree,
336  get_dpo_vector_fe_wedge_dgp(degree),
337  FiniteElementData<dim>::L2)
338 {}
339 
340 
341 
342 template <int dim, int spacedim>
343 std::unique_ptr<FiniteElement<dim, spacedim>>
345 {
346  return std::make_unique<FE_WedgeDGP<dim, spacedim>>(*this);
347 }
348 
349 
350 
351 template <int dim, int spacedim>
352 std::string
354 {
355  std::ostringstream namebuf;
356  namebuf << "FE_WedgeDGP<" << dim << ">(" << this->degree << ")";
357 
358  return namebuf.str();
359 }
360 
361 // explicit instantiations
362 #include "fe_wedge_p.inst"
363 
Definition: fe_q.h:551
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_wedge_p.cc:344
std::string get_name() const override
Definition: fe_wedge_p.cc:353
FE_WedgeDGP(const unsigned int degree)
Definition: fe_wedge_p.cc:334
std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const override
Definition: fe_wedge_p.cc:304
FE_WedgeP(const unsigned int degree)
Definition: fe_wedge_p.cc:165
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_wedge_p.cc:267
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_wedge_p.cc:283
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
Definition: fe_wedge_p.cc:196
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_wedge_p.cc:175
std::string get_name() const override
Definition: fe_wedge_p.cc:184
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_wedge_p.cc:145
FE_WedgePoly(const unsigned int degree, const internal::GenericDoFsPerObject &dpos, const typename FiniteElementData< dim >::Conformity conformity)
Definition: fe_wedge_p.cc:86
const unsigned int degree
Definition: fe_data.h:449
unsigned int n_dofs_per_cell() const
ReferenceCell reference_cell() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition: fe.h:2446
const std::vector< Point< dim > > & get_unit_support_points() const
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2439
unsigned int n_faces() const
ReferenceCell face_reference_cell(const unsigned int face_no) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices() const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Triangle
static constexpr const ::ndarray< unsigned int, 18, 2 > wedge_table_2
internal::GenericDoFsPerObject expand(const unsigned int dim, const std::vector< unsigned int > &dofs_per_object, const ::ReferenceCell reference_cell)
Definition: fe_data.cc:25
static constexpr const ::ndarray< unsigned int, 6, 2 > wedge_table_1
std::vector< std::vector< unsigned int > > object_index
Definition: fe_data.h:195
std::vector< std::vector< unsigned int > > first_object_index_on_face
Definition: fe_data.h:200
std::vector< std::vector< unsigned int > > dofs_per_object_inclusive
Definition: fe_data.h:190
std::vector< std::vector< unsigned int > > dofs_per_object_exclusive
Definition: fe_data.h:185