/********************************************************************/ /* */ /* This file is part of the VARKON Geometry Library. */ /* URL: http://www.varkon.com */ /* */ /* This library is free software; you can redistribute it and/or */ /* modify it under the terms of the GNU Library General Public */ /* License as published by the Free Software Foundation; either */ /* version 2 of the License, or (at your option) any later */ /* version. */ /* */ /* This library is distributed in the hope that it will be */ /* useful, but WITHOUT ANY WARRANTY; without even the implied */ /* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR */ /* PURPOSE. See the GNU Library General Public License for more */ /* details. */ /* */ /* You should have received a copy of the GNU Library General */ /* Public License along with this library; if not, write to the */ /* Free Software Foundation, Inc., 675 Mass Ave, Cambridge, */ /* MA 02139, USA. */ /* */ /* (C)Microform AB 1984-1999, Gunnar Liden, gunnar@microform.se */ /* */ /********************************************************************/ #include "../../DB/include/DB.h" #include "../include/GE.h" /********************************************************************/ /*! */ /* Function: varkon_sur_num2 File: sur902.c */ /* ============================================================= */ /* */ /* Purpose */ /* ------- */ /* */ /* The function is a general function for the calculation */ /* of points on a surface, i.e a function which numerically */ /* finds the solution to a two variable problem. */ /* */ /* Input is a pointer to a structure which defines the */ /* calculation criterion for the points (the objective */ /* function), the method for the numerical calculation, */ /* the number of restarts, the tolerances, etc. */ /* */ /* The problem is solved as a non-linear optimization problem. */ /* The gradient method, with parabola interpolation for the */ /* one-dimensional search, is used. */ /* The convergence is normally slow using the gradient and */ /* the search direction is modified with the Hessian matrix. */ /* */ /* Author: Gunnar Liden */ /* */ /* Revisions */ /* */ /* 1994-06-09 Originally written */ /* 1995-06-07 LINTSUR - lintsur added */ /* 1995-08-27 ntol -> TOL7 for LINTSUR istart=5,6,7,8,9 */ /* core for too many solutions !!!! ????? */ /* 1995-09-16 pat_tra for LINTSUR */ /* 1995-09-30 Bug: hstep not initialized for _newst */ /* Minimum step TOL1/100 --> TOL1/10000 */ /* 1995-10-23 Tests: lintsur X value for angle computation */ /* 1995-12-31 stepmin, lintsur objf_crit and tolerances */ /* 1997-02-27 Bug for boundary */ /* 1999-11-27 Free source code modifications */ /* */ /* !*/ /********************************************************************/ /* ------------- Short description of function -----------------*/ /* */ /*sdescr varkon_sur_num2 Surface two param. calculation */ /* */ /*------------------------------------------------------------- */ /*!New-Page--------------------------------------------------------!*/ /*!-------------- Function calls (internal) ------------------------*/ /* */ #ifdef DEBUG static short initial(); /* Check and initiations For Debug On */ #endif static short retrieve(); /* Retrieve input to local variables */ static short newstart(); /* Calculate a new UV start point */ static short sdirect(); /* One-dimensional search direction */ static short stepleng(); /* Step length for one-dimen. search */ static short nearpt(); /* Objective, gradient.. for CLOSEPT */ static short lintsur(); /* Objective, gradient.. for LINTSUR */ static short parab3p(); /* Parabola interp.: 3 pts */ static short addsol(); /* Add a non-existent solution. */ static DBfloat stepmin(); /* Minimum U,V step length */ /* */ /*-----------------------------------------------------------------!*/ /*!-------------- Flow diagram -------------------------------------*/ /* _________________ */ /* ! ! */ /* ! varkon_sur_num2 ! */ /* ! (sur902) ! */ /* !_________________! */ /* ____________________!__________________ */ /* _!_ _!_ _!_ _!_ _!_ _!_ _!_ */ /* ! ! ! ! ! ! ! ! ! ! ! ! ! ! */ /* ! 1 ! ! 2 ! L ! 3 ! L ! 4 ! ! 5 ! ! 6 ! ! 7 ! */ /* ! ! ! ! a ! ! a ! ! ! ! ! ! ! ! */ /* ! i ! ! r ! b ! n ! b ! O ! ! a ! ! s ! ! S ! */ /* ! n ! ! e ! e ! e ! e ! b ! ! d ! ! d ! ! t ! */ /* ! i ! ! t ! l ! w ! l ! j ! ! d ! ! i ! ! e ! */ /* ! t ! ! r ! ! s ! ! e ! ! s ! ! r ! ! p ! */ /* ! i ! ! i ! n ! t ! n ! c ! ! o ! ! e ! ! ! */ /* ! a ! ! e ! e ! a ! o ! t ! ! l ! ! c ! !___! */ /* ! l ! ! v ! w ! r ! p ! i ! !___! ! t ! ! */ /* !___! ! e ! s ! t ! t ! v ! !___! ! */ /* !___! t !___! ! e ! ! */ /* !___! ! */ /* __________! !_____________ */ /* ___!____ ___!____ ___!______ ____!____ */ /* ! nearpt ! !lintsur ! ! stepleng ! ! parab3p ! */ /* !________! !________! !__________! !_________! !*/ /* */ /*------------------------------------------------------------------*/ /*!New-Page--------------------------------------------------------!*/ /*!-- Static (common) variables for the functions in this file -----*/ /* */ static DBint ctype; /* Case of computation */ static DBint method; /* Method which shall be used */ static DBint nstart; /* The number of restarts */ static DBint istart; /* The current restart number: */ /* Eq. -1: No restart */ /* Eq. 0: Start with input s_uvpt*/ /* Eq. 1: Start with mid pt */ /* Eq. 2: Start with corner pt 1 */ /* (Refer to newstart .. ) */ static DBint maxiter; /* Maximum number of iterations */ static DBTmat pat_tra; /* Transform. matrix for patches */ static DBint nopt; /* Number of optimiz. iterations */ static DBVector extpt; /* External point */ static DBVector v_proj; /* Projection vector (LINTSUR) */ static DBVector s_uvpt; /* Start UV point */ static DBint uvstart; /* Eq. 1: s_uvpt is defined */ static DBfloat us,vs,ue,ve; /* The search area (a patch) */ static DBfloat u_pat,v_pat; /* Current U,V point on surface */ static DBfloat objf; /* Objective function */ static DBfloat objf_crit; /* Objective fctn end criterion */ static DBfloat c_zero; /* The zero (0.0) end criterion */ static DBfloat c_zero2; /* Optional c_zero */ static DBfloat c_infi; /* Infinite criterion */ static DBfloat c_denom; /* Criterion for zero denominator */ static DBfloat ctol; /* Criterion for identical points */ static DBfloat comptol; /* Computer tolerance (accuracy) */ static DBfloat dobjdu,dobjdv; /* Gradient */ static DBfloat d2objdu2,d2objdv2;/* Elements in the Hessian matrix */ static DBfloat d2objdudv; /* Element in the Hessian matrix */ static DBfloat dmdu,dmdv; /* Search direction */ static DBfloat hstep; /* Step U,V for one-dimen. search */ static DBfloat hstep_r3min; /* Minimum R*3 hstep */ static DBfloat gradl; /* Length of gradient vector */ static DBfloat deterh; /* Determinant of Hessian matrix */ static DBfloat obj_1; /* Objective function value 1 */ static DBfloat obj_2; /* Objective function value 2 */ static DBfloat obj_3; /* Objective function value 3 */ static DBint numsol; /* Number of solutions (=*p_npts) */ static DBint icase; /* Case for varkon_sur_eval */ static EVALS xyz; /* Coordinates and derivatives */ /* for a point on a surface */ /*-----------------------------------------------------------------!*/ static DBint n_iter; /* For the maxiter check !Temp. static */ /*!-------------- Function calls (external) ------------------------*/ /* */ /* varkon_sur_eval * Surface evaluation routine */ /* varkon_evals_transf * Transformate EVALS */ /* varkon_angd * Angle in deg. between vectors */ /* varkon_comptol * Retrieve computer tolerance */ /* varkon_erpush * Error message to terminal */ /* */ /*-----------------------------------------------------------------!*/ /*!------------ Error messages and warnings ------------------------*/ /* */ /* SU2943 = Called function xxxxxx failed in varkon_sur_num2 */ /* SU2983 = sur902 Illegal computation case= for varkon_sur_num2 */ /* SU2973 = Internal function xxxxxx failed in varkon_sur_num2 */ /* SU2993 = Severe program error in varkon_sur_num2 (sur902). */ /* */ /*-----------------------------------------------------------------!*/ /*!****************** Function **************************************/ /* */ DBstatus varkon_sur_num2 ( /*-------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBSurf *p_sur, /* Surface (ptr) */ DBPatch *p_pat, /* Alloc. area for topol. patch data (ptr) */ IRUNON *p_comp, /* Computation data (ptr) */ SPOINT *p_spt, /* Current record in point table (ptr) */ /* Out: */ DBint *p_npts, /* Number of output points (ptr) */ EVALS xyz_a[SMAX] ) /* Array with all solutions (R*3 and UV) */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /* */ /*----------------------------------------------------------------- */ char errbuf[80]; /* String for error message fctn erpush */ DBint status; /* Error code from a called function */ /*--------------end-of-declarations---------------------------------*/ /*!New-Page--------------------------------------------------------!*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ /*! */ /* 1. Check of input data and initializations */ /* __________________________________________ */ /* */ /* Initialize output coordinates and derivatives For DEBUG On. */ /* Call of initial, which also prints out the input data. */ /* !*/ #ifdef DEBUG status=initial(p_spt,p_comp,p_npts); if (status<0) { sprintf(errbuf,"initial%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif /*! */ /* 2. Retrieve the input data */ /* __________________________ */ /* */ /* Retrieve input data to local static variables. */ /* Call of retrieve. */ /* !*/ status=retrieve (p_spt,p_comp); #ifdef DEBUG if (status<0) { sprintf(errbuf,"retrieve%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif /*! */ /* 3. Restart (or start) of calculation */ /* ____________________________________ */ /* */ _newst: /* */ /* Calculate a (new) UV start point. */ /* Call of newstart if istart <= nstart */ /* Goto _exit if istart > nstart */ /* !*/ if ( istart > nstart ) goto _exit; status=newstart (); #ifdef DEBUG if (status<0) { sprintf(errbuf,"newstart%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif /*! */ /* 4. Calculate objective function */ /* _______________________________ */ /* */ n_iter = 0; /* */ _nopt: /* Label: New optimization iteration */ /* !*/ n_iter = n_iter + 1 ; if ( n_iter > maxiter ) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902 Number iterations n_iter= %d > maxiter= %d\n" ,n_iter , maxiter ); fprintf(dbgfil(SURPAC), "sur902 Temporary fix: Make a new start \n"); } #endif sprintf(errbuf,"(maxiter)%%varkon_sur_num2 (sur902)"); goto _newst; /* Temporary fix return(varkon_erpush("SU2993",errbuf)); */ } /*! */ /* Check that point u_pat,v_pat is within patch limit us,vs,ue,ve */ /* Calculate objective function, gradient and gradient length. */ /* Call of nearpt for ctype= CLOSEPT */ /* !*/ if ( u_pat >= us && u_pat <= ue && v_pat >= vs && v_pat <= ve ) { ; } else { u_pat = u_pat-hstep*dmdu; v_pat = v_pat-hstep*dmdv; nopt= -10; /* Calculate a new direction */ } if ( ctype == CLOSEPT ) { status= nearpt (p_sur,p_pat); #ifdef DEBUG if (status<0) { sprintf(errbuf,"nearpt%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif } else if ( ctype == LINTSUR ) { status= lintsur (p_sur,p_pat); #ifdef DEBUG if (status<0) { sprintf(errbuf,"lintsur%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif } #ifdef DEBUG if ( objf_crit < 0.1 ) { if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902 n_iter= %d u_pat %7.4f v_pat %7.4f objf_crit %15.12f nopt %d\n" ,n_iter,u_pat , v_pat , objf_crit , nopt ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902 n_iter= %d dobjdu %f dobjdv %f\n" ,n_iter, dobjdu , dobjdv ); } } #endif /*! */ /* 5. Optimal point ? */ /* ___________________ */ /* */ /* There is an optimal point if objf_crit < c_zero : */ /* i. Add the solution. Call of addsol. */ /* ii. Make a new start. Goto _newst. */ /* !*/ if ( objf_crit <= c_zero ) { status= addsol (p_sur,p_pat,xyz_a); #ifdef DEBUG if (status<0) { sprintf(errbuf,"addsol%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif goto _newst; } /*! */ /* 6. Search direction */ /* ___________________ */ /* */ /* Calculate the one-dimensional search direction */ /* Call of sdirect if nopt < 0. */ /* !*/ if ( nopt < 0 ) status=sdirect(); #ifdef DEBUG if (status<0) { sprintf(errbuf,"sdirect%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif /*!New-Page--------------------------------------------------------!*/ /*! */ /* 7. Step length for the one-dimensional search */ /* ______________________________________________ */ /* */ /* Calculate or adjust the step length hstep. */ /* Call of stepleng if nopt < 0. */ /* Solution is outside patch <==> hstep= 0 ==> Goto _newst. */ /* (Search direction in the current point) */ /* !*/ if ( nopt < 0 ) status=stepleng (); #ifdef DEBUG if (status<0) { sprintf(errbuf,"stepleng%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } #endif #ifdef DEBUG if ( objf_crit < 0.1 && dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902 n_iter= %d hstep %18.12f dmdu %f dmdv %f\n" ,n_iter, hstep, dmdu , dmdv ); } if ( dbglev(SURPAC) == 2 && fabs(hstep) < 100.0*stepmin() ) { fprintf(dbgfil(SURPAC), "sur902 !hstep! %18.15f stepmin %18.15f \n" , fabs(hstep), stepmin() ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } if ( dbglev(SURPAC) == 1 && fabs(hstep) < stepmin() ) { fprintf(dbgfil(SURPAC), "sur902 !hstep! %18.15f < stepmin %18.15f ==> Goto _newst\n" , fabs(hstep), stepmin() ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif if ( fabs(hstep) < stepmin() ) { hstep = -c_infi; goto _newst; } /*! */ /* If nopt < 0 (search direction and step length has been calc.) */ /* i. Let obj_1= objf */ /* ii. Let nopt= 1 */ /* iii. Let u_pat= u_pat+hstep*dmdv v_pat= v_pat+hstep*dmdv */ /* iv. Goto _nopt; */ /* !*/ if ( nopt < 0 ) { nopt = 1; obj_1 = objf; u_pat = u_pat+hstep*dmdu; v_pat = v_pat+hstep*dmdv; goto _nopt; } /*! */ /* If nopt = 1 and objf < obj_1 */ /* i. Let obj_2= objf */ /* ii. Let nopt= 2 */ /* iii. Let u_pat= u_pat+hstep*dmdv v_pat= v_pat+hstep*dmdv */ /* iv. Goto _nopt; */ /* Else if nopt = 1 and objf >= obj_1 */ /* i. hstep = hstep/2 (if hstep > 0.1*stepmin) */ /* ii. Goto _nopt; */ /* !*/ if ( nopt == 1 && objf < obj_1 ) { nopt = 2; obj_2 = objf; u_pat = u_pat+hstep*dmdu; v_pat = v_pat+hstep*dmdv; goto _nopt; } else if ( nopt == 1 && objf >= obj_1 ) { hstep = hstep/2; #ifdef DEBUG if ( dbglev(SURPAC) == 1 && fabs(hstep) < 100.0*stepmin() ) { fprintf(dbgfil(SURPAC), "sur902 hstep %18.15f divided by two (2) and nopt = -10\n", fabs(hstep) ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif /* Ugly programmed fix that only partly solves a problem */ u_pat = u_pat-hstep*dmdu; v_pat = v_pat-hstep*dmdv; nopt = -10; if ( fabs(hstep) < 0.1*stepmin() ) { sprintf(errbuf,"(hstep/2)%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } goto _nopt; } /*! */ /* If nopt >= 2 and objf < obj_2 */ /* i. Let obj_1= obj_2 and obj_2= objf */ /* ii. Let nopt= nopt + 1 */ /* iii. Let u_pat= u_pat+hstep*dmdv v_pat= v_pat+hstep*dmdv */ /* iv. Goto _nopt; */ /* Else if nopt >= 2 and objf >= obj_2 */ /* i. Calculate hmin. Call of parab3p. */ /* ii. Let nopt = -10 */ /* iii. Goto _nopt. */ /* !*/ if ( nopt >= 2 && objf < obj_2 ) { nopt = nopt + 1; obj_1 = obj_2; obj_2 = objf; u_pat = u_pat+hstep*dmdu; v_pat = v_pat+hstep*dmdv; goto _nopt; } else if ( nopt >= 2 && objf >= obj_2 ) { obj_3 = objf; status=parab3p(); if (status<0) { sprintf(errbuf,"parab3p%%varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2973",errbuf)); } nopt = -10; goto _nopt; } /*! */ _exit: /* Last restart */ /* !*/ *p_npts = numsol; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902 Exit*varkon_sur_num2 us %6.2f ue %6.2f vs %6.2f ve %6.2f n %d\n" ,p_spt->us , p_spt->ue , p_spt->vs , p_spt->ve , *p_npts); } fflush(dbgfil(SURPAC)); #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ #ifdef DEBUG /*!********* Internal ** function ** initial **For*Debug*On**********/ /* */ /* Purpose */ /* -------- */ /* */ /* The function checks the input data and it initializes the */ /* output variables and the static (common) variables to the */ /* values 1.23456789 and 123456789. */ static short initial(p_spt,p_comp,p_npts) /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ SPOINT *p_spt; /* Current record in point table (ptr) */ IRUNON *p_comp; /* Computation data (ptr) */ DBint *p_npts; /* Number of output points (ptr) */ /* Out: */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ char errbuf[80]; /* String for error message fctn erpush */ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* 1. Printout of input variables */ /* ============================== */ /* !*/ if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "Enter ***** varkon_sur_num2 (sur902) ** initial ******\n"); fprintf(dbgfil(SURPAC), "sur902*initial Comp. case %d Method %d Restrts %d Max iter %d\n", p_comp->pro_un.ctype, p_comp->pro_un.method, p_comp->pro_un.nstart, p_comp->pro_un.maxiter); if ( p_comp->pro_un.ctype == CLOSEPT ) { fprintf(dbgfil(SURPAC), "sur902*initial Case CLOSEPT: External point= %10.2f %10.2f %10.2f \n", p_spt->extpt.x_gm, p_spt->extpt.y_gm, p_spt->extpt.z_gm); } if ( p_comp->pro_un.ctype == LINTSUR ) { fprintf(dbgfil(SURPAC), "sur902*initial Case LINTSUR: External point= %10.2f %10.2f %10.2f \n", p_spt->extpt.x_gm, p_spt->extpt.y_gm, p_spt->extpt.z_gm); fprintf(dbgfil(SURPAC), "sur902*initial Case LINTSUR: Projection vec= %10.2f %10.2f %10.2f \n", p_spt->proj.x_gm, p_spt->proj.y_gm, p_spt->proj.z_gm); } fprintf(dbgfil(SURPAC), "sur902*initial startuvpt %8.4f %8.4f %8.4f\n", p_spt->startuvpt.x_gm , p_spt->startuvpt.y_gm , p_spt->startuvpt.z_gm ); fprintf(dbgfil(SURPAC), "sur902*initial Patch us %8.4f vs %8.4f ue %8.4f ve %8.4f\n", p_spt->us , p_spt->vs , p_spt->ue , p_spt->ve ); } fflush(dbgfil(SURPAC)); /* To file from buffer */ /*! */ /* 2. Check of input data */ /* ====================== */ /* */ /* Check the UV area. Mininum rectangle sides are 10*TOL4. */ /* !*/ if ( fabs(p_spt->ue - p_spt->us) < 10.0*TOL4 ) { sprintf(errbuf,"(ue-us<10*TOL4)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } else if ( fabs(p_spt->ve - p_spt->vs) < 10.0*TOL4 ) { sprintf(errbuf,"(ve-vs<10*TOL4)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* The number of maxiter must be between 10 and 100. */ /* The number of restarts must be between 2 and maxiter. */ /* !*/ if (p_comp->pro_un.maxiter>= 10 || p_comp->pro_un.maxiter<= 100 ) { ; } else { sprintf(errbuf,"(maxiter)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } if (p_comp->pro_un.nstart > 1 || p_comp->pro_un.nstart <= p_comp->pro_un.maxiter ) { ; } else { sprintf(errbuf,"(nstart)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } if ( p_comp->pro_un.ctype == CLOSEPT ) { ; } else if ( p_comp->pro_un.ctype == LINTSUR ) { ; } else { sprintf(errbuf,"(ctype)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* 3. Initiate output variables and static (internal) variables */ /* ============================================================ */ /* !*/ *p_npts = -123456789; us = -1.23456789; ue = -1.23456789; vs = -1.23456789; ve = -1.23456789; ctype = -12345; method = -12345; nstart = -12345; maxiter = -12345; obj_1 = -1.23456789; /* Objective function value 1 */ obj_2 = -1.23456789; /* Objective function value 2 */ obj_3 = -1.23456789; /* Objective function value 3 */ return(SUCCED); } /* End of function */ /********************************************************************/ #endif /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** nearpt ********************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calculates the function f_cur and dfdt */ /* for calculation case: CLOSEPT */ /* */ static short nearpt(p_sur,p_pat) /* */ /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: extpt, u_pat,v_pat */ DBSurf *p_sur; /* Surface (ptr) */ DBPatch *p_pat; /* Alloc. area for topol. patch data (ptr) */ /* Out: Objective function objf, gradient and Hessian matrix */ /* */ /*-----------------------------------------------------------------!*/ /*!------------------- Theory --------------------------------------*/ /* */ /* Objective function objf: */ /* ------------------------ */ /* */ /* The objective function is the distance from an external */ /* point extpt to a point on the surface. The u,v point which */ /* gives the minimum distance is the solution to the problem. */ /* In this point is the surface normal parallell with the */ /* vector from the surface point to the external point extpt. */ /* */ /* objf(u,v) = !r(u,v)-extpt!**2 = */ /* ( xyz.r_x-extpt.x_gm)**2 + */ /* ( xyz.r_y-extpt.y_gm)**2 + */ /* ( xyz.r_z-extpt.z_gm)**2 */ /* */ /* */ /* dobjdu = d(objf)/du = */ /* 2*( xyz.r_x-extpt.x_gm)* xyz.u_x + */ /* 2*( xyz.r_y-extpt.y_gm)* xyz.u_y + */ /* 2*( xyz.r_z-extpt.z_gm)* xyz.u_z */ /* */ /* */ /* dobjdv = d(objf)/dv = */ /* 2*( xyz.r_x-extpt.x_gm)* xyz.v_x + */ /* 2*( xyz.r_y-extpt.y_gm)* xyz.v_y + */ /* 2*( xyz.r_z-extpt.z_gm)* xyz.v_z */ /* */ /* */ /* */ /* d2objdu2 = d2(objf)/du2 = */ /* 2*xyz.u_x* xyz.u_x + */ /* 2*( xyz.r_x-extpt.x_gm)* xyz.u2_x + */ /* 2*xyz.u_y* xyz.u_y + */ /* 2*( xyz.r_y-extpt.y_gm)* xyz.u2_y + */ /* 2*xyz.u_z* xyz.u_z + */ /* 2*( xyz.r_z-extpt.z_gm)* xyz.u2_z */ /* !*/ /*!New-Page--------------------------------------------------------!*/ /*! */ /* d2objdv2 = d2(objf)/dv2 = */ /* 2*xyz.v_x* xyz.v_x + */ /* 2*( xyz.r_x-extpt.x_gm)* xyz.v2_x + */ /* 2*xyz.v_y* xyz.v_y + */ /* 2*( xyz.r_y-extpt.y_gm)* xyz.v2_y + */ /* 2*xyz.v_z* xyz.v_z + */ /* 2*( xyz.r_z-extpt.z_gm)* xyz.v2_z */ /* */ /* */ /* */ /* d2objdudv = d2(objf)/dudv = */ /* 2*xyz.v_x* xyz.u_x + */ /* 2*( xyz.r_x-extpt.x_gm)* xyz.uv_x + */ /* 2*xyz.v_y* xyz.u_y + */ /* 2*( xyz.r_y-extpt.y_gm)* xyz.uv_y + */ /* 2*xyz.v_z* xyz.u_z + */ /* 2*( xyz.r_z-extpt.z_gm)* xyz.uv_z */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ DBfloat snorm[3]; /* Surface normal */ DBfloat vextp[3]; /* Vector to external point */ DBfloat angle; /* Angle between vectors */ DBint status; /* Error code from a called function */ char errbuf[80]; /* String for error message fctn erpush */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* 1. Coordinates and derivatives */ /* ============================== */ /* */ /* Evaluate coordinates and derivatives for u_pat,v_pat */ /* Call of varkon_sur_eval (sur210). */ /* !*/ status=varkon_sur_eval (p_sur,p_pat,icase,u_pat,v_pat,&xyz); #ifdef DEBUG if (status<0) { sprintf(errbuf,"varkon_sur_eval%%varkon_sur_num2 (nearpt)"); return(varkon_erpush("SU2943",errbuf)); } #endif /*! */ /* Objective function value objf. */ /* !*/ objf= ( xyz.r_x-extpt.x_gm )*( xyz.r_x-extpt.x_gm ) + ( xyz.r_y-extpt.y_gm )*( xyz.r_y-extpt.y_gm ) + ( xyz.r_z-extpt.z_gm )*( xyz.r_z-extpt.z_gm ); /*! */ /* Gradient= (dobjdu,dobjdv) */ /* !*/ dobjdu= 2.0*( xyz.r_x-extpt.x_gm )*xyz.u_x + 2.0*( xyz.r_y-extpt.y_gm )*xyz.u_y + 2.0*( xyz.r_z-extpt.z_gm )*xyz.u_z; dobjdv= 2.0*( xyz.r_x-extpt.x_gm )*xyz.v_x + 2.0*( xyz.r_y-extpt.y_gm )*xyz.v_y + 2.0*( xyz.r_z-extpt.z_gm )*xyz.v_z; /*! */ /* ! d2objdu2 d2objdudv ! */ /* Hessian= ! ! */ /* ! d2objdudv d2objdv2 ! */ /* !*/ d2objdu2 = 2.0* xyz.u_x* xyz.u_x + 2.0*( xyz.r_x-extpt.x_gm)* xyz.u2_x + 2.0* xyz.u_y* xyz.u_y + 2.0*( xyz.r_y-extpt.y_gm)* xyz.u2_y + 2.0* xyz.u_z* xyz.u_z + 2.0*( xyz.r_z-extpt.z_gm)* xyz.u2_z; d2objdv2 = 2.0* xyz.v_x* xyz.v_x + 2.0*( xyz.r_x-extpt.x_gm)* xyz.v2_x + 2.0* xyz.v_y* xyz.v_y + 2.0*( xyz.r_y-extpt.y_gm)* xyz.v2_y + 2.0* xyz.v_z* xyz.v_z + 2.0*( xyz.r_z-extpt.z_gm)* xyz.v2_z; d2objdudv = 2.0* xyz.v_x* xyz.u_x + 2.0*( xyz.r_x-extpt.x_gm)* xyz.uv_x + 2.0* xyz.v_y* xyz.u_y + 2.0*( xyz.r_y-extpt.y_gm)* xyz.uv_y + 2.0* xyz.v_z* xyz.u_z + 2.0*( xyz.r_z-extpt.z_gm)* xyz.uv_z; /*! */ /* Gradient vector length gradl. */ /* !*/ gradl= SQRT(dobjdu*dobjdu + dobjdv*dobjdv); /*! */ /* Angle between surface normal and vector to extpt. */ /* Call of varkon_angd (angd). */ /* Let objf_crit= angle between the two vectors. */ /* !*/ snorm[0] = xyz.n_x; snorm[1] = xyz.n_y; snorm[2] = xyz.n_z; vextp[0] = extpt.x_gm - xyz.r_x; vextp[1] = extpt.y_gm - xyz.r_y; vextp[2] = extpt.z_gm - xyz.r_z; if ( objf > 2*ctol ) { status=varkon_angd (snorm , vextp, &angle ); if (status<0) { sprintf(errbuf,"angd%%varkon_sur_num2 (nearpt)"); return(varkon_erpush("SU2943",errbuf)); } if ( angle > 90.0 ) { objf_crit = 180.0 - angle; } else { objf_crit = angle; } } else /* Vector does not exist Added 1995-08-26 */ /* Hopefully will gradient length be zero !!!! */ { objf_crit = SQRT(objf); c_zero = c_zero2; /* Change to coordinate criterion */ } #ifdef DEBUG /* if ( dbglev(SURPAC) == 1 && objf_crit < 1.0 ) ODIN */ if ( dbglev(SURPAC) == 1 ) { if ( angle > 90.0 ) { fprintf(dbgfil(SURPAC), "sur902*nearpt Dist %8.1f Angle %f u_pat %8.4f v_pat %8.4f \n" , -SQRT(objf) , objf_crit , u_pat , v_pat ); } else { fprintf(dbgfil(SURPAC), "sur902*nearpt Dist %8.1f Angle %f u_pat %8.4f v_pat %8.4f \n" , SQRT(objf) , objf_crit , u_pat , v_pat ); } } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*nearpt objf %f dobjdu %f dobjdv %f\n" , objf , dobjdu, dobjdv ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*nearpt objf %10.1f gradl %f u_pat %8.4f v_pat %8.4f \n" , objf , gradl , u_pat , v_pat ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*nearpt d2objdu2 %f d2objdv2 %f d2objdudv %f \n" , d2objdu2, d2objdv2 , d2objdudv ); } fflush(dbgfil(SURPAC)); /* To file from buffer */ #endif return(SUCCED); } /* End of function */ /*!****************************************************************!*/ /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** lintsur *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calculates the function f_cur and dfdt */ /* for calculation case: LINTSUR */ /* */ static short lintsur (p_sur,p_pat) /* */ /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: extpt, (v_proj), u_pat,v_pat */ DBSurf *p_sur; /* Surface (ptr) */ DBPatch *p_pat; /* Alloc. area for topol. patch data (ptr) */ /* Out: Objective function objf, gradient and Hessian matrix */ /* */ /*-----------------------------------------------------------------!*/ /*!------------------- Theory --------------------------------------*/ /* */ /* Objective function objf: */ /* ------------------------ */ /* */ /* The objective function is the distance in the YZ plane */ /* to the point (0,0,0). The surface has been transformated... */ /* The u,v point which gives the minimum distance is the solution */ /* to the problem. This distance is zero ..... */ /* In this point is the X axis parallell with the */ /* vector from the surface point to the external point (0,0,0). */ /* */ /* objf(u,v) = !r(u,v)!**2 in YZ plane = */ /* ( xyz.r_y)**2 + */ /* ( xyz.r_z)**2 */ /* */ /* */ /* dobjdu = d(objf)/du = */ /* 2*( xyz.r_y)* xyz.u_y + */ /* 2*( xyz.r_z)* xyz.u_z */ /* */ /* */ /* dobjdv = d(objf)/dv = */ /* 2*( xyz.r_y)* xyz.v_y + */ /* 2*( xyz.r_z)* xyz.v_z */ /* */ /* */ /* */ /* d2objdu2 = d2(objf)/du2 = */ /* 2*xyz.u_y* xyz.u_y + */ /* 2*( xyz.r_y)* xyz.u2_y + */ /* 2*xyz.u_z* xyz.u_z + */ /* 2*( xyz.r_z)* xyz.u2_z */ /* !*/ /*!New-Page--------------------------------------------------------!*/ /*! */ /* d2objdv2 = d2(objf)/dv2 = */ /* 2*xyz.v_y* xyz.v_y + */ /* 2*( xyz.r_y)* xyz.v2_y + */ /* 2*xyz.v_z* xyz.v_z + */ /* 2*( xyz.r_z)* xyz.v2_z */ /* */ /* */ /* */ /* d2objdudv = d2(objf)/dudv = */ /* 2*xyz.v_y* xyz.u_y + */ /* 2*( xyz.r_y)* xyz.uv_y + */ /* 2*xyz.v_z* xyz.u_z + */ /* 2*( xyz.r_z)* xyz.uv_z */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ EVALS xyz_tra; /* Transformed oordinates and derivatives */ /* for a point on a surface */ /*----------------------------------------------------------------- */ DBint status; /* Error code from a called function */ char errbuf[80]; /* String for error message fctn erpush */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* 1. Coordinates and derivatives */ /* ============================== */ /* */ /* Evaluate coordinates and derivatives for u_pat,v_pat */ /* Call of varkon_sur_eval (sur210). */ /* !*/ status=varkon_sur_eval (p_sur,p_pat,icase,u_pat,v_pat,&xyz); #ifdef DEBUG if (status<0) { sprintf(errbuf,"varkon_sur_eval%%varkon_sur_num2 (lintsur)"); return(varkon_erpush("SU2943",errbuf)); } #endif /*! */ /* Transformate coordinates and derivatives to "xaxis=line system" */ /* Call of varkon_evals_transf (sur640). */ /* !*/ status=varkon_evals_transf ( &xyz, &pat_tra, &xyz_tra ); if (status<0) { sprintf(errbuf,"varkon_evals_transf%%varkon_sur_num2 (lintsur)"); return(varkon_erpush("SU2943",errbuf)); } /*! */ /* Objective function value objf. */ /* !*/ objf= ( xyz_tra.r_y)*( xyz_tra.r_y) + ( xyz_tra.r_z)*( xyz_tra.r_z); /*! */ /* Gradient= (dobjdu,dobjdv) */ /* !*/ dobjdu= 2.0*( xyz_tra.r_y)*xyz_tra.u_y + 2.0*( xyz_tra.r_z)*xyz_tra.u_z; dobjdv= 2.0*( xyz_tra.r_y)*xyz_tra.v_y + 2.0*( xyz_tra.r_z)*xyz_tra.v_z; /*! */ /* ! d2objdu2 d2objdudv ! */ /* Hessian= ! ! */ /* ! d2objdudv d2objdv2 ! */ /* !*/ d2objdu2 = 2.0* xyz_tra.u_y* xyz_tra.u_y + 2.0*( xyz_tra.r_y)* xyz_tra.u2_y + 2.0* xyz_tra.u_z* xyz_tra.u_z + 2.0*( xyz_tra.r_z)* xyz_tra.u2_z; d2objdv2 = 2.0* xyz_tra.v_y* xyz_tra.v_y + 2.0*( xyz_tra.r_y)* xyz_tra.v2_y + 2.0* xyz_tra.v_z* xyz_tra.v_z + 2.0*( xyz_tra.r_z)* xyz_tra.v2_z; d2objdudv = 2.0* xyz_tra.v_y* xyz_tra.u_y + 2.0*( xyz_tra.r_y)* xyz_tra.uv_y + 2.0* xyz_tra.v_z* xyz_tra.u_z + 2.0*( xyz_tra.r_z)* xyz_tra.uv_z; /*! */ /* Gradient vector length gradl. */ /* !*/ gradl= SQRT(dobjdu*dobjdu + dobjdv*dobjdv); objf_crit = SQRT(objf); #ifdef DEBUG if ( dbglev(SURPAC) == 1 && objf_crit < 5.0 ) { fprintf(dbgfil(SURPAC), "sur902*lintsur Distance (objf_crit) %15.10f u_pat %8.4f v_pat %8.4f \n" , objf_crit , u_pat , v_pat ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*lintsur objf %f dobjdu %f dobjdv %f\n" , objf , dobjdu, dobjdv ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*lintsur objf %10.1f gradl %f u_pat %8.4f v_pat %8.4f \n" , objf , gradl , u_pat , v_pat ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*lintsur d2objdu2 %f d2objdv2 %f d2objdudv %f \n" , d2objdu2, d2objdv2 , d2objdudv ); } fflush(dbgfil(SURPAC)); #endif return(SUCCED); } /* End of function */ /*!****************************************************************!*/ /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** parab3p *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function interpolates a start point, a mid point and */ /* an end point with a parabola and calculates hmin for the */ /* minimum point. */ /* */ static short parab3p() /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In and out: Global (static) variables: */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat a0,a1,a2; /* Parabola coefficients */ DBfloat hmin; /* Minimum for X=hmin */ DBfloat minm; /* Minimum minm = a1*hmin+a2*hmin**2 */ /*-----------------------------------------------------------------!*/ char errbuf[80]; /* String for error message fctn erpush */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*parab3p obj_1 %f obj_2 %f obj_3 %f hstep %f\n", obj_1 , obj_2 , obj_3 , hstep ); } #endif /*! */ /* 1. Coefficients for the parabola */ /* -------------------------------- */ /* Parabola interpolation of the three points */ /* Y = a0 + a1*X + a2*X**2 (parabola) */ /* !*/ a0 = obj_1; a2 = (obj_1/2.0 - obj_2 + obj_3/2.0)/hstep/hstep; a1 = (obj_2 - obj_1)/hstep - a2*hstep; /* Minimum point */ if ( ABS(a2) > c_denom ) { hmin = - a1/2.0/a2; minm = a0 + a1*hmin + a2*hmin*hmin; } else { hmin = hstep; minm = obj_2; } if ( minm <= obj_1 && minm <= obj_2 && minm <= obj_3 ) { /* Go (back) to hmin */ u_pat = u_pat - (2.0*hstep - hmin)*dmdu; v_pat = v_pat - (2.0*hstep - hmin)*dmdv; } else { sprintf(errbuf,"(minm)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 && n_iter > maxiter-5 ) { fprintf(dbgfil(SURPAC), "sur902*parab3p minm %f hmin %f u_pat %f v_pat %f\n", minm , hmin , u_pat , v_pat ); } fflush(dbgfil(SURPAC)); #endif return(SUCCED); } /* End of function */ /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** addsol ********************/ /* */ /* Purpose */ /* -------- */ /* */ /* Add a solution if it is a new solution (not already existing). */ /* */ /* Note that the identical points criterion (ctol) is used for */ /* all types of interrogation (c_zero is not used). */ /* */ static short addsol(p_sur,p_pat,xyz_a) /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBSurf *p_sur; /* Surface (ptr) */ DBPatch *p_pat; /* Alloc. area for topol. patch data (ptr) */ EVALS xyz_a[SMAX]; /* Array with all solutions (R*3 and UV) */ /* Global (static) variables: */ /* Points and derivatives xyz from varkon_sur_eval */ /* ..... */ /* Out: Function */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBint i_sol; /* Loop index for previous solutions */ DBint icase_l; /* Calculation case for varkon_sur_eval */ /*-----------------------------------------------------------------!*/ DBfloat dist; /* Distance between solution points (R*3) */ DBint status; /* Error code from a called function */ char errbuf[80]; /* String for error message fctn erpush */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*addsol ctol %f numsol %d \n" , ctol , numsol ); } #endif /*! */ /* 1. Check if solution already exists */ /* __________________________________ */ /* If numsol = 0 add solution (goto adds). */ /* !*/ if ( numsol == 0 ) goto adds; /*! */ /* Start loop previous solutions i_sol= 1,2, .. ,numsol. */ /* !*/ for (i_sol = 1; i_sol <= numsol ; i_sol = i_sol+1) { /*! */ /* Solution exists if the R*3 distance between xyz and xyz_a */ /* is less than ctol. Goto sex if this is the case. */ /* !*/ dist = SQRT ( (xyz.r_x-xyz_a[i_sol-1].r_x)*(xyz.r_x-xyz_a[i_sol-1].r_x) + (xyz.r_y-xyz_a[i_sol-1].r_y)*(xyz.r_y-xyz_a[i_sol-1].r_y) + (xyz.r_z-xyz_a[i_sol-1].r_z)*(xyz.r_z-xyz_a[i_sol-1].r_z) ); #ifdef DEBUG /* 1995-08-27 if ( dbglev(SURPAC) == 1 && dist > 10*ctol ) */ if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*addsol u %10.5f v %10.5f dist to prev. solut. pt %25.10f\n", xyz_a[i_sol-1].u , xyz_a[i_sol-1].v , dist ); } #endif if ( dist <= 10.0*ctol ) goto sex; } /* End loop i_sol= 1,2,..., numsol */ /*! */ /* 2. Add solution */ /* _______________ */ /* */ adds: /* Label: Add solution */ /* */ /* Increase number of solutions. */ numsol = numsol + 1; /* Check that numsol is less or equal SMAX */ /* !*/ if ( numsol > SMAX ) { sprintf(errbuf,"(SMAX)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* Recompute xyz with icase_l= 3 for xxx??? */ /* !*/ if ( ctype == SURPLAN ) { icase_l = 3; status=varkon_sur_eval (p_sur,p_pat,icase_l,xyz.u,xyz.v,&xyz); if (status<0) { sprintf(errbuf,"varkon_sur_eval%%varkon_sur_num2 (addsol)"); return(varkon_erpush("SU2943",errbuf)); } } /*! */ /* Points and coordinates to xyz_a. */ /* !*/ xyz_a[numsol-1]=xyz; /* The R*3 solution point */ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*addsol numsol %d u %8.5f v %8.5f n_iter %d objf_crit %f\n", numsol , xyz_a[numsol-1].u , xyz_a[numsol-1].v , n_iter, objf_crit ); } fflush(dbgfil(SURPAC)); #endif sex: /*! Label sex: Solution already exists !*/ return(SUCCED); } /* End of function */ /*!****************************************************************!*/ /*!New-Page--------------------------------------------------------!*/ /*!********* Internal ** function ** retrieve ***********************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function retrieves input data to local static variables. */ static short retrieve(p_spt,p_comp) /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ SPOINT *p_spt; /* Current record in point table (ptr) */ IRUNON *p_comp; /* Computation data (ptr) */ /* Out: */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* 1. Computation data. */ /* ==================== */ /* */ /* Case of computation, method and number of restarts, etc */ /* Retrieve ctype, method, nstart and maxiter from p_comp. */ /* !*/ ctype = p_comp->pro_un.ctype; method = p_comp->pro_un.method; nstart = p_comp->pro_un.nstart; maxiter = p_comp->pro_un.maxiter; pat_tra = p_comp->pro_un.pat_tra; /* 1995-09-16 */ /*! */ /* Let istart= -1 (start value for number of restarts) */ /* Let hstep = -c_infi (start value for search step) */ /* Let numsol= 0 (start value for number of solutions) */ /* !*/ c_infi = 50000.0; istart = -1; hstep = -c_infi; numsol = 0; /*! */ /* Computer tolerance. Call of varkon_comptol (sur753). */ /* !*/ comptol=varkon_comptol(); c_denom = comptol*1000.0; /*! */ /* 2. Point data. */ /* ============== */ /* Retrieve extpt for ctype= CLOSEPT from SPOINT record p_spt */ /* Let c_zero= ntol and icase= 3 for ctype= CLOSEPT */ /* Retrieve extpt, proj for ctype= LINTSUR from SPOINT record p_spt */ /* Let c_zero= ctol and icase= 3 for ctype= LINTSUR */ /* !*/ if ( ctype == CLOSEPT ) { extpt.x_gm = p_spt->extpt.x_gm ; extpt.y_gm = p_spt->extpt.y_gm ; extpt.z_gm = p_spt->extpt.z_gm ; v_proj.x_gm = 0.0; v_proj.y_gm = 0.0; v_proj.z_gm = 0.0; icase = 3; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ c_zero = p_comp->pro_un.ntol; c_zero2 = p_comp->pro_un.ctol; hstep_r3min = p_comp->pro_un.ctol/100.0; ctol = p_comp->pro_un.ctol/10.0; /* Previous value (/10)... */ } else if ( ctype == LINTSUR ) { extpt.x_gm = p_spt->extpt.x_gm ; extpt.y_gm = p_spt->extpt.y_gm ; extpt.z_gm = p_spt->extpt.z_gm ; v_proj.x_gm = p_spt->proj.x_gm ; v_proj.y_gm = p_spt->proj.y_gm ; v_proj.z_gm = p_spt->proj.z_gm ; icase = 3; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ c_zero = p_comp->pro_un.ctol; c_zero2 = c_zero; hstep_r3min = p_comp->pro_un.ctol/100.0; ctol = p_comp->pro_un.ctol/10.0; /* Previous value ... (/10) */ } /*! */ /* Retrieve search area us,vs,ue,ve from SPOINT record p_spt */ /* !*/ us = p_spt->us; vs = p_spt->vs; ue = p_spt->ue; ve = p_spt->ve; /*! */ /* If the start UV point is defined: */ /* i. Retrieve s_uvpt from SPOINT record p_spt */ /* ii. Let flag uvstart= 1 */ /* If the start UV point not is defined: */ /* i. Let flag uvstart= -1 */ /* !*/ if ( p_spt->startuvpt.z_gm >= 0.0 ) { s_uvpt = p_spt->startuvpt; uvstart = 1; } else { uvstart = -1; } #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*retrieve ctype %d method %d nstart %d maxiter %d istart %d\n", ctype, method, nstart, maxiter , istart ); if ( ctype == CLOSEPT ) { fprintf(dbgfil(SURPAC), "sur902*retrieve Case CLOSEPT: extpt %10.2f %10.2f %10.2f\n", extpt.x_gm, extpt.y_gm, extpt.z_gm ); } if ( ctype == LINTSUR ) { fprintf(dbgfil(SURPAC), "sur902*retrieve Case LINTSUR: extpt %10.2f %10.2f %10.2f\n", extpt.x_gm, extpt.y_gm, extpt.z_gm ); } fprintf(dbgfil(SURPAC), "sur902*retrieve s_uvpt %8.4f %8.4f %8.4f uvstart %d\n", s_uvpt.x_gm , s_uvpt.y_gm , s_uvpt.z_gm , uvstart ); fprintf(dbgfil(SURPAC), "sur902*retrieve Patch us %8.4f vs %8.4f ue %8.4f ve %8.4f\n", us , vs , ue , ve ); fprintf(dbgfil(SURPAC), "sur902*retrieve icase %d c_zero %f ctol %f hstep %f\n", icase , c_zero , ctol , hstep); } /* End dbglev */ #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*!********* Internal ** function ** newstart ***********************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calculates a new UV start value */ static short newstart() /* !*/ /********************************************************************/ /* ------------- Argument declarations -----------------------------*/ /* */ /* In: */ /* Out: */ /*----------------------------------------------------------------- */ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ char errbuf[80]; /* String for error message fctn erpush */ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* 1. Calculate restart value */ /* =========================== */ /* */ /* Start value for the calculation if a start UV point is defined. */ /* If istart = -1 and uvstart = 1: */ /* i. Let current UV point be u_pat,v_pat= startuvpt */ /* ii. Let istart = 0 */ /* !*/ if ( istart == -1 && uvstart == 1 ) { u_pat = s_uvpt.x_gm; v_pat = s_uvpt.y_gm; istart = 0; } /*! */ /* Start value for the calculation if no start UV point is defined. */ /* Else if istart = -1 and uvstart = -1: */ /* i. Let current UV point be u_pat,v_pat be patch mid point */ /* ii. Let istart = 1 */ /* !*/ else if ( istart == -1 && uvstart == -1 ) { u_pat = us + 0.5*(ue-us); v_pat = vs + 0.5*(ve-vs); istart = 1; } /*! */ /* Start value when the start UV point s_uvpt is used */ /* Else if istart = 0 */ /* i. Let current UV point be u_pat,v_pat be patch mid point */ /* ii. Let istart = 1 */ /* !*/ else if ( istart == 0 ) { u_pat = us + 0.5*(ue-us); v_pat = vs + 0.5*(ve-vs); istart = 1; } /*! */ /* Else if istart = 1 */ /* i. Let current UV point be u_pat= us and v_pat= vs */ /* ii. Let istart = 2 */ /* !*/ else if ( istart == 1 ) { u_pat = us; v_pat = vs; istart = 2; } /*! */ /* Else if istart = 2 */ /* i. Let current UV point be u_pat= us and v_pat= ve */ /* ii. Let istart = 3 */ /* !*/ else if ( istart == 2 ) { u_pat = us; v_pat = ve; istart = 3; } /*! */ /* Else if istart = 3 */ /* i. Let current UV point be u_pat= ue and v_pat= ve */ /* ii. Let istart = 4 */ /* !*/ else if ( istart == 3 ) { u_pat = ue; v_pat = ve; istart = 4; } /*! */ /* Else if istart = 4 */ /* i. Let current UV point be u_pat= ue and v_pat= vs */ /* ii. Let istart = 5 */ /* !*/ else if ( istart == 4 ) { u_pat = ue; v_pat = vs; istart = 5; } /*! */ /* Else if istart = 5 */ /* i. Let current UV point be u_pat= (ue-us)/4 and */ /* v_pat= (ve-vs)/2 */ /* ii. Let istart = 6 */ /* !*/ else if ( istart == 5 ) { u_pat = us + 0.25*(ue-us); v_pat = vs + 0.50*(ve-vs); istart = 6; } /*! */ /* Else if istart = 6 */ /* i. Let current UV point be u_pat= 3*(ue-us)/4 and */ /* v_pat= (ve-vs)/2 */ /* ii. Let istart = 7 */ /* !*/ else if ( istart == 6 ) { u_pat = us + 0.75*(ue-us); v_pat = vs + 0.50*(ve-vs); istart = 7; } /*! */ /* Else if istart = 7 */ /* i. Let current UV point be u_pat= (ue-us)/2 and */ /* v_pat= (ve-vs)/4 */ /* ii. Let istart = 8 */ /* !*/ else if ( istart == 7 ) { u_pat = us + 0.50*(ue-us); v_pat = vs + 0.25*(ve-vs); istart = 8; } /*! */ /* Else if istart = 8 */ /* i. Let current UV point be u_pat= (ue-us)/2 and */ /* v_pat= 3*(ve-vs)/4 */ /* ii. Let istart = 9 */ /* !*/ else if ( istart == 8 ) { u_pat = us + 0.50*(ue-us); v_pat = vs + 0.75*(ve-vs); istart = 9; } /*! */ /* Else : Error */ /* !*/ else { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*newstart Error istart %d uvstart %d\n", istart , uvstart ); } #endif sprintf(errbuf,"(istart)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* 2. Start value loop variables */ /* ============================= */ /* */ /* Start value for the the number of optimization iterations: */ nopt = -10; /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*newstart u_pat %f v_pat %f istart %d nopt %d hstep %18.12f\n", u_pat , v_pat , istart , nopt , hstep); } /* End dbglev */ #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** stepleng ******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calculates the step length hstep for the */ /* one-dimensional search. */ /* */ static short stepleng() /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: gradient, .. */ /* Out: Step length hstep */ /* */ /*-----------------------------------------------------------------!*/ /* ------------------- Theory --------------------------------------*/ /* */ /*----------------------------------------------------------------- */ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat ustep; /* Step in the U direction */ DBfloat vstep; /* Step in the V direction */ DBfloat gradm; /* Length of (dmdu,dmdv) search direction */ /*-----------------------------------------------------------------!*/ #ifdef DEBUG char errbuf[80]; /* String for error message fctn erpush */ #endif /*------------end-of-declarations-----------------------------------*/ #ifdef DEBUG if ( dbglev(SURPAC) == 2 && fabs(hstep) < 10.0*stepmin() ) { fprintf(dbgfil(SURPAC), "sur902*stepleng Input hstep %18.15f\n" , hstep ); fflush(dbgfil(SURPAC)); } #endif /*! */ /* Algorithm */ /* ========= */ /* */ /* 1. Step length */ /* ============== */ /* */ /* Search vector length gradm */ /* !*/ gradm= SQRT(dmdu*dmdu + dmdv*dmdv); if ( gradm < c_denom ) goto noadj; /*! */ /* Maximum step in U direction ustep (minimum 3 steps) */ /* !*/ if ( fabs(dmdu) < c_denom ) { ustep = c_infi; } else if ( dmdu > 0.0 ) { /* ustep = (u_pat-us)/3.0/dmdu; 1997-02-27 */ ustep = (u_pat-us-comptol)/3.0/dmdu; if ( ustep < 0 ) ustep = c_infi; } else if ( dmdu < 0.0 ) { /* ustep = -(ue-u_pat)/3.0/dmdu; 1997-02-27 */ ustep = -(ue+comptol-u_pat)/3.0/dmdu; if ( ustep < 0 ) ustep = c_infi; } /*! */ /* Maximum step in V direction vstep (minimum 3 steps) */ /* !*/ if ( fabs(dmdv) < c_denom ) { vstep = c_infi; } else if ( dmdv > 0.0 ) { /* vstep = (v_pat-vs)/3.0/dmdv; 1997-02-27 */ vstep = (v_pat-vs-comptol)/3.0/dmdv; if ( vstep < 0 ) vstep = c_infi; } else if ( dmdv < 0.0 ) { /* vstep = -(ve-v_pat)/3.0/dmdv; 1997-02-27 */ vstep = -(ve+comptol-v_pat)/3.0/dmdv; if ( vstep < 0 ) vstep = c_infi; } /*! */ /* Choose the minimum value of (input) hstep, ustep and vstep. */ /* If hstep has the minimum value, modify hstep to exactly */ /* "hit" a patch border. */ /* !*/ #ifdef DEBUG if(ustep < 0.0 || vstep < 0.0 || hstep > 0.0 ) { sprintf(errbuf,"(step sign)%% varkon_sur_num2 (sur902)"); return(varkon_erpush("SU2993",errbuf)); } #endif if ( ustep < -hstep && ustep < vstep ) { hstep = -ustep; } else if ( vstep < -hstep && vstep < ustep ) { hstep = -vstep; } else /* hstep has the minimum value */ { if ( ustep < vstep ) { while ( ustep > -hstep ) ustep= ustep/2.0; hstep = -ustep; } else { while ( vstep > -hstep ) vstep= vstep/2.0; hstep = -vstep; } } #ifdef DEBUG if ( dbglev(SURPAC) == 1 && objf_crit < 1000.0*c_zero ) { fprintf(dbgfil(SURPAC), "sur902*stepleng ustep %12.9f vstep %12.9f hstep %12.9f (out)\n", ustep,vstep,hstep ); fflush(dbgfil(SURPAC)); } #endif /*! */ /*! */ /* If the first point is outside us,ue,vs,ve is the current */ /* point on the border and search will be outside area. */ /* Let hstep= 0 in this case. */ /* !*/ if ( u_pat + hstep*dmdu >= us && u_pat + hstep*dmdu <= ue && v_pat + hstep*dmdv >= vs && v_pat + hstep*dmdv <= ve ) { ; } else { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*stepleng Outside u_pat %f v_pat %f\n" , u_pat ,v_pat ); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*stepleng Outside u_pat+hstep*dmdu %f v_pat+hstep*dmdv %f\n" , u_pat+hstep*dmdu,v_pat+hstep*dmdv ); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*stepleng ustep %f vstep %f hstep %f \n" , ustep, vstep, hstep ); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*stepleng dmdu %f dmdv %f gradm %f \n" , dmdu , dmdv , gradm ); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur902*stepleng us %6.3f ue %6.3f vs %6.3f ve %6.3f\n" , us, ue, vs, ve ); } #endif hstep = 0.0; } noadj: /* No adjustment of the step length (gradm= 0) */ #ifdef DEBUG if ( dbglev(SURPAC) == 2 && (n_iter > maxiter-5 || fabs(hstep) < 10.0*stepmin() )) { fprintf(dbgfil(SURPAC), "sur902*stepleng ue-u_pat %f ve-v_pat %f\n" , ue-u_pat,ve-v_pat ); } if ( dbglev(SURPAC) == 2 && (n_iter > maxiter-5 || fabs(hstep) < 10.0*stepmin() )) { fprintf(dbgfil(SURPAC), "sur902*stepleng ustep %f vstep %f hstep %f \n" , ustep, vstep, hstep ); } if ( dbglev(SURPAC) == 2 && (n_iter > maxiter-5 || fabs(hstep) < 10.0*stepmin() )) { fprintf(dbgfil(SURPAC), "sur902*stepleng dmdu %f dmdv %f gradm %f hstep %18.15f\n" , dmdu , dmdv , gradm, hstep ); fflush(dbgfil(SURPAC)); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** sdirect *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calculates the direction for the */ /* one-dimensional search (dmdu,dmdv). */ /* */ static short sdirect () /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: gradient dmdu,dmdv */ /* Out: The search direction dmdu,dmdv */ /* */ /*-----------------------------------------------------------------!*/ /* ------------------- Theory --------------------------------------*/ /* */ /*----------------------------------------------------------------- */ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat gradm; /* Length of (dmdu,dmdv) search direction */ /*-----------------------------------------------------------------!*/ DBfloat ellip; /* For the analysis of the Hessian matrix */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* 1. Hessian adjustment */ /* ====================== */ /* */ /* If .. ? .. adjustment of the search direction */ /* */ /* Calculation of the Hessian determinant deterh */ /* and the new search direction (dmdu,dmdv). */ /* !*/ deterh = d2objdu2*d2objdv2 - d2objdudv*d2objdudv; /* Multiplication with the inverse (with deterh ! ) */ dmdu = d2objdv2*dobjdu - d2objdudv*dobjdv; dmdv = - d2objdudv*dobjdu + d2objdu2*dobjdv; if( ABS(deterh) > c_denom ) { dmdu = dmdu/deterh; dmdv = dmdv/deterh; } /*! */ /* Analysis of the Hessian determinant deterh */ /* Let hstep= -0.8*gradm and goto _ellip for an elliptic point. */ /* (Note that the length of (dmdu,dmdv) is to the minimum pt) . */ /* Goto _saddl for a saddle point. */ /* !*/ /* For the analysation of an elliptic point */ ellip = d2objdu2*u_pat*u_pat + d2objdv2*v_pat*v_pat + 2.0*d2objdudv*u_pat*v_pat; if ( ellip > 0.0 && deterh > c_denom ) { /* Start elliptic point */ gradm= SQRT(dmdu*dmdu + dmdv*dmdv); hstep = -0.80*gradm; goto _ellip; } if ( deterh < -c_denom ) { #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*sdirect Saddle pt deterh %f u_pat %8.4f v_pat %8.4f \n" , deterh , u_pat, v_pat ); } #endif goto _saddl; } /*! */ /* 3. Search direction equal to gradient */ /* ===================================== */ /* */ _saddl: /* */ /* Let (dmdu,dmdv)= (dobjdu,dobjdv). */ /* !*/ dmdu = dobjdu; dmdv = dobjdv; /*! */ /* 4. Normalise search direction vector */ /* ===================================== */ /* */ _ellip: /* */ /* Let (dmdu,dmdv)= (dmdu/gradm,dmdv/gradm) */ /* !*/ gradm= SQRT(dmdu*dmdu + dmdv*dmdv); if ( gradm > c_denom ) { dmdu = dmdu/gradm; dmdv = dmdv/gradm; } #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*sdirect deterh %f \n" , deterh ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*sdirect deterh %f ellip %f dmdu %8.4f dmdv %8.4f \n" , deterh , ellip , dmdu, dmdv ); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*!************** Internal ** function ** stepmin ******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calculates the U,V minimum step length for hstep */ /* */ /* Note that this internal function is DBfloat (and not short)! */ /* */ /* */ static DBfloat stepmin () /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: hstep_r3min, xyz and dmdu,dmdv */ /* Out: The function value is the minimum step */ /* */ /*-----------------------------------------------------------------!*/ /*!------------------- Theory --------------------------------------*/ /* */ /* Reference: Faux & Pratt p 107-108 */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat hstep_uvmin; /* Minimum U,V length for hstep */ DBfloat tang_x; /* Search direction tangent, X component */ DBfloat tang_y; /* Search direction tangent, Y component */ DBfloat tang_z; /* Search direction tangent, Z component */ DBfloat tang_length; /* Search direction tangent length */ /*-----------------------------------------------------------------!*/ DBfloat dmdu_n; /* dmdu normalized U component */ DBfloat dmdv_n; /* dmdv normalized V component */ DBfloat v_length; /* Length of a vector */ /*------------end-of-declarations-----------------------------------*/ #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur902*stepmin hstep_r3min %18.15f\n" , hstep_r3min ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif /*! */ /* Algorithm */ /* ========= */ /* !*/ /* 0. Initializations */ hstep_uvmin = 12345.6789; dmdu_n = 0.0; dmdv_n = 0.0; tang_x = 0.0; tang_y = 0.0; tang_z = 0.0; tang_length = -0.123456789; /*! */ /* 1. Tangent in the search direction */ /* */ /* !*/ v_length= SQRT(dmdu*dmdu + dmdv*dmdv); if ( v_length < c_denom ) { hstep_uvmin = c_infi; goto undef; } dmdu_n = dmdu/v_length; dmdv_n = dmdv/v_length; tang_x = xyz.u_x*dmdu_n + xyz.v_x*dmdv_n; tang_y = xyz.u_y*dmdu_n + xyz.v_y*dmdv_n; tang_z = xyz.u_z*dmdu_n + xyz.v_z*dmdv_n; /*! */ /* 2. Minimum U,V step = Minimum R*3 step / Search tangent length */ /* */ /* !*/ tang_length = SQRT(tang_x*tang_x + tang_y*tang_y + tang_z*tang_z); if ( tang_length < c_denom ) { hstep_uvmin = c_infi; goto undef; } hstep_uvmin = hstep_r3min/tang_length; undef: /* Surface or search direction is undefined */ #ifdef DEBUG if ( dbglev(SURPAC) == 1 && hstep_uvmin > 100.0 ) { fprintf(dbgfil(SURPAC), "sur902*stepmin hstep_r3min %18.15f\n" , hstep_r3min ); fprintf(dbgfil(SURPAC), "sur902*stepmin dmdu %25.15f dmdv %25.15f\n" , dmdu , dmdv ); fprintf(dbgfil(SURPAC), "sur902*stepmin dmdu_n %25.15f dmdv_n %25.15f\n" , dmdu_n , dmdv_n ); fprintf(dbgfil(SURPAC), "sur902*stepmin dr/du %10.5f %10.5f %10.5f\n", xyz.u_x, xyz.u_y, xyz.u_z ); fprintf(dbgfil(SURPAC), "sur902*stepmin dr/dv %10.5f %10.5f %10.5f\n", xyz.v_x, xyz.v_y, xyz.v_z ); fprintf(dbgfil(SURPAC), "sur902*stepmin Search direction tangent %10.5f %10.5f %10.5f\n", tang_x,tang_y,tang_z ); fprintf(dbgfil(SURPAC), "sur902*stepmin Search direction tangent lenght %25.15f \n" , tang_length ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 && hstep_uvmin > 100.0 ) { fprintf(dbgfil(SURPAC), "sur902*stepmin hstep %25.15f hstep_uvmin %25.15f\n" , hstep, hstep_uvmin ); fflush(dbgfil(SURPAC)); } #endif return(hstep_uvmin); } /* End of function */ /********************************************************************/