-
Notifications
You must be signed in to change notification settings - Fork 5
/
element.f90
171 lines (156 loc) · 6.69 KB
/
element.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
module element
use connectivity
use initialization
implicit none
real(kind=8),dimension(1,2)::gauss
real(kind=8),dimension(4,2)::Gauss_4
real(kind=8),dimension(2,2)::jaco,inv_J
real(kind=8),dimension(1,4)::N
real(kind=8),dimension(2,4)::N_ab,N_xy
real(kind=8),dimension(3,8)::B
real(kind=8),dimension(8,3)::B_trans
real(kind=8)::det_jacobian,area
real(kind=8),dimension(:),allocatable::Mass
real(kind=8),dimension(8,1)::U_el,F_el
real(kind=8),dimension(8,8)::K_el_int, K_el
contains
function Gauss_integration_points() result(Gauss_4)
real(kind=8),dimension(4,2):: Gauss_4
! integer::i,j
Gauss_4 = reshape((/-1,1,1,-1,-1,-1,1,1/),shape(Gauss_4))
! do i=1,size(Gauss_4,1)
! print*, (Gauss_4(i,j),j=1,size(Gauss_4,2))
! enddo
Gauss_4 = 1/sqrt(3.0)*Gauss_4
end function Gauss_integration_points
subroutine get_deformed_node_coord_parallel(node,en,U_dt)
real(kind=8),dimension(4,2),intent(inout)::node
type(disp_type),intent(in)::U_dt
integer,dimension(1,4),intent(in)::en
integer,dimension(4)::buffer
buffer=en(1,:)
node(:,1) = node(:,1) + U_dt%values(buffer,1)
node(:,2) = node(:,2) + U_dt%values(buffer,2)
end subroutine get_deformed_node_coord_parallel
subroutine get_deformed_node_coord(node,en,U_dt)
real(kind=8),dimension(4,2),intent(inout)::node
type(disp_type),intent(in)::U_dt
integer,dimension(1,4),intent(in)::en
integer,dimension(4)::buffer
buffer=en(1,:)
node(:,1) = node(:,1) + U_dt%values(buffer,1)
node(:,2) = node(:,2) + U_dt%values(buffer,2)
end subroutine get_deformed_node_coord
subroutine shape_function(gauss,N)
real(kind=8),dimension(1,2),intent(in)::gauss
real(kind=8),dimension(1,4),intent(out)::N
real(kind=8)::a,b
real(kind=8),parameter::one=1.0d0
a = gauss(1,1)
b = gauss(1,2)
N(1,1) = 0.25*(one-a)*(one-b)
N(1,2) = 0.25*(one+a)*(one-b)
N(1,3) = 0.25*(one+a)*(one+b)
N(1,4) = 0.25*(one-a)*(one+b)
end subroutine shape_function
function element_area(node) result(area)
real(kind=8),dimension(4,2),intent(in)::node
real(kind=8)::area
real(kind=8),dimension(1,2)::gauss
real(kind=8),dimension(4,2)::gauss_4
real(kind=8)::det_jacobian
integer::intp
gauss_4 = Gauss_integration_points()
area = 0.0
do intp = 1,4
gauss(1,1:2) = gauss_4(intp,1:2)
call Jacobian(gauss,node,jaco,det_jacobian)
!print*, det_jacobian
area = area + det_jacobian
enddo
end function element_area
subroutine create_mass_matrix(element_data,node_data, &
node_nums, element_nums,rho,Mass)
integer::el_id
integer,intent(in)::element_nums,node_nums
type(node_type),intent(in)::node_data
type(element_type),intent(in)::element_data
integer,dimension(1,4) ::en
integer,dimension(4)::en2
real(kind=8),dimension(4,2)::node
real(kind=8)::area, mass_per_node
real(kind=8),intent(in)::rho
real(kind=8),dimension(:),allocatable,intent(out)::Mass
allocate(Mass(2*node_nums))
Mass(:) = 0.d0
!print*,size(Mass)
do el_id=1,element_nums
call element_info(el_id,element_data,node_data,en,node)
area = element_area(node)
mass_per_node = area*rho/4.d0
en2 = en(1,1:4)
Mass(2*en2-1) = Mass(2*en2-1) + mass_per_node
Mass(2*en2) = Mass(2*en2) + mass_per_node
enddo
end subroutine create_mass_matrix
subroutine Jacobian(gauss,node,jaco,det_jacobian)
real(kind=8),dimension(1,2),intent(in)::gauss
real(kind=8),dimension(4,2),intent(in)::node
real(kind=8),dimension(2,2),intent(out)::jaco
real(kind=8),intent(out)::det_jacobian
real(kind=8),dimension(2,4)::dphi
real(kind=8)::a,b
real(kind=8),parameter::one=1.0d0
! integer::i,j,k
a = gauss(1,1)
b = gauss(1,2)
dphi(1,1:4) = 0.25*(/b-one,one-b,one+b,-one-b/)
dphi(2,1:4) = 0.25*(/a-one,-one-a,one+a,one-a/)
jaco=0.0
call dgemm('n','n',2,2,4,1.d0,dphi,2,node,4,0.d0,jaco,2)
det_jacobian = jaco(1,1)*jaco(2,2)-jaco(1,2)*jaco(2,1)
! print *, det_jacobian
end subroutine Jacobian
function get_N_ab(gauss) result(dphi)
real(kind=8),dimension(1,2),intent(in)::gauss
real(kind=8),dimension(2,4)::dphi
real(kind=8)::a,b
real(kind=8),parameter::one=1.0d0
a = gauss(1,1)
b = gauss(1,2)
dphi(1,1:4) = 0.25*(/b-one,one-b,one+b,-one-b/)
dphi(2,1:4) = 0.25*(/a-one,-one-a,one+a,one-a/)
end function get_N_ab
subroutine get_B(jacobian,N_ab,B,inv_J)
real(kind=8),dimension(2,2),intent(in)::jacobian
real(kind=8),dimension(2,4),intent(in)::N_ab
real(kind=8),dimension(2,2),intent(out)::inv_J
real(kind=8),dimension(2,4)::N_xy
real(kind=8),dimension(3,8),intent(out)::B
real(kind=8)::det_jacobian
integer::i,j,k
N_xy(:,:) = 0.0
!-- First need to inverse Jacobian
det_jacobian = jacobian(1,1)*jacobian(2,2)-&
jacobian(1,2)*jacobian(2,1)
! do i=1,2
! write(*,'(2F16.10)') (jacobian(i,j),j=1,2)
! enddo
inv_J(1,1) = jacobian(2,2)/det_jacobian
inv_J(2,2) = jacobian(1,1)/det_jacobian
inv_J(1,2) = -jacobian(1,2)/det_jacobian
inv_J(2,1) = -jacobian(2,1)/det_jacobian
!-- Then multiply inv_J with N_ab
do j=1,4
do i=1,2
do k=1,2
N_xy(i,j) = N_xy(i,j) + inv_J(i,k)*N_ab(k,j)
enddo
enddo
enddo
B(1,1:7:2) = N_xy(1,1:4)
B(2,2:8:2) = N_xy(2,1:4)
B(3,1:7:2) = N_xy(2,1:4)
B(3,2:8:2) = N_xy(1,1:4)
end subroutine get_B
end module element