!オンサイト木橋の3点曲げのモデルを作るプログラム
!局部座屈をみれるかどうか
implicit real*8(a-h,o-z)
dimension x(200000),y(200000),z(200000)&
&,narax(200000),naray(200000),naraz(200000)&
&,e66(6,6),b(3,3)&
&,meshmoku(400000,8),meshkou(400000,8)
!材料の寸法の設定
b1=840.0d-3!木材の幅
b2=9.0d-3!鋼板の幅
h1=10.0d-3!押し板までの高さ
h2=100.0d-3!押し板の高さ
h3=h1!押し板上から木材までの高さ
h4=40.0d-3
h5=180.0d-3
hall=h1*4.0d0+h2*2.0d0+h4*2.0d0+h5
ellitahan=50.0d-3!スパンの半分の長さ
ellhaji=150.0d-3!出ている部分の長さ
ellhaji2=165.0d-3
ellana=265.0d-3
ellaida=ellana*2.0d0+ellitahan*2.0d0
p=-1.000d-1!荷重
! x,y,z方向の要素分割数
nx1=3!木材の分割数
nx2=2!鋼板の分割数
nx=nx1+nx2*2
ny1=1!押し板までの分割数
ny2=1!押し板の分割数
ny3=ny1!押し板上から木材終わりまでの分割数
ny4=2!間の鋼板の分割数
ny5=4!間の鋼板の分割数
nzitahan=1!軸方向の半分の分割数
nzana=2
nzaida=nzana*2+nzitahan*2
nzhaji=1!出ている部分の分割数
nzhaji2=3!出ている部分の分割数
nitimen=(nx+1)*(ny1*2+ny2+1)*2+(nx2+1)*(ny4*2+ny5-1)*2
nitimenana=(nx+1)*(ny1*2+ny2+1)*2+(nx2+1)*ny4*2*2
nosi=10!押し板の数　プレストレスの載荷位置の数(z座標)上下なのでこれの2倍
nza=nzhaji2*2+nzhaji*2+nzitahan*nosi*2+nzaida*(nosi-1)
nzananasi=nzhaji*2+nzhaji2*2+(nzitahan*2+1)*(nosi*2-1)
nzanaari=(nzana-1)*2*(nosi-1)
nsetten=nitimen*nzananasi+nitimenana*nzanaari
ellall=ellhaji*2.0d0+ellhaji2*2.0d0+ellaida*(nosi-1)+ellitahan*2.0d0*nosi
!材料の数
!ヤング率の入力
ezz=7.694d3!z軸方向のヤング率
exx=ezz/25.0d0!x軸方向のヤング率
eyy=ezz/25.0d0!y軸方向のヤング率
!ポアソン比の入力
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!mk=2!鋼板の数
! 1要素のx,y,z方向の長さ：xn,yn,zn
xn1=b1/nx1!木材の1要素あたりの幅方向長さ
xn2=b2/nx2!鋼板の1要素あたりの幅方向長さ
yn1=h1/ny1!押し板までの1要素あたりの高さ方向の長さ 
yn2=h2/ny2!押し板の1要素あたりの高さ方向の長さ
yn3=h3/ny3!押し板上から木材終わりまでの1要素あたりの高さ方向の長さ
yn4=h4/ny4!間の鋼板の1要素あたりの高さ方向の長さ
yn5=h5/ny5!160-340
znhaji2=ellhaji2/nzhaji2!支点間の1要素あたりの長さ
znhaji=ellhaji/nzhaji!出ている部分の1要素あたりの長さ
znana=ellana/nzana
znitahan=ellitahan/nzitahan
znaida=ellaida/nzaida
!print *,nza
!print *,nzananasi+nzanaari
!節点の設定
!下と上部分の節点(軸方向の端っこ)
do l=1,2
do k=1,(nzhaji+nzhaji2)*2
do j=1,ny1*2+ny2+1
do i=1,nx+1
if (k>=1.and.k<=nzhaji+nzhaji2) then
nset=i+(j-1)*(nx+1)+(k-1)*nitimen&
&+(l-1)*(nitimen-(nx+1)*(ny1*2+ny2+1))
naraz(nset)=k
else
nset=i+(j-1)*(nx+1)+(k-nzhaji-nzhaji2)*nitimen&
&+(l-1)*(nitimen-(nx+1)*(ny1*2+ny2+1))&
&+nitimen*(nzananasi-nzhaji-nzhaji2-1)&
&+nitimenana*nzanaari
naraz(nset)=(k-nzhaji-nzhaji2)+nza+1-nzhaji-nzhaji2
end if
narax(nset)=i
naray(nset)=j+(l-1)*(ny1*2+ny2+ny4*2+ny5)
if (i>=1.and.i<=nx2+1) then
x(nset)=xn2*(i-1)
else if (i>=nx2+2.and.i<=nx2+nx1) then
x(nset)=xn1*(i-nx2-2)+xn1+b2
else
x(nset)=xn2*(i-nx2-nx1-1)+b2+b1
end if
if (j>=1.and.j<=ny1+1) then
y(nset)=yn1*(j-1)+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
else if (j>=ny1+2.and.j<=ny1+ny2+1) then
y(nset)=yn2*(j-ny1-2)+h1+yn2+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
else 
y(nset)=yn1*(j-ny1-ny2-2)+h1+h2+yn1+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
end if
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2) then
z(nset)=znhaji2*(k-nzhaji-1)+ellhaji
else if (k>=nzhaji+nzhaji2+1.and.k<=nzhaji+nzhaji2*2) then
z(nset)=znhaji2*(k-nzhaji-nzhaji2)+ellall-ellhaji2-ellhaji
else
z(nset)=znhaji*(k-nzhaji-nzhaji2*2)+ellall-ellhaji
end if
end do
end do
end do
end do
!下と上部分の節点(軸方向穴なし)
do m=1,nosi*2-1
do l=1,2
do k=1,nzitahan*2+1
do j=1,ny1*2+ny2+1
do i=1,nx+1
nset=i+(j-1)*(nx+1)+(k-1)*nitimen&
&+(l-1)*(nitimen-(nx+1)*(ny1*2+ny2+1))&
&+(m-1)*(nitimen*(nzitahan*2+1)+&
&nitimenana*(nzana-1))&
&+nitimen*(nzhaji+nzhaji2)
naraz(nset)=nzhaji+nzhaji2+k+(m-1)*(nzitahan*2+nzana)
narax(nset)=i
naray(nset)=j+(l-1)*(ny1*2+ny2+ny4*2+ny5)
if (i>=1.and.i<=nx2+1) then
x(nset)=xn2*(i-1)
else if (i>=nx2+2.and.i<=nx2+nx1) then
x(nset)=xn1*(i-nx2-2)+xn1+b2
else
x(nset)=xn2*(i-nx2-nx1-1)+b2+b1
end if
if (j>=1.and.j<=ny1+1) then
y(nset)=yn1*(j-1)+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
else if (j>=ny1+2.and.j<=ny1+ny2+1) then
y(nset)=yn2*(j-ny1-2)+h1+yn2+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
else 
y(nset)=yn1*(j-ny1-ny2-2)+h1+h2+yn1+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
end if
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+&
&(m-1)*(ellitahan*2+ellana)
end do
end do
end do
end do
end do
!下と上部分の節点(軸方向穴あり)
do m=1,(nosi-1)*2
do l=1,2
do k=1,nzana-1
do j=1,ny1*2+ny2+1
do i=1,nx+1
nset=i+(j-1)*(nx+1)+(k-1)*nitimenana&
&+(l-1)*(nitimenana-(nx+1)*(ny1*2+ny2+1))&
&+(m-1)*(nitimen*(nzitahan*2+1)+&
&nitimenana*(nzana-1))&
&+nitimen*(nzhaji+nzhaji2+nzitahan*2+1)
naraz(nset)=nzhaji+nzhaji2+nzitahan*2+1+k+(m-1)*(nzitahan*2+nzana)
narax(nset)=i
naray(nset)=j+(l-1)*(ny1*2+ny2+ny4*2+ny5)
if (i>=1.and.i<=nx2+1) then
x(nset)=xn2*(i-1)
else if (i>=nx2+2.and.i<=nx2+nx1) then
x(nset)=xn1*(i-nx2-2)+xn1+b2
else
x(nset)=xn2*(i-nx2-nx1-1)+b2+b1
end if
if (j>=1.and.j<=ny1+1) then
y(nset)=yn1*(j-1)+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
else if (j>=ny1+2.and.j<=ny1+ny2+1) then
y(nset)=yn2*(j-ny1-2)+h1+yn2+(l-1)*(h1*2.0d0+h2+h4*2.0d0+h5)
else 
y(nset)=yn1*(j-ny1-ny2-2)+h1+h2+yn1+(l-1)*(h1*2.0d0+h4*2.0d0+h2+h5)
end if
z(nset)=k*znana+ellhaji+ellhaji2+&
&ellitahan*2.0d0+&
&(m-1)*(ellitahan*2+ellana)
end do
end do
end do
end do
end do
!
!間の鋼板の節点(軸方向の端っこ)
do l=1,2
do k=1,(nzhaji+nzhaji2)*2
do j=1,ny4*2+ny5-1
do i=1,nx2+1
if (k>=1.and.k<=nzhaji+nzhaji2) then
nset=i+(j-1)*(nx2+1)*2+(k-1)*nitimen&
&+(l-1)*(nx2+1)&
&+(nx+1)*(ny1*2+ny2+1)
naraz(nset)=k
else
nset=i+(j-1)*(nx2+1)*2+(k-nzhaji-nzhaji2)*nitimen&
&+(l-1)*(nx2+1)&
&+nitimen*(nzananasi-nzhaji-nzhaji2-1)&
&+nitimenana*nzanaari&
&+(nx+1)*(ny1*2+ny2+1)
naraz(nset)=k-nzhaji-nzhaji2+(nza+1-nzhaji-nzhaji2)
end if
narax(nset)=i+(l-1)*(nx2+nx1)
naray(nset)=j+ny1*2+ny2+1
x(nset)=(i-1)*xn2+(l-1)*(b1+b2)
if (j>=1.and.j<=ny4) then
y(nset)=yn4*j&
&+h1*2.0d0+h2
else if (j>=ny4+1.and.j<=ny4+ny5) then
y(nset)=yn5*(j-ny4)&
&+h1*2.0d0+h2+h4
else 
y(nset)=yn4*(j-ny4-ny5)&
&+h1*2.0d0+h2+h4+h5
end if
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2) then
z(nset)=znhaji2*(k-nzhaji-1)+ellhaji
else if (k>=nzhaji+nzhaji2+1.and.k<=nzhaji+nzhaji2*2) then
z(nset)=znhaji2*(k-nzhaji-nzhaji2)+ellall-ellhaji2-ellhaji
else
z(nset)=znhaji*(k-nzhaji-nzhaji2*2)+ellall-ellhaji
end if
end do
end do
end do
end do
!間の鋼板の節点(軸方向穴なし)
do m=1,nosi*2-1
do l=1,2
do k=1,nzitahan*2+1
do j=1,ny4*2+ny5-1
do i=1,nx2+1
nset=i+(j-1)*(nx2+1)*2+(k-1)*nitimen&
&+(l-1)*(nx2+1)&
&+(m-1)*(nitimen*(nzitahan*2+1)+&
&nitimenana*(nzana-1))&
&+nitimen*(nzhaji+nzhaji2)&
&+(nx+1)*(ny1*2+ny2+1)
naraz(nset)=nzhaji+nzhaji2+k+(m-1)*(nzitahan*2+nzana)
narax(nset)=i+(l-1)*(nx2+nx1)
naray(nset)=j+ny1*2+ny2+1
x(nset)=(i-1)*xn2+(l-1)*(b1+b2)
if (j>=1.and.j<=ny4) then
y(nset)=yn4*j&
&+h1*2.0d0+h2
else if (j>=ny4+1.and.j<=ny4+ny5) then
y(nset)=yn5*(j-ny4)&
&+h1*2.0d0+h2+h4
else 
y(nset)=yn4*(j-ny4-ny5)&
&+h1*2.0d0+h2+h4+h5
end if
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+&
&(m-1)*(ellitahan*2+ellana)
end do
end do
end do
end do
end do
!
!print'("*NODE, NSET=Nall")'
!do i=1,nsetten
!print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
!&i,',',x(i),',',y(i),',',z(i)
!end do
!間の鋼板の節点(軸方向穴あり)
do m=1,(nosi-1)*2
do l=1,2
do k=1,nzana-1
do j=1,ny4*2
do i=1,nx2+1
nset=i+(j-1)*(nx2+1)*2+(k-1)*nitimenana&
&+(l-1)*(nx2+1)&
&+(m-1)*(nitimen*(nzitahan*2+1)+&
&nitimenana*(nzana-1))&
&+nitimen*(nzhaji+nzhaji2+nzitahan*2+1)&
&+(nx+1)*(ny1*2+ny2+1)
naraz(nset)=nzhaji+nzhaji2+nzitahan*2+1+k+(m-1)*(nzitahan*2+nzana)
narax(nset)=i+(l-1)*(nx2+nx1)
x(nset)=(i-1)*xn2+(l-1)*(b1+b2)
if (j>=1.and.j<=ny4) then
y(nset)=yn4*j&
&+h1*2.0d0+h2
naray(nset)=j+ny1*2+ny2+1
else
y(nset)=yn4*(j-ny4-1)&
&+h1*2.0d0+h2+h4+h5
naray(nset)=(j-ny4)+ny1*2+ny2+ny4+ny5
end if
z(nset)=k*znana+ellhaji+ellhaji2+&
&ellitahan*2.0d0+&
&(m-1)*(ellitahan*2+ellana)
end do
end do
end do
end do
end do
!
print'("*NODE, NSET=Nall")'
do i=1,nsetten
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
&i,',',x(i),',',y(i),',',z(i)
end do
!要素の設定
!木材部分
nyoum=0
do j=1,nx1
do i=1,nsetten
if (x(i)<=b2+0.000001d0+(j-1)*xn1.and.&
&x(i)>=b2-0.0000001d0+(j-1)*xn1&
&.and.y(i)<=h1*2.0d0+h2-0.00001d0&
&.and.z(i)<=ellall-znhaji+0.0000001d0) then
nyoum=nyoum+1
meshmoku(nyoum,1)=i
else if (x(i)<=b2+0.000001d0+(j-1)*xn1.and.&
&x(i)>=b2-0.0000001d0+(j-1)*xn1&
&.and.y(i)<=hall-0.00001d0&
&.and.y(i)>=h1*2.0d0+h2+h4*2.0d0+h5-0.000001d0&
&.and.z(i)<=ellall-znhaji+0.0000001d0) then
nyoum=nyoum+1
meshmoku(nyoum,1)=i
end if
end do
end do
!
do j=1,nsetten
do i=1,nyoum
if (naraz(j)==(naraz(meshmoku(i,1))+1)&
&.and.y(j)>=y(meshmoku(i,1))-0.000001d0&
&.and.y(j)<=y(meshmoku(i,1))+0.000001d0&
&.and.x(j)<=x(meshmoku(i,1))+0.000001d0&
&.and.x(j)>=x(meshmoku(i,1))-0.000001d0) then
meshmoku(i,4)=j
end if
end do
end do
!
!
!do i=1,nsetten
!print *,i,naraz(i)
!end do
!
do i=1,nyoum
meshmoku(i,2)=meshmoku(i,1)+nx+1
meshmoku(i,3)=meshmoku(i,4)+nx+1
meshmoku(i,5)=meshmoku(i,1)+1
meshmoku(i,6)=meshmoku(i,2)+1
meshmoku(i,7)=meshmoku(i,3)+1
meshmoku(i,8)=meshmoku(i,4)+1
end do
!
!要素の設定
!鋼板部分
nyouk=0
do m=1,2
do j=1,nx2
do i=1,nsetten
if (narax(i)==j+(m-1)*(nx1+nx2)&
&.and.naray(i)<(ny1*2+ny2+1)&
&.and.naraz(i)<nza+1) then
nyouk=nyouk+1
meshkou(nyouk,1)=i
end if
end do
end do
end do
!
do m=1,2
do j=1,nx2
do i=1,nsetten
if (narax(i)==j+(m-1)*(nx1+nx2)&
&.and.naray(i)>(ny1*2+ny2+ny4*2+ny5)&
&.and.naray(i)<(ny1*2+ny2+1)*2+ny4*2+ny5-1&
&.and.naraz(i)<nza+1) then
nyouk=nyouk+1
meshkou(nyouk,1)=i
end if
end do
end do
end do
!
!
do m=1,2
do j=1,nx2
do i=1,nsetten
if (narax(i)==j+(m-1)*(nx1+nx2)&
&.and.naray(i)<(ny1*2+ny2+ny4*2+ny5+1)&
&.and.naray(i)>(ny1*2+ny2)&
&.and.naraz(i)<nzhaji+nzhaji2+1) then
nyouk=nyouk+1
meshkou(nyouk,1)=i
end if
end do
end do
end do
!
!
do m=1,2
do j=1,nx2
do i=1,nsetten
if (narax(i)==j+(m-1)*(nx1+nx2)&
&.and.naray(i)<(ny1*2+ny2+ny4*2+ny5+1)&
&.and.naray(i)>(ny1*2+ny2)&
&.and.naraz(i)>=nza-nzhaji-nzhaji2&
&.and.naraz(i)<nza+1) then
nyouk=nyouk+1
meshkou(nyouk,1)=i
end if
end do
end do
end do
!
!
do l=1,nosi*2-1
do m=1,2
do j=1,nx2
do i=1,nsetten
if (narax(i)==j+(m-1)*(nx1+nx2)&
&.and.naray(i)<(ny1*2+ny2+ny4*2+ny5+1)&
&.and.naray(i)>(ny1*2+ny2)&
&.and.naraz(i)>nzhaji+nzhaji2+(l-1)*(nzitahan*2&
&+nzana)&
&.and.naraz(i)<nzhaji+nzhaji2+nzitahan*2&
&+(l-1)*(nzitahan*2+nzana)) then
nyouk=nyouk+1
meshkou(nyouk,1)=i
end if
end do
end do
end do
end do
!
do j=1,nsetten
do i=1,nyouk
if (naraz(j)==naraz(meshkou(i,1))+1&
&.and.narax(j)==narax(meshkou(i,1))&
&.and.naray(j)==naray(meshkou(i,1))) then
meshkou(i,4)=j
else if (naraz(j)==naraz(meshkou(i,1))&
&.and.narax(j)==narax(meshkou(i,1))&
&.and.naray(j)==naray(meshkou(i,1))+1) then
meshkou(i,2)=j
else if (naraz(j)==naraz(meshkou(i,1))+1&
&.and.narax(j)==narax(meshkou(i,1))&
&.and.naray(j)==naray(meshkou(i,1))+1) then
meshkou(i,3)=j
end if
end do
end do
!
do i=1,nyouk
meshkou(i,5)=meshkou(i,1)+1
meshkou(i,6)=meshkou(i,2)+1
meshkou(i,7)=meshkou(i,3)+1
meshkou(i,8)=meshkou(i,4)+1
end do
!
!
!print'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")'
!do i=1,nyoum
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&i,',',&
!&meshmoku(i,1),',',&
!&meshmoku(i,2),',',&
!&meshmoku(i,3),',',&
!&meshmoku(i,4),',',&
!&meshmoku(i,5),',',&
!&meshmoku(i,6),',',&
!&meshmoku(i,7),',',&
!&meshmoku(i,8)
!end do
!print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")'
!do i=1,nyouk
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&i+nyoum,',',&
!&meshkou(i,1),',',&
!&meshkou(i,2),',',&
!&meshkou(i,3),',',&
!&meshkou(i,4),',',&
!&meshkou(i,5),',',&
!&meshkou(i,6),',',&
!&meshkou(i,7),',',&
!&meshkou(i,8)
!end do
!共有していない節点をたして要素のなかを変える
do i=1,nsetten
narax(i)=-1
naray(i)=-1
naraz(i)=-1
end do
iset=nsetten
do j=1,nosi
do i=1,nsetten
if&
&(z(i)>=ellhaji+ellhaji2+(ellitahan*2.0d0+ellaida)*(j-1)-0.00001d0&
&.and.z(i)<=ellhaji+ellhaji2+ellitahan*2.0d0&
&+(ellitahan*2.0d0+ellaida)*(j-1)+0.00001d0)&
then
naraz(i)=-2
end if 
end do
end do
!
do i=1,nsetten
if (y(i)<=h1+h2+0.00001d0&
&.and.y(i)>=h1-0.00001d0&
&)then
naray(i)=-3
else if (y(i)<=h1+h2+h1+h5+h1+h2+0.0001d0&
&.and.y(i)>=h1+h2+h1+h5+h1-0.0001d0&
&) then
naray(i)=-3
end if
if (y(i)<=h1+h2+0.00001d0&
&.and.y(i)>=h1-0.00001d0&
&)then
naray(i)=-3
else if (y(i)<=h1+h2+h1+h5+h1+h2+0.0001d0&
&.and.y(i)>=h1+h2+h1+h5+h1-0.0001d0&
&) then
naray(i)=-3
end if
!
if (x(i)<=b2+0.00001d0&
&.and.x(i)>=b2-0.00001d0&
&)then
narax(i)=-4
else if (x(i)<=b1+b2+0.00001d0&
&.and.x(i)>=b1+b2-0.00001d0&
&)then
narax(i)=-5
end if
end do
!
!do i=1,nsetten
!if (narax(i)==-4) then
!print*,i,narax(i)
!end if
!end do
!
do i=1,nsetten
if (naray(i)==-3.and.naraz(i)==-2.and.narax(i)==-4) then
narax(i)=0
else if (naray(i)==-3.and.naraz(i)==-2.and.narax(i)==-5) then
narax(i)=-2
end if
end do
!
!print'("*NODE, NSET=Nall")'
!do i=1,nsetten
!if (narax(i)==-4) then
!print*,i,narax(i)
!end if
!end do
!print'("*NODE, NSET=Nall")'
!do i=1,nsetten
!if (narax(i)==-5) then
!print*,i,narax(i)
!end if
!end do
!print'("*NODE, NSET=Nall")'
!do i=1,nsetten
!if (narax(i)==0) then
!print*,i,narax(i)
!end if
!end do
!do i=1,nsetten
!if (narax(i)==-4) then
!print*,i,narax(i)
!end if
!end do
!
iset=nsetten
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,1))==-4) then
iset=iset+1
narax(meshmoku(i,1))=iset
x(iset)=b2
y(iset)=y(meshmoku(i,1))
z(iset)=z(meshmoku(i,1))
meshmoku(i,1)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,2))==-4) then
iset=iset+1
narax(meshmoku(i,2))=iset
x(iset)=b2
y(iset)=y(meshmoku(i,2))
z(iset)=z(meshmoku(i,2))
meshmoku(i,2)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,3))==-4) then
iset=iset+1
narax(meshmoku(i,3))=iset
x(iset)=b2
y(iset)=y(meshmoku(i,3))
z(iset)=z(meshmoku(i,3))
meshmoku(i,3)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,4))==-4) then
iset=iset+1
narax(meshmoku(i,4))=iset
x(iset)=b2
y(iset)=y(meshmoku(i,4))
z(iset)=z(meshmoku(i,4))
meshmoku(i,4)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,5))==-5) then
iset=iset+1
narax(meshmoku(i,5))=iset
x(iset)=b2+b1
y(iset)=y(meshmoku(i,5))
z(iset)=z(meshmoku(i,5))
meshmoku(i,5)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,6))==-5) then
iset=iset+1
narax(meshmoku(i,6))=iset
y(iset)=y(meshmoku(i,6))
z(iset)=z(meshmoku(i,6))
x(iset)=b2+b1
meshmoku(i,6)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,7))==-5) then
iset=iset+1
narax(meshmoku(i,7))=iset
x(iset)=b2+b1
y(iset)=y(meshmoku(i,7))
z(iset)=z(meshmoku(i,7))
meshmoku(i,7)=iset
end if
end do
do i=1,nx1*(ny1*2+ny2)*nza*2
if (narax(meshmoku(i,8))==-5) then
iset=iset+1
narax(meshmoku(i,8))=iset
x(iset)=b2+b1
y(iset)=y(meshmoku(i,8))
z(iset)=z(meshmoku(i,8))
meshmoku(i,8)=iset
end if
end do
!do i=1,nsetten
!if (narax(i)==-4) then
!print*,i,narax(i)
!end if
!end do
!do i=1,nsetten
!if (narax(i)==-5) then
!print*,i,narax(i)
!end if
!end do
!節点の書く
!鋼板節点
print'("*NODE, NSET=Nall")'
do i=1,iset
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
&i,',',x(i),',',y(i),',',z(i)
end do
!要素を書く
print'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")'
do i=1,nyoum
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i,',',&
&meshmoku(i,1),',',&
&meshmoku(i,2),',',&
&meshmoku(i,3),',',&
&meshmoku(i,4),',',&
&meshmoku(i,5),',',&
&meshmoku(i,6),',',&
&meshmoku(i,7),',',&
&meshmoku(i,8)
end do

