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: s1388.c,v 1.3 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1388
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1388(SISLSurf * ps1,double * gcoons[],int * jnumb1,int * jnumb2,int * jdim,int * jstat)55 s1388(SISLSurf *ps1,double *gcoons[],int *jnumb1,int *jnumb2,int *jdim,int *jstat)
56 #else
57 void s1388(ps1,gcoons,jnumb1,jnumb2,jdim,jstat)
58 SISLSurf *ps1;
59 double *gcoons[];
60 int *jnumb1;
61 int *jnumb2;
62 int *jdim;
63 int *jstat;
64 #endif
65 /*
66 *********************************************************************
67 *
68 *********************************************************************
69 *
70 * PURPOSE : Convert a B-spline surface of order up to 4 in both
71 * direction to a mesh of Coons patches with uniform
72 * parametrization. The function assumes that the B-spline
73 * surface is C1.
74 *
75 *
76 * INPUT : ps1 - Pointer to the surface to be converted
77 *
78 *
79 * OUTPUT : gcoons - Array containing the sequence of Coons patches.
80 * The total number of patches is jnumb1*jnumb2.
81 * The patches are stored in sequence with jdim*16
82 * doubles for each patch. For each corner of the patch
83 * we store in sequence position, derivative in first
84 * direction, derivative in second direction and twist.
85 * This array is allocated inside the routine and must
86 * be released by the calling routine.
87 *
88 * jstat - status messages
89 * = 1 : Orders too high surface
90 * interpolated
91 * = 0 : ok
92 * < 0 : error
93 *
94 *
95 * METHOD : For each polynomial patch the corner position, (1,0)-
96 * (0,1)- and (1,1)-derivative is calculated.
97 *
98 * USE: int knumb1,knumb2,kdim,kstat;
99 * double *scoons;
100 * SISLSurf *qs1;
101 * .
102 * .
103 * s1388(qs1,&scoons,&knumb1,&knumb2,&kdim,&kstat);
104 * .
105 *
106 * If one of the orders of the surface is greater than four
107 * (*jstat==1) then degree reduction (order reduction) should
108 * be applied before using this routine to get a satisfactory
109 * representation of the surface by Coons patches.
110 * The degree reduction routine is s1348.
111 *
112 * REFERENCES :
113 *
114 *-
115 * CALLS : s1424, s6err
116 *
117 * WRITTEN BY : Tor Dokken, SI, Norway, 1988-11
118 *
119 *********************************************************************
120 */
121 {
122 int kstat=0; /* Local status variable. */
123 int kpos=0; /* Position of error. */
124 int kdim; /* Dimension of the space in which the surface lies. */
125 int kder=1; /* Calculate all first derivatives */
126 int klfs=0; /* Pointer into first knot vector */
127 int klft=0; /* Pointer into second knot vector */
128 int kdumlfs; /* Temporary pointer into first knot vector */
129 int kdumlft; /* Temporary pointer into second knot vector */
130 int ksize; /* Number of doubles to store a Coons patch */
131
132 int kj; /* Control variables in for loop */
133 int kn1; /* Number of vertices in first parameter direction */
134 int kn2; /* Number of vertices in first parameter direction */
135 int kk1; /* Order in first parameter direction */
136 int kk2; /* Order in first parameter direction */
137 double *st1; /* Knots in first parameter direction */
138 double *st2; /* Knots in second parameter direction */
139 double spar[2]; /* Current parameter value */
140 double sparx[2]; /* Temporary parameter value */
141 double tdiff1; /* Length of parameter interval in direction 1 */
142 double tdiff2; /* Length of parameter interval in direction 2 */
143 double *scorn1; /* Pointer to corner 1 of current Coons patch */
144 double *scorn2; /* Pointer to corner 2 of current Coons patch */
145 double *scorn3; /* Pointer to corner 3 of current Coons patch */
146 double *scorn4; /* Pointer to corner 4 of current Coons patch */
147 double tdum; /* Temporary variable */
148
149 kn1 = ps1->in1;
150 kn2 = ps1->in2;
151 kk1 = ps1->ik1;
152 kk2 = ps1->ik2;
153 kdim = ps1 ->idim;
154 st1 = ps1 -> et1;
155 st2 = ps1 -> et2;
156
157 /* Calculate number of doubles to store a Coons patch */
158
159 ksize = kdim*16;
160
161
162 /* Allocate array for storage of the coefficients */
163
164 *gcoons = newarray((kn1*kn2*16*kdim),DOUBLE);
165 if (*gcoons == SISL_NULL) goto err101;
166
167 klft = kk2 - 1;
168
169 *jnumb2 = 0;
170
171 scorn1 = *gcoons;
172
173 while (klft < kn2)
174 {
175 *jnumb1 = 0;
176 klfs = kk1 - 1;
177 while (klfs < kn1)
178 {
179
180 /* Set pointers to the corners */
181
182 scorn2 = scorn1 + 4*kdim;
183 scorn3 = scorn2 + 4*kdim;
184 scorn4 = scorn3 + 4*kdim;
185
186 spar[0] = st1[klfs];
187 spar[1] = st2[klft];
188
189 /* The parameter pair spar describes the lower left corner of the
190 current polynomial patch. By evaluating at spar we get
191 pointers into the knot vectors telling where we are. These are
192 klfs and klft t. We can calulate in the other corners by
193 just adding one to these pointers and find their parameter value.
194 */
195
196 /* Calulate first corner of patch */
197
198 s1424(ps1,kder,kder,spar,&klfs,&klft,scorn1,&kstat);
199 if (kstat<0) goto error;
200
201 /* Find length of parameter intervals */
202
203 tdiff1 = st1[klfs+1] - st1[klfs];
204 tdiff2 = st2[klft+1] - st2[klft];
205
206 /* Calculate second corner of patch */
207
208 sparx[0] = st1[klfs+1];
209 sparx[1] = spar[1];
210 kdumlfs = klfs;
211 kdumlft = klft;
212
213 s1424(ps1,kder,kder,sparx,&kdumlfs,&kdumlft,scorn2,&kstat);
214 if (kstat<0) goto error;
215
216
217 /* Calculate third corner of patch */
218
219 sparx[0] = spar[0];
220 sparx[1] = st2[klft+1];
221 kdumlfs = klfs;
222 kdumlft = klft;
223
224 s1424(ps1,kder,kder,sparx,&kdumlfs,&kdumlft,scorn3,&kstat);
225 if (kstat<0) goto error;
226
227 /* Calculate third corner of patch */
228
229 sparx[0] = st1[klfs+1];
230 sparx[1] = st2[klft+1];
231 kdumlfs = klfs;
232 kdumlft = klft;
233
234 s1424(ps1,kder,kder,sparx,&kdumlfs,&kdumlft,scorn4,&kstat);
235 if (kstat<0) goto error;
236
237 /* Scale derivatives to match uniform parametrization */
238
239
240 /* Scale derivatives in first direction */
241
242 for (kj=kdim;kj<2*kdim;kj++)
243 {
244 scorn1[kj] *= tdiff1;
245 scorn2[kj] *= tdiff1;
246 scorn3[kj] *= tdiff1;
247 scorn4[kj] *= tdiff1;
248 }
249
250 /* Scale derivatives in second direction */
251
252 for (kj=2*kdim;kj<3*kdim;kj++)
253 {
254 scorn1[kj] *= tdiff2;
255 scorn2[kj] *= tdiff2;
256 scorn3[kj] *= tdiff2;
257 scorn4[kj] *= tdiff2;
258 }
259
260 /* Scale twist */
261
262 tdum = tdiff1*tdiff2;
263
264 for (kj=3*kdim;kj<4*kdim;kj++)
265 {
266 scorn1[kj] *= tdum;
267 scorn2[kj] *= tdum;
268 scorn3[kj] *= tdum;
269 scorn4[kj] *= tdum;
270 }
271
272 /* Update position of first point of Coons patch */
273
274 scorn1 += ksize;
275
276 klfs += 1;
277 *jnumb1 +=1;
278 }
279
280 klft += 1;
281 *jnumb2 +=1;
282 }
283
284 /* The array is probably too big for the Coons patches, decrease the
285 array size */
286
287
288 /* Allocate array for storage of the coefficients */
289
290 *gcoons = increasearray(*gcoons,((*jnumb1)*(*jnumb2)*16*kdim),DOUBLE);
291 if (*gcoons == SISL_NULL) goto err101;
292
293
294 /* Test if order to high */
295
296 *jdim = kdim;
297
298 if (kk1>4 || kk2 > 4) goto war01;
299
300 *jstat = 0;
301
302 goto out;
303
304 /* Orders too high */
305
306 war01:*jstat=1;
307 goto out;
308
309
310
311
312 /* Error in scratch allocation */
313
314 err101: *jstat = -101;
315 s6err("s1388",*jstat,kpos);
316 goto out;
317
318 /* Error in lower level routine. */
319
320 error : *jstat = kstat;
321 s6err("s1388",*jstat,kpos);
322 goto freeout;
323
324 /* Some error has occured free allocated space */
325
326 freeout:
327 if (*gcoons != SISL_NULL) freearray(*gcoons);
328 goto out;
329
330 out:
331 return;
332 }
333
334