1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1330.c,v 1.3 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1330
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1330(double epar11[],double epar12[],double epar21[],double epar22[],double eval11[],double eval12[],double eval21[],double eval22[],int * jbound,double gpar1[],double gpar2[],int * jstat)55 s1330(double epar11[],double epar12[],double epar21[],double epar22[],
56 double eval11[],double eval12[],double eval21[],double eval22[],
57 int *jbound,double gpar1[],double gpar2[],int *jstat)
58 #else
59 void s1330(epar11,epar12,epar21,epar22,eval11,eval12,eval21,eval22,
60 jbound,gpar1,gpar2,jstat)
61 double epar11[];
62 double epar12[];
63 double epar21[];
64 double epar22[];
65 double eval11[];
66 double eval12[];
67 double eval21[];
68 double eval22[];
69 int *jbound;
70 double gpar1[];
71 double gpar2[];
72 int *jstat;
73 #endif
74 /*
75 *********************************************************************
76 *
77 * PURPOSE : To find if there is an intersection between epar1 and
78 * epar2 with the 4-D SISLbox desribed by eval11[0]:eval11[1]
79 * in the first parameter direction and eval12[0]:eval12[1]
80 * in the second parameter direction in the first patch, and
81 * eval21[0]:eval21[1] in the first parameter direction and
82 * eval22[0]:eval22[1] in the second parameter direction in
83 * the second patch. If there is an intersection find the
84 * intersection closest to epar1.
85 *
86 * INPUT : epar11 - Start of line in first surface
87 * epar12 - Start of line in second surface
88 * epar21 - End of line in first surface
89 * epar22 - End of line in second surface
90 * eval11 - Interval in first parameter direction in patch 1
91 * eval12 - Interval in second parameter direction in patch 1
92 * eval21 - Interval in first parameter direction in patch 2
93 * eval22 - Interval in second parameter direction in patch 2
94 *
95 *
96 * OUTPUT : gpar1 - Parameter pair of intersection in first surface
97 * : gpar2 - Parameter pair of intersection in second surface
98 * jbound - Indicator telling along which boundary
99 * we have an intersection
100 * = 0 : no intersection
101 * = 1 : intersection along u=eval11[0]
102 * = 2 : intersection along v=eval12[1]
103 * = 3 : intersection along u=eval11[1]
104 * = 4 : intersection along v=eval12[0]
105 * = 5 : intersection along s=eval21[0]
106 * = 6 : intersection along t=eval22[1]
107 * = 7 : intersection along s=eval21[1]
108 * = 8 : intersection along t=eval22[0]
109 * jstat - status messages
110 * = 0 : Line outside no intersection
111 * = 1 : Line inside no intersection
112 * = 2 : epar2 outside, epar1 inside, step out
113 * = 3 : epar1 outside, epar2 inside, step in
114 * = 4 : epar2 outside, epar1 on boundary
115 * < 0 : error
116 *
117 *
118 * METHOD :
119 *
120 * REFERENCES :
121 *
122 *-
123 * CALLS :
124 *
125 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway, 16-August-1988
126 * Revised by : Tor Dokken, SI, Oslo, Norway, 3-april-1989
127 * Correction of long steps
128 *
129 *********************************************************************
130 */
131 {
132 int kstat1=0,kstat2=0; /* Local status variable */
133 int kstat=0; /* Local status variable */
134 int kpos=0; /* Position of error */
135 int kins1; /* epar1 inside/outside SISLbox */
136 int kins2; /* epar2 inside/outside SISLbox */
137 int kbound1; /*Intersection indicator along boundary of surface 1*/
138 int kbound2; /*Intersection indicator along boundary of surface 2*/
139
140 double tdom; /* Denominator of last intersection point */
141 double tfak1,tfak2; /* Distance to straight line */
142 double spar11[2]; /*Candidate intersection point two first coordinates*/
143 double spar12[2]; /*Candidate intersection point two last coordinates */
144 double spar21[2]; /*Candidate intersection point two first coordinates*/
145 double spar22[2]; /*Candidate intersection point two last coordinates */
146
147 *jbound = 0;
148
149 /* Test if both ends are inside */
150
151 kins1 = kins2 = 0;
152
153 if (eval11[0] <= epar11[0]+REL_PAR_RES && epar11[0] <= eval11[1]+REL_PAR_RES &&
154 eval12[0] <= epar11[1]+REL_PAR_RES && epar11[1] <= eval12[1]+REL_PAR_RES &&
155 eval21[0] <= epar12[0]+REL_PAR_RES && epar12[0] <= eval21[1]+REL_PAR_RES &&
156 eval22[0] <= epar12[1]+REL_PAR_RES && epar12[1] <= eval22[1]+REL_PAR_RES)
157 kins1 = 1;
158
159 if (eval11[0] <= epar21[0]+REL_PAR_RES && epar21[0] <= eval11[1]+REL_PAR_RES &&
160 eval12[0] <= epar21[1]+REL_PAR_RES && epar21[1] <= eval12[1]+REL_PAR_RES &&
161 eval21[0] <= epar22[0]+REL_PAR_RES && epar22[0] <= eval21[1]+REL_PAR_RES &&
162 eval22[0] <= epar22[1]+REL_PAR_RES && epar22[1] <= eval22[1]+REL_PAR_RES)
163 kins2 = 1;
164
165
166 /* Test if we step from the boundary and out */
167
168 /* if ((eval11[0] == epar11[0] && epar21[0] < eval11[0]) || */
169 /* (epar11[0] == eval11[1] && eval11[1] < epar21[0]) || */
170 /* (eval12[0] == epar11[1] && epar21[1] < eval12[0]) || */
171 /* (epar11[1] == eval12[1] && eval12[1] < epar21[1]) || */
172 /* (eval21[0] == epar12[0] && epar22[0] < eval21[0]) || */
173 /* (epar12[0] == eval21[1] && eval21[1] < epar22[0]) || */
174 /* (eval22[0] == epar12[1] && epar22[1] < eval22[0]) || */
175 /* (epar12[1] == eval22[1] && eval22[1] < epar22[1])) goto war04; */
176 if ((DEQUAL(eval11[0],epar11[0]) && epar21[0] < eval11[0]) ||
177 (DEQUAL(epar11[0],eval11[1]) && eval11[1] < epar21[0]) ||
178 (DEQUAL(eval12[0],epar11[1]) && epar21[1] < eval12[0]) ||
179 (DEQUAL(epar11[1],eval12[1]) && eval12[1] < epar21[1]) ||
180 (DEQUAL(eval21[0],epar12[0]) && epar22[0] < eval21[0]) ||
181 (DEQUAL(epar12[0],eval21[1]) && eval21[1] < epar22[0]) ||
182 (DEQUAL(eval22[0],epar12[1]) && epar22[1] < eval22[0]) ||
183 (DEQUAL(epar12[1],eval22[1]) && eval22[1] < epar22[1]))
184 goto war04;
185
186 if (kins1==1 && kins2==1)
187 goto war01;
188
189 /* Test if both ends are to the left, right, below or above */
190
191 if ( (epar11[0] < eval11[0] && epar21[0] < eval11[0]) ||
192 (eval11[1] < epar11[0] && eval11[1] < epar21[0] ) ||
193 (epar11[1] < eval12[0] && epar21[1] < eval12[0]) ||
194 (eval12[1] < epar11[1] && eval12[1] < epar21[1] ) ||
195 (epar12[0] < eval21[0] && epar22[0] < eval21[0]) ||
196 (eval21[1] < epar12[0] && eval21[1] < epar22[0] ) ||
197 (epar12[1] < eval22[0] && epar22[1] < eval22[0]) ||
198 (eval22[1] < epar12[1] && eval22[1] < epar22[1] ) )
199 goto war00;
200
201
202
203 /* Check if intersection in first two dimensions */
204
205 s1305(epar11,epar21,eval11,eval12,&kbound1,spar11,&kstat);
206
207 if (kstat<0) goto error;
208 kstat1 = kstat;
209 if (kstat1==0) goto war00;
210
211 /* Calculate two last coefficients */
212
213 if (kstat1==2 || kstat1==3)
214 {
215 tfak1 = fabs(spar11[0]-epar11[0]) + fabs(spar11[1]-epar11[1]);
216 tfak2 = fabs(epar21[0]-spar11[0]) + fabs(epar21[1]-spar11[1]);
217 tdom = tfak1 + tfak2;
218 if (DNEQUAL(tdom,DZERO))
219 {
220 spar12[0] = (tfak2*epar12[0] + tfak1*epar22[0])/tdom;
221 spar12[1] = (tfak2*epar12[1] + tfak1*epar22[1])/tdom;
222
223 /* If the two last coefficients are zero, then then this intersection
224 point must be discarded */
225
226 if (spar12[0]<eval21[0] || eval21[1]<spar12[0] ||
227 spar12[1]<eval22[0] || eval22[1]<spar12[1])
228 {
229 /* Intersection point outside */
230
231 kbound1 = 0;
232 }
233 }
234 else
235 {
236 /* epar1, spar and epar2 has equal first coordinates, since there
237 is an intersection all must lie on the boundary, thus all are
238 inside */
239
240 kbound1 = 0;
241 }
242 }
243 else if (kstat1==4 && kins1==1)
244 {
245 /* On boundary and stepping out */
246 goto war04;
247 }
248
249 /* Check if intersection in last two dimensions */
250
251 s1305(epar12,epar22,eval21,eval22,&kbound2,spar22,&kstat);
252
253 if (kstat<0)
254 goto error;
255 kstat2 = kstat;
256 if (kstat2==0)
257 goto war00;
258
259 if (kstat1==1 && kstat2==1)
260 goto war01;
261
262
263 /* Calculate two last coefficients */
264
265 if (kstat2==2 || kstat2==3)
266 {
267 tfak1 = fabs(spar22[0]-epar12[0]) + fabs(spar22[1]-epar12[1]);
268 tfak2 = fabs(epar22[0]-spar22[0]) + fabs(epar22[1]-spar22[1]);
269 tdom = tfak1 + tfak2;
270 if (DNEQUAL(tdom,DZERO))
271 {
272 spar21[0] = (tfak2*epar11[0] + tfak1*epar21[0])/tdom;
273 spar21[1] = (tfak2*epar11[1] + tfak1*epar21[1])/tdom;
274
275 /* If the two last coefficients are zero, then then this intersection
276 point must be discarded */
277
278 if (spar21[0]<eval11[0] || eval11[1]<spar21[0] ||
279 spar21[1]<eval12[0] || eval12[1]<spar21[1])
280 {
281 /* Intersection point outside */
282
283 kbound2 = 0;
284 }
285 }
286 else
287 {
288 /* epar1, spar and epar2 has equal last coordinates, since there
289 is an intersection all must lie on the boundary, thus all are
290 inside */
291
292 kbound2 = 0;
293 }
294 }
295 else if (kstat2==4 && kins1==1)
296 {
297 /* On boundary and stepping out */
298 goto war04;
299 }
300
301 /* kbound1 and kbound2 tells if we have got and intersection with the
302 boundary */
303
304
305 /* If intersections along both boundaries then find which intersection
306 is closest to epar1 */
307
308 if (kbound1!=0 && kbound2!=0)
309 {
310 /*guen int t1,t2,t3,t4;*/ /* temporary varuiables */
311 /*guen changed to */
312 double t1,t2,t3,t4; /* temporary variables */
313
314 t1 = s6dist(spar11,epar11,2);
315 t2 = s6dist(spar12,epar12,2);
316 t3 = s6dist(spar21,epar11,2);
317 t4 = s6dist(spar22,epar12,2);
318
319 if ((t1*t1+t2*t2) < (t3*t3+t4*t4) )
320 kbound2 = 0;
321 else
322 kbound1 = 0;
323 }
324
325 if (kbound1==0 && kbound2 ==0)
326 {
327 /* No intersection */
328 goto war00;
329 }
330 else if (kbound1!=0 && kbound2==0)
331 {
332 /* Intersection with boundary of first patch */
333 memcopy(gpar1,spar11,2,DOUBLE);
334 memcopy(gpar2,spar12,2,DOUBLE);
335 *jbound = kbound1;
336 }
337 else if (kbound1==0 && kbound2!=0)
338 {
339 /* Intersection with boundary of second patch */
340 memcopy(gpar1,spar21,2,DOUBLE);
341 memcopy(gpar2,spar22,2,DOUBLE);
342 *jbound = kbound2+4;
343 }
344
345 if (kins1 == 1)
346 {
347 if (eval11[0] == epar11[0] || epar11[0] == eval11[1] ||
348 eval12[0] == epar11[1] || epar11[1] == eval12[1] ||
349 eval21[0] == epar12[0] || epar12[0] == eval21[1] ||
350 eval22[0] == epar12[1] || epar12[1] == eval22[1])
351 {
352 goto war04;
353 }
354 else
355 {
356 goto war02;
357 }
358 }
359
360 if (kins2 == 1)
361 goto war03;
362
363 goto war05;
364
365 /* Line outside */
366
367 war00:
368 *jstat = 0;
369 goto out;
370
371 /* Line inside */
372
373 war01:
374 *jstat = 1;
375 goto out;
376
377 /* epar1 inside epar2 outside */
378 war02:
379 *jstat = 2;
380 goto out;
381
382 /* epar2 inside epar1 outside */
383 war03:
384 *jstat = 3;
385 goto out;
386
387 /* epar1 on boundary, epar2 outside */
388 war04:
389 *jstat = 4;
390 goto out;
391
392 /* Special error */
393 war05:
394 *jstat = 5;
395 goto out;
396
397 /* Error in lower leve function */
398 error:
399 *jstat = kstat;
400 s6err("s1330",*jstat,kpos);
401 goto out;
402
403
404 out:
405 return;
406 }
407