diff options
Diffstat (limited to 'Dragon/src/NXTITA.f')
| -rw-r--r-- | Dragon/src/NXTITA.f | 308 |
1 files changed, 308 insertions, 0 deletions
diff --git a/Dragon/src/NXTITA.f b/Dragon/src/NXTITA.f new file mode 100644 index 0000000..235ef98 --- /dev/null +++ b/Dragon/src/NXTITA.f @@ -0,0 +1,308 @@ +*DECK NXTITA + FUNCTION NXTITA(POSTRI ,PINPOS,VOLINT) +* +*---------- +* +*Purpose: +* Compute the volume of intersection between +* a 2--D triangle and an annular pin. +* +*Copyright: +* Copyright (C) 2010 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version. +* +*Author(s): G. Marleau. +* +*Parameters: input +* POSTRI spatial description of the triangle with: +* POSTRI(1,J) $X$ location of corner $j$ of triangle; +* POSTRI(2,J) $Y$ location of corner $j$ of triangle. +* PINPOS spatial description of the annular pin region with: +* PINPOS(0) the radius of the annular pin; +* PINPOS(1) the $X$ position of the annular pin center; +* PINPOS(2) the $Y$ position of the annular pin center. +* +*Parameters: output +* NXTITA type of intersection between haxagon and annular pin, where: +* = 0 means that there is no intersection +* between the two regions; +* = 1 means that the hexagon +* is all located inside the annular pin; +* = 2 means that the annular pin +* is all located inside the hexagon; +* =-1 means that the intersection between +* the hexagon and the annular pin is partial. +* VOLINT 2-D volume of intersection (area) between hexagon and +* annular pin. +* +*Reference: +* G. Marleau, +* New Geometries Processing in DRAGON: The NXT: Module, +* Report IGE-260, Polytechnique Montreal, +* Montreal, 2005. +* +*---- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER NXTITA + DOUBLE PRECISION POSTRI(2,3),PINPOS(0:2) + DOUBLE PRECISION VOLINT +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='NXTITA') + INTEGER IPRINT + PARAMETER (IPRINT=100) + DOUBLE PRECISION DCUTOF + PARAMETER (DCUTOF=1.0D-8) + DOUBLE PRECISION DZERO,DONE,DHALF,DSQ3O2 + PARAMETER (DZERO=0.0D0,DONE=1.0D0, + > DHALF=0.5D0,DSQ3O2=0.86602540378444D0) +*---- +* Functions +*---- + DOUBLE PRECISION XDRCST,PI +*---- +* Local variables +*---- + INTEGER IDIR,IFACE,NFPINS,NCIN,NKPINS,ICOR,ICORIN(3),IC + DOUBLE PRECISION VOLOUT,RADP2,X1,Y1,X2,Y2,VOLTRI,ALPHA + DOUBLE PRECISION CPPP(2),CORPOS(2),DIRFAC(2,3),DIRLIN(2,3),RADC, + > DIRN,DISTM,DISTC,DNOR(3),DTAN(3), + > VT,AFPINT(3),XFPINT(3), + > YFPINT(3),TC1(2),TC2(2),DISTS,DISTF +*---- +* Print header if required +*---- + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6000) NAMSBR + WRITE(IOUT,6010) (POSTRI(1,ICOR),POSTRI(2,ICOR),ICOR=1,3) + WRITE(IOUT,6011) (PINPOS(IFACE),IFACE=0,2) + ENDIF +*---- +* Initialize PI, NXTITA and VOLINT +*---- + PI=XDRCST('Pi',' ') + NXTITA=0 + VOLINT=DZERO + VOLOUT=DZERO +*---- +* Evaluate distance from FACES to pin center +*---- + RADP2=PINPOS(0)**2 + NFPINS=0 + NCIN=0 + NKPINS=0 +*---- +* Compute volume of triangle +*---- + X1=POSTRI(1,2)-POSTRI(1,1) + X2=POSTRI(1,3)-POSTRI(1,1) + Y1=POSTRI(2,2)-POSTRI(2,1) + Y2=POSTRI(2,3)-POSTRI(2,1) + VOLTRI=ABS(X1*Y2-X2*Y1)/2.0D0 +*---- +* Analyze each face +*---- + CPPP(1)=POSTRI(1,2) + CPPP(2)=POSTRI(2,2) + CORPOS(1)=POSTRI(1,3) + CORPOS(2)=POSTRI(2,3) + DO IFACE=1,3 + ICOR=IFACE + ICORIN(ICOR)=0 +* write(6,'(A6,5x,2i5)') 'FACE ',IFACE,ICOR + AFPINT(IFACE)=DZERO + XFPINT(IFACE)=DZERO + YFPINT(IFACE)=DZERO + DNOR(IFACE)=DZERO + DTAN(IFACE)=DZERO + RADC=DZERO + DISTM=DZERO +*---- +* Find direction of face and its normal directed inward +* the triangle +*---- + DIRLIN(1,IFACE)=(POSTRI(1,IFACE)-CORPOS(1)) + DIRLIN(2,IFACE)=(POSTRI(2,IFACE)-CORPOS(2)) + DIRN=SQRT(DIRLIN(1,IFACE)**2+DIRLIN(2,IFACE)**2) + DIRLIN(1,IFACE)=DIRLIN(1,IFACE)/DIRN + DIRLIN(2,IFACE)=DIRLIN(2,IFACE)/DIRN + DIRFAC(1,IFACE)=-DIRLIN(2,IFACE) + DIRFAC(2,IFACE)=DIRLIN(1,IFACE) +* write(6,*) 'Corner and pin position with respect to corner' +* write(6,'(4F20.10)') +* > CORPOS(1),CORPOS(2),PINPOS(1)-CORPOS(1), +* > PINPOS(2)-CORPOS(2) +* write(6,*) 'Face tangent and normal' +* write(6,'(4F20.10)') +* > DIRLIN(1,IFACE),DIRLIN(2,IFACE),DIRFAC(1,IFACE),DIRFAC(2,IFACE) + DO IDIR=1,2 + DISTC=PINPOS(IDIR)-CORPOS(IDIR) + DISTM=DISTM+(CPPP(IDIR)-CORPOS(IDIR))*DIRFAC(IDIR,IFACE) + DNOR(IFACE)=DNOR(IFACE)+DISTC*DIRFAC(IDIR,IFACE) + DTAN(IFACE)=DTAN(IFACE)+DISTC*DIRLIN(IDIR,IFACE) + RADC=RADC+DISTC**2 + ENDDO +* write(6,*) 'Distance of center to face',DNOR(IFACE) + IF(DNOR(IFACE) .GT. DISTM+PINPOS(0)) THEN + NXTITA=0 + VOLINT=DZERO + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6012) NAMSBR,NXTITA,VOLINT + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN + ELSE IF(DNOR(IFACE) .LT. -PINPOS(0)) THEN + NXTITA=0 + VOLINT=DZERO + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6012) NAMSBR,NXTITA,VOLINT + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN + ELSE IF(DNOR(IFACE) .GE. PINPOS(0)) THEN + NFPINS=NFPINS+1 + NKPINS=NKPINS+1 + ELSE +*---- +* Point of intersection of current face with pin +*---- + XFPINT(IFACE)=DNOR(IFACE) + YFPINT(IFACE)=SQRT(RADP2-XFPINT(IFACE)**2) + AFPINT(IFACE)=ACOS(XFPINT(IFACE)/PINPOS(0)) + VT=XFPINT(IFACE)*YFPINT(IFACE) + DISTS=DTAN(IFACE)-YFPINT(IFACE) + DISTF=DTAN(IFACE)+YFPINT(IFACE) +* write(6,*) 'DISTS/DISTF/DIRN=',DISTS,DISTF,DIRN + IF(DISTS .LE. DIRN .AND. DISTF .GT. DZERO) THEN +* write(6,'(A9,2X,7F20.10)') +* > 'Volout 1=',XFPINT(IFACE),YFPINT(IFACE), +* > AFPINT(IFACE),VT,VOLOUT, +* > RADP2*AFPINT(IFACE)-VT,VOLOUT+RADP2*AFPINT(IFACE)-VT + VOLOUT=VOLOUT+RADP2*AFPINT(IFACE)-VT + ELSE + NKPINS=NKPINS+1 +* write(6,'(A9,2X,7F20.10)') +* > 'Volout 3=',XFPINT(IFACE),YFPINT(IFACE), +* > AFPINT(IFACE),VT,VOLOUT, +* > RADP2*AFPINT(IFACE)-VT,VOLOUT + VOLOUT=VOLOUT + ENDIF + ENDIF + IF(RADC .LT. RADP2) THEN +*---- +* Identify corners in pin +*---- + ICORIN(ICOR)=1 + NCIN=NCIN+1 + ENDIF + CPPP(1)=CORPOS(1) + CPPP(2)=CORPOS(2) + CORPOS(1)=POSTRI(1,IFACE) + CORPOS(2)=POSTRI(2,IFACE) + ENDDO + IF(NFPINS .EQ. 3) THEN +*---- +* Pin completely inside triangle +*---- + NXTITA=2 + VOLINT=PI*RADP2 + ELSE IF(NCIN .EQ. 3) THEN +*---- +* Triangle completely inside pin +*---- + NXTITA=1 + VOLINT=VOLTRI + ELSE IF(NKPINS .EQ. 3) THEN +*---- +* No intersection between triangle and pin +*---- + NXTITA=0 + VOLINT=DZERO + ELSE +*---- +* For corners inside pin, find intersection of outside surfaces +* and remove from VOLOUT +*---- + NXTITA=-1 + DO ICOR=1,3 +* write(6,*) 'ICORIN=',ICOR,ICORIN(ICOR) + IF(ICORIN(ICOR) .EQ. 1) THEN +*---- +* Point of intersection of previous face with pin in the positive +* direction +*---- + IFACE=ICOR-1 + IF(IFACE .LE. 0) IFACE=3+IFACE + IC=ICOR-2 + IF(IC .LE. 0) IC=3+IC + DO IDIR=1,2 + TC2(IDIR)=POSTRI(IDIR,IC) + > +(DTAN(IFACE)+YFPINT(IFACE))*DIRLIN(IDIR,IFACE) + ENDDO +* write(6,*) 'TC2=',TC2(1),TC2(2) +*---- +* Point of intersection of current face with pin in the negative +* direction +*---- + IFACE=ICOR + IC=ICOR-1 + IF(IC .LE. 0) IC=3+IC + DO IDIR=1,2 + TC1(IDIR)=POSTRI(IDIR,IC) + > +(DTAN(IFACE)-YFPINT(IFACE))*DIRLIN(IDIR,IFACE) + ENDDO +* write(6,*) 'TC1=',TC1(1),TC1(2) +*---- +* Triangle outside is identified by CORPOS(*,IC),TC1,TC2 +* Compute its volume +*---- + X1=TC1(1)-POSTRI(1,IC) + X2=TC2(1)-POSTRI(1,IC) + Y1=TC1(2)-POSTRI(2,IC) + Y2=TC2(2)-POSTRI(2,IC) + VOLTRI=ABS(X1*Y2-X2*Y1)/2.0D0 +* write(6,*) 'Triangle 1=',X1,X2,Y1,Y2,VOLTRI +*---- +* Add contribution from annular sector between T1 and T2 +* Compute distance of T1 to T2 +*---- + Y2=DZERO + DO IDIR=1,2 + Y2=Y2+(TC1(IDIR)-TC2(IDIR))**2 + ENDDO + Y2=Y2/4.0D0 + X2=SQRT(RADP2-Y2) + Y2=SQRT(Y2) + ALPHA=ACOS(X2/PINPOS(0)) +* write(6,*) 'Triangle 2=',X2,Y2,ALPHA, +* > VOLOUT,RADP2*ALPHA,X2*Y2,RADP2*ALPHA-X2*Y2+VOLTRI + VOLOUT=VOLOUT-RADP2*ALPHA+X2*Y2-VOLTRI + ENDIF + ENDDO + VOLINT=PI*RADP2-VOLOUT +* write(6,*) 'Triangle 3=',VOLINT,PI*RADP2,VOLOUT + ENDIF + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6012) NAMSBR,NXTITA,VOLINT + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6010 FORMAT('POSTRI={',5(F20.10,','),F20.10,'};') + 6011 FORMAT('PINPOS={',2(F20.10,','),F20.10,'};') + 6012 FORMAT(A6,'={',I5,',',F20.10,'};') + END |