/********************************************************************/ /* */ /* 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_pat_nurbs File: sur242.c */ /* ============================================================= */ /* */ /* Purpose */ /* ------- */ /* */ /* The function calculates coordinates and derivatives for */ /* a given parameter ( u,v ) point on a NURBS surface. */ /* */ /* Author: Gunnar Liden !*/ /* */ /* Revisions */ /* */ /* 1997-12-09 Originally written */ /* 1998-01-28 d2rdudv for a rational surface */ /* 1999-11-21 Free source code modifications */ /* */ /********************************************************************/ /* ------------- Short description of function -----------------*/ /* */ /*sdescr varkon_pat_nurbs NURBS evaluation function */ /* */ /*------------------------------------------------------------- */ /*!---------------------- Theory -----------------------------------*/ /* */ /* Evaluation of coordinates and derivatives */ /* for a NURBS surface */ /* */ /*-----------------------------------------------------------------!*/ /*!-------------- Function calls (internal) ------------------------*/ /* */ static short findspan(); /* Determine the knot span index */ static short derbafun(); /* Basis functions with derivatives */ static short poideriv(); /* Calculate point and derivatives */ static short basisfun(); /* Nonvanishing basis functions */ static short point(); /* Calculate point coordinates */ /* */ #ifdef DEBUG static short initial(); /* Initialization of variables */ #endif /*-----------------------------------------------------------------!*/ /*!-- Static (common) variables for the functions in this file -----*/ /* */ static DBint order_u; /* Order for surface in U direction */ static DBint order_v; /* Order for surface in V direction */ static DBint deg_u; /* Degree for surface in U direction */ static DBint deg_v; /* Degree for surface in V direction */ static DBint nk_u; /* Number of nodes in U direction */ static DBint nk_v; /* Number of nodes in V direction */ static DBint span_u; /* Knot span index in U direction */ static DBint span_v; /* Knot span index in V direction */ static DBfloat b_u[MAX_NURBD];/* Basis functions in U direction */ static DBfloat b_v[MAX_NURBD];/* Basis functions in V direction */ static DBfloat comptol; /* Computer tolerance (accuracy) */ static DBfloat ctol; /* Coordinate tolerance */ /*-----------------------------------------------------------------!*/ /*!-------------- Function calls (external) ------------------------*/ /* */ /* varkon_poi_nurbs * Get address to a NURBS node */ /* varkon_comptol * Retrieve computer tolerance */ /* varkon_ctol * Coordinate tolerance */ /* varkon_erinit * Initial. of error messages */ /* varkon_erpush * Error message to terminal */ /* */ /*-----------------------------------------------------------------!*/ /*!------------ Error messages and warnings ------------------------*/ /* */ /* SU2973 = Internal function xxxxxx failed in varkon_pat_nurbs */ /* SU2993 = Severe program error in varkon_pat_nurbs (sur242). */ /* */ /*-----------------------------------------------------------------!*/ /*!****************** Function **************************************/ /* */ DBstatus varkon_pat_nurbs ( /*-------------- Argument declarations -----------------------------*/ /* */ /* In: */ GMPATNU *p_patnu, /* NURBS patch (ptr) */ DBint icase, /* Calculation case: */ /* Eq. 0: Only coordinates */ /* Eq. 1: Coordinates and dr/du derivative */ /* Eq. 2: Coordinates and dr/dv derivative */ /* Eq. : All derivatives */ DBfloat u_glob, /* Global U value */ DBfloat v_glob, /* Global V value */ EVALS *p_xyz ) /* Coordinates and derivatives (ptr) */ /* Out: */ /* Data to p_xyz */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat d_u[MAX_NURBD][MAX_NURBD]; /* U Basis and derivatives */ DBfloat d_v[MAX_NURBD][MAX_NURBD]; /* V Basis and derivatives */ DBint n_der_out; /* Requested no. of derivatives*/ /* */ /*-----------------------------------------------------------------!*/ DBfloat *p_knot; /* Knot (U or V) vector (ptr) */ DBfloat t_val; /* Current (U or V) parameter value */ DBint i_span; /* Span */ DBint n_der; /* Requested number of derivatives */ DBint p_degree; /* The degree of the basis functions */ char errbuf[80]; /* String for error message fctn erpush */ short status; /* Error code from a called function */ /* ----------------- Theory ----------------------------------------*/ /* */ /* The coordinates and derivatives ........................ */ /* */ /*----------------------------------------------------------------- */ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 Enter**varkon_pat_nurbs p_patnu %d u %f v %f icase %d\n", (int)p_patnu,u_glob,v_glob,(short)icase); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif /*! */ /* 1. Check of input data and initializations */ /* __________________________________________ */ /* */ /* Initialize output coordinates and derivatives and internal */ /* variables for DEBUG on. */ /* Call of initial. */ /* !*/ #ifdef DEBUG status=initial(p_patnu,u_glob,v_glob,p_xyz); if (status<0) { sprintf(errbuf,"initial%%varkon_pat_nurbs (sur242)"); return(varkon_erpush("SU2973",errbuf)); } #endif /*! */ /* Computer tolerance. Call of varkon_comptol (sur753). */ /* !*/ comptol=varkon_comptol(); /*! */ /* Retrieve coordinate tolerance. Call of varkon_ctol (sur751) */ /* !*/ ctol = varkon_ctol(); /*! */ /* Requested number of derivatives n_der_out defined by icase */ /* !*/ n_der_out = 2; /* Requested no. of derivatives */ /* Check that degree is .. to be added */ /*! */ /* 2. Get number of nodes, order (degree) of surface */ /* ____________________________ */ /* */ /* !*/ /* NURBS: Order and degree in U and V direction */ order_u = p_patnu->order_u; order_v = p_patnu->order_v; deg_u = order_u - 1; deg_v = order_v - 1; /* NURBS: Number of nodes in U and V direction */ nk_u = p_patnu->nk_u; nk_v = p_patnu->nk_v; if ( deg_u >= MAX_NURBD || deg_v >= MAX_NURBD ) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 deg_u %d or deg_v %d >= MAX_NURBD %d\n", (int)deg_u, (int)deg_v, MAX_NURBD ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"MAX_NURBD%%sur242"); return(varkon_erpush("SU2993",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 deg_u %d or deg_v %d \n", (int)deg_u, (int)deg_v); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 nk_u %d or nk_v %d \n", (int) nk_u, (int) nk_v ); fflush(dbgfil(SURPAC)); } #endif /*! */ /* 3. Detetermine the knot span indeces span_u,span_v */ /* __________________________________________________ */ /* */ /* Call of internal function findspan. */ /* !*/ status= findspan (p_patnu, u_glob,v_glob); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 findspan failed\n" ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"findspan%%sur242"); return(varkon_erpush("SU2973",errbuf)); } /*! */ /* 4. Compute the nonvanishing basis functions b_u and b_v */ /* _______________________________________________________ */ /* */ /* Call of internal function basisfun. */ /* !*/ status= basisfun (p_patnu, u_glob,v_glob); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 basisfun failed\n" ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"basisfun%%sur242"); return(varkon_erpush("SU2973",errbuf)); } /*! */ /* 5. Calculate coordinates */ /* ________________________ */ /* */ /* Call of internal function point. */ /* !*/ status= point ( p_patnu, p_xyz ); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 point failed\n" ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"point%%sur242"); return(varkon_erpush("SU2973",errbuf)); } /*! */ /* 6. Calculate basis functions and derivatives */ /* ____________________________________________ */ /* */ /* Calls of internal function derbafun */ /* !*/ p_knot = p_patnu->kvec_u; i_span = span_u; t_val = u_glob; p_degree = deg_u; n_der = deg_u; /* Max */ n_der = n_der_out; status=derbafun(p_knot,t_val,i_span,n_der,p_degree,d_u); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 derbafun U failed\n" ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"derbafun%%sur242"); return(varkon_erpush("SU2973",errbuf)); } p_knot = p_patnu->kvec_v; i_span = span_v; t_val = v_glob; p_degree = deg_v; n_der = deg_v; /* Max */ n_der = n_der_out; status=derbafun(p_knot,t_val,i_span,n_der,p_degree,d_v); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 derbafun V failed\n" ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"derbafun%%sur242"); return(varkon_erpush("SU2973",errbuf)); } /*! */ /* 7. Calculate point and derivatives */ /* __________________________________ */ /* */ /* Call of internal function poideriv */ /* !*/ status=poideriv(p_patnu,d_u,d_v,n_der_out,p_xyz ); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 poideriv failed\n" ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"poideriv%%sur242"); return(varkon_erpush("SU2973",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 r_x %f r_y %f r_z %f \n", p_xyz->r_x,p_xyz->r_y,p_xyz->r_z); fprintf(dbgfil(SURPAC), "sur242 u_x %f u_y %f u_z %f \n", p_xyz->u_x,p_xyz->u_y,p_xyz->u_z); fprintf(dbgfil(SURPAC), "sur242 v_x %f v_y %f v_z %f \n", p_xyz->v_x,p_xyz->v_y,p_xyz->v_z); fprintf(dbgfil(SURPAC), "sur242 u2_x %f u2_y %f u2_z %f \n", p_xyz->u2_x,p_xyz->u2_y,p_xyz->u2_z); fprintf(dbgfil(SURPAC), "sur242 v2_x %f v2_y %f v2_z %f \n", p_xyz->v2_x,p_xyz->v2_y,p_xyz->v2_z); fprintf(dbgfil(SURPAC), "sur242 uv_x %f uv_y %f uv_z %f \n", p_xyz->uv_x,p_xyz->uv_y,p_xyz->uv_z); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242 Exit *** varkon_pat_nurbs x= %8.2f y= %8.2f z= %8.2f \n", p_xyz->r_x,p_xyz->r_y,p_xyz->r_z); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif return(SUCCED); } /* End of function */ /*!****************************************************************!*/ #ifdef DEBUG /*!********* Internal ** function **Defined*Only*For*Debug*On********/ /* */ /* Purpose */ /* -------- */ /* */ /* The function initializes the output variables and the */ /* static (common) variables to 1.23456789 and 123456789. */ /* Printout also of the conic lofting patch data. */ static short initial(p_patnu,u_glob,v_glob,p_xyz) /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ GMPATNU *p_patnu; /* NURBS patch (ptr) */ DBfloat u_glob; /* Global U value */ DBfloat v_glob; /* Global V value */ EVALS *p_xyz; /* Pointer to coordinates and derivatives */ /* for a point on a surface */ /* Out: */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBint i_b; /* Loop index basis functions */ /*-----------------------------------------------------------------!*/ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* !*/ if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242 *** initial: p_xyz= %d\n", (int)p_xyz); } /*! */ /* 1. Initialize output variables EVALS and static variables */ /* !*/ /* Coordinates r(u) */ p_xyz->r_x= 1.23456789; p_xyz->r_y= 1.23456789; p_xyz->r_z= 1.23456789; /* Tangent dr/du */ p_xyz->u_x= -1.23456789; p_xyz->u_y= 1.23456789; p_xyz->u_z= 1.23456789; /* Tangent dr/dv */ p_xyz->v_x= 1.23456789; p_xyz->v_y= 1.23456789; p_xyz->v_z= 1.23456789; /* Second derivative d2r/du2 */ p_xyz->u2_x= 1.23456789; p_xyz->u2_y= 1.23456789; p_xyz->u2_z= 1.23456789; /* Second derivative d2r/dv2 */ p_xyz->v2_x= 1.23456789; p_xyz->v2_y= 1.23456789; p_xyz->v2_z= 1.23456789; /* Twist vector d2r/dudv */ p_xyz->uv_x= 1.23456789; p_xyz->uv_y= 1.23456789; p_xyz->uv_z= 1.23456789; /* Surface normal r(u) */ p_xyz->n_x= 1.23456789; p_xyz->n_y= 1.23456789; p_xyz->n_z= 1.23456789; /* Normal derivative dn/du */ p_xyz->nu_x= 1.23456789; p_xyz->nu_y= 1.23456789; p_xyz->nu_z= 1.23456789; /* Normal derivative dn/dv */ p_xyz->nv_x= 1.23456789; p_xyz->nv_y= 1.23456789; p_xyz->nv_z= 1.23456789; /* Maximum and minimum principal curvature */ p_xyz->kmax= 1.23456789; p_xyz->kmin= 1.23456789; /* Direction in R*3 for kmax */ p_xyz->kmax_x= 1.23456789; p_xyz->kmax_y= 1.23456789; p_xyz->kmax_z= 1.23456789; /* Direction in R*3 for kmin */ p_xyz->kmin_x= 1.23456789; p_xyz->kmin_y= 1.23456789; p_xyz->kmin_z= 1.23456789; /* Direction in u,v for kmax */ p_xyz->kmax_u= 1.23456789; p_xyz->kmax_v= 1.23456789; /* Direction in u,v for kmin */ p_xyz->kmin_u= 1.23456789; p_xyz->kmin_v= 1.23456789; /* Internal varibles */ order_u = I_UNDEF; order_v = I_UNDEF; nk_u = I_UNDEF; nk_v = I_UNDEF; span_u = I_UNDEF; span_v = I_UNDEF; deg_u = I_UNDEF; deg_v = I_UNDEF; for ( i_b = 0; i_b < MAX_NURBD; i_b = i_b + 1) { b_u[i_b] = F_UNDEF; b_v[i_b] = F_UNDEF; } return(SUCCED); } /* End of function */ /*!****************************************************************!*/ #endif /*!********* Internal ** function ***********************************/ /* */ /* Purpose */ /* -------- */ /* */ /* Determine the knot span indeces span_u and span_v */ /* The knot span indeces are integers with values */ /* .. I am not sure that the below code is right !!!! */ /* u_span = u_deg-1, ... , nk_u-1 */ /* v_span = v_deg-1, ... , nk_v-1 */ static short findspan (p_patnu, u_glob, v_glob ) /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ GMPATNU *p_patnu; /* NURBS patch (ptr) */ DBfloat u_glob; /* Global U value */ DBfloat v_glob; /* Global U value */ /* Out: */ /* The value of the function */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat u_knot; /* U knot value */ DBfloat v_knot; /* V knot value */ DBint low; /* Index lower limit in binary search */ DBint high; /* Index upper limit in binary search */ DBint mid; /* Index mid point in binary search */ /*-----------------------------------------------------------------!*/ DBfloat *p_uval; /* U knot value (ptr) */ DBfloat *p_vval; /* V knot value (ptr) */ DBint k_iu; /* Loop index U knot values */ DBint n_iter; /* Number of iterations */ char errbuf[80]; /* String for error message fctn erpush */ #ifdef DEBUG DBint k_iv; /* Loop index V knot values */ #endif /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* !*/ /*! */ /* Check that input parameters have valid values for Debug On */ /* !*/ #ifdef DEBUG k_iu = 0; p_uval = p_patnu->kvec_u + k_iu; u_knot = *p_uval; if ( u_glob < u_knot - comptol ) { sprintf(errbuf,"u_glob < Start knot%%sur242*findspan"); return(varkon_erpush("SU2993",errbuf)); } k_iu = nk_u; p_uval = p_patnu->kvec_u + k_iu - 1; u_knot = *p_uval; if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan u_glob %f End value u_knot= %8.2f \n", u_glob, u_knot); fflush(dbgfil(SURPAC)); /* To file from buffer */ } if ( u_glob > u_knot + comptol ) { sprintf(errbuf,"u_glob > End knot%%sur242*findspan"); return(varkon_erpush("SU2993",errbuf)); } k_iv = 1; p_vval = p_patnu->kvec_v + k_iv - 1; v_knot = *p_vval; if ( v_glob < v_knot - comptol ) { sprintf(errbuf,"v_glob < Start knot%%sur242*findspan"); return(varkon_erpush("SU2993",errbuf)); } k_iv = nk_v; p_vval = p_patnu->kvec_v + k_iv - 1; v_knot = *p_vval; if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan v_glob %f End value v_knot= %8.2f \n", v_glob, v_knot); fflush(dbgfil(SURPAC)); /* To file from buffer */ } if ( v_glob > v_knot + comptol ) { sprintf(errbuf,"v_glob > End knot%%sur242*findspan"); return(varkon_erpush("SU2993",errbuf)); } #endif /* Initialization of local variables for Debug On */ #ifdef DEBUG low = I_UNDEF; high = I_UNDEF; mid = I_UNDEF; #endif n_iter = 0; #ifdef DEBUG for ( k_iu = 1; k_iu <= nk_u; ++k_iu ) { p_uval = p_patnu->kvec_u + k_iu - 1; u_knot = *p_uval; if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan Index (span) %4d U knot value u_knot= %8.2f \n", (int)k_iu-1, u_knot); fflush(dbgfil(SURPAC)); /* To file from buffer */ } } for ( k_iv = 1; k_iv <= nk_v; ++k_iv ) { p_vval = p_patnu->kvec_v + k_iv - 1; v_knot = *p_vval; if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan Index (span) %4d V knot value v_knot= %8.2f \n", (int)k_iv-1, v_knot); fflush(dbgfil(SURPAC)); /* To file from buffer */ } } #endif /*! */ /* 1. Determine the U knot span index span_u */ /* ------------------------------------------ */ /* */ /* !*/ /* Check if the input parameter value is equal to the start knot */ p_uval = p_patnu->kvec_u + deg_u; u_knot = *p_uval; if ( fabs(u_glob - u_knot) < comptol ) { span_u = deg_u; goto u_end; } /* Check if the input parameter value is equal to the end knot */ k_iu = nk_u; p_uval = p_patnu->kvec_u + nk_u - deg_u - 1; u_knot = *p_uval; if ( fabs(u_glob - u_knot) < comptol ) { span_u = nk_u - 1 - deg_u - 1; goto u_end; } /* Initialize lower and upper limits for the binary search */ /* The lower limit is the degree of the NURBS surface and */ /* the upper limit is the number of knot values. */ low = deg_u; high = nk_u - deg_u - 1; /* The mid point is the truncated value (low+high)/2 */ mid = (DBint)((low+high)/2); #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*findspan low %4d high %4d mid %4d Start values\n", (int)low, (int)high, (int)mid ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan U[low] %4.2f U[high] %4.2f U[mid] %4.2f U[mid+1] %4.2f\n", *(p_patnu->kvec_u + low) , *(p_patnu->kvec_u + high), *(p_patnu->kvec_u + mid) , *(p_patnu->kvec_u + mid + 1 )); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif /* Do the binary search */ while ( u_glob < *(p_patnu->kvec_u + mid) || u_glob >= *(p_patnu->kvec_u + mid + 1) ) { n_iter = n_iter + 1; if ( n_iter > 1000 ) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242**findspan n_iter >1000 u_glob %f\n", u_glob ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"n_iter U%%sur242*findspan"); return(varkon_erpush("SU2993",errbuf)); } if ( u_glob < *(p_patnu->kvec_u + mid) ) high = mid; else low = mid; /* The new mid point is the truncated value (low+high)/2 */ mid = (DBint)((low+high)/2); #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*findspan low %4d high %4d mid %4d u_glob %f\n", (int)low, (int)high, (int)mid, u_glob ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*findspan U[low] %4.2f U[high] %4.2f U[mid] %4.2f U[mid+1] %4.2f\n", *(p_patnu->kvec_u + low) , *(p_patnu->kvec_u + high), *(p_patnu->kvec_u + mid) , *(p_patnu->kvec_u + mid + 1) ); fflush(dbgfil(SURPAC)); } #endif } /* End loop binary search */ /* The output knot span index is found */ span_u = mid; u_end:; /* Label: Input U is the end knot U value */ #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*findspan u_glob= %6.4f span_u= %d U[] %f U[+1] %f\n", u_glob, (int)span_u, *(p_patnu->kvec_u + span_u ) , *(p_patnu->kvec_u + span_u +1)); fflush(dbgfil(SURPAC)); } #endif /*! */ /* 2. Determine the V knot span index span_v */ /* ------------------------------------------ */ /* */ /* !*/ /* Check if the input parameter value is equal to the start knot */ p_vval = p_patnu->kvec_v + deg_v; v_knot = *p_vval; if ( fabs(v_glob - v_knot) < comptol ) { span_v = deg_v; goto v_end; } /* Check if the input parameter value is equal to the end knot */ p_vval = p_patnu->kvec_v + nk_v - deg_v - 1; v_knot = *p_vval; if ( fabs(v_glob - v_knot) < comptol ) { span_v = nk_v - 1 - deg_v - 1; goto v_end; } /* Initialize lower and upper limits for the binary search */ /* The lower limit is the degree of the NURBS surface and */ /* the upper limit is the number of knot values. */ low = deg_v; high = nk_v - deg_v - 1; /* The mid point is the truncated value (low+high)/2 */ mid = (DBint)((low+high)/2); #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*findspan low %4d high %4d mid %4d Start values\n", (int)low, (int)high, (int)mid ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*findspan V[low] %4.2f V[high] %4.2f V[mid] %4.2f V[mid+1] %4.2f\n", *(p_patnu->kvec_v + low) , *(p_patnu->kvec_v + high), *(p_patnu->kvec_v + mid) , *(p_patnu->kvec_v + mid + 1) ); fflush(dbgfil(SURPAC)); } #endif /* Do the binary search */ while ( v_glob < *(p_patnu->kvec_v + mid) || v_glob >= *(p_patnu->kvec_v + mid + 1) ) { n_iter = n_iter + 1; if ( n_iter > 1000 ) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242**findspan n_iter >1000 v_glob %f\n", v_glob ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"n_iter V%%sur242*findspan"); return(varkon_erpush("SU2993",errbuf)); } if ( v_glob < *(p_patnu->kvec_v + mid) ) high = mid; else low = mid; /* The new mid point is the truncated value (low+high)/2 */ mid = (DBint)((low+high)/2); #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan low %4d high %4d mid %4d v_glob %f\n", (int)low, (int)high, (int)mid, v_glob ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan V[low] %4.2f V[high] %4.2f V[mid] %4.2f V[mid+1] %4.2f\n", *(p_patnu->kvec_v + low) , *(p_patnu->kvec_v + high), *(p_patnu->kvec_v + mid) , *(p_patnu->kvec_v + mid + 1) ); fflush(dbgfil(SURPAC)); } #endif } /* End loop binary search span_v */ /* The output knot span index is found */ span_v = mid; v_end:; /* Label: Input V is the end knot V value */ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan u_glob= %6.4f span_u= %d U[] %f U[+1] %f\n", u_glob, (int)span_u, *(p_patnu->kvec_u + span_u ) , *(p_patnu->kvec_u + span_u +1)); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*findspan v_glob= %6.4f span_v= %d V[] %f V[+1] %f\n", v_glob, (int)span_v, *(p_patnu->kvec_v + span_v ) , *(p_patnu->kvec_v + span_v +1)); fflush(dbgfil(SURPAC)); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!********* Internal ** function ***********************************/ /* */ /* Purpose */ /* -------- */ /* */ /* Compute the nonvanishing basis functions b_u and b_v */ static short basisfun (p_patnu, u_glob, v_glob ) /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ GMPATNU *p_patnu; /* NURBS patch (ptr) */ DBfloat u_glob; /* Global U value */ DBfloat v_glob; /* Global U value */ /* Out: */ /* The value of the function */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat u_knot1; /* U knot value */ DBfloat u_knot2; /* U knot value */ DBfloat v_knot1; /* V knot value */ DBfloat v_knot2; /* V knot value */ DBfloat temp; /* Temporary basis function value */ DBfloat saved; /* Saved basis function value */ DBfloat left [MAX_NURBD];/* Parameter interval left */ DBfloat right[MAX_NURBD];/* Parameter interval right */ DBint j_n; /* Loop index */ DBint r_n; /* Loop index */ /*-----------------------------------------------------------------!*/ DBfloat denom; /* Denominator */ DBfloat *p_uval; /* U knot value (ptr) */ DBfloat *p_vval; /* V knot value (ptr) */ char errbuf[80]; /* String for error message fctn erpush */ #ifdef DEBUG DBfloat u_sum; /* Sum of basis functions U */ DBfloat v_sum; /* Sum of basis functions V */ #endif /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* !*/ /* Initialization of basis function values for Debug On */ #ifdef DEBUG for ( j_n = 0; j_n <= deg_u; j_n++) { left [j_n] = F_UNDEF; right[j_n] = F_UNDEF; } #endif /*! */ /* 1. Basis functions for U */ /* ------------------------- */ /* */ /* !*/ /*! */ /* First basis function */ /* !*/ b_u[0] = 1.0; for ( j_n = 1; j_n <= deg_u; j_n++) { p_uval = p_patnu->kvec_u + span_u + 1 - j_n; u_knot1 = *p_uval; left [j_n] = u_glob - u_knot1; p_uval = p_patnu->kvec_u + span_u + j_n; u_knot2 = *p_uval; right[j_n] = u_knot2 - u_glob; saved = 0.0; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun j_n= %4d left[%4d]= %f right[%4d] %f\n", (int)j_n, (int)j_n, left[j_n], (int)j_n, right[j_n] ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun u_knot1 %f u_knot2 %f \n", u_knot1, u_knot2 ); fflush(dbgfil(SURPAC)); /* To file from buffer */ } #endif for ( r_n = 0; r_n < j_n; r_n++) { denom = right[r_n+1]+left[j_n-r_n]; if ( fabs(denom) < 10.0*comptol) { sprintf(errbuf,"denom U = 0%%sur242*basisfun"); return(varkon_erpush("SU2993",errbuf)); } temp = b_u[r_n] / denom; b_u[r_n] = saved + right[r_n+1]*temp; saved = left[j_n-r_n]*temp; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun j_n= %4d r_n %d b_u[%4d] %f left[ ]= %f right[ ] %f\n", (int)j_n, (int)r_n, (int)r_n, b_u[r_n], left[j_n-r_n], right[r_n+1] ); fflush(dbgfil(SURPAC)); } #endif } /* End loop r_n */ b_u[j_n] = saved; } /* End loop j_n */ /*! */ /* 2. Basis functions for V */ /* ------------------------- */ /* */ /* !*/ /*! */ /* First basis function */ /* !*/ b_v[0] = 1.0; for ( j_n = 1; j_n <= deg_v; j_n++) { p_vval = p_patnu->kvec_v + span_v + 1 - j_n; v_knot1 = *p_vval; left [j_n] = v_glob - v_knot1; p_vval = p_patnu->kvec_v + span_v + j_n; v_knot2 = *p_vval; right[j_n] = v_knot2 - v_glob; saved = 0.0; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun j_n= %4d left[%4d]= %f right[%4d] %f\n", (int)j_n, (int)j_n, left[j_n], (int)j_n, right[j_n] ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun v_knot1 %f v_knot2 %f \n", v_knot1, v_knot2 ); fflush(dbgfil(SURPAC)); } #endif for ( r_n = 0; r_n < j_n; r_n++) { denom = right[r_n+1]+left[j_n-r_n]; if ( fabs(denom) < 10.0*comptol) { sprintf(errbuf,"denom V = 0%%sur242*basisfun"); return(varkon_erpush("SU2993",errbuf)); } temp = b_v[r_n] / denom; b_v[r_n] = saved + right[r_n+1]*temp; saved = left[j_n-r_n]*temp; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun j_n= %4d r_n %d b_v[%4d] %f left[ ]= %f right[ ] %f\n", (int)j_n, (int)r_n, (int)r_n, b_v[r_n], left[j_n-r_n], right[r_n+1] ); fflush(dbgfil(SURPAC)); } #endif } /* End loop r_n */ b_v[j_n] = saved; } /* End loop j_n */ /*! */ /* 3. Check result for Debug On */ /* ----------------------------- */ /* */ /* !*/ #ifdef DEBUG u_sum = 0.0; for ( j_n = 0; j_n <= deg_u; j_n++) { u_sum = b_u[j_n] + u_sum; } v_sum = 0.0; for ( j_n = 0; j_n <= deg_v; j_n++) { v_sum = b_v[j_n] + v_sum; } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun u_sum= %f v_sum= %f\n", u_sum, v_sum ); fflush(dbgfil(SURPAC)); } if ( fabs(u_sum-1.0) > comptol ) { sprintf(errbuf,"u_sum not 1%%sur242*basisfun"); return(varkon_erpush("SU2993",errbuf)); } if ( fabs(v_sum-1.0) > comptol ) { sprintf(errbuf,"v_sum not 1%%sur242*basisfun"); return(varkon_erpush("SU2993",errbuf)); } #endif #ifdef DEBUG for ( j_n = 0; j_n <= deg_u; j_n++) { if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun b_u[%d]= %f \n", (int)j_n, b_u[j_n] ); fflush(dbgfil(SURPAC)); } } for ( j_n = 0; j_n <= deg_v; j_n++) { if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*basisfun b_v[%d]= %f \n", (int)j_n, b_v[j_n] ); fflush(dbgfil(SURPAC)); } } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!********* Internal ** function ***********************************/ /* */ /* Purpose */ /* -------- */ /* */ /* Compute point coordinates */ static short point (p_patnu, p_xyz ) /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ GMPATNU *p_patnu; /* NURBS patch (ptr) */ EVALS *p_xyz; /* Coordinates and derivatives (ptr) */ /* Out: */ /* The value of the function */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat tempx[MAX_NURBD];/* Temporary vector to store ... */ DBfloat tempy[MAX_NURBD];/* Temporary vector to store ... */ DBfloat tempz[MAX_NURBD];/* Temporary vector to store ... */ DBfloat tempw[MAX_NURBD];/* Temporary vector to store ... */ DBfloat polyx; /* Polygon point X value */ DBfloat polyy; /* Polygon point Y value */ DBfloat polyz; /* Polygon point Z value */ DBfloat polyw; /* Polygon point w value */ DBfloat s_x; /* Surface point X coordinate */ DBfloat s_y; /* Surface point Y coordinate */ DBfloat s_z; /* Surface point Z coordinate */ DBfloat s_w; /* Surface point w coordinate */ DBint l_n; /* Loop index */ DBint k_n; /* Loop index */ /* For varkon_poi_nurbs (sur243) */ DBint rwcase; /* Read/write case: */ /* Eq. 0: Read node */ /* Eq. 1: Write node */ /* Eq. 2: Initialize values */ /* Eq. 3: Printout of polygon points */ /* Eq. 4: Check NURBS data */ DBint i_up; /* Point address in U direction */ DBint i_vp; /* Point address in V direction */ DBHvector node; /* Polygon point */ /*-----------------------------------------------------------------!*/ DBint status; /* Error code from a called function */ char errbuf[80]; /* String for error message fctn erpush */ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* !*/ /* Initialization of local variables for Debug On */ #ifdef DEBUG i_up = I_UNDEF; i_vp = I_UNDEF; polyx = F_UNDEF; polyy = F_UNDEF; polyz = F_UNDEF; polyw = F_UNDEF; s_x = F_UNDEF; s_y = F_UNDEF; s_z = F_UNDEF; s_w = F_UNDEF; for ( l_n = 0; l_n < MAX_NURBD; l_n++) { tempx[l_n] = F_UNDEF; tempy[l_n] = F_UNDEF; tempz[l_n] = F_UNDEF; tempw[l_n] = F_UNDEF; } #endif /*! */ /* Printout of polygon points for Debug On */ /* Call of varkon_poi_nurbs (sur243). */ /* !*/ #ifdef DEBUG rwcase = 3; i_up = I_UNDEF; i_vp = I_UNDEF; status= varkon_poi_nurbs (p_patnu, rwcase, i_up, i_vp, &node); if ( status < 0 ) { sprintf(errbuf, "sur243 List%%sur242*point"); return(varkon_erpush("SU2943",errbuf)); } #endif for ( l_n = 0; l_n <= deg_v; l_n++) { tempx[l_n] = 0.0; tempy[l_n] = 0.0; tempz[l_n] = 0.0; tempw[l_n] = 0.0; for ( k_n = 0; k_n <= deg_u; k_n++) { i_up = span_u - deg_u + k_n; i_vp = span_v - deg_v + l_n; rwcase = 0; status= varkon_poi_nurbs (p_patnu, rwcase, i_up, i_vp, &node); if ( status < 0 ) { sprintf(errbuf, "sur243 Read%%sur242*point"); return(varkon_erpush("SU2943",errbuf)); } polyx = node.x_gm; polyy = node.y_gm; polyz = node.z_gm; polyw = node.w_gm; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "i_up %3d i_vp %3d x %8.2f y %8.2f z %8.2f w %8.4f\n", (int)i_up,(int)i_vp, polyx, polyy, polyz, polyw ); fflush(dbgfil(SURPAC)); } #endif tempx[l_n] = tempx[l_n] + b_u[k_n]*polyx; tempy[l_n] = tempy[l_n] + b_u[k_n]*polyy; tempz[l_n] = tempz[l_n] + b_u[k_n]*polyz; tempw[l_n] = tempw[l_n] + b_u[k_n]*polyw; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "l_n %3d k_n %3d polyx %8.4f b_u[k_n] %8.4f tempx[l_n] %8.2f\n", (int)l_n ,(int)k_n, polyx, b_u[k_n],tempx[l_n]); fflush(dbgfil(SURPAC)); } #endif } } s_x = 0.0; s_y = 0.0; s_z = 0.0; s_w = 0.0; for ( l_n = 0; l_n <= deg_v; l_n++) { s_x = s_x + b_v[l_n]*tempx[l_n]; s_y = s_y + b_v[l_n]*tempy[l_n]; s_z = s_z + b_v[l_n]*tempz[l_n]; s_w = s_w + b_v[l_n]*tempw[l_n]; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "l_n %3d b_v[l_n] %8.4f tempx[l_n] %8.2f s_x %f\n", (int)l_n , b_v[l_n],tempx[l_n], s_x); fflush(dbgfil(SURPAC)); } #endif } if ( fabs(s_w) < 100.0*comptol ) { sprintf(errbuf,"Rational denominator = 0%%sur242*point"); return(erpush("SU2993",errbuf)); } p_xyz->r_x = s_x/s_w; p_xyz->r_y = s_y/s_w; p_xyz->r_z = s_z/s_w; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*POINT X= %8.2f Y= %8.2f Z= %8.2f w= %8.4f\n", s_x,s_y,s_z,s_w); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*point x= %8.2f y= %8.2f z= %8.2f \n", p_xyz->r_x,p_xyz->r_y,p_xyz->r_z); fflush(dbgfil(SURPAC)); } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!********* Internal ** function ***********************************/ /* */ /* Purpose */ /* -------- */ /* */ /* Compute nonvanishing basis functions and their derivatives */ static short derbafun (p_knot,t_val,i_span,n_der,p_degree,ders ) /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBfloat *p_knot; /* Knot (U or V) vector (ptr) */ DBfloat t_val; /* Current (U or V) parameter value */ DBint i_span; /* Span */ DBint n_der; /* Requested number of derivatives */ DBint p_degree; /* The degree of the basis functions */ /* Out: */ DBfloat ders[MAX_NURBD][MAX_NURBD]; /* Basis ftns & derivatives */ /* */ /* */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat k_value_1; /* Knot value 1 from input vector */ DBfloat k_value_2; /* Knot value 2 from input vector */ DBfloat temp; /* Temporary basis function value */ DBfloat saved; /* Saved basis function value */ DBfloat left [MAX_NURBD];/* Parameter interval left */ DBfloat right[MAX_NURBD];/* Parameter interval right */ DBint j_n; /* Loop index */ DBint r_n; /* Loop index */ DBint k_n; /* Loop index */ DBint j1; /* Loop end value */ DBint j2; /* Loop end value */ DBint s1; /* Index in array */ DBint s2; /* Index in array */ DBint rk; /* Index in array */ DBint pk; /* Index in array */ DBfloat deriv; /* Derivative */ /*-----------------------------------------------------------------!*/ /* Basis functions and knot differences */ DBfloat ndu[MAX_NURBD][MAX_NURBD]; /* To store intermediate results */ DBfloat a[2][MAX_NURBD]; DBfloat denom; /* Denominator */ DBfloat *p_tval; /* Knot value (ptr) */ char errbuf[80]; /* String for error message fctn erpush */ #ifdef DEBUG DBfloat u_sum; /* Sum of basis functions */ #endif /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*derbafun Enter t_val %f i_span %d n_der %d p_degree %d\n", t_val,(int)i_span, (int)n_der, (int)p_degree ); fflush(dbgfil(SURPAC)); } #endif /* Initialization of basis function values for Debug On */ #ifdef DEBUG for ( j_n = 0; j_n < MAX_NURBD; j_n++) { left [j_n] = F_UNDEF; right[j_n] = F_UNDEF; a[0][j_n] = F_UNDEF; a[1][j_n] = F_UNDEF; } for ( j_n = 0; j_n < MAX_NURBD; j_n++) { for ( r_n = 0; r_n < MAX_NURBD; r_n++) { ndu[j_n][r_n] = F_UNDEF; ders[j_n][r_n] = F_UNDEF; } } denom = F_UNDEF; temp = F_UNDEF; saved = F_UNDEF; k_value_1 = F_UNDEF; k_value_2 = F_UNDEF; u_sum = F_UNDEF; #endif /* First basis function */ ndu[0][0] = 1.0; for ( j_n = 1; j_n <= p_degree; j_n++) { /* Knot value */ p_tval = p_knot + i_span + 1 - j_n; k_value_1 = *p_tval; left [j_n] = t_val - k_value_1; p_tval = p_knot + i_span + j_n; k_value_2 = *p_tval; right[j_n] = k_value_2 - t_val; saved = 0.0; for ( r_n = 0; r_n < j_n; r_n++) { ndu[j_n][r_n] = right[r_n+1] + left[j_n-r_n]; denom = ndu[j_n][r_n]; if ( fabs(denom) < 10.0*comptol ) { sprintf(errbuf,"denom= 0%%sur242*derbafun"); return(varkon_erpush("SU2993",errbuf)); } temp = ndu[r_n][j_n-1]/denom; ndu[r_n][j_n] = saved + right[r_n+1]*temp; saved = left[j_n-r_n]*temp; } /* End loop r_n */ ndu[j_n][j_n] = saved; } /* End loop j_n */ /* Load the basis functions */ for ( j_n = 0; j_n <= p_degree; j_n++) { ders[0][j_n] = ndu[j_n][p_degree]; } /* End loop j_n */ #ifdef DEBUG for ( j_n = 0; j_n <= p_degree; j_n++) { if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*derbafun ders[%d][%d]= %f Prior to derivative calc.\n", (int)0 , (int)j_n, ders[0 ][j_n] ); fflush(dbgfil(SURPAC)); } } #endif /* Calculation of derivatives */ /* Loop over function index */ for ( r_n = 0; r_n <= p_degree; r_n++) { s1 = 0; s2 = 1; a[0][0] = 1.0; /* Loop to compute the kth derivative */ for ( k_n = 1; k_n <= n_der; k_n++) { deriv = 0.0; rk = r_n - k_n; pk = p_degree - k_n; if ( r_n >= k_n ) { denom = ndu[pk+1][rk]; if ( fabs(denom) < 10.0*comptol ) { sprintf(errbuf,"denom= 0 (2)%%sur242*derbafun"); return(varkon_erpush("SU2993",errbuf)); } a[s2][0] = a[s1][0]/ndu[pk+1][rk]; deriv = a[s2][0]*ndu[rk][pk]; } if (rk >= -1 ) j1 = 1; else j1 = -rk; if ( r_n-1 <= pk ) j2 = k_n - 1; else j2 = p_degree - r_n; for ( j_n = j1; j_n <= j2; j_n++) { denom = ndu[pk+1][rk+j_n]; if ( fabs(denom) < 10.0*comptol ) { sprintf(errbuf,"denom= 0 (3)%%sur242*derbafun"); return(varkon_erpush("SU2993",errbuf)); } a[s2][j_n] = (a[s1][j_n]-a[s1][j_n-1])/denom; deriv = deriv + a[s2][j_n]*ndu[rk+j_n][pk]; } /* End loop j_n */ if ( r_n <= pk ) { denom = ndu[pk+1][r_n]; if ( fabs(denom) < 10.0*comptol ) { sprintf(errbuf,"denom= 0 (4)%%sur242*derbafun"); return(varkon_erpush("SU2993",errbuf)); } a[s2][k_n] = -a[s1][k_n-1]/denom; deriv = deriv + a[s2][k_n]*ndu[r_n][pk]; } /* End r_n <= pk */ ders[k_n][r_n] = deriv; /* Switch rows */ j_n = s1; s1 = s2; s2 = j_n; } /* End loop k_n */ } /* End loop r_n */ /* Multiply with correct factors */ r_n = p_degree; for ( k_n = 1; k_n <= n_der; k_n++) { for ( j_n = 0; j_n <= p_degree; j_n++) { ders[k_n][j_n] = ders[k_n][j_n]*r_n; } /* End loop j_n */ r_n = r_n*(p_degree-k_n); } /* End loop k_n */ #ifdef DEBUG for ( j_n = 0; j_n <= p_degree; j_n++) { for ( r_n = 0; r_n <= p_degree; r_n++) { if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*derbafun ders[%d][%d]= %f \n", (int)j_n, (int)r_n, ders[j_n][r_n] ); fflush(dbgfil(SURPAC)); } } } #endif return(SUCCED); } /* End of function */ /********************************************************************/ /*!********* Internal ** function ***********************************/ /* */ /* Purpose */ /* -------- */ /* */ /* Compute point coordinates */ static short poideriv (p_patnu,d_u,d_v, n_der, p_xyz ) /* */ /******************************************************************!*/ /*!------------- Argument declarations -----------------------------*/ /* */ /* In: */ /* Static variables */ /* deg_u Degree in U */ /* deg_v Degree in V */ /* span_u Span in U */ /* span_v Span in V */ /* */ /* */ GMPATNU *p_patnu; /* NURBS patch (ptr) */ DBfloat d_u[MAX_NURBD][MAX_NURBD];/* U Basis and derivatives */ DBfloat d_v[MAX_NURBD][MAX_NURBD];/* V Basis and derivatives */ DBint n_der; /* Requested no. of derivatives */ EVALS *p_xyz; /* Coord's and derivatives (ptr)*/ /* Out: */ /* The value of the function */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ DBfloat SKL[MAX_NURBD][MAX_NURBD][4]; /* Point and derivatives */ DBfloat tempx[MAX_NURBD]; /* Temporary X vector */ DBfloat tempy[MAX_NURBD]; /* Temporary Y vector */ DBfloat tempz[MAX_NURBD]; /* Temporary Z vector */ DBfloat tempw[MAX_NURBD]; /* Temporary w vector */ DBfloat polyx; /* Polygon point X value */ DBfloat polyy; /* Polygon point Y value */ DBfloat polyz; /* Polygon point Z value */ DBfloat polyw; /* Polygon point w value */ /*-----------------------------------------------------------------!*/ DBfloat h_x; /* Surface point X coordinate, homogenous */ DBfloat h_y; /* Surface point Y coordinate, homogenous */ DBfloat h_z; /* Surface point Z coordinate, homogenous */ DBfloat h_w; /* Surface point w coordinate, homogenous */ DBfloat hu_x ; /* Surface dR/du homogenous */ DBfloat hu_y ; /* Surface dR/du homogenous */ DBfloat hu_z ; /* Surface dR/du homogenous */ DBfloat hu_w ; /* Surface dR/du homogenous */ DBfloat hv_x ; /* Surface dR/dv homogenous */ DBfloat hv_y ; /* Surface dR/dv homogenous */ DBfloat hv_z ; /* Surface dR/dv homogenous */ DBfloat hv_w ; /* Surface dR/dv homogenous */ DBfloat hu2_x ; /* Surface d2R/du2 homogenous */ DBfloat hu2_y ; /* Surface d2R/du2 homogenous */ DBfloat hu2_z ; /* Surface d2R/du2 homogenous */ DBfloat hu2_w ; /* Surface d2R/du2 homogenous */ DBfloat hv2_x ; /* Surface d2R/dv2 homogenous */ DBfloat hv2_y ; /* Surface d2R/dv2 homogenous */ DBfloat hv2_z ; /* Surface d2R/dv2 homogenous */ DBfloat hv2_w ; /* Surface d2R/dv2 homogenous */ DBfloat huv_x ; /* Surface d2R/dudv homogenous */ DBfloat huv_y ; /* Surface d2R/dudv homogenous */ DBfloat huv_z ; /* Surface d2R/dudv homogenous */ DBfloat huv_w ; /* Surface d2R/dudv homogenous */ DBfloat h_d2xdudv; /* Surface d2R/dudv homogenous */ DBfloat h_d2ydudv; /* Surface d2R/dudv homogenous */ DBfloat h_d2zdudv; /* Surface d2R/dudv homogenous */ DBfloat h_d2wdudv; /* Surface d2R/dudv homogenous */ /* For varkon_poi_nurbs (sur243) */ DBint rwcase; /* Read/write case: */ /* Eq. 0: Read node */ /* Eq. 1: Write node */ /* Eq. 2: Initialize values */ /* Eq. 3: Printout of polygon points */ /* Eq. 4: Check NURBS data */ DBint i_up; /* Point address in U direction */ DBint i_vp; /* Point address in V direction */ DBHvector node; /* Polygon point */ DBint l_n; /* Loop index */ DBint k_n; /* Loop index */ DBint s_n; /* Loop index */ DBint r_n; /* Loop index */ DBint min_du; /* Minimum of n_der and deg_u */ DBint min_dv; /* Minimum of n_der and deg_v */ DBint min_dd; /* Minimum of d_n-k_n and min_dv */ DBint status; /* Error code from a called function */ char errbuf[80]; /* String for error message fctn erpush */ /*--------------end-of-declarations---------------------------------*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* !*/ /* Initialization of local variables for Debug On */ #ifdef DEBUG for ( l_n = 0; l_n < MAX_NURBD; l_n++) { for ( k_n = 0; k_n < MAX_NURBD; k_n++) { SKL[l_n][k_n][0] = F_UNDEF; SKL[l_n][k_n][1] = F_UNDEF; SKL[l_n][k_n][2] = F_UNDEF; SKL[l_n][k_n][3] = F_UNDEF; } } for ( l_n = 0; l_n < MAX_NURBD; l_n++) { tempx[l_n] = F_UNDEF; tempy[l_n] = F_UNDEF; tempz[l_n] = F_UNDEF; tempw[l_n] = F_UNDEF; } l_n = I_UNDEF; k_n = I_UNDEF; r_n = I_UNDEF; s_n = I_UNDEF; min_du = I_UNDEF; min_dv = I_UNDEF; min_dd = I_UNDEF; status = I_UNDEF; rwcase = I_UNDEF; i_up = I_UNDEF; i_vp = I_UNDEF; node.x_gm = F_UNDEF; node.y_gm = F_UNDEF; node.z_gm = F_UNDEF; node.w_gm = F_UNDEF; polyx = F_UNDEF; polyy = F_UNDEF; polyz = F_UNDEF; polyw = F_UNDEF; h_x = F_UNDEF; h_y = F_UNDEF; h_z = F_UNDEF; h_w = F_UNDEF; hu_x = F_UNDEF; hu_y = F_UNDEF; hu_z = F_UNDEF; hu_w = F_UNDEF; hv_x = F_UNDEF; hv_y = F_UNDEF; hv_z = F_UNDEF; hv_w = F_UNDEF; hu2_x = F_UNDEF; hu2_y = F_UNDEF; hu2_z = F_UNDEF; hu2_w = F_UNDEF; hv2_x = F_UNDEF; hv2_y = F_UNDEF; hv2_z = F_UNDEF; hv2_w = F_UNDEF; huv_x = F_UNDEF; huv_y = F_UNDEF; huv_z = F_UNDEF; huv_w = F_UNDEF; h_d2xdudv = F_UNDEF; h_d2ydudv = F_UNDEF; h_d2zdudv = F_UNDEF; h_d2wdudv = F_UNDEF; #endif /*! */ /* Printout of polygon points for Debug On */ /* Call of varkon_poi_nurbs (sur243). */ /* !*/ #ifdef DEBUG rwcase = 3; i_up = I_UNDEF; i_vp = I_UNDEF; status= varkon_poi_nurbs (p_patnu, rwcase, i_up, i_vp, &node); if ( status < 0 ) { sprintf(errbuf, "sur243 List%%sur242*poideriv"); return(varkon_erpush("SU2943",errbuf)); } #endif /*! */ /* Initialization of point and derivative array SKL */ /* !*/ if ( n_der < deg_u ) min_du = n_der; else min_du = deg_u; if ( n_der < deg_v ) min_dv = n_der; else min_dv = deg_v; for ( k_n = deg_u+1; k_n <= n_der; k_n++) { for ( l_n = 0; l_n <= n_der-k_n; l_n++) { SKL[k_n][l_n][0] = 0.0; SKL[k_n][l_n][1] = 0.0; SKL[k_n][l_n][2] = 0.0; SKL[k_n][l_n][3] = 0.0; } } for ( l_n = deg_v+1; l_n <= n_der; l_n++) { for ( k_n = 0; k_n <= n_der-l_n; k_n++) { SKL[k_n][l_n][0] = 0.0; SKL[k_n][l_n][1] = 0.0; SKL[k_n][l_n][2] = 0.0; SKL[k_n][l_n][3] = 0.0; } } /*! */ /* Calculate point and derivatives */ /* !*/ for ( k_n = 0; k_n <= min_du; k_n++) { for ( s_n = 0; s_n <= deg_v; s_n++) { tempx[s_n] = 0.0; tempy[s_n] = 0.0; tempz[s_n] = 0.0; tempw[s_n] = 0.0; for ( r_n = 0; r_n <= deg_u; r_n++) { /* Get node point */ i_up = span_u - deg_u + r_n; i_vp = span_v - deg_v + s_n; rwcase = 0; status= varkon_poi_nurbs (p_patnu, rwcase, i_up, i_vp, &node); if ( status < 0 ) { sprintf(errbuf, "sur243 Read%%sur242*poideriv"); return(varkon_erpush("SU2943",errbuf)); } polyx = node.x_gm; polyy = node.y_gm; polyz = node.z_gm; polyw = node.w_gm; tempx[s_n] = tempx[s_n] + d_u[k_n][r_n]*polyx; tempy[s_n] = tempy[s_n] + d_u[k_n][r_n]*polyy; tempz[s_n] = tempz[s_n] + d_u[k_n][r_n]*polyz; tempw[s_n] = tempw[s_n] + d_u[k_n][r_n]*polyw; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "i_up %3d i_vp %3d x %8.2f y %8.2f z %8.2f \n", (int)i_up,(int)i_vp, polyx, polyy, polyz); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "s_n %3d k_n %3d r_n %3d temp(s_n)= %8.2f %8.2f %8.2f\n", (int)s_n ,(int)k_n ,(int)r_n, tempx[s_n], tempy[s_n], tempz[s_n] ); fflush(dbgfil(SURPAC)); } #endif } /* End loop r_n */ } /* End loop s_n */ if ( n_der-k_n < min_dv ) min_dd = n_der-k_n; else min_dd = min_dv; for ( l_n = 0; l_n <= min_dd; l_n++) { SKL[k_n][l_n][0] = 0.0; SKL[k_n][l_n][1] = 0.0; SKL[k_n][l_n][2] = 0.0; SKL[k_n][l_n][3] = 0.0; for ( s_n = 0; s_n <= deg_v; s_n++) { SKL[k_n][l_n][0] = SKL[k_n][l_n][0]+d_v[l_n][s_n]*tempx[s_n]; SKL[k_n][l_n][1] = SKL[k_n][l_n][1]+d_v[l_n][s_n]*tempy[s_n]; SKL[k_n][l_n][2] = SKL[k_n][l_n][2]+d_v[l_n][s_n]*tempz[s_n]; SKL[k_n][l_n][3] = SKL[k_n][l_n][3]+d_v[l_n][s_n]*tempw[s_n]; } /* End loop s_n */ } /* End loop l_n */ } /* End loop k_n */ /*! */ /* Coordinates and derivatives to output point p_xyz */ /* !*/ /* Homogenous coordinates */ h_x = SKL[0][0][0]; h_y = SKL[0][0][1]; h_z = SKL[0][0][2]; h_w = SKL[0][0][3]; if ( fabs(h_w) < 100.0*comptol ) { sprintf(errbuf,"Rational denominator = 0%%sur242*poideriv"); return(varkon_erpush("SU2993",errbuf)); } /* Cartesian coordinates */ p_xyz->r_x = h_x/h_w; p_xyz->r_y = h_y/h_w; p_xyz->r_z = h_z/h_w; /* First derivative */ if ( n_der < 1 ) goto all_deriv; hu_x = SKL[1][0][0]; hu_y = SKL[1][0][1]; hu_z = SKL[1][0][2]; hu_w = SKL[1][0][3]; /*! */ /* Cartesian first derivative u_x, u_y, u_z */ /* !*/ p_xyz->u_x = ( h_w*hu_x - hu_w*h_x )/ h_w/ h_w; p_xyz->u_y = ( h_w*hu_y - hu_w*h_y )/ h_w/ h_w; p_xyz->u_z = ( h_w*hu_z - hu_w*h_z )/ h_w/ h_w; hv_x = SKL[0][1][0]; hv_y = SKL[0][1][1]; hv_z = SKL[0][1][2]; hv_w = SKL[0][1][3]; /*! */ /* Cartesian first derivative v_x, v_y, v_z */ /* !*/ p_xyz->v_x = ( h_w*hv_x - hv_w*h_x )/ h_w/ h_w; p_xyz->v_y = ( h_w*hv_y - hv_w*h_y )/ h_w/ h_w; p_xyz->v_z = ( h_w*hv_z - hv_w*h_z )/ h_w/ h_w; /* Second derivatives */ if ( n_der < 2 ) goto all_deriv; hu2_x = SKL[2][0][0]; hu2_y = SKL[2][0][1]; hu2_z = SKL[2][0][2]; hu2_w = SKL[2][0][3]; hv2_x = SKL[0][2][0]; hv2_y = SKL[0][2][1]; hv2_z = SKL[0][2][2]; hv2_w = SKL[0][2][3]; h_d2xdudv = SKL[1][1][0]; h_d2ydudv = SKL[1][1][1]; h_d2zdudv = SKL[1][1][2]; h_d2wdudv = SKL[1][1][3]; huv_x = SKL[1][1][0]; huv_y = SKL[1][1][1]; huv_z = SKL[1][1][2]; huv_w = SKL[1][1][3]; /*! */ /* Cartesian second derivative u2_x, u2_y, u2_z */ /* !*/ p_xyz->u2_x = ( h_w*hu2_x - hu_w*hu_x)/ h_w/ h_w - ((hu2_w*h_x + hu_w*hu_x)* h_w* h_w - 2.0* h_w*h_x *hu_w*hu_w) / h_w/ h_w/ h_w/ h_w; p_xyz->u2_y = ( h_w*hu2_y - hu_w*hu_y)/ h_w/ h_w - ((hu2_w*h_y + hu_w*hu_y)* h_w* h_w - 2.0* h_w*h_y *hu_w*hu_w) / h_w/ h_w/ h_w/ h_w; p_xyz->u2_z = ( h_w*hu2_z - hu_w*hu_z)/ h_w/ h_w - ((hu2_w*h_z + hu_w*hu_z)* h_w* h_w - 2.0* h_w*h_z *hu_w*hu_w) / h_w/ h_w/ h_w/ h_w; /*! */ /* Cartesian second derivative v2_x, v2_y, v2_z */ /* !*/ p_xyz->v2_x = ( h_w*hv2_x - hv_w*hv_x)/ h_w/ h_w - ((hv2_w*h_x + hv_w*hv_x)* h_w* h_w - 2.0* h_w*h_x *hv_w*hv_w) / h_w/ h_w/ h_w/ h_w; p_xyz->v2_y = ( h_w*hv2_y - hv_w*hv_y)/ h_w/ h_w - ((hv2_w*h_y + hv_w*hv_y)* h_w* h_w - 2.0* h_w*h_y *hv_w*hv_w) / h_w/ h_w/ h_w/ h_w; p_xyz->v2_z = ( h_w*hv2_z - hv_w*hv_z)/ h_w/ h_w - ((hv2_w*h_z + hv_w*hv_z)* h_w* h_w - 2.0* h_w*h_z *hv_w*hv_w) / h_w/ h_w/ h_w/ h_w; /*! !! Derivative d2r/dudv= uv_x,uv_y,uv_z !!!! Not OK !!!! !*/ /* uv_x= -1.23456789; uv_y= -1.23456789; uv_z= -1.23456789; Also sur224 must be updated if .. */ /* Add for rational derivatives ............. */ #ifdef REMOVE_980128 p_xyz->uv_x = h_d2xdudv; p_xyz->uv_y = h_d2ydudv; p_xyz->uv_z = h_d2zdudv; #endif /* REMOVE_980128 */ p_xyz->uv_x = (h_w*h_w*(hv_w*hu_x+h_w*huv_x-huv_w*h_x-hu_w*hv_x)+ 2.0*h_w*hv_w*(h_w*hu_x-hu_w*h_x)) / h_w / h_w / h_w / h_w; p_xyz->uv_y = (h_w*h_w*(hv_w*hu_y+h_w*huv_y-huv_w*h_y-hu_w*hv_y)+ 2.0*h_w*hv_w*(h_w*hu_y-hu_w*h_y)) / h_w / h_w / h_w / h_w; p_xyz->uv_z = (h_w*h_w*(hv_w*hu_z+h_w*huv_z-huv_w*h_z-hu_w*hv_z)+ 2.0*h_w*hv_w*(h_w*hu_z-hu_w*h_z)) / h_w / h_w / h_w / h_w; all_deriv:; /* Label: All requested derivatives are calculated */ #ifdef DEBUG for ( l_n = 0; l_n <= n_der; l_n++) { for ( k_n = 0; k_n <= n_der; k_n++) { if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*poideriv SKL[%3d][%3d][0-3]= %8.2f %8.2f %8.2f %5.3f\n", (int)l_n, (int)k_n, SKL[l_n][k_n][0] , SKL[l_n][k_n][1] , SKL[l_n][k_n][2] , SKL[l_n][k_n][3] ); } } } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*poideriv h_w %5.4f hu_w %7.4f hv_w %7.4f hu2_w %7.4f hv2_w %7.4f\n", h_w, hu_w, hv_w, hu2_w, hv2_w); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*poideriv x= %8.2f y= %8.2f z= %8.2f \n", p_xyz->r_x,p_xyz->r_y,p_xyz->r_z); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur242*poideriv drdu %8.2f %8.2f %8.2f \n", p_xyz->u_x,p_xyz->u_y,p_xyz->u_z); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur242*poideriv drdv %8.2f %8.2f %8.2f \n", p_xyz->v_x,p_xyz->v_y,p_xyz->v_z); fflush(dbgfil(SURPAC)); } #endif return(SUCCED); } /* End of function */ /********************************************************************/