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: s1731.c,v 1.3 2001-03-19 15:58:52 afr Exp $
45 *
46 */
47
48
49 #define S1731
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1731(SISLSurf * ps,SISLSurf ** rsnew,int * jstat)55 s1731(SISLSurf *ps,SISLSurf **rsnew,int *jstat)
56 #else
57 void s1731(ps,rsnew,jstat)
58 SISLSurf *ps;
59 SISLSurf **rsnew;
60 int *jstat;
61 #endif
62 /*
63 *********************************************************************
64 *
65 *********************************************************************
66 *
67 * PURPOSE : To convert a B-spline surface to Bezier surfaces.
68 *
69 *
70 * INPUT : ps - Surface to convert.
71 *
72 *
73 *
74 * OUTPUT : rsnew - The new Bezier represented surface.
75 * jstat - status messages
76 * > 0 : warning
77 * = 0 : ok
78 * < 0 : error
79 *
80 *
81 * METHOD : Inserting knots until all knots
82 * have multiplisity pc->ik.
83 *
84 *
85 * REFERENCES :
86 *
87 *-
88 * CALLS : newarray - Allocate space for array of given type.
89 * new0array - Allocate space whith zero values.
90 * freearray - Free space occupied by given array.
91 * newSurf - Allocate space for a new surf-object.
92 * freeSurf - Free space occupied by given surf-object.
93 * S1701.C - Making the knot-inserten-transformation matrix.
94 * make_sf_kreg - Ensure that the input surface is k-regular.
95 *
96 * WRITTEN BY : Arne Laksaa, SI, 88-11.
97 * REVISED BY : Johannes Kaasa, SI, May 1992 (Introduced NURBS).
98 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, Oct. 1994. Moved free'ing
99 * of 'qkreg' from 'outfree:' to 'out:' to remove memory leak.
100 *
101 **********************************************************************/
102 {
103 register int ki,ki1,ki2; /* Control variable in loop. */
104 register int kj,kj1,kj2,kj3; /* Control variable in loop. */
105 int kstat; /* Local status variable. */
106 int kpos=0; /* Position of error. */
107 int kmy; /* An index to the knot-vector. */
108 int kpl,kfi,kla; /* To posisjon elements in trans.-matrix. */
109 int kk1=ps->ik1; /* Order one of the input surface. */
110 int kk2=ps->ik2; /* Order two of the input surface. */
111 int kn=ps->in1; /* Number of the vertices in input curves. */
112 int kdim=ps->idim; /* Dimensjon of the space in whice surf lies. */
113 int kn1,kn2; /* Number of vertices in the new surface. */
114 double *s1,*s2,*s3; /* Pointers used in loop. */
115 double *st1=SISL_NULL; /* The new knot-vector. */
116 double *st2=SISL_NULL; /* The new knot-vector. */
117 double *sp=SISL_NULL; /* To use in s1701.c */
118 double *salfa=SISL_NULL; /* A line of the trans.-matrix. */
119 double *scoef=SISL_NULL; /* The new vertice. */
120 double *scoefh=SISL_NULL; /* A new vertice for help. */
121 SISLSurf *q1=SISL_NULL; /* Pointer to new surf-object. */
122
123 double *rcoef; /* Potential rational vertices. */
124 int rdim; /* Potential rational dimension. */
125 SISLSurf *qkreg=SISL_NULL; /* Input surface made k-regular. */
126
127 /* Check that we have a surface to treat. */
128
129 if (!ps) goto err150;
130
131 /* Make sure that the surface is k-regular. */
132
133 if (ps->cuopen_1 == SISL_SURF_PERIODIC ||
134 ps->cuopen_2 == SISL_SURF_PERIODIC)
135 {
136 make_sf_kreg(ps,&qkreg,&kstat);
137 if (kstat < 0) goto err153;
138 }
139 else qkreg = ps;
140
141
142 /* Check if the surface is rational. */
143
144 if (qkreg->ikind == 2 || qkreg->ikind == 4)
145 {
146 rcoef = qkreg->rcoef;
147 rdim = kdim + 1;
148 }
149 else
150 {
151 rcoef = qkreg->ecoef;
152 rdim = kdim;
153 }
154
155 /* Allocate space for the kk elements which may not be zero in eache
156 line of the basic transformation matrix, and space for new knots
157 to use in s1701.c */
158
159 if ((salfa=newarray(kk1+kk2,double))==SISL_NULL) goto err101;
160 if ((sp=newarray(kk1+kk2,double))==SISL_NULL) goto err101;
161
162 /* Find the number of vertices in the first direction
163 in the new surface. */
164
165 for(ki=0,kn1=0;ki<kn+kk1;ki+=kj,kn1+=kk1)
166 for(kj=1;ki+kj<kn+kk1 && (qkreg->et1[ki] == qkreg->et1[ki+kj]);kj++);
167 kn1 -= kk1;
168
169 /* Find the number of vertices in the second direction
170 in the new surface. */
171
172 for(kn=qkreg->in2,ki=0,kn2=0;ki<kn+kk2;ki+=kj,kn2+=kk2)
173 for(kj=1;ki+kj<kn+kk2 && (qkreg->et2[ki] == qkreg->et2[ki+kj]);kj++);
174 kn2 -= kk2;
175
176 /* Allocating the new arrays to the new surface. */
177
178 if ((st1=newarray(kn1+kk1,double))==SISL_NULL) goto err101;
179 if ((st2=newarray(kn2+kk2,double))==SISL_NULL) goto err101;
180 if ((scoefh=new0array(kn1*kn*rdim,double))==SISL_NULL) goto err101;
181 if ((scoef=new0array(kn1*kn2*rdim,double))==SISL_NULL) goto err101;
182
183 /* Making the new knotvectors in the first direction */
184
185 for(kn=qkreg->in1,ki=0,ki1=0;ki<kn+kk1;ki+=kj)
186 {
187 for(kj=1;ki+kj<kn+kk1 && (qkreg->et1[ki] == qkreg->et1[ki+kj]);kj++);
188 for(kj1=0;kj1<kk1;kj1++,ki1++) st1[ki1] = qkreg->et1[ki];
189 }
190
191 /* Making the new knotvectors in the second direction. */
192
193 for(kn=qkreg->in2,ki=0,ki1=0;ki<kn+kk2;ki+=kj)
194 {
195 for(kj=1;ki+kj<kn+kk2 && (qkreg->et2[ki] == qkreg->et2[ki+kj]);kj++);
196 for(kj1=0;kj1<kk2;kj1++,ki1++) st2[ki1] = qkreg->et2[ki];
197 }
198
199 /* Updating the coefisientvector to the new surface in the first
200 direction.*/
201
202 for(s1=scoefh,ki2=kn1*kn*rdim,ki=0,kmy=0;ki<kn1;ki++)
203 {
204 /* Here we compute a new line with line number ki of
205 the knot inserten matrix. */
206
207 while(qkreg->et1[kmy+1] <= st1[ki]) kmy++;
208 s1701(ki,kmy,kk1,qkreg->in1,&kpl,&kfi,&kla,st1,qkreg->et1,sp,salfa,&kstat);
209 if (kstat) goto err153;
210
211 /* Compute the kn2*rdim vertices with the same "index". */
212
213 for (kj=0; kj<rdim; kj++,s1++)
214 for (s2=s1,s3=s2+ki2,kj3=kj;s2<s3;s2+=rdim*kn1,kj3+=rdim*qkreg->in1)
215 for (*s2=0,kj1=kfi,kj2=kfi+kpl; kj1<=kla; kj1++,kj2++)
216 *s2 += salfa[kj2] * rcoef[kj1*rdim+kj3];
217 }
218
219 /* Updating the coefisientvector to the new surface in the second
220 direction.*/
221
222 for(s1=scoef,ki2=kn1*rdim,ki=0,kmy=0;ki<kn2;ki++,s1+=kn1*rdim)
223 {
224 /* Here we compute a new line with line number ki of
225 the knot inserten matrix. */
226
227 while(qkreg->et2[kmy+1] <= st2[ki]) kmy++;
228 s1701(ki,kmy,kk2,kn,&kpl,&kfi,&kla,st2,qkreg->et2,sp,salfa,&kstat);
229 if (kstat) goto err153;
230
231 /* Compute the kn1*rdim vertices with the same "index". */
232
233 for (kj=0; kj<rdim; kj++)
234 for (s2=s1+kj,s3=s2+ki2,kj3=kj;s2<s3;s2+=rdim,kj3+=rdim)
235 for (*s2=0,kj1=kfi,kj2=kfi+kpl; kj1<=kla; kj1++,kj2++)
236 *s2 += salfa[kj2] * scoefh[kj1*kn1*rdim+kj3];
237 }
238
239
240 /* Allocating new surface-objects.*/
241
242 if ((q1=newSurf(kn1,kn2,kk1,kk2,st1,st2,scoef,qkreg->ikind,kdim,2)) == SISL_NULL)
243 goto err101;
244
245 /* Updating output. */
246
247 *rsnew = q1;
248 *jstat = 0;
249 goto out;
250
251 /* Error. Subrutine error. */
252
253 err153: *jstat = kstat;
254 goto outfree;
255
256 /* Error. No surface to treat. */
257
258 err150: *jstat = -150;
259 s6err("s1731",*jstat,kpos);
260 goto out;
261
262 /* Error. Allocation error, not enough memory. */
263
264 err101: *jstat = -101;
265 s6err("s1731",*jstat,kpos);
266 goto outfree;
267
268 outfree:
269 if(q1) freeSurf(q1);
270 else
271 {
272 if (st1) freearray(st1);
273 if (st2) freearray(st2);
274 if (scoef) freearray(scoef);
275 }
276
277 /* Free local used memory. */
278
279 out:
280 if (qkreg != SISL_NULL && qkreg != ps) freeSurf(qkreg);
281 if (salfa) freearray(salfa);
282 if (sp) freearray(sp);
283 if (scoefh) freearray(scoefh);
284 }
285