 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "moab/LocalDiscretization/QuadraticHex.hpp"
2 #include "moab/Forward.hpp"
3
4 namespace moab
5 {
6
7 // those are not just the corners, but for simplicity, keep this name
8 //
9 const int QuadraticHex::corner[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, // corner nodes: 0-7
10 { -1, 1, -1 }, // mid-edge nodes: 8-19
11 { -1, -1, 1 }, // center-face nodes 20-25 center node 26
12 { 1, -1, 1 }, //
13 { 1, 1, 1 }, { -1, 1, 1 }, // 4 ----- 19 ----- 7
14 { 0, -1, -1 }, // . | . |
15 { 1, 0, -1 }, // 16 25 18 |
16 { 0, 1, -1 }, // . | . |
17 { -1, 0, -1 }, // 5 ----- 17 ----- 6 |
18 { -1, -1, 0 }, // | 12 | 23 15
19 { 1, -1, 0 }, // | | |
20 { 1, 1, 0 }, // | 20 | 26 | 22 |
21 { -1, 1, 0 }, // | | |
22 { 0, -1, 1 }, // 13 21 | 14 |
23 { 1, 0, 1 }, // | 0 ----- 11 ----- 3
24 { 0, 1, 1 }, // | . | .
25 { -1, 0, 1 }, // | 8 24 | 10
26 { 0, -1, 0 }, // | . | .
27 { 1, 0, 0 }, // 1 ----- 9 ----- 2
28 { 0, 1, 0 }, //
29 { -1, 0, 0 }, { 0, 0, -1 }, { 0, 0, 1 }, { 0, 0, 0 } };
30
31 double QuadraticHex::SH( const int i, const double params )
32 {
33 switch( i )
34 {
35 case -1:
36 return ( params * params - params ) / 2;
37 case 0:
38 return 1 - params * params;
39 case 1:
40 return ( params * params + params ) / 2;
41 default:
42 return 0.;
43 }
44 }
45 double QuadraticHex::DSH( const int i, const double params )
46 {
47 switch( i )
48 {
49 case -1:
50 return params - 0.5;
51 case 0:
52 return -2 * params;
53 case 1:
54 return params + 0.5;
55 default:
56 return 0.;
57 }
58 }
59
60 ErrorCode QuadraticHex::evalFcn( const double* params,
61 const double* field,
62 const int /*ndim*/,
63 const int num_tuples,
64 double* /*work*/,
65 double* result )
66 {
67 assert( params && field && num_tuples > 0 );
68 std::fill( result, result + num_tuples, 0.0 );
69 for( int i = 0; i < 27; i++ )
70 {
71 const double sh = SH( corner[i][0], params[0] ) * SH( corner[i][1], params[1] ) * SH( corner[i][2], params[2] );
72 for( int j = 0; j < num_tuples; j++ )
73 result[j] += sh * field[num_tuples * i + j];
74 }
75
76 return MB_SUCCESS;
77 }
78
79 ErrorCode QuadraticHex::jacobianFcn( const double* params,
80 const double* verts,
81 const int nverts,
82 const int ndim,
83 double* /*work*/,
84 double* result )
85 {
86 assert( 27 == nverts && params && verts );
87 if( 27 != nverts ) return MB_FAILURE;
88 Matrix3* J = reinterpret_cast< Matrix3* >( result );
89 for( int i = 0; i < 27; i++ )
90 {
91 const double sh[3] = { SH( corner[i][0], params[0] ), SH( corner[i][1], params[1] ),
92 SH( corner[i][2], params[2] ) };
93 const double dsh[3] = { DSH( corner[i][0], params[0] ), DSH( corner[i][1], params[1] ),
94 DSH( corner[i][2], params[2] ) };
95
96 for( int j = 0; j < 3; j++ )
97 {
98 ( *J )( j, 0 ) += dsh[0] * sh[1] * sh[2] * verts[ndim * i + j]; // dxj/dr first column
99 ( *J )( j, 1 ) += sh[0] * dsh[1] * sh[2] * verts[ndim * i + j]; // dxj/ds
100 ( *J )( j, 2 ) += sh[0] * sh[1] * dsh[2] * verts[ndim * i + j]; // dxj/dt
101 }
102 }
103
104 return MB_SUCCESS;
105 }
106
107 ErrorCode QuadraticHex::integrateFcn( const double* /*field*/,
108 const double* /*verts*/,
109 const int /*nverts*/,
110 const int /*ndim*/,
111 const int /*num_tuples*/,
112 double* /*work*/,
113 double* /*result*/ )
114 {
115 return MB_NOT_IMPLEMENTED;
116 }
117
118 ErrorCode QuadraticHex::reverseEvalFcn( EvalFcn eval,
119 JacobianFcn jacob,
120 InsideFcn ins,
121 const double* posn,
122 const double* verts,
123 const int nverts,
124 const int ndim,
125 const double iter_tol,
126 const double inside_tol,
127 double* work,
128 double* params,
129 int* is_inside )
130 {
131 assert( posn && verts );
132 return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
133 is_inside );
134 }
135
136 int QuadraticHex::insideFcn( const double* params, const int ndim, const double tol )
137 {
138 return EvalSet::inside_function( params, ndim, tol );
139 }
140
141 ErrorCode QuadraticHex::normalFcn( const int /*ientDim*/,
142 const int /*facet*/,
143 const int /*nverts*/,
144 const double* /*verts*/,
145 double* /*normal[3]*/ )
146 {
147 return MB_NOT_IMPLEMENTED;
148 }
149 } // namespace moab
1.9.1.