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: s1378.c,v 1.3 2001-03-19 15:58:48 afr Exp $
45 *
46 */
47
48
49 #define S1378
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1378(SISLSurf * psurf,double econic[],int ideg,int idim,SISLSurf ** rsurf,int * jstat)55 s1378 (SISLSurf * psurf, double econic[], int ideg, int idim,
56 SISLSurf ** rsurf, int *jstat)
57 #else
58 void
59 s1378 (psurf, econic, ideg, idim, rsurf, jstat)
60 SISLSurf *psurf;
61 double econic[];
62 int ideg;
63 int idim;
64 SISLSurf **rsurf;
65 int *jstat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 * PURPOSE : To put a surface description into the descripiton of
71 * a torus surface described by the input array econic.
72 *
73 * INPUT : psurf - Pointer to input surface
74 * econic - Description of torus
75 * ideg - Type of conic: torus: ideg=1001
76 * idim - Dimension of object space
77 *
78 * OUTPUT : rsurf - The resulting surface
79 * jstat - status messages
80 * > 0 : warning
81 * = 0 : ok
82 * < 0 : error
83 *
84 * METHOD : We first make the appropriate knot vector, then we calulate
85 * parametervalues for the interpolation, then the appropriate
86 * values of the surface put into the conic equation are found,
87 * and at last the surface is interpolated.
88 *
89 * REFERENCES :
90 *
91 * CALLS : s1376,s1890,s1424,s6scpr,s1891,s6err.
92 *
93 * WRITTEN BY : Tor Dokken, SI, 1988-11
94 * REVISED BY : Mike Floater, SI, 1991-01 for a rational surface.
95 *
96 *********************************************************************
97 */
98 {
99 int ikind; /* type of surface psurf is */
100 int kn1; /* Number of vertices of psurf in first par.dir */
101 int kk1; /* Order in psurf in first par.dir */
102 int kn2; /* Number of vertices of psurf in second par.dir */
103 int kk2; /* Order in psurf in second par.dir */
104 int kjkk1; /* Order of interpolated basis in first par.dir */
105 int kjkn1; /* Number of vertices in interpolated basis first.dr*/
106 int kjkk2; /* Order of interpolated basis in first par SISLdir */
107 int kjkn2; /* Number of vertices in interpolated basis secnd.dr*/
108 int kdim; /* Number of dimesions in psurf */
109 int kstat; /* Local status variable */
110 int kpos = 0; /* Position indicator for errors */
111 int kzero = 0; /* Value 0 needed in call s1891 */
112 int kone = 1; /* Value 1 needed in call s1891 */
113 int cuopen; /* Open/Closed flag */
114 int ki, kj, kl; /* Loop control variable */
115 int kp; /* Index of points put into conic equation */
116 int klfs = 0; /* Pointer into knot vector */
117 int klft = 0; /* Pointer into knot vector */
118 double *st1 = SISL_NULL; /* First knot vector is psurf */
119 double *st2 = SISL_NULL; /* Second knot vector is psurf */
120 double *scentr = econic; /* Center of torus */
121 double *saxis = econic + 3; /* Axis of torus */
122 double tbigr = *(econic + 6); /* Big radius of torus */
123 double tsmalr = *(econic + 7);/* Small radius of torus */
124 double tbigr2 = tbigr * tbigr;/* Square of big radius */
125 double tdiffr2 = tbigr2 - tsmalr * tsmalr; /* Difference of square of radia*/
126 double *sval1 = SISL_NULL; /* Array of values of surface put into torus eq. */
127 double *sval2 = SISL_NULL;
128 double *sval3 = SISL_NULL;
129 double *sgt1 = SISL_NULL; /* Knot vector in first parameter direction of
130 surface put into torus equation */
131 double *sgt2 = SISL_NULL; /* Knot vector in second parameter direction of
132 surface put into torus equation */
133 double sy[3]; /* Difference between point and torus center */
134 double tzn; /* Projection of sy onto torus axis */
135 double tyy; /* Square of length of sy */
136 double tzz; /* Square of length of sz */
137 double ty; /* Component of sy */
138 double tz; /* Component of sz */
139 double sder[4]; /* SISLPoint on the surface */
140 double spar[2]; /* Current parameter pair */
141 double ww; /* the weight of sder squared if psurf is rational */
142 double *par1 = SISL_NULL; /* Parameter vaues in direction 1. */
143 double *par2 = SISL_NULL; /* Parameter vaues in direction 2. */
144 int *der1 = SISL_NULL; /* Derivative indicators in direction 1. */
145 int *der2 = SISL_NULL; /* Derivative indicators in direction 2. */
146 SISLSurf *tempsurf = SISL_NULL; /* only used for rational surfaces */
147
148 *jstat = 0;
149
150
151 /* Test if torus. */
152
153 if (ideg != 1001)
154 goto err180;
155
156 if (idim != psurf->idim)
157 goto err104;
158
159 /* Make local pointers. */
160
161 kn1 = psurf->in1;
162 kk1 = psurf->ik1;
163 kn2 = psurf->in2;
164 kk2 = psurf->ik2;
165 kdim = psurf->idim;
166 st1 = psurf->et1;
167 st2 = psurf->et2;
168 ikind = psurf->ikind;
169
170 if (ikind == 2 || ikind == 4)
171 {
172 tempsurf = newSurf (kn1, kn2, kk1, kk2, st1, st2,
173 psurf->rcoef, ikind - 1, kdim + 1, 0);
174 if (tempsurf == SISL_NULL)
175 goto err171;
176 tempsurf->cuopen_1 = psurf->cuopen_1;
177 tempsurf->cuopen_2 = psurf->cuopen_2;
178 }
179 else
180 {
181 tempsurf = psurf;
182 }
183
184 /* Test input. */
185
186 if (kdim != 3)
187 goto err104;
188
189
190 /* Make description of knot array for interpolation in first parameter
191 direction. */
192
193 s1376 (st1, kn1, kk1, &sgt1, &kjkn1, &kjkk1, &kstat);
194 if (kstat < 0)
195 goto error;
196
197
198 /* Make parameter values and derivative indicators. */
199
200 s1890 (sgt1, kjkk1, kjkn1, &par1, &der1, &kstat);
201 if (kstat < 0)
202 goto error;
203
204
205 /* Make description of knot array for interpolation in second parameter
206 direction. */
207
208 s1376 (st2, kn2, kk2, &sgt2, &kjkn2, &kjkk2, &kstat);
209 if (kstat < 0)
210 goto error;
211
212
213 /* Make parameter values and derivative indicators. */
214
215 s1890 (sgt2, kjkk2, kjkn2, &par2, &der2, &kstat);
216 if (kstat < 0)
217 goto error;
218
219
220 /* Allocate array for values of surface put into torus equation. */
221
222 sval1 = newarray (kjkn1 * kjkn2, DOUBLE);
223 if (sval1 == SISL_NULL)
224 goto err101;
225
226
227 /* Calculate values to be interpolated. */
228
229 /* Index of point to be stored. */
230
231 kp = 0;
232
233 for (kj = 0; kj < kjkn2; kj++)
234 {
235
236 spar[1] = par2[kj];
237
238 for (ki = 0; ki < kjkn1; ki++)
239 {
240 /* Calculate values on 3-D surface */
241
242 spar[0] = par1[ki];
243
244 s1424 (tempsurf, 0, 0, spar, &klfs, &klft, sder, &kstat);
245 if (kstat < 0)
246 goto error;
247
248 /*
249 * The calculation of a point on the torus surface
250 * can be done in the following way.
251 *
252 * y = p - scentr
253 * z = y - (y saxis) saxis
254 *
255 * The equation of the torus can be written
256 *
257 * 2 2
258 * (y - R z/sqrt(z z) ) - r = 0
259 *
260 *
261 * or by elliminating the square root:
262 *
263 * f =
264 *
265 * 2 2 2 2 2 2 2
266 * (yy) + 2 (yy)(R -r ) - 4R zz + (R -r ) = 0
267 *
268 * or in 4-D homogeneous coordinates:
269 *
270 * 4
271 * f =
272 *
273 * 2 2 2 2 2 2 4 2 2 2
274 * (yy) + 2 w (yy)(R -r ) - 4w R zz + w (R -r ) = 0
275 *
276 * where Y = T - w*scentr, p+T/w
277 *
278 * We thus need to calculate yy and zz:
279 */
280
281 if (ikind == 2 || ikind == 4)
282 {
283 for (kl = 0; kl < 3; kl++)
284 sy[kl] = sder[kl] - sder[3] * scentr[kl];
285 ww = sder[3] * sder[3];
286 }
287 else
288 {
289 for (kl = 0; kl < 3; kl++)
290 sy[kl] = sder[kl] - scentr[kl];
291 ww = (double) 1.0;
292 }
293
294 tzn = s6scpr (sy, saxis, 3);
295
296 tyy = (double) 0.0;
297 tzz = (double) 0.0;
298
299 /* Make z and necessary derivatives of z */
300
301 for (kl = 0; kl < 3; kl++)
302 {
303 ty = sy[kl];
304 tz = ty - tzn * saxis[kl];
305 tyy += ty * ty;
306 tzz += tz * tz;
307 }
308
309 /* 2 2 2
310 Now tyy = yy and tzz = zz, tbigr2 = R ,tdiffr2 = R - r */
311
312 sval1[kp++] = tyy * tyy + ((double) 2.0 * ww * tyy + ww * ww * tdiffr2) * tdiffr2
313 - (double) 4.0 *ww * tbigr2 * tzz;
314 }
315 }
316
317 cuopen = TRUE;
318
319 /* Interpolate in second parameter direction, the first parameter direction
320 is treated as a point of dimension kjkn1 */
321
322 s1891 (par2, sval1, kjkn1, kjkn2, kone, der2, cuopen, sgt2, &sval2,
323 &kjkn2, kjkk2, kzero, kzero, &kstat);
324 if (kstat < 0)
325 goto error;
326
327
328 /* Interpolate in first parameter direction, perform kjkn2 interpolations
329 of one dimensional data */
330
331 s1891 (par1, sval2, kone, kjkn1, kjkn2, der1, cuopen, sgt1, &sval3,
332 &kjkn1, kjkk1, kzero, kzero, &kstat);
333 if (kstat < 0)
334 goto error;
335
336 *rsurf = SISL_NULL;
337 *rsurf = newSurf (kjkn1, kjkn2, kjkk1, kjkk2, sgt1, sgt2, sval3, 1, 1, 1);
338 if (*rsurf == SISL_NULL)
339 goto err171;
340 (*rsurf)->cuopen_1 = psurf->cuopen_1;
341 (*rsurf)->cuopen_2 = psurf->cuopen_2;
342
343 /* Ok ! */
344
345 goto out;
346
347
348 /* Error in lower level function */
349
350 error:
351 *jstat = kstat;
352 s6err ("s1378", *jstat, kpos);
353 goto out;
354
355 /* Error in space allocation */
356
357 err101:
358 *jstat = -101;
359 s6err ("s1378", *jstat, kpos);
360 goto out;
361
362 /* Dimension not equal to 3 or confliciting dim */
363
364 err104:
365 *jstat = -104;
366 s6err ("s1378", *jstat, kpos);
367 goto out;
368
369 /* Could not create surface. */
370
371 err171:
372 *jstat = -171;
373 s6err ("s1378", *jstat, kpos);
374 goto out;
375
376 /* Error in torus description */
377
378 err180:
379 *jstat = -180;
380 s6err ("s1378", *jstat, kpos);
381 goto out;
382
383 out:
384
385 /* Release allocated arrays */
386
387 if (sgt1 != SISL_NULL)
388 freearray (sgt1);
389 if (sgt2 != SISL_NULL)
390 freearray (sgt2);
391 if (sval1 != SISL_NULL)
392 freearray (sval1);
393 if (sval2 != SISL_NULL)
394 freearray (sval2);
395 if (sval3 != SISL_NULL)
396 freearray (sval3);
397 if (par1 != SISL_NULL)
398 freearray(par1);
399 if (par2 != SISL_NULL)
400 freearray(par2);
401 if (der1 != SISL_NULL)
402 freearray(der1);
403 if (der2 != SISL_NULL)
404 freearray(der2);
405 if ((ikind == 2 || ikind == 4) && (tempsurf != SISL_NULL))
406 freeSurf (tempsurf);
407
408
409 return;
410 }
411