heat_mat_ass_bc_DFLUX.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_heat_mat_ass_bc_dflux
  use m_fstr
 
  implicit none
contains
  !C
  !C***
  !C*** MAT_ASS_DFLUX
  !C***
  !C
  subroutine heat_mat_ass_bc_dflux( hecMESH, hecMAT, fstrHEAT, CTIME, DTIME, beta )
    use m_heat_get_amplitude
    use m_heat_lib_dflux
    use m_static_lib_3d
    integer(kind=kint) :: k, icel, ic_type, isect, isuf, iamp, nn, is, j
    real(kind=kreal)   :: ctime, dtime, qq, val, asect, thick, vol, beta
    type(fstr_heat)          :: fstrheat
    type(hecmwst_matrix)     :: hecMAT
    type(hecmwst_local_mesh) :: hecMESH
    real(kind=kreal)   :: xx(20), yy(20), zz(20)
    real(kind=kreal)   :: vect(20), ss(2000)
    integer(kind=kint) :: ig0, ig, iS0, iE0, nodLocal(20)
    real(kind=kreal), allocatable :: bbak(:)
 
    !C
    !$omp parallel default(none), &
      !$omp&  private(k,icel,ic_type,isect,isuf,iamp,QQ,val,nn,is,j,nodLOCAL,xx,yy,zz,asect,vect,thick), &
      !$omp&  shared(fstrHEAT,CTIME,hecMAT,hecMESH)
    !$omp do
    do k = 1, fstrheat%Q_SUF_tot
 
      icel    = fstrheat%Q_SUF_elem(k)
      ic_type = hecmesh%elem_type(icel)
      isect   = hecmesh%section_ID(icel)
      isuf    = fstrheat%Q_SUF_surf(k)
      iamp    = fstrheat%Q_SUF_ampl(k)
 
      call heat_get_amplitude( fstrheat,iamp,ctime,qq )
      val     = fstrheat%Q_SUF_val (k) * qq
      if( dabs(val) < 1.d-16 ) cycle
      !C**
      nn = hecmw_get_max_node(ic_type)
      !C**
      is = hecmesh%elem_node_index(icel-1)
      do j = 1, nn
        nodlocal(j) = hecmesh%elem_node_item(is+j)
        xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
        yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
        zz(j) = hecmesh%node( 3*nodlocal(j)   )
      enddo
      !C**
      if    ( ic_type.eq.111 ) then
        is = hecmesh%section%sect_R_index(isect)
        asect = hecmesh%section%sect_R_item(is)
        call heat_dflux_111(nn,xx,yy,zz,asect,isuf,val,vect)
 
      elseif( ic_type.eq.231 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        call heat_dflux_231(nn,xx,yy,zz,thick,isuf,val,vect)
 
      elseif( ic_type.eq.232 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        call heat_dflux_232(nn,xx,yy,zz,thick,isuf,val,vect)
 
      elseif( ic_type.eq.241 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        call heat_dflux_241(nn,xx,yy,zz,thick,isuf,val,vect)
 
      elseif( ic_type.eq.242 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        call heat_dflux_242(nn,xx,yy,zz,thick,isuf,val,vect)
 
      elseif( ic_type.eq.341 ) then
        call heat_dflux_341(nn,xx,yy,zz,isuf,val,vect)
 
      elseif( ic_type.eq.342 ) then
        call heat_dflux_342(nn,xx,yy,zz,isuf,val,vect)
 
      elseif( ic_type.eq.351 ) then
        call heat_dflux_351(nn,xx,yy,zz,isuf,val,vect)
 
      elseif( ic_type.eq.352 ) then
        call heat_dflux_352(nn,xx,yy,zz,isuf,val,vect)
 
      elseif( ic_type.eq.361 ) then
        call heat_dflux_361(nn,xx,yy,zz,isuf,val,vect)
 
      elseif( ic_type.eq.362 ) then
        call heat_dflux_362(nn,xx,yy,zz,isuf,val,vect)
 
      elseif( ic_type.eq.731 ) then
        is    = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        call heat_dflux_731(nn,xx,yy,zz,thick,isuf,val,vect)
 
      elseif( ic_type.eq.741 ) then
        is    = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        call heat_dflux_741(nn,xx,yy,zz,thick,isuf,val,vect)
 
      endif
 
      do j = 1, nn
        !$omp atomic
        hecmat%B( nodlocal(j) ) = hecmat%B( nodlocal(j) ) - vect(j)
      enddo
    enddo
    !$omp end do
    !$omp end parallel
 
    if( fstrheat%WL_tot>0 ) then
      allocate( bbak(size(hecmat%B)) )
 
      do ig0 = 1, fstrheat%WL_tot
        ig = fstrheat%weldline(ig0)%egrpid
        is0= hecmesh%elem_group%grp_index(ig-1) + 1
        ie0= hecmesh%elem_group%grp_index(ig  )
        vol =0.d0; bbak=0.d0
        do k=is0, ie0
          icel   = hecmesh%elem_group%grp_item(k)
          ic_type = hecmesh%elem_type(icel)
          isect   = hecmesh%section_ID(icel)
          isuf    = 0
 
          nn = hecmw_get_max_node(ic_type)
          is = hecmesh%elem_node_index(icel-1)
          do j = 1, nn
            nodlocal(j) = hecmesh%elem_node_item(is+j)
            xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
            yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
            zz(j) = hecmesh%node( 3*nodlocal(j)   )
          enddo
 
          if( .not. isweldactive(nn, xx,yy,zz, fstrheat%weldline(ig0), ctime-0.5d0*dtime) ) cycle
 
          vol = vol + volume_c3(ic_type,nn,xx,yy,zz)
          val = fstrheat%weldline(ig0)%I * fstrheat%weldline(ig0)%U *   &
            fstrheat%weldline(ig0)%coe
          write(idbg,*) "Element", hecmesh%global_elem_id(icel),"with dflux", val, vol
 
          if( ic_type.eq.341 ) then
            call heat_dflux_341(nn,xx,yy,zz,isuf,val,vect)
 
          elseif( ic_type.eq.342 ) then
            call heat_dflux_342(nn,xx,yy,zz,isuf,val,vect)
 
          elseif( ic_type.eq.351 ) then
            call heat_dflux_351(nn,xx,yy,zz,isuf,val,vect)
 
          elseif( ic_type.eq.352 ) then
            call heat_dflux_352(nn,xx,yy,zz,isuf,val,vect)
 
          elseif( ic_type.eq.361 ) then
            call heat_dflux_361(nn,xx,yy,zz,isuf,val,vect)
 
          elseif( ic_type.eq.362 ) then
            call heat_dflux_362(nn,xx,yy,zz,isuf,val,vect)
 
 
          else
            write(*,*) '###ERROR### : Element type not supported for heat analysis'
            write(*,*) ' ic_type = ', ic_type
            call hecmw_abort(hecmw_comm_get_comm())
 
          endif
 
          do j = 1, nn
            bbak( nodlocal(j) ) = bbak( nodlocal(j) ) - vect(j)
          enddo
        enddo
 
        if( vol>0 ) then
          bbak = bbak/vol
          hecmat%B = hecmat%B + bbak
        endif
      enddo
 
      deallocate(bbak)
    endif
 
  end subroutine heat_mat_ass_bc_dflux
 
  logical function isweldactive(nn, xx,yy,zz, weldline, ctime)
    integer, intent(in)          ::  nn
    real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:)
    type(tweldline), intent(in)  :: weldline
    real(kind=kreal), intent(in) :: ctime
 
    integer :: i
    real(kind=kreal) :: cpos, wpos, tend
 
    isweldactive = .false.
    tend = weldline%tstart+ (weldline%n2-weldline%n1)/weldline%v
    if( ctime<weldline%tstart .or. ctime>tend ) return
    cpos = 0
    do i=1,nn
      if( weldline%xyz==1 ) then
        cpos = cpos+xx(i)
      elseif( weldline%xyz==2 ) then
        cpos = cpos+yy(i)
      else
        cpos = cpos+zz(i)
      endif
    enddo
    cpos = cpos/nn
    wpos = weldline%n1 + weldline%v*(ctime-weldline%tstart)
    if( dabs(cpos-wpos)<weldline%distol ) isweldactive = .true.
  end function
end module m_heat_mat_ass_bc_dflux