fstr_AddBC.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_addbc
  implicit none
contains
 
  !------------------------------------------------------------------------------------------*
  subroutine fstr_addbc(cstep,hecMESH,hecMAT,fstrSOLID,fstrPARAM,hecLagMAT,iter,conMAT,RHSvector)
    !------------------------------------------------------------------------------------------*
    use m_fstr
    use mcontact
    use m_static_lib_1d
    use m_utilities
    use m_fstr_timeinc
    integer, intent(in)                  :: cstep
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: hecMAT
    type(fstr_solid)                     :: fstrSOLID
    type(fstr_param)                     :: fstrPARAM
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    integer(kind=kint)                   :: iter
    type(hecmwst_matrix), optional       :: conMAT
    real(kind=kreal), optional           :: rhsvector(:) 
 
    integer(kind=kint) :: ig0, ig, ityp, idofS, idofE, idof, iS0, iE0, ik, in
    real(kind=kreal)   :: rhs0, rhs, factor, factor0
    integer(kind=kint) :: ndof, grpid, istot
 
    !for rotation
    integer(kind=kint) :: n_rot, rid, jj_n_amp
    type(trotinfo)     :: rinfo
    real(kind=kreal)   :: 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)
 
    !   ----- Prescibed displacement Boundary Conditions
    do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( iter>1 ) then
        factor=0.d0
      else
        jj_n_amp = fstrsolid%BOUNDARY_ngrp_amp(ig0)
        if( jj_n_amp <= 0 ) then  ! Amplitude not defined
          factor0 = fstrsolid%FACTOR(1)
          factor = fstrsolid%FACTOR(2)
        else
          call table_amp(hecmesh,fstrsolid,cstep,jj_n_amp,fstr_get_time(),factor0)
          call table_amp(hecmesh,fstrsolid,cstep,jj_n_amp,fstr_get_time()+fstr_get_timeinc(),factor)
        endif
        factor = factor - factor0
        if(fstrsolid%step_ctrl(cstep)%solution==stepvisco)then
          factor = 0.d0
          if(factor0 < 1.d-10) factor = 1.d0
        endif
      endif
 
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
      ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
      rhs0 = fstrsolid%BOUNDARY_ngrp_val(ig0)
      !
      ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
      idofs = ityp/10
      idofe = ityp - idofs*10
      !
      istot = fstrsolid%BOUNDARY_ngrp_istot(ig0)
      !
      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
        if( istot == 0 ) then
          rhs= rhs0*factor
        else
          ! not implemented yet
          write(*,*) 'Error: rotational boundary cannot be specified with total value'
          call hecmw_abort( hecmw_comm_get_comm() )
        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( istot == 0 ) then
            rhs = rhs0*factor
          else
            rhs = (rhs0 - fstrsolid%unode_bak(ndof*(in-1)+idof))*factor
          endif
          if(present(rhsvector)) then
            rhsvector(ndof*(in-1)+idof) = rhs
            ! write(6,*) 'BC: ', ndof*(in-1)+idof, RHS
            cycle
          endif
          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
 
    !Apply rotational boundary condition
    !need to fix!!
    do rid = 1, n_rot
      if( .not. rinfo%conds(rid)%active ) cycle
      cdiff = 0.d0
      cdiff0 = 0.d0
      cddisp = 0.d0
 
      if( 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( 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(rhsvector)) then
            rhsvector(ndof*(in-1)+idof) = rhs
            ! write(6,*) 'BC(rot): ', ndof*(in-1)+idof, RHS
            cycle
          endif
          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
    if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
 
    !
    !   ------ Truss element Diagonal Modification
    call truss_diag_modify(hecmat,hecmesh)
 
    !
    !   ------ Equation boundary conditions
    do ig0=1,fstrsolid%n_fix_mpc
      if( fstrsolid%mpc_const(ig0) == 0.d0 ) cycle
      ! we need to confirm if it is active in curr step here
      rhs = fstrsolid%mpc_const(ig0)*factor
      hecmesh%mpc%mpc_const(ig0) = rhs
    enddo
 
  end subroutine fstr_addbc
 
end module m_fstr_addbc