fstr_CreateMatrix.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_creatematrix_and_dampingforce
  use m_fstr
  implicit none
 
  private
  public :: fstr_creatematrix_and_dampingforce
 
contains
 
  !---------------------------------------------------------------------*
  subroutine fstr_creatematrix_and_dampingforce( hecMESH, hecMAT, fstrSOLID, time, tincr, fstrDYNAMIC, coef )
  !---------------------------------------------------------------------*
    use m_static_lib
    use mmechgauss
    use m_dynamic_mass
    use m_elemact
    use m_fstr_nodalkinematics, only: fstr_ensure_finite_rotation_state
 
    type(hecmwst_local_mesh)              :: hecmesh
    type(hecmwst_matrix)                  :: hecmat
    type(fstr_solid)                      :: fstrsolid
    real(kind=kreal), intent(in)          :: time         
    real(kind=kreal), intent(in)          :: tincr        
    type(fstr_dynamic), intent(in), optional :: fstrdynamic
    real(kind=kreal), intent(in), optional :: coef(6)     
 
    integer(kind=kint) :: ndof, itype, is, ie, ic_type, nn, icel, iis, i, j, in, jn
    integer(kind=kint) :: nodlocal(fstrsolid%max_ncon)
    real(kind=kreal)   :: stiff_mat(fstrsolid%max_ncon_stf*6, fstrsolid%max_ncon_stf*6)
    real(kind=kreal)   :: mass_mat(20*6, 20*6), damp_mat(20*6, 20*6)
    real(kind=kreal)   :: mat(20*6, 20*6), lumped(20*6)
    real(kind=kreal)   :: tt(fstrsolid%max_ncon), ecoord(3,fstrsolid%max_ncon)
    real(kind=kreal)   :: u(6,fstrsolid%max_ncon), du(6*20), u_prev(6,fstrsolid%max_ncon)
    real(kind=kreal)   :: veca(20*6), vecb(20*6)
    real(kind=kreal)   :: acc(20*6), vec(20*6), df(20*6)
    real(kind=kreal)   :: a1, a2, a3, b1, b2, b3
    type(tmaterial), pointer :: material
    logical            :: is_dynamic
 
    is_dynamic = present(fstrdynamic) .and. present(coef)
 
    call hecmw_mat_clear( hecmat )
    if( is_dynamic ) then
      fstrsolid%DFORCE = 0.0d0
      a1 = coef(1); a2 = coef(2); a3 = coef(3)
      b1 = coef(4); b2 = coef(5); b3 = coef(6)
    endif
 
    ndof = hecmat%NDOF
    call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
 
    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)
 
      ! ----- Ignore link and patch elements
      if (hecmw_is_etype_link(ic_type)) cycle
      if (hecmw_is_etype_patch(ic_type)) cycle
 
      nn = hecmw_get_max_node(ic_type)
 
      !$omp parallel default(none), &
        !$omp&  private(icel,iiS,nn,j,nodLOCAL,i,in,jn,ecoord,du,u,u_prev,tt, &
        !$omp&          material,stiff_mat,mass_mat,damp_mat, &
        !$omp&          lumped,mat,df,vec,acc,vecA,vecB), &
        !$omp&  shared(iS,iE,hecMESH,ndof,fstrSOLID,ic_type,hecMAT,time,tincr,fstrDYNAMIC, &
        !$omp&         is_dynamic,a1,a2,a3,b1,b2,b3)
      !$omp do
      do icel = is, ie
 
        ! ----- nodal coordinate & displacement
        iis = hecmesh%elem_node_index(icel-1)
        nn  = hecmesh%elem_node_index(icel) - iis
        tt(:) = 0.d0
        do j = 1, nn
          in = hecmesh%elem_node_item(iis+j)
          nodlocal(j) = in
          do i = 1, 3
            ecoord(i,j) = hecmesh%node(3*(in-1)+i)
          enddo
          do i = 1, ndof
            du(ndof*(j-1)+i) = fstrsolid%dunode(ndof*(in-1)+i)
            u(i,j) = fstrsolid%unode(ndof*(in-1)+i) + du(ndof*(j-1)+i)
            u_prev(i,j) = fstrsolid%unode(ndof*(in-1)+i)
          enddo
          if( is_dynamic ) then
            do i = 1, ndof
              vec(ndof*(j-1)+i) = fstrdynamic%VEL(ndof*(in-1)+i,1)
              acc(ndof*(j-1)+i) = fstrdynamic%ACC(ndof*(in-1)+i,1)
            enddo
          endif
          if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
            tt(j) = fstrsolid%temperature( nodlocal(j) )
          endif
        enddo
 
        ! ----- inactive element : assemble dummy stiffness and skip the rest
        if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) then
          stiff_mat = 0.0d0
          call stf_dummy( ndof, nn, ecoord(:,1:nn), u(1:3,1:nn), &
            &  stiff_mat(1:nn*ndof, 1:nn*ndof), fstrsolid%elements(icel) )
          call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiff_mat)
          cycle
        endif
 
        ! ----- element K, M, C (section / material context resolved inside)
        stiff_mat = 0.0d0
        mass_mat  = 0.0d0
        damp_mat  = 0.0d0
 
        call calc_stiff_and_mass_elem( ic_type, nn, ndof, ecoord, u, u_prev, tt, &
          time, tincr, is_dynamic, fstrsolid, hecmesh, icel, nodlocal, &
          stiff_mat, mass_mat, damp_mat, lumped )
 
        ! ----- assemble system matrix (and the dynamic damping force)
        if( is_dynamic ) then
          ! Newmark intermediate vectors
          do i = 1, nn*ndof
            veca(i) = -a3*du(i) + a2*vec(i) + a1*acc(i)
            vecb(i) = -b3*du(i) + b2*vec(i) + b1*acc(i)
          enddo
 
          ! effective dynamic system matrix A = c1*K + c2*M (+ b3*C) and damping force df
          material => fstrsolid%elements(icel)%gausses(1)%pMaterial
          call calc_damping_mat_and_force_elem( ic_type, nn, ndof, material, fstrdynamic, &
            a3, b3, veca, vecb, stiff_mat, mass_mat, damp_mat, mat, df )
 
          call hecmw_mat_ass_elem(hecmat, nn, nodlocal, mat)
 
          do j = 1, nn
            do i = 1, ndof
              !$omp atomic
              fstrsolid%DFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%DFORCE(ndof*(nodlocal(j)-1)+i)+df(ndof*(j-1)+i)
            enddo
          enddo
        else
          ! static / eigen / frequency : assemble K only
          call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiff_mat)
        endif
 
      enddo      ! icel
      !$omp end do
      !$omp end parallel
    enddo        ! itype
 
  end subroutine fstr_creatematrix_and_dampingforce
 
  !---------------------------------------------------------------------*
  subroutine calc_stiff_and_mass_elem( ic_type, nn, ndof, ecoord, u, u_prev, tt, &
      time, tincr, is_dynamic, fstrSOLID, hecMESH, icel, nodLOCAL, &
      stiff_mat, mass_mat, damp_mat, lumped )
  !---------------------------------------------------------------------*
    use m_static_lib
    use m_static_lib_shell, only: shellcomposenodaldisplacement
    use mmechgauss
    use m_dynamic_mass
    use m_fstr_nodalkinematics, only: fstr_get_shell_trial_directors, fstr_get_shell_reference_directors
    use m_fstr_finiterotationkinematics, only: fstr_uses_finite_rotation_kinematics
 
    integer(kind=kint), intent(in)    :: ic_type, ndof
    integer(kind=kint), intent(inout) :: nn  ! some legacy STF_* routines mutate this
    real(kind=kreal),   intent(in)    :: ecoord(:,:), u(:,:), u_prev(:,:), tt(:)
    real(kind=kreal),   intent(in)    :: time, tincr
    logical,            intent(in)    :: is_dynamic
    type(fstr_solid),               intent(inout) :: fstrSOLID
    type(hecmwst_local_mesh),       intent(in)    :: hecMESH
    integer(kind=kint),             intent(in)    :: icel
    integer(kind=kint),             intent(inout) :: nodLOCAL(:)
    real(kind=kreal),   intent(inout) :: stiff_mat(:,:), mass_mat(:,:), damp_mat(:,:)
    real(kind=kreal),   intent(inout) :: lumped(:)
 
    type(tmaterial), pointer :: material
    integer(kind=kint) :: isect, ihead, cdsys_id, sec_opt, i, j
    real(kind=kreal) :: coords(3,3), thick
    real(kind=kreal) :: rho, length, surf
    real(kind=kreal) :: shell_nddisp(6,fstrsolid%max_ncon)
    real(kind=kreal) :: shell_du(6,fstrsolid%max_ncon)
    real(kind=kreal) :: shell_director(3,fstrsolid%max_ncon)
    real(kind=kreal) :: shell_ref_director(3,fstrsolid%max_ncon)
 
    ! ----- section / material context
    isect = hecmesh%section_ID(icel)
    ihead = hecmesh%section%sect_R_index(isect-1)
    cdsys_id = hecmesh%section%sect_orien_ID(isect)
    sec_opt = hecmesh%section%sect_opt(isect)
    coords = 0.0d0
    if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
 
    material => fstrsolid%elements(icel)%gausses(1)%pMaterial
    thick = hecmesh%section%sect_R_item(ihead+1)
 
    if( ic_type==241 .or. ic_type==242 .or. ic_type==231 .or. ic_type==232 .or. ic_type==2322) then
      if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
      call stf_c2( ic_type, nn, ecoord(1:2,1:nn), fstrsolid%elements(icel)%gausses(:), thick, &
        stiff_mat(1:nn*ndof,1:nn*ndof), fstrsolid%elements(icel)%iset, u(1:2,1:nn) )
 
      if( is_dynamic ) call mass_c2(ic_type, nn, ecoord(1:2,1:nn), fstrsolid%elements(icel)%gausses, &
        sec_opt, thick, mass_mat, lumped)
 
    elseif( ic_type==301 ) then
      call stf_c1( ic_type, nn, ecoord(:,1:nn), thick, fstrsolid%elements(icel)%gausses(:), &
        stiff_mat(1:nn*ndof,1:nn*ndof), u(1:3,1:nn) )
 
    elseif( ic_type==361 ) then
      if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
        call stf_c3( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
          stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
      else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then ! B-bar element
        call stf_c3d8bbar( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
          stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
      else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then ! incompatible element
        call stf_c3d8ic( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
          stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), &
          fstrsolid%elements(icel)%aux, tt(1:nn) )
      else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then ! F-bar element
        call stf_c3d8fbar( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
          stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
      endif
 
      if( is_dynamic ) call mass_c3(ic_type, nn, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses, mass_mat, lumped)
 
    elseif( ic_type==341 .or. ic_type==351 .or. ic_type==342 .or. ic_type==352 .or. ic_type==362 ) then
      ! SESNS option: stiffness is assembled separately via ic_type==881/891;
      ! return stiff_mat = 0 here and still compute the element mass below.
      if( .not. (ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 == kel341sesns) ) then
        call stf_c3( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
          stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
      endif
 
      if( is_dynamic ) call mass_c3(ic_type, nn, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses, mass_mat, lumped)
 
    else if( ic_type == 511 ) then
      call stf_connector( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
        stiff_mat(1:nn*ndof,1:nn*ndof), u(1:3,1:nn), tt(1:nn))
 
      if( is_dynamic ) call dmp_connector( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
        damp_mat(1:nn*ndof,1:nn*ndof), u(1:3,1:nn), tt(1:nn))
 
    else if( ic_type == 611 ) then
      if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
      call stf_beam(ic_type, nn, ecoord, hecmesh%section%sect_R_item(ihead+1:), &
        &   material%variables(m_youngs), material%variables(m_poisson), stiff_mat(1:nn*ndof,1:nn*ndof))
 
      if( is_dynamic ) then
        surf = hecmesh%section%sect_R_item(ihead+4)
        length = get_length(ecoord(1:3,1:nn))
        rho = material%variables(m_density)
        call mass_beam(surf, length, rho, mass_mat)
      endif
 
    else if( ic_type == 641 ) then
      if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
      call stf_beam_641(ic_type, nn, ecoord, fstrsolid%elements(icel)%gausses(:), &
        &            hecmesh%section%sect_R_item(ihead+1:), stiff_mat(1:nn*ndof,1:nn*ndof))
 
      if( is_dynamic ) then
        surf = hecmesh%section%sect_R_item(ihead+4)
        length = get_length(ecoord(1:3,1:nn))
        rho = material%variables(m_density)
        call mass_beam_33(surf, length, rho, mass_mat)
      endif
 
    else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
      if( fstr_uses_finite_rotation_kinematics( ic_type, nn, material ) ) then
        shell_du(:, :) = 0.0d0
        do j = 1, nn
          do i = 1, min(ndof, 6)
            shell_du(i,j) = fstrsolid%dunode(ndof*(nodlocal(j)-1)+i)
          enddo
        enddo
        call shellcomposenodaldisplacement( ndof, nn, u_prev(1:6,1:nn), shell_du(1:6,1:nn), &
          shell_nddisp(1:6,1:nn) )
        call fstr_get_shell_trial_directors( fstrsolid, thick, nn, nodlocal(1:nn), shell_director(1:3,1:nn) )
        call fstr_get_shell_reference_directors( fstrsolid, thick, nn, nodlocal(1:nn), &
          shell_ref_director(1:3,1:nn) )
        call stf_shell_mitc(ic_type, nn, ndof, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses(:), &
          &              stiff_mat(1:nn*ndof,1:nn*ndof), thick, 0, nddisp=shell_nddisp(1:6,1:nn), &
          &              element=fstrsolid%elements(icel), nddirector=shell_director(1:3,1:nn), &
          &              ndrefdirector=shell_ref_director(1:3,1:nn))
      else
        if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
        call stf_shell_mitc(ic_type, nn, ndof, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses(:), &
          &              stiff_mat(1:nn*ndof,1:nn*ndof), thick, 0)
      endif
 
      if( is_dynamic ) then
        rho = material%variables(m_density)
        call mass_shell(ic_type, nn, ecoord(1:3,1:nn), rho, thick, fstrsolid%elements(icel)%gausses, mass_mat, lumped)
      endif
 
    else if( ic_type == 761 ) then   ! for shell-solid mixed analysis
      if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
      call stf_shell_mitc(731, 3, 6, ecoord(1:3,1:3), fstrsolid%elements(icel)%gausses(:), &
        &              stiff_mat(1:nn*ndof,1:nn*ndof), thick, 2)
 
      if( is_dynamic ) then
        surf = get_face3(ecoord(1:3,1:nn))
        rho = material%variables(m_density)
        call mass_s3(surf, thick, rho, mass_mat)
      endif
 
    else if( ic_type == 781 ) then   ! for shell-solid mixed analysis
      if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
      call stf_shell_mitc(741, 4, 6, ecoord(1:3,1:4), fstrsolid%elements(icel)%gausses(:), &
        &              stiff_mat(1:nn*ndof,1:nn*ndof), thick, 1)
 
      if( is_dynamic ) then
        surf = get_face4(ecoord(1:3,1:nn))
        rho = material%variables(m_density)
        call mass_s4(surf, thick, rho, mass_mat)
      endif
 
    elseif( ic_type==3414 ) then
      if( material%mtype /= incomp_newtonian ) call createmat_abort( ic_type, 3, material%mtype )
      call stf_c3_vp( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
        stiff_mat(1:nn*ndof,1:nn*ndof), tincr, u_prev(1:4,1:nn) )
 
    else if( ic_type == 881 .or. ic_type == 891 ) then  ! for selective es/ns smoothed fem
      call stf_c3d4_sesns( ic_type, nn, nodlocal, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
        stiff_mat, cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
 
    else
      call createmat_abort( ic_type, 1 )
    endif
 
  end subroutine calc_stiff_and_mass_elem
 
  !---------------------------------------------------------------------*
  subroutine calc_damping_mat_and_force_elem( ic_type, nn, ndof, material, fstrDYNAMIC, &
      a3, b3, vecA, vecB, stiff_mat, mass_mat, damp_mat, mat, df )
  !---------------------------------------------------------------------*
    integer(kind=kint), intent(in) :: ic_type, nn, ndof
    type(tmaterial), pointer, intent(in)    :: material
    type(fstr_dynamic),       intent(in)    :: fstrDYNAMIC
    real(kind=kreal),         intent(in)    :: a3, b3
    real(kind=kreal),         intent(in)    :: veca(:), vecb(:)
    real(kind=kreal),         intent(in)    :: stiff_mat(:,:), mass_mat(:,:), damp_mat(:,:)
    real(kind=kreal),         intent(out)   :: mat(:,:), df(:)
 
    integer(kind=kint) :: i, j
    real(kind=kreal)   :: ray_m, ray_k, c1, c2
    real(kind=kreal)   :: vecc(20*6), kb(20*6)
 
    if( material%is_elem_Rayleigh_damping ) then
      ray_m = material%variables(m_damping_rm)
      ray_k = material%variables(m_damping_rk)
    else
      ray_m = fstrdynamic%ray_m
      ray_k = fstrdynamic%ray_k
    endif
 
    ! ----- LHS: A = c1*K + c2*M (+ b3*C for connector)
    c1 = 1.d0 + ray_k*b3
    c2 = a3   + ray_m*b3
    do i = 1, nn*ndof
      do j = 1, nn*ndof
        mat(j,i) = c1*stiff_mat(j,i) + c2*mass_mat(j,i)
      enddo
    enddo
    if( ic_type == 511 ) then
      do i = 1, nn*ndof
        do j = 1, nn*ndof
          mat(j,i) = mat(j,i) + b3*damp_mat(j,i)
        enddo
      enddo
    endif
 
    ! ----- RHS: df = M*(vecA + ray_m*vecB) + ray_k*K*vecB (+ C*vecB for connector)
    do i = 1, nn*ndof
      vecc(i) = veca(i) + ray_m*vecb(i)
    enddo
 
    kb = 0.0d0
    do i = 1, nn*ndof
      do j = 1, nn*ndof
        kb(i) = kb(i) + stiff_mat(i,j)*vecb(j)
      enddo
    enddo
 
    df = 0.0d0
    df(1:nn*ndof) = matmul(mass_mat(1:nn*ndof,1:nn*ndof), vecc(1:nn*ndof))
    do i = 1, nn*ndof
      df(i) = df(i) + ray_k*kb(i)
    enddo
    if( ic_type == 511 ) then
      do i = 1, nn*ndof
        do j = 1, nn*ndof
          df(i) = df(i) + damp_mat(i,j)*vecb(j)
        enddo
      enddo
    endif
 
  end subroutine calc_damping_mat_and_force_elem
 
  subroutine createmat_abort( ic_type, flag, mtype )
    integer(kind=kint), intent(in)           :: ic_type
    integer(kind=kint), intent(in)           :: flag
    integer(kind=kint), intent(in), optional :: mtype
 
    if( flag == 1 ) then
      write(*,*) '###ERROR### : Element type not supported for static analysis'
    else if( flag == 2 ) then
      write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
    else if( flag == 3 ) then
      write(*,*) '###ERROR### : This element is not supported for this material'
    endif
    write(*,*) ' ic_type = ', ic_type
    if( present(mtype) ) write(*,*) ' mtype = ', mtype
    call hecmw_abort(hecmw_comm_get_comm())
  end subroutine createmat_abort
 
end module m_fstr_creatematrix_and_dampingforce