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: sh1856.c,v 1.2 2001-03-19 15:59:06 afr Exp $
45  *
46  */
47 
48 
49 #define SH1856
50 
51 #include "sislP.h"
52 
53 #if defined(SISLNEEDPROTOTYPES)
sh1856(SISLSurf * ps1,double epoint[],double edir[],int idim,double aepsco,double aepsge,int trackflag,int * jtrack,SISLTrack *** wtrack,int * jpt,double ** gpar,int ** pretop,int * jcrv,SISLIntcurve *** wcurve,int * jstat)54 void sh1856(SISLSurf *ps1,double epoint[],double edir[],int idim,
55 	    double aepsco,double aepsge,
56 	    int trackflag, int *jtrack, SISLTrack *** wtrack,
57 	    int *jpt,double **gpar,int **pretop,int *jcrv,
58 	    SISLIntcurve ***wcurve,int *jstat)
59 #else
60 void sh1856(ps1,epoint,edir,idim,aepsco,aepsge,
61 	    trackflag,jtrack,wtrack,jpt,gpar,pretop,jcrv,wcurve,jstat)
62      SISLSurf     *ps1;
63      double   epoint[];
64      double   edir[];
65      int      idim;
66      double   aepsco;
67      double   aepsge;
68      int       trackflag;
69      int       *jtrack;
70      SISLTrack ***wtrack;
71      int      *jpt;
72      double   **gpar;
73      int      **pretop;
74      int      *jcrv;
75      SISLIntcurve ***wcurve;
76      int      *jstat;
77 #endif
78 /*
79 *********************************************************************
80 *
81 *********************************************************************
82 *
83 * PURPOSE    : Find all intersections between a tensor-product surface
84 *              and an infinite straight line.
85 *
86 *
87 *
88 * INPUT      : ps1    - Pointer to surface.
89 *              epoint - SISLPoint on the line.
90 *              edir   - Direction vector of the line.
91 *              idim   - Dimension of the space in which the line lies.
92 *              aepsco - Computational resolution.
93 *              aepsge - Geometry resolution.
94 *              trackflag - If true, create tracks.
95 *
96 *
97 *
98 * OUTPUT     : jtrack - Number of tracks created
99 *              wtrack - Array of pointers to tracks
100 *              jpt    - Number of single intersection points.
101 *              gpar   - Array containing the parameter values of the
102 *                       single intersection points in the parameter
103 *                       plane of the surface. The points lie continuous.
104 *                       Intersection curves are stored in wcurve.
105 *              pretop - Topology info. for single intersection points.
106 *              *jcrv  - Number of intersection curves.
107 *              wcurve  - Array containing descriptions of the intersection
108 *                       curves. The curves are only described by points
109 *                       in the parameter plane. The curve-pointers points
110 *                       to nothing. (See description of Intcurve
111 *                       in intcurve.dcl).
112 *              jstat  - status messages
113 *                                         > 0      : warning
114 *                                         = 0      : ok
115 *                                         < 0      : error
116 *
117 *
118 * METHOD     : The line is described as the intersection between two
119 *              planes. The vertices of the surface are put into the equation
120 *              of this planes achieving a surface in the two-dimentional
121 *              space. Then the zeroes of this surface is found.
122 *
123 *
124 * REFERENCES :
125 *
126 *-
127 * CALLS      : sh1761 - Perform point object-intersection.
128 *              s1328 - Equation of surface into equations of two planes.
129 *              s1329 - Equation of surface into equation of plane.
130 *              make_sf_kreg   - Ensure k-regularity of surface.
131 *              hp_s1880 - Put intersections on output format.
132 *              s6twonorm - Make two vectors of length one that are normal
133 *                          to a 3d input vector.
134 *              newPoint    - Create new point.
135 *              newObject - Create new object.
136 *              freeObject - Free space occupied by an object.
137 *              freeIntdat  - Free space occupied by an intersection data.
138 *
139 * WRITTEN BY : Vibeke Skytt, SI, 88-06.
140 * REWRITTEN BY : Bjoern Olav Hoset, SI, 89-06.
141 *
142 *********************************************************************
143 */
144 {
145   double *nullp = SISL_NULL;
146   int kstat = 0;           /* Local status varible.                        */
147   int kpos = 0;            /* Position of error.                           */
148   int kdim;                /* Dimension of space in which the point in the
149 			      intersect point and surface problem lies.    */
150   double *spar = SISL_NULL;     /* Dummy array containing parameter values of
151 			      second object of single intersection points. */
152   double spoint[2];        /* SISLPoint to intersect with object.              */
153   double *snorm1 = SISL_NULL;   /* Normal to direction vector of line.          */
154   double *snorm2 = SISL_NULL;   /* Normal to direction vector of line and snorm1.*/
155   SISLSurf *qs = SISL_NULL;         /* Pointer to surface in
156 			      surface/point intersection.*/
157   SISLPoint *qp = SISL_NULL;        /* Pointer to point in
158 			      surface/point intersection.  */
159   SISLObject *qo1 = SISL_NULL;      /* Pointer to surface in
160 			      object/point intersection. */
161   SISLObject *qo2 = SISL_NULL;      /* Pointer to point in
162 			      object/point intersection    */
163   SISLIntdat *qintdat = SISL_NULL;  /* Intersection result */
164   int      ksurf=0;         /* Dummy number of Intsurfs. */
165   SISLIntsurf **wsurf=SISL_NULL;    /* Dummy array of Intsurfs. */
166   int      kdeg=2000;       /* input to int_join_per. */
167   SISLObject *track_obj=SISL_NULL;
168   SISLSurf *qkreg=SISL_NULL; /* Input surface ensured k-regularity. */
169 
170   /* -------------------------------------------------------- */
171 
172   if (ps1->cuopen_1 == SISL_SURF_PERIODIC ||
173       ps1->cuopen_2 == SISL_SURF_PERIODIC)
174   {
175      /* Cyclic surface. */
176 
177      make_sf_kreg(ps1,&qkreg,&kstat);
178      if (kstat < 0) goto error;
179    }
180   else
181     qkreg = ps1;
182 
183   /*
184   * Create new object and connect surface to object.
185   * ------------------------------------------------
186   */
187 
188   if (!(track_obj = newObject (SISLSURFACE)))
189     goto err101;
190   track_obj->s1 = ps1;
191 
192   /*
193    * Check dimension.
194    * ----------------
195    */
196 
197   *jpt  = 0;
198   *jcrv = 0;
199   *jtrack = 0;
200 
201   if (idim != 2 && idim != 3) goto err105;
202   if (idim != qkreg -> idim) goto err106;
203 
204   /*
205    * Allocate space for normal vectors.
206    * ----------------------------------
207    */
208 
209   snorm1 = newarray(idim,double);
210   snorm2 = newarray(idim,double);
211   if (snorm1 == SISL_NULL || snorm2 == SISL_NULL) goto err101;
212 
213   if (idim == 3)
214     {
215 
216       /*
217        * Find two planes that intersect in the given line.
218        * -------------------------------------------------
219        */
220 
221       s6twonorm(edir,snorm1,snorm2,&kstat);
222       if (kstat < 0) goto error;
223 
224       /*
225        * Put the surface into the plane equations.
226        * -----------------------------------------
227        */
228 
229       s1328(qkreg,epoint,snorm1,snorm2,idim,&qs,&kstat);
230       if (kstat < 0) goto error;
231 
232       /*
233        * Create new object and connect point to object.
234        * ----------------------------------------------
235        */
236 
237       kdim      = 2;
238       spoint[0] = spoint[1] = DZERO;
239       if (!(qo2  = newObject(SISLPOINT))) goto err101;
240       if (!(qp   = newPoint(spoint,kdim,1))) goto err101;
241       qo2 -> p1 = qp;
242     }
243   else if (idim == 2)
244     {
245 
246       /*
247        * Find normal vector of line.
248        * ---------------------------
249        */
250 
251       snorm1[0] = edir[1];
252       snorm1[1] = (-1)*edir[0];
253 
254       /*
255        * Put surface into line-equation.
256        * -------------------------------
257        */
258 
259       s1329(qkreg,epoint,snorm1,idim,&qs,&kstat);
260       if (kstat < 0) goto error;
261 
262       /*
263        * Create new object and connect point to object.
264        * ----------------------------------------------
265        */
266 
267       kdim      = 1;
268       spoint[0] = DZERO;
269       if (!(qo2  = newObject(SISLPOINT))) goto err101;
270       if (!(qp   = newPoint(spoint,kdim,1))) goto err101;
271       qo2 -> p1 = qp;
272     }
273 
274   /*
275    * Create new object and connect surface to object.
276    * ------------------------------------------------
277    */
278 
279   if(!(qo1 = newObject(SISLSURFACE))) goto err101;
280   qo1 -> s1 = qs;
281   qo1 -> o1 = qo1;
282 
283   /*
284    * Find intersections.
285    * -------------------
286    */
287 
288   sh1761(qo1,qo2,aepsge,&qintdat,&kstat);
289   if (kstat < 0) goto error;
290 
291   /* Represent degenerated intersection curves as one point.  */
292 
293   sh6degen(track_obj,track_obj,&qintdat,aepsge,&kstat);
294   if (kstat < 0) goto error;
295 
296   /* Join periodic curves */
297   int_join_per( &qintdat,track_obj,track_obj,nullp,kdeg,aepsge,&kstat);
298   if (kstat < 0)
299     goto error;
300 
301   /* Create tracks */
302   if (trackflag && qintdat)
303     {
304       make_tracks (qo1, qo2, 0, nullp,
305 		   qintdat->ilist, qintdat->vlist,
306 		   jtrack, wtrack, aepsge, &kstat);
307       if (kstat < 0)
308 	goto error;
309     }
310 
311   /*
312    * Express intersections on output format.
313    * ---------------------------------------
314    */
315 
316   if (qintdat)/* Only if there were intersections found */
317     {
318       hp_s1880(track_obj, track_obj, kdeg,
319 	       2,0,qintdat,jpt,gpar,&spar,pretop,jcrv,wcurve,&ksurf,&wsurf,&kstat);
320       if (kstat < 0) goto error;
321     }
322 
323   /*
324    * Intersections found.
325    * --------------------
326    */
327 
328   *jstat = 0;
329   goto out;
330 
331   /* Error in space allocation.  */
332 
333  err101: *jstat = -101;
334         s6err("sh1856",*jstat,kpos);
335         goto out;
336 
337   /* Error in input. Dimension different from two or three.  */
338 
339  err105: *jstat = -105;
340         s6err("sh1856",*jstat,kpos);
341         goto out;
342 
343   /* Dimensions conflicting.  */
344 
345  err106: *jstat = -106;
346         s6err("sh1856",*jstat,kpos);
347         goto out;
348 
349   /* Error in lower level routine.  */
350 
351   error : *jstat = kstat;
352         s6err("sh1856",*jstat,kpos);
353         goto out;
354 
355  out:
356 
357   /* Free allocated space.  */
358 
359   if (snorm1)  freearray(snorm1);
360   if (snorm2)  freearray(snorm2);
361   if (spar)    freearray(spar);
362   if (qo1)     freeObject(qo1);
363   if (qo2)     freeObject(qo2);
364   if (qintdat) freeIntdat(qintdat);
365   if (track_obj)
366     {
367        track_obj->s1 = SISL_NULL;
368        freeObject(track_obj);
369     }
370 
371   /* Free local surface.  */
372     if (qkreg != SISL_NULL && qkreg != ps1) freeSurf(qkreg);
373 
374 return;
375 }
376 
377