src/service/itm_grid_simplex.f90

Go to the documentation of this file.

module itm_grid_simplex

  !> Service routines for grids consisting out of simplexes (nodes, triangles, tetrahedra...)

  use itm_types
  use itm_assert
  use itm_combinations

  use itm_grid_common
  use itm_grid_object
  use itm_grid_access
  use itm_grid_subgrid

  use itm_grid_simplex_dict

  implicit none

  !> Definition of default subgrids
  integer, parameter :: GRID_SIMPLEX_SUBGRID_0D = 1 ! all 0d objects
  integer, parameter :: GRID_SIMPLEX_SUBGRID_1D = 2 ! all 1d objects
  integer, parameter :: GRID_SIMPLEX_SUBGRID_2D = 3 ! all 2d objects
  integer, parameter :: GRID_SIMPLEX_SUBGRID_3D = 4 ! all 3d objects

  ! Same as above, with human-readable names
  integer, parameter :: GRID_SIMPLEX_NODES     = GRID_SIMPLEX_SUBGRID_0D ! all 0d objects
  integer, parameter :: GRID_SIMPLEX_EDGES     = GRID_SIMPLEX_SUBGRID_1D ! all 1d objects
  integer, parameter :: GRID_SIMPLEX_TRIANGLES = GRID_SIMPLEX_SUBGRID_2D ! all 2d objects  
  integer, parameter :: GRID_SIMPLEX_TETRAEDER = GRID_SIMPLEX_SUBGRID_3D ! all 3d objects 

