有限要素法　　二次元非定常熱伝導解析

（このプログラムは2次元の熱伝導と熱伝達の解析が可能）

（ここでは、差分法と同じ問題を解きます）

 （表面熱伝達の計算も可能です）

1  solver       

     A.解析例（差分法と同じもの、温度空間分布）

2   FEM-hca-pre-processor-FORTRAN

　B.解析例（差分法と同じもの、時系列データ)

3   FEM-hca-pre-post-processor-Excel-VBA

　  input-data.txt     ****-post.txtとも開いてコピーします。

　コピーはCTRL+A、CTRL+Cです。

　そしてSheet3 を作り[A2]セルに貼り付け、

　コンマのあるデータでセルに格納します。

    マクロを実行すれば、図が描けます。

　animationだけは、Sheet1～Sheet10まで白紙にして

　用意して、Sheet11の[A2]セルに****-post.txtの

　ファイルの中身を貼り付けます。

　実行すると、Sheet1～Sheet10にコンタ図が描かれます。

　PowerPointに貼り付けて、スライドショウで

　1枚あたりの表示時間を0.5秒程度にすると

　アニメーションになります。　

　ここで用語の解説、

　マクロとVBAは同じものを指しています。

　マクロ(macro)の言語はVBA（Visual Basic Application)

　です。

 3-1.Excel-VBA-mesh

 3-2.Excel-VBA-node-number

 3-3.Excel-VBA-element-number

 3-4.Excel-VBA-number

 3-5.Excel-VBA-conta-animation

4 .animation 例（PowerPoint使用、3－4と組み合わせて）

5.FEM理論（Ritz法）

グラフFo＝0.00125、0.0125、0.0625、0.125、0.25の

作成方法は、****-post.txtのファイル（自分で名前をつける）

のファイルの先頭からnode,　elem,　nodeと来たら、

node番号1から101までをコピーしてExcelに貼り付け、

理論解をつける場合はマクロを自分で書き、

散布図を描けば、出来上がり。

グラフtime-seriesの作成方法は、

ファイル****-time-series.txt（自分で実行時名前をつける）

これを、開き、CTRL+A,CTRL+Cですべてコピーして

Excelに貼り付け、散布図を描くと出来上がりです。

ここに出ているプログラムをコピー＆ペイスト

すると、1行おきにブランク行が出来るが、

構わず、******.forに保存して、コンパイル、

実行出来ることを確認しました。

VBAもExcelのファイルを****.xlsmにして、

マクロのところに貼り付けて実行できます。

【1　solver】

　（以下、Cの書いてある所から下をコピーしてTeraPadに

　　　貼り付け、FEM-hca-solver.forという名称で保存）

C

C FEM unsteady heat conduction analysis

C

C-----------------------------------------------------------

C  program file name :  FEM-hca-solver.for

C

C  compile方法

C  >gfortran FEM-hca-solver.for -o FEM-hca-solver.exe

C

C  入力データの作成方法

C　2  pre-processor-program-FORTRANに書いてあります。

C

C　実行方法

C  >FEM-hca-solver

C   input file name ?

C   FEM-hca-input-data.txt

C   ・

C   ・

C

C

C

C----------------------------------------------------------

C 注意書き

C  このプログラムを使用して製品を製作しても、計算結果を使用 

C　し、製品に不備が生じても当方は責任を負いません。

C　このプログラムは、無料公開するものですが、これで、商売を

C　行う事は、謹んで頂きたい。

C

C  このプログラムは、学習用です。

C

C----------------------------------------------------------

C

C

C

C 

C--------------------------------------------------------------------

C     variable    :  mean

C

C      nall       :　ALL NODE NUMBER

C      ieall      :  ALL ELEMENT NUMBER

C      nq         :　ALL Heating Value NODE NUMBER

C      x(i)       :　NODE NUMBER i x coordinate

C      y(i)       :　NODE NUMBER i y coordinate

C      nelem(i,j) :　ELEMENT Compornent NODE NUMBER

C                         （i：ELEMENT NUMBER　j=1,2,3）

C      ramda(I)   :　material NUMBER I heat conductivity

C      h(I)       :　ELEMENT NUMBER I thickness

C      number_of_boundary      :　Temperature Fix ALL NODE NUMBER

C      number_ob_bound(I)      :　Temperature fix node number

C      bound(I)   :　node number I temperature fix temperature

C                                     （boundary condition)

C      iden       :  Heat transfer boundary number

C      ialpha(I,2):  Heat transfer boundary node number(I,1)(I,2)

C      alpha(I)   :  Heat transfer rate of I

C      troom      :  Ambient temperature

C      tinit      :  Initial Temperature

C      qt(I)      :　NODE NUMBER I Heating value

C      all(I,J)   :　ALL stiffness Matrix

C      ek(I,J)    :　Element stiffness Matrix

C      sk3(I,J)   :  element heat transfer matrix

C      c(I)       :  material NUMBER I specific heat (Ws/g℃）

C      roh(I)     :  material NUMBER I density (g/mm3）

C      CC(I,J)    :　Thermal capacity Matrix

C      AKCTP(I,J) :  Unsteady  θ(t+Δt)  calculation  Matrix

C      AKCT(I,J)  :  Unsteady  multiply θ(t) Matrix 

C      TEMP(I)    :  NODE NUMBER I Temperature:θ(t) (℃)

C      dt         :  Δt(sec)

C      time       :  Time(sec)

C      tend       :  END Time(sec)

C      JQ(I)      :　Heating NODE NUMBER（以下K)

C      QTIM(K,J)　:　NODE NUMBER K  Time J   （sec）

C      QT(K,J)    :  NODE NUMBER K   （W)（polygonal line）

