module single_K
    contains
   
subroutine solve(n_num,coods,E,nu,t,K)
implicit none
integer::n_num(3)       !单元节点号,i,j,m
real::coods(6)          !单元三个节点的坐标,i,j,m
real::E,nu,t              !弹性模量和泊松比
real::K(6,6)            !单元刚度矩阵
real::A                 !单元面积
real::B(3,6)            !应变转换矩阵
real::D(3,3)            !弹性矩阵
real::S(3,6)            !应力转换矩阵
integer::i,j

call A_area(coods,A)
call B_maxtri(coods,A,B)
call D_maxtri(E,nu,D)
call S_maxtri(D,B,S)
call K_maxtri(B,S,K)

do i=1,6,1
    do j=1,6,1
        K(i,j)=t*A*K(i,j)
    end do
end do

write(*,*) 'BT='
do i=1,6,1
    do j=1,6,1
        write(*,"(F15.5)",ADVANCE='NO') K(i,j)
    end do
    write(*,"(/)",ADVANCE='NO')
end do
end subroutine

subroutine A_area(coods,A)
implicit none
real::coods(6)
real::A
A=coods(3)*coods(6)+coods(1)*coods(4)+coods(5)*coods(2)-coods(3)*coods(2)-coods(5)*coods(4)-coods(1)*coods(6)
A=A/2.0
write(*,"('Area=',F20.5)") A
end subroutine

subroutine B_maxtri(coods,A,B)
implicit none
real::coods(6)
real::B(3,6)
real::B_ijm(3,2)        !B矩阵的分量
real::temp_i(2),temp_j(2),temp_m(2)
real::ai,bi,ci
real::A
integer::i,j

do i=1,3,1
    do j=1,2,1
        B_ijm(i,j)=0
    end do
end do

do i=1,3,1
    ai=coods(3)*coods(6)-coods(5)*coods(3)
    bi=coods(4)-coods(6)
    ci=coods(5)-coods(3)
   
    temp_i=coods(1:2)
    temp_j=coods(3:4)
    temp_m=coods(5:6)
   
    coods(1:2)=temp_j
    coods(3:4)=temp_m
    coods(5:6)=temp_i
   
    B_ijm(1,1)=bi/2.0/A
    B_ijm(2,2)=ci/2.0/A
    B_ijm(3,1)=ci/2.0/A
    B_ijm(3,2)=bi/2.0/A
   
    B(1:3,(2*i-1):(2*i))=B_ijm
end do

do i=1,3,1
    do j=1,6,1
        write(*,"(F15.5)",ADVANCE='NO') B(i,j)
    end do
    write(*,"(/)",ADVANCE='NO')
end do

end subroutine

subroutine D_maxtri(E,nu,D)
implicit none
real::E
real::nu
real::D(3,3)
integer::i,j

do i=1,3,1
    do j=1,3,1
        if(j==i) then
            D(i,j)=1.0
        else
            D(i,j)=0
        endif
    end do
end do

D(1,2)=nu
D(2,1)=nu
D(3,3)=(1-nu)/2.0

do i=1,3,1
    do j=1,3,1
        D(i,j)=E*D(i,j)/(1-nu**2)
    end do
end do

write(*,*) 'D='
do i=1,3,1
    do j=1,3,1
        write(*,"(F15.5)",ADVANCE='NO') D(i,j)
    end do
    write(*,"(/)",ADVANCE='NO')
end do

end subroutine

subroutine S_maxtri(D,B,S)
implicit none
real::D(3,3)
real::B(3,6)
real::S(3,6)
integer::i,j,m

do i=1,3,1
    do j=1,6,1
        S(i,j)=0
        do m=1,3,1
            S(i,j)=S(i,j)+D(i,m)*B(m,j)
        end do
    end do
end do

write(*,*) 'S='
do i=1,3,1
    do j=1,6,1
        write(*,"(F15.5)",ADVANCE='NO') S(i,j)
    end do
    write(*,"(/)",ADVANCE='NO')
end do

end subroutine

subroutine K_maxtri(B,S,K)
implicit none
real::B(3,6)
real::S(3,6)
real::K(6,6)
real::BT(6,3)
real::i,j,m

do i=1,3,1
    do j=1,6,1
        BT(j,i)=B(i,j)
    end do
end do

write(*,*) 'BT='

do i=1,6,1
    do j=1,3,1
        write(*,"(F15.5)",ADVANCE='NO') BT(i,j)
    end do
    write(*,"(/)",ADVANCE='NO')
end do

do i=1,6,1
    do j=1,6,1
        K(i,j)=0
        do m=1,3,1
            K(i,j)=K(i,j)+BT(i,m)*S(m,j)
        end do
    end do
end do

end subroutine
end module

program main
use single_K
implicit none
integer::n_num(3)=(/1,2,3/)
real::coods(6)=(/0,2000,0,1000,1000,1000/)
real::E=2.0e5
real::nu=0
real::t=10.0
real::K(6,6)

call solve(n_num,coods,E,nu,t,K)
end program