1 #ifndef __MAPP_TETRA_3D_EXT_HH__
2 #define __MAPP_TETRA_3D_EXT_HH__
3 
4 // BSGridALUVecType is defined in BSGrid interface to dune or in gitter_sti.hh
5 class BSGridLinearSurfaceMapping
6 {
7 protected:
8   double _n [3] ;
9   const double (&_p0)[3], (&_p1)[3], (&_p2)[3] ;
10 public:
11   inline BSGridLinearSurfaceMapping
12     (const double (&)[3], const double (&)[3], const double (&)[3]);
13 
14   // same as method normal of LinearSurfaceMapping, just for Dune Vecs
15   inline void normal(double * normal) const;
16 
17   // same as method normal of LinearSurfaceMapping, just for Dune Vecs
18   template <class aluvec_t>
19   inline void normal(aluvec_t & normal) const;
20 };
21 
22 class BSGridBilinearSurfaceMapping {
23 protected:
24   double _n[3][3];
25   double _b[4][3];
26   const double (&_p0)[3];
27   const double (&_p1)[3];
28   const double (&_p2)[3];
29   const double (&_p3)[3];
30 public:
31   inline BSGridBilinearSurfaceMapping(const double(&)[3], const double(&)[3],
32                                       const double(&)[3], const double(&)[3]);
33 
34   inline void normal(const double* local, double* normal) const;
35 
36   template <class aluvec_t>
37   inline void normal(const aluvec_t& local, aluvec_t& normal) const;
38 };
39 
40 inline BSGridLinearSurfaceMapping ::
BSGridLinearSurfaceMapping(const double (& x0)[3],const double (& x1)[3],const double (& x2)[3])41 BSGridLinearSurfaceMapping (const double (&x0)[3],
42                             const double (&x1)[3], const double (&x2)[3])
43   //: LinearSurfaceMapping (x0,x1,x2)
44   : _p0 (x0), _p1 (x1), _p2 (x2)
45 {
46   // copied from LinearSurfaceMapping, ist a bit faster
47   _n[0] = -0.5 * ((_p1[1]-_p0[1]) *(_p2[2]-_p1[2]) - (_p2[1]-_p1[1]) *(_p1[2]-_p0[2])) ;
48   _n[1] = -0.5 * ((_p1[2]-_p0[2]) *(_p2[0]-_p1[0]) - (_p2[2]-_p1[2]) *(_p1[0]-_p0[0])) ;
49   _n[2] = -0.5 * ((_p1[0]-_p0[0]) *(_p2[1]-_p1[1]) - (_p2[0]-_p1[0]) *(_p1[1]-_p0[1])) ;
50 }
51 
normal(double * normal)52 inline void BSGridLinearSurfaceMapping :: normal (double * normal) const
53 {
54   normal[0] = this->_n[0];
55   normal[1] = this->_n[1];
56   normal[2] = this->_n[2];
57   return ;
58 }
59 
60 template <class aluvec_t>
normal(aluvec_t & normal)61 inline void BSGridLinearSurfaceMapping :: normal (aluvec_t & normal) const
62 {
63   normal[0] = this->_n[0];
64   normal[1] = this->_n[1];
65   normal[2] = this->_n[2];
66   return ;
67 }
68 
69 inline BSGridBilinearSurfaceMapping::
BSGridBilinearSurfaceMapping(const double (& x0)[3],const double (& x1)[3],const double (& x2)[3],const double (& x3)[3])70 BSGridBilinearSurfaceMapping(const double(&x0)[3], const double(&x1)[3],
71                              const double(&x2)[3], const double(&x3)[3]) :
72   _p0 (x0), _p1 (x1), _p2 (x2), _p3 (x3) {
73   _b [0][0] = _p0 [0] ;
74   _b [0][1] = _p0 [1] ;
75   _b [0][2] = _p0 [2] ;
76   _b [1][0] = _p3 [0] - _p0 [0] ;
77   _b [1][1] = _p3 [1] - _p0 [1] ;
78   _b [1][2] = _p3 [2] - _p0 [2] ;
79   _b [2][0] = _p1 [0] - _p0 [0] ;
80   _b [2][1] = _p1 [1] - _p0 [1] ;
81   _b [2][2] = _p1 [2] - _p0 [2] ;
82   _b [3][0] = _p2 [0] - _p1 [0] - _b [1][0] ;
83   _b [3][1] = _p2 [1] - _p1 [1] - _b [1][1] ;
84   _b [3][2] = _p2 [2] - _p1 [2] - _b [1][2] ;
85   _n [0][0] = _b [1][1] * _b [2][2] - _b [1][2] * _b [2][1] ;
86   _n [0][1] = _b [1][2] * _b [2][0] - _b [1][0] * _b [2][2] ;
87   _n [0][2] = _b [1][0] * _b [2][1] - _b [1][1] * _b [2][0] ;
88   _n [1][0] = _b [1][1] * _b [3][2] - _b [1][2] * _b [3][1] ;
89   _n [1][1] = _b [1][2] * _b [3][0] - _b [1][0] * _b [3][2] ;
90   _n [1][2] = _b [1][0] * _b [3][1] - _b [1][1] * _b [3][0] ;
91   _n [2][0] = _b [3][1] * _b [2][2] - _b [3][2] * _b [2][1] ;
92   _n [2][1] = _b [3][2] * _b [2][0] - _b [3][0] * _b [2][2] ;
93   _n [2][2] = _b [3][0] * _b [2][1] - _b [3][1] * _b [2][0] ;
94   return ;
95 }
96 
97 inline void
normal(const double * local,double * normal)98 BSGridBilinearSurfaceMapping::normal(const double* local, double* normal) const {
99   double x = local [0];
100   double y = local [1];
101   normal [0] = -( _n [0][0] + _n [1][0] * x + _n [2][0] * y) ;
102   normal [1] = -( _n [0][1] + _n [1][1] * x + _n [2][1] * y) ;
103   normal [2] = -( _n [0][2] + _n [1][2] * x + _n [2][2] * y) ;
104   return ;
105 }
106 
107 template <class aluvec_t>
normal(const aluvec_t & local,aluvec_t & normal)108 inline void BSGridBilinearSurfaceMapping::normal(const aluvec_t& local, aluvec_t& normal) const {
109   double x = local [0];
110   double y = local [1];
111   normal [0] = -( _n [0][0] + _n [1][0] * x + _n [2][0] * y) ;
112   normal [1] = -( _n [0][1] + _n [1][1] * x + _n [2][1] * y) ;
113   normal [2] = -( _n [0][2] + _n [1][2] * x + _n [2][2] * y) ;
114   return ;
115 }
116 
117 
118 
119 #endif
120