This program implements the Euler-Kirchhoff model of a thin rod with circular cross section, subject to external forces derivable from a potential energy function. Though originally designed for insertion into program BEND for use in the winding of superconducting coils, it can also be used to route hoses or optical fibers so as to minimize stress while avoiding obstacles or, in the case of fibers, to minimize optical loss caused by bending (oral remark of Rodger Bossert). Unlike BEND it finds a minimal-energy configuration directly, minimizing the energy from among a sequence of trial configurations under the guidance of one of the many publicly available numerical optimization programs.
Below is a listing of the source:
PROGRAM EKRod
* April 8, 2004
*
* Using its Euler-Kirchhoff mathematical model (see Chapters 18 and 19
* in reference [2]), find the minimal-energy configuration of a thin
* rod. The energy is the sum of the stored elastic strain energy and
* that of an external potential. The rod is constrained at both ends.
* The equations for the centerline are solved by a Direct Method from
* the Calculus of Variations, containing an Indirect Method for solution
* of the equations of twist about that centerline.
*
* At present the energy-evaluation subroutine is restricted to rods with
* a circular cross section.
*
* This program is based on ideas which were developed at Fermi National
* Accelerator Laboratory as further development of program LEAD.
*
* For further information on and development of this program call
* Joe Cook at Fermi National Accelerator Laboratory, (630) 840-2821,
* jmcjmc@fnal.gov, http://tdpc02.fnal.gov/cook/; or at B Associates,
* Inc., (630) 964-2220, jmcook@ieee.org.
IMPLICIT NONE
INTEGER MaxDim, MPMax, LevMax, MaxMax
DOUBLE PRECISION FTOL
PARAMETER (LevMax=3, MaxDim=2+3*(2**LevMax-2),
& MPMax=MaxDim+1, MaxMax=50, FTOL=1.0D-6)
INTEGER I, J, ITER, NDIM, MP, Binom(0:2**LevMax+1),
& Degree, Level, IMax, idum
DOUBLE PRECISION Vertex(MaxDim+1,MaxDim), X(MaxDim), Y(MPMax),
& T0(3), T1(3), Mark0(3), Mark1(3), MinBest,
& Aflex, Atwist, Aelong, L0, Length,
& P(0:2**LevMax+1,3), PBest(0:2**LevMax+1,3),
& C(MaxMax,3), PI, INPROD, NORM, SCALE,
& TEMP1, TEMP2, VECTOR(3), VEC1(3),
& Bezier, UTotal
CHoratiu
DOUBLE PRECISION mass,acc_g
REAL ran1
EXTERNAL Bezier, UTotal
CHARACTER TITLE*72, INNAM*30, OUTNAM*30, Char72*72
& , Char1
LOGICAL DEBUG
COMMON/BEZ/ P, Binom, Level, IMax
COMMON/UTOTL/ Aflex, Atwist, Aelong, L0, C, Length, mass,
& acc_g
COMMON/EndPnt/ T0, T1, Mark0, Mark1
* Glossary
*
* Bezier Curve generator, an external function
* Binom() Binomial coeficients for the Bezier polynomial
* C(i,j) jth coordinate of the ith grid point on the
* centerline of the rod
* Degree Degree of the Bezier polynomial and therefore one
* less than the number of control points, so it
* equals 2**Level + 1.
* I, J Counters
* idum Negative integer to initialize ran1
* L0 Prescribed length of the rod.
* Level Counter for the outer loop, the loop in which
* the number of control points increases.
* LevMax Maximum allowable value for Level
* Mark0,Mark1 Marks at the initial and final positions to show
* where a specified ruling of the cylindrical surface
* of the straight rod is to be placed when it is twisted.
* They do not determine which of the twists, differing
* by 360 degrees, will rotate the final end of the rod
* to line it up with Mark1.
* MaxDim Maximum allowable value for NDIM
* MinBest The lowest value of the objective function currently
* found by the optimizer. It is associated with PBest.
* MP Number of vertices in the simplex
* MPMax Maximum allowed number of vertices in the simplex
* NDIM The number of parameters varied by Nelder-Mead.
* It equals two plus, three times four less than
* the number of control points; which also equals
* 2 + 3*(2**Level - 2).
* P(i,j) jth coordinate of the ith control point, i = 0, ...
* ..., Degree.
* PBest(i,j) The P(i,j) associated with the lowest value, MinBest,
* of the objective function currently found by the optimizer.
* PI Pi.
* ran1 Random number generator for the interval (0,1).
* See page 271 in reference [5].
* SCALE Initial estimate of the size of the convex set
* containing the configured rod.
* TEMPi Temporary scalar variables
* T0() Unit vector tangent to the centerline of the rod,
* pointing into the rod from its initial endpoint.
* T1() Unit vector tangent to the centerline of the rod,
* pointing INTO the rod from its final endpoint.
* UTotal total stored energy, an external function
* VECi() Temporary vector variables
* VECTOR() Temporary vector variable
* Vertex(i,j) jth coordinate of the ith vertex of the simplex
* X(j) jth coordinate of a vertex
* Y() Values of the stored energy at the NDIM+1 vertices
* of the Nelder-Mead simplex
* References
*
* [1] Louis Brand, "Vector and Tensor Analysis", John Wiley & Sons,
* Inc., New York (1947).
* [2] A. E. H. Love, "A Treatise on the Mathematical Theory of
* Elasticity" Dover Publications, New York.
* [3] James L. Kuester and Joe H. Mize, "Optimization Techniques
* with Fortran", McGraw-Hill Book Company, New York (1973).
* [4] W. Murray (editor), "Numerical Methods for Unconstrained
* Optimization", Academic press, New York (1972).
* [5] W. H. Press, B. P. Flannerey, S. A. Teukolsky and W. T.
* Vetterling, Numerical Recipes in FORTRAN, Second Edition,
* Cambridge University Press, New York (1992).
* [6] Gerald Farin, "Curves and Surfaces for Computer Aided Geometric
* Design", 2nd Edition, Academic press, Inc., New York (1990).
* [7] F. S. Hill Jr., "Computer Graphics", Macmillan Publishing Co.,
* New York (1990).
* [8] Su Bu-qing and Liu Ding-yuan, "Computational Geometry, Curve
* and Surface Modeling", Academic Press, Inc., New York (1989).
* [9] I. D. Faux and M. J. Pratt, "Computational Geometry for Design
* and Manufacture", John Wiley & Sons, New York (1979).
* [10] T. v. Kármán and M. Biot, "Mathematical Methods in Engineering",
* McGraw-Hill Book Company, Inc., New York (1940).
* [11] H. Goldstein, "Classical Mechanics", Addison-Wesley Press, Inc.,
* Cambridge (1950).
* [12] J. P. Den Hartog, "Strength of Materials", Dover Publications, Inc.,
* New York (1961).
* [13] H. Cramér, "Mathematical Methods of Statistics", Princeton University
* Press, Princeton (1951).
*
DEBUG = .TRUE.
PI = ACOS(-1.D0)
IMax = 31
idum = -3
WRITE(*,9034)
9034 FORMAT(/' Please enter the complete file name for the disk'
& ,' file from which input is to be'
& /' read. Include the file type'
& /' with separating punctuation ".".')
READ(*,9800) INNAM
OPEN( 35, FILE=INNAM, STATUS='OLD')
READ( 35, 9800) TITLE
WRITE(*,*)
WRITE(*,*) ' The title of the input file is:'
WRITE(*,*)
WRITE(*,*) TITLE
READ( 35, 9800) Char72
READ( 35, 9800) Char72
READ(35,*) P(0,1), P(0,2), P(0,3)
READ(35,*) P(1,1), P(1,2), P(1,3)
READ(35,*) P(2,1), P(2,2), P(2,3)
READ(35,*) P(3,1), P(3,2), P(3,3)
READ(35,*) L0
READ(35,*) Mark0
READ(35,*) Mark1
READ(35,*) Aflex
READ(35,*) Atwist
READ(35,*) Aelong
READ(35,*) mass
READ(35,*) acc_g
CLOSE(35)
TEMP1 = SQRT( (P(0,1)-P(3,1))**2 + (P(0,2)-P(3,2))**2 +
& (P(0,3)-P(3,3))**2)
IF ( TEMP1 .GT. L0*1.001) THEN
WRITE(*,*)
WRITE(*,*)'The proposed initial and final points'
WRITE(*,*)'of the centerline are further apart'
WRITE(*,*)'than its proposed length.'
STOP 'Please resubmit your data.'
ELSE IF ( TEMP1 .GT. L0*.999) THEN
WRITE(*,*)
WRITE(*,*)'The minimum energy configuration is a straight line'
DO I=1,IMax
C(I,1) = P(0,1) + ( P(3,1) - P(0,1))*(I-1)/DBLE(IMax-1)
C(I,2) = P(0,2) + ( P(3,2) - P(0,2))*(I-1)/DBLE(IMax-1)
C(I,3) = P(0,3) + ( P(3,3) - P(0,3))*(I-1)/DBLE(IMax-1)
END DO
* OUTPUT ...
GO TO 9999
END IF
* Compute the unit tangent vectors at the end-points of the
* centerline of the rod.
DO J=1,3
T0(J) = P(1,J) - P(0,J)
T1(J) = P(2,J) - P(3,J)
END DO
TEMP1 = NORM(T0)
IF (TEMP1 .LE. FTOL) THEN
WRITE(*,*)'Your first two points are too close together.'
WRITE(*,*)'Please reenter them.'
* GOTO ...
END IF
TEMP2 = NORM(T1)
IF (TEMP2 .LE. FTOL) THEN
WRITE(*,*)'Your last two points are too close together.'
WRITE(*,*)'Please reenter them.'
* GOTO ...
END IF
T0 = T0/TEMP1
T1 = T1/TEMP2
* Initialize the simplex.
Vertex(1,1) = TEMP1
Vertex(1,2) = TEMP2
Vertex(2,1) = TEMP1/2
Vertex(2,2) = TEMP2
Vertex(3,1) = TEMP1
Vertex(3,2) = TEMP2/2
SCALE = ( TEMP1 + TEMP2)/2
* Orthonormalize Mark0 with respect to T0.
TEMP1 = INPROD( Mark0, T0)
Mark0 = Mark0 - TEMP1 * T0
TEMP1 = NORM( Mark0)
IF ( TEMP1 .LT. 1.D-14) THEN
PRINT*
PRINT*,' INPUT ERROR: The initial mark has no component'
PRINT*,' perpendicular to the initial tangent.'
STOP' '
END IF
Mark0 = Mark0/TEMP1
* Orthonormalize MARK with respect to T1.
TEMP1 = INPROD( Mark1, T1)
Mark1 = Mark1 - TEMP1 * T1
TEMP1 = NORM( Mark1)
IF ( TEMP1 .LT. 1.D-14) THEN
PRINT*
PRINT*,' INPUT ERROR: The final mark has no component'
PRINT*,' perpendicular to the final tangent.'
STOP' '
END IF
Mark1 = Mark1/TEMP1
VECTOR(1) = P(3,1) - P(0,1)
VECTOR(2) = P(3,2) - P(0,2)
VECTOR(3) = P(3,3) - P(0,3)
CALL XPROD( T0, VECTOR, VEC1)
TEMP1 = NORM( VEC1)
CALL XPROD( T1, VECTOR, VEC1)
TEMP2 = NORM( VEC1)
IF ( NORM(VECTOR) .GT. 1.D-14 .AND.
& TEMP1+TEMP2 .LT. 1.D-14 .AND.
& ( ( mass .LT. 1.D-14).OR.
& ( ABS( VECTOR(1))+ABS( VECTOR(2)) .LT. 1.D-14)
& )
& ) THEN
WRITE(*,*)
WRITE(*,*)'Bifurcation: Solutions to the optimization problem,'
WRITE(*,*)'with the given endpoint conditions, exist but are'
WRITE(*,*)'not unique. They should be found by using the'
WRITE(*,*)'analysis of the buckling of columns under axial'
WRITE(*,*)'load (see Section VII.11 in reference [10])'
WRITE(*,*)'instead of using this program. Here LEVEL 1 may'
WRITE(*,*)'output an unstable, highly compressed, straight'
WRITE(*,*)'configuration of the rod. If you wish nonetheless'
WRITE(*,*)'to continue, depress the ENTER key.'
READ*
END IF
* Initialize Binom.
Binom(0) = 1
DO I=1,2**LevMax+1
Binom(I) = 0
END DO
* Initialize the binomial coefficients
* for the cubic Bezier polynomial.
DO J=1,3
DO I=J,1,-1
Binom(I) = Binom(I) + Binom(I-1)
END DO
END DO
***** Enter the first level for a two-parameter minimization of the
* stored energy.
Level = 1
NDIM = 2 + 3*(2**Level - 2)
MP = NDIM + 1
Degree = 2**Level + 1
* Evaluate the stored energy at each vertex of the initial simplex.
DO I=1,MP
DO J=1,NDIM
X(J) = Vertex(I,J)
END DO
Y(I) = Bezier( X, NDIM)
END DO
* Start the optimization by a two-parameter approximation to the
* rod-configuration which minimizes its stored energy.
CALL NELMED( Vertex, Y, MP, MaxDim, NDIM, FTOL, Bezier, ITER)
WRITE(*,*)
WRITE(*,'(1x,a,i5)') 'When the number of parameters was', NDIM
WRITE(*,'(1x,a,i5)') 'the number of iterations was ', iter
WRITE(*,'(1X,A,1X,E12.5)') 'the length of the rod was ', Length
WRITE(*,'(1x,a,2(1X,F12.6))') 'and the minimum was', Y(1)
WRITE(*,*)'Grid points were at'
DO I=1,(IMax+2)/3
WRITE(*,9500) 3*(I-1)+1,
& C(3*(I-1)+1,1), C(3*(I-1)+1,2), C(3*(I-1)+1,3)
END DO
WRITE(*,*)
WRITE(*,*) ' Press ENTER key to continue.'
READ*
IF ( Level .NE. LevMax) THEN
* Invoke Bezier one more time at this level
* to leave the control points in their optimal positions.
DO J=1,NDIM
X(J) = Vertex(1,J)
END DO
TEMP1 = Bezier( X, NDIM)
END IF
MinBest = Y(1)
* CALL PLTCRV( C, IMax)
***** Enter the Level loop.
DO 6000 Level=2,LevMax
NDIM = 2 + 3*(2**Level - 2)
MP = NDIM + 1
Degree = 2**Level + 1
* Advance the binomial coefficients to
* those for the Degreeic Bezier polynomial.
DO J=2**(Level-1)+2,Degree
DO I=J,1,-1
Binom(I) = Binom(I) + Binom(I-1)
END DO
END DO
* Split the middle control points from the previous level
* to form the middle control points for this level.
DO J=1,3
* Move the last control point up in the array
* to make room for the other new control points.
P(Degree,J) = P(2**(Level-1)+1,J)
* Store the next to the last old control point.
VECTOR(J) = P(2**(Level-1),J)
* Move it up the array, and then move it half-way towards
* the last new control point.
P(Degree-1,J) = ( P(Degree,J) + VECTOR(J))/2
* Split the other varying old control points, using I to count
* the old line-segments connecting adjacent pairs of such
* control points.
DO I=2**(Level-1),2,-1
P(Degree - 2*(2**(Level-1) - I + 1), J) =
& ( 3*VECTOR(J) + P(I-1,J))/4
P(Degree - 2*(2**(Level-1) - I + 1) - 1,J) =
& ( VECTOR(J) + 3*P(I-1,J))/4
VECTOR(J) = P(I-1,J)
END DO
* Now move the old P(1,*) half-way towards P(0,*).
P(1,J) = (VECTOR(J) + P(0,J))/2
END DO
PBest = P
* Initialize the simplex as one with its first
* vertex at the best vertex from the previous Level
* by transforming from the control point variables
* to the optimizer variables.
60 CONTINUE ! Target of GOTOs for Level-repeat.
P = PBest
* First take care of the first vertex, Vertex(1,*).
* Set Vertex(1,1) equal to the inner product of
* T0(*) with P(1,*) - P(0,*).
Vertex(1,1) = 0.0D0
DO J=1,3
Vertex(1,1) = Vertex(1,1) + T0(J)*(P(1,J) - P(0,J))
END DO
* Initialize the remaining coordinates of the first vertex
* except for its last.
DO J=1,3
DO I=2,Degree-2
Vertex(1, 1 + 3*(I-2) + J) = P(I,J)
END DO
END DO
* Initialize the last coordinate of the first vertex.
* Set Vertex(1,NDIM) equal to the inner product
* of T1(*) with P(Degree-1,*) - P(Degree,*).
Vertex(1,NDIM) = 0.0D0
DO J=1,3
Vertex(1,NDIM) = Vertex(1,NDIM) +
& T1(J)*(P(Degree-1,J) - P(Degree,J))
END DO
* Now that the first vertex has been taken care of
* initialize the remaining vertices of the simplex.
* The NDIM line segments joining each to the first
* vertex are of random length and are all perpendicular
* to each other.
DO I=2,MP
DO J=1,NDIM
Vertex(I,J) = Vertex(1,J)
END DO
END DO
DO I=2,MP
* Construct an approximately normal probability distribution
* with mean zero and variance one (see page 245 in section
* 19.1 of [13]).
TEMP1 = SQRT(2.E0)*( ran1(idum) + ran1(idum) + ran1(idum) +
& ran1(idum) + ran1(idum) + ran1(idum) - 3)
* Split this distribution away from zero.
IF ( TEMP1 .GE. 0.D0) THEN
TEMP1 = TEMP1 + 1.D0
ELSE
TEMP1 = TEMP1 - 1.D0
END IF
* Scale and add it to the (I-1)th component of the Ith vertex.
Vertex(I,I-1) = Vertex(I,I-1) + TEMP1*SCALE/(4*Level)
END DO
DO I=1,MP
DO J=1,NDIM
X(J) = Vertex(I,J)
END DO
Y(I) = Bezier( X, NDIM)
END DO
CALL NELMED( Vertex, Y, MP, MaxDim, NDIM, FTOL, Bezier, iter)
WRITE(*,*)
WRITE(*,'(1X,A,I5)') 'When the number of parameters was', NDIM
WRITE(*,'(1X,A,I5)') 'the number of iterations was ', iter
WRITE(*,'(1X,A,1X,E12.5)') 'the length of the rod was ', Length
WRITE(*,'(1X,A,2(1X,F12.6))') 'and the minimum was ',
& Y(1)
IF ( MinBest .LT. Y(1)) THEN
WRITE(*,'(1X,A,2(1X,F12.6))') 'whereas the best minimum was',
& MinBest
ELSE
MinBest = Y(1)
PBest = P
END IF
WRITE(*,*)'Grid points were at'
DO I=1,(IMax+2)/3
WRITE(*,9500) 3*(I-1)+1,
& C(3*(I-1)+1,1), C(3*(I-1)+1,2), C(3*(I-1)+1,3)
END DO
WRITE(*,*)
WRITE(*,*) ' Press ENTER key to continue.'
READ*
* Invoke Bezier one more time at this level
* to leave the control points in their optimal positions.
DO J=1,NDIM
X(J) = Vertex(1,J)
END DO
TEMP1 = Bezier( X, NDIM)
* CALL PLTCRV( C, IMax)
PRINT*,' Do you want to try this level again? (N/Y)?'
READ(*,9800) Char1
IF ( Char1.EQ.'Y' .OR. Char1.EQ.'y') GO TO 60
PRINT*,' Do you want to uniformize the control points? (N/Y)?'
READ(*,9800) Char1
IF ( Char1.EQ.'Y' .OR. Char1.EQ.'y') THEN
CALL Uniform( P, Level)
* CALL PLTCRV( C, IMax)
GO TO 60
END IF
IF ( Level .EQ. LevMax) THEN
WRITE(*,9036)
9036 FORMAT(/' Please enter the complete file name for the disk'
& ,' file to which output is to be'
& /' written. Include the file type'
& /' with separating punctuation ".".')
READ(*,9800) OUTNAM
OPEN( 36, FILE=OUTNAM, STATUS='NEW')
WRITE( 36, 9800) TITLE
WRITE( 36, 9550) IMax, ' = IMax'
WRITE( 36, 9800)
&' R(I,1) R(I,2) R(I,3) I'
WRITE( 36, 9800)
&' ------------ ------------ ------------ ---'
DO I=1,IMax
WRITE( 36, 9600) C(I,1), C(I,2), C(I,3), I
END DO
CLOSE(36)
* WRITE(*,*)'Do you want to go to the next point in a parameter v
* &ariation? (Y/N)'
* READ(*,9800) Char1
* IF ( Char1.NE.'N' .AND. Char1.NE.'n') THEN
* CALL HOMOTOP( P, Mark1, L0)
* GO TO 60
* END IF
END IF
6000 CONTINUE
***** Exit from the Level loop.
9999 WRITE(*,*)
STOP ' ***** Normal End of Program EKRod *****'
9500 FORMAT(I5,1X,3(1X,D18.6))
9550 FORMAT(1X,I4,A)
9600 FORMAT(1X,3(1X,D18.6),1X,I5)
9800 FORMAT(A)
END
************************************************************************
SUBROUTINE NELMED( Vertex, y, mp, MaxDim, NDIM, ftol, FUNC, iter)
* Find a local minimum of the given function FUNC using the simplex
* method of Nelder and Mead. See Chapter 9.I, pages 298-308 in
* reference [3], Chapter 2.11, pages 26-28 in [4], or Chapter 10.4,
* pages 402-406 in [5]. This method does not use derivatives of FUNC
* and therefore converges more slowly than those that do but it is of
* wider applicability.
* Glossary
*
* MaxDim Maximum allowable value for NDIM
* MP Number of vertices in the simplex
* NDIM Number of parameters varied by the optimizer.
* Vertex(i,j) jth coordinate of the ith simplex.
* Y(i) Value of the objective function at the vertex
* of the ith simplex.
*
IMPLICIT NONE
INTEGER iter, mp, NDIM, MaxDim, NMAX, ITMAX
DOUBLE PRECISION ftol, Vertex(MaxDim+1,MaxDim), y(mp), FUNC
PARAMETER (NMAX=45, ITMAX=40000)
INTEGER i, ihi, ilo, inhi, j, m,n
EXTERNAL FUNC
DOUBLE PRECISION rtol, sum, swap, ysave, ytry, psum(NMAX), amotry
iter = 0
1 DO n=1,NDIM
sum = 0.
DO m=1,NDIM+1
sum = sum+Vertex(m,n)
END DO
psum(n) = sum
END DO
2 ilo = 1
if (y(1).gt.y(2)) then
ihi = 1
inhi = 2
ELSE
ihi = 2
inhi = 1
END IF
DO i=1,NDIM+1
if(y(i).le.y(ilo)) ilo=i
if(y(i).gt.y(ihi)) then
inhi = ihi
ihi = i
ELSE if(y(i).gt.y(inhi)) then
if(i.ne.ihi) inhi = i
END IF
END DO
* Test to see if NELMED has converged.
rtol = 2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
IF (rtol.lt.ftol) THEN
* NELMED has converged.
swap = y(1)
y(1) = y(ilo)
y(ilo) = swap
DO n=1,NDIM
swap = Vertex(1,n)
Vertex(1,n) = Vertex(ilo,n)
Vertex(ilo,n) = swap
END DO
* Exit to calling program.
RETURN
END IF
IF (iter .GE. ITMAX) THEN
WRITE(*,*) ' ITMAX exceeded in NELMED'
READ*
END IF
iter = iter+2
ytry = amotry(Vertex,y,psum,mp,MaxDim,NDIM,FUNC,ihi,-1.0D0)
if (ytry.le.y(ilo)) then
ytry = amotry(Vertex,y,psum,mp,MaxDim,NDIM,FUNC,ihi,2.0D0)
ELSE if (ytry.ge.y(inhi)) then
ysave = y(ihi)
ytry = amotry(Vertex,y,psum,mp,MaxDim,NDIM,FUNC,ihi,0.5D0)
if (ytry.ge.ysave) then
DO i=1,NDIM+1
if(i.ne.ilo)then
DO j=1,NDIM
psum(j) = 0.5*(Vertex(i,j)+Vertex(ilo,j))
Vertex(i,j) = psum(j)
END DO
y(i) = FUNC( psum, NDIM)
END IF
END DO
iter = iter+NDIM
GOTO 1
END IF
ELSE
iter = iter-1
END IF
GOTO 2
END
************************************************************************
FUNCTION amotry(Vertex, y, psum, mp, MaxDim, NDIM, FUNC, ihi, fac)
INTEGER ihi, mp, NDIM, MaxDim, NMAX
DOUBLE PRECISION amotry, fac, Vertex(MaxDim+1,MaxDim),
& psum(NDIM), y(mp), FUNC
PARAMETER (NMAX=45)
INTEGER j
DOUBLE PRECISION fac1, fac2, ytry, ptry(NMAX)
EXTERNAL FUNC
fac1 = (1.-fac)/NDIM
fac2 = fac1-fac
DO j=1,NDIM
ptry(j) = psum(j)*fac1-Vertex(ihi,j)*fac2
END DO
ytry = FUNC( ptry, NDIM)
if (ytry.lt.y(ihi)) then
y(ihi) = ytry
DO j=1,NDIM
psum(j) = psum(j)-Vertex(ihi,j)+ptry(j)
Vertex(ihi,j) = ptry(j)
END DO
END IF
amotry = ytry
RETURN
END
************************************************************************
DOUBLE PRECISION FUNCTION Bezier( Lambda, NDIM)
* Given the independent variables Lambda, pass the corresponding
* Bezier curve to UTotal for a calculation of the total elastic strain
* and external potential energy of an Euler-Kirchhoff elastica having
* that curve as its centerline.
*
* For information on Bezier curves, see Chapter 3 in reference
* [6]; sections 14.2 and 14.3 in [7]; Chapter 4 in [8]; and Chapter
* 5 in [9].
*
************************************************************************
IMPLICIT NONE
INTEGER NDIM, MaxDim, LevMax, I, II, J, IMax, MaxMax
PARAMETER (LevMax=4, MaxDim=2+3*(2**LevMax-2), MaxMax=50)
INTEGER Binom(0:2**LevMax+1), Degree, Level
DOUBLE PRECISION T, C(3), B(0:2**LevMax+1), R(MaxMax,3), UTotal,
& Lambda( MaxDim), P(0:2**LevMax+1,3), DELTAT,
& T0(3), T1(3), Mark0(3), Mark1(3)
EXTERNAL UTotal
COMMON/BEZ/ P, Binom, Level, IMax
COMMON/EndPnt/ T0, T1, Mark0, Mark1
Degree = 2**Level + 1
* Transform from the optimizer variables
* to the control point variables.
DO J=1,3
P(1,J) = P(0,J) + Lambda(1)*T0(J)
IF ( Lambda(1) .LT. 1.D-15) THEN
Bezier = 1.D20
RETURN
END IF
DO I=2,Degree-2
P(I,J) = Lambda(1 + 3*(I-2) + J)
END DO
P(Degree-1,J) = P(Degree,J) + Lambda(NDIM)*T1(J)
IF ( Lambda(NDIM) .LT. 1.D-15) THEN
Bezier = 1.D30
RETURN
END IF
END DO
DELTAT = 1.D0/(IMax - 1)
DO II=1,IMax
T = (II-1)*DELTAT
* Evaluate the Bernstein polynomials at T.
DO I=0,Degree
B(I) = Binom(I) * (T**I)*((1-T)**(Degree-I))
END DO
* Multiply the Degree+1 vectors P(i,*) representing the control
* points by these Degree+1 scalar functions of T and add them
* together.
C(1) = 0
C(2) = 0
C(3) = 0
DO I=0,Degree
C(1) = C(1) + B(I)*P(I,1)
C(2) = C(2) + B(I)*P(I,2)
C(3) = C(3) + B(I)*P(I,3)
END DO
R(II,1) = C(1)
R(II,2) = C(2)
R(II,3) = C(3)
END DO
Bezier = UTotal( R, IMax)
END
************************************************************************
DOUBLE PRECISION FUNCTION UTotal( R, IMax)
* Compute the length and strain-energy of the curve R as an elastica.
* At present the rod is required to have a circular cross section.
IMPLICIT NONE
INTEGER I, IMax, MaxMax, J
PARAMETER ( MaxMax=50)
DOUBLE PRECISION Length, ENERGY, DSNEW, DSOLD, ANGLE, BOld(3),
& R(MaxMax,3), C(MaxMax,3), XNEW(3), XOld(3),
& T0(3), T1(3), T(3), Twist,
& Aflex, Atwist, Aelong, L0, TauDS, dRotate(3),
& NORM, INPROD, TEMP1, BNew(3), Torse, Beta,
& N1(3), NoTwst(3), Mark0(3), Mark1(3), Pi
CHoratiu
DOUBLE PRECISION mass,acc_g,z
COMMON/UTOTL/ Aflex, Atwist, Aelong, L0, C, Length, mass,
& acc_g
COMMON/EndPnt/ T0, T1, Mark0, Mark1
* Beta Angle between Mark0 and the binormal in the initial
* Frenet frame. It is measured in the sense given by
* the right-hand rule with respect to advancing arc
* length along R.
* BOld() Binormal to the curve at the (I-1)st grid point.
* BNew() Binormal to the curve at the Ith grid point.
* C(i,j) Duplicate of R(i,j) in order to avoid conflict of
* dummy argument with COMMON
* dRotate() Increment of rotation of the Frenet frame.
* Mark0,Mark1 Marks at the initial and final positions to show
* where a specified ruling of the cylindrical surface
* of the straight rod is to be placed when it is twisted.
* They do not determine which of the twists, differing
* by 360 degrees, will rotate the final end of the rod
* to line it up with Mark1.
* NoTwst() When the rod is in its unstrained state its surface
* can be represented by a right-circular cylinder of
* radius a. In order to specify the configuration of
* the rod after it has been twisted, the user picks
* a particular ruling of this cylinder. The initial
* end of this curve-segment is to be placed at R(1,)+
* a*Mark0() and parallel to T0. Its other end is eventually
* to be placed at R(IMax,)+a*Mark1(), after it has been
* twisted about R(,), but first the position must be found
* of other end BEFORE it has been twisted about R(,) so
* that we will have a null position from which to measure
* the final twist. This null position is given by NoTwst.
* It is defined in terms of the of the Frenet normal and
* binormal as basis vectors at R(IMax,) but these vectors
* have rotated about R(I,), introducing an irrelevant twist
* as I advanced so this rotation must be canceled out. The
* rate of change of this irrelevant twist is called "torsion",
* an unfortunate terminology leading to some confusion with
* the twisting of the rod-material which also rotates about
* the centerline. The torsion depends only on the configura-
* tion of the centerline, not on that of the rod-material
* which surrounds it, so it can be computed. Being a rate
* we then integrate it to find the final twist of the Frenet
* frame, which is called here "torse", and which is then
* subtracted to give the configuration of the UNTWISTED
* material. What has been found is a curve-segment whose
* terminal point is R(Imax,)+a*NoTwst() and whose initial
* point is R(1,)+a*Mark0(), and which is everywhere parallel
* to R(,). (See the Example on pages 94-95 of [1] for a
* definition of parallelism of curves. If a plane is passed
* perpendicularly though any point of R(,) it will
* intersect this curve-segment at a point where the
* tangents to the two curves are parallel.)
* NoTwst would make an angle of -torse with the binormal
* at I=IMax if the curve-segment started out at the
* binormal at I=1, but instead it starts out at Mark0 so
* that angle is Beta-torse.
* N1() Frenet normal curvature vector at I=IMax. There the
* triple [XNew,N1,BNew] defines a right-handed
* coordinate system, one in which the binormal vector
* rotates AWAY from N1 when it is advancing with a
* right-handed helical motion.
* Pi Pi
* R(i,j) jth coordinate of the ith grid point on the centerline
* of the rod
* T() Unit tangent vector to the centerline.
* TauDS The increment of torsion
* TEMPi Temporary variables
* T0 Unit tangent vector at the initial point of the
* centerline.
* T1 The NEGATIVE of the unit tangent vector at the final
* point of the centerline.
* Torse Accumulated torsion
* Twist The total twist of the rod from one end to the other.
* Because the rate of twist is a cyclic variable (see
* page 165 in [11]) it is constant. Alternatively,
* because this variable separates out in the Hamiltonian,
* the corresponding energy can be minimized independently
* of the that of the others, and the solution of its
* Euler-Lagrange equations is a constant.
Pi = ACOS(-1.D0)
* I=1:
Torse = 0
XOld = T0
XNEW(1) = R(2,1) - R(1,1)
XNEW(2) = R(2,2) - R(1,2)
XNEW(3) = R(2,3) - R(1,3)
z = (R(1,3) + R(2,3))/2.d0
DSNEW = NORM(XNEW)
IF ( DSNEW .LT. 1.D-14) THEN
WRITE(*,*)
WRITE(*,*)'The first two grid points on the curve'
WRITE(*,*)'are coincident.'
* Return a large energy so that it will be rejected
* by the minimizer.
UTotal = 1.D20
RETURN
END IF
T = XNEW/DSNEW
CALL XPROD( XOld, XNew, BOld)
TEMP1 = NORM( BOld)
IF ( TEMP1 .LT. 1.D-14) THEN
WRITE(*,*)
WRITE(*,*)'The first two grid points on the curve'
WRITE(*,*)'are parallel to the initial tangent.'
STOP' '
END IF
BOld = BOld/TEMP1
CALL XPROD( BOld, Mark0, dRotate)
Beta = ASIN( INPROD( T0, dRotate))
IF ( INPROD( Mark0, BOld) .LT. 0.D0) THEN
IF ( Beta .GE. 0.D0 ) THEN
Beta = + Pi - Beta
ELSE
Beta = - Pi - Beta
END IF
END IF
Length = DSNEW
ANGLE = ACOS( MIN( INPROD( T0, XNew)/DSNEW, 1.D0)) / 2
CHoratiu added the gravitational energy for the first segment:
ENERGY = .5 * Aflex * ANGLE**2 / DSNEW - mass*length*acc_g*z
* The movement to the second grid point has been completed.
DO 1000 I=2,IMax-1
DSOLD = DSNEW
XOLD = XNEW
XNEW(1) = R(I+1,1) - R(I,1)
XNEW(2) = R(I+1,2) - R(I,2)
XNEW(3) = R(I+1,3) - R(I,3)
z = (R(I+1,3) + R(I,3))/2.d0
DSNEW = NORM(XNEW)
Length = Length + DSNEW
ANGLE = ACOS(MIN( INPROD( XOLD, XNEW) / (DSOLD*DSNEW), 1.D0))
CHoratiu added the gravitational energy for each segment:
ENERGY = ENERGY + Aflex * ANGLE**2 / ( DSOLD + DSNEW)
& - mass*dsnew*acc_g*z
CALL XPROD( XOld, XNEW, BNew)
TEMP1 = NORM( BNew)
IF ( TEMP1 .LT. 1.D-14) THEN
WRITE(*,*)
WRITE(*,*)'Colinear grid points at', I
GO TO 900
END IF
BNew = BNew/TEMP1
CALL XPROD( BOld, BNew, dRotate)
TauDS = ASIN( INPROD( T, dRotate))
IF ( INPROD( BOld, BNew) .LT. 0.D0) THEN
IF ( TauDS .GE. 0.D0 ) THEN
TauDS = + Pi - TauDS
ELSE
TauDS = - Pi - TauDS
END IF
END IF
Torse = Torse + TauDS
* Prepare for the next grid point.
BOld = BNew
* CALL XPROD( BOld, XNew, N1)
T = XNEW/DSNEW
900 CONTINUE
1000 CONTINUE
* Store R-values in C in order to avoid conflict of dummy argument
* with COMMON.
DO I=1,IMax
DO J=1,3
C(I,J) = R(I,J)
END DO
END DO
* I = IMax:
ANGLE = ACOS( MIN( -INPROD( T1, XNEW) / DSNEW, 1.D0)) / 2
CALL XPROD( T1, XNEW, BNew) ! Remember: T1 points backwards.
TEMP1 = NORM( BNew)
IF ( TEMP1 .LT. 1.D-14) THEN
WRITE(*,*)
WRITE(*,*)'The last two grid points on the curve'
WRITE(*,*)'are parallel to the final tangent.'
STOP' '
END IF
BNew = BNew/TEMP1
* BNew is now the final (i.e., at I=IMax) binormal.
CALL XPROD( BOld, BNew, dRotate)
TauDS = ASIN( INPROD( T, dRotate))
IF ( INPROD( BOld, BNew) .LT. 0.D0) THEN
IF ( TauDS .GE. 0.D0 ) THEN
TauDS = + Pi - TauDS
ELSE
TauDS = - Pi - TauDS
END IF
END IF
Torse = Torse + TauDS
CALL XPROD( T1, BNew, N1) ! Remember: T1 points backwards
TEMP1 = NORM( N1)
IF ( TEMP1 .LT. 1.D-14) THEN
WRITE(*,*)
WRITE(*,*)'The curve makes a 90 degree turn here.'
* Return a large energy so that it will be rejected
* by the minimizer.
UTotal = 1.D20
RETURN
END IF
N1 = N1/TEMP1
* Rotate BNew to remove the accumulated torsion, and then to
* put the origin of the untwisted coordinate system at Mark0.
* (Remember: A right-handed rotation rotates BNew towards -N1.)
NoTwst = COS(Beta-Torse)*BNew - SIN(Beta-Torse)*N1
CALL XPROD( NoTwst, Mark1, dRotate)
Twist = ASIN( -INPROD( T1, dRotate))
IF ( INPROD( NoTwst, Mark1) .LT. 0.D0) THEN
IF ( Twist .GE. 0.D0 ) THEN
Twist = + Pi - Twist
ELSE
Twist = - Pi - Twist
END IF
END IF
IF ( Length .GT. 1.1 * L0) THEN
UTotal = 1.D20 * ( Length - 1.1 * L0)
RETURN
END IF
UTotal = ENERGY + .5 * Aflex * ANGLE**2 / DSNEW
& + Aelong * (Length - L0)**2 / 2
& + Atwist * Twist**2 / ( 2 * Length)
* See the bottom of page 164 in [12] for the reasons why these three
* contributions to the total energy are present as separate summands.
CHoratiu added the gravitational energy for the last segment
& - mass*dsnew*acc_g*z
9500 FORMAT( I5, 1X, 3( 4X, E9.2))
END
************************************************************************
DOUBLE PRECISION FUNCTION INPROD(X,Y)
* Compute the inner product of X and Y.
DOUBLE PRECISION X(3), Y(3)
INPROD = X(1) * Y(1) + X(2) * Y(2) + X(3) * Y(3)
END
************************************************************************
DOUBLE PRECISION FUNCTION NORM(X)
* Compute the length of X.
DOUBLE PRECISION X(3)
NORM = SQRT( X(1)**2 + X(2)**2 + X(3)**2)
END
************************************************************************
SUBROUTINE XPROD( X1, X2, X)
* Compute the cross product X of the two vectors X1 and X2.
DOUBLE PRECISION X1(3), X2(3), X(3)
X(1) = X1(2) * X2(3) - X1(3) * X2(2)
X(2) = X1(3) * X2(1) - X1(1) * X2(3)
X(3) = X1(1) * X2(2) - X1(2) * X2(1)
END
************************************************************************
SUBROUTINE Uniform( P, Level)
* February 4, 2004
*
* Report on the degree of uniformity of a distribution of a sequence of
* Degree+1 points P in Euclidean 3-space, along the broken line segments
* connecting them in order. (A uniform sequence of control points is defined
* to be one which is equidistributed along a line obtained by straightening
* out the broken line.)
*
* At the user's option, remedy the nonuniformity.
IMPLICIT NONE
INTEGER LevMax, Level, Degree, I, J, K
PARAMETER (LevMax=4)
DOUBLE PRECISION P(0:2**LevMax+1,3), P_Old(0:2**LevMax+1,3),
& Vector(3), LengthP(2**LevMax), Mean, Norm,
& DeltaP(2**LevMax-1), Skewness, DeltaUniform
* & ,R2(2**LevMax+2,2)
CHARACTER Char1*1
* Glossary
*
* Degree One less than the number of control points.
* DeltaUniform Distance between each of the Degree-1 new intermediate
* control points.
* I,J,K Counters.
* LengthP(i) Distance from P(1,) to P(i,), as measured along the
* sequence of line segments connecting these two points.
* Mean The average of all of the LengthP(i)s.
* P(i,j) The jth component of the ith control point, counting
* from i=0.
* P_Old(i,j) Copy of the input P(i,j).
* Skewness An indication of the deviation of the mean of the
* Degree-1 intermediate control points, as measured
* along the broken line segment connecting them. It
* is normalized to lie in the interval [-1.+1].
* Vector() Dummy vector.
Degree = 2**Level + 1
* Compute the lengths DeltaP(I) of each of the straight line segments
* connecting, in sequence, the Degree-1 intermediate control points,
* and also their cummulative sums, LengthP(I).
LengthP(1) = 0
DO J=1,3
Vector(J) = P(2,J) - P(1,J)
END DO
DeltaP(1) = Norm( Vector)
LengthP(2) = DeltaP(1)
DO I=2,Degree-2
DO J=1,3
Vector(J) = P(I+1,J) - P(I,J)
END DO
DeltaP(I) = Norm( Vector)
LengthP(I+1) = LengthP(I) + DeltaP(I)
END DO
* Compute the mean value of all of the LengthP(i)s.
Mean = 0
DO I=2,Degree-2
Mean = Mean + LengthP(I)
END DO
Mean = Mean/(Degree - 3)
* Scale to the interval [-1,+1].
Skewness = 2 * Mean/LengthP(Degree-1) - 1
* IF ( ABS( Mean/LengthP(Degree)) .LT. .2D0 ) RETURN
PRINT*
PRINT*,' The skewness of the sequence of control points, as measur
&ed'
PRINT*,' along the broken line connecting them and normalized to t
&he'
WRITE(*,9810)'interval between -1 and +1, is', Skewness, '.'
PRINT*,' Do you want to uniformize these control points? (Y/N)'
READ(*,9800) Char1
IF ( Char1.EQ.'N' .OR. Char1.EQ.'n') RETURN
* ELSE uniformize.
DeltaUniform = LengthP(Degree-1)/(Degree-2)
* Store the input control points to make room for the uniformized
* control points.
DO I=0,Degree
DO J=1,3
P_Old(I,J) = P(I,J)
END DO
END DO
* Arrange the new control points uniformly along the same broken line
* that contained the old control points.
I = 1
DO K=1,Degree-3
100 I = I + 1
IF ( K * DeltaUniform .LE. LengthP(I)) THEN
DO J=1,3
P(K+1,J) = P_Old(I-1,J) + ( P_Old(I,J) - P_Old(I-1,J)) *
& ( K * DeltaUniform - LengthP(I-1)) / DeltaP(I-1)
END DO
ELSE
GO TO 100
END IF
I = I - 1
END DO
9800 FORMAT(A)
9810 Format(2X,A,F6.3,A)
END
************************************************************************
FUNCTION ran1(idum)
* Extracted from the "Numerical Recipes" CDROM. See page 271 of
* "Numerical Recipes in FORTRAN", Second Edition, by W. H. Press,
* S. A. Teukolsky, W. T. Vetterling and Brian P. Flannery, Cambridge
* University Press, Cambridge, (1992).
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
REAL ran1,AM,EPS,RNMX
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
*NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER j,k,iv(NTAB),iy
SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
if (idum.le.0.or.iy.eq.0) then
idum=max(-idum,1)
do j=NTAB+8,1,-1
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
if (j.le.NTAB) iv(j)=idum
end do
iy=iv(1)
endif
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
j=1+iy/NDIV
iy=iv(j)
iv(j)=idum
ran1=min(AM*iy,RNMX)
return
END
***********************************************************************
Legal Notices