C JWISIA.FOR H18.5.5 BY S.I（ソース公開ver.1） C C PARAMETER文での配列宣言用定数の説明 C NELM;要素（囲み）数、NODE;節点（交点）数 C IHENGD;入力線数、IHENTD;辺数、ITENSD;１直線上の交点数 C NDP;節点からの辺数、NLN;１要素の節点数 C IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) PARAMETER ( NELM = 3000, NODE = 3000 ) PARAMETER ( IHENGD = 5000, IHENTD = 8000 ) PARAMETER ( ITENSD = 1000, NDP = 50, NLN = 500,NLN1=NLN+1 ) CHARACTER LAYHA2*1 DIMENSION SEN1(2,IHENGD), SEN2(2,IHENGD), ISENC(IHENGD) DIMENSION TEN(2,NODE), IHEN(2,IHENTD), IHENC(IHENTD) DIMENSION ITENS(ITENSD), XX(ITENSD), NN(ITENSD) DIMENSION YX(2,NODE), JUN(NODE) DIMENSION NODK(NODE), NODP(NDP,NODE), NOD(0:NLN,NELM) DIMENSION XYG(2,NELM), YXG(2,NELM), N(0:NLN), M(NLN) DIMENSION XEL(NLN1), YEL(NLN1), JUN2(NELM), IHEC0(IHENTD) C OPEN ( 2, FILE = 'JWC_TEMP.TXT' ) OPEN ( 7, FILE = 'WORK1.TXT' ) OPEN ( 8, FILE = 'AREPRT.TXT' ) CALL SCALIN ( SCAL ) C 図寸1mmの1/10倍（図寸0.1mm）を判定誤差としている。 GOSA = SCAL*0.1 IKETA = KETKET( GOSA ) C WRITE(*,*) '番号の付け方；1か2を入力して下さい。' WRITE(*,*) ' =1;左下から右上、=2；左上から右下 ' WRITE(*,*) ' ' READ(*,*) INB WRITE(*,*) ' ' C WRITE(*,*) '判定誤差（o）＝',GOSA WRITE(*,*) ' ' WRITE(*,*) '現在計算中' C C 直線の読み込み、及び辺、節点ﾃﾞｰﾀの作成 CALL SENDAT( SEN1, SEN2, ISENC, TEN, IHEN, IHENC + ,ITENS, XX, NN, IHENGD, IHENTD, ITENSD, IT, IH, IHENK + ,NODK, NODP, NDP, NODE, XXRR, YYUU, GOSA ) C statim = timef() C C 囲みの探索（節点と辺のデータから要素データを作成する） CALL KAKOMI( TEN, IHEN, IHENC, IHENTD, IT, IH, IHENK + , NODK, NODP, NDP, NODE, NLN, N, M, XEL, YEL, NLN1 + , NELM, IE, ITES, XYG, NOD, IHEC0, GOSA ) C write(*,*) ' 主な計算時間(秒)=', timef() write(*,*) ' 要素数=',ie, ' 主な計算回数=',ites C WRITE(8,*) ' 　主な計算回　＝　',ITES WRITE(8,*) ' 節点（交点）数＝　',IT WRITE(8,*) ' 要素（囲み）数＝　',IE IT2= IT*2 CALL DECOPY ( TEN, YX, IT2 ) C C 節点YXの順を右へ及び上へ変更しTENに、順をJUNに出力 C JUN(旧節点番号)＝新節点番号 CALL NODCHG ( YX, TEN, JUN, IT, SCAL, 1, INB ) C WRITE(8,*) ' 節点ＮＯ．　X(mm）、　y(mm)' WRITE(8,'(I4,2E15.6)') ( I,( TEN(J,I),J=1,2),I=1,IT) C C 要素重心XYGの順を右へ及び上へ変更しYXGに、順をJUN2に出力 C JUN2(新要素番号)＝旧要素番号 CALL NODCHG ( XYG, YXG, JUN2, IE, SCAL, 2, INB ) C C 以下は８番ファイルへの結果出力と7番ファイルへの作図出力 C 　HHH3、HHH4；節点番号、要素番号文字高さ(mm) C 　XXRR、YYUU；図形の右端、上端座標(mm) C WRITE(8,*) ' 要素NO.　面積(mm2）、重心（Ｘ、Ｙ）、点数、　節点NO.' XXRR = XXRR + SCAL*20 XS = 0. YS = 0. HHH3 = 3.0 HHH4 = 4.0 ANG = 0.0 WRITE(7,'(3Hlyf)') C DO 70 I = 1, IT DATDAT = FLOAT(I) XXGG = TEN(1,I)-HHH3*SCAL*0.5*( 3.5-0.5*AINT(LOG10(DATDAT)) ) YYGG = TEN(2,I)-HHH3*SCAL*0.5 CALL NUMBRS ( XXGG, YYGG, HHH3, DATDAT, 4, -1, ANG, 0 ) 70 CONTINUE C WRITE(7,'(3Hcn5)') WRITE(7,1001) XXRR+0.0*SCAL, YYUU+10.*SCAL, 10., 0.0 WRITE(7,1002) XXRR+15.*SCAL, YYUU+10.*SCAL, 10., 0.0 1001 FORMAT(2Hch,4E14.6,2H ",3HNO.) 1002 FORMAT(2Hch,4E14.6,2H ",10H面積（u）) WRITE(7,'(3Hlc5)') WRITE(7,'(3Hlt9)') ARET = 0.0 DO 80 II = 1, IE I = JUN2(II) ARE = AREA ( NOD, N, YX, NLN, NODE, I, NELM ) WRITE(8,'(I4,3E15.6,I4,20I5,/(23X,20I5))') + II, ARE, YXG(1,II), YXG(2,II), + NOD(0,I), ( JUN(N(J)), J = 1, NOD(0,I) ) WRITE(7,'(3Hly0)') DATDAT = FLOAT(II ) XXGG = YXG(1,II)-HHH4*SCAL*0.5*( 3.5-0.5*AINT(LOG10(DATDAT)) ) YYGG = YXG(2,II)-HHH4*SCAL*0.5 XS = XXRR + FLOAT(II/100)*40.*SCAL YS = YYUU - FLOAT(MOD(II,100))*5.*SCAL CALL NUMBRS ( XS, YS,HHH3, DATDAT, 4, -1, ANG, 0 ) ARE = ARE / 1.D6 ARET = ARET + ARE CALL NUMBRS + ( XS+10.*SCAL, YS, HHH3, ARE, 10, IKETA, ANG, 0 ) WRITE(7,'(2Hly,A1)') LAYHA2(MOD(II-1,14)+1) WRITE(7,'(3Hlc3)') WRITE(7,'(3Hlt9)') CALL NUMBRS ( XXGG, YYGG, HHH4, DATDAT, 4, -1, ANG, 0 ) WRITE(7,'(4E14.6)') + YXG(1,II)-SCAL*3.0, YXG(2,II)-SCAL*3.0 + , YXG(1,II)+SCAL*3.0, YXG(2,II)+SCAL*3.0 WRITE(7,'(4E14.6)') + YXG(1,II)-SCAL*3.0, YXG(2,II)+SCAL*3.0 + , YXG(1,II)+SCAL*3.0, YXG(2,II)-SCAL*3.0 WRITE(7,'(3Hlc5)') WRITE(7,'(3Hlt9)') DO 81 JJ = 1, NOD(0,I) WRITE(7,'(4E14.6)') YX(1,N(JJ)), YX(2,N(JJ)) + ,YX(1,N(JJ+1)), YX(2,N(JJ+1)) 81 CONTINUE 80 CONTINUE WRITE(7,'(3Hly0)') WRITE(7,'(3Hcn5)') WRITE(7,1003) XXRR+0.*SCAL, YYUU+25.*SCAL, 10., 0.0 CALL NUMBRS +(XXRR+7.*5.*SCAL,YYUU+26.*SCAL,HHH4,ARET,13,IKETA,ANG, 0 ) 1003 FORMAT(2Hch,4E14.6,2H ",12H総面積（u）) C C 7番(WORK1.TXT)へ出力した作図データをJWC_TEMP.TXTへ転送 CALL JWCOUT C STOP END C C 点１と点２の間の線分長を与える FUNCTION SENL ( X1, Y1, X2, Y2 ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) X12 = ABS(X1-X2) Y12 = ABS(Y1-Y2) SENL = X12**2 + Y12**2 IF ( SENL. EQ. 0. ) RETURN SENL = SQRT( SENL ) RETURN END C C 囲みで点Ｋを既に通っていないかﾁｪｯｸする C =1なら既点有り、=0なら既点なし FUNCTION NKHAN( N, K ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION N(0:K) NKHAN = 1 DO 10 I = 1, K-2 IF ( N(I). EQ. N(K) ) RETURN 10 CONTINUE NKHAN = 0 RETURN END C C 点１と２がある辺の使用可能回数IHENCを出力する FUNCTION IHEMC( N1, N2, IHEN, IHENC, IH ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION IHEN(2,IH), IHENC(IH) I1 = MAX(N1,N2) I2 = MIN(N1,N2) DO 10 I = 1, IH J1 = IHEN(1,I) J2 = IHEN(2,I) IF ( J1. EQ. I1. AND. J2. EQ. I2 ) THEN IHEMC = IHENC(I) RETURN ENDIF 10 CONTINUE WRITE(*,*) ' ERROR IN IHEMC ' STOP RETURN END C C 点１と２がある辺の使用可能回数IHENCを１減らす SUBROUTINE IHELC( N1, N2, IHEN, IHENC, IHENK, IH ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION IHEN(2,IH), IHENC(IH) I1 = MAX(N1,N2) I2 = MIN(N1,N2) DO 10 I = 1, IH J1 = IHEN(1,I) J2 = IHEN(2,I) IF ( J1. EQ. I1. AND. J2. EQ. I2 ) THEN IHENC(I) = IHENC(I)-1 IHENK = IHENK - 1 RETURN ENDIF 10 CONTINUE WRITE(*,*) ' ERROR IN IHELC ' RETURN END C C 点配列XYを左から右へDET毎上へ並び替えYXとして出力する。 C DETは同じ標高の段とするかしないかの標高差(mm)用判定子 C 新旧順番の対応はJUNに以下のように出力する。 C IHA=1ならJUN(旧NO)=新NO　←節点用 C IHA=2ならJUN(新NO)=旧NO　←要素用 SUBROUTINE NODCHG( XY, YX, JUN, NNP, DET, IHA, INB ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION XY(2,NNP), JUN(NNP), X(3000), Y(3000), IYS(3000) DIMENSION YX(2,NNP), N(3000), NN(3000), XX(3000), YY(3000) DIMENSION OS(2000) C 番号の付け方；INB=1;下から上、INB=2;上から下 C DATA INB / 1 / DO 1 I = 1, NNP X(I) = XY(1,I) Y(I) = XY(2,I) 1 CONTINUE C IF ( INB. EQ. 1 ) THEN CALL MINJUN ( Y, N, NNP ) ELSE CALL MAXJUN ( Y, N, NNP ) ENDIF C IY = 1 IYS(1) = 1 IF ( INB. EQ. 1 ) THEN Y0 = -1.E20 ELSE Y0 = 1.E20 ENDIF DO 10 I = 1, NNP Y1 = Y(N(I)) IF ( INB. EQ. 1. AND. Y1. GT. Y0+DET +. OR.INB. NE. 1. AND. Y1. LT. Y0-DET ) THEN IY = IY + 1 IYS(IY) = I ENDIF Y0 = Y1 10 CONTINUE IYK = IY IYS(IYK+1) = NNP+1 C DO 20 I = 1, IYK II = 0 DO 21 J = IYS(I), IYS(I+1)-1 II = II + 1 XX(II) = X(N(J)) YY(II) = Y(N(J)) OS(II) = N(J) 21 CONTINUE CALL MINJUN ( XX, NN, IYS(I+1)-IYS(I) ) DO 22 J = IYS(I), IYS(I+1)-1 III = NN(J-IYS(I)+1) IF ( IHA. EQ. 1 ) JUN(OS(III)) = J IF ( IHA. EQ. 2 ) JUN(J) = OS(III) YX(1,J) = XY(1,OS(III)) YX(2,J) = XY(2,OS(III)) 22 CONTINUE 20 CONTINUE C RETURN END C C Ｘ値の昇順をＮに返す SUBROUTINE MINJUN ( X, N, K ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) C DIMENSION X(K), N(K) C DO 10 I = 1, K N(I) = I 10 CONTINUE C IF ( K. LE. 1 ) RETURN C DO 20 I = 1, K-1 SM = 1.E10 L = I DO 21 J = I, K IF ( X( N(J) ). LT. SM ) THEN SM = X( N(J) ) L = J ENDIF 21 CONTINUE L0 = N(I) N(I) = N(L) N(L) = L0 20 CONTINUE C RETURN END C C Ｘ値の降順をＮに返す SUBROUTINE MAXJUN ( X, N, K ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) C DIMENSION X(K), N(K) C DO 10 I = 1, K N(I) = I 10 CONTINUE C IF ( K. LE. 1 ) RETURN C DO 20 I = 1, K-1 SM = -1.E10 L = I DO 21 J = I, K IF ( X( N(J) ). GT. SM ) THEN SM = X( N(J) ) L = J ENDIF 21 CONTINUE L0 = N(I) N(I) = N(L) N(L) = L0 20 CONTINUE C RETURN END C C 16進数のﾚｰﾔｰNOを10進数のNOに変える FUNCTION LAYHAN(LY) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) CHARACTER LY*1, CLY(16)*1 DATA CLY / '0','1','2','3','4','5','6','7' + ,'8','9','a','b','c','d','e','f' / DO 10 I = 1, 16 IF ( LY. EQ. CLY(I) ) THEN LAYHAN = I RETURN ENDIF 10 CONTINUE WRITE(*,*) ' ERROR IN LAYHAN ' RETURN END C C 10進数を16進数のﾚｰﾔｰNOに変える FUNCTION LAYHA2(LY) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) CHARACTER LAYHA2*1, CLY(0:15)*1 DATA CLY / '0','1','2','3','4','5','6','7' + ,'8','9','a','b','c','d','e','f' / LL = MOD(LY,16) LAYHA2 = CLY(LL) RETURN END C C TXTDATの文字数を返す FUNCTION NUMCHA ( TXTDAT ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) CHARACTER TXTDAT*500 NUMCHA = 0 DO 10 I = 500, 1, -1 IF ( TXTDAT(I:I) . NE. ' ' ) THEN NUMCHA = I RETURN ENDIF 10 CONTINUE RETURN END C C 閉合点列の面積算定と点列並びを反時計回りに変更する。 FUNCTION AREA ( NOD, N, TEN, NLN, NODE, III,NELM ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION NOD(0:NLN,NELM), N(0:NLN), TEN(2,NODE) DIMENSION XP(3000), YP(3000) C 閉合点列の面積を求める。 K = NOD(0,III)+1 DO 10 I = 1, K-1 N(I) = NOD(I,III) XP(I) = TEN(1,N(I)) YP(I) = TEN(2,N(I)) 10 CONTINUE N(K) = N(1) XP(K) = XP(1) YP(K) = YP(1) C AREA = 0. DO 30 I = 1, K-1 X1 = XP(I) Y1 = YP(I) X2 = XP(I+1) Y2 = YP(I+1) AREA = AREA + ( X2-X1 ) * ( Y2 + Y1 ) * 0.5 30 CONTINUE C C 面積が＋なら点列は時計回りなので反時計回りに変更する IF ( AREA. GT. 0. ) THEN DO 40 I = 2, K-1 N(I) = NOD(K-1-I+2,III) 40 CONTINUE ENDIF C AREA = ABS ( AREA ) RETURN END C C WORK1.TXT(ﾒｲﾝでは7番)へ出力した作図データをJWC_TEMP.TXTへ転送 SUBROUTINE JWCOUT IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) CHARACTER TXTDAT*500, FMT*6 CLOSE(2) CLOSE(7) OPEN ( 2, FILE = 'JWC_TEMP.TXT' ) OPEN ( 7, FILE = 'WORK2.TXT') OPEN ( 9, FILE = 'WORK1.TXT' ) C JWC_TEMP.TXTの諸条件ﾃﾞｰﾀ2行〜13行をWORK2.TXTへ転送 READ(2,'(A500)') TXTDAT DO 1 I = 1, 12 READ(2,'(A500)',END=2) TXTDAT WRITE(FMT,'(2H(A,I3,1H))') NUMCHA(TXTDAT) WRITE(7,FMT) TXTDAT(1:NUMCHA(TXTDAT)) 1 CONTINUE 2 CONTINUE C 続けてWORK1.TXTの作図ﾃﾞｰﾀをWORK2.TXTへ書き込む 3 READ(9,'(A500)',END=4) TXTDAT WRITE(FMT,'(2H(A,I3,1H))') NUMCHA(TXTDAT) WRITE(7,FMT) TXTDAT(1:NUMCHA(TXTDAT)) GO TO 3 4 CONTINUE C WORK2.TXTを巻戻し、JWC_TEMP.TXTを開きなおす REWIND 7 CLOSE (2) OPEN ( 2, FILE = 'JWC_TEMP.TXT' ) C WORK2.TXTをJWC_TEMP.TXTをへ転送 5 READ( 7,'(A500)',END=6) TXTDAT WRITE(FMT,'(2H(A,I3,1H))') NUMCHA(TXTDAT) WRITE(2,FMT) TXTDAT(1:NUMCHA(TXTDAT)) GO TO 5 6 CONTINUE RETURN END C C 実数DATAを半角右詰めで描画する。 C XS,YS：描画始点（左下）、H：文字高、DATA：数値（実数） C I1：描画全文字数、I2：小数点以下描画ケタ数（整数で描くなら-1） C 例えば、123.123を123.12と描くならI1>=8、I2=2 C AN：文字のX軸からの角度（水平なら0.0，垂直上向きなら90.0） C I：=0なら+なし、=1なら+描く SUBROUTINE NUMBRS ( XS, YS, H, DAT, I1, I2, AN, I ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) CHARACTER FMT*30, HANNUM*40 DATA XMOD, YMOD, XCMO, YCMO / 0., 0., 1., 1./ FMT(1:30) = ' ' HANNUM(1:40) = ' ' ICN = NINT ( H*XCMO ) IF ( H*XCMO. LT. 2.25 ) ICN = 1 IF ( H*XCMO. GE. 2.25. AND. H*XCMO. LT. 2.75 ) ICN = 2 IF ( ICN. GT. 10 ) ICN = 10 WRITE(7,'(2Hcn,I1)') ICN IF ( I2. GE. 0 ) THEN WRITE(FMT,'(2H(F,I2,1H.,I2,1H))') I1,I2 WRITE(HANNUM,FMT) DAT IF ( HANNUM(I1-I2-1:I1-I2). EQ. ' .' ) + HANNUM(I1-I2-1:I1-I2) = '0.' IF ( HANNUM(I1-I2-1:I1-I2). EQ. '-.' ) + HANNUM(I1-I2-2:I1-I2) = '-0.' ELSE WRITE(FMT,'(2H(I,I2,1H))') I1 WRITE(HANNUM,FMT) NINT(DAT) ENDIF IF ( I. EQ. 1. AND. DAT. GT. 0. ) THEN KETA = LOG10(DAT) IF ( KETA. LT. 0 ) KETA = 0 HANNUM(I1-I2-1-KETA-1:I1-I2-1-KETA-1) = '+' ENDIF WRITE(FMT,'(20H(3Hch ,4E14.6,2H ",A,I2,1H))') I1 WRITE(7,FMT) + XS*XCMO+XMOD, YS*YCMO+YMOD, + 10.*COS(AN/57.296), 10.*SIN(AN/57.296), HANNUM RETURN END C C Ｄ(mm)正方をm2単位で表示するに適す小数点以下桁数を返す FUNCTION KETKET(D) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) A = ANINT ( LOG10( (D/1000.)**2 ) ) I = -1 IF ( A. LT. 1. ) I = 1 KETKET = I * NINT ( ABS(A) ) IF ( KETKET. LT. -1 ) KETKET = -1 IF ( KETKET. GT. 7 ) KETKET = 7 RETURN END C C JW_TEMP.TXTから書き込みﾚｰﾔｰのｽｹｰﾙを読む SUBROUTINE SCALIN ( SCALE ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER MOJ1*1, TXTDAT*500 DIMENSION SCAL(16) READ(2,'(A500)') ( TXTDAT, I = 1, 2 ) SCAL(1) = READAT(2,3) DO 10 I = 2, 16 SCAL(I) = READAT(0,-1) 10 CONTINUE DO 1002 I = 1, 5 READ(2,'(A500)') TXTDAT 1002 CONTINUE READ(2,'(2X,A1)') MOJ1 SCALE = SCAL(LAYHAN(MOJ1)) RETURN END C C MT番ﾌｧｲﾙから書式不明ﾃﾞｰﾀを読み込む FUNCTION READAT ( MT, M1 ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER TXTDAT*500, FMT*7 C MT；装置番号（前回と同じ行なら0にする） C M1；開始ｶﾗﾑ、=-1なら前回ｶﾗﾑ+1 READAT = 1.E20 IF ( M1. EQ. -1. AND. I2. EQ. 0 ) RETURN IF ( MT. NE. 0 ) READ(MT,'(A500)') TXTDAT IF ( M1. EQ. -1 ) M1 = I2+1 I1 = M1 I2 = 0 DO 10 I = I1, 500 IF ( TXTDAT(I:I). NE. ' ' ) I2 = I IF ( I2. EQ. 0 ) GO TO 10 IF ( TXTDAT(I:I). EQ. ' ' ) THEN I2 = I-1 GO TO 11 ENDIF 10 CONTINUE 11 CONTINUE IF ( I2. EQ. 0 ) RETURN WRITE(FMT,'(2H(F,I2,3H.0))') I2-I1+1 READ(TXTDAT(I1:I2),FMT) READAT RETURN END C C 辺の角度（ﾗｼﾞｱﾝ）を求める。 FUNCTION RADTAN ( X1, Y1, X2, Y2 ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PAI = 4.0*ATAN(1.0) IF ( X2-X1. NE. 0. ) THEN RADTAN = ATAN( (Y2-Y1)/(X2-X1) ) IF ( X2. LT. X1 ) RADTAN = RADTAN + PAI ELSE IF ( Y2-Y1. GT. 0. ) RADTAN = 0.5*PAI IF ( Y2-Y1. LT. 0. ) RADTAN = 1.5*PAI IF ( Y2-Y1. EQ. 0. ) RADTAN = 0. ENDIF IF ( RADTAN. LT. 0. ) RADTAN = RADTAN + 2.*PAI RETURN END C C 直線の読み込み、及び辺、節点ﾃﾞｰﾀの作成 SUBROUTINE SENDAT( SEN1, SEN2, ISENC, TEN, IHEN, IHENC + ,ITENS, XX, NN, IHENGD, IHENTD, ITENSD, IT, IH, IHENK + ,NODK, NODP, NDP, NODE, XXRR, YYUU, GOSA ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) CHARACTER MOJ2*2 DIMENSION SEN1(2,IHENGD), SEN2(2,IHENGD), ISENC(IHENGD) DIMENSION TEN(2,NODE), IHEN(2,IHENTD), IHENC(IHENTD) DIMENSION ITENS(ITENSD), XX(ITENSD), NN(ITENSD) DIMENSION NODK(NODE), NODP(NDP,NODE) C C JWWINの作図ﾃﾞｰﾀ読み込み C XXRR =-1.E10 YYUU =-1.E10 IK = 0 LC = 1 1 READ(2,'(A2)',END=2) MOJ2 IF ( MOJ2. EQ. 'lc' ) THEN BACKSPACE 2 READ(2,'(A2,I1)') MOJ2, LC GO TO 1 ENDIF IF ( MOJ2(1:1). NE. ' ' ) GO TO 1 BACKSPACE 2 IK = IK + 1 READ(2,*) SEN1(1,IK), SEN1(2,IK), SEN2(1,IK), SEN2(2,IK) XXRR = MAX( SEN1(1,IK), SEN2(1,IK), XXRR ) YYUU = MAX( SEN1(2,IK), SEN2(2,IK), YYUU ) ISENC(IK) = LC GO TO 1 2 CONTINUE KSEN = IK C C @基準線をＸ軸へ座標変換して各線との交点(Y=0)を求める C A交点をＸ座標昇順に並べ替え元座標へ変換して節点配列TENへ格納 C 　（同じ座標の節点がある場合は既往節点番号をITENSに記録） C B辺データを作成 C IT = 0 IH = 0 CALL INTSET ( IHENC, 1, IHENTD, 2 ) C @基準線をＸ軸へ座標変換して各線との交点(Y=0)を求める DO 10 I = 1, KSEN X0 = SEN1(1,I) Y0 = SEN1(2,I) AN = RADTAN( X0,Y0,SEN2(1,I),SEN2(2,I) ) C = COS(0.0-AN) S = SIN(0.0-AN) XR = (SEN2(1,I)-X0)*C - (SEN2(2,I)-Y0)*S II = 0 DO 11 J = 1, KSEN IF ( J. EQ. I ) GO TO 11 X1 = (SEN1(1,J)-X0)*C - (SEN1(2,J)-Y0)*S Y1 = (SEN1(1,J)-X0)*S + (SEN1(2,J)-Y0)*C X2 = (SEN2(1,J)-X0)*C - (SEN2(2,J)-Y0)*S Y2 = (SEN2(1,J)-X0)*S + (SEN2(2,J)-Y0)*C IF ( MAX(X1,X2). LT. 0.-GOSA ) GO TO 11 IF ( MIN(X1,X2). GT. XR+GOSA ) GO TO 11 IF ( MAX(Y1,Y2). LT. 0.-GOSA ) GO TO 11 IF ( MIN(Y1,Y2). GT. 0.+GOSA ) GO TO 11 IF ( ABS(Y2-Y1). LT. GOSA ) GO TO 11 XXX = X1 + (0.-Y1)*(X2-X1)/(Y2-Y1) IF ( XXX. LT. 0.-GOSA. OR. XXX. GT. XR+GOSA ) GO TO 11 II = II + 1 IF ( ABS(XXX). LT. GOSA ) XXX =0. IF ( ABS(XXX-XR). LT. GOSA ) XXX =XR XX(II) = XXX 11 CONTINUE C A交点をＸ座標昇順に並べ替え元座標へ変換して節点配列TENへ格納 CALL MINJUN ( XX, NN, II ) C = COS(AN) S = SIN(AN) XXK0 = -1.E10 IL = 0 IHH = IH DO 12 K = 1, II IF ( XX(NN(K))-XXK0. LT. GOSA ) GO TO 12 XXK0 = XX(NN(K)) IL = IL+1 X1 = X0 + XX(NN(K))*C Y1 = Y0 + XX(NN(K))*S NCK = NODCHK ( X1,Y1,TEN,IT, GOSA ) IF ( NCK. EQ. 0 ) THEN IT = IT+1 ITENS(IL)=IT TEN(1,IT)=X1 TEN(2,IT)=Y1 ELSE ITENS(IL)=NCK ENDIF C B辺データを作成 IF ( IL. EQ. 1 ) GO TO 12 IH = IH + 1 J1 = ITENS(IL-1) J2 = ITENS(IL ) IHEN(1,IH) = MAX ( J1, J2 ) IHEN(2,IH) = MIN ( J1, J2 ) DO 15 IIH = 1, IHH IF ( IHEN(1,IIH). EQ. IHEN(1,IH). AND. + IHEN(2,IIH). EQ. IHEN(2,IH) ) THEN IF ( ISENC(I). EQ. 6 ) IHENC(IIH) = 1 IH = IH-1 GO TO 12 ENDIF 15 CONTINUE IF ( ISENC(I). EQ. 6 ) IHENC(IH) = 1 12 CONTINUE 10 CONTINUE C C 節点から分岐する辺の数と、分岐先節点番号を登録 C INTSET;　配列NODKの1〜IT番に0を入れる C NODK(IHEN(1,I));I番目の辺の片方の節点から分岐するヶ数 C NODP(II,JJ);JJ番の節点からII番目に分岐する先の節点番号 C IHENK；辺の総使用可能回数 C CALL INTSET ( NODK, 1, IT, 0 ) IHENK = 0 DO 20 I = 1, IH NODK(IHEN(1,I)) = NODK(IHEN(1,I)) + 1 NODP( NODK(IHEN(1,I)) , IHEN(1,I) ) = IHEN(2,I) NODK(IHEN(2,I)) = NODK(IHEN(2,I)) + 1 NODP( NODK(IHEN(2,I)) , IHEN(2,I) ) = IHEN(1,I) IHENK = IHENK + IHENC(I) 20 CONTINUE RETURN END C C Ｘ、Ｙと同じ座標の節点がなければNODECHK=0、あれば=節点番号 FUNCTION NODCHK ( X,Y,TEN,IT, GOSA ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION TEN(2,IT) NODCHK=0 DO 10 I = 1, IT IF ( SENL(X,Y,TEN(1,I),TEN(2,I)). LT. GOSA ) THEN NODCHK=I RETURN ENDIF 10 CONTINUE RETURN END C C 整数配列のN1番以後をＫに初期化する SUBROUTINE INTSET ( M, N1, N, K ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION M(N) DO 10 I = N1, N M(I) = K 10 CONTINUE RETURN END C C 閉合点列の面積と重心を求める SUBROUTINE GYUSIN (X, Y, N, GX, GY, A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(1),Y(1) AX=0. AY=0. RX=0. RY=0. DO 10 I=1,N-1 AX = AX+(X(I)+X(I+1))*(Y(I+1)-Y(I))/2 RX = RX+X(I)*(Y(I+1)-Y(I))*X(I)/2 + +(X(I+1)-X(I))*(Y(I+1)-Y(I))/2*(X(I)+(X(I+1)-X(I))/3) AY = AY+(Y(I)+Y(I+1))*(X(I+1)-X(I))/2 RY = RY+Y(I)*(X(I+1)-X(I))*Y(I)/2 + +(Y(I+1)-Y(I))*(X(I+1)-X(I))/2*(Y(I)+(Y(I+1)-Y(I))/3) 10 CONTINUE A = ABS(AX) IF ( A. LT. 1.E-20 ) THEN GX = X(1) GY = Y(1) RETURN ENDIF GX = RX/AX GY = RY/AY RETURN END C C 点が閉合点列の内(=1)か外(=0)か線上(=-1)かを判定 C （点から水平に左側に引いた直線と点列の交点数が奇数なら内） C （点と同じ高さに水平な辺がある場合は回転してやり直す） FUNCTION INOUTO ( X , Y , N , XPP , YPP, GS ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION X(N), Y(N) C T = -1.484 DO 100 III = 1, 5 INOUTO = 0 T = T + 1.484 CALL KAITEN ( XPP, YPP, T, XP, YP ) DO 10 I = 1, N-1 J = I + 1 K = J + 1 IF ( K. GT. N ) K = K - N + 1 C CALL KAITEN ( X(I) , Y(I), T, XI, YI ) CALL KAITEN ( X(J) , Y(J), T, XJ, YJ ) CALL KAITEN ( X(K) , Y(K), T, XK, YK ) C DX = ABS( XJ - XI ) DY = ABS( YJ - YI ) XMI = MIN( XI, XJ ) IF ( DX + DY. LE. GS ) GO TO 10 IF ( XMI. GT. XP ) GO TO 10 C IF ( DY. LE. GS ) THEN IF ( ZEROMI(YI,YP,GS). NE. 0.0 ) GO TO 10 IF ( ZEROMI(XP,XI,GS)*ZEROMI(XP,XJ,GS).LE.0. ) THEN INOUTO = -1 RETURN ELSE GO TO 100 ENDIF ENDIF C PAY = ZEROMI(YP,YI,GS)*ZEROMI(YP,YJ,GS) IF ( PAY.GT.0. ) GO TO 10 IF ( PAY.LT.0. ) THEN XX = XI+(YP-YI)*(XJ-XI)/(YJ-YI) PAX = ZEROMI(XX,XP,GS) IF ( PAX. GT. 0. ) GO TO 10 IF ( PAX. LT. 0. ) THEN INOUTO = INOUTO + 1 GO TO 10 ELSE INOUTO = -1 RETURN ENDIF ENDIF IF ( PAY.EQ.0. AND. ZEROMI(YP,YJ,GS). EQ. 0. ) THEN IF (ZEROMI(YP,YI,GS)*ZEROMI(YP,YK,GS). LT. 0. ) THEN INOUTO = INOUTO + 1 ENDIF ENDIF 10 CONTINUE INOUTO = MOD ( INOUTO, 2 ) RETURN 100 CONTINUE RETURN END C C Ｔラジアンで座標を回転させる SUBROUTINE KAITEN ( X, Y, T, XX, YY ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) C = COS ( T ) S = SIN ( T ) XX = X * C + Y * S YY = Y * C - X * S RETURN END C C Ｘ１−Ｘ２を返すが、２点が誤差ＧＳ以下なら０とする。 FUNCTION ZEROMI( X1, X2, GS ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) ZEROMI = 0.0 A = ABS(X1-X2) IF ( A. LE. GS ) RETURN ZEROMI = X1- X2 RETURN END C C 囲みの探索（節点と辺のデータから要素データを作成する） SUBROUTINE KAKOMI( TEN, IHEN, IHENC, IHENTD, IT, IH, IHENK + , NODK, NODP, NDP, NODE, NLN, N, M, XEL, YEL, NLN1 + , NELM, IE, ITES, XYG, NOD, IHEC0, GOSA ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION TEN(2,NODE), IHEN(2,IHENTD), IHENC(IHENTD) DIMENSION NODK(NODE), NODP(NDP,NODE), NOD(0:NLN,NELM) DIMENSION XYG(2,NELM), N(0:NLN), M(NLN) DIMENSION XEL(NLN1), YEL(NLN1), IHEC0(IHENTD) C ITESKAI=500000000までの計算時間は数分程度 DATA ITESKAI / 500000000 / C C TEN；節点XY座標、IHEN；辺の両端節点番号、IHENC；辺の使用可能回数 C IT；節点数、IH；辺数、IHEN；辺の総使用可能回数 C IE；要素数、ITES；主な計算回数、XYG；要素重心座標 C NOD；要素構成節点数と番号、GOSA；同じ点と判定する誤差(mm) C C DO 10 三角形から順にNKAI角形までの囲みを探索する C DO 11 節点１から最終節点ITまで順に始点（位K=0)をずらす C DO 12 位をK+1からNKAKU（最高位）まで進める C DO 13 各位毎、前回分岐順+1の分岐順（辺）となる節点を設定 C N(IK)は、IK位の節点番号 C M(IK)は、IK位の節点が前の位の節点から何番目の分岐順かを示す C NODK(N(IK-1))は、前の位の節点から出る辺数（分岐数） C NODP(II、N(K-1))）は、前の位の節点からII番目の分岐先節点番号 C NKHAN(N,IK)は検討中の点列内に同じ節点がないか調べている C IHEMCは検討中の辺の使用可能回数がゼロでないかを調べている C CALL INCOPY ( IHENC, IHEC0, IH ) NKAI = NLN IF ( NLN. GT. IT ) NKAI = IT ITES = 0 IE = 0 DO 10 NKAKU = 3, NKAI DO 11 IP = 1, IT K = 0 N(K) = IP CALL INTSET( M, 1, NKAI, 0 ) 111 CONTINUE DO 12 IK = K+1, NKAKU MIK = M(IK) DO 13 II = MIK+1, NODK(N(IK-1)) ITES = ITES+1 IF ( ITES. GT. ITESKAI ) THEN WRITE(*,*) ' ERROR ITES=',ITES, ' PLEASE INPUT SOMEVALUE' READ(*,*) AIAI STOP ENDIF M(IK) = II N(IK) = NODP(II,N(IK-1)) IF ( NKHAN(N,IK). NE. 0 ) GO TO 13 IF ( IHEMC(N(IK),N(IK-1),IHEN,IHENC,IH). EQ. 0 ) GO TO 13 IF ( IK. NE. NKAKU ) THEN GO TO 12 ELSE IF ( N(IK). NE. IP ) THEN GO TO 13 ELSE GO TO 112 ENDIF ENDIF 13 CONTINUE IF ( IK. EQ. 1 ) THEN GO TO 11 ELSE K = IK-2 M(IK) = 0 GO TO 111 ENDIF 12 CONTINUE C 112 CONTINUE C 以下は囲み(NKAKU角形)が閉じた時の処理 C （囲みの中に辺や節点がないか判定した後、要素の登録を行う） C INOUTOは、その点が囲み内なら=1、囲み外なら=0 C IHELCは、囲みの各辺の使用可能回数を１減らしている IE = IE + 1 NOD(0,IE) = NKAKU XL = 1.E10 XR = -1.E10 YD = 1.E10 YU = -1.E10 IHSM = 0 DO 51 IEE = 1, NKAKU NOD(IEE,IE) = N(IEE-1) XEL(IEE) = TEN(1,NOD(IEE,IE)) YEL(IEE) = TEN(2,NOD(IEE,IE)) IF ( XEL(IEE). LT. XL ) XL = XEL(IEE) IF ( XEL(IEE). GT. XR ) XR = XEL(IEE) IF ( YEL(IEE). LT. YD ) YD = YEL(IEE) IF ( YEL(IEE). GT. YU ) YU = YEL(IEE) IHSM = IHSM + IHEMC(N(IEE),N(IEE-1),IHEN,IHEC0,IH) 51 CONTINUE C 囲みの辺が全て線色６（空洞用）でないか判定 IF ( IHSM. LE. NKAKU ) THEN IE = IE -1 M(NKAKU) = 0 K = NKAKU-2 GO TO 111 ENDIF C 囲みの中に既往の要素の辺の中点がないか判定 XEL(NKAKU+1) = XEL(1) YEL(NKAKU+1) = YEL(1) CALL GYUSIN ( XEL,YEL,NKAKU+1,XYG(1,IE),XYG(2,IE), AAA) DO 52 IEE = 1, IE-1 IF ( XYG(1,IEE). LT. XL. OR. XYG(1,IEE). GT. XR ) GO TO 52 IF ( XYG(2,IEE). LT. YD. OR. XYG(2,IEE). GT. YU ) GO TO 52 XX0 = TEN(1,NOD(NOD(0,IEE),IEE)) YY0 = TEN(2,NOD(NOD(0,IEE),IEE)) DO 53 IKK = 1, NOD(0,IEE) XX1 = TEN(1,NOD(IKK,IEE)) YY1 = TEN(2,NOD(IKK,IEE)) IHANT = INOUTO(XEL,YEL,NKAKU+1,(XX0+XX1)/2.,(YY0+YY1)/2.,GOSA ) IF ( IHANT. EQ. 0 ) GO TO 151 XX0 = XX1 YY0 = YY1 53 CONTINUE IE = IE -1 M(NKAKU) = 0 K = NKAKU-2 GO TO 111 151 CONTINUE 52 CONTINUE C 囲みの中に節点（交点）がないか判定 DO 54 ITT = 1, IT XX1 = TEN(1,ITT) YY1 = TEN(2,ITT) IF ( XX1. LT. XL. OR. XX1. GT. XR ) GO TO 54 IF ( YY1. LT. YD. OR. YY1. GT. YU ) GO TO 54 DO 55 IKK = 1, NOD(0,IE) IF ( ITT. EQ. NOD(IKK,IE) ) GO TO 54 55 CONTINUE IHANT = INOUTO (XEL, YEL, NKAKU+1, XX1, YY1, GOSA ) IF ( IHANT. EQ. 1 ) THEN IE = IE -1 M(NKAKU) = 0 K = NKAKU-2 GO TO 111 ENDIF 54 CONTINUE DO 56 IEE = 1, NKAKU CALL IHELC ( N(IEE-1),N(IEE),IHEN,IHENC,IHENK,IH ) 56 CONTINUE C 使用可能な辺がもう無い場合終了 IF ( IHENK. EQ. 0 ) GO TO 777 C １位(K=1)または２位(K=2)の節点から再実施 IF ( IHEMC(N(0),N(1),IHEN,IHENC,IH). EQ. 0 ) THEN K = 0 ELSE K = 1 ENDIF CALL INTSET( M, K+2, NKAKU, 0 ) GO TO 111 11 CONTINUE 10 CONTINUE C 777 CONTINUE C RETURN END C C 実数配列Ｘ１の値をＸ２にコピーする。 SUBROUTINE DECOPY ( X1, X2, NDM ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION X1(1), X2(1) DO 10 I = 1, NDM X2(I) = X1(I) 10 CONTINUE RETURN END C C 整数配列Ｉ１の値をＩ２にコピーする。 SUBROUTINE INCOPY ( I1, I2, NDM ) IMPLICIT DOUBLE PRECISION ( A-H, O-Z ) DIMENSION I1(1), I2(1) DO 10 I = 1, NDM I2(I) = I1(I) 10 CONTINUE RETURN END