13472J/1.128J/2158J/16940J COMPUTATIONAL GEOMETRY Lecture 20 Pro R: Ki. Pat ricalakis Copyright( 2003 Massachusetts Institute of Technology Contents 20 Advanced topics in differential geometry eslcs 20.1.1 Motivation 20.1.2 Definition 20.1.3 Governing equations 20.1.4 Two-point boundary value problem 20.1.5 Example 20.2 Developable surface 20.2.1 Motivation 22223580002 20.2.2 Definition 20.2.3 Developable surface in terms of Bezier surface 20.2.4 Development of developable surface(flattening) 20.3 Umbilics 20.3.1 Motivation 20.3.2 Definition 20.3.3 Computation of umbilical points 20.3. 4 Classification 20.3.5 Characteristic lines 20.4 Parabolic, ridge and sub-parabolic points 20.4.2 Focal surfaces 20.4.3 Parabolic points 22 20.4. 4 Ridge points 0.4.5 Sub-parabolic points Bibliography Reading in the Textbook Geodesics: Chapter 10, pp. 265-291 Umbilics: Chapter 9, pp. 231-264
13.472J/1.128J/2.158J/16.940J COMPUTATIONAL GEOMETRY Lecture 20 Dr. K. H. Ko Prof. N. M. Patrikalakis Copyright ≤c 2003 Massachusetts Institute of Technology Contents 20 Advanced topics in differential geometry 2 20.1 Geodesics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 20.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 20.1.2 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 20.1.3 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 20.1.4 Two-point boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . 5 20.1.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 20.2 Developable surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 20.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 20.2.2 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 20.2.3 Developable surface in terms of B´ezier surface . . . . . . . . . . . . . . . . . . . 12 20.2.4 Development of developable surface (flattening) . . . . . . . . . . . . . . . . . . 13 20.3 Umbilics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.2 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.3 Computation of umbilical points . . . . . . . . . . . . . . . . . . . . . . . . . . 15 20.3.4 Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 20.3.5 Characteristic lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 20.4 Parabolic, ridge and sub-parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . 21 20.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 20.4.2 Focal surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 20.4.3 Parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 20.4.4 Ridge points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 20.4.5 Sub-parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Bibliography 25 Reading in the Textbook • Geodesics : Chapter 10, pp.265 - 291 • Umbilics : Chapter 9, pp.231 - 264 1
Lecture 20 Advanced topics in differential geometry 20.1 Geodesics In this section we study the computation of shortest path between two points on free-form surfaces [14, 11 20.1.1 Motivation ship design robot motion planning · terrain navigation installation of underwater cable 20.1.2 Definition . t: Unit tangent vector of C at P n: Unit normal vector of c at p N: Unit surface normal vector of s at P u: Unit vector perpendicular to t in the tangent plane defined by n x t
Lecture 20 Advanced topics in differential geometry 20.1 Geodesics In this section we study the computation of shortest path between two points on free-form surfaces [14, 11]. 20.1.1 Motivation • ship design • robot motion planning • terrain navigation • installation of underwater cable 20.1.2 Definition • t: Unit tangent vector of C at P • n: Unit normal vector of C at P • N: Unit surface normal vector of S at P • u: Unit vector perpendicular to t in the tangent plane defined by N × t. 2
Figure 20.1: Definition of geodesic curvature We can decompose the curvature vector k of C into n component kn, which is called normal curvature vector, and u component ka, which is called geodesic curvature vector k=kn+ko=-knn+kou 20.1 Here kn and Kg are the normal and geodesic curvatures, respectively and defined as follows (20.3) · Consequently, ds Geodesic paths are sometimes defined as shortest path between points on a surface however this is not always a satisfactory definition Definition: Geodesics are curves of zero geodesic curvature 24 20.1.3 Governing equations The unit tangent vector of the curve C on the surface r is given by dr(u(s), v(s)) du ds Hence using chain rules ds =ruu)+2rw'ds ds +Iw(ds)+ruds2krd-u dt dudu du
t N u n k kg kn S P C Figure 20.1: Definition of geodesic curvature. • We can decompose the curvature vector k of C into N component kn, which is called normal curvature vector, and u component kg, which is called geodesic curvature vector k = kn + kg = −ρnN + ρgu (20.1) Here ρn and ρg are the normal and geodesic curvatures, respectively and defined as follows: ρn = −k N· (20.2) ρg = k u· (20.3) • Consequently, dt ρg = ds · (N × t) (20.4) • Geodesic paths are sometimes defined as shortest path between points on a surface, however this is not always a satisfactory definition. Definition: Geodesics are curves of zero geodesic curvature [24]. 20.1.3 Governing equations • The unit tangent vector of the curve C on the surface r is given by dr(u(s), v(s)) du dv t = = ru + rv (20.5) ds ds ds • Hence using chain rules ⎤ ⎦2 ⎤ ⎦2 dt du du dv dv d2u d2v = ruu + 2ruv + rvv + ru + rv (20.6) ds ds ds ds ds ds2 ds2 3
Consider that the surface normal n has the direction of normal of the geodesic line tn Since kn= g, equation(20. 7)can be rewritten as dt dt r,= (20.8) By substituting equation(20.6)into equations(20.8) we obtain du d 2u ds d +(re·ra ds2+ha2=0(20.9) dudu rua‘r ds )+2(rwu ro)ds ds +(roo r s)+ 0(20.10 By eliminating d! from equation(20.)using equation(2010), and eliminating d! from equation(20.10)using equation(20.9)and employing the Christoffel symbols, we obtain du du ds2 dsds (20.11) ds2 +Ii 12 dsds ds (20.12) Where Tik(i, 3, k= 1, 2)are the Christoffel symbols defined as follows TI GE-BC),l≈2EF1-EE+FE GE -2FF, +FE 2(E EG.-FE 2(EG-F2) 20.13) 2G,-GGu+ FG EGn-2FF+FG 2(EG-F2) 2(EG-F2) These two second order differential equations can be rewritten as a system of four first order differential equations 6 20.14 d q 中一山一 (20.16) (20.17
• Consider that the surface normal N has the direction of normal of the geodesic line ±n n r· u = 0, n r· v = 0. (20.7) dt • Since kn = ds , equation (20.7) can be rewritten as dt dt ru = 0, rv = 0 (20.8) ds · ds · • By substituting equation (20.6) into equations (20.8) we obtain ⎤ ⎦2 ⎤ ⎦2 du du dv dv d2u d2v (ruu · ru) + 2(ruv · ru) + (rvv · ru) + E + F = 0 (20.9) ds ds ds ⎤ ⎦2 du du dv ds ds2 ds2 ⎤ ⎦2 dv d2u d2v (ruu · rv) + 2(ruv · rv) + (rvv · rv) + F + G = 0 (20.10) ds ds ds ds ds2 ds2 v By eliminating d2 from equation (20.9) using equation (20.10), and eliminating d2u • ds2 ds2 from equation (20.10) using equation (20.9) and employing the Christoffel symbols, we obtain ⎤ ⎤ ⎦2 d2u du⎦2 du dv dv + �1 + 2�1 + �1 = 0 (20.11) ds2 11 22 ds ⎤ d2v du⎦2 12 ds ds ds ⎤ ⎦2 du dv dv + �2 + 2�2 + �2 = 0 (20.12) ds2 11 22 ds 12 ds ds ds • jk Where � (i, j, k = 1, 2) are the Christoffel symbols defined as follows: i �1 GEu − 2FFu + FEv �2 2EFu − EEv + FEu 11 = 11 = � 2(EG − F2) , 2(EG − F2) 1 GEv − FGu �2 EGu − FEv 12 = 12 = � 2(EG − F2) , 2(EG − F2) (20.13) 1 2GFv − GGu + FGv �2 EGv − 2FFv + FGu 22 = 22 = 2(EG − F2) , 2(EG − F2) • These two second order differential equations can be rewritten as a system of four first order differential equations [6]. du = p (20.14) ds dv = q (20.15) ds dp 2 2 = −�1 − 2�1 11p 12pq − �1 (20.16) 22q ds dq 2 2 = −�2 − 2�2 11p 12pq − �2 (20.17) 22q ds 4
Euler Lagrange Equations(Calculus of Variations We can also find this result by means of the general rules of the calculus of variations We want to minimize d f(u, v, i)du (20.18) f(u,u,d)=VE+2Fi+Gi2, Udu 20.1 This leads to the condition af d af 0 (20.20 from which we can derive the same differential equation for geodesics 20. 1.4 Two-point boundary value problem We can solve a system of four first order ordinary differential equations(20. 14)to(20.17) Initial-value problem(IVP), where all four boundary conditions are given at one point, or as Boundary-value problem(BVP), where four boundary conditions are specified at two distinct points The first order differential equation for a boundary value problem can be written in vector form as where, ds=g(s, y), y(A)=(uA, UA, PA, qA), y(B)=(uB,UB, PB, qB)? (20.21) where pA, PB, gA and gB are unknowns (u,0,P,q) g=(P,q,-Ip2-212p-r292,-I1p2-212p-r22)2(20.23) There are two commonly used approaches to the numerical solution of BVP Shooting method: easy to implement but unstable 2. Relaxation method: more sophisticated but stable · Shooting method We assume values at s= A, which are not given as boundary conditions at s= A and compute the solution of the resulting iVp to s=B
• Euler Lagrange Equations (Calculus of Variations) We can also find this result by means of the general rules of the calculus of variations. We want to minimize � b � b I = ds = f(u, v, v˙)du (20.18) a a where dv 2 f(u, v, v˙) = � E + 2F v˙ + Gv˙ , v˙ = (20.19) du This leads to the condition πf d πf (20.20) πv − du πv˙ = 0 from which we can derive the same differential equation for geodesics. 20.1.4 Two-point boundary value problem • We can solve a system of four first order ordinary differential equations (20.14) to (20.17) as – Initial-value problem (IVP), where all four boundary conditions are given at one point, or as – Boundary-value problem (BVP), where four boundary conditions are specified at two distinct points. • The first order differential equation for a boundary value problem can be written in vector form as: dy = g(s, y), y(A) = (uA, vA, pA, qA) T , y(B) = (uB, vB, pB, qB) T (20.21) ds where pA, pB, qA and qB are unknowns, y = (u, v, p, q) T (20.22) 2 12pq − �1 2 2 22q2 ) T g = (p, q, −�1 − 2�1 22q , −�2 − 2�2 12pq − �2 (20.23) 11p 11p • There are two commonly used approaches to the numerical solution of BVP. 1. Shooting method: easy to implement but unstable. 2. Relaxation method: more sophisticated but stable. • Shooting method – We assume values at s = A, which are not given as boundary conditions at s = A and compute the solution of the resulting IVP to s = B. 5
The computed values of y(B)will not, in general, agree with the corresponding boundary condition at s= B Consequently, we need to adjust the initial values and try again The process is repeated until the computed values at the final point agree with the boundary conditions and referred as shooting method Formulation: Using the first fundamental form, given pa we can obtain gA from FpA±VFmA-G(EnA-1) Thus we assume pA and solve the differential equation as IvP using, say runge- Kutta method. Here we also have to assume the entire arc length of the geodesic path s to stop the integration. Thus the unknowns can be considered as pA and s If we denote the computed value of(uB, UB) as(uB, UB), the difference can be given as (ub -uB, UB We need to adjust pa and s to make the difference zero This can be done by employing the Newtons method B n+1 dub uB(PA+△pA)-ubB(PA) dus uB(s+△s)-uB(s) ap B UR(PA+ ApA)-uR(pa) auk u*(s+As)-uB(s) e Relaxation method The second method is based on a finite difference approximation to a on a mesh of points in the interval [A, B This method starts with an initial guess and improves the solution iteratively and referred as. direct method. relaxation method or finite difference method The shooting method is often very sensitive to the unknown initial angles at point A and unless a good initial guess is provided, the integrated path will never reach the other point B, while the relaxation method starts with two end points fixed and relaxes to the true solution and hence it is much more stable Let us consider a mesh of points satisfying A=S1<S2 <sm=B. We approximate the n first order differential equations by the trapezoidal rule [8 -11 =Gk+Gk-l],k=2,3 (20.24)
� � � (�) (�) – The computed values of y(B) will not, in general, agree with the corresponding boundary condition at s = B. – Consequently, we need to adjust the initial values and try again. – The process is repeated until the computed values at the final point agree with the boundary conditions and referred as shooting method. – Formulation: Using the first fundamental form, given pA we can obtain qA from − 2 F pA ± F2pA − G(Ep2 q A − 1) A = . G Thus we assume pA and solve the differential equation as IVP using, say RungeKutta method. Here we also have to assume the entire arc length of the geodesic path s to stop the integration. Thus the unknowns can be considered as pA and s. � B), the difference can be given � If we denote the computed value of (uB, vB) as (uB, v as (u This can be done by employing the Newton’s method � T B − vB) � B − uB, v . We need to adjust pA and s to make the difference zero. ⎣ � �−1 B �s � B ⎤ ⎦ ⎤ ⎦ �u �u ⎤ ⎦ pA pA �pA uB − uB = − � B �s � B �v �v s n+1 s �p vB − vB n A n where � B � B(pA) � B(pA + �pA) − u � B � B(s) � πu u πu uB(s + �s) − u = , = πPA �pA πs �s � B v � B(pA) � B(pA + �pA) − v � B � B(s) � πv πv vB(s + �s) − v = , = πPA �pA πs �s • Relaxation method dy – The second method is based on a finite difference approximation to ds on a mesh of points in the interval [A, B]. – This method starts with an initial guess and improves the solution iteratively and referred as, direct method, relaxation method or finite difference method. – The shooting method is often very sensitive to the unknown initial angles at point A and unless a good initial guess is provided, the integrated path will never reach the other point B, while the relaxation method starts with two end points fixed and relaxes to the true solution and hence it is much more stable. – Let us consider a mesh of points satisfying A = s1 < s2 < . . . < sm = B. We approximate the n first order differential equations by the trapezoidal rule [8]. Yk − Yk−1 1 = [Gk + Gk−1], k = 2, 3, . . . , m (20.24) sk − sk−1 2 6
where the n-vectors Yk, Gk are meant to approximate y(sk)and g(sk) with bound ary conditions (uA, UA, P1, q1), Ym=B=(u Yi has n1=2 known components, while Ym has n2=n-n1=4-2=2 known components Boundary Conditions Unknown S1 4 SM-1 SM Figure 20.2: Mesh points This discrete approximation will be accurate to the order of h2(h= maxk sk Sk-1)). Equation(20.24)forms a system of(m-1)n nonlinear algebraic equations with(m-1)n unknowns Let us refer to equation(20. 24)as FR=(FLk, F2 F and the boundary conditions(20.25)as F1=(F1 F m+1 (F1 12.n Ym-B then we h F
������������������������������������� where the n-vectors Yk, Gk are meant to approximate y(sk) and g(sk) with boundary conditions Y1 = β = (uA, vA, p1, q1) T , Ym = � = (uB, vB, pm, qm) T (20.25) – Y1 has n1 = 2 known components, while Ym has n2 = n − n1 = 4 − 2 = 2 known components. �������������������������������������✥ u1 u2 u3 u4 uM−2 uM−1 uM �������������������������������������✥ Boundary Conditions Unknowns ������������������������������������� v1 v2 v3 v4 vM−2 vM−1 vM �������������������������������������✥ ������������������������������������� p1 p2 p3 p4 pM−2 pM−1 pM �������������������������������������✥ ������������������������������������� q1 q2 q3 q4 . . . . . . . . . qM−2 qM−1 qM s1 s2 s3 s4 sM−2 sM−1 sM Figure 20.2: Mesh points. s – This discrete approximation will be accurate to the order of h2 (h = maxk{sk − k−1}). Equation (20.24) forms a system of (m − 1)n nonlinear algebraic equations with (m − 1)n unknowns. – Let us refer to equation (20.24) as 1 Fk = (F1,k, F2,k, . . . , Fn,k) T = Yk − Yk−1 [Gk + Gk−1] = 0 (20.26) sk − sk−1 − 2 – and the boundary conditions (20.25) as F F1 = (F1,1, F2,1, . . . , Fn1,1) T = Y1 − β = 0, (20.27) m+1 = (F1,m+1, F2,m+1, . . . , Fn2,m+1) T = Ym − � = 0 – then we have mn nonlinear algebraic equations F = (FT 2 , . . . , FT 1 , FT m+1) T = 0 (20.28) 7
oints Tolerance Iterations Geodesic Distance 0E310110|1:6611.8651.661 Table 20.1: Numerical conditions and results for the computation of the geodesic path between corner points of the wave-like surface and can be computed by quadratically convergent Newton iteration, if a sufficiently accurate starting vector Y(O)=(YT, Y2, .. Ym) is provided. The Newton iter- ation scheme is given by (+1) Y+△Y Jo]△Y0 (20.30) where jJo] is the mn by mn Jacobian matrix of FO) with respect to Yo Figure 20.3: Geodesic paths on the wave-like bicubic B-spline surface between points of two corners 20.1.5 Example Bilinear surface r(u, v)=(u, v, uv) Eu=0, F E= 2 G2=0
Points Tolerance Iterations Geodesic Distance m κN L M R L M R 101 1.0E-3 10 1 10 1.661 1.865 1.661 Table 20.1: Numerical conditions and results for the computation of the geodesic path between corner points of the wave-like surface. – and can be computed by quadratically convergent Newton iteration, if a sufficiently accurate starting vector Y(0) = (YT 2 , . . . , YT 1 , YT m)T is provided. The Newton iteration scheme is given by Y(i+1) [J = Y(i) + �Y(i) (20.29) (i) ]�Y(i) = −F(i) (20.30) where [J(i) ] is the mn by mn Jacobian matrix of F(i) with respect to Y(i) . x y z Figure 20.3: Geodesic paths on the wave-like bicubic B-spline surface between points of two corners. 20.1.5 Example Bilinear surface r(u, v) = (u, v, uv). E E E = 1 + v2, F = uv, G = 1 + u2 u = 0, Fu = v, Gu = 2u v = 2v, Fv = u, Gv = 0 8
科人 Figure 20.4: Convergence of the right geodesic path in Figure 20.3 2()v+ut:(2 2(1+02)(1+2)-2=0 2 u2+u2+1 Finally the differential equations are given by 山一凼 d山 ds +u2+1 dq
x y z Figure 20.4: Convergence of the right geodesic path in Figure 20.3. �1 −2(uv) · v + uv · (2v) = = 0 11 2 2[(1 + v2)(1 + u2) − u v2] �2 = �1 22 = �22 = 0 2 � 11 1 v 12 = � u2 + v2 + 1 2 u 12 = u2 + v2 + 1 Finally the differential equations are given by du = p ds dv = q ds dp −2u = ds u2 + v2 + 1 pq dq −2v = ds u2 + v2 + 1 pq. 9
20.2 Developable surface Developable surfaces are a special class of surfaces that can be developed or unfolded onto a plane without stretching or tearing 1, 9, 5, 16, 15 20.2.1 Motivation In shipbuilding, doubly curved plates are manufactured by roller and line heating pro- cesses, while the singly curved plates(developable surface) are manufactured by roller The use of developable surfaces has several advantages such as lower labor costs in construction, smaller capital investment in equipment, ease of repair and simple tools for construction In automobile production, body panels, upholstery and window glass are developable surfaces 20.2.2 Definition a ruled surface is defined as a surface generated by the motion of a straight line referred as a generator or ruling 24] The mathematical representation of a ruled surface is given by r(u, )=rA(u)+vo(u) 0.3 where r(a) is a 3d curve called the directrix or base curve of the ruled surface and a u) is a unit vector which gives the direction of the ruling at each point on the directrix see Figure 20.5 r(u, v)=ra(u)+ va(u) ruling or directrix Figure 20.5: Definition of ruled surface 10
20.2 Developable surface Developable surfaces are a special class of surfaces that can be developed or unfolded onto a plane without stretching or tearing [1, 9, 5, 16, 15]. 20.2.1 Motivation • In shipbuilding, doubly curved plates are manufactured by roller and line heating processes, while the singly curved plates (developable surface) are manufactured by roller only. • The use of developable surfaces has several advantages such as lower labor costs in construction, smaller capital investment in equipment, ease of repair and simple tools for construction. • In automobile production, body panels, upholstery and window glass are developable surfaces. 20.2.2 Definition • A ruled surface is defined as a surface generated by the motion of a straight line referred as a generator or ruling [24]. • The mathematical representation of a ruled surface is given by r(u, v) = rA(u) + vβ(u) (20.31) where r(u) is a 3D curve called the directrix or base curve of the ruled surface and β(u) is a unit vector which gives the direction of the ruling at each point on the directrix see Figure 20.5. ruling or generator directrix α(u) r(u,v)= + rA(u) rA(u) vα(u) Figure 20.5: Definition of ruled surface 10