5
7
* \brief The UGGridGeometry class and its specializations
8
#include <dune/common/array.hh>
9
10
#include <dune/common/fmatrix.hh>
11
#include <dune/geometry/genericgeometry/geometry.hh>
13
#include "uggridrenumberer.hh"
12
#include <dune/geometry/multilineargeometry.hh>
18
/** \brief Defines the geometry part of a mesh entity.
17
/** \brief Defines the geometry part of a mesh entity.
21
\tparam mydim Dimension of the corresponding reference element
22
\tparam coorddim Each corner is a point with coorddim coordinates.
20
\tparam mydim Dimension of the corresponding reference element
21
\tparam coorddim Each corner is a point with coorddim coordinates.
24
This version is actually used only for mydim==coorddim. The mydim<coorddim cases
25
are in specializations below.
27
template<int mydim, int coorddim, class GridImp>
28
class UGGridGeometry :
29
public GeometryDefaultImplementation <mydim, coorddim, GridImp, UGGridGeometry>
23
This version is actually used only for mydim==coorddim. The mydim<coorddim cases
24
are in specializations below.
26
template<int mydim, int coorddim, class GridImp>
27
class UGGridGeometry :
28
public GeometryDefaultImplementation <mydim, coorddim, GridImp, UGGridGeometry>
31
30
typedef typename GridImp::ctype UGCtype;
33
32
template <int codim_, int dim_, class GridImp_>
34
33
friend class UGGridEntity;
38
37
/** \brief Default constructor
42
jacobianIsUpToDate_ = false;
43
jacobianInverseIsUpToDate_ = false;
41
jacobianIsUpToDate_ = false;
42
jacobianInverseIsUpToDate_ = false;
46
/** \brief Return the element type identifier
45
/** \brief Return the element type identifier
48
47
* UGGrid supports triangles and quadrilaterals in 2D, and
49
48
* tetrahedra, pyramids, prisms, and hexahedra in 3D.
56
55
//! return the number of corners of this element.
57
56
int corners () const {
58
return UG_NS<coorddim>::Corners_Of_Elem(target_);
57
return UG_NS<coorddim>::Corners_Of_Elem(target_);
61
//! access to coordinates of corners. Index is the number of the corner
60
//! access to coordinates of corners. Index is the number of the corner
62
61
FieldVector<UGCtype, coorddim> corner (int i) const;
64
/** \brief Maps a local coordinate within reference element to
63
/** \brief Maps a local coordinate within reference element to
65
64
* global coordinate in element */
66
65
FieldVector<UGCtype, coorddim> global (const FieldVector<UGCtype, mydim>& local) const;
68
/** \brief Maps a global coordinate within the element to a
67
/** \brief Maps a global coordinate within the element to a
69
68
* local coordinate in its reference element */
70
69
FieldVector<UGCtype, mydim> local (const FieldVector<UGCtype, coorddim>& global) const;
73
Integration over a general element is done by integrating over the reference element
74
and using the transformation from the reference element to the global element as follows:
75
\f[\int\limits_{\Omega_e} f(x) dx = \int\limits_{\Omega_{ref}} f(g(l)) A(l) dl \f] where
76
\f$g\f$ is the local to global mapping and \f$A(l)\f$ is the integration element.
78
For a general map \f$g(l)\f$ involves partial derivatives of the map (surface element of
79
the first kind if \f$d=2,w=3\f$, determinant of the Jacobian of the transformation for
80
\f$d=w\f$, \f$\|dg/dl\|\f$ for \f$d=1\f$).
82
For linear elements, the derivatives of the map with respect to local coordinates
83
do not depend on the local coordinates and are the same over the whole element.
85
For a structured mesh where all edges are parallel to the coordinate axes, the
86
computation is the length, area or volume of the element is very simple to compute.
88
Each grid module implements the integration element with optimal efficiency. This
89
will directly translate in substantial savings in the computation of finite element
72
Integration over a general element is done by integrating over the reference element
73
and using the transformation from the reference element to the global element as follows:
74
\f[\int\limits_{\Omega_e} f(x) dx = \int\limits_{\Omega_{ref}} f(g(l)) A(l) dl \f] where
75
\f$g\f$ is the local to global mapping and \f$A(l)\f$ is the integration element.
77
For a general map \f$g(l)\f$ involves partial derivatives of the map (surface element of
78
the first kind if \f$d=2,w=3\f$, determinant of the Jacobian of the transformation for
79
\f$d=w\f$, \f$\|dg/dl\|\f$ for \f$d=1\f$).
81
For linear elements, the derivatives of the map with respect to local coordinates
82
do not depend on the local coordinates and are the same over the whole element.
84
For a structured mesh where all edges are parallel to the coordinate axes, the
85
computation is the length, area or volume of the element is very simple to compute.
87
Each grid module implements the integration element with optimal efficiency. This
88
will directly translate in substantial savings in the computation of finite element
92
91
UGCtype integrationElement (const FieldVector<UGCtype, mydim>& local) const;
94
93
UGCtype volume() const {
99
// coorddim*coorddim is an upper bound for the number of vertices
100
UGCtype* cornerCoords[coorddim*coorddim];
101
UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
102
return UG_NS<coorddim>::Area_Of_Element(corners(),
103
const_cast<const double**>(cornerCoords));
98
// coorddim*coorddim is an upper bound for the number of vertices
99
UGCtype* cornerCoords[coorddim*coorddim];
100
UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
101
return UG_NS<coorddim>::Area_Of_Element(corners(),
102
const_cast<const double**>(cornerCoords));
106
105
//! The inverse transpose of the Jacobian matrix of the mapping from the reference element to this element
109
108
const FieldMatrix<UGCtype, mydim,coorddim>& jacobianTransposed (const FieldVector<UGCtype, mydim>& local) const;
114
/** \brief Init the element with a given UG element */
115
void setToTarget(typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target)
118
jacobianIsUpToDate_ = false;
119
jacobianInverseIsUpToDate_ = false;
122
//! The jacobian inverse transposed
123
mutable FieldMatrix<UGCtype,coorddim,mydim> jac_inverse_;
124
//! The jacobian transposed
125
mutable FieldMatrix<UGCtype,mydim,coorddim> jac_;
113
/** \brief Init the element with a given UG element */
114
void setToTarget(typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target)
117
jacobianIsUpToDate_ = false;
118
jacobianInverseIsUpToDate_ = false;
121
//! The jacobian inverse transposed
122
mutable FieldMatrix<UGCtype,coorddim,mydim> jac_inverse_;
123
//! The jacobian transposed
124
mutable FieldMatrix<UGCtype,mydim,coorddim> jac_;
127
126
/** \brief For simplices the Jacobian matrix is a constant, hence it needs
128
127
to be computed only once for each new element. This can save some
130
129
mutable bool jacobianInverseIsUpToDate_;
131
130
mutable bool jacobianIsUpToDate_;
133
// in element mode this points to the element we map to
134
// in coord_mode this is the element whose reference element is mapped into the father's one
135
typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target_;
141
/****************************************************************/
143
/* Specialization for faces in 3d */
145
/****************************************************************/
147
template<class GridImp>
148
class UGGridGeometry<2, 3, GridImp> :
149
public GenericGeometry::BasicGeometry<2, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,2,3> >
152
template <int codim_, int dim_, class GridImp_>
153
friend class UGGridEntity;
155
template <class GridImp_>
156
friend class UGGridIntersectionIterator;
158
typedef typename GenericGeometry::BasicGeometry<2, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,2,3> > Base;
161
static const int mydimension = 2;
162
static const int coorddimension = 3;
164
/** \brief Setup method with a geometry type and a set of corners
165
\param coordinates The corner coordinates in DUNE numbering
167
void setup(const GeometryType& type, const std::vector<FieldVector<typename GridImp::ctype,3> >& coordinates)
170
// Yes, a strange way, but the only way, as BasicGeometry doesn't have a setup method
171
static_cast< Base & >( *this ) = Base( type, coordinates );
177
/****************************************************************/
179
/* Specialization for faces in 2d */
181
/****************************************************************/
183
template<class GridImp>
184
class UGGridGeometry <1, 2, GridImp> :
185
public GenericGeometry::BasicGeometry<1, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,1,2> >
188
template <int codim_, int dim_, class GridImp_>
189
friend class UGGridEntity;
191
template <class GridImp_>
192
friend class UGGridIntersectionIterator;
194
typedef typename GenericGeometry::BasicGeometry<1, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,1,2> > Base;
197
static const int mydimension = 1;
198
static const int coorddimension = 2;
200
/** \brief Constructor with a geometry type and a set of corners */
201
void setup(const GeometryType& type, const std::vector<FieldVector<typename GridImp::ctype,2> >& coordinates)
204
// Yes, a strange way, but the only way, as BasicGeometry doesn't have a setup method
205
static_cast< Base & >( *this ) = Base( type, coordinates );
132
// in element mode this points to the element we map to
133
// in coord_mode this is the element whose reference element is mapped into the father's one
134
typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target_;
140
/****************************************************************/
142
/* Specialization for faces in 3d */
144
/****************************************************************/
146
template<class GridImp>
147
class UGGridGeometry<2, 3, GridImp> :
148
public CachedMultiLinearGeometry<typename GridImp::ctype, 2, 3>
152
/** \brief Constructor from a given geometry type and a vector of corner coordinates */
153
UGGridGeometry(const GeometryType& type, const std::vector<FieldVector<typename GridImp::ctype,3> >& coordinates)
154
: CachedMultiLinearGeometry<typename GridImp::ctype, 2, 3>(type, coordinates)
160
/****************************************************************/
162
/* Specialization for faces in 2d */
164
/****************************************************************/
166
template<class GridImp>
167
class UGGridGeometry <1, 2, GridImp> :
168
public CachedMultiLinearGeometry<typename GridImp::ctype,1,2>
172
/** \brief Constructor from a given geometry type and a vector of corner coordinates */
173
UGGridGeometry(const GeometryType& type, const std::vector<FieldVector<typename GridImp::ctype,2> >& coordinates)
174
: CachedMultiLinearGeometry<typename GridImp::ctype, 1, 2>(type, coordinates)
210
179
namespace FacadeOptions
added
removed
dune/grid/uggrid/uggridgeometry.hh
