Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
LinearTet.cpp
Go to the documentation of this file.
1 #include "moab/LocalDiscretization/LinearTet.hpp"
2 #include "moab/Forward.hpp"
3 #include <algorithm>
4 #include <cmath>
5 #include <limits>
6 
7 namespace moab
8 {
9 
10 const double LinearTet::corner[4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
11 
12 ErrorCode LinearTet::initFcn( const double* verts, const int nverts, double*& work )
13 {
14  // allocate work array as:
15  // work[0..8] = T
16  // work[9..17] = Tinv
17  // work[18] = detT
18  // work[19] = detTinv
19  if( nverts != 4 )
20  {
21  std::cout << "Invalid Tetrahedron. Expected 4 vertices.\n";
22  return MB_FAILURE;
23  }
24 
25  assert( verts );
26 
27  Matrix3 J( verts[1 * 3 + 0] - verts[0 * 3 + 0], verts[2 * 3 + 0] - verts[0 * 3 + 0],
28  verts[3 * 3 + 0] - verts[0 * 3 + 0], verts[1 * 3 + 1] - verts[0 * 3 + 1],
29  verts[2 * 3 + 1] - verts[0 * 3 + 1], verts[3 * 3 + 1] - verts[0 * 3 + 1],
30  verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2],
31  verts[3 * 3 + 2] - verts[0 * 3 + 2] );
32 
33  // Update the work array
34  if( !work ) work = new double[20];
35 
36  J.copyto( work );
37  J.inverse().copyto( work + Matrix3::size );
38  work[18] = J.determinant();
39  work[19] = ( work[18] < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / work[18] );
40 
41  return MB_SUCCESS;
42 }
43 
44 ErrorCode LinearTet::evalFcn( const double* params,
45  const double* field,
46  const int /*ndim*/,
47  const int num_tuples,
48  double* /*work*/,
49  double* result )
50 {
51  assert( params && field && num_tuples > 0 );
52  std::vector< double > f0( num_tuples );
53  std::copy( field, field + num_tuples, f0.begin() );
54  std::copy( field, field + num_tuples, result );
55 
56  for( unsigned i = 1; i < 4; ++i )
57  {
58  double p = 0.5 * ( params[i - 1] + 1 ); // transform from -1 <= p <= 1 to 0 <= p <= 1
59  for( int j = 0; j < num_tuples; j++ )
60  result[j] += ( field[i * num_tuples + j] - f0[j] ) * p;
61  }
62 
63  return MB_SUCCESS;
64 }
65 
66 ErrorCode LinearTet::integrateFcn( const double* field,
67  const double* /*verts*/,
68  const int nverts,
69  const int /*ndim*/,
70  const int num_tuples,
71  double* work,
72  double* result )
73 {
74  assert( field && num_tuples > 0 );
75  std::fill( result, result + num_tuples, 0.0 );
76  for( int i = 0; i < nverts; ++i )
77  {
78  for( int j = 0; j < num_tuples; j++ )
79  result[j] += field[i * num_tuples + j];
80  }
81  double tmp = work[18] / 24.0;
82  for( int i = 0; i < num_tuples; i++ )
83  result[i] *= tmp;
84 
85  return MB_SUCCESS;
86 }
87 
88 ErrorCode LinearTet::jacobianFcn( const double*, const double*, const int, const int, double* work, double* result )
89 {
90  // jacobian is cached in work array
91  assert( work );
92  std::copy( work, work + 9, result );
93  return MB_SUCCESS;
94 }
95 
96 ErrorCode LinearTet::reverseEvalFcn( EvalFcn eval,
97  JacobianFcn jacob,
98  InsideFcn ins,
99  const double* posn,
100  const double* verts,
101  const int nverts,
102  const int ndim,
103  const double iter_tol,
104  const double inside_tol,
105  double* work,
106  double* params,
107  int* is_inside )
108 {
109  assert( posn && verts );
110  return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
111  is_inside );
112 }
113 
114 int LinearTet::insideFcn( const double* params, const int, const double tol )
115 {
116  return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[2] >= -1.0 - tol &&
117  params[0] + params[1] + params[2] <= 1.0 + tol );
118 }
119 
120 ErrorCode LinearTet::evaluate_reverse( EvalFcn eval,
121  JacobianFcn jacob,
122  InsideFcn inside_f,
123  const double* posn,
124  const double* verts,
125  const int nverts,
126  const int ndim,
127  const double iter_tol,
128  const double inside_tol,
129  double* work,
130  double* params,
131  int* inside )
132 {
133  // TODO: should differentiate between epsilons used for
134  // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
135  // right now, fix the tolerance used for NR
136  const double error_tol_sqr = iter_tol * iter_tol;
137  CartVect* cvparams = reinterpret_cast< CartVect* >( params );
138  const CartVect* cvposn = reinterpret_cast< const CartVect* >( posn );
139 
140  // find best initial guess to improve convergence
141  CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ),
142  CartVect( -1, -1, 1 ) };
143  double resl = std::numeric_limits< double >::max();
144  CartVect new_pos, tmp_pos;
145  ErrorCode rval;
146  for( unsigned int i = 0; i < 4; i++ )
147  {
148  rval = ( *eval )( tmp_params[i].array(), verts, ndim, ndim, work, tmp_pos.array() );
149  if( MB_SUCCESS != rval ) return rval;
150  double tmp_resl = ( tmp_pos - *cvposn ).length_squared();
151  if( tmp_resl < resl )
152  {
153  *cvparams = tmp_params[i];
154  new_pos = tmp_pos;
155  resl = tmp_resl;
156  }
157  }
158 
159  // residual is diff between old and new pos; need to minimize that
160  CartVect res = new_pos - *cvposn;
161  Matrix3 J;
162  rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() );
163 #ifndef NDEBUG
164  double det = J.determinant();
165  assert( det > std::numeric_limits< double >::epsilon() );
166 #endif
167  Matrix3 Ji = J.inverse();
168 
169  int iters = 0;
170  // while |res| larger than tol
171  int dum, *tmp_inside = ( inside ? inside : &dum );
172  while( res % res > error_tol_sqr )
173  {
174  if( ++iters > 25 )
175  {
176  // if we haven't converged but we're outside, that's defined as success
177  *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
178  if( !( *tmp_inside ) )
179  return MB_SUCCESS;
180  else
181  return MB_INDEX_OUT_OF_RANGE;
182  }
183 
184  // new params tries to eliminate residual
185  *cvparams -= Ji * res;
186 
187  // get the new forward-evaluated position, and its difference from the target pt
188  rval = ( *eval )( params, verts, ndim, ndim, work, new_pos.array() );
189  if( MB_SUCCESS != rval ) return rval;
190  res = new_pos - *cvposn;
191  }
192 
193  if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );
194 
195  return MB_SUCCESS;
196 } // Map::evaluate_reverse()
197 
198 ErrorCode LinearTet::normalFcn( const int ientDim,
199  const int facet,
200  const int nverts,
201  const double* verts,
202  double normal[3] )
203 {
204  // assert(facet < 4 && ientDim == 2 && nverts == 4);
205  if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed tet :: expected value = 4 " );
206  if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " );
207  if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-3" );
208 
209  int id0 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][0];
210  int id1 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][1];
211  int id2 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][2];
212 
213  double x0[3], x1[3];
214 
215  for( int i = 0; i < 3; i++ )
216  {
217  x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
218  x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i];
219  }
220 
221  double a = x0[1] * x1[2] - x1[1] * x0[2];
222  double b = x1[0] * x0[2] - x0[0] * x1[2];
223  double c = x0[0] * x1[1] - x1[0] * x0[1];
224  double nrm = sqrt( a * a + b * b + c * c );
225 
226  if( nrm > std::numeric_limits< double >::epsilon() )
227  {
228  normal[0] = a / nrm;
229  normal[1] = b / nrm;
230  normal[2] = c / nrm;
231  }
232  return MB_SUCCESS;
233 }
234 
235 } // namespace moab