/********************************************************************/ /* */ /* 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_num1 File: sur900.c */ /* ============================================================= */ /* */ /* Purpose */ /* ------- */ /* */ /* The function is a general function for the calculation */ /* of points on an iso-parametric curve on a surface. */ /* */ /* 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 methods available are Newton-Rhapson, parabola inter- */ /* polation, .... */ /* */ /* Author: Gunnar Liden */ /* !*/ /* Revisions */ /* */ /* 1994-06-03 Originally written */ /* 1996-02-23 icase= 4 for I_SILHOU and F_SILHOU, unused variabl */ /* 1996-04-29 analyze added */ /* 1996-05-02 analyze trimming for silhouettes */ /* 1998-09-26 analyze Tolerance idpoint for intersects, */ /* TOL2->idpoint, ana2pts added */ /* 1999-11-24 Free source code modifications */ /* */ /********************************************************************/ /* ------------- Short description of function -----------------*/ /* */ /*sdescr varkon_sur_num1 Surface one param. calculation */ /* */ /*------------------------------------------------------------- */ /*!-------------- Function calls (internal) ------------------------*/ /* */ #ifdef DEBUG static short initial(); /* Initialization of variables */ #endif static short fsilhou(); /* Function f_c and dfdu for F_SILH */ static short isilhou(); /* Function f_c and dfdu for I_SILH */ static short ssilhou(); /* Function f_c and dfdu for S_SILH */ static short intplan(); /* Function f_c and dfdu for SURPLAN */ static short addsol(); /* Add a non-existent solution. */ static short analyze(); /* Analyze boundary */ static short ana2pts(); /* Additional comparison of solutions */ /* */ /*-----------------------------------------------------------------!*/ #ifdef NOT_YET_IMPLEMENTED static short newton1(); /* Newton Rhapson method */ static short fisopho(); /* Function f_c and dfdu for F_ISOP */ static short iisopho(); /* Function f_c and dfdu for I_ISOP */ #endif /* NOT_YET_IMPLEMENTED */ #ifdef NOT_USED_FOR_THE_MOMENT static DBfloat tglobal_pre; /* Previous parameter value */ static short parab2d(); /* Parabola interp.: 2 pts+derivative */ #endif /* NOT_USED_FOR_THE_MOMENT */ /*!-- Static (common) variables for the functions in this file -----*/ /* */ static DBint ctype; /* Case of computation */ static DBfloat c_zero; /* The zero (0.0) end criterion */ static DBfloat idpoint; /* Identical points criterion */ static DBfloat i_zero; /* Criterion for isoparameter */ static DBint method; /* Method which shall be used */ static DBfloat t_s,t_e; /* Param. interv. for calcul. */ static DBfloat t_c; /* Constant parameter value */ static DBint iso_p; /* Eq. 1: Const U Eq. 2: Const V */ static DBfloat tglobal; /* Current parameter value */ static DBfloat f_cur; /* Current (objective) fctn value */ static DBfloat f_pre; /* Previous function value */ static DBfloat dfdt; /* Function derivative w.r.t t */ static DBint numsol; /* Number of solutions (=*p_npts) */ static DBfloat ctol; /* Coordinate tolerance */ static DBfloat ntol; /* Surface normal tolerance */ static DBfloat comptol; /* Computer tolerance */ static DBfloat pi_x,pi_y,pi_z,pi_d;/* Intersecting plane */ static DBfloat eye_x,eye_y,eye_z; /* Eye point for reflection crvs */ static DBfloat eyev_x,eyev_y,eyev_z;/* Eye vect. for reflection crvs */ static DBVector sdir; /* Direction for spine silhouette */ static DBfloat angle; /* Angle in degrees for isophotes */ static EVALS xyz; /* Coordinates and derivatives */ /* for a point on a surface */ /*-----------------------------------------------------------------!*/ static DBint i_rg; /* Temporary !!!! */ /*!-------------- Function calls (external) ------------------------*/ /* */ /* varkon_comptol * Retrieve computer tolerance */ /* varkon_sur_eval * Surface evaluation routine */ /* varkon_vec_projpla * Project vector onto plane */ /* varkon_ini_evals * Initiate EVALS */ /* varkon_erpush * Error message to terminal */ /* */ /*-----------------------------------------------------------------!*/ /*!------------ Error messages and warnings ------------------------*/ /* */ /* SU2943 = Called function xxxxxx failed in varkon_sur_num1 */ /* SU2983 = sur900 Illegal computation case= for varkon_sur_num1 */ /* SU2973 = Internal function xxxxxx failed in varkon_sur_num1 */ /* SU2993 = Severe program error in varkon_sur_num1 (sur900). */ /* */ /*-----------------------------------------------------------------!*/ /*!****************** Function **************************************/ /* */ DBstatus varkon_sur_num1 ( /*-------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBSurf *p_sur, /* Surface (ptr) */ DBPatch *p_pat, /* Alloc. area for topol. patch data (ptr) */ DBfloat t_in[2][2], /* Iso-parameter curve limits */ /* (0,0): Start U (0,1): Start V */ /* (1,0): End U (1,1): End V */ IRUNON *p_comp, /* Computation data (ptr) */ /* Out: */ DBint *p_npts, /* Number of output points (ptr) */ EVALS xyz_a[SMAX], /* Array with all solutions (R*3 and UV) */ DBint *p_b_case ) /* Boundary point case: */ /* Eq. 0: No points */ /* Eq. 1: Normal points */ /* Eq. 2: Curve coincides w. isoparameter */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBint nstart; /* The number of restarts */ DBint maxiter; /* Maximum number of iterations */ DBfloat deltat; /* The parameter t step for the restarts */ DBint i_r; /* Loop index for the restarts */ DBint j_n; /* Loop index for the numerical solution */ DBfloat h; /* Step determined by Newton !!! Temp. !! */ DBint icase; /* Calculation case for varkon_sur_eval */ /* */ /*-----------------------------------------------------------------!*/ DBint f_a; /* Flag for analyzis */ /* Eq. 0: No point (numsol= 0) */ /* Eq. 1: One end point (numsol= 0) */ /* Eq. 2: Two points (numsol= 2) */ char errbuf[80]; /* String for error message fctn erpush */ DBfloat dum1; /* Dummy float */ DBint status; /* Error code from a called function */ #ifdef DEBUG DBint ix1; /* Temporarely used loop index */ #endif /*--------------end-of-declarations---------------------------------*/ /*!New-Page--------------------------------------------------------!*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*Enter varkon_sur_num1: us %f vs %f ue %f ve %f\n", t_in[0][0],t_in[0][1],t_in[1][0],t_in[1][1]); } #endif /*! */ /* 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. */ /* !*/ /* Printout of input data is in function initial */ icase = I_UNDEF; #ifdef DEBUG status=initial(p_sur,p_pat,t_in,p_comp,p_npts,xyz_a); if (status<0) { sprintf(errbuf,"initial%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } #endif comptol = varkon_comptol(); *p_b_case = I_UNDEF; /*! */ /* Case of computation, method, idpoint and number of restarts */ ctype = p_comp->ipl_un.ctype; method = p_comp->ipl_un.method; nstart = p_comp->ipl_un.nstart; maxiter = p_comp->ipl_un.maxiter; /* 1998-09-26 TODO TOL2= 0.001 and idpoint= 0.01 in SURPAC.TOL */ /* The default tolerance idpoint should be changed to 0.001 */ /* and /10.0 should be deleted. */ idpoint = p_comp->ipl_un.idpoint/10.0; /* Determine with t_in if U (iso_p=1) or V (iso_p=2) is constant */ /* !*/ if ( fabs(t_in[0][0]-t_in[1][0]) < TOL4 ) iso_p = 1; else if ( fabs(t_in[0][1]-t_in[1][1]) < TOL4 ) iso_p = 2; else { sprintf(errbuf,"(non-isop)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* Retrieve computation data for the input computation case */ /* (for instance the intersection plane for ctype= SURPLAN) */ /* and definition of the end criterion for c_zero */ /* !*/ if ( ctype == F_SILH ) { eye_x = p_comp->sil_un.eye_x ; eye_y = p_comp->sil_un.eye_y ; eye_z = p_comp->sil_un.eye_z ; icase = 4; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ c_zero = p_comp->sil_un.ntol; i_zero = 50.0*c_zero; /* TODO use idpoint, but test first */ } else if ( ctype == I_SILH ) { eyev_x = p_comp->sil_un.eyev_x ; eyev_y = p_comp->sil_un.eyev_y ; eyev_z = p_comp->sil_un.eyev_z ; icase = 4; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ ntol = p_comp->sil_un.ntol; c_zero = ntol; i_zero = 50.0*c_zero; /* TODO use idpoint, but test first */ } else if ( ctype == S_SILH ) { sdir.x_gm = p_comp->sil_un.eyev_x ; sdir.y_gm = p_comp->sil_un.eyev_y ; sdir.z_gm = p_comp->sil_un.eyev_z ; icase = 2; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ ntol = p_comp->sil_un.ntol; c_zero = ntol; i_zero = 50.0*c_zero; /* TODO use idangle, but test first */ } else if ( ctype == F_ISOP ) { eye_x = p_comp->iso_un.eye_x ; eye_y = p_comp->iso_un.eye_y ; eye_z = p_comp->iso_un.eye_z ; angle = p_comp->iso_un.angle ; icase = 4; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ ntol = p_comp->sil_un.ntol; c_zero = ntol; i_zero = 50.0*c_zero; /* TODO use idangle, but test first */ } else if ( ctype == I_ISOP ) { eye_x = p_comp->iso_un.eye_x ; eye_y = p_comp->iso_un.eye_y ; eye_z = p_comp->iso_un.eye_z ; angle = p_comp->iso_un.angle ; icase = 4; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ ntol = p_comp->sil_un.ntol; c_zero = ntol; i_zero = 50.0*c_zero; /* TODO use idangle, but test first */ } else if ( ctype == SURPLAN ) { pi_x = p_comp->ipl_un.in_x; pi_y = p_comp->ipl_un.in_y; pi_z = p_comp->ipl_un.in_z; pi_d = p_comp->ipl_un.in_d; if ( iso_p == 1 ) icase = 2; if ( iso_p == 2 ) icase = 1; ctol = p_comp->ipl_un.ctol; /* End criterion for the numerical solution ( !!! geo995 !!! ) */ c_zero = ctol/10.0; i_zero = p_comp->ipl_un.idpoint; } else { sprintf(errbuf,"(ctype )%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* Parameter interval t_s and t_e for the numerical solution */ /* !*/ if ( iso_p == 1) { t_s = t_in[0][1]; t_e = t_in[1][1]; t_c = t_in[0][0]; } else if ( iso_p == 2) { t_s = t_in[0][0]; t_e = t_in[1][0]; t_c = t_in[0][1]; } if ( t_s > t_e ) { dum1 = t_s; t_s = t_e; t_e = dum1; } /*! */ /* The parameter step for the restarts */ deltat = (t_e-t_s-TOL4)/((DBfloat)nstart-1.0); /* The number of output points */ numsol = 0; /* !*/ /*! */ /* Analyze if whole, or part of boundary is the solution */ /* Call of internal function analyze. */ /* !*/ analyze (p_sur,p_pat,nstart,icase,xyz_a,&f_a); if ( f_a >= 1 ) goto boundsol; /*! */ /* 2. Numerical solution for a one parameter function */ /* __________________________________________________ */ /* */ /* Start loop: Restarts i_r=1,2, .. nstart */ /* !*/ i_rg = 0; /* Temporary !!!!! */ for (i_r = 1; i_r <= nstart ; i_r = i_r+1) { /*! */ /* New parameter start value */ tglobal = t_s + deltat*((DBfloat)i_r-1.0); /* Start loop: Numerical i_n=1,2, .. */ /* !*/ for(j_n = 1; j_n <= maxiter ; j_n = j_n+1) { /*! */ /* Call of varkon_sur_eval (sur210) for all surface types */ /* (Error SU2943 if called function fails) */ /* !*/ i_rg = i_rg+1; /* Temporary !!!!! */ if (iso_p == 2) status= varkon_sur_eval (p_sur,p_pat,icase,tglobal,t_c,&xyz); else status= varkon_sur_eval (p_sur,p_pat,icase,t_c,tglobal,&xyz); if (status<0) { sprintf(errbuf,"varkon_sur_eval%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2943",errbuf)); } /*! */ /* Calculation of function f_cur and derivative dfdu. */ /* Call of fsilhou for ctype = F_SILH */ /* Call of isilhou for ctype = I_SILH */ /* Call of ssilhou for ctype = S_SILH */ /* Call of fisopho for ctype = F_ISOP */ /* Call of iisopho for ctype = I_ISOP */ /* Call of intplan for ctype = SURPLAN */ /* (Error SU2973 if internal function fails) */ /* !*/ f_pre = f_cur; /* Temporary !!! */ if ( ctype == F_SILH ) status=fsilhou(); if (status<0) { sprintf(errbuf,"fsilhou%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } else if ( ctype == I_SILH ) status=isilhou(); if (status<0) { sprintf(errbuf,"isilhou%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } else if ( ctype == S_SILH ) status=ssilhou(); if (status<0) { sprintf(errbuf,"ssilhou%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } else if ( ctype == F_ISOP ) status=intplan(); if (status<0) { sprintf(errbuf,"fisopho%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } else if ( ctype == I_ISOP ) status=intplan(); if (status<0) { sprintf(errbuf,"iisopho%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } else if ( ctype == SURPLAN ) status=intplan(); if (status<0) { sprintf(errbuf,"intplan%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } /* !!!! ?????? else { sprintf(errbuf,"(ctype 2)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } !!!! ???? */ /*! */ /* Is function value f_cur <= c_zero ? */ /* !*/ if( fabs(f_cur) <= c_zero ) { /*! */ /* Yes, f_cur= 0 for curve parameter value tglobal: */ /* Increase the number of solutions numsol and add */ /* the solution (R*3 and U,V) to output array xyz_a */ /* if the parameter tglobal is between t_s and t_e */ /* and if the solution not exists (call of addsol). */ /* !*/ /* ?????? Note ?: tglobal= t_s and t_e-TOL4 will always be */ /* since the minimum number of restarts is 2 */ /* TODO !!!! Problem to have the solution as a u val */ /* Accepting +/- TOL4 will give the wrong segment */ /* adress */ if ( tglobal >= t_s ) /* Parameter is greater or */ { /* equal interval start t_s */ if ( tglobal < t_e ) /* Parameter value is less */ { /* than interval end t_e */ /* 1998-09-26 addsol(p_sur,p_pat,xyz_a);*//* Add non-existent solution */ /* Add non-existent solution */ status = addsol(p_sur,p_pat,xyz_a); if ( status < 0 ) return(status); } /* End tglobal < t_e */ } /* End tglobal >= t_s */ j_n = maxiter; /* Refer to Note 2 above */ } /* End if f_cur <= c_zero */ /*! */ /* Is derivative dfdu too small ? */ /* !*/ if( fabs(dfdt) <= TOL3 ) /* Derivative <= TOL3 ? */ { /* */ /*! */ /* Yes, dfdu is less or equal TOL3: */ j_n = maxiter; /* ==> end of loop j_n */ h = TOL6; /* ??? Newton only ??? */ /* !*/ } /* End dfdu<= TOL3 */ else /* Derivative > TOL3 */ { /* */ h = -f_cur/dfdt; /* Newton step ... */ } /* End if: dfdu <= TOL3 */ /* Solution outside 0 < u < 1 ? */ if (fabs(tglobal-t_s) < TOL4 ) { if ( h < 0 ) j_n = maxiter; } if ( fabs(tglobal-t_e) < TOL4 ) { if ( h > 0 ) j_n = maxiter; } if ( fabs(h) > 2.0 ) j_n = maxiter; /* Compute the next parameter value */ tglobal = tglobal + h; if ( tglobal < t_s - TOL4 ) tglobal = t_s + 0.1*TOL4; if ( tglobal > t_e + TOL4 ) tglobal = t_e - 0.1*TOL4; #ifdef NOT_USED_FOR_THE_MOMENT status=parab2d(); if (status<0) { sprintf(errbuf,"parab2d%%varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2973",errbuf)); } #endif /* NOT_USED_FOR_THE_MOMENT */ if ( j_n > maxiter ) { sprintf(errbuf,"(maxiter)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } } /* end loop j_n : Numerical solution */ } /* end loop i_r : Restarts with deltat */ boundsol:; /*! */ /* 3. Exit */ /* ________ */ /* */ /* Number of solutions to output variable */ *p_npts = numsol; /* !*/ if ( f_a == 2 ) *p_b_case = 2; else if ( f_a != 2 && *p_npts > 0 ) *p_b_case = 1; else *p_b_case = 0; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { if (*p_npts> 0 ) { for (ix1 = 1; ix1 <= *p_npts; ix1 = ix1+1) { fprintf(dbgfil(SURPAC), "sur900 Point %d U %9.6f V %9.6f X %9.2f Y %9.2f Z %9.2f\n", (int)ix1 , xyz_a[ix1-1].u ,xyz_a[ix1-1].v , xyz_a[ix1-1].r_x,xyz_a[ix1-1].r_y,xyz_a[ix1-1].r_z); } } } #endif return(SUCCED); } /* End of function */ /********************************************************************/ #ifdef DEBUG /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*!********* Internal ** function ** initial ************************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function initializes the output variables and the */ /* static (common) variables to 1.23456789 and 123456789. */ static short initial(p_sur,p_pat,t_in,p_comp,p_npts,xyz_a) /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBSurf *p_sur; /* Surface (ptr) */ DBPatch *p_pat; /* Alloc. area for topol. patch data (ptr) */ DBfloat t_in[2][2]; /* Iso-parameter curve limits */ /* (0,0): Start U (0,1): Start V */ /* (1,0): End U (1,1): End V */ IRUNON *p_comp; /* Computation data (ptr) */ DBint *p_npts; /* Number of output points (ptr) */ EVALS xyz_a[SMAX]; /* Array with all solutions (R*3 and UV) */ /* Out: */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ DBint iso; /* Eq. 1: Const U Eq. 2: Const V */ /*----------------------------------------------------------------- */ EVALS xyz_d; /* Temporary used variable */ DBint i_p; /* Loop index */ 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_num1 (sur900) ** initial ******\n"); fprintf(dbgfil(SURPAC), "sur900 Start t_in(0,0-1) %f %f End t_in(1,0-1) %f %f \n", t_in[0][0],t_in[0][1],t_in[1][0],t_in[1][1]); fprintf(dbgfil(SURPAC), "sur900 Case of comp. %d Method %d No restarts %d Max iter %d\n", (int)p_comp->ipl_un.ctype, (int)p_comp->ipl_un.method, (int)p_comp->ipl_un.nstart, (int)p_comp->ipl_un.maxiter); if ( p_comp->ipl_un.ctype == I_SILH) { fprintf(dbgfil(SURPAC), "sur900 Case I_SILH: Eye vector= %f %f %f \n", p_comp->sil_un.eyev_x, p_comp->sil_un.eyev_y, p_comp->sil_un.eyev_z); } if ( p_comp->ipl_un.ctype == S_SILH) { fprintf(dbgfil(SURPAC), "sur900 Case S_SILH: U crv tangent %f %f %f\n", p_comp->sil_un.eyev_x, p_comp->sil_un.eyev_y, p_comp->sil_un.eyev_z); } if ( p_comp->ipl_un.ctype == F_SILH) { fprintf(dbgfil(SURPAC), "sur900 Case F_SILH: Eye point= %f %f %f \n", p_comp->sil_un.eye_x, p_comp->sil_un.eye_y, p_comp->sil_un.eye_z); } if ( p_comp->ipl_un.ctype == SURPLAN) { fprintf(dbgfil(SURPAC), "sur900 Case SURPLAN: Intersect plane: %f %f %f %f\n", p_comp->ipl_un.in_x, p_comp->ipl_un.in_y, p_comp->ipl_un.in_z, p_comp->ipl_un.in_d); } } fflush(dbgfil(SURPAC)); /* To file from buffer */ if ( fabs(t_in[0][0]-t_in[1][0]) < TOL4 ) iso = 1; else if ( fabs(t_in[0][1]-t_in[1][1]) < TOL4 ) iso = 2; else { sprintf(errbuf,"(non-isop)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } /*! Parameter values must be in same patch Error SU2993 if not OK. !*/ if (iso == 2 ) { /* Not complete !!!! */ if ( fabs(t_in[0][0]-t_in[1][0]) > 1.0 ) { sprintf(errbuf,"t_in(0-1)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } if ( fabs(t_in[0][0]-t_in[1][0]) < TOL4 ) { sprintf(errbuf,"t_in(0-1)=0%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } } if (iso == 1 ) { /* Not complete !!!! */ if ( fabs(t_in[0][1]-t_in[1][1]) > 1.0 ) { sprintf(errbuf,"t_in(2-3)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } if ( fabs(t_in[0][1]-t_in[1][1]) < TOL4 ) { sprintf(errbuf,"t_in(2-3)=0%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } } /*! The number of restarts must be between 2 and maxiter */ /* Error SU2993 if not OK. !*/ if (p_comp->ipl_un.nstart > 1 || p_comp->ipl_un.nstart <= p_comp->ipl_un.maxiter ) { ; } else { sprintf(errbuf,"(nstart)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* 2. Initiate output variables and static (internal) variables */ /* !*/ *p_npts = -123456789; eyev_x = -1.23456789; eyev_y = -1.23456789; eyev_z = -1.23456789; eye_x = -1.23456789; eye_y = -1.23456789; eye_z = -1.23456789; angle = -1.23456789; pi_x = -1.23456789; pi_y = -1.23456789; pi_z = -1.23456789; pi_d = -1.23456789; for (i_p=1; i_p<=SMAX; i_p= i_p+1) /* Start loop */ { /* Initialize variabel EVALS */ /* Call of varkon_ini_evals (sur170). */ /* */ varkon_ini_evals (&xyz_d); xyz_a[i_p-1] = xyz_d; } /* End loop */ return(SUCCED); } /* End of function */ /********************************************************************/ #endif /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** fsilhou *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calcultes the function f_cur and dfdt */ /* for calculation case: F_SILH */ /* */ static short fsilhou() /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: */ /* Points and derivatives xyz from varkon_sur_eval */ /* The eye point eye_x, eye_y, eye_z and flag iso_p */ /* Out: Function f_cur and derivative dfdt */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* Function value: */ f_cur = xyz.n_x*xyz.r_x + xyz.n_y*xyz.r_y + xyz.n_z*xyz.r_z - eye_x*xyz.n_x - eye_y*xyz.n_y - eye_z*xyz.n_z; /* !*/ if ( iso_p == 2 ) { /*! */ /* Function derivative value dfdt for a constant V line (iso_p=2): */ dfdt = xyz.nu_x*xyz.r_x + xyz.nu_y*xyz.r_y + xyz.nu_z*xyz.r_z + xyz.n_x*xyz.u_x + xyz.n_y*xyz.u_y + xyz.n_z*xyz.u_z - xyz.nu_x*eye_x - xyz.nu_y*eye_y - xyz.nu_z*eye_z; /* !*/ } else { /*! */ /* Function derivative value dfdt for a constant U line (iso_p=1): */ dfdt = xyz.nv_x*xyz.r_x + xyz.nv_y*xyz.r_y + xyz.nv_z*xyz.r_z + xyz.n_x*xyz.v_x + xyz.n_y*xyz.v_y + xyz.n_z*xyz.v_z - xyz.nv_x*eye_x - xyz.nv_y*eye_y - xyz.nv_z*eye_z; /* !*/ } #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur900*fsilhou t_c %f tglobal %f f_cur %f dfdt %f\n" , t_c , tglobal, f_cur, dfdt ); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** isilhou *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calcultes the function f_cur and dfdt */ /* for calculation case: I_SILH */ /* */ static short isilhou() /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: */ /* Points and derivatives xyz from varkon_sur_eval */ /* The eye vector eyev_x, eyev_y, eyev_z and flag iso_p */ /* Out: Function f_cur and derivative dfdt */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* Function value: */ f_cur = xyz.n_x*eyev_x + xyz.n_y*eyev_y + xyz.n_z*eyev_z; /* !*/ if ( iso_p == 2 ) /* 1 eller 2 ????? */ { /*! */ /* Function derivative value dfdt for a constant V line (iso_p=2): */ dfdt = xyz.nu_x*eyev_x + xyz.nu_y*eyev_y + xyz.nu_z*eyev_z; /* !*/ } else { /*! */ /* Function derivative value dfdt for a constant U line (iso_p=1): */ dfdt = xyz.nv_x*eyev_x + xyz.nv_y*eyev_y + xyz.nv_z*eyev_z; /* !*/ } #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { /* fprintf(dbgfil(SURPAC), "POI_FREE(#%d,VEC(%f,%f)); ! t xyz.n_y xyz.nu_y= %f \n", (int)i_rg, tglobal, xyz.n_y, xyz.nu_y); fprintf(dbgfil(SURPAC), "POI_FREE(#%d,VEC(%f,%f)); ! t xyz.n_z xyz.nu_z= %f \n", (int)i_rg+100, tglobal, xyz.n_z, xyz.nu_z); */ fprintf(dbgfil(SURPAC), "sur900*isilhou t_c %f tglobal %f f_cur %f dfdt %f\n" , t_c , tglobal, f_cur, dfdt ); fprintf(dbgfil(SURPAC), "sur900*isilhou xyz.nu_x %f xyz.nu_y %f xyz.nu_z %f\n" , xyz.nu_x,xyz.nu_y,xyz.nu_z ); fprintf(dbgfil(SURPAC), "sur900*isilhou xyz.nv_x %f xyz.nv_y %f xyz.nv_z %f\n" , xyz.nv_x,xyz.nv_y,xyz.nv_z ); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** ssilhou *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calcultes the function f_cur and dfdt */ /* for calculation case: S_SILH */ /* */ static short ssilhou() /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: */ /* Points and derivatives xyz from varkon_sur_eval */ /* The vector sdir */ /* Out: Function f_cur and derivative dfdt */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ DBVector plane; /* Plane normal */ DBint ocase; /* Output case. Eq. 1: Normalised vector */ DBVector psdir; /* U crv tangent projected in spine plane */ /*----------------------------------------------------------------- */ DBint status; /* Error code from a called function */ char errbuf[80]; /* String for error message fctn erpush */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ plane.x_gm= xyz.spt_x; plane.y_gm= xyz.spt_y; plane.z_gm= xyz.spt_z; ocase= 1; status= varkon_vec_projpla ( sdir, plane, ocase, &psdir ); if (status<0) { sprintf(errbuf,"varkon_sur_eval%%sur900*ssilhou"); return(varkon_erpush("SU2943",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*ssilhou spt_x %f spt_y %f spt_z %f\n" , xyz.spt_x, xyz.spt_y, xyz.spt_z ); } #endif /* Function value: */ f_cur= psdir.y_gm*xyz.v_z+psdir.x_gm*xyz.v_y+psdir.z_gm*xyz.v_x -psdir.y_gm*xyz.v_x-psdir.x_gm*xyz.v_z-psdir.z_gm*xyz.v_y; /* !*/ if ( iso_p == 1 ) { /*! */ /* Function derivative value dfdt for a constant V line (iso_p=1): */ dfdt=psdir.y_gm*xyz.v2_z+psdir.x_gm*xyz.v2_y+psdir.z_gm*xyz.v2_x -psdir.y_gm*xyz.v2_x-psdir.x_gm*xyz.v2_z-psdir.z_gm*xyz.v2_y; /* !*/ } else { sprintf(errbuf,"(ssilhou)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*ssilhou t_c %f tglobal %f f_cur %f dfdt %f\n" , t_c , tglobal, f_cur, dfdt ); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** intplan *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function calcultes the function f_cur and dfdt */ /* for calculation case: SURPLAN */ /* */ static short intplan() /* */ /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: Global (static) variables: */ /* Points and derivatives xyz from varkon_sur_eval */ /* Plane pi_x, pi_y, pi_z, pi_d and flag iso_p */ /* Out: Function f_cur and derivative dfdt */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /* --------------- Internal variables ------------------------------*/ /* */ /*----------------------------------------------------------------- */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* Function value: */ f_cur = pi_x*xyz.r_x + pi_y*xyz.r_y + pi_z*xyz.r_z - pi_d; /* !*/ if ( iso_p == 2 ) { /*! */ /* Function derivative value dfdt for a constant V line (iso_p=2): */ dfdt = pi_x*xyz.u_x + pi_y*xyz.u_y + pi_z*xyz.u_z; /* !*/ } else { /*! */ /* Function derivative value dfdt for a constant U line (iso_p=1): */ dfdt = pi_x*xyz.v_x + pi_y*xyz.v_y + pi_z*xyz.v_z; /* !*/ } #ifdef DEBUG if ( dbglev(SURPAC) == 2 && fabs(f_cur) < 0.1 ) { fprintf(dbgfil(SURPAC), "sur900*intplan* tglobal= %15.10f f_cur= %15.10f\n" , tglobal, f_cur); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur900*intplan* tglobal= %f f_cur= %f dfdt= %f\n" , tglobal, f_cur, dfdt ); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*! !*/ #ifdef NOT_USED_FOR_THE_MOMENT /*!************** Internal ** function ** parab2d *******************/ /* */ /* Purpose */ /* -------- */ /* */ /* The function interpolates a start point, an end point and */ /* and an end derivative with a parabola and calculates X */ /* for F(X)=0 */ /* */ static short parab2d() /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In and out: Global (static) variables: */ /* tglobal, tglobal_pre, f_cur, f_pre */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat u1; /* Local coordinate for tglobal_pre */ DBfloat v1; /* Local coordinate for f_pre */ DBfloat b1,b2; /* Coefficients for parabola */ DBfloat v; /* Local coordinate corresp. to parab=0 */ DBfloat u_r1,u_r2; /* Two solutions (roots) */ /*-----------------------------------------------------------------!*/ DBfloat dum1; /* Dummy DBfloat */ /*------------end-of-declarations-----------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "Enter ***** varkon_sur_num1 (sur900) ** parab2d ******\n"); fprintf(dbgfil(SURPAC), "sur900 tglobal_pre= %f f_pre= %f (dfdt= %f)\n", tglobal_pre,f_pre,dfdt); fprintf(dbgfil(SURPAC), "sur900 tglobal= %f f_cur= %f \n",tglobal ,f_cur); } #endif /*! */ /* 1. Check of input data and initiations */ /* -------------------------------------- */ /* */ /* Determine case. Goto _a for tglobal_pre <= tglobal */ /* Goto _b for tglobal_pre > tglobal */ /* !*/ if ( tglobal_pre > tglobal ) goto _b; /* else ( tglobal_pre <= tglobal ) goto _a; */ /*_a: Label for tglobal_pre <= tglobal */ /*! 2.a Create parabola !*/ /* Input points to local coordinates */ u1 = tglobal - tglobal_pre; v1 = f_cur- f_pre ; /*! Coefficients for the parabola !*/ b1 = 2.0*v1 - u1*dfdt; b2 = v1 - b1; /*! */ /* 3.a Calculate new tglobal for parabola= 0 */ /* !*/ v = 0.0 - f_pre; if ( fabs(b2) > TOL1 ) { u_r1 = -b1*u1/(2.0*b2) + SQRT( (b1*u1/(2.0*b2))*(b1*u1/(2.0*b2)) + u1*u1/b2*v ); u_r2 = -b1*u1/(2.0*b2) - SQRT( (b1*u1/(2.0*b2))*(b1*u1/(2.0*b2)) + u1*u1/b2*v ); } else { if ( fabs(b1) > TOL1 ) { u_r1 = v*u1/b1; u_r2 = u_r1; } else { /* Error */ u_r1 = u1; u_r2 = u_r1; } } /* Choose root nearest current u value (tglobal) in the loop */ if ( fabs(u_r1 - u1) < fabs(u_r2-u1) ) { dum1 = tglobal; tglobal = u_r1 + tglobal_pre; tglobal_pre = dum1; f_pre = f_cur; } else { dum1 = tglobal; tglobal = u_r2 + tglobal_pre; tglobal_pre = dum1; f_pre = f_cur; } #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "Exit *** varkon_sur_num1 *** parab2d *** tglobal= %f \n" ,tglobal); /* fprintf(dbgfil(SURPAC), "sur900 Out (_a) tglobal_pre= %f f_pre= %f \n", tglobal_pre,f_pre); fprintf(dbgfil(SURPAC), "sur900 tglobal= %f \n",tglobal ); */ } #endif return(SUCCED); /*! */ _b: /* Label for tglobal_pre > tglobal */ /* 2.b Create parabola */ /* Input points to local coordinates */ /* !*/ u1 = tglobal_pre - tglobal; v1 = f_pre - f_cur ; /*! */ /* Coefficients for the parabola */ /* !*/ b1 = u1*dfdt; b2 = v1 - b1; /*! */ /* 3.b Calculate new tglobal for parabola= 0 */ /* !*/ /* Parabola= 0 in local coordinates */ v = 0.0 - f_cur; if ( fabs(b2) > TOL1 ) { u_r1 = -b1*u1/(2.0*b2) + SQRT( (b1*u1/(2.0*b2))*(b1*u1/(2.0*b2)) + u1*u1/b2*v ); u_r2 = -b1*u1/(2.0*b2) - SQRT( (b1*u1/(2.0*b2))*(b1*u1/(2.0*b2)) + u1*u1/b2*v ); } else { if ( fabs(b1) > TOL1 ) { u_r1 = v*u1/b1; u_r2 = u_r1; } else { /* Error */ u_r1 = u1; u_r2 = u_r1; } } /* Choose root nearest current u value (tglobal) in the loop */ if ( fabs(u_r1 - u1) < fabs(u_r2-u1) ) { tglobal_pre = tglobal; tglobal = u_r1 + tglobal; f_pre = f_cur; } else { tglobal_pre = tglobal; tglobal = u_r2 + tglobal; f_pre = f_cur; } #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "Exit *** varkon_sur_num1 *** parab2d *** tglobal= %f \n" ,tglobal); /* fprintf(dbgfil(SURPAC), "sur900 Out (_b) tglobal_pre= %f f_pre= %f \n", tglobal_pre,f_pre); fprintf(dbgfil(SURPAC), "sur900 tglobal= %f \n",tglobal ); */ } #endif return(SUCCED); } /* End of function */ #endif /* NOT_USED_FOR_THE_MOMENT */ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** addsol ********************/ /* */ /* Purpose */ /* -------- */ /* */ /* Add a solution if it is a new solution (not already existing). */ /* */ /* Note that the identical points criterion (idpoint) 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 */ /* The eye point eye_x, eye_y, eye_z and flag iso_p */ /* 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 f_e; /* Flag equal points: 1=Equal 0=Not equal */ short 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), "sur900*addsol idpoint %f numsol %d \n" , idpoint , (int)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 idpoint. Goto sex if this is the case. */ /* Make additional comparison (call ana2pts) if points are */ /* nearly equal. */ /* !*/ 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 if ( dbglev(SURPAC) == 1 && dist > idpoint ) { fprintf(dbgfil(SURPAC), "sur900*addsol u %10.5f v %10.5f dist %f idpoint %f \n", xyz_a[i_sol-1].u , xyz_a[i_sol-1].v , dist , idpoint ); } #endif /*QQ 1998-09-26 if ( dist <= TOL2 ) goto sex; */ if ( dist <= idpoint ) goto sex; else if ( dist <= 100.0*idpoint ) { f_e = I_UNDEF; status = ana2pts (p_sur,p_pat,xyz_a,i_sol, &f_e); if ( status < 0 ) { sprintf(errbuf,"ana2pts%%(sur900*addsol)"); return(varkon_erpush("SU2973",errbuf)); } /* The additional comparison result ........ */ if (1 == f_e) 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_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } /*! */ /* Recompute xyz with icase_l= 3 for SURPLAN or S_SILH */ /* !*/ if ( ctype == SURPLAN || ctype == S_SILH ) { #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur900*addsol Recompute _sur_eval u %10.5f v %10.5f\n", xyz.u , xyz.v ); fflush(dbgfil(SURPAC)); } #endif icase_l = 3; status= varkon_sur_eval (p_sur,p_pat,icase_l,xyz.u,xyz.v,&xyz); if (status<0) { sprintf(errbuf,"sur210%%(sur900*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) == 2 ) { fprintf(dbgfil(SURPAC), "sur900*addsol x %10.2f y %10.2f z %10.2f\n", xyz.r_x , xyz.r_y , xyz.r_z); fflush(dbgfil(SURPAC)); } #endif sex: /*! Label sex: Solution already exists !*/ return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** analyze********************/ /* */ /* Purpose */ /* -------- */ /* */ /* Determine if the whole, or part of the boundary curve is the */ /* solution. */ /* */ static short analyze (p_sur,p_pat,nstart,icase,xyz_a,p_f_a) /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBSurf *p_sur; /* Surface (ptr) */ DBPatch *p_pat; /* Alloc. area for topol. patch data (ptr) */ DBint nstart; /* The number of restarts */ DBint icase; /* Calculation case for varkon_sur_eval */ /* Out: */ /* Global static variable numsol. */ EVALS xyz_a[SMAX]; /* Array with all solutions (R*3 and UV) */ DBint *p_f_a; /* Flag for analyzis (ptr) */ /* Eq. 0: No point (numsol= 0) */ /* Eq. 1: One end point (numsol= 0) */ /* Eq. 2: Two points (numsol= 2) */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat deltap; /* The parameter t step for points */ DBint i_poi; /* Loop index point along the boundary */ DBint n_poi; /* Number of points along the boundary */ DBint i_start; /* Start point inside c_zero */ DBint i_end; /* End point inside c_zero */ /*-----------------------------------------------------------------!*/ EVALS xyz_start; /* Coordinates and derivatives for i_start */ EVALS xyz_end; /* Coordinates and derivatives for i_end */ 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) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*analyze i_zero %f numsol %d \n" , i_zero , (int)numsol ); } #endif /*! */ /* 1. Check of input data for Debug On */ /* */ /* !*/ #ifdef DEBUG if ( numsol != 0 ) { sprintf(errbuf,"(numsol)%% varkon_sur_num1*analyze"); return(varkon_erpush("SU2993",errbuf)); } if ( nstart < 1 ) { sprintf(errbuf,"(nstart)%% varkon_sur_num1*analyze"); return(varkon_erpush("SU2993",errbuf)); } #endif /*! */ /* Parameter step for along the boundary curve for n_poi= 2*nstart */ /* .... more programming for number of of points (R*3 distance) .. */ /* .... more programming use dfdt, curve tangents ................ */ /* !*/ n_poi = 2*nstart; deltap = (t_e-t_s-comptol)/((DBfloat)n_poi-1.0); /*! */ /* Start loop i_poi= 1, 2, 3, .... n_poi */ /* !*/ *p_f_a = 0; i_start = -1; i_end = -1; for (i_poi = 1; i_poi <= n_poi; i_poi = i_poi+1) { /*! */ /* Current parameter value tglobal */ tglobal = t_s + deltap*((DBfloat)i_poi-1.0); /*! */ /* Call of varkon_sur_eval (sur210) for all surface types */ /* !*/ if (iso_p == 2) status= varkon_sur_eval (p_sur,p_pat,icase,tglobal,t_c,&xyz); else status= varkon_sur_eval (p_sur,p_pat,icase,t_c,tglobal,&xyz); if (status<0) { sprintf(errbuf,"varkon_sur_eval%%varkon_sur_num1*analyze"); return(varkon_erpush("SU2943",errbuf)); } /*! */ /* Calculation of function f_cur and derivative dfdu. */ /* Call of fsilhou for ctype = F_SILH */ /* Call of isilhou for ctype = I_SILH */ /* Call of ssilhou for ctype = S_SILH */ /* Call of fisopho for ctype = F_ISOP */ /* Call of iisopho for ctype = I_ISOP */ /* Call of intplan for ctype = SURPLAN */ /* (Error SU2973 if internal function fails) */ /* !*/ if ( ctype == F_SILH ) status=fsilhou(); else if ( ctype == I_SILH ) status=isilhou(); else if ( ctype == S_SILH ) status=ssilhou(); else if ( ctype == F_ISOP ) status=intplan(); else if ( ctype == I_ISOP ) status=intplan(); else if ( ctype == SURPLAN ) status=intplan(); else { sprintf(errbuf,"(ctype 2)%% varkon_sur_num1 (sur900)"); return(varkon_erpush("SU2993",errbuf)); } if (status<0) { sprintf(errbuf,"f, dfdt%%varkon_sur_num1*analyze"); return(varkon_erpush("SU2973",errbuf)); } /* 1998-09-26 if ( fabs(f_cur) < 50.0*c_zero && i_start < 0 ) */ if ( fabs(f_cur) < i_zero && i_start < 0 ) { xyz_start = xyz; i_start = i_poi; } /*1998-09-26 else if ( fabs(f_cur) < 50.0*c_zero && i_start > 0 ) */ else if ( fabs(f_cur) < i_zero && i_start > 0 ) { xyz_end = xyz; i_end = i_poi; } #ifdef DEBUG if ( dbglev(SURPAC) == 1 && fabs(f_cur) < 100.0*i_zero && fabs(f_cur) >= i_zero ) { fprintf(dbgfil(SURPAC), "sur900*analyze i_start %d i_end %d i_poi %d t %f f_cur %f Outside\n" , (int)i_start, (int)i_end, (int)i_poi, tglobal,f_cur ); } if ( dbglev(SURPAC) == 1 && fabs(f_cur) < i_zero ) { fprintf(dbgfil(SURPAC), "sur900*analyze i_start %d i_end %d i_poi %d t %f f_cur %f Inside\n" , (int)i_start, (int)i_end, (int)i_poi, tglobal,f_cur ); } #endif } /*! */ /* End loop i_poi= 1, 2, 3, .... n_poi */ /* !*/ if ( i_start == 1 && i_end == n_poi ) { tglobal = t_s; xyz = xyz_start; addsol(p_sur,p_pat,xyz_a); /* Add non-existent solution */ tglobal = t_e; xyz = xyz_end; addsol(p_sur,p_pat,xyz_a); /* Add non-existent solution */ *p_f_a = 2; } else if ( i_start == 1 && i_end < 0 ) { *p_f_a = 1; } else if ( i_start == n_poi && i_end < 0 ) { *p_f_a = 1; } else if ( i_start >= 1 && i_end >= n_poi ) { tglobal = t_s + deltap*((DBfloat)i_start-1.0); xyz = xyz_start; addsol(p_sur,p_pat,xyz_a); /* Add non-existent solution */ tglobal = t_s + deltap*((DBfloat)i_end -1.0); xyz = xyz_end; addsol(p_sur,p_pat,xyz_a); /* Add non-existent solution */ *p_f_a = 2; } #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*analyze i_start %d i_end %d n_poi %d *p_f_a %d numsol %d\n" , (int)i_start, (int)i_end, (int)n_poi, (int)*p_f_a, (int)numsol ); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur900*analyze Output numsol %d \n" , (int)numsol ); } #endif /*! */ return(SUCCED); } /* End of function */ /********************************************************************/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*! !*/ /*!************** Internal ** function ** ana2pts********************/ /* */ /* Purpose */ /* -------- */ /* */ /* Additional analysis in order to find out if two solutions are */ /* identical. */ /* */ /* Input points is one existing solution in the output array xyz_a */ /* and the global (static) variable xyz. */ /* This function is used to find a more exact solution when, */ /* for the plane/surface intersect, the plane is nearly tangential */ /* to a patch boundary. The numerical method will in these cases */ /* give multiple solutions which can be at a greater distance */ /* idpoint even if the end criterion ctol (c_zero) is many times */ /* less than idpoint. */ /* An alternative solution could be to run this function (sur900) */ /* with an end criterion c_zero = c_zero/10.0 or c_zero/100.0 */ /* */ /* This function assumes that function analyze finds an interval */ /* of solutions (more programming is actually needed in analyze ..)*/ /* */ /* Note that the input (existing) solution point will be adjusted */ /* if the soltions are equal */ /* */ static short ana2pts (p_sur,p_pat,xyz_a,i_sol,p_f_e) /* */ /* !*/ /********************************************************************/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBSurf *p_sur; /* Surface (ptr) */ DBPatch *p_pat; /* Alloc. area for topol. patch data (ptr) */ /* Points and derivatives xyz from varkon_sur_eval */ /* Out: */ /* Global static variable numsol. */ EVALS xyz_a[SMAX]; /* Array with all solutions (R*3 and UV) */ DBint i_sol; /* Index to comparison point */ DBint *p_f_e; /* Flag equal points: 1=Equal 0=Not equal */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat deltap; /* The parameter t step for points */ DBint i_poi; /* Loop index point along the boundary */ DBint n_poi; /* Number of points along the boundary */ DBfloat t_s_ana; /* Start parameter value for analysis */ DBfloat t_e_ana; /* End parameter value for analysis */ /*-----------------------------------------------------------------!*/ DBint i_start; /* Start point inside c_zero */ DBint i_end; /* End point inside c_zero */ DBint i_minim; /* Start point inside c_zero */ DBfloat min_dist; /* Start parameter value for analysis */ DBint icase_l; /* Calculation case for varkon_sur_eval */ EVALS xyz_input; /* Input coordinates and derivatives */ EVALS xyz_minim; /* Coordinates and derivatives for minimum */ 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) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts idpoint %f numsol %d \n" , idpoint , (int)numsol ); } #endif /*! */ /* 1. Check of input data for Debug On */ /* */ /* !*/ /*! */ /* Initialize output flag to 1 <==> Points are equal */ /* !*/ *p_f_e = 1; /*! */ /* Input point xyz to local parameter (xyz used for the calulation */ /* of the objective function. Necessary to reset xyz at exit of */ /* of this function. Refer to addsol, which calls this function)... */ /* !*/ xyz_input = xyz; /* TODO Use V3MOME */ /* Initialize minimum value TODO: Call initial instead ........ */ xyz_minim = xyz; /* TODO Use V3MOME */ /*! */ /* Parameter step for points between the two solutions */ /* .... more programming for number of of points (R*3 distance) .. */ /* .... more programming use dfdt, curve tangents ................ */ /* !*/ /* Isoparameter curve is U or V */ if ( 2 == iso_p ) { t_s_ana = xyz_a[i_sol-1].u; t_e_ana = xyz.u; } else { t_s_ana = xyz_a[i_sol-1].v; t_e_ana = xyz.v; } n_poi = 5; deltap = (t_e_ana-t_s_ana)/((DBfloat)n_poi-1.0); #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts t_s_ana %f t_e_ana %f deltap %f n_poi %d \n" , t_s_ana, t_e_ana, deltap, (int)n_poi); } #endif /*! */ /* Start loop i_poi= 1, 2, 3, .... n_poi */ /* !*/ i_minim = -1; min_dist = 50000.0; i_start = -1; i_end = -1; /* Calculate all derivatives, since solution may be modified */ /* Loss of performance, not necessary.. Refer to addsol */ icase_l = 3; for (i_poi = 1; i_poi <= n_poi; i_poi = i_poi+1) { /*! */ /* Current parameter value tglobal */ tglobal = t_s_ana + deltap*((DBfloat)i_poi-1.0); /*! */ /* Call of varkon_sur_eval (sur210) for all surface types */ /* !*/ if (iso_p == 2) status= varkon_sur_eval (p_sur,p_pat,icase_l,tglobal,t_c,&xyz); else status= varkon_sur_eval (p_sur,p_pat,icase_l,t_c,tglobal,&xyz); #ifdef DEBUG if ( dbglev(SURPAC) == 1 && status < 0 ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts sur210 failed status= %d \n" , (int)status ); } #endif if (status<0) { sprintf(errbuf,"sur210%%sur900*ana2pts"); return(varkon_erpush("SU2943",errbuf)); } /*! */ /* Calculation of function f_cur and derivative dfdu. */ /* Call of fsilhou for ctype = F_SILH */ /* Call of isilhou for ctype = I_SILH */ /* Call of ssilhou for ctype = S_SILH */ /* Call of fisopho for ctype = F_ISOP */ /* Call of iisopho for ctype = I_ISOP */ /* Call of intplan for ctype = SURPLAN */ /* (Error SU2973 if internal function fails) */ /* !*/ if ( ctype == F_SILH ) status=fsilhou(); else if ( ctype == I_SILH ) status=isilhou(); else if ( ctype == S_SILH ) status=ssilhou(); else if ( ctype == F_ISOP ) status=intplan(); else if ( ctype == I_ISOP ) status=intplan(); else if ( ctype == SURPLAN ) status=intplan(); else { sprintf(errbuf,"(ctype)%% (sur900*ana2pts)"); return(varkon_erpush("SU2993",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 && status < 0 ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts Calculation f, dfdu failed status= %d \n" , (int)status ); } #endif if (status<0) { sprintf(errbuf,"f, dfdt%%varkon_sur_num1*ana2pts"); return(varkon_erpush("SU2973",errbuf)); } /*! */ /* Compare minimum distance */ /* !*/ if ( fabs(f_cur) < min_dist ) { min_dist = fabs(f_cur); xyz_minim = xyz; i_minim = i_poi; } /*! */ /* Change flag to not equal points if objective function value */ /* is greater than criterion i_zero. */ /* !*/ if ( fabs(f_cur) >= i_zero ) *p_f_e = 0; #ifdef DEBUG if ( dbglev(SURPAC) == 1 && fabs(f_cur) >= i_zero ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts *p_f_e %d i_poi %d f_cur %f Outside\n" , (int)*p_f_e , (int)i_poi, f_cur ); } if ( dbglev(SURPAC) == 1 && fabs(f_cur) < i_zero ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts *p_f_e %d i_poi %d f_cur %15.12f dfdt %f\n" , (int)*p_f_e , (int)i_poi, f_cur, dfdt ); } #endif } /*! */ /* End loop i_poi= 1, 2, 3, .... n_poi */ /* !*/ /*! */ /* Modify existing solution if minimum point */ /* !*/ if ( 1 == *p_f_e ) xyz_a[i_sol-1] = xyz_minim; /*! */ /* Reset point xyz to the input value */ /* !*/ xyz = xyz_input; /* TODO Use V3MOME */ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur900*ana2pts n_poi %d *p_f_e %d i_minim %d min_dist %15.10f\n" , (int)n_poi, (int)*p_f_e, (int)i_minim, min_dist ); } #endif return(SUCCED); } /* End of function */ /********************************************************************/