fstr_contact_output.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_contact_output
  use hecmw
  use m_fstr
  use mcontactdef
  use m_fstr_contact_geom
  implicit none
 
  public :: initialize_contact_output_vectors
  public :: initialize_embed_vectors
  public :: setup_contact_elesurf_for_area
  public :: calc_contact_area
  public :: fstr_setup_parancon_contactvalue
  public :: update_contact_state_vectors
 
  private :: calc_nodalarea_surfelement
 
contains
 
  subroutine initialize_contact_output_vectors(fstrSOLID,hecMAT)
    type(fstr_solid)       :: fstrsolid
    type(hecmwst_matrix)   :: hecmat
 
    if( .not. associated(fstrsolid%CONT_NFORCE) ) then
      allocate( fstrsolid%CONT_NFORCE(hecmat%NP*3) )
      fstrsolid%CONT_NFORCE(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_FRIC) ) then
      allocate( fstrsolid%CONT_FRIC(hecmat%NP*3) )
      fstrsolid%CONT_FRIC(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_RELVEL) ) then
      allocate( fstrsolid%CONT_RELVEL(hecmat%NP*3) )
      fstrsolid%CONT_RELVEL(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_STATE) ) then
      allocate( fstrsolid%CONT_STATE(hecmat%NP*1) )
      fstrsolid%CONT_STATE(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_AREA) ) then
      allocate( fstrsolid%CONT_AREA(hecmat%NP) )
      fstrsolid%CONT_AREA(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_NTRAC) ) then
      allocate( fstrsolid%CONT_NTRAC(hecmat%NP*3) )
      fstrsolid%CONT_NTRAC(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_FTRAC) ) then
      allocate( fstrsolid%CONT_FTRAC(hecmat%NP*3) )
      fstrsolid%CONT_FTRAC(:) = 0.d0
    end if
 
  end subroutine
 
  subroutine initialize_embed_vectors(fstrSOLID,hecMAT)
    type(fstr_solid)       :: fstrsolid
    type(hecmwst_matrix)   :: hecmat
 
    if( .not. associated(fstrsolid%EMBED_NFORCE) ) then
      allocate( fstrsolid%EMBED_NFORCE(hecmat%NP*3) )
      fstrsolid%EMBED_NFORCE(:) = 0.d0
    end if
  end subroutine
 
  subroutine setup_contact_elesurf_for_area( hecMESH, fstrSOLID )
    type( hecmwst_local_mesh ), intent(in)       :: hecmesh
    type(fstr_solid), intent(inout)              :: fstrsolid
 
    integer(kind=kint) :: i, j, k, sgrp_id, is, ie, eid, sid, n_cdsurfs
    logical, pointer :: cdef_surf(:,:)
 
    if( associated(fstrsolid%CONT_SGRP_ID) ) deallocate(fstrsolid%CONT_SGRP_ID)
 
    allocate(cdef_surf(l_max_elem_surf,hecmesh%n_elem))
    cdef_surf(:,:) = .false.
 
    ! label contact defined surfaces
    n_cdsurfs = 0
    do i=1, fstrsolid%n_contacts
      !grpid = fstrSOLID%contacts(i)%group
      !if( .not. fstr_isContactActive( fstrSOLID, grpid, cstep ) ) then
      !  call clear_contact_state(fstrSOLID%contacts(i));  cycle
      !endif
 
      do k=1,2 !slave,master
        if( k==1 ) then !slave
          sgrp_id = fstrsolid%contacts(i)%surf_id1_sgrp
        else if( k==2 ) then !master
          sgrp_id = fstrsolid%contacts(i)%surf_id2
        end if
 
        if( sgrp_id <= 0 ) cycle
 
        is = hecmesh%surf_group%grp_index(sgrp_id-1) + 1
        ie = hecmesh%surf_group%grp_index(sgrp_id  )
        do j=is,ie
          eid = hecmesh%surf_group%grp_item(2*j-1)
          sid = hecmesh%surf_group%grp_item(2*j)
          ! only internal and boundary element should be added
          if( .not. cdef_surf(sid,eid) ) n_cdsurfs = n_cdsurfs + 1
          cdef_surf(sid,eid) = .true.
        enddo
      end do
    enddo
 
    !gather info of contact defined surfaces
    allocate(fstrsolid%CONT_SGRP_ID(2*n_cdsurfs))
    n_cdsurfs = 0
    do i=1,hecmesh%n_elem
      do j=1,l_max_elem_surf
        if( cdef_surf(j,i) ) then
          n_cdsurfs = n_cdsurfs + 1
          fstrsolid%CONT_SGRP_ID(2*n_cdsurfs-1) = i
          fstrsolid%CONT_SGRP_ID(2*n_cdsurfs  ) = j
        endif
      end do
    end do
    deallocate(cdef_surf)
 
  end subroutine
 
  subroutine calc_contact_area( hecMESH, fstrSOLID, flag )
    type( hecmwst_local_mesh ), intent(in)       :: hecmesh
    type(fstr_solid), intent(inout)              :: fstrsolid
    integer(kind=kint), intent(in)               :: flag
 
    integer(kind=kint), parameter :: ndof=3
    integer(kind=kint) :: i, isuf, icel, sid, etype, nn, is, idx
    integer(kind=kint) :: ndlocal(l_max_elem_node)
    real(kind=kreal), allocatable :: coord(:)
    real(kind=kreal)   :: ecoord(ndof,l_max_elem_node), vect(l_max_elem_node)
 
    fstrsolid%CONT_AREA(:) = 0.d0
 
    if( .not. associated(fstrsolid%CONT_SGRP_ID) ) return
 
    allocate(coord(ndof*hecmesh%n_node))
    do i=1,ndof*hecmesh%n_node
      coord(i) = hecmesh%node(i)+fstrsolid%unode(i)
    end do
    if( flag == 1 ) then
      do i=1,ndof*hecmesh%n_node
        coord(i) = coord(i)+fstrsolid%dunode(i)
      end do
    end if
 
    do isuf=1,size(fstrsolid%CONT_SGRP_ID)/2
      icel = fstrsolid%CONT_SGRP_ID(2*isuf-1)
      sid  = fstrsolid%CONT_SGRP_ID(2*isuf  )
 
      etype = hecmesh%elem_type(icel)
      nn = hecmw_get_max_node(etype)
      is = hecmesh%elem_node_index(icel-1)
      ndlocal(1:nn) = hecmesh%elem_node_item (is+1:is+nn)
 
      do i=1,nn
        idx = ndof*(ndlocal(i)-1)
        ecoord(1:ndof,i) = coord(idx+1:idx+ndof)
      end do
 
      call calc_nodalarea_surfelement( etype, ecoord, sid, vect )
 
      do i=1,nn
        idx = ndlocal(i)
        fstrsolid%CONT_AREA(idx) = fstrsolid%CONT_AREA(idx) + vect(i)
      end do
 
    end do
 
    deallocate(coord)
  end subroutine
 
  subroutine calc_nodalarea_surfelement( etype, ecoord, sid, vect )
    integer(kind=kint), intent(in)   :: etype
    real(kind=kreal), intent(in)     :: ecoord(:,:)
    integer(kind=kint), intent(in)   :: sid
    real(kind=kreal), intent(out)    :: vect(:)
 
    integer(kind=kint), parameter :: ndof=3
    integer(kind=kint) :: nod(l_max_surface_node)
    integer(kind=kint) :: nsur, stype, ig0, i
    real(kind=kreal)   :: localcoord(2), normal(3), area, wg
    real(kind=kreal)   :: scoord(ndof,l_max_surface_node), h(l_max_surface_node)
 
    vect(:) = 0.d0
 
    call getsubface( etype, sid, stype, nod )
    nsur = getnumberofnodes( stype )
    do i=1,nsur
      scoord(1:ndof,i) = ecoord(1:ndof,nod(i))
    end do
 
    area = 0.d0
    do ig0=1,numofquadpoints( stype )
      call getquadpoint( stype, ig0, localcoord(1:2) )
      call getshapefunc( stype, localcoord(1:2), h(1:nsur) )
 
      wg=getweight( stype, ig0 )
      ! normal = dx/dr_1 \times dx/dr_2
      normal(1:3) = surfacenormal( stype, nsur, localcoord(1:2), scoord(1:ndof,1:nsur) )
      area = area + dsqrt(dot_product(normal,normal))*wg
    enddo
    area = area/dble(nsur)
    do i=1,nsur
      vect(nod(i)) = area
    end do
 
  end subroutine
 
  subroutine fstr_setup_parancon_contactvalue(hecMESH,ndof,vec,vtype)
    use m_fstr
    implicit none
    type(hecmwst_local_mesh), intent(in)      :: hecmesh
    integer(kind=kint), intent(in)            :: ndof
    real(kind=kreal), pointer, intent(inout)  :: vec(:)
    integer(kind=kint), intent(in)            :: vtype !1:value, 2:state
    !
    integer(kind=kint) ::  i,n,i0,n_loc
    integer(kind=kint) :: offset, pid, lid
    integer(kind=kint), allocatable :: displs(:)
    real(kind=kreal), allocatable   :: vec_all(:)
 
    !
    n_loc = hecmesh%nn_internal
    allocate(displs(0:nprocs))
    displs(:) = 0
    displs(myrank+1) = n_loc
    call hecmw_allreduce_i(hecmesh, displs, nprocs+1, hecmw_sum)
    do i=1,nprocs
      displs(i) = displs(i-1) + displs(i)
    end do
    offset = displs(myrank)
    n = displs(nprocs)
 
    allocate(vec_all(ndof*n))
 
    if( vtype == 1 ) then
      vec_all(:) = 0.d0
      do i= 1, hecmesh%n_node
        pid = hecmesh%node_ID(i*2)
        lid = hecmesh%node_ID(i*2-1)
        i0 = (displs(pid) + (lid-1))*ndof
        vec_all(i0+1:i0+ndof) = vec((i-1)*ndof+1:i*ndof)
        vec((i-1)*ndof+1:i*ndof) = 0.d0
      enddo
 
      call hecmw_allreduce_r(hecmesh, vec_all, n*ndof, hecmw_sum)
      do i= 1, hecmesh%n_node
        pid = hecmesh%node_ID(i*2)
        lid = hecmesh%node_ID(i*2-1)
        i0 = (displs(pid) + (lid-1))*ndof
        if (dot_product(vec_all(i0+1:i0+ndof),vec_all(i0+1:i0+ndof)) == 0.d0) cycle
        vec((i-1)*ndof+1:i*ndof) = vec_all(i0+1:i0+ndof)
      enddo
      do i=1,ndof*n_loc
        vec(i) = vec_all(offset*ndof+i)
      end do
    else if( vtype == 2 ) then
      vec_all(:) = -1000.d0
      do i= 1, hecmesh%n_node
        if( vec(i) == 0.d0 ) cycle
        pid = hecmesh%node_ID(i*2)
        lid = hecmesh%node_ID(i*2-1)
        i0 = displs(pid) + lid
        vec_all(i0) = vec(i)
      enddo
 
      call hecmw_allreduce_r(hecmesh, vec_all, n, hecmw_max)
      do i= 1, hecmesh%n_node
        pid = hecmesh%node_ID(i*2)
        lid = hecmesh%node_ID(i*2-1)
        i0 = displs(pid) + lid
        if(vec_all(i0) == -1000.d0) cycle
        vec(i) = vec_all(i0)
      enddo
    end if
 
    deallocate(displs,vec_all)
  end subroutine
 
  subroutine update_contact_state_vectors( contact, dt, relvel_vec, state_vec )
    type( tcontact ), intent(in)      :: contact
    real(kind=kreal), intent(in)      :: dt             
    real(kind=kreal), intent(inout)   :: relvel_vec(:) 
    real(kind=kreal), intent(inout)   :: state_vec(:)  
 
    integer(kind=kint)  :: i, j, slave
 
    if( contact%method == contacts2s ) then
      ! MORTAR=YES: the contact state lives per integration point in slave_surf
      ! (contact%slave is unused). Aggregate each segment's state onto its slave nodes
      ! for nodal output, using the same convention as the NODE-SURF path (free nodes
      ! show CONTACTFREE, active nodes show CONTACTSTICK/SLIP; max over segments
      ! sharing a node).
      do i = 1, size(contact%slave_surf)
        do j = 1, size(contact%slave_surf(i)%nodes)
          slave = contact%slave_surf(i)%nodes(j)
          if( state_vec(slave) < 0.1d0 .or. contact%slave_surf(i)%state > 0 ) &
          &  state_vec(slave) = dble(contact%slave_surf(i)%state)
        enddo
      enddo
      return
    endif
 
    do i= 1, size(contact%slave)
      slave = contact%slave(i)
      if( state_vec(slave) < 0.1d0 .or. contact%states(i)%state > 0 ) &
      &  state_vec(slave) = dble(contact%states(i)%state)
 
      if( is_contact_free(contact%states(i)%state) ) cycle   ! not in contact or near
      if( dt < 1.d-16 ) cycle ! too small delta t
      relvel_vec(3*slave-2:3*slave) = contact%states(i)%reldisp(1:3)/dt
    enddo
 
  end subroutine update_contact_state_vectors
 
end module m_fstr_contact_output