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: sh1855.c,v 1.2 2001-03-19 15:59:06 afr Exp $
45  *
46  */
47 
48 
49 #define SH1855
50 
51 #include "sislP.h"
52 
53 #if defined(SISLNEEDPROTOTYPES)
54 void
sh1855(SISLSurf * ps1,double ecentr[],double aradius,double enorm[],int idim,double aepsco,double aepsge,int trackflag,int * jtrack,SISLTrack *** wtrack,int * jpt,double ** gpar,int ** pretop,int * jcrv,SISLIntcurve *** wcurve,int * jstat)55 sh1855(SISLSurf *ps1,double ecentr[],double aradius,
56 	   double enorm[],int idim,double aepsco,double aepsge,
57 	int trackflag, int *jtrack, SISLTrack *** wtrack,
58 	   int *jpt,double **gpar,int **pretop,int *jcrv,SISLIntcurve ***wcurve,int *jstat)
59 #else
60 void sh1855(ps1,ecentr,aradius,enorm,idim,aepsco,aepsge,
61 	trackflag,jtrack,wtrack,jpt,gpar,pretop,jcrv,wcurve,jstat)
62      SISLSurf     *ps1;
63      double   ecentr[];
64      double   aradius;
65      double   enorm[];
66      int      idim;
67      double   aepsco;
68      double   aepsge;
69      int       trackflag;
70      int       *jtrack;
71      SISLTrack ***wtrack;
72      int      *jpt;
73      double   **gpar;
74      int      **pretop;
75      int      *jcrv;
76      SISLIntcurve ***wcurve;
77      int      *jstat;
78 #endif
79 /*
80 *********************************************************************
81 *
82 *********************************************************************
83 *
84 * PURPOSE    : Find all intersections between a tensor-product surface
85 *              and a circle.
86 *
87 *
88 *
89 * INPUT      : ps1      - Pointer to surface.
90 *              ecentr   - Center of the circle.
91 *              aradius  - Radius of the circle.
92 *              enorm    - Normal vector to the plane in which the circle
93 *                         lies.
94 *              idim     - Dimension of the space in which the circle lies.
95 *              aepsco   - Computational resolution.
96 *              aepsge   - Geometry resolution.
97 *              trackflag - If true, create tracks.
98 *
99 *
100 *
101 * OUTPUT     : jtrack - Number of tracks created
102 *              wtrack - Array of pointers to tracks
103 *              jpt    - Number of single intersection points.
104 *              gpar   - Array containing the parameter values of the
105 *                       single intersection points in the parameter
106 *                       plane of the surface. The points lie continuous.
107 *                       Intersection curves are stored in wcurve.
108 *              pretop - Topology info. for single intersection points.
109 *              *jcrv  - Number of intersection curves.
110 *              wcurve  - Array containing descriptions of the intersection
111 *                       curves. The curves are only described by points
112 *                       in the parameter plane. The curve-pointers points
113 *                       to nothing. (See description of Intcurve
114 *                       in intcurve.dcl).
115 *              jstat  - status messages
116 *                                         > 0      : warning
117 *                                         = 0      : ok
118 *                                         < 0      : error
119 *
120 *
121 * METHOD     : The vertices of the surface are put into the equation of
122 *              the circle achieving a surface in the two-dimentional space.
123 *              Then the zeroes of this surface is found.
124 *
125 *
126 * REFERENCES :
127 *
128 *-
129 * CALLS      : sh1761 - Perform point object-intersection.
130 *              s1320 - Put equation of surface into equation of implicit
131 *                      surface.
132 *              s1324 - Represent circel as two implicit functions.
133 *              make_sf_kreg   - Ensure k-regularity of surface.
134 *              hp_s1880 - Put intersections on output format.
135 *              newObject - Create new object.
136 *              newPoint    - Create new point.
137 *              freeObject - Free space occupied by an object.
138 *              freeIntdat  - Free space occupied by an intersection data.
139 *
140 * WRITTEN BY : Vibeke Skytt, SI, 88-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 = 2;               /* Dimension of space in which the point in the
149 				 intersect point/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 *scirc = SISL_NULL;       /* Description of a circel as two implicit
154 				 functions.                                */
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   /*
194    * Check dimension.
195    * ----------------
196    */
197 
198   *jpt  = 0;
199   *jcrv = 0;
200   *jtrack = 0;
201 
202   if (idim != qkreg -> idim) goto err106;
203 
204   /*
205    * Allocate space for matrix describing a cone.
206    * --------------------------------------------
207    */
208 
209   if ((scirc = newarray(2*(idim+1)*(idim+1),double)) == SISL_NULL) goto err101;
210 
211   /*
212    * Make a matrix of dimension (idim+1)x(idim+1) describing a
213    * cone as an implicit function.
214    * ---------------------------------------------------------
215    */
216 
217   s1324(ecentr,aradius,enorm,idim,scirc,&kstat);
218   if (kstat < 0) goto error;
219 
220   /*
221    * Put the description of the input surface into the implicit
222    * equation for the cone.
223    * ----------------------------------------------------------
224    */
225 
226   s1320(qkreg,scirc,2,0,&qs,&kstat);
227   if (kstat < 0) goto error;
228 
229   /*
230    * Create new object and connect surface to object.
231    * ------------------------------------------------
232    */
233 
234   if(!(qo1 = newObject(SISLSURFACE))) goto err101;
235   qo1 -> s1 = qs;
236   qo1 -> o1 = qo1;
237 
238   /*
239    * Create new object and connect point to object.
240    * ----------------------------------------------
241    */
242 
243   if(!(qo2 = newObject(SISLPOINT))) goto err101;
244   spoint[0] = spoint[1] = DZERO;
245   if(!(qp = newPoint(spoint,kdim,1))) goto err101;
246   qo2 -> p1 = qp;
247 
248   /*
249    * Find intersections.
250    * -------------------
251    */
252 
253   sh1761(qo1,qo2,aepsge,&qintdat,&kstat);
254   if (kstat < 0) goto error;
255 
256   /* Represent degenerated intersection curves as one point.  */
257 
258   sh6degen(track_obj,track_obj,&qintdat,aepsge,&kstat);
259   if (kstat < 0) goto error;
260 
261   /* Join periodic curves */
262   int_join_per( &qintdat,track_obj,track_obj,nullp,kdeg,aepsge,&kstat);
263   if (kstat < 0)
264     goto error;
265 
266   /* Create tracks */
267   if (trackflag && qintdat)
268     {
269       make_tracks (qo1, qo2, 0, nullp,
270 		   qintdat->ilist, qintdat->vlist,
271 		   jtrack, wtrack, aepsge, &kstat);
272       if (kstat < 0)
273 	goto error;
274     }
275 
276   /*
277    * Express intersections on output format.
278    * ---------------------------------------
279    */
280 
281   if (qintdat)/* Only if there were intersections found */
282     {
283       hp_s1880(track_obj, track_obj, kdeg,
284 	       2,0,qintdat,jpt,gpar,&spar,pretop,jcrv,wcurve,&ksurf,&wsurf,&kstat);
285       if (kstat < 0) goto error;
286     }
287 
288   /*
289    * Intersections found.
290    * --------------------
291    */
292 
293   *jstat = 0;
294   goto out;
295 
296   /* Error in space allocation.  */
297 
298  err101: *jstat = -101;
299         s6err("sh1855",*jstat,kpos);
300         goto out;
301 
302   /* Dimensions conflicting.  */
303 
304  err106: *jstat = -106;
305         s6err("sh1855",*jstat,kpos);
306         goto out;
307 
308   /* Error in lower level routine.  */
309 
310   error : *jstat = kstat;
311         s6err("sh1855",*jstat,kpos);
312         goto out;
313 
314  out:
315 
316   /* Free allocated space.  */
317 
318   if (spar)    freearray(spar);
319   if (scirc)   freearray(scirc);
320   if (qo1)     freeObject(qo1);
321   if (qo2)     freeObject(qo2);
322   if (qintdat) freeIntdat(qintdat);
323   if (track_obj)
324     {
325        track_obj->s1 = SISL_NULL;
326        freeObject(track_obj);
327     }
328 
329   /* Free local surface.  */
330     if (qkreg != SISL_NULL && qkreg != ps1) freeSurf(qkreg);
331 
332 return;
333 }
334 
335