SUBROUTINE TRIANINT(NMAP,M,ISWT,ITIME) USE BLKMAP USE BLK1MOD SAVE ! INCLUDE 'BLK1.COM' DIMENSION WGT(8) REAL*8 XMINL,YMINL,XMAXL,YMAXL ! data itime/0/ ! LOOK FOR MATCHING POINTS ITIMESKP=0 DO K=1,MAXPTS DISQ=(XUSR(M)-XMAP(K))**2+(YUSR(M)-YMAP(K))**2 IF(DISQ .LT. 1.) THEN WD(M)=VAL(K) FPN = WD(M)*10. X = CORD(M,1) Y = CORD(M,2) - .11 IF(X .GT. 0. .AND. X .LT. HSIZE .AND. & & Y .GT. 0. .AND. Y .LT. 7.5) THEN CALL RRED CALL NUMBR(X,Y,0.1,FPN,0.0,-1) endif ITIMESKP=1 GO TO 300 ENDIF ENDDO ! Search for element that has circumcircle around the node IF(ISWT .NE. 0) THEN IF(ITIME .EQ. 0) NSTART=1 ELSE NSTART=1 ENDIF DO N=NSTART,NELTS IF(NOPEL(N,1) .EQ. 0) GO TO 200 if(RADS(N) .eq. 0.) then CALL CCENTRE(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3)) & &,YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3)) & &,XCEN(N),YCEN(N),RADS(N)) endif IF(RADS(N)+XCEN(N) .GE. XUSR(M)) THEN NSTART=N GO TO 210 ENDIF 200 CONTINUE ENDDO 210 CONTINUE 220 continue ! WRITE(155,*) M,NSTART DO N=NSTART,NELTS IF(NOPEL(N,1) .EQ. 0) GO TO 250 if(RADS(N) .eq. 0.) then CALL CCENTRE(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3)) & &,YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3)) & &,XCEN(N),YCEN(N),RADS(N)) endif xminl=min(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3))) xmaxl=max(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3))) yminl=min(YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3))) ymaxl=max(YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3))) ! IF(M .EQ. 6316) THEN ! WRITE(156,'(2I6,6F15.2)') M,N,XUSR(M),XMINL,XMAXL,YUSR(M),YMINL,YMAXL ! ENDIF if(xusr(m) .lt. xminl-0.01 .or. xusr(m) .gt. xmaxl+0.01) then go to 250 elseif(yusr(m) .lt. yminl-0.01 .or. yusr(m) .gt. ymaxl+0.01) then go to 250 endif ! IF(M .EQ. 6316) WRITE(156,*) 'PASSED X AND Y TEST',N DISQ=(XUSR(M)-XCEN(N))**2+(YUSR(M)-YCEN(N))**2 IF(DISQ .LE. RADS(N)**2*1.0001) THEN ! IF(M .EQ. 6316) write(156,*) m,n,disq,rads(n)**2,xusr(m),xcen(n) ! We have a candidate CALL GETWT(N,XUSR(M),YUSR(M),WGT,1) DO K=1,3 IF(WGT(K) .LT. -1E-4 .OR. WGT(K) .GT. 1.0001) THEN WRITE(142,*) 'REJECT',m,n,disq,rads(n)**2,wgt(1),wgt(2),wgt(3) if(nstart .gt. 1) then nstart=1 go to 220 endif GO TO 250 ENDIF ENDDO WD(M)=WGT(1)*VAL(NOPEL(N,1))+WGT(2)*VAL(NOPEL(N,2))+WGT(3)*VAL(NOPEL(N,3)) FPN = WD(M)*10. X = CORD(M,1) Y = CORD(M,2) - .11 IF(X .GT. 0. .AND. X .LT. HSIZE .AND. & & Y .GT. 0. .AND. Y .LT. 7.5) THEN CALL RRED CALL NUMBR(X,Y,0.1,FPN,0.0,-1) endif GO TO 300 ENDIF 250 CONTINUE ENDDO 300 CONTINUE IF(ITIMESKP .EQ. 0) ITIME=1 RETURN END SUBROUTINE GETWT(N,XSW,YSW,WGT,ISWT) !- !......SUBROUTINE TO EVALUATE FUNCTION AT GRID POINTS !- !- N = ELEMENT NUMBER !_ XSW = X COORDINATE OF DESIRED POINT !_ YSW = Y COORDINATE OF DESIRED POINT ! WGT(8) = ARRAY OF WEIGHTING FUNCTIONS ! ISWT = SWITCH FOR CHOICE BETWEEN LINEAR AND QUADRATIC WEIGHTING ! = 1 FOR LINEAR ! = 2 FOR QUADRATIC ! FROM COMMON ! NOP = LIST OF NODAL CONNECTIONS AROUND AN ELEMET ! CORD = REAL*8 ARRAY OF NODAL COORDINATES ! USE BLKMAP USE BLK1MOD REAL*8 XN,DNX,DNY,XSW,YSW DOUBLE PRECISION XG,YG,XK,YK,XP,YP ! INCLUDE 'BLK1.COM' !- DIMENSION X(9),Y(9),WGT(8) !- DATA TOL/0.01/ !- !- !......DETERMINE ELEMENT TYPE !- !IPKOCT93 ADD if(n .eq. 1910) then aaa=0 endif NCN=6 IT=2 !- !......ESTABLISH LOCAL COORDINATES FOR EACH NODE POINT OF ELEMENT !- K1=NOPEL(N,1) X(1)=0. Y(1)=0. DO 300 K=3,NCN,2 K2=NOPEL(N,K/2+1) X(K)=XMAP(K2)-XMAP(K1) Y(K)=YMAP(K2)-YMAP(K1) 300 END DO X(2)=X(3)/2. Y(2)=Y(3)/2. X(4)=(X(3)+X(5))/2. Y(4)=(Y(3)+Y(5))/2. X(6)=X(5)/2. Y(6)=Y(5)/2. xminl=min(x(1),x(3),x(5)) yminl=min(y(1),y(3),y(5)) xmaxl=max(x(1),x(3),x(5)) ymaxl=max(y(1),y(3),y(5)) !- !......ESTABLISH LOCAL COORDINATES OF DESIRED POINT !- XP=XSW-XMAP(K1) YP=YSW-YMAP(K1) if(xp .lt. xminl .or. xp .gt. xmaxl) then wgt(1)=2.0 return elseif(yp .lt. yminl .or. yp .gt. ymaxl) then wgt(1)=2.0 return endif XG=0. YG=0. !- !......ITERATE TO FIND LOCAL COORDINATE !- DO 400 ITER=1,10 DXKDX=0. DXKDY=0. DYKDX=0. DYKDY=0. XK=-XP YK=-YP DO 350 K=2,NCN XK=XK+XN(IT,K,XG,YG)*X(K) YK=YK+XN(IT,K,XG,YG)*Y(K) DXKDX=DXKDX+DNX(IT,K,XG,YG)*X(K) DYKDX=DYKDX+DNX(IT,K,XG,YG)*Y(K) DXKDY=DXKDY+DNY(IT,K,XG,YG)*X(K) DYKDY=DYKDY+DNY(IT,K,XG,YG)*Y(K) 350 END DO DET=DXKDX*DYKDY-DXKDY*DYKDX DX=(-DYKDY*XK+DXKDY*YK)/DET DY=( DYKDX*XK-DXKDX*YK)/DET XG=XG+DX YG=YG+DY IF(ABS(DX).LT.TOL .AND. ABS(DY).LT.TOL) GO TO 420 400 END DO !- !......NOW GET WEIGHTING FUNCTIONS FOR QUAD FUNCTION !- 420 CONTINUE DO K=1,NCN WGT(K)=XN(IT,K,XG,YG) END DO IF(ISWT .EQ. 1) THEN !- !- REDUCE TO LINEAR FUNCTION BY ADDING TERMS !- DO K=2,NCN,2 WGT(K-1)=WGT(K-1)+WGT(K)/2. IF(K .LT. NCN) THEN WGT(K+1)=WGT(K+1)+WGT(K)/2. ELSE WGT(1)=WGT(1)+WGT(K)/2. ENDIF ENDDO !- !- THEN COMPACT ARRAY !- DO K=1,NCN/2 WGT(K)=WGT(2*K-1) ENDDO ENDIF RETURN END