dune-fem  2.8-git
geometrygridpart/intersection.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_GRIDPART_GEOMETRYGRIDPART_INTERSECTION_HH
2 #define DUNE_FEM_GRIDPART_GEOMETRYGRIDPART_INTERSECTION_HH
3 
4 #include <type_traits>
5 
6 #include <dune/common/version.hh>
8 
9 #include <dune/geometry/affinegeometry.hh>
11 
12 namespace Dune
13 {
14 
15  namespace Fem
16  {
17 
18  // GeometryGridPartIntersection
19  // ----------------------------
20 
21  template< class GridFamily >
23  {
25  typedef typename std::remove_const< GridFamily >::type::Traits Traits;
26 
27  public:
28  typedef typename std::remove_const< GridFamily >::type::ctype ctype;
29 
30  static const int dimension = std::remove_const< GridFamily >::type::dimension;
31  static const int dimensionworld = std::remove_const< GridFamily >::type::dimensionworld;
32 
33  typedef typename Traits::template Codim< 0 >::Entity Entity;
34  typedef typename Traits::template Codim< 0 >::Geometry ElementGeometry;
35  typedef typename Traits::template Codim< 1 >::Geometry Geometry;
36  typedef typename Traits::template Codim< 1 >::LocalGeometry LocalGeometry;
37 
38  typedef typename Traits::GridFunctionType GridFunctionType;
39 
40  typedef FieldVector< ctype, dimensionworld > GlobalCoordinate;
41  typedef FieldVector< ctype, dimension-1 > LocalCoordinate;
42 
43  private:
44  typedef typename Geometry::Implementation GeometryImplType;
45 
46  typedef typename Traits::HostGridPartType HostGridPartType;
47  typedef typename HostGridPartType::IntersectionType HostIntersectionType;
48 
49  public:
51  GeometryGridPartIntersection ( const GridFunctionType &gridFunction, const typename ElementGeometry::Implementation &insideGeo, HostIntersectionType hostIntersection )
52  : hostIntersection_( std::move( hostIntersection ) ), gridFunction_( &gridFunction ), insideGeo_( insideGeo )
53  {}
54 
55  operator bool () const { return bool( hostIntersection_ ); }
56 
57  Entity inside () const { return typename Entity::Implementation( gridFunction(), hostIntersection().inside() ); }
58  Entity outside () const { return typename Entity::Implementation( gridFunction(), hostIntersection().outside() ); }
59 
60  int boundaryId () const {
62  }
63 
64  bool boundary () const { return hostIntersection().boundary(); }
65 
66  bool conforming () const { return hostIntersection().conforming(); }
67 
68  int twistInSelf() const { return hostIntersection().impl().twistInSelf(); }
69 
70  int twistInNeighbor() const { return hostIntersection().impl().twistInNeighbor(); }
71 
72  bool neighbor () const { return hostIntersection().neighbor(); }
73 
74  std::size_t boundarySegmentIndex () const { return hostIntersection().boundarySegmentIndex(); }
75 
76  LocalGeometry geometryInInside () const { return hostIntersection().geometryInInside(); }
77  LocalGeometry geometryInOutside () const { return hostIntersection().geometryInOutside(); }
78 
79  Geometry geometry () const
80  {
81  typedef typename Geometry::Implementation Impl;
82  return Geometry( Impl( insideGeo_, geometryInInside(), 2*insideGeo_.impl().localFunction().order()+1 ) );
83  }
84 
85  bool equals ( const This &other ) const { return hostIntersection() == other.hostIntersection(); }
86 
87  GeometryType type () const { return hostIntersection().type(); }
88 
89  int indexInInside () const { return hostIntersection().indexInInside(); }
90  int indexInOutside () const { return hostIntersection().indexInOutside(); }
91 
93  {
94  const auto &refElement = ReferenceElements< ctype, dimension >::general( insideGeo_.type() );
95  const auto &refNormal = refElement.integrationOuterNormal( indexInInside() );
96 
97  const auto jit = insideGeo_.jacobianInverseTransposed( geometryInInside().global( local ) );
98 
99  GlobalCoordinate normal;
100  jit.mv( refNormal, normal );
101  // double det = std::sqrt( GeometryGridPartGeometryType::MatrixHelper::template detATA<dimensionworld,dimension>( jit ) );
102  // return normal *= ctype( 1 ) / sqrt(det);
103  return normal *= geometry().integrationElement( local ) / normal.two_norm();
104 
105 #if 0
107 FieldVector< ctype, dimension > x( geometryInInside().global( local ) );
108 FieldVector< ctype, dimensionworld > normal, tau, nu;
109 FieldMatrix< ctype, 1, dimensionworld > tauT = geometry().jacobianTransposed(local);
110 
111 for (int i = 0; i != dimensionworld; ++i)
112  tau[i] = tauT[0][i];
113 
114 tau /= tau.two_norm();
115 
116 FieldMatrix< ctype, dimensionworld, dimensionworld > localFnGrad, localFnGradT;
117 gridFunction().localFunction(*hostIntersection().inside()).jacobian(x, localFnGrad);
118 
119 for (int i = 0; i != dimensionworld; ++i)
120  for (int j = 0; j != dimensionworld; ++j )
121  localFnGradT[i][j] = localFnGrad[j][i];
122 
123 crossProduct(localFnGradT[0], localFnGradT[1], nu );
124 nu /= nu.two_norm();
125 
126 /*
127 nu = insideGeo_.global(x);
128 nu /= nu.two_norm();
129 */
130 
131 crossProduct(nu, tau, normal);
132 
133 // get conormal in same direction as discrete conormal
134 if (normal*hostIntersection().integrationOuterNormal(local) < 0)
135  normal *= -1;
136 
137 // normal /= normal.two_norm();
138 // normal *= GeometryImplType(insideGeo_, gridFunction_, hostIntersection_, affineGeometry).integrationElement(local);
139 normal *= geometry().integrationElement(local)/normal.two_norm();
140  // return hostIntersection().integrationOuterNormal(local);
141  return normal;
142 #endif
143  }
144 
145 #if 0
146  void crossProduct(const FieldVector< ctype, dimensionworld > &vec1, const FieldVector< ctype, dimensionworld > &vec2, FieldVector< ctype, dimensionworld > &ret) const
147  {
148  assert( dimensionworld == 3 );
149 
150  ret[0] = vec1[1]*vec2[2] - vec1[2]*vec2[1];
151  ret[1] = vec1[2]*vec2[0] - vec1[0]*vec2[2];
152  ret[2] = vec1[0]*vec2[1] - vec1[1]*vec2[0];
153  }
154 #endif
155 
157  {
158  const auto &refElement = Dune::ReferenceElements< ctype, dimension >::general( insideGeo_.type() );
159  const auto &refNormal = refElement.integrationOuterNormal( indexInInside() );
160 
161  GlobalCoordinate normal;
162  insideGeo_.jacobianInverseTransposed( geometryInInside().global( local ) ).mv( refNormal, normal );
163  return normal;
164  }
165 
167  {
168  GlobalCoordinate normal = outerNormal( local );
169  return normal *= (ctype( 1 ) / normal.two_norm());
170  }
171 
173  {
174  const auto &refElement = Dune::ReferenceElements< ctype, dimension >::general( insideGeo_.type() );
175  const auto &refNormal = refElement.integrationOuterNormal( indexInInside() );
176 
177  GlobalCoordinate normal;
178  insideGeo_.jacobianInverseTransposed( geometryInInside().center() ).mv( refNormal, normal );
179  return normal *= (ctype( 1 ) / normal.two_norm());
180  }
181 
182  const HostIntersectionType &hostIntersection () const { return hostIntersection_; }
183 
185  {
186  assert( gridFunction_ );
187  return *gridFunction_;
188  }
189 
190  private:
191  HostIntersectionType hostIntersection_;
192  const GridFunctionType *gridFunction_ = nullptr;
193  typename ElementGeometry::Implementation insideGeo_;
194  };
195 
196  } // namespace Fem
197 
198 } // namespace Dune
199 
200 #endif // #ifndef DUNE_FEM_GRIDPART_GEOMETRYGRIDPART_INTERSECTION_HH
Definition: bindguard.hh:11
Definition: geometrygridpart/intersection.hh:23
Geometry geometry() const
Definition: geometrygridpart/intersection.hh:79
GeometryGridPartIntersection(const GridFunctionType &gridFunction, const typename ElementGeometry::Implementation &insideGeo, HostIntersectionType hostIntersection)
Definition: geometrygridpart/intersection.hh:51
const GridFunctionType & gridFunction() const
Definition: geometrygridpart/intersection.hh:184
int indexInOutside() const
Definition: geometrygridpart/intersection.hh:90
bool boundary() const
Definition: geometrygridpart/intersection.hh:64
FieldVector< ctype, dimensionworld > GlobalCoordinate
Definition: geometrygridpart/intersection.hh:40
FieldVector< ctype, dimension-1 > LocalCoordinate
Definition: geometrygridpart/intersection.hh:41
int twistInNeighbor() const
Definition: geometrygridpart/intersection.hh:70
static const int dimension
Definition: geometrygridpart/intersection.hh:30
LocalGeometry geometryInOutside() const
Definition: geometrygridpart/intersection.hh:77
int indexInInside() const
Definition: geometrygridpart/intersection.hh:89
std::remove_const< GridFamily >::type::ctype ctype
Definition: geometrygridpart/intersection.hh:28
GlobalCoordinate integrationOuterNormal(const LocalCoordinate &local) const
Definition: geometrygridpart/intersection.hh:92
Traits::template Codim< 0 >::Geometry ElementGeometry
Definition: geometrygridpart/intersection.hh:34
Traits::template Codim< 1 >::Geometry Geometry
Definition: geometrygridpart/intersection.hh:35
static const int dimensionworld
Definition: geometrygridpart/intersection.hh:31
GlobalCoordinate centerUnitOuterNormal() const
Definition: geometrygridpart/intersection.hh:172
bool equals(const This &other) const
Definition: geometrygridpart/intersection.hh:85
LocalGeometry geometryInInside() const
Definition: geometrygridpart/intersection.hh:76
Traits::template Codim< 0 >::Entity Entity
Definition: geometrygridpart/intersection.hh:33
const HostIntersectionType & hostIntersection() const
Definition: geometrygridpart/intersection.hh:182
Traits::template Codim< 1 >::LocalGeometry LocalGeometry
Definition: geometrygridpart/intersection.hh:36
bool conforming() const
Definition: geometrygridpart/intersection.hh:66
std::size_t boundarySegmentIndex() const
Definition: geometrygridpart/intersection.hh:74
Entity inside() const
Definition: geometrygridpart/intersection.hh:57
GlobalCoordinate unitOuterNormal(const LocalCoordinate &local) const
Definition: geometrygridpart/intersection.hh:166
GeometryType type() const
Definition: geometrygridpart/intersection.hh:87
GlobalCoordinate outerNormal(const LocalCoordinate &local) const
Definition: geometrygridpart/intersection.hh:156
Entity outside() const
Definition: geometrygridpart/intersection.hh:58
Traits::GridFunctionType GridFunctionType
Definition: geometrygridpart/intersection.hh:38
bool neighbor() const
Definition: geometrygridpart/intersection.hh:72
int twistInSelf() const
Definition: geometrygridpart/intersection.hh:68
int boundaryId() const
Definition: geometrygridpart/intersection.hh:60
Definition: boundaryidprovider.hh:36