hecmw_es_mesh_connectivity.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 hecmw_es_mesh_connectivity
  use hecmw_util
  use hecmw_varray_int
  use hecmw_etype
  implicit none
 
  private
  public :: hecmw_create_smoothing_element_connectivity
 
contains
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  common functions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  subroutine concat_int_list(inlist, n_in, outlist)
    integer, allocatable, intent(in) :: inlist(:)
    integer, intent(in)              :: n_in
    integer, pointer, intent(inout)  :: outlist(:)
 
    integer :: i, n_ori, n_new, lbd
    integer, allocatable :: tmp(:)
 
    n_ori = size(outlist)
    lbd = lbound(outlist,1)
    allocate(tmp(n_ori))
 
    do i=1,n_ori
      tmp(i) = outlist(lbd+i-1)
    end do
 
    n_new = n_ori + n_in
 
    deallocate(outlist)
    allocate(outlist(lbd:lbd+n_new-1))
    do i=1,n_ori
      outlist(lbd+i-1) = tmp(i)
    end do
    do i=1,n_in
      outlist(lbd+i-1+n_ori) = inlist(i)
    end do
 
    deallocate(tmp)
  end subroutine
 
  subroutine add_newelements_to_hecmesh( etype, elems, sectionID, hecMESH )
    integer, intent(in)                       :: etype
    type( hecmwST_varray_int ), allocatable, intent(in) :: elems(:)
    integer, allocatable, intent(in)          :: sectionID(:)
    type(hecmwST_local_mesh), intent(inout)   :: hecMESH
 
    integer(kind=kint) :: i
    integer(kind=kint) :: n_newelem, n_newelemnodes, nitem
    integer(kind=kint), allocatable :: iwork(:)
 
    n_newelem = size(elems)
    n_newelemnodes = 0
    do i=1,n_newelem
      n_newelemnodes = n_newelemnodes + hecmw_varray_int_get_nitem(elems(i))
    end do
    allocate(iwork(0:n_newelemnodes))
 
    ! elem_type_index(:)
    iwork(1) = hecmesh%elem_type_index(hecmesh%n_elem_type)+n_newelem
    call concat_int_list(iwork, 1, hecmesh%elem_type_index)
 
    ! elem_type_item(:)
    iwork(1) = etype
    call concat_int_list(iwork, 1, hecmesh%elem_type_item)
 
    ! elem_type(:)
    iwork(1:n_newelem) = etype
    call concat_int_list(iwork, n_newelem, hecmesh%elem_type)
 
    ! section_ID(:)
    call concat_int_list(sectionid, n_newelem, hecmesh%section_ID)
 
    ! elem_mat_ID_index(:)
    iwork(0) = hecmesh%elem_mat_ID_index(ubound(hecmesh%elem_mat_ID_index,1))
    do i=1,n_newelem
      iwork(i) = iwork(i-1) + 1
    end do
    call concat_int_list(iwork, n_newelem, hecmesh%elem_mat_ID_index)
 
    ! elem_mat_ID_item(:)
 
    ! elem_node_index(:)
    iwork(0) = hecmesh%elem_node_index(ubound(hecmesh%elem_node_index,1))
    do i=1,n_newelem
      iwork(i) = iwork(i-1) + hecmw_varray_int_get_nitem(elems(i))
    end do
    call concat_int_list(iwork, n_newelem, hecmesh%elem_node_index)
 
    ! elem_node_item(:)
    n_newelemnodes = 0
    do i=1,n_newelem
      nitem = hecmw_varray_int_get_nitem(elems(i))
      call hecmw_varray_int_get_item_all( elems(i), iwork(n_newelemnodes+1:n_newelemnodes+nitem) )
      n_newelemnodes = n_newelemnodes + nitem
    end do
    call concat_int_list(iwork, n_newelemnodes, hecmesh%elem_node_item)
 
    ! elem_ID(:)
    do i=1,n_newelem
      iwork(2*i-1) = hecmesh%n_elem + 1
      iwork(2*i  ) = 0
    end do
    call concat_int_list(iwork, 2*n_newelem, hecmesh%elem_ID)
 
    ! global_elem_ID(:)
    do i=1,n_newelem
      iwork(i) = 0
    end do
    call concat_int_list(iwork, n_newelem, hecmesh%global_elem_ID)
 
    ! elem_internal_list(:)
    ! elem_mat_int_index(:)
    ! elem_mat_int_val(:)
    ! elem_val_index(:)
    ! elem_val_item(:)
 
    ! n_elem
    !hecMESH%n_elem = hecMESH%n_elem + n_newelem
 
    ! n_elem_gross
    ! ne_internal
    !hecMESH%ne_internal = hecMESH%ne_internal + n_newelem
 
    ! n_elem_type
    hecmesh%n_elem_type = hecmesh%n_elem_type + 1
 
    ! n_elem_mat_ID
 
  end subroutine
 
  subroutine reorder_tet1( nid, nodlocal )
    integer, intent(in)    :: nid
    integer, intent(inout) :: nodlocal(4)
 
    integer :: i, idx
    integer :: nodlocal_ori(4)
 
    nodlocal_ori(1:4) = nodlocal(1:4)
    idx = 0
    do i=1,4
      if( nodlocal(i) == nid ) then
        idx = i
        exit
      end if
    end do
 
    if( idx == 1 ) then
      return
    else if( idx == 2 ) then
      nodlocal(1) = nodlocal_ori(2)
      nodlocal(2) = nodlocal_ori(3)
      nodlocal(3) = nodlocal_ori(1)
    else if( idx == 3 ) then
      nodlocal(1) = nodlocal_ori(3)
      nodlocal(2) = nodlocal_ori(1)
      nodlocal(3) = nodlocal_ori(2)
    else if( idx == 4 ) then
      nodlocal(1) = nodlocal_ori(4)
      nodlocal(2) = nodlocal_ori(2)
      nodlocal(3) = nodlocal_ori(1)
      nodlocal(4) = nodlocal_ori(3)
    end if
 
  end subroutine
 
  subroutine reorder_tet1_edge( nid1, nid2, nodlocal )
    integer, intent(in)    :: nid1
    integer, intent(in)    :: nid2
    integer, intent(inout) :: nodlocal(4)
 
    integer :: i, idx
    integer :: nodlocal_ori(4)
 
    call reorder_tet1( nid1, nodlocal )
 
    nodlocal_ori(1:4) = nodlocal(1:4)
    idx = 1
    do i=2,4
      if( nodlocal(i) == nid2 ) then
        idx = i
        exit
      end if
    end do
 
    if( idx == 2 ) then
      return
    else if( idx == 3 ) then
      nodlocal(2) = nodlocal_ori(3)
      nodlocal(3) = nodlocal_ori(4)
      nodlocal(4) = nodlocal_ori(2)
    else if( idx == 4 ) then
      nodlocal(2) = nodlocal_ori(4)
      nodlocal(3) = nodlocal_ori(2)
      nodlocal(4) = nodlocal_ori(3)
    else
      stop "error in reorder_tet1_edge"
    end if
 
  end subroutine
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  create node-smoothed elements
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  subroutine initialize_node_elem_list( hecMESH, n_elem_ori, is_selem_list, node_elem_list )
    type(hecmwST_local_mesh),target              :: hecMESH
    integer(kind=kint), intent(in)               :: n_elem_ori
    logical, allocatable, intent(in)             :: is_selem_list(:)
    type( hecmwST_varray_int ), allocatable, intent(inout) :: node_elem_list(:)
 
    integer(kind=kint) :: i, icel, iS, iE, nid
 
    call hecmw_varray_int_initialize_all( node_elem_list, hecmesh%n_node, 8 )
 
    do icel=1,n_elem_ori
      if( .not. is_selem_list(icel) ) cycle
 
      is = hecmesh%elem_node_index(icel-1)+1
      ie = hecmesh%elem_node_index(icel)
 
      do i=is,ie
        nid = hecmesh%elem_node_item(i)
        call hecmw_varray_int_add( node_elem_list(nid), icel )
      end do
    end do
 
  end subroutine
 
  subroutine create_node_smoothed_elements( hecMESH, node_elem_list, nselems, sectionID )
    type(hecmwST_local_mesh), intent(in)         :: hecMESH
    type( hecmwST_varray_int ), allocatable, intent(in)    :: node_elem_list(:)
    type( hecmwST_varray_int ), allocatable, intent(inout) :: nselems(:)
    integer, allocatable, intent(inout)          :: sectionID(:)
 
    integer(kind=kint) :: n_nselems
    type( hecmwST_varray_int ), allocatable :: n_sections(:)
    integer(kind=kint) :: i, j, k, icel, isect, iS, iE
    integer(kind=kint) :: nodlocal(4)
 
    !get num of node-smoothed elements
    call hecmw_varray_int_initialize_all( n_sections, hecmesh%n_node, 2 )
 
    n_nselems = 0
    do i=1,hecmesh%n_node
      do j=1,hecmw_varray_int_get_nitem(node_elem_list(i))
        icel = hecmw_varray_int_get_item( node_elem_list(i), j )
        isect = hecmesh%section_ID(icel)
        call hecmw_varray_int_add_if_not_exits( n_sections(i), isect )
      enddo
      n_nselems = n_nselems + hecmw_varray_int_get_nitem(n_sections(i))
    end do
 
    allocate(sectionid(n_nselems))
 
    !create node-smoothed elements
    !elements which belongs to different sections must not be connected
    call hecmw_varray_int_initialize_all( nselems, n_nselems, 16 )
 
    n_nselems = 0
    do i=1,hecmesh%n_node
      do j=1,hecmw_varray_int_get_nitem( n_sections(i) )
        isect = hecmw_varray_int_get_item( n_sections(i), j )
        n_nselems = n_nselems + 1
        sectionid(n_nselems) = isect
        do k=1,hecmw_varray_int_get_nitem(node_elem_list(i))
          icel = hecmw_varray_int_get_item( node_elem_list(i), k )
          if ( hecmesh%section_ID(icel) == isect ) then
            is = hecmesh%elem_node_index(icel-1)+1
            ie = hecmesh%elem_node_index(icel)
            if( ie-is+1 /= 4 ) stop "error in add_elemments_smoothed_by_node(1)"
            nodlocal(1:4) = hecmesh%elem_node_item(is:ie)
            call reorder_tet1( i, nodlocal )
            if( hecmw_varray_int_get_nitem(nselems(n_nselems)) == 0 ) call hecmw_varray_int_add( nselems(n_nselems), nodlocal(1) )
            call hecmw_varray_int_expand( nselems(n_nselems), 3, nodlocal(2:4) )
          end if
        end do
      enddo
    end do
 
    !finalize n_sections
    call hecmw_varray_int_finalize_all( n_sections )
 
  end subroutine
 
  subroutine add_elemments_smoothed_by_node( hecMESH, n_elem_ori, is_selem_list )
    type(hecmwST_local_mesh),target  :: hecMESH
    integer(kind=kint), intent(in)   :: n_elem_ori
    logical, allocatable, intent(in) :: is_selem_list(:)
 
    type( hecmwST_varray_int ), allocatable :: node_elem_list(:)
    type( hecmwST_varray_int ), allocatable :: nselems(:)
    integer, allocatable :: sectionID(:)
 
 
    !initialize node_elem_list
    call initialize_node_elem_list( hecmesh, n_elem_ori, is_selem_list, node_elem_list )
 
    !create node-smoothed elements and store them in list data
    call create_node_smoothed_elements( hecmesh, node_elem_list, nselems, sectionid )
 
    !finalize node_elem_list
    call hecmw_varray_int_finalize_all( node_elem_list )
 
    !add node-smoothed elements to hecMESH
    call add_newelements_to_hecmesh( 881, nselems, sectionid, hecmesh )
 
    !finalize nselems
    call hecmw_varray_int_finalize_all( nselems )
 
  end subroutine
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  create edge-smoothed elements
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  subroutine initialize_edge_elem_list( hecMESH, n_elem_ori, is_selem_list, edges, edge_elem_list )
    type(hecmwST_local_mesh),target              :: hecMESH
    integer(kind=kint), intent(in)               :: n_elem_ori
    logical, allocatable, intent(in)             :: is_selem_list(:)
    integer, allocatable, intent(inout)          :: edges(:)
    type( hecmwST_varray_int ), allocatable, intent(inout) :: edge_elem_list(:)
 
    type( hecmwST_varray_int ), allocatable :: edge_con(:)
    type( hecmwST_varray_int ), allocatable :: edge_id(:)
    integer(kind=kint) :: ndlocal(4)
    integer(kind=kint), parameter :: leid(2,6) = reshape((/1,2,2,3,3,1,1,4,2,4,3,4/),shape(leid))
 
    integer(kind=kint) :: i, j, icel, iS, iE, nid1, nid2
    integer(kind=kint) :: n_edges, edid
 
    !create edge connectivity list in node num order
    !edge_con(i)%items : list of nodes connecting to i by edge
    call hecmw_varray_int_initialize_all( edge_con, hecmesh%n_node, 8 )
 
    do icel=1,n_elem_ori
      if( .not. is_selem_list(icel) ) cycle
 
      is = hecmesh%elem_node_index(icel-1)+1
      ie = hecmesh%elem_node_index(icel)
 
      ndlocal(1:4) = hecmesh%elem_node_item(is:ie)
      do i=1,6
        nid1 = ndlocal(leid(1,i))
        nid2 = ndlocal(leid(2,i))
        if( nid1 < nid2 ) then
          call hecmw_varray_int_add_if_not_exits( edge_con(nid1), nid2 )
        else
          call hecmw_varray_int_add_if_not_exits( edge_con(nid2), nid1 )
        end if
      end do
    end do
 
    !create edges and edge_id
    n_edges = 0
    do i=1,hecmesh%n_node
      n_edges = n_edges + hecmw_varray_int_get_nitem(edge_con(i))
    end do
 
    allocate(edges(2*n_edges))
    call hecmw_varray_int_initialize_all( edge_id, hecmesh%n_node, 8 )
    n_edges = 0
    do i=1,hecmesh%n_node
      do j=1,hecmw_varray_int_get_nitem(edge_con(i))
        n_edges = n_edges + 1
        edges(2*n_edges-1) = i
        edges(2*n_edges  ) = hecmw_varray_int_get_item( edge_con(i), j )
        call hecmw_varray_int_add( edge_id(i), n_edges )
      end do
    end do
 
    !create edge-elem list
    call hecmw_varray_int_initialize_all( edge_elem_list, n_edges, 8 )
 
    do icel=1,n_elem_ori
      if( .not. is_selem_list(icel) ) cycle
 
      is = hecmesh%elem_node_index(icel-1)+1
      ie = hecmesh%elem_node_index(icel)
 
      ndlocal(1:4) = hecmesh%elem_node_item(is:ie)
      do i=1,6
        nid1 = ndlocal(leid(1,i))
        nid2 = ndlocal(leid(2,i))
        if( nid1 < nid2 ) then
          edid = hecmw_varray_int_get_item( edge_id(nid1), hecmw_varray_int_find( edge_con(nid1), nid2 ) )
        else
          edid = hecmw_varray_int_get_item( edge_id(nid2), hecmw_varray_int_find( edge_con(nid2), nid1 ) )
        end if
        call hecmw_varray_int_add( edge_elem_list(edid), icel )
      end do
    end do
 
    call hecmw_varray_int_finalize_all( edge_con )
    call hecmw_varray_int_finalize_all( edge_id )
 
  end subroutine
 
  subroutine create_edge_smoothed_elements( hecMESH, edges, edge_elem_list, eselems, sectionID )
    type(hecmwST_local_mesh), intent(in)         :: hecMESH
    integer(kind=kint), allocatable, intent(in)  :: edges(:)
    type( hecmwST_varray_int ), allocatable, intent(in)    :: edge_elem_list(:)
    type( hecmwST_varray_int ), allocatable, intent(inout) :: eselems(:)
    integer, allocatable, intent(inout)          :: sectionID(:)
 
    integer(kind=kint) :: n_edges, n_eselems
    type( hecmwST_varray_int ), allocatable :: n_sections(:)
    integer(kind=kint) :: i, j, k, icel, isect, iS, iE
    integer(kind=kint) :: nodlocal(4)
 
    n_edges = size(edge_elem_list)
 
    !get num of edge-smoothed elements
    call hecmw_varray_int_initialize_all( n_sections, n_edges, 2 )
 
    n_eselems = 0
    do i=1,n_edges
      do j=1,hecmw_varray_int_get_nitem(edge_elem_list(i))
        icel = hecmw_varray_int_get_item( edge_elem_list(i), j )
        isect = hecmesh%section_ID(icel)
        call hecmw_varray_int_add_if_not_exits( n_sections(i), isect )
      enddo
      n_eselems = n_eselems + hecmw_varray_int_get_nitem(n_sections(i))
    end do
 
    allocate(sectionid(n_eselems))
 
    !create edge-smoothed elements
    !elements which belongs to different sections must not be connected
    call hecmw_varray_int_initialize_all( eselems, n_eselems, 16 )
 
    n_eselems = 0
    do i=1,n_edges
      do j=1,hecmw_varray_int_get_nitem(n_sections(i))
        isect = hecmw_varray_int_get_item( n_sections(i), j )
        n_eselems = n_eselems + 1
        sectionid(n_eselems) = isect
        do k=1,hecmw_varray_int_get_nitem(edge_elem_list(i))
          icel = hecmw_varray_int_get_item( edge_elem_list(i), k )
          if ( hecmesh%section_ID(icel) == isect ) then
            is = hecmesh%elem_node_index(icel-1)+1
            ie = hecmesh%elem_node_index(icel)
            if( ie-is+1 /= 4 ) stop "error in add_elemments_smoothed_by_node(1)"
            nodlocal(1:4) = hecmesh%elem_node_item(is:ie)
            call reorder_tet1_edge( edges(2*i-1), edges(2*i), nodlocal )
            if( hecmw_varray_int_get_nitem(eselems(n_eselems)) == 0 ) &
              &  call hecmw_varray_int_expand( eselems(n_eselems), 2, nodlocal(1:2) )
            call hecmw_varray_int_expand( eselems(n_eselems), 2, nodlocal(3:4) )
          end if
        end do
      enddo
    end do
 
    !finalize n_sections
    call hecmw_varray_int_finalize_all( n_sections )
 
  end subroutine
 
  subroutine add_elemments_smoothed_by_edge( hecMESH, n_elem_ori, is_selem_list )
    type(hecmwST_local_mesh),target :: hecMESH
    integer(kind=kint), intent(in)  :: n_elem_ori
    logical, allocatable, intent(in) :: is_selem_list(:)
 
    integer, allocatable          :: edges(:)
    type( hecmwST_varray_int ), allocatable :: edge_elem_list(:)
    type( hecmwST_varray_int ), allocatable :: eselems(:)
    integer, allocatable :: sectionID(:)
 
 
    !initialize edge_elem_list
    call initialize_edge_elem_list( hecmesh, n_elem_ori, is_selem_list, edges, edge_elem_list )
 
    !create edge-smoothed elements and store them in list data
    call create_edge_smoothed_elements( hecmesh, edges, edge_elem_list, eselems, sectionid )
 
    !finalize edges and edge_elem_list
    deallocate(edges)
    call hecmw_varray_int_finalize_all( edge_elem_list )
 
    !add edge-smoothed elements to hecMESH
    call add_newelements_to_hecmesh( 891, eselems, sectionid, hecmesh )
 
    !finalize eselems
    call hecmw_varray_int_finalize_all( eselems )
 
  end subroutine
 
  subroutine hecmw_create_smoothing_element_connectivity( hecMESH, is_selem_list )
    type(hecmwst_local_mesh),target :: hecmesh
    logical, allocatable, intent(in) :: is_selem_list(:)
 
    integer(kind=kint) :: n_elem_ori
 
    n_elem_ori = hecmesh%n_elem
 
    call add_elemments_smoothed_by_node( hecmesh, n_elem_ori, is_selem_list )
 
    call add_elemments_smoothed_by_edge( hecmesh, n_elem_ori, is_selem_list )
 
  end subroutine
 
end module