1 #ifndef DUNE_FEM_GRIDPART_GEOGRIDPART_INTERSECTION_HH
2 #define DUNE_FEM_GRIDPART_GEOGRIDPART_INTERSECTION_HH
3 
4 #include <type_traits>
5 #include <utility>
6 
7 #include <dune/fem/gridpart/geogridpart/cornerstorage.hh>
8 
9 namespace Dune
10 {
11 
12   namespace Fem
13   {
14 
15     // GeoIntersection
16     // --------------
17 
18     template< class GridFamily >
19     class GeoIntersection
20     {
21       typedef typename std::remove_const< GridFamily >::type::Traits Traits;
22 
23     public:
24       typedef typename std::remove_const< GridFamily >::type::ctype ctype;
25 
26       static const int dimension = std::remove_const< GridFamily >::type::dimension;
27       static const int dimensionworld = std::remove_const< GridFamily >::type::dimensionworld;
28 
29       typedef typename Traits::template Codim< 0 >::Entity Entity;
30       typedef typename Traits::template Codim< 0 >::Geometry ElementGeometry;
31       typedef typename Traits::template Codim< 1 >::Geometry Geometry;
32       typedef typename Traits::template Codim< 1 >::LocalGeometry LocalGeometry;
33 
34       typedef typename Traits::CoordFunctionType CoordFunctionType;
35 
36     private:
37       typedef typename Entity::Implementation EntityImplType;
38       typedef typename ElementGeometry::Implementation ElementGeometryImplType;
39       typedef typename Geometry::Implementation GeometryImplType;
40 
41       typedef typename Traits::HostGridPartType HostGridPartType;
42       typedef typename HostGridPartType::IntersectionType HostIntersectionType;
43 
44       typedef GeoIntersectionCoordVector< GridFamily > CoordVectorType;
45 
46     public:
GeoIntersection(const CoordFunctionType & coordFunction,const ElementGeometry & insideGeo,HostIntersectionType hostIntersection)47       GeoIntersection ( const CoordFunctionType &coordFunction, const ElementGeometry &insideGeo, HostIntersectionType hostIntersection )
48       : coordFunction_( &coordFunction ),
49         insideGeo_( insideGeo.impl() ),
50         hostIntersection_( std::move( hostIntersection ) )
51       {}
52 
inside() const53       Entity inside () const
54       {
55         return EntityImplType( coordFunction(), hostIntersection().inside() );
56       }
57 
outside() const58       Entity outside () const
59       {
60         return EntityImplType( coordFunction(), hostIntersection().outside() );
61       }
62 
boundary() const63       bool boundary () const
64       {
65         return hostIntersection().boundary();
66       }
67 
conforming() const68       bool conforming () const
69       {
70         return hostIntersection().conforming();
71       }
72 
neighbor() const73       bool neighbor () const
74       {
75         return hostIntersection().neighbor();
76       }
77 
boundaryId() const78       int boundaryId () const
79       {
80         return hostIntersection().boundaryId();
81       }
82 
boundarySegmentIndex() const83       size_t boundarySegmentIndex () const
84       {
85         return hostIntersection().boundarySegmentIndex();
86       }
87 
geometryInInside() const88       LocalGeometry geometryInInside () const
89       {
90         return hostIntersection().geometryInInside();
91       }
92 
geometryInOutside() const93       LocalGeometry geometryInOutside () const
94       {
95         return hostIntersection().geometryInOutside();
96       }
97 
geometry() const98       Geometry geometry () const
99       {
100         const LocalGeometry &localGeo = geometryInInside();
101         CoordVectorType coords( insideGeo_, localGeo );
102         return Geometry( GeometryImplType( type(), coords ) );
103       }
104 
type() const105       GeometryType type () const
106       {
107         return hostIntersection().type();
108       }
109 
indexInInside() const110       int indexInInside () const
111       {
112         return hostIntersection().indexInInside();
113       }
114 
indexInOutside() const115       int indexInOutside () const
116       {
117         return hostIntersection().indexInOutside();
118       }
119 
120       FieldVector< ctype, dimensionworld >
integrationOuterNormal(const FieldVector<ctype,dimension-1> & local) const121       integrationOuterNormal ( const FieldVector< ctype, dimension-1 > &local ) const
122       {
123         auto refElement = referenceElement< ctype, dimension>( insideGeo_.type() );
124 
125         FieldVector< ctype, dimension > x( geometryInInside().global( local ) );
126         typedef typename ElementGeometryImplType::JacobianInverseTransposed JacobianInverseTransposed;
127         const JacobianInverseTransposed &jit = insideGeo_.jacobianInverseTransposed( x );
128         const FieldVector< ctype, dimension > &refNormal = refElement.integrationOuterNormal( indexInInside() );
129 
130         FieldVector< ctype, dimensionworld > normal;
131         jit.mv( refNormal, normal );
132         normal *= ctype( 1 ) / jit.det();
133         return normal;
134       }
135 
136       FieldVector< ctype, dimensionworld >
outerNormal(const FieldVector<ctype,dimension-1> & local) const137       outerNormal ( const FieldVector< ctype, dimension-1 > &local ) const
138       {
139         auto refElement = referenceElement< ctype, dimension>( insideGeo_.type() );
140 
141         FieldVector< ctype, dimension > x( geometryInInside().global( local ) );
142         typedef typename ElementGeometryImplType::JacobianInverseTransposed JacobianInverseTransposed;
143         const JacobianInverseTransposed &jit = insideGeo_.jacobianInverseTransposed( x );
144         const FieldVector< ctype, dimension > &refNormal = refElement.integrationOuterNormal( indexInInside() );
145 
146         FieldVector< ctype, dimensionworld > normal;
147         jit.mv( refNormal, normal );
148         return normal;
149       }
150 
151       FieldVector< ctype, dimensionworld >
unitOuterNormal(const FieldVector<ctype,dimension-1> & local) const152       unitOuterNormal ( const FieldVector< ctype, dimension-1 > &local ) const
153       {
154         FieldVector< ctype, dimensionworld > normal = outerNormal( local );
155         normal *= (ctype( 1 ) / normal.two_norm());
156         return normal;
157       }
158 
centerUnitOuterNormal() const159       FieldVector< ctype, dimensionworld > centerUnitOuterNormal () const
160       {
161         auto refFace = referenceElement< ctype, dimension-1 >( type() );
162         return unitOuterNormal( refFace.position( 0, 0 ) );
163       }
164 
coordFunction() const165       const CoordFunctionType &coordFunction () const
166       {
167         assert( coordFunction_ );
168         return *coordFunction_;
169       }
170 
hostIntersection() const171       const HostIntersectionType &hostIntersection () const
172       {
173         return hostIntersection_;
174       }
175 
176     private:
177       const CoordFunctionType *coordFunction_ = nullptr;
178       ElementGeometryImplType insideGeo_;
179       HostIntersectionType hostIntersection_;
180     };
181 
182   } // namespace Fem
183 
184 } // namespace Dune
185 
186 #endif // #ifndef DUNE_FEM_GRIDPART_GEOGRIDPART_INTERSECTION_HH
187