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:
45 *
46 */
47 #define S1505
48
49 #include "sislP.h"
50
51 #if defined(SISLNEEDPROTOTYPES)
52 void
s1505(const SISLSurf * ps1,int ider,int m1,int m2,double * ebder1,double * ebder2,int * ileft1,int * ileft2,double eder[],double norm[],int * jstat)53 s1505(const SISLSurf *ps1,int ider,int m1, int m2,double *ebder1,double
54 *ebder2, int *ileft1,int *ileft2,double eder[],double norm[],int *jstat)
55 #else
56 void s1505(ps1,ider,m1,m2,ebder1,ebder2,ileft1,ileft2,eder,norm,jstat)
57 SISLSurf *ps1;
58 int ider;
59 int m1;
60 int m2;
61 double *ebder1;
62 double *ebder2;
63 int *ileft1;
64 int *ileft2;
65 double eder[];
66 double norm[];
67 int *jstat;
68 #endif
69 /*
70 *********************************************************************
71 *
72 *********************************************************************
73 *
74 * PURPOSE : Evaluate the surface pointed at by ps1 over an m1 * m2
75 * grid, assuming that the B-splines have been
76 * pre-evaluated, by s1504, over that grid.
77 * The knots et1 and et2 and grid points (x[i],y[j]) are not needed.
78 * Compute ider derivatives.
79 *
80 * INPUT : ps1 - Pointer to the surface to evaluate.
81 * ider - Number of derivatives to calculate.
82 * < 0 : No derivative calculated.
83 * = 0 : Position calculated.
84 * = 1 : Position and first derivative calculated.
85 * etc.
86 * m1 - Number of grid points in first direction.
87 * m2 - Number of grid points in first direction.
88 * ileft1 - Array of indexes to the intervals in the knotvector
89 * in the first parameter direction where each subsequence
90 * of k1 non-zero B-splines are located.
91 * ileft2 - Array of indexes to the intervals in the knotvector
92 * in the second parameter direction where each subsequence
93 * of k2 non-zero B-splines are located.
94 * ebder1 - Triple array of dimension [(ider+1)*k1*m1] containing
95 * values of the k1 nonzero B-splines and their
96 * derivatives at the points x[0],...,x[m1-1]
97 * (i.e. pre-evaluated B-splines).
98 * These numbers are stored in the following order:
99 * First the (ider+1) derivatives of the first nonzero
100 * B-spline at x[0]. Then the (ider+1) derivatives of
101 * the second nonzero B-spline at x[0], etc.
102 * Later we repeat for x[1],... etc.
103 *
104 * ebder2 - Triple array of dimension [(ider+1)*k2*m2] containing
105 * values of the k2 nonzero B-splines and their
106 * derivatives at the points y[0],...,y[m2-1]
107 *
108 * OUTPUT : eder - Array where the derivatives of the surface
109 * are placed, dimension
110 * idim * ((ider+1)(ider+2) / 2) * m1 * m2.
111 * The sequence is position,
112 * first derivative in first parameter direction,
113 * first derivative in second parameter direction,
114 * (2,0) derivative, (1,1) derivative, (0,2)
115 * derivative, etc. at point (x[0],y[0]),
116 * followed by the same information at (x[1],y[0]),
117 * etc.
118 * norm - Normals of surface. Is calculated if ider >= 1.
119 * Dimension is idim*m1*m2.
120 * The normals are not normalized.
121 * jstat - status messages
122 * = 2 : Surface is degenerate
123 * at some point, normal
124 * has zero length.
125 * = 1 : Surface is close to
126 * degenerate at some point
127 * Angle between tangents,
128 * less than angular tolerance.
129 * = 0 : ok
130 * < 0 : error
131 *
132 * METHOD : The code and method is similar to that of s1421 except that
133 * the B-splines and their derivatives have already been
134 * calculated (and are given in ebder1 and ebder2) and
135 * we evaluate over a whole grid rather than at one point.
136 * The method is to find the control points and control derivatives
137 * of each isocurve in x (fixed y value). We then
138 * evaluate at each x value along the isocurve.
139 *
140 * CALLS : s6err - Error handling routine
141 * s6strider - Make derivative of rational expression
142 *
143 * WRITTEN BY : Michael Floater, SINTEF, May 1998.
144 *********************************************************************
145 */
146 {
147 int kstat=0; /* Local status variable. */
148 int kpos=0; /* The position of error. */
149 int i1,i2; /* Loop variables. */
150 int i1pos,i2pos; /* Offset indexes. */
151 int kn1,kn2; /* The number of B-splines accociated with the knot
152 vectors st1 and st2. */
153 int kk1,kk2; /* The polynomial order of the surface in the two
154 directions. */
155 int kdim; /* The dimension of the space in which the surface
156 lies. Equivalently, the number of components
157 of each B-spline coefficient. */
158 int kleft1,kleft2; /* Local versions of knot intervals. */
159 int ki,kx,kjh; /* Control variables in for loops and for stepping
160 through arrays. */
161 int kih2; /* Index for stepping through ebder2. */
162 int ky,kl,kl1,kl2; /* Control variables in for loops and for stepping
163 through arrays. */
164 double *scoef; /* The B-spline coefficients of the surface.
165 This is an array of dimension [kn2*kn1*kdim]. */
166 double tt; /* Dummy variable used for holding an array element
167 in a for loop. */
168 double *ew=SISL_NULL; /* Pointer to an array of dimension [kn1*(ider+1)*kdim]
169 which will be used to store the result of the first
170 matrix multiplication. */
171 double *sder=SISL_NULL; /* Pointer to array used for storage of points, if
172 rational has room for homogenous coordinates. */
173 double *enorm=SISL_NULL; /* Array for surface normal. */
174 int size; /* Space occupied by points and derivs at one eval. */
175 int sizeh; /* Space occupied by homogeneous points and derivs . */
176 int size1,size2; /* Useful variables. */
177 int ederpos; /* Index of position in eder. */
178 int normpos; /* Index of position in norm. */
179
180 int knumb2; /* Necessary size of ew */
181
182 int tot,temp; /* Temporary variables. */
183
184 kn1 = ps1 -> in1;
185 kn2 = ps1 -> in2;
186 kk1 = ps1 -> ik1;
187 kk2 = ps1 -> ik2;
188 kdim = ps1 -> idim;
189 /* Check the input. */
190
191 if (kdim < 1) goto err102;
192 if (kk1 < 1) goto err115;
193 if (kn1 < kk1 || kn2 < kk2) goto err116;
194 if (ider < 0) goto err178;
195
196 *jstat = 0;
197
198 if (ps1->ikind == 2 || ps1->ikind == 4)
199 {
200 scoef = ps1 -> rcoef;
201 kdim +=1;
202 }
203 else
204 {
205 scoef = ps1 -> ecoef;
206 }
207 sizeh = kdim*(ider+1)*(ider+2)/2;
208 if((sder=newarray(sizeh,DOUBLE)) == SISL_NULL) goto err101;
209 if((enorm=newarray(ps1->idim,DOUBLE)) == SISL_NULL) goto err101;
210
211 size = ps1->idim*(ider+1)*(ider+2)/2;
212 size1 = (ider+1)*kk1;
213 size2 = (ider+1)*kk2;
214
215 /* Allocate space for B-spline values and derivatives and one work array. */
216
217 knumb2 = kn1*(ider+1)*kdim;
218 if((ew=newarray(knumb2,DOUBLE)) == SISL_NULL) goto err101;
219
220 ederpos = 0;
221 normpos = 0;
222
223 /* Run through grid points in the y direction. */
224 for(i2=0, i2pos=0; i2<m2; i2++, i2pos += size2)
225 {
226 kleft2 = ileft2[i2];
227
228 /* Compute the control points (and control derivatives
229 of the v = x[i2] isocurve. */
230
231 /* Set all the elements of ew to 0. */
232 for(ki=0; ki<knumb2; ki++) ew[ki] = DZERO;
233
234 /* ki steps through the appropriate kk2 rows of B-spline coefficients
235 while kih2 steps through the B-spline value and derivatives for the
236 B-spline given by ki. */
237
238 kih2 = i2pos;
239 for (ki=kleft2-kk2+1; ki<=kleft2; ki++)
240 {
241 /* kx counts through the ider+1 derivatives to be computed.
242 kjh steps through ew once for each ki to accumulate the
243 contribution from the different B-splines.
244 kl1 points to the first component of the first B-spline coefficient
245 in row no. ki of the B-spline coefficient matrix.
246 */
247
248 kjh = 0; kl1 = kdim*ki*kn1;
249 for (kx=0; kx<=ider; kx++)
250 {
251 /* The value of the B-spline derivative is stored in tt while
252 kl2 steps through the kdim components of all the B-spline
253 coefficients that multiply nonzero B-splines along st1.
254 */
255
256 tt = ebder2[kih2++]; kl2 = kl1;
257 for (kl=0; kl<kdim*kn1; kl++,kjh++,kl2++)
258 {
259 ew[kjh] += scoef[kl2]*tt;
260 }
261 }
262 }
263
264 /* Run through grid points in the x direction
265 evaluating along the iso-curve defined by y = y[i2]. */
266 for(i1=0, i1pos=0; i1<m1; i1++, i1pos += size1)
267 {
268 kleft1 = ileft1[i1];
269
270 /* Set all the elements of sder to 0. */
271 for(ki=0; ki<sizeh; ki++) sder[ki] = DZERO;
272
273 for(ky=0; ky<=ider; ky++)
274 {
275 kl1 = kdim * (ky * kn1 + kleft1 - kk1 + 1);
276 for(kx=0; kx<=ider-ky; kx++)
277 {
278 tot = kx + ky;
279 temp = ((tot * (tot+1)) >> 1) + ky;
280 kjh = temp * kdim;
281
282 for(ki=0; ki<kk1; ki++)
283 {
284 tt = ebder1[i1pos + (ider+1) * ki + kx];
285 for(kl=0; kl<kdim; kl++)
286 {
287 sder[kjh+kl] += ew[kl1 + kdim * ki + kl] * tt;
288 }
289 }
290 }
291 }
292
293 /* If rational surface calculate the derivatives based on
294 derivatives in homogenous coordinates */
295
296 if (ps1->ikind == 2 || ps1->ikind == 4)
297 {
298 s6strider(sder,ps1->idim,ider,eder+ederpos,&kstat);
299 if (kstat<0) goto error;
300 }
301 else
302 {
303 for(ki=0; ki<sizeh; ki++) eder[ederpos+ki] = sder[ki];
304 }
305
306 /* Calculate normal if idim==3 and ider>0. */
307
308 if (ider>0 && ps1->idim ==3)
309 {
310
311 s6crss(eder+ederpos+ps1->idim,eder+ederpos+2*ps1->idim,enorm);
312
313 s6norm(enorm,3,norm+normpos,&kstat);
314 }
315
316 ederpos += size;
317 normpos += 3;
318 }
319 }
320
321 goto out;
322
323 /* Not enough memory. */
324 err101: *jstat = -101;
325 s6err("s1505",*jstat,kpos);
326 goto out;
327
328 /* kdim less than 1. */
329 err102: *jstat = -102;
330 s6err("s1505",*jstat,kpos);
331 goto out;
332
333 /* Polynomial order less than 1. */
334 err115: *jstat = -115;
335 s6err("s1505",*jstat,kpos);
336 goto out;
337
338 /* Fewer B-splines than the order. */
339 err116: *jstat = -116;
340 s6err("s1505",*jstat,kpos);
341 goto out;
342
343 /* Illegal derivative requested. */
344 err178: *jstat = -178;
345 s6err("s1221",*jstat,kpos);
346 goto out;
347
348 /* Error in lower level routine. */
349
350 error: *jstat = kstat;
351 s6err("s1505",*jstat,kpos);
352 goto out;
353
354 out:
355 /* Free memory. */
356 if(sder != SISL_NULL) freearray(sder);
357 if(enorm != SISL_NULL) freearray(enorm);
358 if (ew != SISL_NULL) freearray(ew);
359
360 return;
361 }
362