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: s1711.c,v 1.4 2001-03-19 15:58:52 afr Exp $
45 *
46 */
47
48
49 #define S1711
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1711(SISLSurf * ps,int ipar,double apar,SISLSurf ** rsnew1,SISLSurf ** rsnew2,int * jstat)55 s1711(SISLSurf *ps,int ipar,double apar,SISLSurf **rsnew1,SISLSurf **rsnew2,int *jstat)
56 #else
57 void s1711(ps,ipar,apar,rsnew1,rsnew2,jstat)
58 SISLSurf *ps;
59 int ipar;
60 double apar;
61 SISLSurf **rsnew1;
62 SISLSurf **rsnew2;
63 int *jstat;
64 #endif
65 /*
66 ********************************************************************
67 *
68 *********************************************************************
69 *
70 * PURPOSE : Subdivide a B-spline surface at a given parameter-value
71 * and a given direction.
72 *
73 *
74 *
75 * INPUT : ps - Surface to subdivide.
76 * ipar - 1 means first direction, 2 means second direction.
77 * apar - Parameter-value at which to subdivide.
78 *
79 *
80 *
81 * OUTPUT : rsnew1 - First part of the subdivided surface.
82 * rsnew2 - Second part of the subdivided surface.
83 * jstat - status messages
84 * > 0 : warning
85 * = 0 : ok
86 * < 0 : error
87 *
88 *
89 * METHOD : Inserting knots at the subdividing-point until
90 * we have a ktuple-knot. Then we may separate the
91 * surface into two parts.
92 * The problem is to treat a three dimensional array that
93 * is locatet in a one dimensional array (coeffisients).
94 * To solve this problem we use constants to marsch in different
95 * direction in the array, and to mark end of lines or columns.
96 *
97 *
98 * REFERENCES :
99 *
100 *-
101 * CALLS : newSurf - Allocate space for a new curve-object.
102 * freeSurf - Free space occupied by given curve-object.
103 * S1700.C - Making the knot-inserten-transformation matrix.
104 *
105 * WRITTEN BY : Arne Laksaa, SI, 88-06.
106 * MODIFIED BY : Mike Floater, SI, 90-12. Subdivide rational surfaces.
107 * MODIFIED BY : Arne Laksaa, SI, 92-09. Move apar to closest knot
108 * if it is close to a knot. Using D(N)EQUAL().
109 * MODIFIED BY : Paal Fugelli, SINTEF, Oslo 15/07-1994. Changed from icopy==0
110 * to icopy==2 in creation of output surface - memory leak.
111 *
112 **********************************************************************/
113 {
114 int kstat; /* Local status variable. */
115 int kpos=0; /* Position of error. */
116 int kmy; /* An index to the knot-vector. */
117 int kv,kv1; /* Number of knots we have to insert. */
118 int kpl,kfi,kla; /* To posisjon elements in trans.-matrix.*/
119 int kk,kksec; /* Order of the input surface. */
120 int kn,knsec; /* Number of the vertices in input curves.*/
121 int kdim=ps->idim; /* Dimensjon of the space in whice surface
122 lies. */
123 int kind=ps->ikind; /* Type of surface ps is. */
124 int kn1,kn2; /* Number of vertices in the new surfaces.*/
125 int knum; /* Number of knots less and equal than
126 the intersection point. */
127 int ki,ki1,ki2; /* Control variable in loop. */
128 int kj,kj1,kj2; /* Control variable in loop. */
129 int k1m,k2m,k3m,k4m; /* Variables to mark directons in array.*/
130 int newkind=1; /* Type of surface subsurfaces are. */
131 double *s1,*s2,*s3,*s4;/* Pointers used in loop. */
132 double *st,*stsec; /* The old knot-vectors. */
133 double *st1=SISL_NULL; /* The first first new knot-vector. */
134 double *st1sec=SISL_NULL; /* The first second new knot-vector. */
135 double *st2=SISL_NULL; /* The second first new knot-vector. */
136 double *st2sec=SISL_NULL; /* The second second new knot-vector. */
137 double *salfa=SISL_NULL; /* A line of the trans.-matrix. */
138 double *scoef; /* Pointer to vertices. */
139 double *scoef1=SISL_NULL; /* The first new vertice. */
140 double *scoef2=SISL_NULL; /* The second new vertice. */
141 SISLSurf *q1=SISL_NULL; /* Pointer to new surface-object. */
142 SISLSurf *q2=SISL_NULL; /* Pointer to new surface-object. */
143 double salfa_local[5];/* Local help array. */
144
145 /* if ps is rational, do subdivision in homogeneous coordinates */
146 /* just need to set up correct dim and kind for the new surfaces at end of routine */
147 if(kind == 2 || kind == 4)
148 {
149 scoef = ps->rcoef;
150 kdim++;
151 newkind++;
152 }
153 else
154 {
155 scoef = ps->ecoef;
156 }
157
158 /* Check that we have a surface to subdivide. */
159
160 if (!ps) goto err150;
161
162 /* Making constants and ponters to mark direction. */
163
164 if (ipar==1)
165 {
166 /* If ipar is 1 we have to split the "three" dimentional
167 coeffisient matrix along a column. In this case k4m is
168 the distance beetween each element in the clumn.
169 For each element in the column we have to treat a part
170 of a line, to march along the line we use k1m.*/
171
172 st = ps->et1;
173 stsec = ps->et2;
174 kn = ps->in1;
175 knsec = ps->in2;
176 kk = ps->ik1;
177 kksec = ps->ik2;
178 k1m = kdim;
179 k4m = kdim*kn;
180 }
181 else
182 {
183 /* If ipar is 2 we have to split the "three" dimentional
184 coeffisient matrix along a line. In this case k4m is
185 the distance beetween each element in the line.
186 For each element in the line we have to treat a part
187 of a column, to march along the column we use k1m.*/
188
189 st = ps->et2;
190 stsec = ps->et1;
191 kn = ps->in2;
192 knsec = ps->in1;
193 kk = ps->ik2;
194 kksec = ps->ik1;
195 k1m = kdim*knsec;
196 k4m = kdim;
197 }
198
199 /* Check that the intersection point is an interior point. */
200
201 if ((apar < *st && DNEQUAL(apar, *st)) ||
202 (apar > st[kn+kk-1] && DNEQUAL(apar, st[kn+kk-1])))
203 goto err158;
204
205 /* Allocate space for the kk elements which may not be zero in eache
206 line of the basic transformation matrix.*/
207
208 if (kk > 5)
209 {
210 if ((salfa = newarray (kk, double)) == SISL_NULL) goto err101;
211 }
212 else salfa = salfa_local;
213
214 /* Find the number of the knots which is smaller or like
215 the intersection point, and how many knots we have to insert.*/
216
217 s1 = st;
218 kv = kk; /* The maximum number of knots we have to insert. */
219
220 if ((apar > s1[0] && DNEQUAL(apar, s1[0])) &&
221 (apar < s1[kn+kk-1] && DNEQUAL(apar, s1[kn+kk-1])))
222 {
223 /* Using binear search*/
224 kj1=0;
225 kj2=kk+kn-1;
226 knum = (kj1+kj2)/2;
227 while (knum != kj1)
228 {
229 if ((s1[knum] < apar ) && DNEQUAL(s1[knum], apar))
230 kj1=knum; else kj2=knum;
231 knum = (kj1+kj2)/2;
232 }
233 knum++; /* The smaller knots.*/
234
235 while (DEQUAL(s1[knum], apar))
236 {
237 apar = s1[knum];
238 knum++;
239 kv--;
240 }
241 /* The knots thats like the */
242 /* intersection point. */
243 }
244 else if (DEQUAL(apar,s1[0]))
245 {
246 apar = s1[0];
247 knum = 0;
248 while (s1[knum] == apar)
249 /* The knots thats like the intersection point. */
250 knum++;
251 }
252 else if (DEQUAL(apar,s1[kn+kk-1]))
253 {
254 apar = s1[kn+kk-1];
255 knum = kn+kk-1;
256 while (s1[knum-1] == apar)
257 /* The knots thats like the intersection point. */
258 knum--;
259 }
260 /* Find the number of vertices in the two new curves. */
261
262 kn1 = knum + kv - kk;
263 kn2 = kn + kk - knum;
264
265 /* Allocating the new arrays to the two new curves. */
266
267 if ((st1=newarray(kn1+kk,double))==SISL_NULL) goto err101;
268 if ((st1sec=newarray(knsec+kksec,double))==SISL_NULL) goto err101;
269 if ((st2=newarray(kn2+kk,double))==SISL_NULL) goto err101;
270 if ((st2sec=newarray(knsec+kksec,double))==SISL_NULL) goto err101;
271 if ((scoef1=newarray(kn1*kdim*knsec,double))==SISL_NULL) goto err101;
272 if ((scoef2=newarray(kn2*kdim*knsec,double))==SISL_NULL) goto err101;
273
274 /* Copying the knotvectors from the old curve to the new curves */
275
276 memcopy(st1,st,kn1,double);
277 memcopy(st2+kk,st+knum,kn2,double);
278 memcopy(st1sec,stsec,knsec+kksec,double);
279 memcopy(st2sec,stsec,knsec+kksec,double);
280
281 /* Updating the knotvectors by inserting the new k-touple knot */
282
283 for(s2=st1+kn1,s3=st2,s4=s3+kk;s3<s4;s2++,s3++) *s2 = *s3 = apar;
284
285 /* Copying the coefisientvectors to the new curves.*/
286
287 if (ipar == 1)
288 for (ki=0; ki<knsec; ki++)
289 {
290 memcopy(scoef1+ki*kdim*kn1,scoef+ki*kdim*kn,
291 kdim*kn1,double);
292 memcopy(scoef2+ki*kdim*kn2,scoef+kdim*(ki*kn+knum-kk),
293 kdim*kn2,double);
294 }
295 else
296 {
297 memcopy(scoef1,scoef,kdim*kn1*knsec,double);
298 memcopy(scoef2,scoef+kdim*(knum-kk)*knsec,
299 kdim*kn2*knsec,double);
300 }
301
302 /* Updating the coefisientvectors to the new surfaces.*/
303
304 /* Updating the first surface. */
305
306 /* If we imagine that the matrix is turned in such a way that we are
307 splitting it along a column, then for each element in the column
308 we have to treat a par of a line, to march along the line
309 in the first new matrix we use k1m, And we use k3m as a mark
310 at the end of the column in this new matrix.*/
311
312 if(ipar==1)
313 {
314 k2m=kdim*kn1;
315 k3m=kdim*kn1*knsec;
316 }
317 else
318 {
319 k2m=kdim;
320 k3m=kdim*knsec;
321 }
322 knum -= kk - 1;
323 for (ki=max(0,knum),kv1=max(0,-knum),s1=scoef1+ki*k1m;ki<kn1;ki++,s1+=k1m)
324 {
325 /* Initialising:
326 knum = knum-kk+1, Index of the first vertice to change.
327 ki = knum, Index of the vertices we are going to
328 change. Starting with knum, but if
329 knum is negativ we start at zero.
330 kv1 = 0, Number if new knots between index ki
331 and ki+kk. We are starting one below
332 becase we are counting up before using
333 it. If knum is negativ we are not
334 starting at zero but at -knum.
335 s1=scoef1+ki*k1m Pointer at the first vertice to
336 change. */
337
338 /* Using the Oslo-algorithm to make a transformation-vector
339 from the old vertices to one new vertice. */
340
341 kmy=ki;
342 s1700(kmy,kk,kn,++kv1,&kpl,&kfi,&kla,st,apar,salfa,&kstat);
343 if (kstat) goto err153;
344
345 /* Compute the knsec*kdim vertices with the "same index". */
346
347 for (s2=s1,s3=s2+k3m,ki2=0; s2<s3; s2+=k2m,ki2+=k4m)
348 for (kj=0,s4=s2; kj<kdim; kj++,s4++)
349 for (*s4=0,kj1=kfi,kj2=kfi+kpl; kj1<=kla;kj1++,kj2++)
350 *s4 += salfa[kj2] * scoef[k1m*kj1+ki2+kj];
351 }
352
353 /* And the second surface. */
354
355 /* If we imagine that the matrix is turned in such a way that we are
356 splitting it along a column, then for each element in the column
357 we have to treat a par of a line, to march along the line
358 in the second new matrix we use k1m, And we use k3m as a mark
359 at the end of the column in this new matrix.*/
360
361 if(ipar==1)
362 {
363 k2m=kdim*kn2;
364 k3m=kdim*kn2*knsec;
365 }
366 else
367 {
368 k2m=kdim;
369 k3m=kdim*knsec;
370 }
371
372 for (ki1=min(kn1+kv-1,kn+kv),s1=scoef2; ki<ki1; ki++,s1+=k1m)
373 {
374 /* Initialising:
375 ki1 = kn1+kv-1, the index of the vertice next to the
376 last vertice we have to change.
377 If we do not have so many vertices,
378 we have to use the index next to the
379 last vertice we have, kn+kv.
380 s1=scoef2 Pointer at the first vertice to
381 change. */
382
383
384 /* Using the Oslo-algorithm to make a transformation-vector
385 from the old vertices to one new vertice. */
386
387 s1700(kmy,kk,kn,kv1--,&kpl,&kfi,&kla,st,apar,salfa,&kstat);
388 if (kstat) goto err153;
389
390
391 /* Compute the knsec*kdim vertices with the "same index". */
392
393 for (s2=s1,s3=s2+k3m,ki2=0; s2<s3; s2+=k2m,ki2+=k4m)
394 for (kj=0,s4=s2; kj<kdim; kj++,s4++)
395 for (*s4=0,kj1=kfi,kj2=kfi+kpl; kj1<=kla;kj1++,kj2++)
396 *s4 += salfa[kj2] * scoef[k1m*kj1+ki2+kj];
397 }
398
399
400 /* Allocating new surface-objects.*/
401 /* use ps->idim rather than kdim in case ps is rational */
402
403
404 if (ipar==1)
405 {
406 if ((q1=newSurf(kn1,knsec,kk,kksec,st1,st1sec, /* PFU 15/07-94 */
407 scoef1,newkind,ps->idim,2)) == SISL_NULL) goto err101;
408 if ((q2=newSurf(kn2,knsec,kk,kksec,st2,st2sec, /* PFU 15/07-94 */
409 scoef2,newkind,ps->idim,2)) == SISL_NULL) goto err101;
410 }
411 else
412 {
413 if ((q1=newSurf(knsec,kn1,kksec,kk,st1sec,st1, /* PFU 15/07-94 */
414 scoef1,newkind,ps->idim,2)) == SISL_NULL) goto err101;
415 if ((q2=newSurf(knsec,kn2,kksec,kk,st2sec,st2, /* PFU 15/07-94 */
416 scoef2,newkind,ps->idim,2)) == SISL_NULL) goto err101;
417 }
418
419
420 /* Updating output. */
421
422 *rsnew1 = q1;
423 *rsnew2 = q2;
424 *jstat = 0;
425 goto out;
426
427
428 /* Error. Error in lower level function. */
429
430 err153: *jstat = kstat;
431 goto outfree;
432
433
434 /* Error. No surface to subdevice. */
435
436 err150: *jstat = -150;
437 s6err("s1711",*jstat,kpos);
438 goto out;
439
440
441 /* Error. The intersection-point is outside the surface. */
442
443 err158: *jstat = -158;
444 s6err("s1711",*jstat,kpos);
445 goto out;
446
447
448 /* Error. Allocation error, not enough memory. */
449
450 err101: *jstat = -101;
451 s6err("s1711",*jstat,kpos);
452 goto outfree;
453
454
455 outfree:
456 if(q1) freeSurf(q1);
457 q1 = SISL_NULL;
458 if(q2) freeSurf(q2);
459 q2 = SISL_NULL;
460
461 /* Free local used memory. */
462
463 out:
464 /* if(!q1) */
465 /* { */
466 /* if (st1) freearray(st1); */
467 /* if (st1sec) freearray(st1sec); */
468 /* if (scoef1) freearray(scoef1); */
469 /* } */
470
471 /* if(!q2) */
472 /* { */
473 /* if (st2) freearray(st2); */
474 /* if (st2sec) freearray(st2sec); */
475 /* if (scoef2) freearray(scoef2); */
476 /* } */
477
478 if (kk > 5 && salfa)
479 freearray (salfa);
480 return;
481 }
482