contains

  !> Write a n-dimensional grid consisting of simplexes (0d: nodes, 1d: edges, 2d: triangles, 3d: tetrahedra)
  !> into a grid descriptor, as well as the default subgrids for objects of all dimensions.
  !>
  !> The dimension n of the grid is taken as size(coordtype).
  !> 
  !> @param grid Grid descriptor to fill
  !> @param coordtype Dimension(n). Defines coordinate types / labels for the 
  !>                  individual spaces axes. See the 
  !>                  constants defined in itm_grid.f90 (COORDTYPE_*)
  !> @param nodes     The positions of the grid nodes in space (with axes as defined in coordtype). 
  !>                  First dimension is node index, second coordinate index, 
  !>                  i.e. nodes(i,:) is the position vector of node i.
  !> @param simplexes Node index lists for the invidiual n-dimensional cells/simplexes in the grid. 
  !>                  First dimension is simplex index, second dimension is node index (as referring to the nodes array),
  !>                  i.e. simplices(i,:) is the vector of node indices for simplex i.
  !>                  Grid node coordinates in the individual dimensions. 
  !> @param id        Name/Identifier of the grid.
  !> @param createSubgrids Optional flag controlling whether default subgrids
  !>                  are created. Default is .true.
  subroutine gridSetupSimplex( grid, coordtype, nodes, simplexes, id, createSubgrids )
    type(type_complexgrid), intent(out) :: grid 
    integer, dimension(:), intent(in) :: coordtype
    real(R8), dimension(:, :), intent(in) :: nodes
    integer, dimension(:, :), intent(in) :: simplexes
    character(*), intent(in), optional :: id
    logical, intent(in), optional :: createSubgrids
    
    ! internal
    integer :: idim, ndim, iSg, iObj
    type(GridObject) :: obj
    integer, allocatable :: classes(:,:)
    logical :: lCreateSubgrids = .true.

    character(9), parameter :: genId(0:3) = (/ 'nodes    ', 'edges    ', 'triangles', 'simplexes' /)

    if (present(createSubgrids)) lCreateSubgrids = createSubgrids

    call assert( size( coordtype ) == size( nodes, 2 ), &
         & "gridSetupSimplex: size of coordtype and nodes don't match" )
    call assert( size(coordtype) == size(simplexes, 2) - 1, &
         & "gridSetupSimplex: shape of simplexes is inconsistent with shape of coordtypes" )

    ndim = size(coordtype)
    allocate( grid % spaces(1) )

    ! ... fill in the grid data
    call gridSetupSimplexSpace( grid % spaces(1), coordtype, nodes, simplexes )

    ! Create subgrids 
    if (lCreateSubgrids) then
        call gridCreateDefaultSubGrids(grid, id)
    end if

  end subroutine gridSetupSimplex


  !> Set up a 1d space consisting of simplices. 
  subroutine gridSetupSimplexSpace( space, coordtype, nodes, simplexes )
    type(type_complexgrid_space), intent(inout) :: space !> The space descriptor to fill
    integer, dimension(:), intent(in) :: coordtype !> The coordinate type of the space
    real(R8), dimension(:, :), intent(in) :: nodes
    integer, dimension(:, :), intent(in) :: simplexes
   
    ! internal
    integer :: ndim, nnode, nsim
    integer :: nparents(size(nodes, 1))
    integer, allocatable :: children(:,:), simplex(:)
    type(SimplexDict) :: dictN, dictC

    integer :: isim, inode, ichild, idim
    integer :: MAX_ENTRIES = 10

    type(SimplexDictIterator) :: iter

    
    ndim = size( coordtype )  ! Simplex of dimension ndim has ndim+1 nodes
    nnode = size( nodes, 1 )
    nsim = size( simplexes, 1 )
    call assert( ndim == size(simplexes, 2) - 1 )

    ! Set coordinate types
    ! TODO: support multiple geometries
    allocate( space % coordtype(ndim, 1) )
    space % coordtype(:, 1) = coordtype

    ! Allocate object definition arrays
    allocate( space % objects(NDIM + 1) )

    ! Set up node information
    ! Geometry: normal geometry representation, third dimension unused.
    ! TODO: support multiple geometries
    allocate( space % objects(1) % geo(size(nodes, 1), ndim, 1, 1) )
    space % objects(1) % geo(:, :, 1, 1) = nodes
    ! space % nodes % xpoints, altgeo, alias: unused for this grid

    ! Do the ndim-simplices as one-offs
    
    ! Populate dictionary of ndim-1-simplices.
    ! We have to set up the entire dictionary first to have the correct global simplex indices.
    call createSimplexDict( dictC, ndim-1, nnode, MAX_ENTRIES )
    allocate(children(ndim+1, ndim))
    do isim = 1, nsim
        children = listChildren(simplexes(isim,:))
        do ichild = 1, size(children, 1)
            call insertSimplex( dictC, sortSimplex(children(ichild, :)) )
        end do
    end do

    ! Fill in object boundaries

    ! Boundary array
    allocate( space % objects(ndim+1) % boundary( nsim, ndim+1 ) )

    ! First set all to undefined
    space % objects(ndim+1) % boundary = GRID_UNDEFINED
    do isim = 1, nsim
        children = listChildren(simplexes(isim,:))
        do ichild = 1, size(children, 1)
            space % objects(ndim+1) % boundary(isim, ichild) &
                & = simplexIndex(dictC, sortSimplex(children(ichild, :)))
        end do
    end do
    deallocate(children)

    ! From here on, do lower-dimensional objects in a standard way
    do idim = ndim - 1, 1 , -1
        write (*,*) "Writing simplexes of dimension", idim
        allocate(simplex(idim+1), children(idim+1, idim))

        dictN = dictC
        call createSimplexDict( dictC, idim-1, nnode, MAX_ENTRIES )
        call resetIterator(iter)
        do 
            call nextSimplex(dictN, iter, simplex)           
            if (simplexUndefined(simplex)) exit

            children = listChildren(simplex)
            do ichild = 1, size(children, 1)
                call insertSimplex( dictC, sortSimplex(children(ichild, :)) )
            end do
        end do

        allocate( space % objects(idim+1) % boundary( countSimplexes(dictN), idim+1 ) )
        call resetIterator(iter)
        isim = 0
        do 
            call nextSimplex(dictN, iter, simplex)
            if (simplexUndefined(simplex)) exit

            isim = isim + 1
            children = listChildren(simplex)
            do ichild = 1, size(children, 1)               
                space % objects(idim+1) % boundary(isim, ichild) &
                    & = simplexIndex(dictC, sortSimplex(children(ichild, :)))                
            end do
        end do    

        deallocate(simplex, children)
    end do

  end subroutine gridSetupSimplexSpace


  subroutine writeSimplexes(simplexes)
    integer, intent(in) :: simplexes(:,:)

    ! internal
    integer :: i

    do i = 1, size(simplexes, 1)
        write (*,*) "(", simplexes(i,:), ")"
    end do

  end subroutine writeSimplexes

  
end module itm_grid_simplex