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: sh1782.c,v 1.3 2005-02-28 09:04:50 afr Exp $
45  *
46  */
47 
48 
49 #define SH1782
50 
51 #include "sislP.h"
52 
53 
54 
55 #if defined(SISLNEEDPROTOTYPES)
56 static void
57 sh1782_s9sf_pt (SISLObject *, SISLObject *, double, SISLIntdat **,
58 		SISLIntpt **, int, int, int *);
59 static void
60 sh1782_s9sf_cu (SISLObject *, SISLObject *, double, SISLIntdat **,
61 		SISLIntpt **, int, int, int *);
62 static void
63 sh1782_s9sf_sf (SISLObject *, SISLObject *, double, SISLIntdat **,
64 		SISLIntpt **, int, int, int *);
65 #else
66 static void sh1782_s9sf_pt ();
67 static void sh1782_s9sf_cu ();
68 static void sh1782_s9sf_sf ();
69 #endif
70 
71 #if defined(SISLNEEDPROTOTYPES)
72 void
sh1782(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat * pintdat,int ipar,double apar,SISLIntdat ** rintdat,int * jnewpt,int * jstat)73 sh1782 (SISLObject * po1, SISLObject * po2, double aepsge,
74 	SISLIntdat * pintdat, int ipar, double apar,
75 	SISLIntdat ** rintdat, int *jnewpt, int *jstat)
76 #else
77 void
78 sh1782 (po1, po2, aepsge, pintdat, ipar, apar, rintdat, jnewpt, jstat)
79      SISLObject *po1;
80      SISLObject *po2;
81      double aepsge;
82      SISLIntdat *pintdat;
83      int ipar;
84      double apar;
85      SISLIntdat **rintdat;
86      int *jnewpt;
87      int *jstat;
88 #endif
89 /*
90 *********************************************************************
91 *
92 *********************************************************************
93 *
94 * PURPOSE    : Insert all points in one intersection data structure
95 *              into another intersection data structure with one more
96 *              parameter direction.
97 *
98 *
99 * INPUT      : po1      - Pointer to the first object in the intersection.
100 *              po2      - Pointer to the second object in the intersection.
101 *              aepsge   - Geometry resolution.
102 *              pintdat  - Intersection data structure of a problem in a
103 *                         reduced parameter space.
104 *              ipar     - Missing parameter direction in the reduced space.
105 *              apar     - Value of the missing parameter for all intersection
106 *                         points.
107 *
108 *
109 * OUTPUT     : rintdat  - Intersection data of object-object intersection.
110 *              jnewpt   - Number of new intersection points.
111 *              jstat    - status messages
112 *                                > 0   : Warning.
113 *                                = 0   : Ok.
114 *                                < 0   : Error.
115 *
116 *
117 * METHOD     :
118 *
119 * CALLS      : sh1782_s9sf_pt  - Set edge info 1D and 2D surface-point.
120 *              sh1782_s9sf_cu  - Set edge info 3D surface-curve.
121 *              sh1782_s9sf_sf  - Set edge info 3D surface-surface.
122 *              sh6idput        - Put point into higher dimensional parameter
123 *                                space.
124 *
125 * REFERENCES :
126 *
127 * WRITTEN BY : Ulf J. Krystad, SI, 06.91
128 *********************************************************************
129 */
130 {
131   int kstat = 0;		/* Status variable.                     */
132   int kdim;			/* Dimension of geometry space.         */
133   int knpoint;			/* Number of int.pt. in array.          */
134   SISLIntpt **uintpt = SISL_NULL;	/* Array storing intersection points.   */
135 
136   *jnewpt = 0;
137 
138   /* Test if an intersection data structure exist.  */
139 
140   if (pintdat == SISL_NULL)
141     {
142       *jstat = 0;
143       goto out;
144     }
145 
146   /* Fetch dimension of geometry space. */
147 
148   if (po1->iobj == SISLPOINT)
149     kdim = po1->p1->idim;
150   else if (po1->iobj == SISLCURVE)
151     kdim = po1->c1->idim;
152   else
153     kdim = po1->s1->idim;
154 
155   /* Insert the intersection points from pintdat into the structure
156      *rintdat belonging to the current problem.  */
157 
158   sh6idput (po1, po2, rintdat, pintdat, ipar, apar, &uintpt, &knpoint, &kstat);
159   if (kstat < 0)
160     goto error;
161 
162   if (knpoint == 0)
163     {
164       *jstat = 0;
165       goto out;
166     }
167 
168 
169   /* Browse on the dimension of geometry space and the type of
170      the input objects.     */
171 
172   if (kdim <= 2 && ((po1->iobj == SISLSURFACE && po2->iobj == SISLPOINT)
173 		    || (po2->iobj == SISLSURFACE && po1->iobj == SISLPOINT)))
174     {
175       /* Compute pre-topology in one-dimensional surface-level value
176          intersection.           */
177 
178       sh1782_s9sf_pt (po1, po2, aepsge, rintdat, uintpt, knpoint,
179 		      ipar, &kstat);
180       if (kstat < 0)
181 	goto error;
182 
183     }
184 
185   else if (kdim == 3 &&
186 	 ((po1->iobj == SISLSURFACE && po2->iobj == SISLCURVE) ||
187 	  (po1->iobj == SISLCURVE && po2->iobj == SISLSURFACE )))
188     {
189       /* Surface-surface intersection in 3-dimensional geometry space. */
190 
191       sh1782_s9sf_cu (po1, po2, aepsge, rintdat, uintpt, knpoint,
192 		      ipar, &kstat);
193       if (kstat < 0)
194 	goto error;
195     }
196 
197   else if (kdim == 3 && po1->iobj == SISLSURFACE
198 	   && po2->iobj == SISLSURFACE)
199     {
200       /* Surface-surface intersection in 3-dimensional geometry space. */
201 
202       sh1782_s9sf_sf (po1, po2, aepsge, rintdat, uintpt, knpoint,
203 		      ipar, &kstat);
204       if (kstat < 0)
205 	goto error;
206     }
207 
208   /* Task performed.  */
209 
210   *jstat = 0;
211   goto out;
212 
213   /* Error in lower level routine.  */
214 
215 error:*jstat = kstat;
216   goto out;
217 
218 out:if (uintpt != SISL_NULL)
219     freearray (uintpt);
220 }
221 
222 
223 
224 
225 #if defined(SISLNEEDPROTOTYPES)
226 static void
sh1782_s9sf_pt(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat ** rintdat,SISLIntpt ** uintpt,int kpoint,int ipar,int * jstat)227 sh1782_s9sf_pt (SISLObject * po1, SISLObject * po2, double aepsge,
228 		SISLIntdat ** rintdat, SISLIntpt ** uintpt, int kpoint,
229 		int ipar, int *jstat)
230 #else
231 static void
232 sh1782_s9sf_pt (po1, po2, aepsge, rintdat, uintpt, kpoint, ipar, jstat)
233      SISLObject *po1;
234      SISLObject *po2;
235      double aepsge;
236      SISLIntdat **rintdat;
237      SISLIntpt **uintpt;
238      int kpoint;
239      int ipar;
240      int *jstat;
241 #endif
242 /*
243 *********************************************************************
244 *
245 *********************************************************************
246 *
247 * PURPOSE    : Set pre-topology data in 1-dimensional surface-point
248 *              intersection.
249 *
250 *
251 * INPUT      : po1      - Pointer to the first object in the intersection.
252 *              po2      - Pointer to the second object in the intersection.
253 *              aepsge   - Geometry resolution.
254 *              rintdat  - Intersection data structure.
255 *              uintpt   - Intersection points.
256 *              knpoint  - No of int points
257 *              ipar     - Parameter direction of missing parameter in the
258 *                         lower dimensional parameter space.
259 *
260 * OUTPUT     : uintpt   - Updated intersection point.
261 *              jstat    - status messages
262 *                                > 0   : Warning.
263 *                                = 0   : Ok.
264 *                                < 0   : Error.
265 *
266 *
267 *
268 * METHOD     :
269 *
270 * CALLS      : sh6getgeom -  Get geometry of intersection point.
271 *              sh6gettop  -  Get topology of intersection point.
272 *              sh6settop  -  Set topology of intersection point.
273 * REFERENCES :
274 *
275 * WRITTEN BY : Ulf J. Krystad, SI, Oslo, Norway. 07.91.
276 *********************************************************************
277 */
278 {
279   int kstat = 0;		/* Status variable.                        */
280   int kdim;			/* Dimension of geometry space.            */
281   double tdum, tsign;		/* Parameters used to browse between cases.*/
282   double *sder;			/* Surface geometry values.                */
283   double *sdum;			/* Normal of surface(dummy)                */
284   int ind;			/* Index for curve's par val.              */
285   int ki, kj, klow, khigh;	/* Help indexes into uintpt                */
286   int index1, index2;		/* Dummy in this context                   */
287   int kant;			/* Edge indicator                          */
288   int *edge_f, *edge_l;		/* Pointers to atribute edge in intpt      */
289   SISLSurf *ps;			/* The 1D Surface                          */
290   /* --------------------------------------------------------------------- */
291 
292   /* Test input */
293   kdim = (po1->iobj == SISLSURFACE) ? po1->s1->idim : po2->s1->idim;
294   if (kdim != 1 && kdim !=2)
295     goto err106;
296 
297   if (ipar == 0)
298     ind = 1;
299   else
300     ind = 0;
301 
302   if (po1->iobj == SISLSURFACE)
303     tsign = (double) 1.0;
304   else
305     tsign = -(double) 1.0;
306 
307   /* Loop for all points */
308   for (ki = 0; ki < kpoint; ki++)
309     for (kj = ki + 1; kj < kpoint; kj++)
310       /* Connected ? */
311       {
312 	sh6getlist (uintpt[ki], uintpt[kj], &index1, &index2, &kstat);
313 	if (kstat < 0)
314 	  goto error;
315 	if (!kstat)
316 	  {
317 	    /* They are connected, sort on edge par val*/
318 	    if (uintpt[ki]->epar[ind] <
319 		uintpt[kj]->epar[ind])
320 	      {
321 		klow = ki;
322 		khigh = kj;
323 	      }
324 	    else
325 	      {
326 		klow = kj;
327 		khigh = ki;
328 	      }
329 
330 	    /* Fetch geometry of surface */
331 	    sh6getgeom ((po1->iobj == SISLSURFACE) ? po1 : po2,
332 			(po1->iobj == SISLSURFACE) ? 1 : 2,
333 			uintpt[klow], &sder, &sdum, aepsge, &kstat);
334 
335 	    if (kstat < 0)
336 	      goto error;
337 
338 	    if (ipar == 0)
339 	      tdum = tsign * sder[1];
340 	    else
341 	      tdum = -tsign * sder[2];
342 
343 	    if (tdum > DZERO)
344 	      sh6setdir (uintpt[klow], uintpt[khigh], &kstat);
345 	    else if (tdum < 0)
346 	      sh6setdir (uintpt[khigh], uintpt[klow], &kstat);
347 
348 	      sh6setcnsdir (uintpt[klow], uintpt[khigh], ipar, &kstat);
349 
350 	    /* Mark edge intersection curve */
351 	    if (sh6ismain (uintpt[ki]) && sh6ismain (uintpt[kj]) &&
352 		po1->o1 == po1 &&
353 		po2->o1 == po2)
354 	      {
355 
356 		ps = po1->s1;
357 		edge_f = &(uintpt[ki]->edge_1);
358 		edge_l = &(uintpt[kj]->edge_1);
359 
360 		if (po2->iobj == SISLSURFACE)
361 		  {
362 		    edge_f = &(uintpt[ki]->edge_2);
363 		    edge_l = &(uintpt[kj]->edge_2);
364 		    ps = po2->s1;
365 		  }
366 
367 		kant = 0;
368 		if (ipar == 0)
369 		  {
370 		    if (DEQUAL (uintpt[ki]->epar[ipar],
371 				ps->et1[ps->ik1 - 1]))
372 		      kant = -1;
373 		    else if (DEQUAL (uintpt[ki]->epar[ipar],
374 				     ps->et1[ps->in1]))
375 		      kant = 1;
376 		  }
377 		else
378 		  {
379 		    if (DEQUAL (uintpt[ki]->epar[ipar], ps->et2[ps->ik2 - 1]))
380 		      kant = 1;
381 		    else if (DEQUAL (uintpt[ki]->epar[ipar],
382 				     ps->et2[ps->in2]))
383 		      kant = -1;
384 		  }
385 
386 		if (kant)
387 		  {
388 		    if ((double) kant * tdum > DZERO)
389 		      *edge_f = *edge_l = SI_RIGHT;
390 		    else if ((double) kant * tdum < DZERO)
391 		      *edge_f = *edge_l = SI_LEFT;
392 
393 		  }
394 
395 	      }
396 
397 	  }
398       }
399 
400   *jstat = 0;
401   goto out;
402 
403   /* Error in input. Incorrect dimension.  */
404 
405 err106:*jstat = -106;
406   goto out;
407 
408   /* Error lower level routine.  */
409 
410 error:*jstat = kstat;
411   goto out;
412 
413 out:
414   return;
415 }
416 
417 
418 #if defined(SISLNEEDPROTOTYPES)
419 static void
sh1782_s9sf_cu(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat ** rintdat,SISLIntpt ** uintpt,int kpoint,int ipar,int * jstat)420 sh1782_s9sf_cu (SISLObject * po1, SISLObject * po2, double aepsge,
421 		SISLIntdat ** rintdat, SISLIntpt ** uintpt, int kpoint,
422 		int ipar, int *jstat)
423 #else
424 static void
425 sh1782_s9sf_cu (po1, po2, aepsge, rintdat, uintpt, kpoint, ipar, jstat)
426      SISLObject *po1;
427      SISLObject *po2;
428      double aepsge;
429      SISLIntdat **rintdat;
430      SISLIntpt **uintpt;
431      int kpoint;
432      int ipar;
433      int *jstat;
434 #endif
435 /*
436 *********************************************************************
437 *
438 *********************************************************************
439 *
440 * PURPOSE    : Set pre-topology data in 3-dimensional surface-curve
441 *              intersection.
442 *
443 *
444 * INPUT      : po1      - Pointer to the first object in the intersection.
445 *              po2      - Pointer to the second object in the intersection.
446 *              aepsge   - Geometry resolution.
447 *              rintdat  - Intersection data structure.
448 *              uintpt   - Intersection points.
449 *              kpoint   - No of int points
450 *              ipar     - Parameter direction of missing parameter in the
451 *                         lower dimensional parameter space.
452 *
453 * OUTPUT     : uintpt   - Updated intersection point.
454 *              jstat    - status messages
455 *                                > 0   : Warning.
456 *                                = 0   : Ok.
457 *                                < 0   : Error.
458 *
459 *
460 * METHOD     :
461 *
462 * CALLS      : s6crss    -  Cross product between two vectors.
463 *              s6scpr    -  Scalar product between two vectors.
464 *              sh6getgeo -  Fetch geometry information of intersection point.
465 *
466 * REFERENCES :
467 *
468 * WRITTEN BY : ALA & VSK okt. 1992.
469 *********************************************************************
470 */
471 {
472   int kstat = 0;		/* Status variable.                        */
473   int kdim;			/* Dimension of geometry space.            */
474   double tdum;			/* Parameter used to browse on the cases.  */
475   double *val_s1;		/* Surface geometry                        */
476   double *val_s2;		/* Surface geometry                        */
477   double *snorm1;		/* Normal vector of surface.               */
478   double *snorm2;		/* Normal vector of surface.               */
479   double *sdir;	        	/* Direction vector of intersection curve. */
480   double *stang;		/* Tangent vector of edge curve.           */
481   int ind;			/* Index for curve's par val.              */
482   int ki, kj, klow, khigh;	/* Help indexes into uintpt                */
483   int index1, index2;		/* Dummy in this context                   */
484   int kcurve;			/* Which object is the curve?		   */
485   /* --------------------------------------------------------------------- */
486 
487   *jstat = 0;
488 
489   /* Test input */
490 
491   if (po1->iobj == SISLCURVE && po2->iobj == SISLSURFACE)
492   {
493      kcurve = 1;
494      kdim = po1->c1->idim;
495      if (kdim != 3)
496 	goto err104;
497      if (kdim != po2->s1->idim)
498 	goto err106;
499   }
500   else if (po1->iobj == SISLSURFACE && po2->iobj == SISLCURVE)
501   {
502      kcurve = 2;
503      kdim = po1->s1->idim;
504      if (kdim != 3)
505 	goto err104;
506      if (kdim != po2->c1->idim)
507 	goto err106;
508   }
509   else
510     goto err104;
511 
512   /* Get index for curve par val  */
513   if (kcurve == 1 && ipar == 1)
514     ind = 2;
515   else  if ((kcurve == 1 && ipar == 2) ||
516 	    (kcurve == 2 && ipar == 0))
517     ind = 1;
518   else  if (kcurve == 2 && ipar == 1)
519     ind = 0;
520   else
521      goto out;
522 
523   /* Loop for all points */
524   for (ki = 0; ki < kpoint; ki++)
525     for (kj = ki + 1; kj < kpoint; kj++)
526       /* Connected ? */
527       {
528 	sh6getlist (uintpt[ki], uintpt[kj], &index1, &index2, &kstat);
529 	if (kstat < 0)
530 	  goto error;
531 	if (!kstat)
532 	  {
533 	    /* They are connected, sort on edge par val*/
534 	    if (uintpt[ki]->epar[ind] <
535 		uintpt[kj]->epar[ind])
536 	      {
537 		klow = ki;
538 		khigh = kj;
539 	      }
540 	    else
541 	      {
542 		klow = kj;
543 		khigh = ki;
544 	      }
545 
546 	    /* Get geometry first object */
547 	    sh6getgeom (po1, 1, uintpt[klow], &val_s1, &snorm1, aepsge, &kstat);
548 	    if (kstat < 0)
549 	      goto error;
550 
551 	    /* Get geometry second object */
552 	    sh6getgeom (po2, 2, uintpt[klow], &val_s2, &snorm2, aepsge, &kstat);
553 	    if (kstat < 0)
554 	      goto error;
555 
556 	    /* Fetch tangent of edge curve.  */
557 	    if (ipar == 0)
558 	      stang = val_s1 + 6;
559 	    else if (ipar == 1 && kcurve == 2)
560 	      stang = val_s1 + 3;
561 	    else if (ipar == 1 && kcurve == 1)
562 	      stang = val_s2 + 6;
563 	    else
564 	      stang = val_s2 + 3;
565 
566 	    /* The direction of the intersection curve is set along
567 	       the direction of the curve in the intersection.  */
568 	    if (kcurve == 1)
569 	       sdir = val_s1 + 3;
570 	    else
571 	       sdir = val_s2 + 3;
572 
573 	    tdum = s6scpr (sdir, stang, kdim);
574 
575 	    if (tdum > 0)
576 	      sh6setdir (uintpt[klow], uintpt[khigh], &kstat);
577 	    else if (tdum < 0)
578 	      sh6setdir (uintpt[khigh], uintpt[klow], &kstat);
579 
580 	      sh6setcnsdir (uintpt[klow], uintpt[khigh], ipar, &kstat);
581 	  }
582       }
583   *jstat = 0;
584   goto out;
585 
586 
587   /* Error in input. Dimension not equal to 3. */
588 
589 err104:*jstat = -104;
590   goto out;
591 
592   /* Error in input. Conflicting dimensions.  */
593 
594 err106:*jstat = -106;
595   goto out;
596 
597   /* Error lower level routine.  */
598 
599 error:*jstat = kstat;
600   goto out;
601 
602 out:
603   return;
604 }
605 
606 
607 
608 #if defined(SISLNEEDPROTOTYPES)
609 static void
sh1782_s9sf_sf(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat ** rintdat,SISLIntpt ** uintpt,int kpoint,int ipar,int * jstat)610 sh1782_s9sf_sf (SISLObject * po1, SISLObject * po2, double aepsge,
611 		SISLIntdat ** rintdat, SISLIntpt ** uintpt, int kpoint,
612 		int ipar, int *jstat)
613 #else
614 static void
615 sh1782_s9sf_sf (po1, po2, aepsge, rintdat, uintpt, kpoint, ipar, jstat)
616      SISLObject *po1;
617      SISLObject *po2;
618      double aepsge;
619      SISLIntdat **rintdat;
620      SISLIntpt **uintpt;
621      int kpoint;
622      int ipar;
623      int *jstat;
624 #endif
625 /*
626 *********************************************************************
627 *
628 *********************************************************************
629 *
630 * PURPOSE    : Set pre-topology data in 3-dimensional surface-surface
631 *              intersection.
632 *
633 *
634 * INPUT      : po1      - Pointer to the first object in the intersection.
635 *              po2      - Pointer to the second object in the intersection.
636 *              aepsge   - Geometry resolution.
637 *              rintdat  - Intersection data structure.
638 *              uintpt   - Intersection points.
639 *              knpoint  - No of int points
640 *              ipar     - Parameter direction of missing parameter in the
641 *                         lower dimensional parameter space.
642 *
643 * OUTPUT     : uintpt   - Updated intersection point.
644 *              jstat    - status messages
645 *                                > 0   : Warning.
646 *                                = 0   : Ok.
647 *                                < 0   : Error.
648 *
649 *
650 * METHOD     :
651 *
652 * CALLS      : s6crss    -  Cross product between two vectors.
653 *              s6scpr    -  Scalar product between two vectors.
654 *              sh6getgeo -  Fetch geometry information of intersection point.
655 *
656 * REFERENCES :
657 *
658 * WRITTEN BY : Ulf J. Krystad, SI, Oslo, Norway. 07.91.
659 *********************************************************************
660 */
661 {
662   int kstat = 0;		/* Status variable.                        */
663   int kdim;			/* Dimension of geometry space.            */
664   double tdum;			/* Parameter used to browse on the cases.  */
665   double *val_s1;		/* Surface geometry                        */
666   double *val_s2;		/* Surface geometry                        */
667   double *snorm1;		/* Normal vector of surface.               */
668   double *snorm2;		/* Normal vector of surface.               */
669   double sdir[3];		/* Direction vector of intersection curve. */
670   double *stang;		/* Tangent vector of edge curve.           */
671   int ind;			/* Index for curve's par val.              */
672   int ki, kj, klow, khigh;	/* Help indexes into uintpt                */
673   int index1, index2;		/* Dummy in this context                   */
674   int *edge_f, *edge_l;		/* Pointers to atribute edge in intpt      */
675   int kant;			/* Edge indicator                          */
676   SISLSurf *ps;			/* The 3D Surface enhanced                 */
677   /* --------------------------------------------------------------------- */
678 
679   /* Test input */
680   kdim = po1->s1->idim;
681   if (kdim != 3)
682     goto err104;
683   if (kdim != po2->s1->idim)
684     goto err106;
685 
686   /* Get index for curve par val  */
687   if (ipar == 0)
688     ind = 1;
689   else if (ipar == 1)
690     ind = 0;
691   else if (ipar == 2)
692     ind = 3;
693   else
694     ind = 2;
695 
696   /* Loop for all points */
697   for (ki = 0; ki < kpoint; ki++)
698     for (kj = ki + 1; kj < kpoint; kj++)
699       /* Connected ? */
700       {
701 	sh6getlist (uintpt[ki], uintpt[kj], &index1, &index2, &kstat);
702 	if (kstat < 0)
703 	  goto error;
704 	if (!kstat)
705 	  {
706 	    /* They are connected, sort on edge par val*/
707 	    if (uintpt[ki]->epar[ind] <
708 		uintpt[kj]->epar[ind])
709 	      {
710 		klow = ki;
711 		khigh = kj;
712 	      }
713 	    else
714 	      {
715 		klow = kj;
716 		khigh = ki;
717 	      }
718 
719 
720 
721 
722 	    /* Get geometry first object */
723 	    sh6getgeom (po1, 1, uintpt[klow], &val_s1, &snorm1, aepsge, &kstat);
724 	    if (kstat < 0)
725 	      goto error;
726 
727 	    /* Get geometry second object */
728 	    sh6getgeom (po2, 2, uintpt[klow], &val_s2, &snorm2, aepsge, &kstat);
729 	    if (kstat < 0)
730 	      goto error;
731 
732 	    /* Get geometry lower level object */
733 	    /* Fetch tangent of edge curve.  */
734 	    if (ipar == 0)
735 	      stang = val_s1 + 6;
736 	    else if (ipar == 1)
737 	      stang = val_s1 + 3;
738 	    else if (ipar == 2)
739 	      stang = val_s2 + 6;
740 	    else
741 	      stang = val_s2 + 3;
742 
743 	    /* Compute direction of intersection curve.  */
744 	    s6crss (snorm1, snorm2, sdir);
745 
746 	    tdum = s6scpr (sdir, stang, kdim);
747 
748 	    if (tdum > 0)
749 	      sh6setdir (uintpt[klow], uintpt[khigh], &kstat);
750 	    else if (tdum < 0)
751 	      sh6setdir (uintpt[khigh], uintpt[klow], &kstat);
752 
753 	      sh6setcnsdir (uintpt[klow], uintpt[khigh], ipar, &kstat);
754 	    /* Mark edge intersection curve */
755 	    if (sh6ismain (uintpt[ki]) && sh6ismain (uintpt[kj]) &&
756 		po1->o1 == po1 &&
757 		po2->o1 == po2)
758 	      {
759 
760 		if (ipar < 2)
761 		  /* First surf */
762 		  {
763 		    edge_f = &(uintpt[ki]->edge_1);
764 		    edge_l = &(uintpt[kj]->edge_1);
765 		    ps = po1->s1;
766 		  }
767 		else
768 		  {
769 		    edge_f = &(uintpt[ki]->edge_2);
770 		    edge_l = &(uintpt[kj]->edge_2);
771 		    ps = po2->s1;
772 		  }
773 
774 		kant = 0;
775 		if (ipar == 0 || ipar == 2)
776 		  {
777 		    if (DEQUAL (uintpt[ki]->epar[ipar],
778 				ps->et1[ps->ik1 - 1]))
779 		      kant = -1;
780 		    else if (DEQUAL (uintpt[ki]->epar[ipar],
781 				     ps->et1[ps->in1]))
782 		      kant = 1;
783 		  }
784 		else
785 		  {
786 		    if (DEQUAL (uintpt[ki]->epar[ipar], ps->et2[ps->ik2 - 1]))
787 		      kant = 1;
788 		    else if (DEQUAL (uintpt[ki]->epar[ipar],
789 				     ps->et2[ps->in2]))
790 		      kant = -1;
791 		  }
792 
793 		if (kant)
794 		  {
795 		    if ((double) kant * tdum > DZERO)
796 		      *edge_f = *edge_l = SI_RIGHT;
797 		    else if ((double) kant * tdum < DZERO)
798 		      *edge_f = *edge_l = SI_LEFT;
799 
800 		  }
801 
802 
803 
804 
805 
806 	      }
807 	  }
808       }
809   *jstat = 0;
810   goto out;
811 
812 
813   /* Error in input. Dimension not equal to 3. */
814 
815 err104:*jstat = -104;
816   goto out;
817 
818   /* Error in input. Conflicting dimensions.  */
819 
820 err106:*jstat = -106;
821   goto out;
822 
823   /* Error lower level routine.  */
824 
825 error:*jstat = kstat;
826   goto out;
827 
828 out:
829   return;
830 }
831