初期ひずみ

どうも応力解析の初期条件として初期ひずみは入れることができないようなのでひずみから応力を計算して、形式を整えて出力するサブルーチンを書いた。
最近、論文を書かなきゃいけないのに付録たるサブルーチンに没頭気味。

***  READ FROM .fil ***
C
      SUBROUTINE HKSMAIN
C
      INCLUDE 'aba_param.inc'
      CHARACTER FNAME*80,FNAMEO*80
C 最大レコード数
      PARAMETER(NMAX=1000000)
C 最大エレメント数
      PARAMETER(MAXE=1000000)
C キーのためのオフセット
      PARAMETER(KOFFSET=2)
C ヤング率[M L-1 T-2]
　　　PARAMETER(YOUNGR=210.0E+7)
C ポアソン比
      PARAMETER(POISON=0.3E+0)
      DIMENSION ARRAY(513),JRRAY(NPRECD,513),LRUNIT(2,1)
      DIMENSION PEX(MAXE),PEY(MAXE),PEZ(MAXE),PEA(MAXE),PES(MAXE)
      DIMENSION NELN(MAXE)
      EQUIVALENCE(ARRAY(1),JRRAY(1,1))
C
C 入力ファイル名と出力ファイル名
C 入力ファイル名では拡張子は省略
C 
      WRITE(6,*)'Enter the name of the input file (w/o extension):'
      READ(5,'(A)') FNAME
      WRITE(6,*)'Enter the name of the output file:'
      READ(5,'(A)') FNAME0      
C
      DO 10 I=1,MAXE
         PEX(I) = 0.D+00
         PEY(I) = 0.D+00
         PEZ(I) = 0.D+00
         PES(I) = 0.D+00
         PEA(I) = 0.D+00
         NELN(I) = 0
   10 CONTINUE
      NTHRES=0
C
*** READ FROM .fil ***
C 初期化
C 読み込む結果ファイルの数
      NRU = 1
C デバイス番号
      LRUNIT(1,1) = 8
C 結果ファイルがバイナリの場合2,アスキーの場合1
   44 CONTINUE
      WRITE(6,*)'Enter the format of input file (1-ASCII, 2-binary):'
      READ(5,*)LRUNIT(2,1)
      IF (LRUNIT(2,1).NE. 1 .AND. LRUNIT(2,1) .NE. 2) THEN
            WRITE(6,*) 'ERROR! This number must be 1 or 2'
            GOTO 44
      END IF
C 書き出し形式、読み込みの場合は何の意味も持たない
      LOUTF = 0
      CALL INITPF(FNAME,NRU,LRUNIT,LOUTF)
C デバイス番号
      JUNIT = 8
      CALL DBRNU(JUNIT)
C
      DO 100 J=1,NMAX
C レコードの順次読み出し
         CALL DBFILE(0,ARRAY,JRCD)
C EOFの検出
         IF (JRCD .NE. 0) GOTO 110
         KEY = JRRAY(1,KOFFSET)
C エレメント番号
         IF (KEY .EQ. 1) THEN
            NELNO = JRRAY(1,KOFFSET+1)
            NELN(NELNO) = NELNO
         END IF
C
C 塑性ひずみPE
         IF (KEY .EQ. 22) THEN
            PEX(NELNO) = ARRAY(KOFFSET+1)
            PEY(NELNO) = ARRAY(KOFFSET+2)
            PEZ(NELNO) = ARRAY(KOFFSET+3)
            PES(NELNO) = ARRAY(KOFFSET+4)
            PEA(NELNO) = ARRAY(KOFFSET+6)
         END IF
C
  100 CONTINUE
  110 CONTINUE
      KEND = J-1
      WRITE(6,*) 'KEND=',KEND
      WRITE(6,*) ' '
C
*** WRITE TO FILE ***
C
      OPEN(51,FILE=FNAMEO)
C
C
      WRITE(51,*)'**Initial conditions,TYPE=STRESS'
      DO K=1,MAXE
         IF (NELN(K) .NE. 0) THEN
            PEAL=PEX(K)+PEY(K)+PEZ(K)
            SX=YOUNGR/(1+POISON)*(PEX(K)-
     1       (POISON/(1.0-2.0*POISON))*(PEAL)))
            SY=YOUNGR/(1+POISON)*(PEY(K)-
     1       (POISON/(1.0-2.0*POISON))*(PEAL)))
            SZ=YOUNGR/(1+POISON)*(PEZ(K)-
     1       (POISON/(1.0-2.0*POISON))*(PEAL)))
            WRITE(51,*)NELN(K),',',SX,',',SY,',',SZ,',0.0,0.0,0.0'
         END IF
      END DO
C
      WRITE(6,*) ' '
      WRITE(51,*) ' '

      STOP
      END

あ、ちなみにここに書いてあるコードをコピペで使うような人はいないと思うけど一応いっておくと、宣伝条項なしの新しいほうのBSDライセンスに準拠ってことで。