surf_ele.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 msurfelement
  use hecmw_util
  use elementinfo
  use m_utilities
  use bucket_search
 
  implicit none
 
  integer(kind=kint), parameter       :: l_max_surface_node =20
  integer(kind=kint), parameter       :: l_max_elem_node = 100
  integer(kind=kint), parameter       :: l_max_elem_surf = 6
 
  integer(kind=kint), parameter       :: n_neighbor_max_init = 8
 
  type tsurfelement
    integer(kind=kint)              :: eid
    integer(kind=kint)              :: etype
    integer(kind=kint), pointer     :: nodes(:)=>null()     
    integer(kind=kint)              :: n_neighbor
    integer(kind=kint)              :: n_neighbor_max
    integer(kind=kint), pointer     :: neighbor(:)=>null()  
    real(kind=kreal)                :: reflen               
    real(kind=kreal)                :: xavg(3)              
    real(kind=kreal)                :: dmax                 
    integer(kind=kint)              :: bktid
    real(kind=kreal), pointer       :: vertex_normals(:,:)=>null()  
    real(kind=kreal), pointer       :: intermediate_points(:,:)=>null()  
  end type tsurfelement
 
  integer(kind=kint), parameter, private :: debug = 0
 
contains
 
  subroutine initialize_surf( eid, etype, nsurf, surf )
    integer(kind=kint), intent(in)    :: eid
    integer(kind=kint), intent(in)    :: etype
    integer(kind=kint), intent(in)    :: nsurf
    type(tsurfelement), intent(inout) :: surf
    integer(kind=kint) :: i, n, outtype, nodes(100)
    surf%eid = eid
 
    if( nsurf > 0 ) then
      call getsubface( etype, nsurf, outtype, nodes )
      surf%etype = outtype
      n=getnumberofnodes( outtype )
    else !treat 3d master element as it is
      surf%etype = etype
      n=getnumberofnodes( etype )
      do i=1,n
        nodes(i) = i
      enddo
    endif
    allocate( surf%nodes(n) )
    surf%nodes(1:n)=nodes(1:n)
    surf%n_neighbor = 0
    surf%n_neighbor_max = 0
    surf%reflen  = -1.d0
    surf%xavg(:) =  0.d0
    surf%dmax    = -1.d0
    surf%bktID   = -1
    allocate( surf%vertex_normals(3,n) )
    surf%vertex_normals(:,:) = 0.d0
  end subroutine
 
  subroutine finalize_surf( surf )
    type(tsurfelement), intent(inout) :: surf
    if( associated(surf%nodes) ) deallocate( surf%nodes )
    if( associated(surf%neighbor) ) deallocate( surf%neighbor )
    if( associated(surf%vertex_normals) ) deallocate( surf%vertex_normals )
    if( associated(surf%intermediate_points) ) deallocate( surf%intermediate_points )
  end subroutine
 
  subroutine write_surf( file, surf )
    integer(kind=kint), intent(in)    :: file
    type(tsurfelement), intent(in)    :: surf
    integer(kind=kint) :: i
    write(file,*) "Element:",surf%eid,"Surface type=",surf%etype
    if( associated( surf%nodes ) )   &
      write(file,*) ( surf%nodes(i),i=1,size(surf%nodes) )
    if( associated( surf%neighbor ) )   &
      write(file,*) ( surf%neighbor(i),i=1,surf%n_neighbor )
  end subroutine
 
  subroutine find_surface_neighbor( surf, bktDB )
    type(tsurfelement), intent(inout) :: surf(:)
    type(bucketdb), intent(in) :: bktDB
    integer(kind=kint) :: i, j, ii,jj, nd1, nd2, nsurf
    integer(kind=kint) :: k, oldsize, newsize
    integer(kind=kint), pointer :: dumarray(:) => null()
    integer(kind=kint) :: bktID, ncand, js
    integer(kind=kint), allocatable :: indexSurf(:)
    if (debug >= 1) write(0,*) 'DEBUG: find_surface_neighbor: start'
 
    nsurf = size(surf)
 
    !$omp parallel do default(none), &
      !$omp&private(i,ii,nd1,j,jj,nd2,oldsize,newsize,dumarray,k,bktID,ncand,indexSurf,js), &
      !$omp&shared(nsurf,surf,bktDB)
    do i=1,nsurf
      bktid = bucketdb_getbucketid(bktdb, surf(i)%xavg)
      ncand = bucketdb_getnumcand(bktdb, bktid)
      if (ncand == 0) cycle
      allocate(indexsurf(ncand))
      call bucketdb_getcand(bktdb, bktid, ncand, indexsurf)
      jloop: do js=1,ncand
        j = indexsurf(js)
        if( i==j ) cycle
        if( associated(surf(i)%neighbor) ) then
          if ( any( surf(i)%neighbor(1:surf(i)%n_neighbor)==j ) ) cycle
        endif
        do ii=1, size(surf(i)%nodes)
          nd1 = surf(i)%nodes(ii)
          do jj=1, size(surf(j)%nodes)
            nd2 = surf(j)%nodes(jj)
            if( nd1==nd2 ) then
              surf(i)%n_neighbor = surf(i)%n_neighbor+1
 
              if( .not. associated(surf(i)%neighbor) ) then
                allocate( surf(i)%neighbor(n_neighbor_max_init) )
                surf(i)%n_neighbor_max = n_neighbor_max_init
                surf(i)%neighbor(1) = j
              else if( surf(i)%n_neighbor > surf(i)%n_neighbor_max ) then
                oldsize = surf(i)%n_neighbor_max
                newsize = oldsize * 2
                dumarray => surf(i)%neighbor
                allocate( surf(i)%neighbor(newsize) )
                surf(i)%n_neighbor_max = newsize
                do k=1,oldsize
                  surf(i)%neighbor(k) = dumarray(k)
                enddo
                surf(i)%neighbor(oldsize+1) = j
                deallocate( dumarray )
              else
                surf(i)%neighbor(surf(i)%n_neighbor) = j
              endif
 
              cycle jloop
            endif
          enddo
        enddo
      enddo jloop
      deallocate(indexsurf)
    enddo
    !$omp end parallel do
 
    if (debug >= 1) write(0,*) 'DEBUG: find_surface_neighbor: end'
  end subroutine
 
  integer(kind=kint) function next_position( surf, cpos )
    type(tsurfelement), intent(in) :: surf
    real(kind=kreal), intent(in)   :: cpos(2)   
 
    integer(kind=kint)          :: i
    real(kind=kreal) :: maxv(3)
    next_position = surf%eid
    if( .not. associated(surf%neighbor) ) return   ! do nothing when not initialized
    maxv(:) = 0.d0
    select case(surf%etype)
      case( fe_tri3n, fe_tri6n )
        if( all(cpos>0.d0) .and. all(cpos<1.d0) ) return
        if( size(surf%neighbor)<3 )  return
        if( all(cpos(:)>0.d0) ) return
        do i=1,3
          if( cpos(i)< 0.d0 ) maxv(i) = -cpos(i)
        enddo
        next_position = maxloc(maxv,1)
      case( fe_quad4n, fe_quad8n )
        if( all(cpos>-1.d0) .and. all(cpos<1.d0) ) return
        if( size(surf%neighbor)<4 ) return
        if( cpos(1)<-1.d0 .and. dabs(cpos(2))<1.d0 ) then
          next_position = 1
        elseif( cpos(1)>1.d0 .and. dabs(cpos(2))<1.d0 ) then
          next_position = 3
        elseif( dabs(cpos(1))<1.d0 .and. cpos(2)<-1.d0 ) then
          next_position = 2
        elseif( dabs(cpos(1))<1.d0 .and. cpos(2)>1.d0 ) then
          next_position = 4
        elseif( cpos(1)<-1.d0 .and. cpos(2)<-1.d0 ) then
          if( cpos(1)>cpos(2) ) then
            next_position = 2
          else
            next_position = 1
          endif
        elseif( cpos(1)<-1.d0 .and. cpos(2)>1.d0 ) then
          if( dabs(cpos(1))>cpos(2) ) then
            next_position = 1
          else
            next_position = 4
          endif
        elseif( cpos(1)>1.d0 .and. cpos(2)<-1.d0 ) then
          if( cpos(1)>dabs(cpos(2)) ) then
            next_position = 3
          else
            next_position = 2
          endif
        elseif( cpos(1)>1.d0 .and. cpos(2)>1.d0 ) then
          if( cpos(1)>cpos(2) ) then
            next_position =  3
          else
            next_position = 4
          endif
        endif
      case default
        stop "type of surface element not defined"
    end select
    next_position = surf%neighbor( next_position )
 
  end function next_position
 
  subroutine update_surface_reflen( surf, coord )
    type(tsurfelement), intent(inout) :: surf(:)
    real(kind=kreal), intent(in)      :: coord(:)  
    real(kind=kreal) :: elem(3, l_max_surface_node), r0(2)
    integer(kind=kint) :: nn, i, j, iss
    !$omp parallel do default(none) private(j,nn,i,iss,elem,r0) shared(surf,coord)
    do j=1,size(surf)
      nn = size(surf(j)%nodes)
      do i=1,nn
        iss = surf(j)%nodes(i)
        elem(1:3,i) = coord(3*iss-2:3*iss)
      enddo
      call getelementcenter( surf(j)%etype, r0 )
      surf(j)%reflen = getreferencelength( surf(j)%etype, nn, r0, elem )
    enddo
    !$omp end parallel do
  end subroutine update_surface_reflen
 
  subroutine update_surface_box_info( surf, currpos )
    type(tsurfelement), intent(inout) :: surf(:)
    real(kind=kreal), intent(in) :: currpos(:)      
    integer(kind=kint) :: nn, i, j, iss
    real(kind=kreal) :: elem(3,l_max_surface_node),xmin(3), xmax(3), xsum(3), lx(3)
    !$omp parallel do default(none) private(i,nn,iss,elem,xmin,xmax,xsum,lx) shared(surf,currpos)
    do j=1, size(surf)
      nn = size(surf(j)%nodes)
      do i=1,nn
        iss = surf(j)%nodes(i)
        elem(1:3,i) = currpos(3*iss-2:3*iss)
      enddo
      do i=1,3
        xmin(i) = minval(elem(i,1:nn))
        xmax(i) = maxval(elem(i,1:nn))
        xsum(i) = sum(elem(i,1:nn))
      enddo
      surf(j)%xavg(:) = xsum(:) / nn
      lx(:) = xmax(:) - xmin(:)
      surf(j)%dmax = maxval(lx) * 0.5d0
    enddo
    !$omp end parallel do
  end subroutine update_surface_box_info
 
  logical function is_in_surface_box(surf, x0, exp_rate)
    type(tsurfelement), intent(inout) :: surf
    real(kind=kreal), intent(in) :: x0(3)      
    real(kind=kreal), intent(in) :: exp_rate   
    real(kind=kreal) :: dist(3), er
    er = max(exp_rate, 1.d0)
    dist(:) = abs(x0(:) - surf%xavg(:))
    if ( maxval(dist(:)) < surf%dmax * er ) then
      is_in_surface_box = .true.
    else
      is_in_surface_box = .false.
    endif
  end function is_in_surface_box
 
  subroutine update_surface_bucket_info(surf, bktDB)
    type(tsurfelement), intent(inout) :: surf(:)
    type(bucketdb), intent(inout) :: bktDB
    real(kind=kreal) :: x_min(3), x_max(3), d_max
    integer(kind=kint) :: nsurf, i, j
    if (debug >= 1) write(0,*) 'DEBUG: update_surface_bucket_info: start'
    nsurf = size(surf)
    if (nsurf == 0) return
    x_min(:) = surf(1)%xavg(:)
    x_max(:) = surf(1)%xavg(:)
    d_max = surf(1)%dmax*2.d0
    do i = 2, nsurf
      do j = 1, 3
        if (surf(i)%xavg(j) < x_min(j)) x_min(j) = surf(i)%xavg(j)
        if (surf(i)%xavg(j) > x_max(j)) x_max(j) = surf(i)%xavg(j)
        if (surf(i)%dmax > d_max) d_max = surf(i)%dmax
      enddo
    enddo
    do j = 1, 3
      x_min(j) = x_min(j) - d_max
      x_max(j) = x_max(j) + d_max
    enddo
    call bucketdb_setup(bktdb, x_min, x_max, d_max, nsurf)
    !$omp parallel do default(none) private(i) shared(nsurf,surf,bktDB)
    do i = 1, nsurf
      surf(i)%bktID = bucketdb_getbucketid(bktdb, surf(i)%xavg)
      call bucketdb_registerpre(bktdb, surf(i)%bktID)
    enddo
    !$omp end parallel do
    call bucketdb_allocate(bktdb)
    do i = 1, nsurf
      call bucketdb_register(bktdb, surf(i)%bktID, i)
    enddo
    if (debug >= 1) write(0,*) 'DEBUG: update_surface_bucket_info: end'
  end subroutine update_surface_bucket_info
 
  subroutine calc_surf_vertex_normals( surf, elem )
    use elementinfo, only: getvertexcoord, surfacenormal
    type(tsurfelement), intent(inout) :: surf
    real(kind=kreal), intent(in)      :: elem(3,*)  
    
    integer(kind=kint) :: j, nn, etype
    real(kind=kreal) :: cnpos(2), normal(3)
    
    etype = surf%etype
    nn = size(surf%nodes)
    
    do j = 1, nn
      call getvertexcoord(etype, j, cnpos)
      normal = surfacenormal(etype, nn, cnpos, elem)
      normal = normal / dsqrt(dot_product(normal, normal))
      surf%vertex_normals(:, j) = normal
    enddo
    
  end subroutine calc_surf_vertex_normals
 
  subroutine calc_all_surf_vertex_normals( surf, currpos )
    type(tsurfelement), intent(inout) :: surf(:)
    real(kind=kreal), intent(in)      :: currpos(:) 
    
    integer(kind=kint) :: i, j, nn, iss
    real(kind=kreal) :: elem(3, l_max_elem_node)
    
    !$omp parallel do default(none) private(i,j,nn,iSS,elem) shared(surf,currpos)
    do i = 1, size(surf)
      nn = size(surf(i)%nodes)
      do j = 1, nn
        iss = surf(i)%nodes(j)
        elem(1:3, j) = currpos(3*iss-2:3*iss)
      enddo
      call calc_surf_vertex_normals(surf(i), elem(1:3,1:nn))
    enddo
    !$omp end parallel do
    
  end subroutine calc_all_surf_vertex_normals
 
end module msurfelement