fstr_ass_load.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_ass_load
  implicit none
contains
  !
  !======================================================================!
 
  subroutine fstr_ass_load(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
    !======================================================================!
    use m_fstr
    use m_static_lib
    use m_fstr_precheck
    use mmechgauss
    use mreadtemp
    use muload
    use m_fstr_spring
    use m_common_struct
    use m_utilities
 
    integer(kind=kint), intent(in)       :: cstep
    real(kind=kreal), intent(in)         :: ctime 
    type(hecmwst_matrix), intent(inout)  :: hecmat
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)      :: fstrSOLID
    type(fstr_param), intent(inout)      :: fstrPARAM
 
    real(kind=kreal)   :: xx(20), yy(20), zz(20)
    real(kind=kreal)   :: params(0:6)
    real(kind=kreal)   :: vect(60)
    integer(kind=kint) :: iwk(60)
    integer(kind=kint) :: nodLocal(20)
    real(kind=kreal)   :: tt(20), tt0(20), coords(3, 3), factor
    integer(kind=kint) :: ndof, ig0, ig, ityp, ltype, iS0, iE0, ik, in, i, j
    integer(kind=kint) :: icel, ic_type, nn, is, isect, id, iset, nsize
    integer(kind=kint) :: itype, iE, ierror, grpid, cdsys_ID, jj_n_amp
    real(kind=kreal)   :: fval, rho, thick, pa1
    logical :: fg_surf
    integer(kind=kint) :: tstep
    type(tmaterial), pointer :: material
    integer(kind=kint) :: ihead
    real(kind=kreal)   :: a
 
    !for torque load
    integer(kind=kint) :: n_rot, rid, n_nodes, idof
    type(trotinfo)     :: rinfo
    real(kind=kreal)   :: tval, normal(3), direc(3), ccoord(3), cdisp(3), cdiff(3)
 
    ndof = hecmat%NDOF
 
    ! -------------------------------------------------------------------
    !  CLOAD
    ! -------------------------------------------------------------------
    n_rot = fstrsolid%CLOAD_ngrp_rot
    if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
 
    fstrsolid%GL(:) = 0.0d0
    fstrsolid%EFORCE(:) = 0.0d0
    do ig0 = 1, fstrsolid%CLOAD_ngrp_tot
      grpid = fstrsolid%CLOAD_ngrp_GRPID(ig0)
      if( .not. fstr_isloadactive( fstrsolid, grpid, cstep ) ) cycle
      jj_n_amp = fstrsolid%CLOAD_ngrp_amp(ig0)
      if( jj_n_amp <= 0 ) then  ! Amplitude not defined
        factor = fstrsolid%FACTOR(2)
      else
        call table_amp(hecmesh,fstrsolid,cstep,jj_n_amp,ctime,factor)
      endif
 
      if( fstr_isloadactive( fstrsolid, grpid, cstep-1 ) ) factor = 1.0d0
      ig = fstrsolid%CLOAD_ngrp_ID(ig0)
      ityp = fstrsolid%CLOAD_ngrp_DOF(ig0)
      fval = fstrsolid%CLOAD_ngrp_val(ig0)
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig  )
 
      if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 ) then ! setup torque load information
        rid = fstrsolid%CLOAD_ngrp_rotID(ig0)
        if( .not. rinfo%conds(rid)%active ) then
          rinfo%conds(rid)%active = .true.
          rinfo%conds(rid)%center_ngrp_id = fstrsolid%CLOAD_ngrp_centerID(ig0)
          rinfo%conds(rid)%torque_ngrp_id = ig
        endif
        if( ityp>ndof ) ityp = ityp-ndof
        rinfo%conds(rid)%vec(ityp) = factor*fval
        cycle
      endif
 
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
        fstrsolid%GL(ndof*(in-1)+ityp) = fstrsolid%GL(ndof*(in-1)+ityp)+factor*fval
      enddo
    enddo
 
    !Add torque load to fstrSOLID%GL
    do rid = 1, n_rot
      if( .not. rinfo%conds(rid)%active ) cycle
      !get number of slave nodes
      n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
 
      !get center node
      ig = rinfo%conds(rid)%center_ngrp_id
      do idof = 1, ndof
        ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
        cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
        cdisp(idof) = cdisp(idof) + hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%dunode)
      enddo
      ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
 
      tval = dsqrt(dot_product(rinfo%conds(rid)%vec(1:ndof),rinfo%conds(rid)%vec(1:ndof)))
      if( tval < 1.d-16 ) then
        write(*,*) '###ERROR### : norm of torque vector must be > 0.0'
        call hecmw_abort( hecmw_comm_get_comm() )
      endif
      normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof)/tval
      tval = tval/dble(n_nodes)
 
      ig = rinfo%conds(rid)%torque_ngrp_id
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig  )
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
        cdiff(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in)+fstrsolid%unode(ndof*(in-1)+1:ndof*in) &
          &  +fstrsolid%dunode(ndof*(in-1)+1:ndof*in)-ccoord(1:ndof)
        call cross_product(normal,cdiff,vect(1:ndof))
        fval = dot_product(vect(1:ndof),vect(1:ndof))
        if( fval < 1.d-16 ) then
          write(*,*) '###ERROR### : torque node is at the same position as that of center node in rotational surface.'
          call hecmw_abort( hecmw_comm_get_comm() )
        endif
        vect(1:ndof) = (tval/fval)*vect(1:ndof)
        fstrsolid%GL(ndof*(in-1)+1:ndof*in) = fstrsolid%GL(ndof*(in-1)+1:ndof*in)+vect(1:ndof)
      enddo
    enddo
    if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
    !
    ! -------------------------------------------------------------------
    !  DLOAD
    ! -------------------------------------------------------------------
    do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
      grpid = fstrsolid%DLOAD_ngrp_GRPID(ig0)
      if( .not. fstr_isloadactive( fstrsolid, grpid, cstep ) ) cycle
      jj_n_amp = fstrsolid%DLOAD_ngrp_amp(ig0)
      if( jj_n_amp <= 0 ) then  ! Amplitude not defined
        factor = fstrsolid%FACTOR(2)
      else
        call table_amp(hecmesh,fstrsolid,cstep,jj_n_amp,ctime,factor)
      endif
 
      if( fstr_isloadactive( fstrsolid, grpid, cstep-1 ) ) factor = 1.0d0
      ig = fstrsolid%DLOAD_ngrp_ID(ig0)
      ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
      do i = 0, 6
        params(i)= fstrsolid%DLOAD_ngrp_params(i,ig0)
      enddo
      ! ----- START & END
      fg_surf = (ltype == 100)
      if( fg_surf ) then                  ! surface group
        is0 = hecmesh%surf_group%grp_index(ig-1) + 1
        ie0 = hecmesh%surf_group%grp_index(ig  )
      else                                ! element group
        is0 = hecmesh%elem_group%grp_index(ig-1) + 1
        ie0 = hecmesh%elem_group%grp_index(ig  )
      endif
      do ik = is0, ie0
        if( fg_surf ) then                ! surface group
          ltype   = hecmesh%surf_group%grp_item(2*ik)*10
          icel    = hecmesh%surf_group%grp_item(2*ik-1)
          ic_type = hecmesh%elem_type(icel)
        else                              ! element group
          icel    = hecmesh%elem_group%grp_item(ik)
          ic_type = hecmesh%elem_type(icel)
        endif
        if( hecmw_is_etype_link(ic_type) ) cycle
        if( hecmw_is_etype_patch(ic_type) ) cycle
        ! if( ic_type==3422 ) ic_type=342
        nn = hecmw_get_max_node(ic_type)
        ! ----- node ID
        is = hecmesh%elem_node_index(icel-1)
        if( fstrsolid%DLOAD_follow == 0 ) then
          do j = 1, nn
            nodlocal(j) = hecmesh%elem_node_item (is+j)
            ! ----- nodal coordinate
            xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
            yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
            zz(j) = hecmesh%node( 3*nodlocal(j)   )
            ! ----- create iwk array ***
            do i = 1, ndof
              iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
            enddo
          enddo
        else
          do j = 1, nn
            nodlocal(j) = hecmesh%elem_node_item (is+j)
            ! ----- nodal coordinate
            if (ndof==2) then
              xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 2*nodlocal(j)-1 )+fstrsolid%dunode( 2*nodlocal(j)-1 )
              yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 2*nodlocal(j)   )+fstrsolid%dunode( 2*nodlocal(j)   )
            else if (ndof==3) then
              xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 3*nodlocal(j)-2 )+fstrsolid%dunode( 3*nodlocal(j)-2 )
              yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 3*nodlocal(j)-1 )+fstrsolid%dunode( 3*nodlocal(j)-1 )
              zz(j) = hecmesh%node( 3*nodlocal(j)   )+fstrsolid%unode( 3*nodlocal(j)   )+fstrsolid%dunode( 3*nodlocal(j)   )
            else if (ndof==6) then
              xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 6*nodlocal(j)-5 )+fstrsolid%dunode( 6*nodlocal(j)-5 )
              yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 6*nodlocal(j)-4 )+fstrsolid%dunode( 6*nodlocal(j)-4 )
              zz(j) = hecmesh%node( 3*nodlocal(j)   )+fstrsolid%unode( 6*nodlocal(j)-3 )+fstrsolid%dunode( 6*nodlocal(j)-3 )
            endif
            ! ----- create iwk array ***
            do i = 1, ndof
              iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
            enddo
          enddo
        end if
        ! ----- section  ID
        isect = hecmesh%section_ID(icel)
        ! ----- Get Properties
        material => fstrsolid%elements(icel)%gausses(1)%pMaterial
        rho = material%variables(m_density)
        call fstr_get_thickness(hecmesh,isect,thick)
        ! ----- Section Data
        if( ndof == 2 ) then
          id=hecmesh%section%sect_opt(isect)
          if( id == 0 ) then
            iset = 1
          else if( id == 1 ) then
            iset = 0
          else if( id == 2) then
            iset = 2
          endif
          pa1=1.d0
        endif
        ! ----- Create local stiffness
        if (ic_type==301)then
          ihead = hecmesh%section%sect_R_index(isect-1)
          call dl_c1(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,thick,ltype,params,vect(1:nn*ndof),nsize)
 
        elseif( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
          call dl_c2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset)
 
        else if ( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .or.   &
            ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
          call dl_c3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize)
 
        else if ( ic_type == 641 ) then
          ihead = hecmesh%section%sect_R_index(isect-1)
          call dl_beam_641(ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), rho, ltype, params, &
            hecmesh%section%sect_R_item(ihead+1:), vect(1:nn*ndof), nsize)
 
        else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
          call dl_shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrsolid%elements(icel)%gausses)
 
        else if( ( ic_type==761 ) .or. ( ic_type==781 ) ) then
          call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
            fstrsolid%elements(icel)%gausses)
 
        else
          nsize = 0
          write(*,*)"### WARNING: DLOAD",ic_type
 
        endif
        ! ----- Add vector
        do j=1,nsize
          fstrsolid%GL( iwk(j) )=fstrsolid%GL( iwk(j) )+factor*vect(j)
        enddo
      enddo
    enddo
 
    ! -----Uload
    call uloading( cstep, factor, fstrsolid%GL )
 
    !C
    !C Update for fstrSOLID%GL
    !C
    call hecmw_update_r (hecmesh,fstrsolid%GL,hecmesh%n_node, hecmesh%n_dof)
 
    call hecmw_mat_clear_b( hecmat )
    do i=1, hecmesh%n_node*  hecmesh%n_dof
      hecmat%B(i)=fstrsolid%GL(i)-fstrsolid%QFORCE(i)
    enddo
 
    do i=1, hecmat%NDOF*hecmat%NP
      !thermal load is not considered
      fstrsolid%EFORCE(i) = fstrsolid%GL(i)
    enddo
 
    ! -------------------------------------------------------------------
    !  TLOAD : THERMAL LOAD USING TEMPERATURE
    ! -------------------------------------------------------------------
    !C
    !C Set Temperature
    !C
    if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
      do ig0 = 1, fstrsolid%TEMP_ngrp_tot
        grpid = fstrsolid%TEMP_ngrp_GRPID(ig0)
        if( .not. fstr_isloadactive( fstrsolid, grpid, cstep ) ) cycle
        factor = fstrsolid%factor(2)
        if( fstr_isloadactive( fstrsolid, grpid, cstep-1 ) ) factor = 1.0d0
        ig = fstrsolid%TEMP_ngrp_ID(ig0)
        fval =fstrsolid%TEMP_ngrp_val(ig0)
        is0 = hecmesh%node_group%grp_index(ig-1)+1
        ie0 = hecmesh%node_group%grp_index(ig  )
        do ik = is0, ie0
          in   = hecmesh%node_group%grp_item(ik)
          pa1 = fstrsolid%temp_bak( in )
          fstrsolid%temperature( in ) = pa1+(fval-pa1)*factor
        enddo
      enddo
 
      if( fstrsolid%TEMP_irres > 0 ) then
        call read_temperature_result(hecmesh, fstrsolid%TEMP_irres, fstrsolid%TEMP_tstep, &
          &  fstrsolid%TEMP_rtype, fstrsolid%TEMP_interval, fstrsolid%TEMP_factor, ctime, &
          &  fstrsolid%temperature, fstrsolid%temp_bak)
      endif
 
      ! ----- element TYPE loop.
      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)
        if( hecmw_is_etype_link(ic_type) ) cycle
        if( hecmw_is_etype_patch(ic_type) ) cycle
        ! ----- Set number of nodes
        nn = hecmw_get_max_node(ic_type)
 
        ! ----- element loop
        do icel = is, ie
 
          ! ----- node ID
          is= hecmesh%elem_node_index(icel-1)
          do j=1,nn
            nodlocal(j)=hecmesh%elem_node_item(is+j)
            ! ----- nodal coordinate
            if (ndof==2) then
              xx(j)=hecmesh%node(3*nodlocal(j)-2)+fstrsolid%unode(ndof*nodlocal(j)-1)
              yy(j)=hecmesh%node(3*nodlocal(j)-1)+fstrsolid%unode(ndof*nodlocal(j)  )
            else if (ndof==3) then
              xx(j)=hecmesh%node(3*nodlocal(j)-2)+fstrsolid%unode(ndof*nodlocal(j)-2)
              yy(j)=hecmesh%node(3*nodlocal(j)-1)+fstrsolid%unode(ndof*nodlocal(j)-1)
              zz(j)=hecmesh%node(3*nodlocal(j)  )+fstrsolid%unode(ndof*nodlocal(j))
            endif
            tt0(j)=fstrsolid%last_temp( nodlocal(j) )
            tt(j) = fstrsolid%temperature( nodlocal(j) )
            ! ----- create iwk array ***
            do i=1,ndof
              iwk(ndof*(j-1)+i)=ndof*(nodlocal(j)-1)+i
            enddo
          enddo
 
          ! ----- section  Data
          isect= hecmesh%section_ID(icel)
          cdsys_id = hecmesh%section%sect_orien_ID(isect)
          call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
 
          if( ndof == 2 ) then
            id=hecmesh%section%sect_opt(isect)
            if( id==0 ) then
              iset = 1
            else if( id == 1 ) then
              iset = 0
            else if( id == 2 ) then
              iset = 2
            endif
            pa1 = 1.0d0
          endif
 
          if( ic_type == 641 ) then
 
            isect= hecmesh%section_ID(icel)
            ihead = hecmesh%section%sect_R_index(isect-1)
 
            call tload_beam_641( ic_type, nn, ndof, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),    &
              fstrsolid%elements(icel)%gausses, hecmesh%section%sect_R_item(ihead+1:), &
              vect(1:nn*ndof) )
 
            do j = 1, ndof*nn
              hecmat%B( iwk(j) ) = hecmat%B( iwk(j) )+vect(j)
            end do
            cycle
          end if
 
          ! ----- Create local stiffness
          if(ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 ) then
            call tload_c2( ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn),      &
              fstrsolid%elements(icel)%gausses,pa1, iset, vect(1:nn*2) )
 
          else if( ic_type == 361 ) then
            if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then
              call tload_c3                                                          &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),       &
                fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
            else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then
              call tload_c3d8bbar                                                          &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),       &
                fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
            else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then
              call tload_c3d8ic                                                            &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),       &
                fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
            else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then
              call tload_c3d8fbar                                                            &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),       &
                fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
            endif
 
          else if( ic_type == 341 .or. ic_type == 351 .or.                       &
              ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
            call tload_c3                                                                &
              ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),       &
              fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
 
          else if( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 ) then
            if( myrank == 0 ) then
              write(imsg,*) '*------------------------', &
                '-------------------*'
              write(imsg,*) ' Thermal loading option for shell elements', &
                'not yet available.'
              write(imsg,*) '*------------------------', &
                '-------------------*'
              call hecmw_abort( hecmw_comm_get_comm())
            endif
 
          endif
 
          ! ----- Add vector
          do j = 1, ndof*nn
            ! fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )+vect(j)
            hecmat%B( iwk(j) ) = hecmat%B( iwk(j) )+vect(j)
          enddo
 
        enddo
      enddo
    endif
 
    ! ----- Spring force
    call fstr_update_ndforce_spring( cstep, hecmesh, fstrsolid, hecmat%B )
 
  end subroutine fstr_ass_load
 
end module m_fstr_ass_load