deal.II version GIT relicensing-2076-g6b43d56a25 2024-11-02 12:30:00+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
util.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 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
16#ifndef dealii_matrix_free_util_h
17#define dealii_matrix_free_util_h
18
19
20#include <deal.II/base/config.h>
21
24
26
28
30
32
33
34namespace internal
35{
36 namespace MatrixFreeFunctions
37 {
43 template <int dim>
44 inline std::pair<ReferenceCell, ::hp::QCollection<dim - 1>>
46 const bool do_assert = true)
47 {
48 if (dim == 2 || dim == 3)
49 {
50 for (unsigned int i = 1; i <= 4; ++i)
51 if (quad == QGaussSimplex<dim>(i))
52 return {ReferenceCells::get_simplex<dim>(),
53 ::hp::QCollection<dim - 1>(
54 QGaussSimplex<dim - 1>(i))};
55
56 for (unsigned int i = 1; i <= 5; ++i)
57 if (quad == QWitherdenVincentSimplex<dim>(i))
58 return {ReferenceCells::get_simplex<dim>(),
59 ::hp::QCollection<dim - 1>(
60 QWitherdenVincentSimplex<dim - 1>(i))};
61
62 for (unsigned int i = 1; i <= 3; ++i)
63 {
64 const FE_SimplexP<dim> fe(i);
65 if (quad == Quadrature<dim>(fe.get_unit_support_points()))
66 return {ReferenceCells::get_simplex<dim>(),
67 ::hp::QCollection<dim - 1>(Quadrature<dim - 1>(
68 fe.get_unit_face_support_points()))};
69 }
70 }
71
72 if (dim == 3)
73 for (unsigned int i = 1; i <= 3; ++i)
74 if (quad == QGaussWedge<dim>(i))
75 {
76 QGauss<dim - 1> quad(i);
77 QGaussSimplex<dim - 1> tri(i);
78
79 return {
81 ::hp::QCollection<dim - 1>(tri, tri, quad, quad, quad)};
82 }
83
84 if (dim == 3)
85 for (unsigned int i = 1; i <= 2; ++i)
86 if (quad == QGaussPyramid<dim>(i))
87 {
88 QGauss<dim - 1> quad(i);
89 QGaussSimplex<dim - 1> tri(i);
90
91 return {
93 ::hp::QCollection<dim - 1>(quad, tri, tri, tri, tri)};
94 }
95
96 // note: handle hypercubes last since normally this function is not
97 // called for hypercubes
98 for (unsigned int i = 1; i <= 5; ++i)
99 if (quad == QGauss<dim>(i))
100 return {ReferenceCells::get_hypercube<dim>(),
101 ::hp::QCollection<dim - 1>(QGauss<dim - 1>(i))};
102
103 if (do_assert)
105
106 return {ReferenceCells::Invalid, ::hp::QCollection<dim - 1>()};
107 }
108
109
110
119 template <int dim>
120 inline std::pair<Quadrature<dim - 1>, Quadrature<dim - 1>>
122 {
124 quad.size() > 0,
126 "There is nothing useful you can do with a MatrixFree/FEEvaluation "
127 "object when using a quadrature formula with zero "
128 "quadrature points!"));
129
130 if (dim == 2 || dim == 3)
131 {
132 for (unsigned int i = 1; i <= 4; ++i)
133 if (quad == QGaussSimplex<dim>(i))
134 {
135 if (dim == 2)
136 return {QGaussSimplex<dim - 1>(i), // line!
138 else
139 return {Quadrature<dim - 1>(), QGaussSimplex<dim - 1>(i)};
140 }
141
142 for (unsigned int i = 1; i <= 5; ++i)
143 if (quad == QWitherdenVincentSimplex<dim>(i))
144 {
145 if (dim == 2)
146 return {QWitherdenVincentSimplex<dim - 1>(i), // line!
148 else
149 return {Quadrature<dim - 1>(),
151 }
152
153 for (unsigned int i = 1; i <= 3; ++i)
154 {
155 const FE_SimplexP<dim> fe(i);
156 if (quad == Quadrature<dim>(fe.get_unit_support_points()))
157 {
158 if (dim == 2)
159 return {Quadrature<dim - 1>(
160 fe.get_unit_face_support_points()), // line!
162 else
163 return {Quadrature<dim - 1>(),
165 fe.get_unit_face_support_points())};
166 }
167 }
168 }
169
170 if (dim == 3)
171 for (unsigned int i = 1; i <= 3; ++i)
172 if (quad == QGaussWedge<dim>(i))
173 return {QGauss<dim - 1>(i), QGaussSimplex<dim - 1>(i)};
174
175 if (dim == 3)
176 for (unsigned int i = 1; i <= 2; ++i)
177 if (quad == QGaussPyramid<dim>(i))
178 return {QGauss<dim - 1>(i), QGaussSimplex<dim - 1>(i)};
179
180 // note: handle hypercubes last since normally this function is not
181 // called for hypercubes
182 for (unsigned int i = 1; i <= 5; ++i)
183 if (quad == QGauss<dim>(i))
184 return {QGauss<dim - 1>(i), Quadrature<dim - 1>()};
185
187
188 return {Quadrature<dim - 1>(), Quadrature<dim - 1>()};
189 }
190
191 inline DEAL_II_ALWAYS_INLINE unsigned int
192 indicate_power_of_two(const unsigned int vectorization_length)
193 {
194 unsigned int vectorization_length_bits = 0;
195 unsigned int my_length = vectorization_length;
196 while (my_length >>= 1)
198 return 1 << vectorization_length_bits;
199 }
200
201 } // end of namespace MatrixFreeFunctions
202} // end of namespace internal
203
205
206#endif
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Invalid
unsigned int indicate_power_of_two(const unsigned int vectorization_length)
Definition util.h:192
std::pair< Quadrature< dim - 1 >, Quadrature< dim - 1 > > get_unique_face_quadratures(const Quadrature< dim > &quad)
Definition util.h:121
std::pair< ReferenceCell, ::hp::QCollection< dim - 1 > > get_face_quadrature_collection(const Quadrature< dim > &quad, const bool do_assert=true)
Definition util.h:45