C      NTE        :  write Time series data    all node number(max40）

C      NTEMP(I)   :　write Time series data   node number

C      IPTIME(I)  :  OUTPUT CONTA FIG IPTIME(I=10)

C                 :         DT*IPTIME(I)=OUTPUT TIME

C

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /GGGGG/ AKCTP(5000,5001),AKCT(5000,5000)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /MMMMM/ CC(5000,5000)

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      common /PPPPP/ number_of_q(5000),q(5000)

C

      call open_file

C

C     add

C     ****************************

      call current_time('FEM unsteady heat conduction analysis start')

C     ****************************

C     add end

C

C

      call file_read(title)

C

      call file_write(title)

C

      call pre_processor

C

C

      write(9,1740)(NTEMP(i),i=1,NTE)

 1740 format('TIME,',100(i4,'-temp,'))

C

C---------------------------------------------

C      write(3,*)'ALL'

C      do 100 i=1,nall

C      write(3,200)(i,j,ALL(i,j),j=1,nall+1)

C  200 format(10('(',I4,',',I4,')=',e10.4))

C  100 continue

C---------------------------------------------

C

C

       do 1450 i=1,nall

      TEMP(i)=tinit

 1450 continue

C

C Calculation_UNSTEADY_START

C

      itime_end=int(tend/dt)+3

      write(6,*)'dt=',dt

      write(6,*)('IP',i,IPTIME(i),i=1,10)

C

C

C     do  start

C

      do 1000 itime=1,itime_end

      time=FLOAT(itime)*dt

      write(6,*)'TIME=',time

      do 1500 i=1,nall

      do 1510 j=1,nall+1

      ALL(i,j)=0.0

 1510 continue

 1500 continue

C

C   Left side substitution /// 左辺代入

C

C

      do 400 i=1,nall

      do 410 j=1,nall

      ALL(i,j)=AKCTP(i,j)

  410 continue

  400 continue

C

C    Left side substitution end /// 左辺代入終

C

C    Right side substitution start   ///  右辺代入始め

C

C

      do 420 i=1,nall

      do 430 j=1,nall

      ALL(i,nall+1)=ALL(i,nall+1)+AKCT(i,j)*TEMP(j)

  430 continue

  420 continue

C

C     Right side substitution end   ///   右辺代入終わり

C

      call right_side_substitution(time)

C

      call set_up_boundary_condition

C

      call gauss(ALL,5000,5001,nall,ierr)

      if(ierr.eq.0)then

      else

      write(3,*)'ierr=',ierr

      write(6,*)'ierr=',ierr

      stop

      endif

C

C

      do 1010 i=1,10

      if(IPTIME(i).eq.itime)then

      do 1020 j=1,nall

      TEMPTI(j,i)=ALL(j,nall+1)

 1020 continue

      else

      endif

 1010 continue

C

      do 1030 i=1,nall

      TEMP(i)=ALL(i,nall+1)

 1030 continue

C

C

      write(9,1730)time,(TEMP(NTEMP(i)),i=1,NTE)

 1730 format(e16.8,100(',',e16.8))

 1000 continue

C

C Calculation_UNSTEADY_end

C

C

      call write_post

C

      call result_by_all

C

      call current_time('FEM usteady heatconduction analysis stop')

C

C

      call close_file

C

      stop

      end

C

C--------------------------------------------------------------------

C

      subroutine open_file

      character file_name*120

C

      write(6,*)'input_file_name'

      read(5,100)file_name

  100 format(A)

      open(UNIT=2,file=file_name,STATUS='OLD')

C

      write(6,*)'output_file_name'

      read(5,100)file_name

      open(UNIT=3,file=file_name,STATUS='UNKNOWN')

C     

      write(6,*)'post_file_name'

      read(5,100)file_name

      open(UNIT=8,file=file_name,STATUS='UNKNOWN')

C

      write(6,*)'time_series_file_name'

      read(5,100)file_name

      open(UNIT=9,file=file_name,STATUS='UNKNOWN')

C

C

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine file_read(title)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /GGGGG/ AKCTP(5000,5001),AKCT(5000,5000)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      common /PPPPP/ number_of_q(5000),q(5000)

      character INDEX*8

      write(6,*)'file_read_IN'

      read(2,*)nall,ieall,number_of_boundary,iden,troom,nq,

     &         dt,iprint,tend,tinit,mat_vareity

      write(6,*)'nall=',nall,'ieall=',ieall

      read(2,100)INDEX

      write(6,100)INDEX

  100 format(A)

      do 200 i=1,nall

      read(2,*)j,x(j),y(j)

  200 continue

      read(2,100)INDEX

      do 210 i=1,ieall

      read(2,*)k,(nelem(i,j),j=1,3),material_number(i),h(i)

  210 continue

      read(2,100)INDEX

      do 500 i=1,mat_vareity

      read(2,*)k,ramda(k),c(k),roh(k)

  500 continue

      if(number_of_boundary.ne.0)then

      read(2,100)INDEX

      do 220 i=1,number_of_boundary

      read(2,*)number_of_bound(i),bound(number_of_bound(i))

  220 continue

      else

      endif

      if(iden.ne.0)then

      read(2,100)INDEX

      do 230 i=1,iden

      read(2,*)ialpha(i,1),ialpha(i,2),alpha(i)

  230 continue

      else

      endif

      if(nq.ne.0)then

      read(2,100)INDEX

C

C

      do 240 L=1,nq

      read(2,*)JQ(L),(QTIM(JQ(L),k),QT(JQ(L),k),k=1,11)

  240 continue

C

C

      else

      endif

      read(2,100)INDEX

      write(6,100)INDEX

      read(2,*)NTE

      read(2,*)(NTEMP(i),i=1,NTE)

      read(2,100)INDEX

      write(6,100)INDEX

      read(2,*)(IPTIME(i),i=1,10)

      read(2,100)title

C

C

      if(nq.ne.0)then

      do 2204 L=1,nq

      do 2208 i=1,9

      T1=QTIM(JQ(L),i)

      T2=QTIM(JQ(L),i+1)

      Q1=QT(JQ(L),i)

      Q2=QT(JQ(L),i+1)

      A1(JQ(L),i)=(Q2-Q1)/(T2-T1)

      if((T2-T1).eq.0.)A1(JQ(L),i)=0.0

      B1(JQ(L),i)=Q1-A1(JQ(L),i)*T1

 2208 continue

 2204 continue

      write(3,*)'Heating Value  Poligonal line'

      do 3602 L=1,nq

      write(3,3600)JQ(L),(A1(JQ(L),i),B1(JQ(L),i),i=1,10)

 3600 format('Heating Value node number=',i4,2x,

     &  10(2x,'A1=',G14.6,1x,'B1=',G14.6))

 3602 continue

      else

      endif

C

C

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine pre_processor

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /GGGGG/ AKCTP(5000,5001),AKCT(5000,5000)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /MMMMM/ CC(5000,5000)

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      dimension ek(3,3)

      do 1000 ielem=1,ieall

      call make_element_stiffness_matrix(ek)

      call make_all_stiffness_matrix(ek)

 1000 continue

      if(iden.ne.0)then

      call set_up_heat_transfer_condition()

      else

      endif

C

      if(nq.ne.0)then

      call heat_value

      else

      endif

      if(number_of_boundary.ne.0)then

      call set_up_boundary_condition

      else

      endif

C

C

C

      call creation_CC

C

C

      do 2000 i=1,5000

      do 2100 j=1,5000

      ak(i,j)=0.0

 2100 continue

 2000 continue

      do 2200 i=1,nall

      do 2300 j=1,nall

      ak(i,j)=ALL(i,j)

 2300 continue

 2200 continue

C    AKCTP          AKCTP=1/2*ak+1/dt*CC

C

      do 2400 i=1,nall

      do 2500 j=1,nall

      AKCTP(i,j)=ak(i,j)/2.0+CC(i,j)/dt

 2500 continue

 2400 continue

C

C    AKCT            AKCT=-1/2*ak+1/dt*CC

C

      do 2410 i=1,nall

      do 2510 j=1,nall

      AKCT(i,j)=-ak(i,j)/2.0+CC(i,j)/dt

 2510 continue

 2410 continue

C

C

      do 400 i=1,5000

      IQT(i)=1

  400 continue

C

C

      return

      end

C

C-----------------------------------------------------------

C

      subroutine creation_CC

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /MMMMM/ CC(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      dimension CD(3,3)

      do 10 i=1,5000

      do 20 j=1,5000

      CC(i,j)=0.0

   20 continue

   10 continue

      do 30 ielem=1,ieall

      do 40 i=1,3

      do 50 j=1,3

      CD(i,j)=0.0

   50 continue

   40 continue

      CD(1,1)=1.0

      CD(2,2)=1.0

      CD(3,3)=1.0

      call calculation_s

      mat_n=material_number(ielem)

      do 60 i=1,3

      do 70 j=1,3

      CD(i,j)=CD(i,j)*c(mat_n)*roh(mat_n)*s*h(ielem)/3.0

   70 continue

   60 continue

C

C

C

      do 80 i=1,3

      do 90 j=1,3

      i1=nelem(ielem,i)

      j1=nelem(ielem,j)

      CC(i1,j1)=CC(i1,j1)+CD(i,j)

   90 continue

   80 continue

   30 continue

      return

      end

C

C-----------------------------------------------------------

C

      subroutine heat_value

      common /FFFFF/ ALL(5000,5001)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /PPPPP/ number_of_q(5000),q(5000)

      if(nq.ne.0)then

      do 100 i=1,nq

      noq=number_of_q(i)

      ALL(noq,nall+1)=ALL(noq,nall+1)+q(noq)

  100 continue

      else

      endif

      return

      end

C

C-----------------------------------------------------------

C

      subroutine set_up_heat_transfer_condition()

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /FFFFF/ all(5000,5001)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      dimension sk3(2,2)

      do 10 k=1,iden

      i=ialpha(k,1)

      j=ialpha(k,2)

      call element_hit(i,j,iielem)

      xa=x(i)-x(j)

      ya=y(i)-y(j)

      sl3=sqrt(xa*xa+ya*ya)

      sk3(1,1)=2.0*alpha(k)*sl3*h(iielem)/6.0

      sk3(1,2)=alpha(k)*sl3*h(iielem)/6.0

      sk3(2,1)=sk3(1,2)

      sk3(2,2)=sk3(1,1)

      all(i,i)=all(i,i)+sk3(1,1)

      all(i,j)=all(i,j)+sk3(1,2)

      all(j,i)=all(j,i)+sk3(2,1)

      all(j,j)=all(j,j)+sk3(2,2)

      all(i,nall+1)=all(i,nall+1)+alpha(k)*troom*sl3*h(iielem)/2.0

      all(j,nall+1)=all(j,nall+1)+alpha(k)*troom*sl3*h(iielem)/2.0

C

C

   10 continue

      return

      end

C

C-----------------------------------------------------------------------

C

      subroutine element_hit(i,j,iielem)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      do 10 ii=1,ieall

      i1=nelem(ii,1)

      i2=nelem(ii,2)

      i3=nelem(ii,3)

      if(i.eq.i1.and.j.eq.i2)go to 20

      if(i.eq.i2.and.j.eq.i3)go to 20

      if(i.eq.i3.and.j.eq.i1)go to 20

      if(i.eq.i2.and.j.eq.i1)go to 20

      if(i.eq.i3.and.j.eq.i2)go to 20

      if(i.eq.i1.and.j.eq.i3)go to 20

   10 continue

      write(3,100)i,j

  100 format('It is not neighbor node number ',i4,2x,i4,

     &'Heat transfer Boundary')

      stop

   20 iielem=ii

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine make_element_stiffness_matrix(ek)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /GGGGG/ AKCTP(5000,5001),AKCT(5000,5000)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      dimension A(3,3),AT(3,3),B(2,3),BT(3,2),PM(3,2),PA(3,3),ek(3,3)

      call set_up_A_matrix(A,3,3)

C

C

      call set_up_B_matrix(B,2,3)

C

C

      call transpose_matrix(A,3,3,AT,3,3)

C

C

      call transpose_matrix(B,2,3,BT,3,2)

C

C

      call multiply_matrix(AT,3,3,BT,3,2,PM,3,2)

C

C

      call multiply_matrix(PM,3,2,B,2,3,PA,3,3)

C

C

      call multiply_matrix(PA,3,3,A,3,3,ek,3,3)

C

C

      call calculation_s

      do 1000 i=1,3

      do 2000 j=1,3

      mat_n=material_number(ielem)

      ek(i,j)=ek(i,j)*ramda(mat_n)*s*h(ielem)

 2000 continue

 1000 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine set_up_A_matrix(A,n1,n2)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /GGGGG/ AKCTP(5000,5001),AKCT(5000,5000)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      dimension A(n1,n2)

      i=nelem(ielem,1)

      j=nelem(ielem,2)

      k=nelem(ielem,3)

C

C

      A(1,1)=x(j)*y(k)-x(k)*y(j)

      A(1,2)=x(k)*y(i)-x(i)*y(k)

      A(1,3)=x(i)*y(j)-x(j)*y(i)

      A(2,1)=y(j)-y(k)

      A(2,2)=y(k)-y(i)

      A(2,3)=y(i)-y(j)

      A(3,1)=x(k)-x(j)

      A(3,2)=x(i)-x(k)

      A(3,3)=x(j)-x(i)

      call calculation_s

      do 10 i=1,3

      do 20 j=1,3

      A(i,j)=A(i,j)*(1.0/(2.0*s))

   20 continue

   10 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine calculation_s

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /GGGGG/ AKCTP(5000,5001),AKCT(5000,5000)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      i=nelem(ielem,1)

      j=nelem(ielem,2)

      k=nelem(ielem,3)

      s=(x(j)-x(i))*(y(k)-y(i))-(x(k)-x(i))*(y(j)-y(i))

      s=abs(s)*0.5

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine set_up_B_matrix(B,n1,n2)

      dimension B(n1,n2)

C

C

      B(1,1)=0.0

      B(1,2)=1.0

      B(1,3)=0.0

      B(2,1)=0.0

      B(2,2)=0.0

      B(2,3)=1.0

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine transpose_matrix(tm1,n1,n2,tm2,m1,m2)

      dimension tm1(n1,n2),tm2(m1,m2)

      do 20 i=1,n1

      do 10 j=1,n2

      tm2(j,i)=tm1(i,j)

   10 continue

   20 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine multiply_matrix(A,n1,n2,B,m1,m2,C,k1,k2)

      dimension A(n1,n2),B(m1,m2),C(k1,k2)

      do 30 i=1,k1

      do 20 j=1,k2

      C(i,j)=0.0

   20 continue

   30 continue

      do 50 i=1,n1

      do 60 j=1,m2

      do 70 k=1,n2

      C(i,j)=C(i,j)+A(i,k)*B(k,j)

   70 continue

   60 continue

   50 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine make_all_stiffness_matrix(ek)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /FFFFF/ ALL(5000,5001)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      dimension ek(3,3)

      do 10 i=1,3

      do 20 j=1,3

      ii=nelem(ielem,i)

      jj=nelem(ielem,j)

      ALL(ii,jj)=ALL(ii,jj)+ek(i,j)

   20 continue

   10 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine set_up_boundary_condition

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /FFFFF/ ALL(5000,5001)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      if(number_of_boundary.ne.0)then

      do 10 i=1,number_of_boundary

      nb=number_of_bound(i)

      do 20 j=1,nall

      ALL(nb,j)=0.0

   20 continue

      ALL(nb,nb)=1.0

      ALL(nb,nall+1)=bound(nb)

   10 continue

      else

      endif

      return

      end

C

C--------------------------------------------------------------------

C

C

C                                       

      subroutine right_side_substitution(time)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /FFFFF/ all(5000,5001)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      dimension sk3(2,2)

      if(iden.ne.0)then

      do 10 k=1,iden

      i=ialpha(k,1)

      j=ialpha(k,2)

      call element_hit(i,j,iielem)

      xa=x(i)-x(j)

      ya=y(i)-y(j)

      sl3=sqrt(xa*xa+ya*ya)

      sk3(1,1)=2.0*alpha(k)*sl3*h(iielem)/6.0

      sk3(1,2)=alpha(k)*sl3*h(iielem)/6.0

      sk3(2,1)=sk3(1,2)

      sk3(2,2)=sk3(1,1)

C      all(i,i)=all(i,i)+sk3(1,1)

C      all(i,j)=all(i,j)+sk3(1,2)

C      all(j,i)=all(j,i)+sk3(2,1)

C      all(j,j)=all(j,j)+sk3(2,2)

      all(i,nall+1)=all(i,nall+1)+alpha(k)*troom*sl3*h(iielem)/2.0

      all(j,nall+1)=all(j,nall+1)+alpha(k)*troom*sl3*h(iielem)/2.0

C

C

   10 continue

      else

      endif

      if(nq.ne.0)then

      do 20 i=1,nq

      IF(time.GE.QTIM(JQ(i),IQT(JQ(i))+1))IQT(JQ(i))=IQT(JQ(i))+1

   20 CONTINUE

      DO 1210 L=1,NQ

      QX=A1(JQ(L),IQT(JQ(L)))*time+B1(JQ(L),IQT(JQ(L)))

C      write(6,*)'nall=',nall,'QX=',QX

C      write(3,*)'time=',time,'QX=',QX,'nall=',nall

      all(JQ(L),nall+1)=all(JQ(L),nall+1)+QX

 1210 CONTINUE

      else

      endif

      return

      end

C

C--------------------------------------------------------------------

C

C

      subroutine file_write(title)

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      common /PPPPP/ number_of_q(5000),q(5000)

C     confirmation input data

C

      write(3,100)title

  100 format(A)

      write(3,300)nall,ieall,number_of_boundary,nq,iden,

     &         mat_vareity,troom

  300 format('numer of nude    =',I5/

     &       'number of element=',I5/

     &       'boundary condition No.=',I5/

     &       'heating node total No.=',I5/

     &       'heat transfer   IDEN  =',I5/

     &       'material vareity      =',I5/

     &       'ambient temperature   =',F8.2)

      write(3,*)'NODE'

      do 320 i=1,nall

      write(3,310)i,x(i),y(i)

  310 format(I5,2G16.8)

  320 continue

      write(3,*)'ELEMENT'

      write(3,330)

  330 format('element number     i     j      k     ',

     &       ' material number   h')

      do 340 i=1,ieall

      write(3,350)i,(nelem(i,j),j=1,3),material_number(i),h(i)

  340 continue

  350 format(8X,I5,3I6,8X,I3,2X,G16.6)

      if(number_of_boundary.ne.0)then

      write(3,*)'Temperature Fix Boundary Condition'

      write(3,*)'NODE number    Fix Temperature'

      do 360 i=1,number_of_boundary

      write(3,370)number_of_bound(i),bound(number_of_bound(i))

  360 continue

  370 format(I10,3X,G14.6)

      else

      endif

      if(iden.ne.0)then

      write(3,380)

  380 format('Heat transfer condition'/

     &       'Heat transfer node number   Heat transfer ratio'/

     &       '        I      J')

      do 390 i=1,iden

      write(3,400)ialpha(i,1),ialpha(i,2),alpha(i)

  400 format(5X,2I5,4X,G16.8)

  390 continue

      else

      endif

      if(nq.ne.0)then

      write(3,*)'Heating node number      heating value'

      do 430 i=1,nq

      L=number_of_q(i)

      write(3,440)L,(QTIM(L,j),QT(L,j),j=1,11)

  430 continue

  440 format(10X,I5,11(3X,G16.8,3X,G16.8))

      else

      endif

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine result_by_all

      common /FFFFF/ ALL(5000,5001)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      write(3,*)'Result_New'

      write(3,100)(i,ALL(i,nall+1),i=1,nall)

  100 format(10(i4,'=',e14.6))

      return

      end

C     

C---------------------------------------------------------------

C     gauss    : 連立方程式を解く / solve simultaneous equations

C     

C     A        : hairetsumei 配列名 / Array name

C     n1,n2    : hairetsu ookisa / Array size

C     n        : michisui no kazu / (変数の数) / Number of variables

C     ierr     : error judgement 0=normal  1 or 2=error

C     A(i,n+1) : answer(i=1 to n)

C

C

C

      subroutine gauss(A,n1,n2,n,ierr)

      dimension A(n1,n2)

      ierr=0

      do 10 k=1,n-1

      Lmax=k

      do 20 i=k,n

      if(abs(A(i,k)).GT.abs(A(Lmax,k)))then

      Lmax=i

      else

      endif

   20 continue

      do 30 i=1,n+1

      b=A(Lmax,i)

      A(Lmax,i)=A(k,i)

      A(k,i)=b

   30 continue

      do 40 i=k+1,n

      if(A(k,k).eq.0.0)then

      ierr=1

      return

      else

      endif

      C=A(i,k)/A(k,k)

      do 50 j=1,n+1

      A(i,j)=A(i,j)-A(k,j)*C

   50 continue

   40 continue

   10 continue

      if(A(n,n).eq.0.0)then

      ierr=2

      return

      else

      endif

      A(n,n+1)=A(n,n+1)/A(n,n)

      do 60 i=2,n

      j=n+1-i

      do 70 k=j+1,n

      A(j,n+1)=A(j,n+1)-A(j,k)*A(k,n+1)

   70 continue

      A(j,n+1)=A(j,n+1)/A(j,j)

   60 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine write_post

      common /AAAAA/ x(5000),y(5000),nelem(5000,3)

      common /BBBBB/ material_number(5000),h(5000)

      common /CCCCC/ bound(5000),number_of_bound(5000)

      common /DDDDD/ ialpha(5000,2),alpha(5000)

      common /EEEEE/ ramda(50),c(50),roh(50)

      common /FFFFF/ ALL(5000,5001)

      common /HHHHH/ A1(5000,11),B1(5000,11),QTIM(5000,11),QT(5000,11)

      common /IIIII/ NTEMP(100),IQT(5000),JQ(5000)

      common /JJJJJ/ TEMP(5000),TEMPTI(5000,10),IPTIME(10)

      common /KKKKK/ ak(5000,5000)

      common /NNNNN/ nall,ieall,number_of_boundary,iden,troom,nq,ielem

      common /OOOOO/ dt,time,iprint,tend,tinit,mat_vareity,NTE,s

      common /PPPPP/ number_of_q(5000),q(5000)

      write(8,3000)nall,ieall,1.0,1.0,1.0

 3000 format(I4,',',I4,3(',',E12.5))

      write(8,*)'NODE'

      do 3010 i=1,nall

      write(8,3020)i,x(i),y(i)

 3020 format(I4,',',E14.6,',',E14.6)

 3010 continue

      write(8,*)'ELEM'

      do 3030 i=1,ieall

      write(8,3040)i,(nelem(i,j),j=1,3)

 3040 format(3(I4,','),I4)

 3030 continue

      write(8,*)'NODE'

      do 3050 i=1,nall

      write(8,3060)i,(TEMPTI(i,j),j=1,10)

 3060 format(I4,10(',',E14.6))

 3050 continue

      return

      end

C

C--------------------------------------------------------------------

C

      subroutine close_file

      close(2)

      close(3)

      close(8)

      close(9)

      return

      end

C------------------------------------------------------------------

C

C

      subroutine current_time(c)

      implicit none

      character(*),intent(in)::c

      integer::ti(1:8)

      character(10)::b(1:3)

      call date_and_time(b(1),b(2),b(3),ti)

      write(3,100)c,ti(1),"/",ti(2),"/",ti(3)," ",ti(5),":",ti(6),":",

     &ti(7),",(yyyy/mm/dd hh:mm:ss)"

  100 format(A20,2X,I4,A1,I2,A1,I2,A1,I2,A1,I2,A1,I2,A22)

      return

      end subroutine current_time

【2  FEM-hca-pre-processor-FORTRAN】

　（以下、Cの書いてある所から下をコピーしてTeraPadに

　　　貼り付け、FEM-hca-preprocessor.forという名称で保存）

C  pre processer

C  file name FEM-hca-preprocessor.for

C-------------------------------------------------------------

C  4-FDM参照のこと。MinGWが必要。

C

C-------------------------------------------------------------

C  compile 方法

C  >gfortran FEM-hca-preprocessor.for -o FEM-hca-preprocessor.exe

C

C-------------------------------------------------------------

C  実行方法（FEM入力データ作成）

C  >FEM-hca-preprocesssor

C

C-------------------------------------------------------------

C  FEM入力データ　　FEM_HCA_input_data.txt

C-------------------------------------------------------------

C

C

C

C FEMunsteady-pre-processer-square-init-temp-100deg-both-sides-0deg.for

C

      CHARACTER FLOUT*80

      DIMENSION X(2500),Y(2500),NELM(2500,3)

      DIMENSION IALUFA(2500,2),ALUFA(2500)

      DIMENSION QTIM(2500,10),QT(2500,10),NNQ(2500)

      DIMENSION NTEMP(40),NB(2500),IPTIME(10),CKKAP(10)

      OPEN(UNIT=2,FILE='FEM_HCA_input_data.txt',

     &STATUS='UNKNOWN')

      K=0

      DO 10 J=1,3

      DO 12 I=1,101

      K=K+1

      X(K)=FLOAT(I-1)*0.001   !NODE NUMBER K X COORDINATE

      Y(K)=FLOAT(J-1)*0.001   !NODE NUMBER K Y COORDINATE

   12 CONTINUE

   10 CONTINUE

      NODE=K

      NELM(1,1)=2          !ELEMENT NUMBER 1 i=1  (NODE NUMBER)

      NELM(1,2)=102          !ELEMENT NUMBER 1 j=103(NODE NUMBER)

      NELM(1,3)=1        !ELEMENT NUMBER 1 k=102(NODE NUMBER)

      NELM(2,1)=2

      NELM(2,2)=103

      NELM(2,3)=102

      DO 30 I=3,200

      DO 32 J=1,3

      NELM(I,J)=NELM(I-2,J)+1 !ELEMENT NUMBER I 

   32 CONTINUE

   30 CONTINUE

      DO 40 I=201,400

      DO 42 J=1,3

      NELM(I,J)=NELM(I-200,J)+101

   42 CONTINUE

   40 CONTINUE

C      NQ=2               ! HEATING NODE TOTAL NUMBER

C      NNQ(1)=931         ! HEATING NODE NUMBER 931

C      NNQ(2)=330         ! HEATING NODE NUMBER 330

C      IALUFA(1,1)=1      ! HEAT TRANSFER BOUNDARY NODE NUMBER 1

C      IALUFA(1,2)=2      ! HEAT TRANSFER BOUNDARY NODE NUMBER 2

C      ALUFA(1)=14.E-6    ! HEAT TRANSFER COEFFICIENT

      DO 50 I=2,30

C      IALUFA(I,1)=IALUFA(I-1,1)+1

C      IALUFA(I,2)=IALUFA(I-1,2)+1

C      ALUFA(I)=14.E-6

   50 CONTINUE

C      IALUFA(31,1)=31

C      IALUFA(31,2)=62

C      ALUFA(31)=14.E-6

      DO 60 I=32,60

C      ALUFA(I)=14.E-6

      DO 62 J=1,2

C      IALUFA(I,J)=IALUFA(I-1,J)+31

   62 CONTINUE

   60 CONTINUE

      CKKAP(1)=0.00125

      CKKAP(2)=0.0125

      CKKAP(3)=0.0625

      CKKAP(4)=0.125

      CKKAP(5)=0.25

      CKKAP(6)=0.25

      CKKAP(7)=0.25

      CKKAP(8)=0.25

      CKKAP(9)=0.25

      CKKAP(10)=0.25

      TROOM=0.0        ! ambient temperature

      IDEN=0           ! heat transfer boundary   total number

      NQ=0             ! heating node total number

      NBOUND=6         ! fix temperature nude total number

      NELEM=400        ! total element number

      RAMDA=48.0       ! heat conductivity(W/mK)

      C=460.54         ! specific heat(J/kgK)

      ROH=7870.0       ! density(Kg/m^3)

      H=0.001          ! thickness (Z coordinate)(m)

      SL=0.1           ! Length(m)

      AKKAP=RAMDA/C/ROH/SL/SL

      TEND=CKKAP(10)/AKKAP+0.5   ! time end

      TINT=100.        ! initial temperature

      DT=0.1           ! dt

      IPRINT=20        ! print out interval of iteration

      AJUDG=1.E-7      ! 

      NB(1)=1          ! fix temperature node number    1

      NB(2)=101        ! fix temperature node number  101

      NB(3)=102        ! fix temperature node number  102

      NB(4)=202        ! fix temperature node number  202

      NB(5)=203        ! fix temperature node number  203

      NB(6)=303        ! fix temperature node number  303

      DO 310 I=1,10

      IPTIME(I)=CKKAP(I)/AKKAP/DT

  310 CONTINUE

      DO 220 I=1,31

C      NB(I)=I

  220 CONTINUE

      DO 230 I=32,64

C      NB(I)=I+930-31

  230 CONTINUE

      IF(NQ.NE.0)THEN

      DO 200 I=1,NQ

      DO 210 J=1,10    ! always J=1,10

C                                     ! heating value line graph

C      QTIM(NNQ(I),J)=FLOAT(J-1)*10.0 ! heating time   time number J

C      QT(NNQ(I),J)=1.0               ! heating value  time number J

  210 CONTINUE

  200 CONTINUE

      ELSE

      ENDIF

      MAT_VARIETY=1

      MAT_NO=1

      WRITE(2,1000)NODE,NELEM,NBOUND,IDEN,TROOM,NQ,DT,IPRINT,

     &TEND,TINT,MAT_VARIETY

 1000 FORMAT(4(I4,','),E16.8,',',I4,',',E16.8,',',I10,',',2(E11.4,',')

     &,I3)

      WRITE(2,300)'NODE'

  300 FORMAT(A4)

      DO 70 I=1,NODE

      WRITE(2,1100)I,X(I),Y(I)

 1100 FORMAT(I4,',',E16.8,',',E16.8)

   70 CONTINUE

      WRITE(2,300)'ELEM'

      DO 80 I=1,NELEM

      WRITE(2,1200)I,(NELM(I,J),J=1,3),MAT_NO,H

 1200 FORMAT(5(I4,','),E16.8)

   80 CONTINUE

      write(2,*)'MATR'

      write(2,1240)1,RAMDA,C,ROH

 1240 format(i3,3(',',e16.8))

      IF(NBOUND.NE.0)THEN

      WRITE(2,300)'BOUN'

      DO 1220 I=1,NBOUND

      WRITE(2,1230)NB(I),0.

 1230 FORMAT(I4,',',E16.8)

 1220 CONTINUE

      ELSE

      ENDIF

      IF(IDEN.NE.0)THEN

      WRITE(2,300)'HTRN'

      DO 90 I=1,IDEN

      WRITE(2,1300)(IALUFA(I,J),J=1,2),ALUFA(I)

 1300 FORMAT(2(I4,','),E16.8)

   90 CONTINUE

      ELSE

      ENDIF

      IF(NQ.NE.0)THEN

      WRITE(2,300)'HEAT'

      DO 2000 I=1,NQ

      WRITE(2,1400)NNQ(I),(QTIM(NNQ(I),J),QT(NNQ(I),J),J=1,10)

 1400 FORMAT(I4,10(',',E16.8,',',E16.8))

 2000 CONTINUE

      ELSE

      ENDIF

      WRITE(2,*)'POST'

      WRITE(2,1500)5   ! out put total node number (Max 40)

 1500 FORMAT(I4)

      WRITE(2,1600)11,21,31,41,51 !out put node number (time series)

 1600 FORMAT(I4,4(',',I4))

C      DO 1800 I=1,10

C      IPTIME(I)=I*11+1

C 1800 CONTINUE

      WRITE(2,*)'CONT'

      WRITE(2,1700)(IPTIME(I),I=1,10)

 1700 FORMAT(9(I15,','),I15)

      WRITE(2,*)'FEM HEAT CONDUCTION VERIFICATION'

      CLOSE(2)

      STOP

      END

【3－1．Excel-VBA-mesh】

Sub draw()

'FEM FIG DRAW

      Dim x(2500) As Single

      Dim y(2500) As Single

      Dim temp(2500) As Single

      Dim i As Integer

      Dim j As Integer

      Dim k As Integer

      Dim n As Integer

      Dim m As Integer

      Dim x1 As Single

      Dim x2 As Single

      Dim y1 As Single

      Dim y2 As Single

      Dim elm(5000, 3) As Integer

      Dim node As Integer

      Dim Nelm As Integer

      Sheets("sheet3").Select

      node = Worksheets("sheet3").Range("A2")

      Nelm = Worksheets("sheet3").Range("B2")

      [A3].Select

      For i = 1 To node

      j = ActiveCell.Offset(i, 0).Value

      x(j) = ActiveCell.Offset(i, 1).Value

      y(j) = ActiveCell.Offset(i, 2).Value

      Next i

      For i = 1 To Nelm

      j = ActiveCell.Offset(i + 1 + node, 0).Value

      elm(j, 1) = ActiveCell.Offset(i + 1 + node, 1).Value

      elm(j, 2) = ActiveCell.Offset(i + 1 + node, 2).Value

      elm(j, 3) = ActiveCell.Offset(i + 1 + node, 3).Value

      Next i

      Dim xmax As Single

      Dim xmin As Single

      Dim ymin As Single

      Dim ymax As Single

      xmax = 0#

      xmin = 999#

      ymin = 999#

      ymax = 0#

      For i = 1 To node

      If (x(i) <= xmin) Then

      xmin = x(i)

      Else

      End If

      If (x(i) >= xmax) Then

      amax = x(i)

      Else

      End If

      If (y(i) <= ymin) Then

      ymin = y(i)

      Else

      End If

      If (y(i) >= ymax) Then

      ymax = y(i)

      Else

      End If

      Next i

      Set RngStart = Worksheets("sheet2").Range("B40")

      xs = RngStart.Left

      ys = RngStart.Top

      Set RngEnd = Worksheets("sheet2").Range("D4")

      ye = RngEnd.Top

      xe = xs + (ys - ye)

      Dim zmax As Single

      Dim zmin As Single

      zmax = xmax

      If (ymax >= zmax) Then

      zmax = ymax

      Else

      End If

      zmin = xmin

      If (ymin <= zmin) Then

      zmin = ymin

      Else

      End If

      dx = (xe - xs) / (zmax - zmin)

      dy = (ye - ys) / (zmax - zmin)

      xg = (xe - xs) * (-zmin / (zmax - zmin)) + xs

      yg = (ye - ys) * (-zmin / (ymax - ymin)) + ys

      Sheets("sheet2").Select

      For n = 1 To Nelm

      i = elm(n, 1)

      j = elm(n, 2)

      k = elm(n, 3)

      x1 = x(i) * dx + xg

      y1 = y(i) * dy + yg

      x2 = x(j) * dx + xg

      y2 = y(j) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(j) * dx + xg

      y1 = y(j) * dy + yg

      x2 = x(k) * dx + xg

      y2 = y(k) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(k) * dx + xg

      y1 = y(k) * dy + yg

      x2 = x(i) * dx + xg

      y2 = y(i) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      Next n

      MsgBox ("end")

End Sub

【3－2．Excel-VBA-node-number】

Sub draw()

'FEM FIG DRAW

      Dim x(2500) As Single

      Dim y(2500) As Single

      Dim fx(2500) As Single

      Dim fy(2500) As Single

      Dim ix0(2500) As Integer

      Dim iy0(2500) As Integer

      Dim i As Integer

      Dim j As Integer

      Dim k As Integer

      Dim n As Integer

      Dim m As Integer

      Dim elm(5000, 3) As Integer

      Dim node As Integer

      Dim Nelm As Integer

      Sheets("sheet3").Select

      node = Worksheets("sheet3").Range("A2")

      Nelm = Worksheets("sheet3").Range("B2")

      [A3].Select

      For i = 1 To node

      j = ActiveCell.Offset(i, 0).Value

      x(j) = ActiveCell.Offset(i, 1).Value

      y(j) = ActiveCell.Offset(i, 2).Value

      Next i

      For i = 1 To Nelm

      j = ActiveCell.Offset(i + 1 + node, 0).Value

      elm(j, 1) = ActiveCell.Offset(i + 1 + node, 1).Value

      elm(j, 2) = ActiveCell.Offset(i + 1 + node, 2).Value

      elm(j, 3) = ActiveCell.Offset(i + 1 + node, 3).Value

      Next i

      Dim xmax As Single

      Dim xmin As Single

      Dim ymin As Single

      Dim ymax As Single

      xmax = 0#

      xmin = 999#

      ymin = 999#

      ymax = 0#

      For i = 1 To node

      If (x(i) <= xmin) Then

      xmin = x(i)

      Else

      End If

      If (x(i) >= xmax) Then

      amax = x(i)

      Else

      End If

      If (y(i) <= ymin) Then

      ymin = y(i)

      Else

      End If

      If (y(i) >= ymax) Then

      ymax = y(i)

      Else

      End If

      Next i

      Set RngStart = Worksheets("sheet2").Range("B40")

      xs = RngStart.Left

      ys = RngStart.Top

      Set RngEnd = Worksheets("sheet2").Range("D4")

      ye = RngEnd.Top

      xe = xs + (ys - ye)

      Dim zmax As Single

      Dim zmin As Single

      zmax = xmax

      If (ymax >= zmax) Then

      zmax = ymax

      Else

      End If

      zmin = xmin

      If (ymin <= zmin) Then

      zmin = ymin

      Else

      End If

      dx = (xe - xs) / (zmax - zmin)

      dy = (ye - ys) / (zmax - zmin)

      xg = (xe - xs) * (-zmin / (zmax - zmin)) + xs

      yg = (ye - ys) * (-zmin / (ymax - ymin)) + ys

      Sheets("sheet2").Select

      For n = 1 To Nelm

      i = elm(n, 1)

      j = elm(n, 2)

      k = elm(n, 3)

      x1 = x(i) * dx + xg

      y1 = y(i) * dy + yg

      x2 = x(j) * dx + xg

      y2 = y(j) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(j) * dx + xg

      y1 = y(j) * dy + yg

      x2 = x(k) * dx + xg

      y2 = y(k) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(k) * dx + xg

      y1 = y(k) * dy + yg

      x2 = x(i) * dx + xg

      y2 = y(i) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      Next n

      'node number draw

      Dim add_str As String

      Sheets("sheet2").Select

      [G40].Select

      c = ActiveSheet.Range("G40").Select

      Dim tmp_str As String

      tmp_str = Selection.Address

      Dim it As Range

      For n = 1 To node

      xa = x(n) * dx + xg

      ya = y(n) * dy + yg

      ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=xa, _

             Top:=ya, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = n

       End With

       Next n

       'element No. draw

 '      For n = 1 To Nelm

 '      i = elm(n, 1)

 '      j = elm(n, 2)

 '      k = elm(n, 3)

 '      xa = ((x(i) + x(j) + x(k)) / 3#) * dx + xg

 '      ya = ((y(i) + y(j) + y(k)) / 3#) * dy + yg

       'element No. WRITE

 '      ActiveSheet.Shapes.AddTextbox( _

 '            Orientation:=msoTextOrientationHorizontal, _

 '            Left:=xa, _

 '            Top:=ya, _

 '            Width:=40#, _

 '            Height:=20#).Select

 '            Selection.Border.LineStyle = xlLineStyleNone

 '      With Selection

 '            .Characters.Text = n

 '      End With

 '      Next n

      MsgBox ("end")

End Sub

【3－3．Excel-VBA-element-number】

Sub draw()

'FEM FIG DRAW

      Dim x(2500) As Single

      Dim y(2500) As Single

      Dim fx(2500) As Single

      Dim fy(2500) As Single

      Dim ix0(2500) As Integer

      Dim iy0(2500) As Integer

      Dim i As Integer

      Dim j As Integer

      Dim k As Integer

      Dim n As Integer

      Dim m As Integer

      Dim elm(5000, 3) As Integer

      Dim node As Integer

      Dim Nelm As Integer

      Sheets("sheet3").Select

      node = Worksheets("sheet3").Range("A2")

      Nelm = Worksheets("sheet3").Range("B2")

      [A3].Select

      For i = 1 To node

      j = ActiveCell.Offset(i, 0).Value

      x(j) = ActiveCell.Offset(i, 1).Value

      y(j) = ActiveCell.Offset(i, 2).Value

      Next i

      For i = 1 To Nelm

      j = ActiveCell.Offset(i + 1 + node, 0).Value

      elm(j, 1) = ActiveCell.Offset(i + 1 + node, 1).Value

      elm(j, 2) = ActiveCell.Offset(i + 1 + node, 2).Value

      elm(j, 3) = ActiveCell.Offset(i + 1 + node, 3).Value

      Next i

      Dim xmax As Single

      Dim xmin As Single

      Dim ymin As Single

      Dim ymax As Single

      xmax = 0#

      xmin = 999#

      ymin = 999#

      ymax = 0#

      For i = 1 To node

      If (x(i) <= xmin) Then

      xmin = x(i)

      Else

      End If

      If (x(i) >= xmax) Then

      amax = x(i)

      Else

      End If

      If (y(i) <= ymin) Then

      ymin = y(i)

      Else

      End If

      If (y(i) >= ymax) Then

      ymax = y(i)

      Else

      End If

      Next i

      Set RngStart = Worksheets("sheet2").Range("B40")

      xs = RngStart.Left

      ys = RngStart.Top

      Set RngEnd = Worksheets("sheet2").Range("D4")

      ye = RngEnd.Top

      xe = xs + (ys - ye)

      Dim zmax As Single

      Dim zmin As Single

      zmax = xmax

      If (ymax >= zmax) Then

      zmax = ymax

      Else

      End If

      zmin = xmin

      If (ymin <= zmin) Then

      zmin = ymin

      Else

      End If

      dx = (xe - xs) / (zmax - zmin)

      dy = (ye - ys) / (zmax - zmin)

      xg = (xe - xs) * (-zmin / (zmax - zmin)) + xs

      yg = (ye - ys) * (-zmin / (ymax - ymin)) + ys

      Sheets("sheet2").Select

      For n = 1 To Nelm

      i = elm(n, 1)

      j = elm(n, 2)

      k = elm(n, 3)

      x1 = x(i) * dx + xg

      y1 = y(i) * dy + yg

      x2 = x(j) * dx + xg

      y2 = y(j) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(j) * dx + xg

      y1 = y(j) * dy + yg

      x2 = x(k) * dx + xg

      y2 = y(k) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(k) * dx + xg

      y1 = y(k) * dy + yg

      x2 = x(i) * dx + xg

      y2 = y(i) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      Next n

      'node number draw

      Dim add_str As String

      Sheets("sheet2").Select

      [G40].Select

      c = ActiveSheet.Range("G40").Select

      Dim tmp_str As String

      tmp_str = Selection.Address

      Dim it As Range

 '     For n = 1 To node

 '

 '     xa = x(n) * dx + xg

 '     ya = y(n) * dy + yg

 '     ActiveSheet.Shapes.AddTextbox( _

 '            Orientation:=msoTextOrientationHorizontal, _

 '            Left:=xa, _

 '            Top:=ya, _

 '            Width:=40#, _

 '            Height:=20#).Select

 '            Selection.Border.LineStyle = xlLineStyleNone

 '      With Selection

 '            .Characters.Text = n

 '      End With

 '      Next n

       'element No. draw

       For n = 1 To Nelm

       i = elm(n, 1)

       j = elm(n, 2)

       k = elm(n, 3)

       xa = ((x(i) + x(j) + x(k)) / 3#) * dx + xg

       ya = ((y(i) + y(j) + y(k)) / 3#) * dy + yg

       'element No. WRITE

       ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=xa, _

             Top:=ya, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = n

       End With

       Next n      

      MsgBox ("end")

End Sub

【3－4．Excel-VBA-number】

Sub draw()

'FEM FIG DRAW

      Dim x(2500) As Single

      Dim y(2500) As Single

      Dim fx(2500) As Single

      Dim fy(2500) As Single

      Dim ix0(2500) As Integer

      Dim iy0(2500) As Integer

      Dim i As Integer

      Dim j As Integer

      Dim k As Integer

      Dim n As Integer

      Dim m As Integer

      Dim elm(5000, 3) As Integer

      Dim node As Integer

      Dim Nelm As Integer

      Sheets("sheet3").Select

      node = Worksheets("sheet3").Range("A2")

      Nelm = Worksheets("sheet3").Range("B2")

      [A3].Select

      For i = 1 To node

      j = ActiveCell.Offset(i, 0).Value

      x(j) = ActiveCell.Offset(i, 1).Value

      y(j) = ActiveCell.Offset(i, 2).Value

      Next i

      For i = 1 To Nelm

      j = ActiveCell.Offset(i + 1 + node, 0).Value

      elm(j, 1) = ActiveCell.Offset(i + 1 + node, 1).Value

      elm(j, 2) = ActiveCell.Offset(i + 1 + node, 2).Value

      elm(j, 3) = ActiveCell.Offset(i + 1 + node, 3).Value

      Next i

      Dim xmax As Single

      Dim xmin As Single

      Dim ymin As Single

      Dim ymax As Single

      xmax = 0#

      xmin = 999#

      ymin = 999#

      ymax = 0#

      For i = 1 To node

      If (x(i) <= xmin) Then

      xmin = x(i)

      Else

      End If

      If (x(i) >= xmax) Then

      amax = x(i)

      Else

      End If

      If (y(i) <= ymin) Then

      ymin = y(i)

      Else

      End If

      If (y(i) >= ymax) Then

      ymax = y(i)

      Else

      End If

      Next i

      Set RngStart = Worksheets("sheet2").Range("B40")

      xs = RngStart.Left

      ys = RngStart.Top

      Set RngEnd = Worksheets("sheet2").Range("D4")

      ye = RngEnd.Top

      xe = xs + (ys - ye)

      Dim zmax As Single

      Dim zmin As Single

      zmax = xmax

      If (ymax >= zmax) Then

      zmax = ymax

      Else

      End If

      zmin = xmin

      If (ymin <= zmin) Then

      zmin = ymin

      Else

      End If

      dx = (xe - xs) / (zmax - zmin)

      dy = (ye - ys) / (zmax - zmin)

      xg = (xe - xs) * (-zmin / (zmax - zmin)) + xs

      yg = (ye - ys) * (-zmin / (ymax - ymin)) + ys

      Sheets("sheet2").Select

      For n = 1 To Nelm

      i = elm(n, 1)

      j = elm(n, 2)

      k = elm(n, 3)

      x1 = x(i) * dx + xg

      y1 = y(i) * dy + yg

      x2 = x(j) * dx + xg

      y2 = y(j) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(j) * dx + xg

      y1 = y(j) * dy + yg

      x2 = x(k) * dx + xg

      y2 = y(k) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      x1 = x(k) * dx + xg

      y1 = y(k) * dy + yg

      x2 = x(i) * dx + xg

      y2 = y(i) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

      Next n

      'node number draw

      Dim add_str As String

      Sheets("sheet2").Select

      [G40].Select

      c = ActiveSheet.Range("G40").Select

      Dim tmp_str As String

      tmp_str = Selection.Address

      Dim it As Range

      For n = 1 To node

      xa = x(n) * dx + xg

      ya = y(n) * dy + yg

      ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=xa, _

             Top:=ya, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = n

       End With

       Next n

       'element No. draw

       For n = 1 To Nelm

       i = elm(n, 1)

       j = elm(n, 2)

       k = elm(n, 3)

       xa = ((x(i) + x(j) + x(k)) / 3#) * dx + xg

       ya = ((y(i) + y(j) + y(k)) / 3#) * dy + yg

       'element No. WRITE

       ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=xa, _

             Top:=ya, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = n

       End With

       Next n      

      MsgBox ("end")

End Sub

animationは、sheet11を作りそこへ、[A2]セルへ、post-ファイルをコピペし、

カンマで区切られたデータとして展開する。

実行でsheet1からsheet10にコンタ図を描く。

そのsheet1からsheet10の図をprint screenでコピーし

PowerPointに1枚ずつの図として貼り付け、

スライドショーの設定で1枚2秒表示で自動で

sheet1からsheet10を表示する。

アニメーションとなる。

【3－5．Excel-VBA-animation】

Sub draw()

'    FEM FIG DRAW

'    Color conta fig

'       written by ARAI JUN

'

      Dim x(2500) As Double

      Dim y(2500) As Double

      Dim temp(2500) As Double

      Dim tempall(2500, 10) As Double

      Dim i As Integer

      Dim j As Integer

      Dim k As Integer

      Dim n As Integer

      Dim m As Integer

      Dim i1 As Integer

      Dim j1 As Integer

      Dim k1 As Integer

      Dim x1 As Double

      Dim x2 As Double

      Dim y1 As Double

      Dim y2 As Double

      Dim elm(5000, 3) As Integer

      Dim node As Integer

      Dim Nelm As Integer

      Dim dt As Double

      Dim iptime(10) As Long

      Sheets("sheet11").Select

      node = Worksheets("sheet11").Range("A2")

      Nelm = Worksheets("sheet11").Range("B2")

      dt = Worksheets("sheet11").Range("C2")

      [A1].Select

      For i = 1 To 10

      iptime(i) = ActiveCell.Offset(1, 2 + i).Value

      Next i

      [A3].Select

      For i = 1 To node

      j = ActiveCell.Offset(i, 0).Value

      x(j) = ActiveCell.Offset(i, 1).Value

      y(j) = ActiveCell.Offset(i, 2).Value

      Next i

      For i = 1 To Nelm

      j = ActiveCell.Offset(i + 1 + node, 0).Value

      elm(j, 1) = ActiveCell.Offset(i + 1 + node, 1).Value

      elm(j, 2) = ActiveCell.Offset(i + 1 + node, 2).Value

      elm(j, 3) = ActiveCell.Offset(i + 1 + node, 3).Value

      Next i

      For i = 1 To node

      j = ActiveCell.Offset(i + 2 + Nelm + node, 0).Value

      For k = 1 To 10

      tempall(j, k) = ActiveCell.Offset(i + 2 + Nelm + node, k).Value

      Next k

      Next i

       tallmin = 9999#

       tallmax = -9999#

       n1 = 26

       For j = 1 To 10

       For i = 1 To node

       If (tempall(i, j) >= tallmax) Then

       tallmax = tempall(i, j)

       Else

       End If

       If (tempall(i, j) <= tallmin) Then

       tallmin = tempall(i, j)

       Else

       End If

       Next i

       Next j

      Dim xmax As Double

      Dim xmin As Double

      Dim ymin As Double

      Dim ymax As Double

      xmax = 0#

      xmin = 999#

      ymin = 999#

      ymax = 0#

      For i = 1 To node

      If (x(i) <= xmin) Then

      xmin = x(i)

      Else

      End If

      If (x(i) >= xmax) Then

      xmax = x(i)

      Else

      End If

      If (y(i) <= ymin) Then

      ymin = y(i)

      Else

      End If

      If (y(i) >= ymax) Then

      ymax = y(i)

      Else

      End If

      Next i

      Set RngStart = Worksheets("sheet1").Range("B40")

      xs = RngStart.Left

      ys = RngStart.Top

      Set RngEnd = Worksheets("sheet1").Range("D4")

      ye = RngEnd.Top

      xe = xs + (ys - ye)

      Dim zmax As Double

      Dim zmin As Double

      zmax = xmax

      If (ymax >= zmax) Then

      zmax = ymax

      Else

      End If

      zmin = xmin

      If (ymin <= zmin) Then

      zmin = ymin

      Else

      End If

      dx = (xe - xs) / (zmax - zmin)

      dy = (ye - ys) / (zmax - zmin)

      xg = (xe - xs) * (-zmin / (zmax - zmin)) + xs

      yg = (ye - ys) * (-zmin / (ymax - ymin)) + ys

      For jj = 1 To 10

      If (jj = 1) Then

      Sheets("sheet1").Select

      Else

      End If

      If (jj = 2) Then

      Sheets("sheet2").Select

      Else

      End If

      If (jj = 3) Then

      Sheets("sheet3").Select

      Else

      End If

      If (jj = 4) Then

      Sheets("sheet4").Select

      Else

      End If

      If (jj = 5) Then

      Sheets("sheet5").Select

      Else

      End If

      If (jj = 6) Then

      Sheets("sheet6").Select

      Else

      End If

      If (jj = 7) Then

      Sheets("sheet7").Select

      Else

      End If

      If (jj = 8) Then

      Sheets("sheet8").Select

      Else

      End If

      If (jj = 9) Then

      Sheets("sheet9").Select

      Else

      End If

      If (jj = 10) Then

      Sheets("sheet10").Select

      Else

      End If

      'outline draw  start

      For n = 1 To Nelm

      icount = 1

      i = elm(n, 1)

      j = elm(n, 2)

      k = elm(n, 3)

      'i-j line draw

      For m = 1 To Nelm

      If (n = m) Then GoTo label4

      i1 = elm(m, 1)

      j1 = elm(m, 2)

      k1 = elm(m, 3)

      If (i = i1 And j = k1) Then GoTo label1

      If (i = j1 And j = i1) Then GoTo label1

      If (i = k1 And j = j1) Then GoTo label1

      GoTo label4

label1:

      icount = icount + 1

label4:

      Next m

      If (icount >= 2) Then GoTo label6

      x1 = x(i) * dx + xg

      y1 = y(i) * dy + yg

      x2 = x(j) * dx + xg

      y2 = y(j) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

label6:

      'j-k line draw

      icount = 1

      For m = 1 To Nelm

      If (n = m) Then GoTo label9

      i1 = elm(m, 1)

      j1 = elm(m, 2)

      k1 = elm(m, 3)

      If (j = i1 And k = k1) Then GoTo label2

      If (j = j1 And k = i1) Then GoTo label2

      If (j = k1 And k = j1) Then GoTo label2

      GoTo label9

label2:

      icount = icount + 1

label9:

      Next m

      If (icount >= 2) Then GoTo label10

      x1 = x(j) * dx + xg

      y1 = y(j) * dy + yg

      x2 = x(k) * dx + xg

      y2 = y(k) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

label10:

      'k-i line draw

      icount = 1

      For m = 1 To Nelm

      If (n = m) Then GoTo label12

      i1 = elm(m, 1)

      j1 = elm(m, 2)

      k1 = elm(m, 3)

      If (k = i1 And i = k1) Then GoTo label3

      If (k = j1 And i = i1) Then GoTo label3

      If (k = k1 And i = j1) Then GoTo label3

      GoTo label12

label3:

      icount = icount + 1

label12:

      Next m

      If (icount >= 2) Then GoTo label13

      x1 = x(k) * dx + xg

      y1 = y(k) * dy + yg

      x2 = x(i) * dx + xg

      y2 = y(i) * dy + yg

      ActiveSheet.Shapes.AddLine x1, y1, x2, y2

label13:

      Next n

      'outline draw end

       'color-conta draw start

       Dim myCht As Chart

       Dim mySts As Series

       Dim Npts As Integer

       Dim myBuilder As FreeformBuilder

       Dim myShape As Shape

       For i = 1 To node

       temp(i) = tempall(i, jj)

       Next i

       tmin = 9999#

       tmax = -9999#

       n1 = 26

       For i = 1 To node

       If (temp(i) >= tmax) Then

       tmax = temp(i)

       Else

       End If

       If (temp(i) <= tmin) Then

       tmin = temp(i)

       Else

       End If

       Next i

       For n = 1 To Nelm

       i = elm(n, 1)

       j = elm(n, 2)

       k = elm(n, 3)

       x1 = x(i)

       y1 = y(i)

       t1 = temp(i)

       x2 = x(j)

       y2 = y(j)

       t2 = temp(j)

       x3 = x(k)

       y3 = y(k)

       t3 = temp(k)

       Call sub_color_conta(x1, y1, t1, x2, y2, t2, x3, y3, t3, xg, yg, dx, dy, tallmin, tallmax, n1)

       Next n

           ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=0, _

             Top:=20, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = "TIME="

'             .Size = 0.2

       End With

      ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=80, _

             Top:=20, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = dt * iptime(jj)

'             .Size = 0.2

       End With

     ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=0, _

             Top:=70, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = "MAX"

'             .Size = 0.2

       End With

      ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=0, _

             Top:=90, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = tmax

'             .Size = 0.2

       End With

      ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=0, _

             Top:=170, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = "MIN"

'             .Size = 0.2

       End With

      ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=0, _

             Top:=190, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = tmin

'             .Size = 0.2

       End With

       If (tmin < 0#) Then

         ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=10, _

             Top:=190, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

         With Selection

             .Characters.Text = Abs(tmin)

'             .Size = 0.2

         End With

       Else

       End If

       Dim temp2(100)

       dtemp = (tallmax - tallmin) / (n1 - 1)

       For i = 1 To n1

       temp2(i) = tallmin + (i - 1) * dtemp

       Next i

       For i = 1 To n1 - 1

       x1 = xmax * dx + xg + 100

       y1 = ys - 24 * i

       x2 = xmax * dx + xg + 150

       y2 = y1

       x3 = x2

       y3 = ys - 24 * (i - 1)

       x4 = x1

       y4 = y3

       Call poly4_scale(x1, y1, x2, y2, x3, y3, x4, y4, i)

       Next i

       For i = 1 To n1

       x1 = xmax * dx + xg + 160

       y1 = ys - 24 * (i - 1) - 10

            ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=x1, _

             Top:=y1, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

       With Selection

             .Characters.Text = temp2(i)

'             .Size = 0.2

       End With

       If (temp2(i) < 0#) Then

        ActiveSheet.Shapes.AddTextbox( _

             Orientation:=msoTextOrientationHorizontal, _

             Left:=x1 + 10#, _

             Top:=y1, _

             Width:=40#, _

             Height:=20#).Select

             Selection.Border.LineStyle = xlLineStyleNone

        With Selection

             .Characters.Text = Abs(temp2(i))

'             .Size = 0.2

        End With

       Else

       End If

       Next i

       Next jj

      'color-conta draw  end      

      MsgBox ("end")

      End Sub

      Sub sub_color_conta(x1, y1, t1, x2, y2, t2, x3, y3, t3, xg, yg, dx, dy, tallmin, tallmax, n1)

      Dim tt(3) As Double

      Dim xx(3) As Double

      Dim yy(3) As Double

      Dim temp1(101) As Double

      Dim xa As Double

      Dim ya As Double

      Dim xb As Double

      Dim yb As Double

      Dim xa1 As Double

      Dim ya1 As Double

      Dim xb1 As Double

      Dim yb1 As Double

      Dim tq As Double

      Dim xq As Double

      Dim yq As Double

      Dim k1 As Integer

      Dim k2 As Integer

      Dim k3 As Integer

      Dim ta As Double

      Dim tb As Double

      Dim dtemp As Double

      dtemp = (tallmax - tallmin) / (n1 - 1)

      For i = 1 To n1

      temp1(i) = tallmin + (i - 1) * dtemp

      Next i

      xx(1) = x1

      yy(1) = y1

      tt(1) = t1

      xx(2) = x2

      yy(2) = y2

      tt(2) = t2

      xx(3) = x3

      yy(3) = y3

      tt(3) = t3

      'quick-sort

      For i = 1 To 3

      For j = 1 To 2

       If (tt(j) >= tt(j + 1)) Then

        tq = tt(j)

        tt(j) = tt(j + 1)

        tt(j + 1) = tq

        xq = xx(j)

        xx(j) = xx(j + 1)

        xx(j + 1) = xq

        yq = yy(j)

        yy(j) = yy(j + 1)

        yy(j + 1) = yq

       Else

       End If

      Next j

      Next i

      x1 = xx(1)

      x2 = xx(2)

      x3 = xx(3)

      y1 = yy(1)

      y2 = yy(2)

      y3 = yy(3)

      t1 = tt(1)

      t2 = tt(2)

      t3 = tt(3)

      j = 0

      For i = 1 To n1 - 1

      If (t1 >= temp1(i) And t1 <= temp1(i + 1)) Then

         If (j = 0) Then

          k1 = i

          j = j + 1

          Else

          End If

         Else

         End If

      Next i

      For i = k1 To n1 - 1

      If (t3 >= temp1(i) And t3 <= temp1(i + 1)) Then

      k3 = i

      Else

      End If

      Next i

      If (k1 = k3) Then

      k2 = k1

      Else

       For i = k1 To k3

        If (t2 >= temp1(i) And t2 <= temp1(i + 1)) Then

         k2 = i

        Else

        End If

       Next i

      End If

      k4 = k3 - k2

      k5 = k2 - k1

      If (k1 = k3) Then GoTo label200

      If (k1 = k2) Then

          'side 1-3  -  side 2-3   line draw  k1 = k2 <> k3

            For i = k2 + 1 To k3

            If (t3 = t1) Then

            ta = 0#

            Else

            ta = (temp1(i) - t1) / (t3 - t1)

            End If

            If (t2 = t3) Then

            tb = 0#

            Else

            tb = (temp1(i) - t2) / (t3 - t2)

            End If

            If (i = 1) Then

            icol = 2

            Else

            icol = i

            End If

            xa = x1 + (x3 - x1) * ta

            ya = y1 + (y3 - y1) * ta

            xb = x2 + (x3 - x2) * tb

            yb = y2 + (y3 - y2) * tb

            If (i = k2 + 1) Then Call poly4(x1, y1, x2, y2, xa, ya, xb, yb, dx, dy, xg, yg, icol - 1)

            If (i > k2 + 1 And i <= k3) Then Call poly4(xa, ya, xb, yb, xa1, ya1, xb1, yb1, dx, dy, xg, yg, icol - 1)

            If (i = k3) Then Call poly3(xa, ya, xb, yb, x3, y3, dx, dy, xg, yg, icol)

            xa1 = xa

            ya1 = ya

            xb1 = xb

            yb1 = yb

            Next i

          Else

          If (k2 = k3) Then

              'side 1-2  - side 1-3     line draw    K1 <> K2 = K3

              For i = k1 + 1 To k2

              If (t3 = t1) Then

              ta = 1#

              Else

              ta = (temp1(i) - t1) / (t3 - t1)

              End If

              If (t1 = t2) Then

              tb = 1#

              Else

              tb = (temp1(i) - t1) / (t2 - t1)

              End If

              If (i = 1) Then

              icol = 2

              Else

              icol = i

              End If

              xa = x1 + (x3 - x1) * ta

              ya = y1 + (y3 - y1) * ta

              xb = x1 + (x2 - x1) * tb

              yb = y1 + (y2 - y1) * tb

              If (i = k1 + 1) Then Call poly3(xa, ya, xb, yb, x1, y1, dx, dy, xg, yg, icol - 1)

              If (i > k1 + 1 And i <= k2) Then Call poly4(xa, ya, xb, yb, xa1, ya1, xb1, yb1, dx, dy, xg, yg, icol - 1)

              If (i = k2) Then Call poly4(xa, ya, xb, yb, x2, y2, x3, y3, dx, dy, xg, yg, icol)

              xa1 = xa

              ya1 = ya

              xb1 = xb

              yb1 = yb

              Next i

            Else

                         'side 1-2  - side 1-3     line draw    k1 <> k2 <> k3

                 For i = k1 + 1 To k2

                 If (t3 = t1) Then

                 ta = 1#

                 Else

                 ta = (temp1(i) - t1) / (t3 - t1)

                 End If

                 If (t1 = t2) Then

                 tb = 1#

                 Else

                 tb = (temp1(i) - t1) / (t2 - t1)

                 End If

                 If (i = 1) Then

                 icol = 2

                 Else

                 icol = i

                 End If

                 xa = x1 + (x3 - x1) * ta

                 ya = y1 + (y3 - y1) * ta

                 xb = x1 + (x2 - x1) * tb

                 yb = y1 + (y2 - y1) * tb

                 If (i = k1 + 1) Then Call poly3(xa, ya, xb, yb, x1, y1, dx, dy, xg, yg, icol - 1)

                 If (i > k1 + 1 And i <= k2) Then Call poly4(xa, ya, xb, yb, xa1, ya1, xb1, yb1, dx, dy, xg, yg, icol - 1)

                 xa1 = xa

                 ya1 = ya

                 xb1 = xb

                 yb1 = yb

                 Next i

           'side 1-3  -  side 2-3   line draw   k1 <> k2 <> k3

                 For i = k2 + 1 To k3

                 If (t3 = t1) Then

                 ta = 0#

                 Else

                 ta = (temp1(i) - t1) / (t3 - t1)

                 End If

                 If (t2 = t3) Then

                 tb = 0#

                 Else

                 tb = (temp1(i) - t2) / (t3 - t2)

                 End If

                 xa = x1 + (x3 - x1) * ta

                 ya = y1 + (y3 - y1) * ta

                 xb = x2 + (x3 - x2) * tb

                 yb = y2 + (y3 - y2) * tb

                 If (i = 1) Then

                 icol = 2

                 Else

                 icol = i

                 End If

                 If (i = k2 + 1) Then Call poly5(xa1, ya1, xb1, yb1, xa, ya, xb, yb, x2, y2, dx, dy, xg, yg, icol - 1)

labelc100:

                 If (i > k2 + 1 And i <= k3) Then Call poly4(xa, ya, xb, yb, xa1, ya1, xb1, yb1, dx, dy, xg, yg, icol - 1)

                 If (i = k3) Then Call poly3(xa, ya, xb, yb, x3, y3, dx, dy, xg, yg, icol)

                 xa1 = xa

                 ya1 = ya

                 xb1 = xb

                 yb1 = yb

                 Next i

            End If

        End If

       GoTo label400

label200:

       If (k1 = 0) Then

       icol = 1

       Else

       icol = k1

       End If

       Call poly3(x1, y1, x2, y2, x3, y3, dx, dy, xg, yg, icol)

label400:

      End Sub

      Sub arctan(xr, yr, phi)

      apai = 3.14159265358979

      If (xr = 0# And yr >= 0#) Then

      phi = apai / 2#

      Else

      End If

      If (xr = 0# And yr <= 0#) Then

      phi = -apai / 2#

      Else

      End If

      If (xr = 0#) Then GoTo label1000

      a = yr / xr

      b = Atn(a)

      phi = b

      If (xr <= 0# And yr <= 0#) Then

      phi = b - apai

      Else

      End If

      If (xr <= 0# And yr >= 0#) Then

      phi = b + apai

      Else

      End If

label1000:

      End Sub

Sub poly3(x1, y1, x2, y2, x3, y3, dx, dy, xg, yg, icol)

      Dim xq3(3) As Double

      Dim yq3(3) As Double

      Dim theta3(3) As Double

      xq3(1) = x1

      xq3(2) = x2

      xq3(3) = x3

      yq3(1) = y1

      yq3(2) = y2

      yq3(3) = y3

      xq = (x1 + x2 + x3) / 3#

      yq = (y1 + y2 + y3) / 3#

      xr = x1 - xq

      yr = y1 - yq

      Call arctan(xr, yr, phi)

      theta3(1) = phi

      xr = x2 - xq

      yr = y2 - yq

      Call arctan(xr, yr, phi)

      theta3(2) = phi

      xr = x3 - xq

      yr = y3 - yq

      Call arctan(xr, yr, phi)

      theta3(3) = phi

      'quick-sort

      For i = 1 To 3

      For j = 1 To 2

      If (theta3(j) <= theta3(j + 1)) Then

      ta = theta3(j)

      theta3(j) = theta3(j + 1)

      theta3(j + 1) = ta

      xq = xq3(j)

      xq3(j) = xq3(j + 1)

      xq3(j + 1) = xq

      yq = yq3(j)

      yq3(j) = yq3(j + 1)

      yq3(j + 1) = yq

      Else

      End If

      Next j

      Next i

       Xnode = xq3(1) * dx + xg

       Ynode = yq3(1) * dy + yg

       Set myBuilder = ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, Xnode, Ynode)

       For i = 2 To 3

       Xnode = xq3(i) * dx + xg

       Ynode = yq3(i) * dy + yg

       myBuilder.AddNodes msoSegmentLine, msoEditingAuto, Xnode, Ynode

       Next i

       Set myShape = myBuilder.ConvertToShape

       With myShape

       .Fill.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.Weight = 1.5

       End With

       End Sub

      Sub poly4(x1, y1, x2, y2, x3, y3, x4, y4, dx, dy, xg, yg, icol)

      Dim xq4(4) As Double

      Dim yq4(4) As Double

      Dim theta4(4) As Double

      xq4(1) = x1

      xq4(2) = x2

      xq4(3) = x3

      yq4(1) = y1

      yq4(2) = y2

      yq4(3) = y3

      xq4(4) = x4

      yq4(4) = y4

      xq = (x1 + x2 + x3 + x4) / 4#

      yq = (y1 + y2 + y3 + y4) / 4#

      xr = x1 - xq

      yr = y1 - yq

      Call arctan(xr, yr, phi)

      theta4(1) = phi

      xr = x2 - xq

      yr = y2 - yq

      Call arctan(xr, yr, phi)

      theta4(2) = phi

      xr = x3 - xq

      yr = y3 - yq

      Call arctan(xr, yr, phi)

      theta4(3) = phi

      xr = x4 - xq

      yr = y4 - yq

      Call arctan(xr, yr, phi)

      theta4(4) = phi

      'quick-sort

      For i = 1 To 4

      For j = 1 To 3

      If (theta4(j) <= theta4(j + 1)) Then

      ta = theta4(j)

      theta4(j) = theta4(j + 1)

      theta4(j + 1) = ta

      xq = xq4(j)

      xq4(j) = xq4(j + 1)

      xq4(j + 1) = xq

      yq = yq4(j)

      yq4(j) = yq4(j + 1)

      yq4(j + 1) = yq

      Else

      End If

      Next j

      Next i

       Xnode = xq4(1) * dx + xg

       Ynode = yq4(1) * dy + yg

       Set myBuilder = ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, Xnode, Ynode)

       For i = 2 To 4

       Xnode = xq4(i) * dx + xg

       Ynode = yq4(i) * dy + yg

       myBuilder.AddNodes msoSegmentLine, msoEditingAuto, Xnode, Ynode

       Next i

       Set myShape = myBuilder.ConvertToShape

       With myShape

       .Fill.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.Weight = 1.5

       End With

      End Sub

Sub poly5(x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, dx, dy, xg, yg, icol)

      Dim xq5(5) As Double

      Dim yq5(5) As Double

      Dim theta5(5) As Double

      xq5(1) = x1

      xq5(2) = x2

      xq5(3) = x3

      yq5(1) = y1

      yq5(2) = y2

      yq5(3) = y3

      xq5(4) = x4

      yq5(4) = y4

      xq5(5) = x5

      yq5(5) = y5

      xq = (x1 + x2 + x3 + x4 + x5) / 5#

      yq = (y1 + y2 + y3 + y4 + y5) / 5#

      xr = x1 - xq

      yr = y1 - yq

      Call arctan(xr, yr, phi)

      theta5(1) = phi

      xr = x2 - xq

      yr = y2 - yq

      Call arctan(xr, yr, phi)

      theta5(2) = phi

      xr = x3 - xq

      yr = y3 - yq

      Call arctan(xr, yr, phi)

      theta5(3) = phi

      xr = x4 - xq

      yr = y4 - yq

      Call arctan(xr, yr, phi)

      theta5(4) = phi

      xr = x5 - xq

      yr = y5 - yq

      Call arctan(xr, yr, phi)

      theta5(5) = phi

      'quick-sort

      For i = 1 To 5

      For j = 1 To 4

      If (theta5(j) <= theta5(j + 1)) Then

      ta = theta5(j)

      theta5(j) = theta5(j + 1)

      theta5(j + 1) = ta

      xq = xq5(j)

      xq5(j) = xq5(j + 1)

      xq5(j + 1) = xq

      yq = yq5(j)

      yq5(j) = yq5(j + 1)

      yq5(j + 1) = yq

      Else

      End If

      Next j

      Next i

       Xnode = xq5(1) * dx + xg

       Ynode = yq5(1) * dy + yg

       Set myBuilder = ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, Xnode, Ynode)

       For i = 2 To 5

       Xnode = xq5(i) * dx + xg

       Ynode = yq5(i) * dy + yg

       myBuilder.AddNodes msoSegmentLine, msoEditingAuto, Xnode, Ynode

       Next i

       Set myShape = myBuilder.ConvertToShape

       With myShape

       .Fill.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.Weight = 1.5

       End With

      End Sub

       Sub poly4_scale(x1, y1, x2, y2, x3, y3, x4, y4, icol)

       Dim xq4_c(4) As Double

       Dim yq4_c(4) As Double

       xq4_c(1) = x1

       xq4_c(2) = x2

       xq4_c(3) = x3

       xq4_c(4) = x4

       yq4_c(1) = y1

       yq4_c(2) = y2

       yq4_c(3) = y3

       yq4_c(4) = y4

       Xnode = xq4_c(1)

       Ynode = yq4_c(1)

       Set myBuilder = ActiveSheet.Shapes.BuildFreeform(msoEditingAuto, Xnode, Ynode)

       For i = 2 To 4

       Xnode = xq4_c(i)

       Ynode = yq4_c(i)

       myBuilder.AddNodes msoSegmentLine, msoEditingAuto, Xnode, Ynode

       Next i

       Set myShape = myBuilder.ConvertToShape

       With myShape

       .Fill.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.ForeColor.RGB = RGB((icol - 1) * 10, 0, (25 - icol) * 10) '

       .Line.Weight = 1.5

       End With

End Sub