/usr/include/dune/grid/uggrid/uggridgeometry.hh is in libdune-grid-dev 2.2.1-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 | #ifndef DUNE_UGGRIDGEOMETRY_HH
#define DUNE_UGGRIDGEOMETRY_HH
/** \file
* \brief The UGGridGeometry class and its specializations
*/
#include <dune/common/array.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/genericgeometry/geometry.hh>
#include "uggridrenumberer.hh"
namespace Dune {
/** \brief Defines the geometry part of a mesh entity.
* \ingroup UGGrid
\tparam mydim Dimension of the corresponding reference element
\tparam coorddim Each corner is a point with coorddim coordinates.
This version is actually used only for mydim==coorddim. The mydim<coorddim cases
are in specializations below.
*/
template<int mydim, int coorddim, class GridImp>
class UGGridGeometry :
public GeometryDefaultImplementation <mydim, coorddim, GridImp, UGGridGeometry>
{
typedef typename GridImp::ctype UGCtype;
template <int codim_, int dim_, class GridImp_>
friend class UGGridEntity;
public:
/** \brief Default constructor
*/
UGGridGeometry()
{
jacobianIsUpToDate_ = false;
jacobianInverseIsUpToDate_ = false;
}
/** \brief Return the element type identifier
*
* UGGrid supports triangles and quadrilaterals in 2D, and
* tetrahedra, pyramids, prisms, and hexahedra in 3D.
*/
GeometryType type () const;
//! returns true if type is simplex, false otherwise
bool affine() const { return type().isSimplex(); }
//! return the number of corners of this element.
int corners () const {
return UG_NS<coorddim>::Corners_Of_Elem(target_);
}
//! access to coordinates of corners. Index is the number of the corner
FieldVector<UGCtype, coorddim> corner (int i) const;
/** \brief Maps a local coordinate within reference element to
* global coordinate in element */
FieldVector<UGCtype, coorddim> global (const FieldVector<UGCtype, mydim>& local) const;
/** \brief Maps a global coordinate within the element to a
* local coordinate in its reference element */
FieldVector<UGCtype, mydim> local (const FieldVector<UGCtype, coorddim>& global) const;
/**
Integration over a general element is done by integrating over the reference element
and using the transformation from the reference element to the global element as follows:
\f[\int\limits_{\Omega_e} f(x) dx = \int\limits_{\Omega_{ref}} f(g(l)) A(l) dl \f] where
\f$g\f$ is the local to global mapping and \f$A(l)\f$ is the integration element.
For a general map \f$g(l)\f$ involves partial derivatives of the map (surface element of
the first kind if \f$d=2,w=3\f$, determinant of the Jacobian of the transformation for
\f$d=w\f$, \f$\|dg/dl\|\f$ for \f$d=1\f$).
For linear elements, the derivatives of the map with respect to local coordinates
do not depend on the local coordinates and are the same over the whole element.
For a structured mesh where all edges are parallel to the coordinate axes, the
computation is the length, area or volume of the element is very simple to compute.
Each grid module implements the integration element with optimal efficiency. This
will directly translate in substantial savings in the computation of finite element
stiffness matrices.
*/
UGCtype integrationElement (const FieldVector<UGCtype, mydim>& local) const;
UGCtype volume() const {
if (mydim==0)
return 1;
// coorddim*coorddim is an upper bound for the number of vertices
UGCtype* cornerCoords[coorddim*coorddim];
UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
return UG_NS<coorddim>::Area_Of_Element(corners(),
const_cast<const double**>(cornerCoords));
}
//! The inverse transpose of the Jacobian matrix of the mapping from the reference element to this element
const FieldMatrix<UGCtype, coorddim,mydim>& jacobianInverseTransposed (const FieldVector<UGCtype, mydim>& local) const;
//! The transpose of the Jacobian matrix of the mapping from the reference element to this element
const FieldMatrix<UGCtype, mydim,coorddim>& jacobianTransposed (const FieldVector<UGCtype, mydim>& local) const;
private:
/** \brief Init the element with a given UG element */
void setToTarget(typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target)
{
target_ = target;
jacobianIsUpToDate_ = false;
jacobianInverseIsUpToDate_ = false;
}
//! The jacobian inverse transposed
mutable FieldMatrix<UGCtype,coorddim,mydim> jac_inverse_;
//! The jacobian transposed
mutable FieldMatrix<UGCtype,mydim,coorddim> jac_;
/** \brief For simplices the Jacobian matrix is a constant, hence it needs
to be computed only once for each new element. This can save some
assembly time. */
mutable bool jacobianInverseIsUpToDate_;
mutable bool jacobianIsUpToDate_;
// in element mode this points to the element we map to
// in coord_mode this is the element whose reference element is mapped into the father's one
typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target_;
};
/****************************************************************/
/* */
/* Specialization for faces in 3d */
/* */
/****************************************************************/
template<class GridImp>
class UGGridGeometry<2, 3, GridImp> :
public GenericGeometry::BasicGeometry<2, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,2,3> >
{
template <int codim_, int dim_, class GridImp_>
friend class UGGridEntity;
template <class GridImp_>
friend class UGGridIntersectionIterator;
typedef typename GenericGeometry::BasicGeometry<2, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,2,3> > Base;
public:
static const int mydimension = 2;
static const int coorddimension = 3;
/** \brief Setup method with a geometry type and a set of corners
\param coordinates The corner coordinates in DUNE numbering
*/
void setup(const GeometryType& type, const std::vector<FieldVector<typename GridImp::ctype,3> >& coordinates)
{
// set up base class
// Yes, a strange way, but the only way, as BasicGeometry doesn't have a setup method
static_cast< Base & >( *this ) = Base( type, coordinates );
}
};
/****************************************************************/
/* */
/* Specialization for faces in 2d */
/* */
/****************************************************************/
template<class GridImp>
class UGGridGeometry <1, 2, GridImp> :
public GenericGeometry::BasicGeometry<1, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,1,2> >
{
template <int codim_, int dim_, class GridImp_>
friend class UGGridEntity;
template <class GridImp_>
friend class UGGridIntersectionIterator;
typedef typename GenericGeometry::BasicGeometry<1, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,1,2> > Base;
public:
static const int mydimension = 1;
static const int coorddimension = 2;
/** \brief Constructor with a geometry type and a set of corners */
void setup(const GeometryType& type, const std::vector<FieldVector<typename GridImp::ctype,2> >& coordinates)
{
// set up base class
// Yes, a strange way, but the only way, as BasicGeometry doesn't have a setup method
static_cast< Base & >( *this ) = Base( type, coordinates );
}
};
namespace FacadeOptions
{
/** \brief Switch on the new implementation for the Geometry interface class
* \deprecated Eventually the new implementation will be hardwired,
* and this switch may disappear without prior warning!
*/
template< int mydim, int cdim, class GridImp>
struct StoreGeometryReference<mydim,cdim,GridImp,UGGridGeometry>
{
static const bool v = false;
};
}
} // namespace Dune
#endif
|