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
get_h_(double xtmp[][4],double u[],double * h)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
length(double x[4][3],double u[])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
intersect(double x[4][3],double u[],double points[4][3])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
u_line(double s,double xc[3],double udir[3],double pnt[3])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
diffVt(double vec1[3],double vec2[3],double ans[3])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
crossProd(double a[3],double b[3],double ans[3])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
dotProd(double a[3],double b[3])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
vecMag(double a[3])196 double vecMag(double a[3]){
197 return sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
198 }
199
vecEqual(double a[3],double b[3])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
isPntInTri(double t_xyz[3][3],double t_normal[3],double p_xyz[3],double tol)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