FrontISTR  5.2.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 !-------------------------------------------------------------------------------
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  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 
270 end 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