!CLT縺ｮ阯､逕ｰ1螻､讀懃ｮ�(CCX)
! 蠑ｾ蝪第ｧ隗｣譫舌∝濠隗｣譫舌�
! 蠅�阜譚｡莉ｶ繧偵�縺｣縺溘ｊ縺ｫ蜷医ｏ縺帙ｋ繧医≧縺ｫ螟峨∴縺溘�
! 螟我ｽ崎ｨ医�菴咲ｽｮ縺ｧ螟我ｽ阪ｒ蜃ｺ蜉帙☆繧九ｈ縺�↓螟峨∴縺溘�(霈芽差菴咲ｽｮ縺ｯ謇句虚縺ｧ螟峨∴繧九∋縺�)
implicit none
!#--蜈･蜉�
real    :: b1=3380.d-3    ! m   :x霆ｸ(讖玖ｻｸ)譁ｹ蜷代�髟ｷ縺�1(H蠖｢驪ｼ縺ｮ髟ｷ縺�)
real    :: b2=310.d-3     ! m   :x霆ｸ(讖玖ｻｸ)譁ｹ蜷代�髟ｷ縺�2(荳｡縺ゅ∪繧�)
real    :: h=240.d-3      ! m   :蜈ｨ螻､縺ｮ鬮倥＆(y霆ｸ)
real    :: ell1=860.d-3   ! m   :z霆ｸ(謾ｯ髢�)譁ｹ蜷代�髟ｷ縺�1
real    :: ell2=60.d-3    ! m   :z霆ｸ(謾ｯ髢�)譁ｹ蜷代�髟ｷ縺�2(繝懊Ν繝�~繧ｴ繝�)
real    :: ell3=80.d-3    ! m   :z霆ｸ(謾ｯ髢�)譁ｹ蜷代�髟ｷ縺�3(遶ｯ~繝懊Ν繝�)
integer :: nsou=2         !     :螻､謨ｰ
integer :: nx1=50         !     :x霆ｸ(讖玖ｻｸ)譁ｹ蜷大�蜑ｲ謨ｰ1
integer :: nx2=3          !     :x霆ｸ(讖玖ｻｸ)譁ｹ蜷大�蜑ｲ謨ｰ2
integer :: ny=12          !     :y霆ｸ(鬮倥＆)譁ｹ蜷大�蜑ｲ謨ｰ
integer :: nz1=50         !     :z霆ｸ(謾ｯ髢�)譁ｹ蜷大�蜑ｲ謨ｰ1
integer :: nz2=3          !     :z霆ｸ(謾ｯ髢�)譁ｹ蜷大�蜑ｲ謨ｰ2
integer :: nz3=3          !     :z霆ｸ(謾ｯ髢�)譁ｹ蜷大�蜑ｲ謨ｰ3
real    :: p=1000.d-3/2   ! MN  :闕ｷ驥�/2
integer :: nstep=100      !     :蜃ｺ蜉帙せ繝�ャ繝玲焚
real    :: saikax=200.d-3 ! m   :霈芽差譚ｿ縺ｮ蟷�(x)
real    :: saikaz=500.d-3 ! m   :霈芽差譚ｿ縺ｮ蟷�(z)
real    :: itix=2000.d-3  ! m   :讖玖ｻｸ譁ｹ蜷�(x)縺ｮ霈芽差譚ｿ荳ｭ蠢�ｽ咲ｽｮ
!real    :: itix=1000.d-3  ! m   :讖玖ｻｸ譁ｹ蜷�(x)縺ｮ霈芽差譚ｿ荳ｭ蠢�ｽ咲ｽｮ
real    :: itiz=0.d-3     ! m   :謾ｯ髢捺婿蜷�(z)縺ｮ霈芽差譚ｿ荳ｭ蠢�ｽ咲ｽｮ
real    :: exx=1.7d3      ! MPa :x霆ｸ(讖玖ｻｸ)譁ｹ蜷代�繝､繝ｳ繧ｰ邇�
real    :: eyy=5.0d3/25   ! MPa :y霆ｸ(鬮倥＆)譁ｹ蜷代�繝､繝ｳ繧ｰ邇�
real    :: ezz=3.5d3      ! MPa :z霆ｸ(謾ｯ髢�)譁ｹ蜷代�繝､繝ｳ繧ｰ邇�
real    :: poizx=0.4d0    !     :繝昴い繧ｽ繝ｳ豈�
!---
real poixy,poixz,poiyx,poiyz,poizy,gxy,gxz,gyz
real :: x(1000000),y(1000000),z(1000000) !遽轤ｹ蠎ｧ讓�
integer :: mesh(1000000,8)               !隕∫ｴ�縺ｨ隕∫ｴ�繧呈ｧ区�縺吶ｋ遽轤ｹ
real :: e66(6,6),b(3,3),dx1,dx2,dy,dz1,dz2,dz3,ts
integer i,j,k,l,nx,nz,nset,nyou,nsaika
!---繝昴い繧ｽ繝ｳ豈�
poizy=poizx
poixy=poizx*exx/ezz
poixz=poixy
poiyx=poizx*eyy/ezz
poiyz=poiyx
!---縺帙ｓ譁ｭ蠑ｾ諤ｧ菫よ焚
gxy=5.0d3/15
gxz=gxy
gyz=gxy
!---隕∫ｴ�髟ｷ縺�
dx1=b1/nx1
dx2=b2/nx2
dy=h/(ny*nsou)
dz1=ell1/nz1
dz2=ell2/nz2
dz3=ell3/nz3
nx=nx1+nx2*2
nz=nz1+nz2+nz3
!---遽轤ｹ蠎ｧ讓吶�險育ｮ�
nset=0
do k=0,nz1
   do j=0,ny*nsou
      do i=0,nx2
         nset=1+nset
         x(nset)=i*dx2
         y(nset)=j*dy
         z(nset)=k*dz1
      enddo
      do i=1,nx1
         nset=1+nset
         x(nset)=b2+i*dx1
         y(nset)=j*dy
         z(nset)=k*dz1
      enddo
      do i=1,nx2
         nset=1+nset
         x(nset)=b1+b2+i*dx2
         y(nset)=j*dy
         z(nset)=k*dz1
      enddo
   enddo
