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