dynamic_mat_ass_bc.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_mat_ass_bc
 
contains
 
 
  subroutine dynamic_mat_ass_bc(cstep, hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, fstrPARAM, hecLagMAT, t_curr, iter, conMAT)
    use m_fstr
    use m_table_dyn
    use mcontact
    use m_utilities
 
    implicit none
    integer, intent(in)                  :: cstep
    type(hecmwst_matrix)                 :: hecMAT
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(fstr_dynamic)                   :: fstrDYNAMIC
    type(fstr_param)                     :: fstrPARAM
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    real(kind=kreal)                     :: t_curr
    integer, optional                    :: iter
    type(hecmwst_matrix), optional       :: conMAT
 
    integer(kind=kint) :: ig0, ig, ityp, NDOF, iS0, iE0, ik, in, idofS, idofE, idof
 
    integer(kind=kint) :: flag_u, grpid
    real(kind=kreal)   :: rhs, f_t, f_t1
 
    !for rotation
    integer(kind=kint) :: n_rot, rid, n_nodes
    type(trotinfo)     :: rinfo
    real(kind=kreal)   :: theta, normal(3), direc(3), ccoord(3), cdiff(3), cdiff0(3)
    real(kind=kreal)   :: cdisp(3), cddisp(3)
    real(kind=kreal)   :: rotation_factor
    !
    ndof = hecmat%NDOF
    n_rot = fstrsolid%BOUNDARY_ngrp_rot
    if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
 
    flag_u = 1
    !C=============================C
    !C-- implicit dynamic analysis
    !C=============================C
    if( fstrdynamic%idx_eqa == 1 ) then
 
      do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
        ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
        grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
        if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
        rhs  = fstrsolid%BOUNDARY_ngrp_val(ig0)
 
        if( present(iter) ) then
          if( iter>1 ) then
            rhs=0.d0
          else
            f_t1 = rhs
            call fstr_get_amplitude_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr-fstrdynamic%t_delta, f_t1, flag_u)
            f_t = rhs
            call fstr_get_amplitude_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, f_t, flag_u)
            rhs = f_t - f_t1
          endif
        else
          call fstr_get_amplitude_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, rhs, flag_u)
        endif
 
        ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
        idofs = ityp/10
        idofe = ityp - idofs*10
 
        is0 = hecmesh%node_group%grp_index(ig-1) + 1
        ie0 = hecmesh%node_group%grp_index(ig  )
 
        if( fstrsolid%BOUNDARY_ngrp_rotID(ig0) > 0 ) then ! setup rotation information
          rid = fstrsolid%BOUNDARY_ngrp_rotID(ig0)
          if( .not. rinfo%conds(rid)%active ) then
            rinfo%conds(rid)%active = .true.
            rinfo%conds(rid)%center_ngrp_id = fstrsolid%BOUNDARY_ngrp_centerID(ig0)
            rinfo%conds(rid)%torque_ngrp_id = ig
          endif
          do idof=idofs,idofe
            if( idof>ndof ) then
              rinfo%conds(rid)%vec(idof-ndof) = rhs
            else
              rinfo%conds(rid)%vec(idof) = rhs
            endif
          enddo
          cycle
        endif
 
        do ik = is0, ie0
          in = hecmesh%node_group%grp_item(ik)
 
          do idof = idofs, idofe
            if(present(conmat)) then
              call hecmw_mat_ass_bc(hecmat, in, idof, rhs, conmat)
            else
              call hecmw_mat_ass_bc(hecmat, in, idof, rhs)
            endif
            if( fstr_is_contact_active() .and. fstrparam%contact_algo == kcaslagrange  &
                .and. fstrparam%nlgeom .and. fstrdynamic%idx_resp == 1 )  then
              if(present(conmat)) then
                call hecmw_mat_ass_bc_contactlag(conmat,heclagmat,in,idof,rhs)
              else
                call hecmw_mat_ass_bc_contactlag(hecmat,heclagmat,in,idof,rhs)
              endif
            endif
          enddo
        enddo
 
      enddo
 
      !Apply rotational boundary condition
      do rid = 1, n_rot
        if( .not. rinfo%conds(rid)%active ) cycle
        cdiff = 0.d0
        cdiff0 = 0.d0
        cddisp = 0.d0
 
        if( present(iter) ) then
          if( iter > 1 ) then
            rotation_factor = 0.d0  ! No additional rotation for subsequent iterations
          else
            rotation_factor = 1.0d0  ! Use the rotation vector as-is (already contains increment)
          endif
        else
          rotation_factor = 1.0d0  ! Use the rotation vector as-is
        endif
 
        if( rotation_factor > 0.d0 ) then
          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)
            cddisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmat%B)
          enddo
          ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
        endif
 
        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)
          if( rotation_factor > 0.d0 ) then
            cdiff0(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in)+fstrsolid%unode(ndof*(in-1)+1:ndof*in)-ccoord(1:ndof)
            cdiff(1:ndof) = cdiff0(1:ndof)
            call rotate_3dvector_by_rodrigues_formula(rinfo%conds(rid)%vec(1:ndof),cdiff(1:ndof))
          endif
          do idof = 1, ndof
            rhs = cdiff(idof)-cdiff0(idof)+cddisp(idof)
            if(present(conmat)) then
              call hecmw_mat_ass_bc(hecmat, in, idof, rhs, conmat)
            else
              call hecmw_mat_ass_bc(hecmat, in, idof, rhs)
            endif
            if( fstr_is_contact_active() .and. fstrparam%solution_type == kststatic   &
                .and. fstrparam%contact_algo == kcaslagrange ) then
              if(present(conmat)) then
                call hecmw_mat_ass_bc_contactlag(conmat,heclagmat,in,idof,rhs)
              else
                call hecmw_mat_ass_bc_contactlag(hecmat,heclagmat,in,idof,rhs)
              endif
            endif
          enddo
        enddo
      enddo
      !C
      !C-- end of implicit dynamic analysis
      !C
 
      !C=============================C
      !C-- explicit dynamic analysis
      !C=============================C
    else if( fstrdynamic%idx_eqa == 11 ) then
      !C
      ndof = hecmat%NDOF
      do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
        ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
        rhs  = fstrsolid%BOUNDARY_ngrp_val(ig0)
 
        call fstr_get_amplitude_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, rhs, flag_u)
 
        ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
 
        is0 = hecmesh%node_group%grp_index(ig-1) + 1
        ie0 = hecmesh%node_group%grp_index(ig  )
        idofs = ityp/10
        idofe = ityp - idofs*10
 
        do ik = is0, ie0
          in = hecmesh%node_group%grp_item(ik)
 
          do idof = idofs, idofe
            hecmat%B        (ndof*in-(ndof-idof)) = rhs
            fstrdynamic%VEC1(ndof*in-(ndof-idof)) = 1.0d0
          end do
        enddo
      enddo
      !C
      !C-- end of explicit dynamic analysis
      !C
    end if
 
    if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
 
  end subroutine dynamic_mat_ass_bc
 
 
  !C***
  !C***
  subroutine dynamic_bc_init(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, t_curr)
    use m_fstr
    use m_table_dyn
 
    implicit none
    type(hecmwst_matrix)     :: hecmat
    type(hecmwst_local_mesh) :: hecMESH
    type(fstr_solid)         :: fstrSOLID
    type(fstr_dynamic)       :: fstrDYNAMIC
    real(kind=kreal)         :: t_curr
 
    integer(kind=kint) :: NDOF, ig0, ig, ityp, iS0, iE0, ik, in, idofS, idofE, idof
    integer(kind=kint) :: flag_u, grpid
    real(kind=kreal)   :: rhs, f_t
 
    flag_u = 1
    ndof = hecmat%NDOF
 
    do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
      ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
      rhs  = fstrsolid%BOUNDARY_ngrp_val(ig0)
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, 1 ) ) cycle
 
      call fstr_get_amplitude_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, rhs, flag_u)
 
      ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
 
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig  )
      idofs = ityp/10
      idofe = ityp - idofs*10
 
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
 
        do idof = idofs, idofe
          fstrdynamic%DISP(ndof*in-(ndof-idof),1) = rhs
        end do
      enddo
    enddo
 
    return
  end subroutine dynamic_bc_init
 
  subroutine dynamic_explicit_ass_bc(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, t_curr, iter)
    use m_fstr
    use m_table_dyn
    use mcontact
    use m_utilities
 
    implicit none
    type(hecmwst_matrix)                 :: hecmat
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(fstr_dynamic)                   :: fstrDYNAMIC
    real(kind=kreal)                     :: t_curr
    integer, optional                    :: iter
 
    integer(kind=kint) :: ig0, ig, ityp, NDOF, iS0, iE0, ik, in, idofS, idofE, idof
 
    integer(kind=kint) :: flag_u
    real(kind=kreal)   :: rhs, f_t, f_t1
 
    !for rotation
    integer(kind=kint) :: n_rot, rid, n_nodes
    type(trotinfo)     :: rinfo
    real(kind=kreal)   :: theta, normal(3), direc(3), ccoord(3), cdiff(3), cdiff0(3)
    real(kind=kreal)   :: cdisp(3), cddisp(3)
    !
    ndof = hecmat%NDOF
    n_rot = fstrsolid%BOUNDARY_ngrp_rot
    if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
 
    flag_u = 1
 
      !C
      ndof = hecmat%NDOF
      do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
        ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
        rhs  = fstrsolid%BOUNDARY_ngrp_val(ig0)
 
        call fstr_get_amplitude_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, rhs, flag_u)
 
        ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
 
        is0 = hecmesh%node_group%grp_index(ig-1) + 1
        ie0 = hecmesh%node_group%grp_index(ig  )
        idofs = ityp/10
        idofe = ityp - idofs*10
 
        if( fstrsolid%BOUNDARY_ngrp_rotID(ig0) > 0 ) then ! setup rotation information
          rid = fstrsolid%BOUNDARY_ngrp_rotID(ig0)
          if( .not. rinfo%conds(rid)%active ) then
            rinfo%conds(rid)%active = .true.
            rinfo%conds(rid)%center_ngrp_id = fstrsolid%BOUNDARY_ngrp_centerID(ig0)
            rinfo%conds(rid)%torque_ngrp_id = ig
          endif
          do idof=idofs,idofe
            if( idof>ndof ) then
              rinfo%conds(rid)%vec(idof-ndof) = rhs
            else
              rinfo%conds(rid)%vec(idof) = rhs
            endif
          enddo
          cycle
        endif
 
        do ik = is0, ie0
          in = hecmesh%node_group%grp_item(ik)
 
          do idof = idofs, idofe
            hecmat%B(ndof*in-(ndof-idof)) = rhs*fstrdynamic%VEC1(ndof*in-(ndof-idof))
        !    fstrDYNAMIC%VEC1(NDOF*in-(NDOF-idof)) = 1.0d0
          end do
        enddo
      enddo
 
      !Apply rotational boundary condition (prescribed absolute displacement)
      !  The explicit scheme stores the absolute prescribed displacement in hecMAT%B
      !  (B = u*VEC1, so that X = B/VEC1 = u). The rotation factor in rinfo%vec is the
      !  total rotation at t_curr (val*f_t), hence the original offset r0 = x0 - c is
      !  rotated and u = (R*r0 - r0) + c_disp.
      do rid = 1, n_rot
        if( .not. rinfo%conds(rid)%active ) cycle
        cdiff = 0.d0
        cdiff0 = 0.d0
        cdisp = 0.d0
 
        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)
        enddo
 
        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)
          cdiff0(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in) - ccoord(1:ndof)
          cdiff(1:ndof) = cdiff0(1:ndof)
          call rotate_3dvector_by_rodrigues_formula(rinfo%conds(rid)%vec(1:ndof),cdiff(1:ndof))
          do idof = 1, ndof
            rhs = cdiff(idof) - cdiff0(idof) + cdisp(idof)
            hecmat%B(ndof*in-(ndof-idof)) = rhs*fstrdynamic%VEC1(ndof*in-(ndof-idof))
          enddo
        enddo
      enddo
 
    if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
 
  end subroutine dynamic_explicit_ass_bc
 
end module m_dynamic_mat_ass_bc