1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GEOGRID_INTERSECTION_HH
4 #define DUNE_GEOGRID_INTERSECTION_HH
5 
6 #include <dune/grid/geometrygrid/declaration.hh>
7 #include <dune/grid/geometrygrid/cornerstorage.hh>
8 
9 namespace Dune
10 {
11 
12   namespace GeoGrid
13   {
14 
15     // Intersection
16     // ------------
17 
18     template< class Grid, class HostIntersection >
19     class Intersection
20     {
21       typedef typename HostIntersection::Geometry HostGeometry;
22       typedef typename HostIntersection::LocalGeometry HostLocalGeometry;
23 
24       typedef typename std::remove_const< Grid >::type::Traits Traits;
25 
26     public:
27       typedef typename Traits::ctype ctype;
28 
29       static const int dimension = Traits::dimension;
30       static const int dimensionworld = Traits::dimensionworld;
31 
32       typedef typename Traits::template Codim< 0 >::Entity Entity;
33       typedef typename Traits::template Codim< 1 >::Geometry Geometry;
34       typedef typename Traits::template Codim< 1 >::LocalGeometry LocalGeometry;
35 
36       typedef typename Traits::template Codim< 0 >::Geometry ElementGeometry;
37 
38     private:
39       typedef GeoGrid::IntersectionCoordVector< Grid > CoordVector;
40 
41       typedef typename Traits::template Codim< 0 >::EntityImpl EntityImpl;
42 
43       typedef typename Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
44       typedef typename Traits::template Codim< 0 >::GeometryImpl ElementGeometryImpl;
45 
46     public:
47 
Intersection()48       Intersection()
49       {}
50 
Intersection(const HostIntersection & hostIntersection,const ElementGeometryImpl & insideGeo)51       explicit Intersection ( const HostIntersection &hostIntersection, const ElementGeometryImpl &insideGeo )
52         : hostIntersection_( hostIntersection )
53         , insideGeo_ ( insideGeo )
54         , geo_( grid() )
55       {}
56 
Intersection(HostIntersection && hostIntersection,const ElementGeometryImpl & insideGeo)57       explicit Intersection ( HostIntersection&& hostIntersection, const ElementGeometryImpl &insideGeo )
58         : hostIntersection_( std::move( hostIntersection ) )
59         , insideGeo_ ( insideGeo )
60         , geo_( grid() )
61       {}
62 
equals(const Intersection & other) const63       bool equals ( const Intersection &other) const
64       {
65         return hostIntersection_ == other.hostIntersection_;
66       }
67 
operator bool() const68       explicit operator bool () const { return bool( hostIntersection_ ); }
69 
inside() const70       Entity inside () const
71       {
72         return EntityImpl( insideGeo_, hostIntersection().inside() );
73       }
74 
outside() const75       Entity outside () const
76       {
77         return EntityImpl( grid(), hostIntersection().outside() );
78       }
79 
boundary() const80       bool boundary () const { return hostIntersection().boundary(); }
81 
conforming() const82       bool conforming () const { return hostIntersection().conforming(); }
83 
neighbor() const84       bool neighbor () const { return hostIntersection().neighbor(); }
85 
boundarySegmentIndex() const86       size_t boundarySegmentIndex () const
87       {
88         return hostIntersection().boundarySegmentIndex();
89       }
90 
geometryInInside() const91       LocalGeometry geometryInInside () const
92       {
93         return hostIntersection().geometryInInside();
94       }
95 
geometryInOutside() const96       LocalGeometry geometryInOutside () const
97       {
98         return hostIntersection().geometryInOutside();
99       }
100 
geometry() const101       Geometry geometry () const
102       {
103         if( !geo_ )
104         {
105           CoordVector coords( insideGeo_, geometryInInside() );
106           geo_ = GeometryImpl( grid(), type(), coords );
107         }
108         return Geometry( geo_ );
109       }
110 
type() const111       GeometryType type () const { return hostIntersection().type(); }
112 
indexInInside() const113       int indexInInside () const
114       {
115         return hostIntersection().indexInInside();
116       }
117 
indexInOutside() const118       int indexInOutside () const
119       {
120         return hostIntersection().indexInOutside();
121       }
122 
123       FieldVector< ctype, dimensionworld >
integrationOuterNormal(const FieldVector<ctype,dimension-1> & local) const124       integrationOuterNormal ( const FieldVector< ctype, dimension-1 > &local ) const
125       {
126         const LocalGeometry geoInInside = geometryInInside();
127         const int idxInInside = indexInInside();
128 
129         auto refElement = referenceElement< ctype, dimension >( insideGeo_.type() );
130 
131         FieldVector< ctype, dimension > x( geoInInside.global( local ) );
132         const typename ElementGeometryImpl::JacobianInverseTransposed &jit = insideGeo_.jacobianInverseTransposed( x );
133         FieldVector< ctype, dimension > refNormal = refElement.integrationOuterNormal( idxInInside );
134 
135         FieldVector< ctype, dimensionworld > normal;
136         jit.mv( refNormal, normal );
137         if( !conforming() )
138           normal *= geoInInside.volume() / refElement.template geometry< 1 >( idxInInside ).volume();
139         normal *= jit.detInv();
140         //normal *= insideGeo_.integrationElement( x );
141         return normal;
142       }
143 
144       FieldVector< ctype, dimensionworld >
outerNormal(const FieldVector<ctype,dimension-1> & local) const145       outerNormal ( const FieldVector< ctype, dimension-1 > &local ) const
146       {
147         auto refElement = referenceElement< ctype, dimension >( insideGeo_.type() );
148 
149         FieldVector< ctype, dimension > x( geometryInInside().global( local ) );
150         const typename ElementGeometryImpl::JacobianInverseTransposed &jit = insideGeo_.jacobianInverseTransposed( x );
151         FieldVector< ctype, dimension > refNormal = refElement.integrationOuterNormal( indexInInside() );
152 
153         FieldVector< ctype, dimensionworld > normal;
154         jit.mv( refNormal, normal );
155         return normal;
156       }
157 
158       FieldVector< ctype, dimensionworld >
unitOuterNormal(const FieldVector<ctype,dimension-1> & local) const159       unitOuterNormal ( const FieldVector< ctype, dimension-1 > &local ) const
160       {
161         FieldVector< ctype, dimensionworld > normal = outerNormal( local );
162         normal *= (ctype( 1 ) / normal.two_norm());
163         return normal;
164       }
165 
centerUnitOuterNormal() const166       FieldVector< ctype, dimensionworld > centerUnitOuterNormal () const
167       {
168         auto refFace = referenceElement< ctype, dimension-1 >( type() );
169         return unitOuterNormal( refFace.position( 0, 0 ) );
170       }
171 
hostIntersection() const172       const HostIntersection &hostIntersection () const
173       {
174         return hostIntersection_;
175       }
176 
grid() const177       const Grid &grid () const { return insideGeo_.grid(); }
178 
179     private:
180       HostIntersection hostIntersection_;
181       ElementGeometryImpl insideGeo_;
182       mutable GeometryImpl geo_;
183     };
184 
185   } // namespace GeoGrid
186 
187 } // namespace Dune
188 
189 #endif // #ifndef DUNE_GEOGRID_INTERSECTION_HH
190