Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
QuadraticHex.cpp
Go to the documentation of this file.
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