enddo
do k=1,nz2
   do j=0,ny*nsou
      do i=0,nx2
         nset=1+nset
         x(nset)=i*dx2
         y(nset)=j*dy
         z(nset)=ell1+k*dz2
      enddo
      do i=1,nx1
         nset=1+nset
         x(nset)=b2+i*dx1
         y(nset)=j*dy
         z(nset)=ell1+k*dz2
      enddo
      do i=1,nx2
         nset=1+nset
         x(nset)=b1+b2+i*dx2
         y(nset)=j*dy
         z(nset)=ell1+k*dz2
      enddo
   enddo
enddo
do k=1,nz3
   do j=0,ny*nsou
      do i=0,nx2
         nset=1+nset
         x(nset)=i*dx2
         y(nset)=j*dy
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx1
         nset=1+nset
         x(nset)=b2+i*dx1
         y(nset)=j*dy
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx2
         nset=1+nset
         x(nset)=b1+b2+i*dx2
         y(nset)=j*dy
         z(nset)=ell1+ell2+k*dz3
      enddo
   enddo
enddo
!---隕∫ｴ�繧呈ｧ区�縺吶ｋ遽轤ｹ縺ｮ險育ｮ�
nyou=0
do k=0,nz-1
   do j=0,ny*nsou-1
      do i=1,nx
         nyou=1+nyou
         mesh(nyou,1)=i+j*(nx+1)+k*((nx+1)*(ny*nsou+1))
         mesh(nyou,2)=mesh(nyou,1)+(nx+1)
         mesh(nyou,3)=mesh(nyou,2)+(nx+1)*(ny*nsou+1)
         mesh(nyou,4)=mesh(nyou,1)+(nx+1)*(ny*nsou+1)
         mesh(nyou,5)=mesh(nyou,1)+1
         mesh(nyou,6)=mesh(nyou,5)+nx+1
         mesh(nyou,7)=mesh(nyou,6)+(nx+1)*(ny*nsou+1)
         mesh(nyou,8)=mesh(nyou,5)+(nx+1)*(ny*nsou+1)
      enddo
   enddo
