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