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