dynamic_mass.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_dynamic_mass
  use m_fstr
 
contains
 
  subroutine mass_c2(etype, nn, ecoord, gausses, sec_opt, thick, mass, lumped, temperature)
    use mmechgauss
    use m_matmatrix
    use elementinfo
    implicit none
    type(tgaussstatus), intent(in) :: gausses(:)
    integer(kind=kint), intent(in) :: etype
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in)  :: ecoord(2,nn)           
    real(kind=kreal), intent(out) :: mass(:,:)              
    real(kind=kreal), intent(out) :: lumped(:)              
    real(kind=kreal), intent(in), optional :: temperature(nn) 
    type(tmaterial), pointer :: matl
    integer(kind=kint), parameter :: ndof = 2
    integer(kind=kint) :: i, j, lx, sec_opt
    real(kind=kreal) :: naturalcoord(2)
    real(kind=kreal) :: func(nn), thick
    real(kind=kreal) :: det, wg, rho
    real(kind=kreal) :: d(2,2), n(2, nn*ndof), dn(2, nn*ndof)
    real(kind=kreal) :: gderiv(nn,2)
    logical :: is_lumped
 
    mass(:,:) = 0.0d0
    lumped = 0.0d0
    matl => gausses(1)%pMaterial
 
    if(sec_opt == 2) thick = 1.0d0
 
    do lx = 1, numofquadpoints(etype)
      call getquadpoint(etype, lx, naturalcoord)
      call getshapefunc(etype, naturalcoord, func)
      call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
 
      if(present(temperature))then
        !ina(1) = temperature
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
        !if(ierr)then
          rho = matl%variables(m_density)
        !else
        !  rho = outa(1)
        !endif
      else
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
        !if(ierr)then
          rho = matl%variables(m_density)
        !else
        !  rho = outa(1)
        !endif
      endif
 
      d = 0.0d0
      d(1,1) = rho*thick
      d(2,2) = rho*thick
 
      wg = getweight(etype, lx)*det
 
      n = 0.0d0
      do i = 1, nn
        n(1,2*i-1) = func(i)
        n(2,2*i  ) = func(i)
      enddo
 
      dn(1:2, 1:nn*ndof) = matmul(d, n(1:2, 1:nn*ndof))
      do  j = 1,nn*ndof
        do i = 1,nn*ndof
          mass(i,j) = mass(i,j) + dot_product(n(:,i), dn(:,j))*wg
        enddo
      enddo
    enddo
 
    is_lumped = .true.
    if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
  end subroutine mass_c2
 
  subroutine mass_c3(etype, nn, ecoord, gausses, mass, lumped, temperature)
    use mmechgauss
    use m_matmatrix
    use elementinfo
    implicit none
    type(tgaussstatus), intent(in) :: gausses(:)
    integer(kind=kint), intent(in) :: etype
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in)  :: ecoord(3,nn)           
    real(kind=kreal), intent(out) :: mass(:,:)              
    real(kind=kreal), intent(out) :: lumped(:)              
    real(kind=kreal), intent(in), optional :: temperature(nn) 
    type(tmaterial), pointer :: matl
    integer(kind=kint), parameter :: ndof = 3
    integer(kind=kint) :: i, j, lx
    real(kind=kreal) :: naturalcoord(3)
    real(kind=kreal) :: func(nn)
    real(kind=kreal) :: det, wg, rho
    real(kind=kreal) :: d(3, 3), n(3, nn*ndof), dn(3, nn*ndof)
    real(kind=kreal) :: gderiv(nn, 3)
    logical :: is_lumped
 
    mass(:,:) = 0.0d0
    lumped = 0.0d0
    matl => gausses(1)%pMaterial
 
    do lx = 1, numofquadpoints(etype)
      call getquadpoint(etype, lx, naturalcoord)
      call getshapefunc(etype, naturalcoord, func)
      call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
 
      if(present(temperature))then
        !ina(1) = temperature
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
        !if(ierr)then
          rho = matl%variables(m_density)
        !else
        !  rho = outa(1)
        !endif
      else
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
        !if(ierr)then
          rho = matl%variables(m_density)
        !else
        !  rho = outa(1)
        !endif
      endif
 
      d = 0.0d0
      d(1,1) = rho
      d(2,2) = rho
      d(3,3) = rho
 
      wg = getweight(etype,lx)*det
 
      n = 0.0d0
      do i = 1, nn
        n(1,3*i-2) = func(i)
        n(2,3*i-1) = func(i)
        n(3,3*i  ) = func(i)
      enddo
 
      dn(1:3, 1:nn*ndof) = matmul(d, n(1:3, 1:nn*ndof))
      do  j = 1,nn*ndof
        do i = 1,nn*ndof
          mass(i,j) = mass(i,j) + dot_product(n(:,i), dn(:,j))*wg
        enddo
      enddo
    enddo
 
    is_lumped = .true.
    if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
  end subroutine mass_c3
 
  subroutine mass_s3(surf, thick, rho, mass)
    implicit none
    real(kind=kreal) :: surf
    real(kind=kreal) :: thick
    real(kind=kreal) :: rho
    real(kind=kreal), intent(out) :: mass(:,:)
    integer(kint) :: i
    mass = 0.0d0
    do i = 1, 9
      mass(i,i) = surf*thick*rho/3.0d0
    enddo
  end subroutine mass_s3
 
  subroutine mass_s4(surf, thick, rho, mass)
    implicit none
    real(kind=kreal) :: surf
    real(kind=kreal) :: thick
    real(kind=kreal) :: rho
    real(kind=kreal), intent(out) :: mass(:,:)
    integer(kint) :: i
    mass = 0.0d0
    do i = 1, 12
      mass(i,i) = surf*thick*rho/4.0d0
    enddo
  end subroutine mass_s4
 
  subroutine mass_beam(surf, length, rho, mass)
    implicit none
    real(kind=kreal) :: surf
    real(kind=kreal) :: length
    real(kind=kreal) :: rho
    real(kind=kreal), intent(out) :: mass(:,:)
    integer(kint) :: i
    mass = 0.0d0
    do i = 1, 12
      mass(i,i) = 0.5d0*surf*length*rho
    enddo
  end subroutine mass_beam
 
  subroutine mass_beam_33(surf, length, rho, mass)
    implicit none
    real(kind=kreal) :: surf
    real(kind=kreal) :: length
    real(kind=kreal) :: rho
    real(kind=kreal), intent(out) :: mass(:,:)
    integer(kint) :: i
    mass = 0.0d0
    do i = 1, 12
      mass(i,i) = 0.5d0*surf*length*rho
    enddo
  end subroutine mass_beam_33
 
  subroutine get_lumped_mass(nn, ndof, mass, lumped)
    use hecmw
    implicit none
    integer(kind=kint) :: i, j, nn, ndof
    real(kind=kreal) :: lumped(:), mass(:,:)
    real(kind=kreal) :: diag_mass, total_mass
 
    total_mass = 0.0d0
    do i = 1, nn*ndof, ndof
      do j = 1, nn*ndof, ndof
        total_mass = total_mass + mass(j,i)
      enddo
    enddo
 
    diag_mass = 0.0d0
    do i = 1, nn*ndof, ndof
      diag_mass = diag_mass + mass(i,i)
    enddo
 
    diag_mass = 1.0d0/diag_mass
    do i = 1, nn*ndof
      lumped(i) = lumped(i) + mass(i,i)*total_mass*diag_mass
    enddo
 
    mass = 0.0d0
    do i = 1, nn*ndof
      mass(i,i) = lumped(i)
    enddo
  end subroutine get_lumped_mass
 
  function get_length(ecoord)
    use hecmw
    implicit none
    real(kind=kreal) :: get_length, ecoord(3,20)
 
    get_length = dsqrt( &
      (ecoord(1,2) - ecoord(1,1))**2 + &
      (ecoord(2,2) - ecoord(2,1))**2 + &
      (ecoord(3,2) - ecoord(3,1))**2 )
  end function get_length
 
  function get_face3(ecoord)
    use hecmw
    implicit none
    real(kind=kreal) :: get_face3, ecoord(3,20)
    real(kind=kreal) :: a1, a2, a3
    real(kind=kreal) :: x(3), y(3), z(3)
 
    x(1) = ecoord(1,1); y(1) = ecoord(2,1); z(1) = ecoord(3,1)
    x(2) = ecoord(1,2); y(2) = ecoord(2,2); z(2) = ecoord(3,2)
    x(3) = ecoord(1,3); y(3) = ecoord(2,3); z(3) = ecoord(3,3)
 
    a1 = (x(2) - x(1))**2 + (y(2) - y(1))**2 + (z(2) - z(1))**2
    a2 = (x(1) - x(3))*(x(2) - x(1)) &
     & + (y(1) - y(3))*(y(2) - y(1)) &
     & + (z(1) - z(3))*(z(2) - z(1))
    a3 = (x(3) - x(1))**2 + (y(3) - y(1))**2 + (z(3) - z(1))**2
 
    get_face3 = 0.5d0*dsqrt(a1*a3 - a2*a2)
  end function get_face3
 
  function get_face4(ecoord)
    use hecmw
    implicit none
    integer(kind=kint) :: lx, ly
    real(kind=kreal) :: get_face4, ecoord(3,20)
    real(kind=kreal) :: xg(2), ri, si, rp, sp, rm, sm, hr(4), hs(4)
    real(kind=kreal) :: xr, xs, yr, ys, zr, zs
    real(kind=kreal) :: x(4), y(4), z(4), det
 
    x(1) = ecoord(1,1); y(1) = ecoord(2,1); z(1) = ecoord(3,1)
    x(2) = ecoord(1,2); y(2) = ecoord(2,2); z(2) = ecoord(3,2)
    x(3) = ecoord(1,3); y(3) = ecoord(2,3); z(3) = ecoord(3,3)
    x(4) = ecoord(1,4); y(4) = ecoord(2,4); z(4) = ecoord(3,4)
 
    xg(1) = -0.5773502691896258d0
    xg(2) = -xg(1)
    get_face4 = 0.0d0
 
    do lx = 1, 2
      ri = xg(lx)
      do ly = 1, 2
        si = xg(ly)
        rp = 1.0d0 + ri
        sp = 1.0d0 + si
        rm = 1.0d0 - ri
        sm = 1.0d0 - si
 
        !C*  FOR R-COORDINATE
        hr(1) =  0.25d0*sp
        hr(2) = -0.25d0*sp
        hr(3) = -0.25d0*sm
        hr(4) =  0.25d0*sm
 
        !C*  FOR S-COORDINATE
        hs(1) =  0.25d0*rp
        hs(2) =  0.25d0*rm
        hs(3) = -0.25d0*rm
        hs(4) = -0.25d0*rp
 
        !C*JACOBI MATRIX
        xr = hr(1)*x(1) + hr(2)*x(2) + hr(3)*x(3) + hr(4)*x(4)
        xs = hs(1)*x(1) + hs(2)*x(2) + hs(3)*x(3) + hs(4)*x(4)
        yr = hr(1)*y(1) + hr(2)*y(2) + hr(3)*y(3) + hr(4)*y(4)
        ys = hs(1)*y(1) + hs(2)*y(2) + hs(3)*y(3) + hs(4)*y(4)
        zr = hr(1)*z(1) + hr(2)*z(2) + hr(3)*z(3) + hr(4)*z(4)
        zs = hs(1)*z(1) + hs(2)*z(2) + hs(3)*z(3) + hs(4)*z(4)
 
        det = (yr*zs - zr*ys)**2 + (zr*xs - xr*zs)**2 + (xr*ys - yr*xs)**2
        det = dsqrt(det)
 
        get_face4 = get_face4 + det
      enddo
    enddo
  end function get_face4
 
end module m_dynamic_mass