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