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