print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")'
do i=1,nyouk
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i+nyoum,',',&
&meshkou(i,1),',',&
&meshkou(i,2),',',&
&meshkou(i,3),',',&
&meshkou(i,4),',',&
&meshkou(i,5),',',&
&meshkou(i,6),',',&
&meshkou(i,7),',',&
&meshkou(i,8)
end do
! 拘束節点の節点番号
print'("*NSET,NSET=kousoku1")'
do i=1,iset
if (z(i)<=ellhaji+0.000001d0.and.z(i)>=ellhaji-0.000001d0&
&.and.y(i)<=0.000001.and.y(i)>=-0.000001) then
print'(i6,a)',i,','
end if
end do
print'("*NSET,NSET=kousoku2")'
do i=1,iset
if (z(i)<=ellall-ellhaji+0.000001d0.and.z(i)>=ellall-ellhaji-0.000001d0&
&.and.y(i)<=0.000001.and.y(i)>=-0.000001) then
print'(i6,a)',i,','
end if
end do

!! 載荷点の節点番号
print'("*NSET,NSET=saika")'
do i=1,iset
if (z(i)<=ellall/2.0d0+0.000001d0.and.z(i)>=ellall/2.0d0-0.000001d0&
&.and.y(i)<=h4+0.000001d0.and.y(i)>=h4-0.000001d0&
&.and.x(i)>=b2+xn1-0.00001d0.and.&
&x(i)<=b2+b1-xn1+0.000001d0) then
print'(i6,a)',i,','
end if 
end do
! 載荷点の下の節点番号
!print'("*NSET,NSET=saikasita")'
print'("*NSET,NSET=saikasita")'
do i=1,iset
if (z(i)<=ellall/2.0d0+0.000001d0.and.z(i)>=ellall/2.0d0-0.000001d0&
&.and.y(i)<=0.000001d0.and.y(i)>=-0.000001d0) then
print'(i6,a)',i,','
end if 
end do
print'("*NSET,NSET=saikamen")'
do i=1,iset
if (z(i)<=ellall/2.0d0+0.000001d0.and.z(i)>=ellall/2.0d0-0.000001d0) then
print'(i6,a)',i,','
end if 
end do
! 境界条件
print'("*BOUNDARY")'
print'("kousoku1,1,3")'
print'("kousoku2,1,2")'
! 材料の設定　ヤング率
! 鋼板の設定
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("206.d3, 0.3")'
print'("*SOLID SECTION,ELSET=Ekou,MATERIAL=EL1")'
! 木材の設定
print'("*MATERIAL,NAME=EL2")'
print'("*ELASTIC,TYPE=ORTHOROPIC")'
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
print'(f12.4,a,f12.4,a,f12.4,a,&
&f12.4,a,f12.4,a,f12.4,a,f12.4,&
& a ,f12.4,a)',&
&e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
&,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
print'(f12.4,a,a)',e66(4,4),',','293'
print'("*SOLID SECTION,ELSET=Emoku,MATERIAL=EL2")'
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*BUCKLE")'
print'("5,0.01")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/(nx1-1)
! datのオプション
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
print'("*EL PRINT,ELSET=Eall")'
print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE,POSITION=AVERAGED AT NODES")'
print'("S")'
print'("*END STEP")'
!
!
!
         end
!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3)
!応力－ひずみ行列：e66
do i=1,6
do j=1,6
e66(i,j)=0.d0
end do
end do
!  
e66(1,1)=1.d0/exx
e66(1,2)=-poixy/exx
e66(1,3)=-poixz/exx
e66(2,1)=-poiyx/eyy
e66(2,2)=1.d0/eyy
e66(2,3)=-poiyz/eyy
e66(3,1)=-poizx/ezz
e66(3,2)=-poizy/ezz
e66(3,3)=1.d0/ezz
!
!****** 逆行列を取って、応力－ひずみ行列に変換する******
det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
&+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
&-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
do i=1,3
do j=1,3
e66(i,j)=b(j,i)/det
end do
end do
!
!****** せん断補正係数を入力 ******
!print*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!print*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!print*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
!
!
return
end
