1 // (c) bernhard schupp 1997 - 1998
2 #include <config.h>
3 
4 #include <dune/alugrid/common/alugrid_assert.hh>
5 #include <cmath>
6 #include <cstdlib>
7 
8 #include "mapp_cube_3d.h"
9 
10 namespace ALUGrid
11 {
12 
13   const alucoord_t TrilinearMapping :: _epsilon = 1.0e-8 ;
14   const alucoord_t QuadraturCube3Dbasis :: _p2 [4][3] = { { .816496580927726, .0,  .5773502691896258},
15                                                       { .0, .816496580927726, -.5773502691896258},
16                                                       { -.816496580927726, .0, .5773502691896258},
17                                                       { .0, -.816496580927726, -.5773502691896258} } ;
18   const alucoord_t QuadraturCube3Dbasis :: _p3 [8][3] = { {  .5773502691896258,  .5773502691896258,  .5773502691896258},
19                                                       { -.5773502691896258,  .5773502691896258,  .5773502691896258},
20                                                       {  .5773502691896258, -.5773502691896258,  .5773502691896258},
21                                                       { -.5773502691896258, -.5773502691896258,  .5773502691896258},
22                                                       {  .5773502691896258,  .5773502691896258, -.5773502691896258},
23                                                       { -.5773502691896258,  .5773502691896258, -.5773502691896258},
24                                                       {  .5773502691896258, -.5773502691896258, -.5773502691896258},
25                                                       { -.5773502691896258, -.5773502691896258, -.5773502691896258} } ;
26 
27   const alucoord_t QuadraturCube2Dbasis :: _p1 [2] = { .0, .0 } ;
28   const alucoord_t QuadraturCube2Dbasis :: _p3 [4][2] = { { .5773502691896258,  .5773502691896258 },
29                                                       {-.5773502691896258,  .5773502691896258 },
30                                                       { .5773502691896258, -.5773502691896258 },
31                                                       {-.5773502691896258, -.5773502691896258 } } ;
32 
33 
linear(const alucoord_t (& p)[3])34   void TrilinearMapping :: linear(const alucoord_t (&p)[3]) {
35     alucoord_t x = .5 * (p[0] + 1.) ;
36     alucoord_t y = .5 * (p[1] + 1.) ;
37     alucoord_t z = .5 * (p[2] + 1.) ;
38     alucoord_t t0 = .5 ;
39     alucoord_t t3 = y * z ;
40     alucoord_t t8 = x * z ;
41     alucoord_t t13 = x * y ;
42     Df[2][0] = t0 * ( a[1][2] + y * a[4][2] + z * a[6][2] + t3 * a[7][2] ) ;
43     Df[2][1] = t0 * ( a[2][2] + x * a[4][2] + z * a[5][2] + t8 * a[7][2] ) ;
44     Df[1][2] = t0 * ( a[3][1] + y * a[5][1] + x * a[6][1] + t13 * a[7][1] ) ;
45     Df[2][2] = t0 * ( a[3][2] + y * a[5][2] + x * a[6][2] + t13 * a[7][2] ) ;
46     Df[0][0] = t0 * ( a[1][0] + y * a[4][0] + z * a[6][0] + t3 * a[7][0] ) ;
47     Df[0][2] = t0 * ( a[3][0] + y * a[5][0] + x * a[6][0] + t13 * a[7][0] ) ;
48     Df[1][0] = t0 * ( a[1][1] + y * a[4][1] + z * a[6][1] + t3 * a[7][1] ) ;
49     Df[0][1] = t0 * ( a[2][0] + x * a[4][0] + z * a[5][0] + t8 * a[7][0] ) ;
50     Df[1][1] = t0 * ( a[2][1] + x * a[4][1] + z * a[5][1] + t8 * a[7][1] ) ;
51 
52   }
53 
det(const alucoord_t (& point)[3])54   alucoord_t TrilinearMapping :: det(const alucoord_t (&point)[3]) {
55           //  Determinante der Abbildung f:[-1,1]^3 -> Hexaeder im Punkt point.
56     linear (point) ;
57     return (DetDf = Df[0][0] * Df[1][1] * Df[2][2] - Df[0][0] * Df[1][2] * Df[2][1] -
58                     Df[1][0] * Df[0][1] * Df[2][2] + Df[1][0] * Df[0][2] * Df[2][1] +
59                     Df[2][0] * Df[0][1] * Df[1][2] - Df[2][0] * Df[0][2] * Df[1][1]) ;
60   }
61 
inverse(const alucoord_t (& p)[3])62   void TrilinearMapping :: inverse(const alucoord_t (&p)[3]) {
63           //  Kramer - Regel, det() rechnet Df und DetDf neu aus.
64     alucoord_t val = 1.0 / det(p) ;
65     Dfi[0][0] = ( Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1] ) * val ;
66     Dfi[0][1] = ( Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2] ) * val ;
67     Dfi[0][2] = ( Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1] ) * val ;
68     Dfi[1][0] = ( Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2] ) * val ;
69     Dfi[1][1] = ( Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0] ) * val ;
70     Dfi[1][2] = ( Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2] ) * val ;
71     Dfi[2][0] = ( Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0] ) * val ;
72     Dfi[2][1] = ( Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1] ) * val ;
73     Dfi[2][2] = ( Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0] ) * val ;
74     return ;
75   }
76 
world2map(const alucoord_t (& wld)[3],alucoord_t (& map)[3])77   void TrilinearMapping :: world2map (const alucoord_t (&wld)[3], alucoord_t (&map)[3]) {
78           //  Newton - Iteration zum Invertieren der Abbildung f.
79     alucoord_t err = 10.0 * _epsilon ;
80 #ifdef ALUGRIDDEBUG
81     int count = 0 ;
82 #endif
83     map [0] = map [1] = map [2] = .0 ;
84     do {
85       alucoord_t upd [3] ;
86       map2world (map, upd) ;
87       inverse (map) ;
88       alucoord_t u0 = upd [0] - wld [0] ;
89       alucoord_t u1 = upd [1] - wld [1] ;
90       alucoord_t u2 = upd [2] - wld [2] ;
91       alucoord_t c0 = Dfi [0][0] * u0 + Dfi [0][1] * u1 + Dfi [0][2] * u2 ;
92       alucoord_t c1 = Dfi [1][0] * u0 + Dfi [1][1] * u1 + Dfi [1][2] * u2 ;
93       alucoord_t c2 = Dfi [2][0] * u0 + Dfi [2][1] * u1 + Dfi [2][2] * u2 ;
94       map [0] -= c0 ;
95       map [1] -= c1 ;
96       map [2] -= c2 ;
97       err = fabs (c0) + fabs (c1) + fabs (c2) ;
98       alugrid_assert (count ++ < 1000) ;
99     } while (err > _epsilon) ;
100   }
101 
102 } // namespace ALUGrid
103