fstr_StiffMatrix.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_stiffmatrix
  use m_fstr
  implicit none
 
  private
  public :: fstr_stiffmatrix
 
contains
 
  !---------------------------------------------------------------------*
  subroutine fstr_stiffmatrix( hecMESH, hecMAT, fstrSOLID, time, tincr)
    !---------------------------------------------------------------------*
    use m_static_lib
    use mmechgauss
 
    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( tmaterial ), pointer :: material
 
    real(kind=kreal)   :: stiffness(fstrsolid%max_ncon_stf*6, fstrsolid%max_ncon_stf*6)
    integer(kind=kint) :: nodlocal(fstrsolid%max_ncon)
    real(kind=kreal)   :: tt(fstrsolid%max_ncon), ecoord(3,fstrsolid%max_ncon)
    real(kind=kreal)   :: thick
    integer(kind=kint) :: ndof, itype, is, ie, ic_type, nn, icel, iis, i, j
    real(kind=kreal)   :: u(6,fstrsolid%max_ncon), du(6,fstrsolid%max_ncon), coords(3,3), u_prev(6,fstrsolid%max_ncon)
    integer            :: isect, ihead, cdsys_id
 
    ! ----- initialize
    call hecmw_mat_clear( hecmat )
 
    ndof = hecmat%NDOF
    tt(:) = 0.d0
 
    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
      ! ----- Set number of nodes
 
      ! ----- element loop
      !$omp parallel default(none), &
        !$omp&  private(icel,iiS,nn,j,nodLOCAL,i,ecoord,du,u,u_prev,tt,cdsys_ID,coords, &
        !$omp&          material,thick,stiffness,isect,ihead), &
        !$omp&  shared(iS,iE,hecMESH,ndof,fstrSOLID,ic_type,hecMAT,time,tincr)
      !$omp do
      do icel= is, ie
 
        ! ----- nodal coordinate & displacement
        iis= hecmesh%elem_node_index(icel-1)
        nn = hecmesh%elem_node_index(icel)-iis
     !   if( nn>150 ) stop "elemental nodes > 150!"
 
        do j=1,nn
          nodlocal(j)= hecmesh%elem_node_item (iis+j)
          do i=1, 3
            ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
          enddo
          do i=1,ndof
            du(i,j) = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
            u(i,j)  = fstrsolid%unode(ndof*nodlocal(j)+i-ndof) + du(i,j)
            u_prev(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
          enddo
          if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 )  &
            tt(j)=fstrsolid%temperature( nodlocal(j) )
        enddo
 
        isect = hecmesh%section_ID(icel)
        ihead = hecmesh%section%sect_R_index(isect-1)
        cdsys_id = hecmesh%section%sect_orien_ID(isect)
        if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
        thick = hecmesh%section%sect_R_item(ihead+1)
        material => fstrsolid%elements(icel)%gausses(1)%pMaterial
 
        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 stiffmat_abort( ic_type, 2 )
          call stf_c2( ic_type,nn,ecoord(1:2,1:nn),fstrsolid%elements(icel)%gausses(:),thick,  &
            stiffness(1:nn*ndof,1:nn*ndof), fstrsolid%elements(icel)%iset,          &
            u(1:2,1:nn) )
 
        elseif ( ic_type==301 ) then
          call stf_c1( ic_type,nn,ecoord(:,1:nn),thick,fstrsolid%elements(icel)%gausses(:),   &
            stiffness(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(:),                &
              stiffness(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(:),                &
              stiffness(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(:), &
              stiffness(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(:),                &
              stiffness(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn), tt(1:nn) )
          endif
 
        elseif (ic_type==341 .or. ic_type==351 .or. ic_type==342 .or. ic_type==352 .or. ic_type==362 ) then
          if( ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 == kel341sesns ) cycle ! skip smoothed fem
          call stf_c3                                                                              &
            ( ic_type, nn, ecoord(:, 1:nn), fstrsolid%elements(icel)%gausses(:),                &
            stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
 
        else if( ic_type == 511) then
          call stf_connector( ic_type,nn,ecoord(:,1:nn),fstrsolid%elements(icel)%gausses(:),   &
            stiffness(1:nn*ndof,1:nn*ndof), u(1:3,1:nn), tincr, tt(1:nn) )
 
        else if( ic_type == 611) then
          if( material%nlgeom_flag /= infinitesimal ) call stiffmat_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), stiffness(1:nn*ndof,1:nn*ndof))
 
        else if( ic_type == 641 ) then
          if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
          call stf_beam_641(ic_type, nn, ecoord, fstrsolid%elements(icel)%gausses(:), &
            &            hecmesh%section%sect_R_item(ihead+1:), stiffness(1:nn*ndof,1:nn*ndof))
 
        else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
          if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
          call stf_shell_mitc(ic_type, nn, ndof, ecoord(1:3, 1:nn), fstrsolid%elements(icel)%gausses(:), &
            &              stiffness(1:nn*ndof, 1:nn*ndof), thick, 0)
 
        else if( ic_type == 761 ) then   !for shell-solid mixed analysis
          if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
          call stf_shell_mitc(731, 3, 6, ecoord(1:3, 1:3), fstrsolid%elements(icel)%gausses(:), &
            &              stiffness(1:nn*ndof, 1:nn*ndof), thick, 2)
 
        else if( ic_type == 781 ) then   !for shell-solid mixed analysis
          if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
          call stf_shell_mitc(741, 4, 6, ecoord(1:3, 1:4), fstrsolid%elements(icel)%gausses(:), &
            &              stiffness(1:nn*ndof, 1:nn*ndof), thick, 1)
 
        elseif ( ic_type==3414 ) then
          if( material%mtype /= incomp_newtonian) call stiffmat_abort( ic_type, 3, material%mtype )
          call stf_c3_vp                                                           &
            ( ic_type, nn, ecoord(:, 1:nn),fstrsolid%elements(icel)%gausses(:), &
            stiffness(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(:),         &
            stiffness, cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
        else
          call stiffmat_abort( ic_type, 1 )
        endif
 
        ! elemact element
        if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) then
          call stf_dummy( ndof, nn, ecoord(:,1:nn), u(1:3,1:nn), &
            &  stiffness(1:nn*ndof, 1:nn*ndof), fstrsolid%elements(icel) )
          !stiffness(:,:) = fstrSOLID%elements(icel)%elemact_coeff*stiffness(:,:)
        end if
        !
        ! ----- CONSTRUCT the GLOBAL MATRIX STARTED
        call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiffness)
 
      enddo      ! icel
      !$omp end do
      !$omp end parallel
    enddo        ! itype
 
  end subroutine fstr_stiffmatrix
 
  subroutine stiffmat_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
 
end module m_fstr_stiffmatrix