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: s1329.c,v 1.3 2001-03-19 15:58:44 afr Exp $
45 *
46 */
47
48
49 #define S1329
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1329(SISLSurf * psold,double epoint[],double enorm[],int idim,SISLSurf ** rsnew,int * jstat)55 s1329(SISLSurf *psold,double epoint[],double enorm[],int idim,
56 SISLSurf **rsnew,int *jstat)
57 #else
58 void s1329(psold,epoint,enorm,idim,rsnew,jstat)
59 SISLSurf *psold;
60 double epoint[];
61 double enorm[];
62 int idim;
63 SISLSurf **rsnew;
64 int *jstat;
65 #endif
66 /*
67 *********************************************************************
68 *
69 *********************************************************************
70 *
71 * PURPOSE : Put the equation of the surface pointed at by psold
72 * into the plane given by the point epoint and the normal
73 * enorm. The result is an equation where the new one-
74 * dimensional surface rsnew is to be equal to zero.
75 *
76 *
77 *
78 * INPUT : psold - Pointer to input surface.
79 * epoint - SISLPoint in the plane.
80 * enorm - Normal to the plane.
81 * idim - Dimension of the space in which the plane lies.
82 *
83 *
84 *
85 * OUTPUT : rsnew - The new one-dimensional surface.
86 * jstat - status messages
87 * > 0 : warning
88 * = 0 : ok
89 * < 0 : error
90 *
91 *
92 * METHOD :
93 *
94 *
95 * REFERENCES :
96 *
97 *-
98 * CALLS : newarray - Allocate space for an array of given type.
99 * freearray - Free space occupied by a given array.
100 * newSurf - Create and initialize new surface.
101 *
102 * WRITTEN BY : Vibeke Skytt, SI, 88-06.
103 * REVISED BY : Mike Floater, SI, 91-04.
104 * CORRECTED BY: Ulf J. Krystad, SI, 91-07.
105 * DEBUGGED BY : Mike Floater, SI, 94-06. Use scSave.
106 *********************************************************************
107 */
108 {
109 int kpos = 0; /* Position of error. */
110 int kdim; /* Dimension of the space in which the output
111 surface lies. */
112 int kn1,kn2; /* Number of coefficients of surface. */
113 int kk1,kk2; /* Order of surface. */
114 int ikind; /* kind of surface psold is */
115 double *scoef = SISL_NULL; /* Coeffecient array of new surface. */
116 double *s1,*s2; /* Pointers used to traverse scoef. */
117 double *sc=SISL_NULL; /* Pointer used to traverse psold->ecoef. */
118 double *scSave=SISL_NULL; /* Pointer to vertices in rational case. */
119 double *rscoef; /* Scaled coefficients if psold is rational */
120 double *s3; /* Stop pointer for each vertex in psold->ecoef. */
121 double *spoint; /* Pointer used to traverse the point epoint. */
122 double *snorm; /* Pointer used to traverse the normal enorm. */
123 double wmin,wmax;/* min and max values of the weights if rational */
124 double scale; /* factor for scaling weights if rational */
125 int i; /* loop variable */
126 int idimp1; /* idim+1 */
127
128 /* Test input. */
129
130 if (idim != psold -> idim) goto err106;
131
132 /* Set simple variables of the new surface. */
133
134 kdim = 1;
135 kn1 = psold -> in1;
136 kn2 = psold -> in2;
137 kk1 = psold -> ik1;
138 kk2 = psold -> ik2;
139 ikind = psold -> ikind;
140
141 /* rational surfaces are a special case */
142 if(ikind == 2 || ikind == 4)
143 {
144 /* scale the coeffs so that min. weight * max. weight = 1 */
145 idimp1=idim+1;
146 rscoef = psold -> rcoef;
147 wmin=rscoef[idim];
148 wmax=rscoef[idim];
149 for(i=idim; i< kn1*kn2*idimp1; i+=idimp1)
150 {
151 if(rscoef[i] < wmin) wmin=rscoef[i];
152 if(rscoef[i] > wmax) wmax=rscoef[i];
153 }
154 scale=1.0/sqrt(wmin*wmax);
155 if ((sc = newarray(idimp1*kn1*kn2,DOUBLE)) == SISL_NULL) goto err101;
156 for(i=0; i< kn1*kn2*idimp1; i++)
157 {
158 sc[i]=rscoef[i]*scale;
159 }
160
161 scSave = sc;
162 }
163 else
164 {
165 sc = psold -> ecoef;
166 }
167
168 /* Allocate space for coeffecient of the new surface. */
169
170 if ((scoef = newarray(kdim*kn1*kn2,DOUBLE)) == SISL_NULL) goto err101;
171
172 /* Compute coefficients of new surface. */
173
174 for (s1=scoef,s2=s1+kn1*kn2; s1<s2; s1++)
175 {
176 *s1 = (double)0.0;
177 if(ikind == 2 || ikind == 4)
178 {
179 /* surface is rational so we're using idim+1 - d homogeneous coords */
180 for (s3=sc+idim,spoint=epoint,snorm=enorm; sc<s3; sc++,spoint++,snorm++)
181 {
182 /* UJK, Turned direction to get right sign in 1D*/
183 /* *s1 += ((*s3)*(*spoint) - *sc)*(*snorm); */
184 *s1 += (*sc - (*s3)*(*spoint))*(*snorm);
185 }
186 sc++;
187 }
188 else
189 {
190 /* surface is not rational so we're using ordinary idim - d coords */
191 for (s3=sc+idim,spoint=epoint,snorm=enorm; sc<s3; sc++,spoint++,snorm++)
192 /* UJK, Turned direction to get right sign in 1D*/
193 /* *s1 += (*spoint - *sc)*(*snorm); */
194 *s1 += (*sc - *spoint)*(*snorm);
195 }
196 }
197
198 if(ikind == 2 || ikind == 4) freearray(scSave);
199
200 /* Create output surface. */
201
202 *rsnew = newSurf(kn1,kn2,kk1,kk2,psold->et1,psold->et2,scoef,1,kdim,1);
203 if (*rsnew == SISL_NULL) goto err101;
204
205 /* Task done. */
206
207 *jstat = 0;
208 goto out;
209
210 /* Error in space allocation. */
211
212 err101:
213 *jstat = -101;
214 s6err("s1329",*jstat,kpos);
215 goto out;
216
217 /* Error in input. Confliction dimensions. */
218
219 err106 :
220 *jstat = -106;
221 s6err("s1329",*jstat,kpos);
222 goto out;
223
224 out:
225 /* Free space allocated for local array. */
226
227 if (scoef != SISL_NULL) freearray(scoef);
228 return;
229 }
230