FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
dynamic_mass.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
7
8contains
9
10 subroutine mass_c2(etype, nn, ecoord, gausses, sec_opt, thick, mass, lumped, temperature)
11 use mmechgauss
12 use m_matmatrix
13 use elementinfo
14 implicit none
15 type(tgaussstatus), intent(in) :: gausses(:)
16 integer(kind=kint), intent(in) :: etype
17 integer(kind=kint), intent(in) :: nn
18 real(kind=kreal), intent(in) :: ecoord(2,nn)
19 real(kind=kreal), intent(out) :: mass(:,:)
20 real(kind=kreal), intent(out) :: lumped(:)
21 real(kind=kreal), intent(in), optional :: temperature(nn)
22 type(tmaterial), pointer :: matl
23 integer(kind=kint), parameter :: ndof = 2
24 integer(kind=kint) :: i, j, lx, sec_opt
25 real(kind=kreal) :: naturalcoord(2)
26 real(kind=kreal) :: func(nn), thick
27 real(kind=kreal) :: det, wg, rho
28 real(kind=kreal) :: d(2,2), n(2, nn*ndof), dn(2, nn*ndof)
29 real(kind=kreal) :: gderiv(nn,2)
30 logical :: is_lumped
31
32 mass(:,:) = 0.0d0
33 lumped = 0.0d0
34 matl => gausses(1)%pMaterial
35
36 if(sec_opt == 2) thick = 1.0d0
37
38 do lx = 1, numofquadpoints(etype)
39 call getquadpoint(etype, lx, naturalcoord)
40 call getshapefunc(etype, naturalcoord, func)
41 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
42
43 if(present(temperature))then
44 !ina(1) = temperature
45 !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
46 !if(ierr)then
47 rho = matl%variables(m_density)
48 !else
49 ! rho = outa(1)
50 !endif
51 else
52 !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
53 !if(ierr)then
54 rho = matl%variables(m_density)
55 !else
56 ! rho = outa(1)
57 !endif
58 endif
59
60 d = 0.0d0
61 d(1,1) = rho*thick
62 d(2,2) = rho*thick
63
64 wg = getweight(etype, lx)*det
65
66 n = 0.0d0
67 do i = 1, nn
68 n(1,2*i-1) = func(i)
69 n(2,2*i ) = func(i)
70 enddo
71
72 dn(1:2, 1:nn*ndof) = matmul(d, n(1:2, 1:nn*ndof))
73 forall(i = 1:nn*ndof, j = 1:nn*ndof)
74 mass(i,j) = mass(i,j) + dot_product(n(:,i), dn(:,j))*wg
75 end forall
76 enddo
77
78 is_lumped = .true.
79 if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
80 end subroutine mass_c2
81
82 subroutine mass_c3(etype, nn, ecoord, gausses, mass, lumped, temperature)
83 use mmechgauss
84 use m_matmatrix
85 use elementinfo
86 implicit none
87 type(tgaussstatus), intent(in) :: gausses(:)
88 integer(kind=kint), intent(in) :: etype
89 integer(kind=kint), intent(in) :: nn
90 real(kind=kreal), intent(in) :: ecoord(3,nn)
91 real(kind=kreal), intent(out) :: mass(:,:)
92 real(kind=kreal), intent(out) :: lumped(:)
93 real(kind=kreal), intent(in), optional :: temperature(nn)
94 type(tmaterial), pointer :: matl
95 integer(kind=kint), parameter :: ndof = 3
96 integer(kind=kint) :: i, j, lx
97 real(kind=kreal) :: naturalcoord(3)
98 real(kind=kreal) :: func(nn)
99 real(kind=kreal) :: det, wg, rho
100 real(kind=kreal) :: d(3, 3), n(3, nn*ndof), dn(3, nn*ndof)
101 real(kind=kreal) :: gderiv(nn, 3)
102 logical :: is_lumped
103
104 mass(:,:) = 0.0d0
105 lumped = 0.0d0
106 matl => gausses(1)%pMaterial
107
108 do lx = 1, numofquadpoints(etype)
109 call getquadpoint(etype, lx, naturalcoord)
110 call getshapefunc(etype, naturalcoord, func)
111 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
112
113 if(present(temperature))then
114 !ina(1) = temperature
115 !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
116 !if(ierr)then
117 rho = matl%variables(m_density)
118 !else
119 ! rho = outa(1)
120 !endif
121 else
122 !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
123 !if(ierr)then
124 rho = matl%variables(m_density)
125 !else
126 ! rho = outa(1)
127 !endif
128 endif
129
130 d = 0.0d0
131 d(1,1) = rho
132 d(2,2) = rho
133 d(3,3) = rho
134
135 wg = getweight(etype,lx)*det
136
137 n = 0.0d0
138 do i = 1, nn
139 n(1,3*i-2) = func(i)
140 n(2,3*i-1) = func(i)
141 n(3,3*i ) = func(i)
142 enddo
143
144 dn(1:3, 1:nn*ndof) = matmul(d, n(1:3, 1:nn*ndof))
145 forall(i = 1:nn*ndof, j = 1:nn*ndof)
146 mass(i,j) = mass(i,j) + dot_product(n(:,i), dn(:,j))*wg
147 end forall
148 enddo
149
150 is_lumped = .true.
151 if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
152 end subroutine mass_c3
153
154 subroutine get_lumped_mass(nn, ndof, mass, lumped)
155 use hecmw
156 implicit none
157 integer(kind=kint) :: i, j, nn, ndof
158 real(kind=kreal) :: lumped(:), mass(:,:)
159 real(kind=kreal) :: diag_mass, total_mass
160
161 total_mass = 0.0d0
162 do i = 1, nn*ndof, ndof
163 do j = 1, nn*ndof, ndof
164 total_mass = total_mass + mass(j,i)
165 enddo
166 enddo
167
168 diag_mass = 0.0d0
169 do i = 1, nn*ndof, ndof
170 diag_mass = diag_mass + mass(i,i)
171 enddo
172
173 diag_mass = 1.0d0/diag_mass
174 do i = 1, nn*ndof
175 lumped(i) = lumped(i) + mass(i,i)*total_mass*diag_mass
176 enddo
177
178 mass = 0.0d0
179 do i = 1, nn*ndof
180 mass(i,i) = lumped(i)
181 enddo
182 end subroutine get_lumped_mass
183
184 function get_length(ecoord)
185 use hecmw
186 implicit none
187 real(kind=kreal) :: get_length, ecoord(3,20)
188
189 get_length = dsqrt( &
190 (ecoord(1,2) - ecoord(1,1))**2 + &
191 (ecoord(2,2) - ecoord(2,1))**2 + &
192 (ecoord(3,2) - ecoord(3,1))**2 )
193 end function get_length
194
195 function get_face3(ecoord)
196 use hecmw
197 implicit none
198 real(kind=kreal) :: get_face3, ecoord(3,20)
199 real(kind=kreal) :: a1, a2, a3
200 real(kind=kreal) :: x(3), y(3), z(3)
201
202 x(1) = ecoord(1,1); y(1) = ecoord(2,1); z(1) = ecoord(3,1)
203 x(2) = ecoord(1,2); y(2) = ecoord(2,2); z(2) = ecoord(3,2)
204 x(3) = ecoord(1,3); y(3) = ecoord(2,3); z(3) = ecoord(3,3)
205
206 a1 = (x(2) - x(1))**2 + (y(2) - y(1))**2 + (z(2) - z(1))**2
207 a2 = (x(1) - x(3))*(x(2) - x(1)) &
208 & + (y(1) - y(3))*(y(2) - y(1)) &
209 & + (z(1) - z(3))*(z(2) - z(1))
210 a3 = (x(3) - x(1))**2 + (y(3) - y(1))**2 + (z(3) - z(1))**2
211
212 get_face3 = 0.5d0*dsqrt(a1*a3 - a2*a2)
213 end function get_face3
214
215 function get_face4(ecoord)
216 use hecmw
217 implicit none
218 integer(kind=kint) :: lx, ly
219 real(kind=kreal) :: get_face4, ecoord(3,20)
220 real(kind=kreal) :: xg(2), ri, si, rp, sp, rm, sm, hr(4), hs(4)
221 real(kind=kreal) :: xr, xs, yr, ys, zr, zs
222 real(kind=kreal) :: x(4), y(4), z(4), det
223
224 x(1) = ecoord(1,1); y(1) = ecoord(2,1); z(1) = ecoord(3,1)
225 x(2) = ecoord(1,2); y(2) = ecoord(2,2); z(2) = ecoord(3,2)
226 x(3) = ecoord(1,3); y(3) = ecoord(2,3); z(3) = ecoord(3,3)
227 x(4) = ecoord(1,4); y(4) = ecoord(2,4); z(4) = ecoord(3,4)
228
229 xg(1) = -0.5773502691896258d0
230 xg(2) = -xg(1)
231 get_face4 = 0.0d0
232
233 do lx = 1, 2
234 ri = xg(lx)
235 do ly = 1, 2
236 si = xg(ly)
237 rp = 1.0d0 + ri
238 sp = 1.0d0 + si
239 rm = 1.0d0 - ri
240 sm = 1.0d0 - si
241
242 !C* FOR R-COORDINATE
243 hr(1) = 0.25d0*sp
244 hr(2) = -0.25d0*sp
245 hr(3) = -0.25d0*sm
246 hr(4) = 0.25d0*sm
247
248 !C* FOR S-COORDINATE
249 hs(1) = 0.25d0*rp
250 hs(2) = 0.25d0*rm
251 hs(3) = -0.25d0*rm
252 hs(4) = -0.25d0*rp
253
254 !C*JACOBI MATRIX
255 xr = hr(1)*x(1) + hr(2)*x(2) + hr(3)*x(3) + hr(4)*x(4)
256 xs = hs(1)*x(1) + hs(2)*x(2) + hs(3)*x(3) + hs(4)*x(4)
257 yr = hr(1)*y(1) + hr(2)*y(2) + hr(3)*y(3) + hr(4)*y(4)
258 ys = hs(1)*y(1) + hs(2)*y(2) + hs(3)*y(3) + hs(4)*y(4)
259 zr = hr(1)*z(1) + hr(2)*z(2) + hr(3)*z(3) + hr(4)*z(4)
260 zs = hs(1)*z(1) + hs(2)*z(2) + hs(3)*z(3) + hs(4)*z(4)
261
262 det = (yr*zs - zr*ys)**2 + (zr*xs - xr*zs)**2 + (xr*ys - yr*xs)**2
263 det = dsqrt(det)
264
265 get_face4 = get_face4 + det
266 enddo
267 enddo
268 end function get_face4
269
270end module m_dynamic_mass
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
This module contains subroutines used in 3d eigen analysis for.
Definition: dynamic_mass.f90:6
real(kind=kreal) function get_length(ecoord)
subroutine mass_c2(etype, nn, ecoord, gausses, sec_opt, thick, mass, lumped, temperature)
real(kind=kreal) function get_face4(ecoord)
real(kind=kreal) function get_face3(ecoord)
subroutine mass_c3(etype, nn, ecoord, gausses, mass, lumped, temperature)
subroutine get_lumped_mass(nn, ndof, mass, lumped)
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13