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: s1370.c,v 1.3 2001-03-19 15:58:47 afr Exp $
45 *
46 */
47
48
49 #define S1370
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1370(SISLCurve * pcurv,double earray[],int idim,int inarr,int ratflag,SISLCurve ** rcurv,int * jstat)55 s1370 (SISLCurve * pcurv, double earray[], int idim, int inarr,
56 int ratflag, SISLCurve ** rcurv, int *jstat)
57 #else
58 void
59 s1370 (pcurv, earray, idim, inarr, ratflag, rcurv, jstat)
60 SISLCurve *pcurv;
61 double earray[];
62 int idim;
63 int inarr;
64 int ratflag;
65 SISLCurve ** rcurv;
66 int *jstat;
67 #endif
68 /*
69 *********************************************************************
70 *
71 * PURPOSE : To put a curve description into the implicit
72 * second order surface described by the input array.
73 *
74 * INPUT : pcurv - Pointer to input curve
75 * earray - The description of the input array
76 * dimension (idim+1)x(idim+1) (xinarr)
77 * idim - Put curve into implicit equation
78 * inarr - Number of parallel matrices in earray.
79 * inarr should be less or equal to 3.
80 * ratflag - If pcurv is nonrational it is ignored.
81 * Otherwise:
82 * If ratflag = 0 rcurv is the nonrational numerator
83 * If ratflag = 1 rcurv is a full rational curve
84 *
85 * OUTPUT : rcurv - The resulting curve
86 * jstat - status messages
87 * > 0 : warning
88 * = 0 : ok
89 * < 0 : error
90 *
91 * METHOD : Dependent on the type of object we make:
92 *
93 * F(S,T) = (P(s,t),1) EARRAY (P(s,t),1)
94 *
95 * by sampling enough point to use interpolation for reproduction.
96 *
97 * REFERENCES :
98 *
99 * CALLS : s1893,s6err.
100 *
101 * WRITTEN BY : Tor Dokken, SI, Oslo , Norway
102 * REVISED BY : Mike Floater, SI, Oslo 11/4/91 for rational curves
103 * REVISED BY : Mike Floater, SI, Oslo 11/9/91 -- ratflag.
104 * REVISED BY : Michael Floater, SI, June 92. The rational stuff
105 * was completely messed up after translation of
106 * the fortran part to c. But it works now.
107 * REVISED BY : Christophe Birkeland, SI, July 1992.
108 * REVISED BY : Christophe Rene Birkeland, SINTEF Oslo, May 1993.
109 * jcurve removed, other minor changes
110 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, September 1994.
111 * Didn't work for rationals - '(*rcurve)->rcoef' returned from
112 * s1893() was SISL_NULL (must copy from ecoef).
113 *
114 *********************************************************************
115 */
116 {
117 int kpos = 0;
118 int kstat = 0;
119 SISLCurve *icurve = SISL_NULL; /* Temporary SISLCurve. */
120 int kn; /* Number of vertices of pcurv */
121 int kk; /* Order in pcurv */
122 int kdim; /* Number of dimesions in pcurv */
123 int kdimp1; /* Dimension of earray should be kdim+1 */
124 double *st = SISL_NULL; /* First knot vector is pcurv */
125 double *scoef = SISL_NULL; /* Vertices of pcurv */
126 int ikind; /* kind of surface pcurv is */
127 double *rscoef = SISL_NULL; /* Scaled coefficients if pcurv is rational */
128 double wmin, wmax; /* min and max values of the weights if rational */
129 double scale; /* factor for scaling weights if rational */
130 int i; /* loop variable */
131 double *sarray = SISL_NULL; /* Array for calculating denominator if used */
132 int knarr; /* Number of parallel arrays to use. */
133 int nkind; /* Kind of output curve (rcurf). */
134
135 *jstat = 0;
136
137 /* Make local pointers. */
138
139 kn = pcurv->in;
140 kk = pcurv->ik;
141 kdim = pcurv->idim;
142 st = pcurv->et;
143 ikind = pcurv->ikind;
144
145 kdimp1 = kdim + 1;
146
147 /* Test input. */
148
149 if (kdim != idim || (kdim != 2 && kdim != 3))
150 goto err104;
151 if (inarr < 1 || 3 < inarr) goto err172;
152
153 /* rational surfaces are a special case. */
154 if (ikind == 2 || ikind == 4)
155 {
156 kdim++;
157
158 /* scale the coeffs so that min. weight * max. weight = 1. */
159
160 rscoef = pcurv->rcoef;
161 wmin = rscoef[kdim-1];
162 wmax = rscoef[kdim-1];
163
164 for (i = 2*kdim-1; i < kn * kdim; i += kdim)
165 {
166 if (rscoef[i] < wmin)
167 wmin = rscoef[i];
168 if (rscoef[i] > wmax)
169 wmax = rscoef[i];
170 }
171 scale = (double) 1.0 / sqrt (wmin * wmax);
172 scoef = newarray (kn * kdim, DOUBLE);
173 if (scoef == SISL_NULL)
174 goto err101;
175
176 for (i = 0; i < kn * kdim; i++)
177 scoef[i] = rscoef[i] * scale;
178 }
179 else
180 scoef = pcurv->ecoef;
181
182 icurve = newCurve (kn, kk, st, scoef, 1, kdim, 1);
183 if (icurve == SISL_NULL)
184 goto err171;
185
186 icurve->cuopen = pcurv->cuopen;
187
188 if ((ikind == 2 || ikind == 4) && ratflag == 1)
189 {
190 /* Output curve will also be rational. */
191
192 nkind = 2;
193
194 /* Add an extra parallel array to pick up the weights
195 of the subsequent homogeneous vertices of rcurv. */
196
197 knarr = inarr + 1;
198 sarray = new0array (kdimp1 * kdimp1 * knarr, DOUBLE);
199 if (sarray == SISL_NULL) goto err101;
200
201 memcopy (sarray, earray, kdimp1 * kdimp1 * inarr, DOUBLE);
202 sarray[kdimp1 * kdimp1 * knarr - 1] = (DOUBLE) 1.0;
203 }
204 else
205 {
206 nkind = 1;
207 knarr = inarr;
208 sarray = earray;
209 }
210
211 /* Put curve into implicit surface. */
212
213 s1893 (icurve, sarray, kdimp1, knarr, 0, 0, rcurv, &kstat);
214 if (kstat < 0) goto error;
215
216 if (*rcurv == SISL_NULL) goto err171;
217
218 if ( ikind == 2 || ikind == 4 )
219 {
220 /* Free arrays. */
221
222 if (scoef) freearray (scoef);
223 if (ratflag && sarray) freearray (sarray);
224
225 if ( ratflag == 1 )
226 {
227 /* Output from s1893 is a dim+1 non-rational curve. */
228 /* Convert homogeneous curve to rational form (rcoef is SISL_NULL here). */
229
230 (*rcurv)->rcoef = newarray((*rcurv)->in * (*rcurv)->idim, DOUBLE);
231 memcopy((*rcurv)->rcoef, (*rcurv)->ecoef,
232 (*rcurv)->in * (*rcurv)->idim, DOUBLE);
233
234 (*rcurv)->idim --; /* Adjust from the homogeneus coordinates. */
235 (*rcurv)->ikind = 2; /* i.e. rational */
236
237 }
238 }
239
240
241 /* Ok ! */
242
243 goto out;
244
245 /* Error in lower level function. */
246
247 error:
248 *jstat = kstat;
249 s6err ("s1370", *jstat, kpos);
250 goto out;
251
252 /* Allocation problems. */
253
254 err101:
255 *jstat = -101;
256 s6err ("s1370", *jstat, kpos);
257 goto out;
258
259 /* Dimension not equal to 3. */
260
261 err104:
262 *jstat = -104;
263 s6err ("s1370", *jstat, kpos);
264 goto out;
265
266 /* Could not create curve */
267
268 err171:
269 *jstat = -171;
270 s6err ("s1370", *jstat, kpos);
271 goto out;
272
273 /* Dimension inarr not equal to 1,2 or 3. */
274
275 err172:
276 *jstat = -172;
277 s6err ("s1370", *jstat, kpos);
278 goto out;
279
280 out:
281 if (icurve != SISL_NULL) freeCurve (icurve);
282 return;
283 }
284