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: s1174.c,v 1.5 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1174
50
51 #include "sislP.h"
52
53 /*
54 * Forward declarations.
55 * ---------------------
56 */
57 #if defined(SISLNEEDPROTOTYPES)
58 static void s1174_s9corr(double [],double,double,double,double,double,double);
59 static void s1174_s9dir(double *,double *,double []);
60 #else
61 static void s1174_s9corr();
62 static void s1174_s9dir();
63 #endif
64
65 #if defined(SISLNEEDPROTOTYPES)
66 void
s1174(SISLSurf * psurf,double estart[],double eend[],double enext[],double gpos[],int * jstat)67 s1174(SISLSurf *psurf,double estart[],
68 double eend[], double enext[], double gpos[],int *jstat)
69 #else
70 void s1174(psurf,estart,eend,enext,gpos,jstat)
71 SISLSurf *psurf;
72 double estart[];
73 double eend[];
74 double enext[];
75 double gpos[];
76 int *jstat;
77 #endif
78 /*
79 *********************************************************************
80 *
81 *********************************************************************
82 *
83 * PURPOSE : Newton iteration on a onedimensional surface.
84 * The function finds a local extremum.
85 *
86 *
87 * INPUT : psurface - Pointer to the surface.
88 * estart - Start values of parameter intervalls.
89 * eend - End value of parameter intervalls.
90 * enext - Parameter start value for iteration.
91 *
92 *
93 * OUTPUT : gpos - Parameter values of the found extremum.
94 * jstat - status messages
95 * = 1 : Extremum found.
96 * = 0 : Extremum NOT found.
97 * < 0 : error.
98 *
99 *
100 * METHOD : Newton iteration in two parameter direction.
101 *
102 *
103 * REFERENCES :
104 *
105 *
106 * WRITTEN BY : Ulf J. Krystad, SI, OCTOBER 1990
107 * Changed by : Per OEyvind Hvidsten, SINTEF, 11-94
108 * Added an initialization of tdist.
109 *
110 *********************************************************************
111 */
112 {
113 int kstat = 0; /* Local status variable. */
114 int kpos = 0; /* Position of error. */
115 int kleft1=0; /* Variables used in the evaluator. */
116 int kleft2=0; /* Variables used in the evaluator. */
117 int kder=2; /* Order of derivatives to be calulated */
118 int knbit; /* Number of iterations */
119 int kdir; /* Changing direction. */
120 double tdelta[2]; /* Parameter intervals of the surface. */
121 double tdist = 0.0; /* Euclidian norm of derivative vector */
122 double tprev; /* Previous Euclidian norm of derivative vector*/
123 double td[2],t1[2],tdn[2];/* Distances between old and new parameter
124 value in the two parameter directions. */
125 double sval[7]; /* Value ,first and second derivatiev of surf. */
126 double *snorm=sval+7; /* Normal vector of the surface, dummy. */
127 double snext[2]; /* Parameter values */
128 double tol = (double)10000.0*REL_COMP_RES; /* Singularity tolerance */
129 /* --------------------------------------------------------------------- */
130
131 /* Test input. */
132 if (psurf->idim != 1) goto err106;
133
134 /* Fetch endpoints and the intervals of parameter interval of curves. */
135
136 tdelta[0] = psurf->et1[psurf->in1] - psurf->et1[psurf->ik1 - 1];
137 tdelta[1] = psurf->et2[psurf->in2] - psurf->et2[psurf->ik2 - 1];
138
139
140 /* Initiate variables. */
141 gpos[0] = enext[0];
142 gpos[1] = enext[1];
143
144 /* Evaluate 0-2.st derivatives of surface */
145 s1421(psurf,kder,gpos,&kleft1,&kleft2,sval,snorm,&kstat);
146 if (kstat < 0) goto error;
147
148 /* Get Euclidian norm of derivative vector */
149 tprev = sqrt(sval[1]*sval[1] + sval[2]*sval[2]);
150
151 /* Compute the Newton stepdistanse vector. */
152 s1174_s9dir(td,td+1,sval);
153
154 if ( (fabs(td[0]/tdelta[0]) <= REL_COMP_RES) &&
155 (fabs(td[1]/tdelta[1]) <= REL_COMP_RES))
156 goto stop_it;
157
158 /* Adjust if we are not inside the parameter intervall. */
159 t1[0] = td[0];
160 t1[1] = td[1];
161 s1174_s9corr(t1,gpos[0],gpos[1],estart[0],eend[0],estart[1],eend[1]);
162
163 /* Iterate to find the intersection point. */
164
165 for (knbit = 0; knbit < 50; knbit++)
166 {
167 /* Evaluate 0-2.st derivatives of surface */
168
169 snext[0] = gpos[0] + t1[0];
170 snext[1] = gpos[1] + t1[1];
171
172 s1421(psurf,kder,snext,&kleft1,&kleft2,sval,snorm,&kstat);
173 if (kstat < 0) goto error;
174
175 /* Get Euclidian norm of derivative vector */
176 tdist = sqrt(sval[1]*sval[1] + sval[2]*sval[2]);
177
178 /* Compute the Newton stepdistanse vector. */
179 s1174_s9dir(tdn,tdn+1,sval);
180
181 /* Check if the direction of the step have change. */
182
183 kdir = (s6scpr(td,tdn,2) >= DZERO); /* 0 if changed. */
184
185 if (tdist <= tprev || kdir)
186 {
187 /* Ordinary converging. */
188
189 gpos[0] += t1[0];
190 gpos[1] += t1[1];
191
192 td[0] = t1[0] = tdn[0];
193 td[1] = t1[1] = tdn[1];
194
195 /* Adjust if we are not inside the parameter intervall. */
196 s1174_s9corr(t1,gpos[0],gpos[1],estart[0],eend[0],estart[1],eend[1]);
197
198
199 if ( (fabs(t1[0]/tdelta[0]) <= REL_COMP_RES) &&
200 (fabs(t1[1]/tdelta[1]) <= REL_COMP_RES))
201 {
202 gpos[0] += t1[0];
203 gpos[1] += t1[1];
204
205 break;
206 }
207
208 tprev = tdist;
209 }
210
211 else
212 {
213 /* Not converging, half step length try again. */
214
215 t1[0] /= (double)2;
216 t1[1] /= (double)2;
217 /* knbit--; */
218 }
219 }
220
221 /* Iteration stopped, test if point is extremum */
222
223 stop_it:
224
225 if (tdist <= tol)
226 *jstat = 1;
227 else
228 *jstat = 0;
229
230
231 /* Test if the iteration is close to a knot */
232 if (fabs(gpos[0] - psurf->et1[kleft1])/tdelta[0] < tol)
233 gpos[0] = psurf->et1[kleft1];
234 else if (fabs(gpos[0] - psurf->et1[kleft1+1])/tdelta[0] < tol)
235 gpos[0] = psurf->et1[kleft1+1];
236
237 if (fabs(gpos[1] - psurf->et2[kleft2])/tdelta[1] < tol)
238 gpos[1] = psurf->et2[kleft2];
239 else if (fabs(gpos[1] - psurf->et2[kleft2+1])/tdelta[1] < tol)
240 gpos[1] = psurf->et2[kleft2+1];
241
242 /* Iteration completed. */
243 goto out;
244
245 /* --------------------------------------------------------------------- */
246 /* Error in input. Dimension not equal to 1 */
247 err106: *jstat = -106;
248 s6err("s1174",*jstat,kpos);
249 goto out;
250
251 /* Error in lower level routine. */
252 error : *jstat = kstat;
253 s6err("s1174",*jstat,kpos);
254 goto out;
255
256 out:;
257 }
258
259 #if defined(SISLNEEDPROTOTYPES)
260 static void
s1174_s9corr(double gd[],double acoef1,double acoef2,double astart1,double aend1,double astart2,double aend2)261 s1174_s9corr(double gd[], double acoef1,double acoef2,double astart1,
262 double aend1,double astart2, double aend2)
263 #else
264 static void s1174_s9corr(gd,acoef1,acoef2,astart1,aend1,astart2,aend2)
265 double gd[];
266 double acoef1;
267 double acoef2;
268 double astart1;
269 double aend1;
270 double astart2;
271 double aend2;
272 #endif
273 /*
274 *********************************************************************
275 *
276 *********************************************************************
277 *
278 * PURPOSE : To be sure that we are inside the boorder of the
279 * parameter plan. If we are outside clipping is used
280 * to corrigate the step value.
281 *
282 *
283 * INPUT : acoef1 - Coeffisient in the first direction.
284 * acoef2 - Coeffisient in the second direction.
285 * astart1 - The lower boorder in first direction.
286 * aend1 - The higher boorder in first direction.
287 * estart2 - The lower boorder in second direction.
288 * eend2 - The higher boorder in second direction.
289 *
290 *
291 *
292 * INPUT/OUTPUT : gd - Old and new step value.
293 *
294 *
295 * METHOD : We are cutting a line inside a rectangle.
296 * In this case we always know that the startpoint of
297 * the line is inside the rectangel, and we may therfor
298 * use a simple kind of clipping.
299 *
300 *
301 * REFERENCES :
302 *
303 *
304 * WRITTEN BY : Arne Laksaa, SI, Feb 1989
305 *
306 *********************************************************************
307 */
308 {
309 if (acoef1 + gd[0] < astart1) gd[0] = astart1 - acoef1;
310 else if (acoef1 + gd[0] > aend1) gd[0] = aend1 - acoef1;
311
312 if (acoef2 + gd[1] < astart2) gd[1] = astart2 - acoef2;
313 else if (acoef2 + gd[1] > aend2) gd[1] = aend2 - acoef2;
314 }
315
316 #if defined(SISLNEEDPROTOTYPES)
317 static void
s1174_s9dir(double * cdiff1,double * cdiff2,double evals[])318 s1174_s9dir(double *cdiff1, double *cdiff2,double evals[])
319 #else
320 static void s1174_s9dir(cdiff1,cdiff2,evals)
321 double *cdiff1;
322 double *cdiff2;
323 double evals[];
324 #endif
325 /*
326 *********************************************************************
327 *
328 *********************************************************************
329 *
330 * PURPOSE : To compute the step according to a Newton scheme
331 * for finding an extremal to a one dimensional surface.
332 *
333 *
334 * INPUT : evals - Value and derivatives in point on surface.
335 *
336 *
337 * OUTPUT : cdiff1 - Parameter increment in first direction.
338 * cdiff2 - Parameter increment in second direction.
339 *
340 *
341 *
342 * METHOD : This is a one dimensional case. We want to find x,y such that
343 *
344 * x,y : Sx(x,y) = 0
345 * Sy(x,y) = 0
346 *
347 * Using Taylor we get:
348 *
349 * x,y : Sx+DxSxx+DySxy =0
350 * Sy+DySyy+DxSxy =0
351 *
352 * x,y: SxxDx + SxyDy = - Sx
353 * SyyDy + SxyDx = - Sy
354 *
355 * The solution of this matrix equation is the
356 * following function.
357 * No effort is done if the matrix is singular,
358 *
359 * REFERENCES :
360 *
361 *
362 * WRITTEN BY : Ulf J. Krystad, SI, OCTOBER 1990
363 *
364 *********************************************************************
365 */
366 {
367 double tdiv; /* Determinant */
368 double ta11,ta12,ta21,ta22; /* The matrix */
369 double tmax; /* The largest value in matrix */
370 double tb1,tb2; /* The right hand side. */
371 double tderx,tderxx; /* Derivatives */
372 double tdery,tderyy;
373 double tderxy;
374 double tdeltax,tdeltay; /* Locals for the step value to be determined. */
375 /* --------------------------------------------------------------------- */
376
377 /* Init */
378 tderx = evals[1];
379 tdery = evals[2];
380 tderxx = evals[3];
381 tderxy = evals[4];
382 tderyy = evals[5];
383 tdeltax = DZERO;
384 tdeltay = DZERO;
385 *cdiff1 = DZERO;
386 *cdiff2 = DZERO;
387
388
389 /* Building the matrix. */
390
391 ta11 = tderxx;
392 ta12 = tderxy;
393 ta21 = tderxy;
394 ta22 = tderyy;
395 tb1 = -tderx;
396 tb2 = -tdery;
397
398 tmax = max(fabs(ta11),max(fabs(ta12),max(fabs(ta21),fabs(ta22))));
399
400 if (DEQUAL(tb1+tmax,tmax) && DEQUAL(tb2+tmax,tmax))
401 {
402 /* Finished, we have found a max. */
403 }
404 else
405 {
406 tdiv = ta11*ta22 - ta21*ta12;
407 if (fabs(tdiv) > MAX(tmax*REL_COMP_RES,REL_COMP_RES))
408 {
409 /* The matrix is ok, solve the system using Cramers rule. */
410 tdeltax = tb1*ta22 - tb2*ta12;
411 tdeltay = ta11*tb2 - ta21*tb1;
412 tdeltax /= tdiv;
413 tdeltay /= tdiv;
414 }
415 else if (max (fabs(ta11),fabs(ta22)) > REL_COMP_RES)
416 {
417 if (fabs(ta11) > fabs(ta22))
418 tdeltax = tb1/ta11;
419 else
420 tdeltay = tb2/ta22;
421 }
422
423 }
424
425 *cdiff1 = tdeltax;
426 *cdiff2 = tdeltay;
427
428 }
429