enddo
!---遽轤ｹ蠎ｧ讓吶�蜃ｺ蜉�
print'("*NODE, NSET=Nall")'
do i=1,nset
   print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
       &i,',',x(i),',',y(i),',',z(i)
end do
!---隕∫ｴ�縺ｮ蜃ｺ蜉�
do k=0,nsou-1
   print'(A,i1)',"*ELEMENT, TYPE=C3D8, ELSET=E",k+1
   do j=0,nz-1
      do i=1,nx*ny
         l=i+j*(nx*ny*nsou)+k*(nx*ny)
         print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
             &l,',',&
             &mesh(l,1),',',&
             &mesh(l,2),',',&
             &mesh(l,3),',',&
             &mesh(l,4),',',&
             &mesh(l,5),',',&
             &mesh(l,6),',',&
             &mesh(l,7),',',&
             &mesh(l,8)
      enddo
   enddo
enddo
!#--蠅�阜譚｡莉ｶ縺ｮ遽轤ｹ鄒､
print'("*NSET,NSET=kotei1")'!H蠖｢驪ｼ縺ｮ蜈ｨ髱｢
do i=nz1+1,nz+1
   do j=nx2+1,nx2+nx1+1
      l=(nx+1)*(ny*nsou+1)*i-(nx+1)+j
      print'(i6,a)',l
   end do
enddo
print'("*NSET,NSET=kotei2")'!H蠖｢驪ｼ縺ｮ蜊企擇
do i=nz1+1,nz1+nz2+1
   do j=nx2+1,nx2+nx1+1
      l=(nx+1)*(ny*nsou+1)*i-(nx+1)+j
      print'(i6,a)',l
   end do
enddo
print'("*NSET,NSET=kotei3")'!繝懊Ν繝医�邱壻ｸ�
i=nz1+nz2+1
do j=nx2+1,nx2+nx1+1
   l=(nx+1)*(ny*nsou+1)*i-(nx+1)+j
   print'(i6,a)',l
end do
print'("*NSET,NSET=koteiz")'!z霆ｸ蟇ｾ遘ｰ髱｢
do i=1,nset
   if (z(i)>-dz1/2.and.z(i)<dz1/2) then
      print'(i6,a)',i
   endif
end do
!---霈芽差譚｡莉ｶ縺ｮ遽轤ｹ鄒､
print'("*NSET,NSET=saika")'
nsaika=0
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>itiz-saikaz/2-dz1/3.and.z(i)<itiz+saikaz/2+dz1/3.and.&
       &x(i)>itix-saikax/2-dx1/3.and.x(i)<itix+saikax/2+dx1/3) then
      nsaika=nsaika+1
      print'(i6,a)',i
   endif
end do
!#--螟我ｽ崎ｨ医�遽轤ｹ鄒､
print'("*NSET,NSET=hen1")'!螟我ｽ崎ｨ�1
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>250.d-3-dx2/2.and.x(i)<250.d-3+dx2/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen2")'!螟我ｽ崎ｨ�2
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>1000.d-3-dx1/2.and.x(i)<1000.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen3")'!螟我ｽ崎ｨ�3
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>1750.d-3-dx1/2.and.x(i)<1750.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen7")'!螟我ｽ崎ｨ�7
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>2000.d-3-dx1/2.and.x(i)<2000.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen4")'!螟我ｽ崎ｨ�4
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>2250.d-3-dx1/2.and.x(i)<2250.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen5")'!螟我ｽ崎ｨ�5
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>3000.d-3-dx1/2.and.x(i)<3000.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen6")'!螟我ｽ崎ｨ�6
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>3750.d-3-dx2/2.and.x(i)<3750.d-3+dx2/2) then
      print'(i6,a)',i
   endif
