fstr_EIG_setMASS.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_fstr_eig_setmass
 
  use m_eigen_lib
  use m_static_get_prop
  use m_dynamic_mass
  use m_fstr
 
contains
 
  subroutine setmass(fstrSOLID, hecMESH, hecMAT, fstrEIG)
    use hecmw_util
    use m_fstr
    use mmaterial
    use m_static_lib_shell
    implicit none
    type(hecmwst_matrix)     :: hecMAT
    type(hecmwst_local_mesh) :: hecMESH
    type(fstr_solid)         :: fstrSOLID
    type(fstr_eigen)         :: fstrEIG
    integer(kind=kint) :: i, iS, iE, ii, nn, j, jn, jS, k
    integer(kind=kint) :: N, NP, NDOF
    integer(kind=kint) :: icel, ic_type, itype, isect, ihead, sec_opt, cid
    integer(kind=kint) :: nodLOCAL(20)
    real(kind=kreal) :: val, surf, chkmass
    real(kind=kreal) :: rho, thick, length
    real(kind=kreal) :: lumped(120), mass(20*6, 20*6), ecoord(3,20)
 
    if(myrank == 0)then
      write(imsg,*)"* mass matrix generation"
    endif
 
    n  = hecmat%N
    np = hecmat%NP
    ndof = hecmesh%n_dof
 
    allocate(fstreig%mass(np*ndof))
    fstreig%mass = 0.0d0
    val    = 0.0d0
 
    do itype = 1, hecmesh%n_elem_type
      is = hecmesh%elem_type_index(itype-1) + 1
      ie = hecmesh%elem_type_index(itype  )
      ic_type = hecmesh%elem_type_item(itype)
      nn = hecmw_get_max_node(ic_type)
 
      if(hecmw_is_etype_link(ic_type)) cycle
      if(hecmw_is_etype_patch(ic_type)) cycle
      if(ic_type == 3414) cycle
      if(ic_type == 881 .or. ic_type == 891) cycle  !skip selective ES/NS smoothed elements (mass uses original 341)
 
      do icel = is, ie
        js = hecmesh%elem_node_index(icel-1)
        do j = 1, nn
          nodlocal(j) = hecmesh%elem_node_item(js+j)
          do i = 1, 3
            ecoord(i,j) = hecmesh%node(3*(nodlocal(j)-1)+i)
          enddo
        enddo
 
        isect = hecmesh%section_ID(icel)
        ihead = hecmesh%section%sect_R_index(isect-1)
        cid   = hecmesh%section%sect_mat_ID_item(isect)
        sec_opt = hecmesh%section%sect_opt(isect)
        rho   = fstrsolid%materials(cid)%variables(m_density)
        thick = fstrsolid%materials(cid)%variables(m_thick)
 
        lumped = 0.0d0
        if(ic_type == 231 .or. ic_type == 232 .or. ic_type == 241 .or. ic_type == 242)then
          call mass_c2(ic_type, nn, ecoord(1:2,1:nn), fstrsolid%elements(icel)%gausses, sec_opt, thick, mass, lumped)
 
        elseif(ic_type == 341 .or. ic_type == 342 .or. ic_type == 351 .or. ic_type == 352 .or. &
             & ic_type == 361 .or. ic_type == 362 )then
          call mass_c3(ic_type, nn, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses, mass, lumped)
 
        elseif(ic_type==731 .or. ic_type==741 .or. ic_type==743) then
          rho = fstrsolid%materials(cid)%variables(m_density)
          thick = fstrsolid%materials(cid)%variables(m_thick)
          call mass_shell(ic_type, nn, ecoord(1:3,1:nn), rho, thick, fstrsolid%elements(icel)%gausses, mass, lumped)
 
        elseif(ic_type == 761)then
          surf = get_face3(ecoord(1:3,1:nn))
          rho = fstrsolid%materials(cid)%variables(m_density)
          thick = fstrsolid%materials(cid)%variables(m_thick)
          val = surf*thick*rho/3.0d0
 
        elseif(ic_type == 781)then
          surf = get_face4(ecoord(1:3,1:nn))
          rho = fstrsolid%materials(cid)%variables(m_density)
          thick = fstrsolid%materials(cid)%variables(m_thick)
          val = surf*thick*rho/4.0d0
 
        elseif(ic_type == 611 .or. ic_type == 641)then
          surf = hecmesh%section%sect_R_item(ihead+4)
          length = get_length(ecoord(1:3,1:nn))
          rho = fstrsolid%materials(cid)%variables(m_density)
          val = 0.5d0*surf*length*rho
 
        elseif(ic_type == 301)then
          surf = hecmesh%section%sect_R_item(ihead+1)
          length = get_length(ecoord(1:3,1:nn))
          rho = fstrsolid%materials(cid)%variables(m_density)
          val = 0.5d0*surf*length*rho
 
        elseif( ic_type/100 == 5 )then !skip interface element
        else
          write(*,*)"** error setMASS"
        endif
 
        do j = 1,nn
          jn = nodlocal(j)
          js = ndof*(jn-1)
 
          if(ic_type == 611)then
            fstreig%mass(js+1) = fstreig%mass(js+1) + val
            fstreig%mass(js+2) = fstreig%mass(js+2) + val
            fstreig%mass(js+3) = fstreig%mass(js+3) + val
            fstreig%mass(js+4) = fstreig%mass(js+4) + 0.0d0
            fstreig%mass(js+5) = fstreig%mass(js+5) + 0.0d0
            fstreig%mass(js+6) = fstreig%mass(js+6) + 0.0d0
 
          elseif(ic_type == 641)then
            if(j == 1 .or. j == 2)then
              fstreig%mass(js+1) = fstreig%mass(js+1) + val
              fstreig%mass(js+2) = fstreig%mass(js+2) + val
              fstreig%mass(js+3) = fstreig%mass(js+3) + val
            elseif(j == 3 .or. j == 4)then
              fstreig%mass(js+1) = fstreig%mass(js+1) + 0.0d0
              fstreig%mass(js+2) = fstreig%mass(js+2) + 0.0d0
              fstreig%mass(js+3) = fstreig%mass(js+3) + 0.0d0
            end if
 
          elseif(ic_type == 761)then
            if(j == 1 .or. j == 2 .or. j == 3)then
              fstreig%mass(js+1) = fstreig%mass(js+1) + val
              fstreig%mass(js+2) = fstreig%mass(js+2) + val
              fstreig%mass(js+3) = fstreig%mass(js+3) + val
            elseif(j == 4 .or. j == 5 .or. j == 6)then
              fstreig%mass(js+1) = fstreig%mass(js+1) + 0.0d0
              fstreig%mass(js+2) = fstreig%mass(js+2) + 0.0d0
              fstreig%mass(js+3) = fstreig%mass(js+3) + 0.0d0
            endif
 
          else if( ic_type == 781 ) then
            if(j == 1 .or. j == 2 .or. j == 3 .or. j == 4)then
              fstreig%mass(js+1) = fstreig%mass(js+1) + val
              fstreig%mass(js+2) = fstreig%mass(js+2) + val
              fstreig%mass(js+3) = fstreig%mass(js+3) + val
            elseif(j == 5 .or. j == 6 .or. j == 7 .or. j == 8)then
              fstreig%mass(js+1) = fstreig%mass(js+1) + 0.0d0
              fstreig%mass(js+2) = fstreig%mass(js+2) + 0.0d0
              fstreig%mass(js+3) = fstreig%mass(js+3) + 0.0d0
            endif
 
          elseif(ic_type == 301) then
            fstreig%mass(js+1) = fstreig%mass(js+1) + val
            fstreig%mass(js+2) = fstreig%mass(js+2) + val
            fstreig%mass(js+3) = fstreig%mass(js+3) + val
 
          else
            do k = 1, ndof
              fstreig%mass(js+k) = fstreig%mass(js+k) + lumped(ndof*(j-1)+k)
            enddo
          endif
        enddo
      enddo
    enddo
 
    call hecmw_update_r(hecmesh, fstreig%mass, np, ndof)
 
    chkmass = 0.0d0
    do i = 1, n
      ii = ndof*(i-1) + 1
      chkmass = chkmass + fstreig%mass(ii)
    end do
    call hecmw_allreduce_r1(hecmesh, chkmass, hecmw_sum)
    fstreig%totalmass = chkmass
 
    if(myrank == 0)then
      write(imsg,"(a,1pe12.5)")"** Total mass: ", chkmass
    endif
  end subroutine setmass
 
end module m_fstr_eig_setmass