Branch data Line data Source code
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 : 538893 : double QuadraticHex::SH( const int i, const double params )
32 : : {
33 [ + + + - ]: 538893 : switch( i )
34 : : {
35 : : case -1:
36 : 179631 : return ( params * params - params ) / 2;
37 : : case 0:
38 : 179631 : return 1 - params * params;
39 : : case 1:
40 : 179631 : return ( params * params + params ) / 2;
41 : : default:
42 : 0 : return 0.;
43 : : }
44 : : }
45 : 215541 : double QuadraticHex::DSH( const int i, const double params )
46 : : {
47 [ + + + - ]: 215541 : switch( i )
48 : : {
49 : : case -1:
50 : 71847 : return params - 0.5;
51 : : case 0:
52 : 71847 : return -2 * params;
53 : : case 1:
54 : 71847 : return params + 0.5;
55 : : default:
56 : 0 : return 0.;
57 : : }
58 : : }
59 : :
60 : 3992 : ErrorCode QuadraticHex::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples,
61 : : double* /*work*/, double* result )
62 : : {
63 [ + - ][ + - ]: 3992 : assert( params && field && num_tuples > 0 );
[ - + ]
64 [ + - ]: 3992 : std::fill( result, result + num_tuples, 0.0 );
65 [ + + ]: 111776 : for( int i = 0; i < 27; i++ )
66 : : {
67 : 107784 : const double sh = SH( corner[i][0], params[0] ) * SH( corner[i][1], params[1] ) * SH( corner[i][2], params[2] );
68 [ + + ]: 431136 : for( int j = 0; j < num_tuples; j++ )
69 : 323352 : result[j] += sh * field[num_tuples * i + j];
70 : : }
71 : :
72 : 3992 : return MB_SUCCESS;
73 : : }
74 : :
75 : 2661 : ErrorCode QuadraticHex::jacobianFcn( const double* params, const double* verts, const int nverts, const int ndim,
76 : : double* /*work*/, double* result )
77 : : {
78 [ + - ][ + - ]: 2661 : assert( 27 == nverts && params && verts );
[ - + ]
79 [ - + ]: 2661 : if( 27 != nverts ) return MB_FAILURE;
80 : 2661 : Matrix3* J = reinterpret_cast< Matrix3* >( result );
81 [ + + ]: 74508 : for( int i = 0; i < 27; i++ )
82 : : {
83 [ + - ][ + - ]: 71847 : const double sh[3] = { SH( corner[i][0], params[0] ), SH( corner[i][1], params[1] ),
84 [ + - ]: 71847 : SH( corner[i][2], params[2] ) };
85 [ + - ][ + - ]: 71847 : const double dsh[3] = { DSH( corner[i][0], params[0] ), DSH( corner[i][1], params[1] ),
86 [ + - ]: 71847 : DSH( corner[i][2], params[2] ) };
87 : :
88 [ + + ]: 287388 : for( int j = 0; j < 3; j++ )
89 : : {
90 [ + - ]: 215541 : ( *J )( j, 0 ) += dsh[0] * sh[1] * sh[2] * verts[ndim * i + j]; // dxj/dr first column
91 [ + - ]: 215541 : ( *J )( j, 1 ) += sh[0] * dsh[1] * sh[2] * verts[ndim * i + j]; // dxj/ds
92 [ + - ]: 215541 : ( *J )( j, 2 ) += sh[0] * sh[1] * dsh[2] * verts[ndim * i + j]; // dxj/dt
93 : : }
94 : : }
95 : :
96 : 2661 : return MB_SUCCESS;
97 : : }
98 : :
99 : 0 : ErrorCode QuadraticHex::integrateFcn( const double* /*field*/, const double* /*verts*/, const int /*nverts*/,
100 : : const int /*ndim*/, const int /*num_tuples*/, double* /*work*/,
101 : : double* /*result*/ )
102 : : {
103 : 0 : return MB_NOT_IMPLEMENTED;
104 : : }
105 : :
106 : 1331 : ErrorCode QuadraticHex::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
107 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
108 : : const double inside_tol, double* work, double* params, int* is_inside )
109 : : {
110 [ + - ][ - + ]: 1331 : assert( posn && verts );
111 : : return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
112 : 1331 : is_inside );
113 : : }
114 : :
115 : 2662 : int QuadraticHex::insideFcn( const double* params, const int ndim, const double tol )
116 : : {
117 : 2662 : return EvalSet::inside_function( params, ndim, tol );
118 : : }
119 : :
120 : 0 : ErrorCode QuadraticHex::normalFcn( const int /*ientDim*/, const int /*facet*/, const int /*nverts*/,
121 : : const double* /*verts*/, double* /*normal[3]*/ )
122 : : {
123 : 0 : return MB_NOT_IMPLEMENTED;
124 : : }
125 [ + - ][ + - ]: 228 : } // namespace moab
|