end do
!---蠅�阜譚｡莉ｶ
print'("*BOUNDARY")'
print'("kotei2,1,2")'
print'("koteiz,3")'!z霆ｸ蟇ｾ遘ｰ髱｢
!#--譚先侭險ｭ螳�
print'("*MATERIAL,NAME=ELt")'!蠑募ｼｵ
print'("*ELASTIC,TYPE=ORTHO")'
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)',e66(4,4)
print'("*PLASTIC")'
print'("54.88,0.0")'
print'("54.88,0.5")'
print'("*MATERIAL,NAME=ELc")'!蝨ｧ邵ｮ
print'("*ELASTIC,TYPE=ORTHO")'
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)',e66(4,4)
print'("*PLASTIC")'
print'("22.54,0.0")'
print'("27.44,0.5")'
!#--螻謇蠎ｧ讓咏ｳｻ縺ｮ險ｭ螳�
print'("*ORIENTATION,NAME=OR1,SYSTEM=RECTANGULAR")'
print'("1.,0.,0.,0.,1.,0.")'!蜈ｨ菴灘ｺｧ讓咏ｳｻ縺ｨ荳閾ｴ
print'("*ORIENTATION,NAME=OR2,SYSTEM=RECTANGULAR")'
print'("0.,0.,-1.,0.,1.,0.")'!y霆ｸ蝗槭ｊ縺ｫ90蠎ｦ蝗櫁ｻ｢
print'("*SOLID SECTION,ELSET=E1,MATERIAL=ELc,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E2,MATERIAL=ELt,ORIENTATION=OR1")'
!print'("*SOLID SECTION,ELSET=E3,MATERIAL=EL1,ORIENTATION=OR1")'
!print'("*SOLID SECTION,ELSET=E4,MATERIAL=EL1,ORIENTATION=OR2")'
!print'("*SOLID SECTION,ELSET=E5,MATERIAL=EL1,ORIENTATION=OR2")'
!print'("*SOLID SECTION,ELSET=E6,MATERIAL=EL1,ORIENTATION=OR1")'
!print'("*SOLID SECTION,ELSET=E7,MATERIAL=EL1,ORIENTATION=OR2")'
!print'("*SOLID SECTION,ELSET=E8,MATERIAL=EL1,ORIENTATION=OR1")'
!---蜃ｺ蜉帙せ繝�ャ繝玲焚
print'("*TIME POINTS,NAME=T1")'
do i=1,nstep
   ts=real(i)/nstep
   print'(f12.4,a)',ts,','
enddo
!---霈芽差譚｡莉ｶ
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nsaika
!print'("*FREQUENCY")'!蝗ｺ譛牙､隗｣譫�
!print'("10")'!謖ｯ蜍輔Δ繝ｼ繝�
!#--dat繝輔ぃ繧､繝ｫ縺ｮ險ｭ螳�
print'("*NODE PRINT,NSET=hen1,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�1
print'("U")'
print'("*NODE PRINT,NSET=hen2,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�2
print'("U")'
print'("*NODE PRINT,NSET=hen3,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�3
print'("U")'
print'("*NODE PRINT,NSET=hen7,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�7
print'("U")'
print'("*NODE PRINT,NSET=hen4,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�4
print'("U")'
print'("*NODE PRINT,NSET=hen5,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�5
print'("U")'
print'("*NODE PRINT,NSET=hen6,TIMEPOINTS=T1")'!螟我ｽ崎ｨ�6
print'("U")'
print'("*NODE PRINT,NSET=saika,TIMEPOINTS=T1")'!蟷ｳ蝮�喧縺励◆縺ｻ縺�′縺�＞��
print'("U")'
!print'("*EL PRINT,ELSET=")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE")'
print'("S")'
print'("*END STEP")'
!
end

!-------- 逶ｴ莠､逡ｰ譁ｹ諤ｧ縺ｮ螳壽焚繧貞�縺吶�繝ｭ繧ｰ繝ｩ繝� --------
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit none
real :: e66(6,6),b(3,3)
real exx,eyy,ezz,poixy,poixz,poiyz,gyz,gxz,gxy,poiyx,poizx,poizy
real det
integer i,j
!implicit real*8(a-h,o-z)
! e66(6,6),b(3,3)
!蠢懷鴨�阪�縺壹∩陦悟��啼66
!do i=1,6
!do j=1,6
!e66(i,j)=0.d0
!end do
!end do
e66=0.d0
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