DUNE-FEM (unstable)

axisalignedcubegeometry.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5
6#ifndef DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
7#define DUNE_GEOMETRY_AXISALIGNED_CUBE_GEOMETRY_HH
8
13#include <bitset>
14
18
19#include <dune/geometry/referenceelements.hh>
20#include <dune/geometry/type.hh>
21
22
23namespace Dune {
24
48 template <class CoordType, unsigned int dim, unsigned int coorddim>
50 {
51
52
53 public:
54
56 constexpr static int mydimension = dim;
57
59 constexpr static int coorddimension = coorddim;
60
62 typedef CoordType ctype;
63
66
69
71 typedef ctype Volume;
72
79 typedef typename std::conditional<dim==coorddim,
82
89 typedef typename std::conditional<dim==coorddim,
92
100 using Jacobian = std::conditional_t<dim==coorddim, DiagonalMatrix<ctype,dim>, FieldMatrix<ctype,coorddim,dim> >;
101
109 using JacobianInverse = std::conditional_t<dim==coorddim, DiagonalMatrix<ctype,dim>, FieldMatrix<ctype,dim,coorddim> >;
110
119
126 : lower_(lower),
127 upper_(upper),
128 axes_()
129 {
130 static_assert(dim==coorddim, "Use this constructor only if dim==coorddim!");
131 // all 'true', but is never actually used
132 axes_ = (1<<coorddim)-1;
133 }
134
144 const std::bitset<coorddim>& axes)
145 : lower_(lower),
146 upper_(upper),
147 axes_(axes)
148 {
149 assert(axes.count()==dim);
150 for (size_t i=0; i<coorddim; i++)
151 if (not axes_[i])
152 upper_[i] = lower_[i];
153 }
154
160 : lower_(lower)
161 {}
162
165 {
166 return GeometryTypes::cube(dim);
167 }
168
171 {
172 GlobalCoordinate result;
173 if (dim == coorddim) { // fast case
174 for (size_t i=0; i<coorddim; i++)
175 result[i] = lower_[i] + local[i]*(upper_[i] - lower_[i]);
176 } else if (dim == 0) { // a vertex -- the other fast case
177 result = lower_; // hope for named-return-type-optimization
178 } else { // slow case
179 size_t lc=0;
180 for (size_t i=0; i<coorddim; i++)
181 result[i] = (axes_[i])
182 ? lower_[i] + local[lc++]*(upper_[i] - lower_[i])
183 : lower_[i];
184 }
185 return result;
186 }
187
190 {
191 LocalCoordinate result;
192 if (dim == coorddim) { // fast case
193 for (size_t i=0; i<dim; i++)
194 result[i] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
195 } else if (dim != 0) { // slow case
196 size_t lc=0;
197 for (size_t i=0; i<coorddim; i++)
198 if (axes_[i])
199 result[lc++] = (global[i] - lower_[i]) / (upper_[i] - lower_[i]);
200 }
201 return result;
202 }
203
206 {
207 JacobianTransposed result;
208
209 // Actually compute the result. Uses different methods depending
210 // on what kind of matrix JacobianTransposed is.
211 jacobianTransposed(result);
212
213 return result;
214 }
215
218 {
220
221 // Actually compute the result. Uses different methods depending
222 // on what kind of matrix JacobianTransposed is.
224
225 return result;
226 }
227
229 Jacobian jacobian([[maybe_unused]] const LocalCoordinate& local) const
230 {
231 return jacobianTransposed(local).transposed();
232 }
233
236 {
237 return jacobianInverseTransposed(local).transposed();
238 }
239
243 Volume integrationElement([[maybe_unused]] const LocalCoordinate& local) const
244 {
245 return volume();
246 }
247
250 {
251 GlobalCoordinate result;
252 if (dim==0)
253 result = lower_;
254 else {
255 // Since lower_==upper_ for unused coordinates, this always does the right thing
256 for (size_t i=0; i<coorddim; i++)
257 result[i] = CoordType(0.5) * (lower_[i] + upper_[i]);
258 }
259 return result;
260 }
261
263 int corners() const
264 {
265 return 1<<dim;
266 }
267
270 {
271 GlobalCoordinate result;
272 if (dim == coorddim) { // fast case
273 for (size_t i=0; i<coorddim; i++)
274 result[i] = (k & (1<<i)) ? upper_[i] : lower_[i];
275 } else if (dim == 0) { // vertex
276 result = lower_; // rely on named return-type optimization
277 } else { // slow case
278 unsigned int mask = 1;
279
280 for (size_t i=0; i<coorddim; i++) {
281 if (not axes_[i])
282 result[i] = lower_[i];
283 else {
284 result[i] = (k & mask) ? upper_[i] : lower_[i];
285 mask = (mask<<1);
286 }
287 }
288 }
289
290
291 return result;
292 }
293
296 {
297 ctype vol = 1;
298 if (dim == coorddim) { // fast case
299 for (size_t i=0; i<dim; i++)
300 vol *= upper_[i] - lower_[i];
301 // do nothing if dim == 0
302 } else if (dim != 0) { // slow case
303 for (size_t i=0; i<coorddim; i++)
304 if (axes_[i])
305 vol *= upper_[i] - lower_[i];
306 }
307 return vol;
308 }
309
311 bool affine() const
312 {
313 return true;
314 }
315
316 friend typename Dune::ReferenceElements< ctype, dim >::ReferenceElement referenceElement ( const AxisAlignedCubeGeometry & /* geometry */ )
317 {
319 }
320
321 private:
322 // jacobianTransposed: fast case --> diagonal matrix
323 void jacobianTransposed ( DiagonalMatrix<ctype,dim> &jacobianTransposed ) const
324 {
325 for (size_t i=0; i<dim; i++)
326 jacobianTransposed.diagonal()[i] = upper_[i] - lower_[i];
327 }
328
329 // jacobianTransposed: slow case --> dense matrix
330 void jacobianTransposed ( FieldMatrix<ctype,dim,coorddim> &jacobianTransposed ) const
331 {
332 if (dim==0)
333 return;
334
335 size_t lc = 0;
336 for (size_t i=0; i<coorddim; i++)
337 if (axes_[i])
338 jacobianTransposed[lc++][i] = upper_[i] - lower_[i];
339 }
340
341 // jacobianInverseTransposed: fast case --> diagonal matrix
342 void jacobianInverseTransposed ( DiagonalMatrix<ctype,dim> &jacobianInverseTransposed ) const
343 {
344 for (size_t i=0; i<dim; i++)
345 jacobianInverseTransposed.diagonal()[i] = CoordType(1.0) / (upper_[i] - lower_[i]);
346 }
347
348 // jacobianInverseTransposed: slow case --> dense matrix
349 void jacobianInverseTransposed ( FieldMatrix<ctype,coorddim,dim> &jacobianInverseTransposed ) const
350 {
351 if (dim==0)
352 return;
353
354 size_t lc = 0;
355 for (size_t i=0; i<coorddim; i++)
356 if (axes_[i])
357 jacobianInverseTransposed[i][lc++] = CoordType(1.0) / (upper_[i] - lower_[i]);
358 }
359
361
363
364 std::bitset<coorddim> axes_;
365 };
366
367} // namespace Dune
368#endif
A geometry implementation for axis-aligned hypercubes.
Definition: axisalignedcubegeometry.hh:50
Volume volume() const
Return the element volume.
Definition: axisalignedcubegeometry.hh:295
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower, const Dune::FieldVector< ctype, coorddim > upper, const std::bitset< coorddim > &axes)
Constructor from a lower left and an upper right corner.
Definition: axisalignedcubegeometry.hh:142
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower, const Dune::FieldVector< ctype, coorddim > upper)
Constructor from a lower left and an upper right corner.
Definition: axisalignedcubegeometry.hh:124
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Inverse Jacobian of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:235
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Inverse Jacobian transposed of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:217
std::conditional< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, dim, coorddim > >::type JacobianTransposed
Return type of jacobianTransposed.
Definition: axisalignedcubegeometry.hh:81
AxisAlignedCubeGeometry(const Dune::FieldVector< ctype, coorddim > lower)
Constructor from a single point only.
Definition: axisalignedcubegeometry.hh:159
static constexpr int mydimension
Dimension of the cube element.
Definition: axisalignedcubegeometry.hh:56
static constexpr int coorddimension
Dimension of the world space that the cube element is embedded in.
Definition: axisalignedcubegeometry.hh:59
GlobalCoordinate corner(int k) const
Return world coordinates of the k-th corner of the element.
Definition: axisalignedcubegeometry.hh:269
ctype Volume
Type used for volume.
Definition: axisalignedcubegeometry.hh:71
FieldVector< ctype, dim > LocalCoordinate
Type used for a vector of element coordinates.
Definition: axisalignedcubegeometry.hh:65
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
Jacobian transposed of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:205
FieldVector< ctype, coorddim > GlobalCoordinate
Type used for a vector of world coordinates.
Definition: axisalignedcubegeometry.hh:68
LocalCoordinate local(const GlobalCoordinate &global) const
Map a point in global (world) coordinates to element coordinates.
Definition: axisalignedcubegeometry.hh:189
CoordType ctype
Type used for single coordinate coefficients.
Definition: axisalignedcubegeometry.hh:62
std::conditional_t< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, coorddim, dim > > Jacobian
Return type of jacobian.
Definition: axisalignedcubegeometry.hh:100
GeometryType type() const
Type of the cube. Here: a hypercube of the correct dimension.
Definition: axisalignedcubegeometry.hh:164
int corners() const
Return the number of corners of the element.
Definition: axisalignedcubegeometry.hh:263
Jacobian jacobian(const LocalCoordinate &local) const
Jacobian of the transformation from local to global coordinates.
Definition: axisalignedcubegeometry.hh:229
std::conditional< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, coorddim, dim > >::type JacobianInverseTransposed
Return type of jacobianInverseTransposed.
Definition: axisalignedcubegeometry.hh:91
Volume integrationElement(const LocalCoordinate &local) const
Return the integration element, i.e., the determinant term in the integral transformation formula.
Definition: axisalignedcubegeometry.hh:243
GlobalCoordinate center() const
Return center of mass of the element.
Definition: axisalignedcubegeometry.hh:249
AxisAlignedCubeGeometry()=default
Constructs an empty geometry.
bool affine() const
Return if the element is affine. Here: yes.
Definition: axisalignedcubegeometry.hh:311
std::conditional_t< dim==coorddim, DiagonalMatrix< ctype, dim >, FieldMatrix< ctype, dim, coorddim > > JacobianInverse
Return type of jacobianInverse.
Definition: axisalignedcubegeometry.hh:109
GlobalCoordinate global(const LocalCoordinate &local) const
Map a point in local (element) coordinates to world coordinates.
Definition: axisalignedcubegeometry.hh:170
A diagonal matrix of static size.
Definition: diagonalmatrix.hh:53
A dense n x m matrix.
Definition: fmatrix.hh:117
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
This file implements a quadratic diagonal matrix of fixed size.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
Mask< V > mask(ADLTag< 0, std::is_same< V, Mask< V > >::value >, const V &v)
implements Simd::mask()
Definition: defaults.hh:153
Dune namespace.
Definition: alignedallocator.hh:13
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition: referenceelements.hh:146
static const ReferenceElement & cube()
get hypercube reference elements
Definition: referenceelements.hh:168
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)