1 /* fortran wrapper for C program length which returns the length 2 of an element in a given direction. 3 */ 4 #include "vec_func.h" 5 6 7 double length(double [][3],double []); 8 int intersect(double x[4][3],double u[],double points[4][3] ); 9 int isPntInTri(double [][3],double [], double [], double ); 10 11 #ifdef sun4_5 12 get_h_(double xtmp[][4], double u[], double *h) 13 #endif 14 #ifdef LINUX 15 get_h_(double xtmp[][4], double u[], double *h) 16 #endif 17 #ifdef ibm 18 get_h(double xtmp[][4], double u[], double *h) 19 #endif 20 #ifdef sgi 21 void get_h(double xtmp[][4], double u[], double *h) 22 #endif 23 #ifdef decalp 24 void get_h(double xtmp[][4], double u[], double *h) 25 #endif 26 #ifdef intel 27 void GET_H(double xtmp[][4], double u[], double *h) 28 #endif 29 { 30 double x[4][3]; 31 int i,j; 32 /*recall that arrays in fortran and C are transposed */ 33 for(i=0; i < 4; i++) 34 for(j=0; j < 3; j++) 35 x[i][j] = xtmp[j][i]; 36 37 *h = length(x,u); 38 } 39 40 double length(double x[4][3],double u[] ){ 41 int nint; 42 double tmp[3],he; 43 double points[4][3]; 44 45 /* compute the intersection points. nint is returned 46 as the number of faces that this line intersected 47 within the tet */ 48 nint = intersect(x,u,points); 49 50 /* if nint is larger than one, then the vector intersected 51 an edge (nint=3) or corner (nint=4) of the tet. In this 52 case, we need to make sure that we are computing the 53 length between two different points. */ 54 55 if(nint > 1){ /* vecEqual returns 1 if 56 these two vectors are the same. 57 */ 58 if( vecEqual(points[0],points[1]) ){ 59 if( vecEqual(points[0],points[2]) ){ 60 diffVt(points[0],points[3],tmp); 61 } 62 else{ 63 diffVt(points[0],points[2],tmp); 64 } 65 } 66 else{ 67 diffVt(points[0],points[1],tmp); 68 } 69 } 70 else{ 71 diffVt(points[0],points[1],tmp); 72 } 73 74 /* compute the length of the element */ 75 he = vecMag(tmp); 76 return he; 77 } 78 79 /* 80 This function returns the intersection points of a line 81 through the centroid of a tet in a given direction. 82 */ 83 84 85 /* u_line is a function that defines the parametric equation 86 of a line given a point on the lineg, a direction, and a 87 parameter value */ 88 89 void u_line(double s, double xc[3], double udir[3], double pnt[3]); 90 91 /* 92 All points which intersect the tet are returned, however 93 only two of them need be used. The value of the function on 94 exit is the number of points of intersection, ie. the number 95 of faces not parallel to udir */ 96 97 int intersect(double x[4][3],double u[],double points[4][3] ){ 98 int pos = 0; 99 double tol = 0.000001; 100 double eps = 0.000001; 101 int face,i,j,k; 102 double xc[3]; 103 double xbar[3][3]; 104 double xref[3],v1[3],v2[3],normal[3]; 105 /* define local vertex numbers of four faces of the tet */ 106 int f_index[4][3] = {{0,1,2},{1,3,2},{0,2,3},{0,3,1}}; 107 double n_dot_u; 108 double xr_xc[3], s, pt[3]; 109 110 /* compute centroid coordinates */ 111 for(i=0; i<3; i++) 112 xc[i] = 0.25*(x[0][i]+x[1][i]+x[2][i]+x[3][i]); 113 114 /* loop over each face */ 115 for( face= 0 ; face< 4; face++ ){ 116 /* compute coordinates of this face */ 117 for( j= 0 ; j<3 ; j++) 118 for( k= 0; k< 3; k++) 119 xbar[j][k] = x[f_index[face][j]][k]; 120 121 /* set reference coordinate */ 122 if( face != 2 ){ 123 for( i= 0; i<3; i++ ) 124 xref[i] = xbar[0][i]; 125 } 126 else{ /* face 2 is the only face that 127 can't use xbar[0] as a reference */ 128 for( i= 0; i<3; i++ ) 129 xref[i] = xbar[1][i]; 130 } 131 132 /* compute two vectors which describe the plane of this 133 face */ 134 for( i=0; i< 3; i++){ 135 v1[i] = xbar[1][i] - xbar[0][i]; 136 v2[i] = xbar[2][i] - xbar[0][i]; 137 } 138 /* get the normal to this plane */ 139 crossProd(v1,v2,normal); 140 141 n_dot_u = dotProd(normal,u); 142 143 /* if u is not parallel to this face, then find 144 the point of intersection */ 145 if( fabs(n_dot_u) > eps ){ 146 diffVt( xref, xc, xr_xc ); 147 148 /* s is the parameter value where the line parallel to 149 to the vector u intersects this face*/ 150 s = dotProd(normal,xr_xc)/n_dot_u; 151 152 u_line(s,xc,u,pt); /* find intersection point */ 153 154 /* if this point lies on the tet, then save it */ 155 if( isPntInTri(xbar,normal,pt,tol) ){ 156 for(i=0; i<3; i++) 157 points[pos][i] = pt[i]; 158 pos += 1; 159 } 160 } 161 else{ /* do nothing */ 162 163 } 164 } 165 166 return pos; 167 } 168 169 /* parametric equation of the line passing through xc 170 in the direction of u. 171 */ 172 173 void u_line(double s, double xc[3], double udir[3], double pnt[3]){ 174 int i; 175 for(i=0; i<3; i++) 176 pnt[i] = xc[i]+s*udir[i]; 177 } 178 179 180 void diffVt(double vec1[3],double vec2[3], double ans[3]){ 181 ans[0] = vec1[0]-vec2[0]; 182 ans[1] = vec1[1]-vec2[1]; 183 ans[2] = vec1[2]-vec2[2]; 184 } 185 186 void crossProd(double a[3], double b[3], double ans[3]){ 187 ans[0] = a[1]*b[2]-b[1]*a[2]; 188 ans[1] = a[2]*b[0]-b[2]*a[0]; 189 ans[2] = a[0]*b[1]-b[0]*a[1]; 190 } 191 192 double dotProd(double a[3], double b[3]){ 193 return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; 194 } 195 196 double vecMag(double a[3]){ 197 return sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]); 198 } 199 200 int vecEqual(double a[3], double b[3]){ 201 return (a[0]==b[0] && a[1]==b[1] && a[2]==b[2]); 202 } 203 204 /* 205 Determines if a point is inside or out of a triangle 206 Pt will be declared in triangle if distance to side < tol 207 if u want the tolerancing to be relative to 't_xyz', 208 pass tol (absolute tolerance) * size of 't_xyz' 209 */ 210 211 int isPntInTri( 212 double t_xyz[3][3], 213 double t_normal[3], /* Normal to triangle (not normed) */ 214 double p_xyz[3], 215 double tol 216 ) 217 218 { 219 220 double cross_p[3]; 221 double d; 222 int side; 223 int status; 224 double cross_p_sqnorm; 225 double dnbr; 226 double dtol; 227 double t_vec[3]; 228 double t_vec_norm; 229 double tp_vec[3]; 230 double tp_vec_sqnorm; 231 232 status= 1; 233 234 for ( side= 0 ; side< 3 ; side++ ) { 235 236 diffVt(t_xyz[(side+1)%3],t_xyz[side],t_vec); 237 crossProd(t_normal,t_vec,cross_p); 238 cross_p_sqnorm= dotProd(cross_p,cross_p); 239 240 diffVt(p_xyz,t_xyz[side],tp_vec); 241 242 d= dotProd(cross_p,tp_vec); 243 dnbr= d * d; 244 dtol= tol*tol*cross_p_sqnorm; 245 246 if ( dnbr < dtol ) { 247 248 /* 249 Pt is considered to be on the side (extended to infinity) 250 */ 251 252 /* 253 Make sure distance to point 0 or 1 less than ... 254 */ 255 256 t_vec_norm= vecMag(t_vec); 257 dtol= t_vec_norm + tol; 258 dtol *= dtol; 259 260 tp_vec_sqnorm= dotProd(tp_vec,tp_vec); 261 if ( tp_vec_sqnorm > dtol ) { 262 status= 0; 263 break; 264 } 265 266 diffVt(p_xyz,t_xyz[(side+1)%3],tp_vec); 267 tp_vec_sqnorm= dotProd(tp_vec,tp_vec); 268 if ( tp_vec_sqnorm > dtol ) { 269 status= 0; 270 break; 271 } 272 273 continue; 274 } 275 276 if ( d > 0.0 ) { 277 278 /* 279 Pt is strictly in the interior of the side 280 */ 281 } 282 else if ( d < 0.0 ) { 283 284 /* 285 Pt is strictly in the outisde of the side 286 */ 287 288 status= 0; 289 break; 290 } 291 292 } 293 294 return status; 295 296 } 297 298 299 300 301 302 303 304 305 306 307 308