GridConnection.f90

Go to the documentation of this file.

!> Refactoring issues towards parallel:
!!
!! * remove camelCase
!<
module gridconnectionmodule
  use kindmodule, only: i4b, dp, lgp
  use simmodule, only: ustop
  use constantsmodule, only: lenmempath, dzero, dpio180, lenmodelname
  use characterstringmodule
  use memorymanagermodule, only: mem_allocate, mem_deallocate
  use memoryhelpermodule, only: create_mem_path
  use listmodule, only: listtype, isequaliface
  use numericalmodelmodule
  use disumodule
  use disconnexchangemodule
  use virtualmodelmodule, only: virtualmodeltype, get_virtual_model, &
                                get_virtual_model_from_list
  use virtualexchangemodule, only: virtualexchangetype, get_virtual_exchange
  use cellwithnbrsmodule
  use connectionsmodule
  use sparsemodule, only: sparsematrix
  use interfacemapmodule
  use basedismodule, only: dis_transform_xy
  use csrutilsmodule
  use stlvecintmodule
  implicit none
  private
 
  ! Initial nr of neighbors for sparse matrix allocation
  integer(I4B), parameter :: initnrneighbors = 7
 
  !> This class is used to construct the connections object for
  !! the interface model's spatial discretization/grid.
  !!
  !! It works as follows:
  !!
  !! 1: construct basic instance, allocate data structures
  !!    based on nr. of primary connections
  !! 2: add primary connections from exchanges
  !! 3: add secondary, tertiary, ... MODEL connections, depending
  !!    on the size of the computational stencil
  !! 4: extend the connection, creating the full data structure
  !!    and relations
  !! 5: build the connections object, from which the grid for
  !!    the interface model can be constructed
  !!
  !! A note on the different indices:
  !!
  !! We have
  !!
  !! - GLOBAL index, which technically labels the row in the solution matrix
  !! - REGIONAL index, running over all models that participate in the interface grid
  !! - LOCAL index, local to each model
  !! - INTERFACE index, numbering the cells in the interface grid
  !<
  type, public :: gridconnectiontype
 
    character(len=LENMEMPATH) :: memorypath
    integer(I4B) :: internalstencildepth !< stencil size for the interior
    integer(I4B) :: exchangestencildepth !< stencil size at the interface
 
    class(numericalmodeltype), pointer :: model => null() !< the model for which this grid connection exists
    class(disconnexchangetype), pointer :: primaryexchange => null() !< pointer to the primary exchange for this interface
 
    integer(I4B), pointer :: nrofboundarycells => null() !< nr of boundary cells with connection to another model
    type(cellwithnbrstype), dimension(:), pointer :: boundarycells => null() !< cells on our side of the primary connections
    type(cellwithnbrstype), dimension(:), pointer :: connectedcells => null() !< cells on the neighbors side of the primary connection
    type(stlvecint), pointer :: haloexchanges !< all exchanges that are potentially part of this interface
 
    integer(I4B), pointer :: nrofcells => null() !< the total number of cells in the interface
    type(globalcelltype), dimension(:), pointer :: idxtoglobal => null() !< a map from interface index to global coordinate
    integer(I4B), dimension(:), pointer, contiguous :: idxtoglobalidx => null() !< a (flat) map from interface index to global index,
                                                                                !! stored in mem. mgr. so can be used for debugging
    type(listtype) :: regionalmodels !< the models participating in the interface
    integer(I4B), dimension(:), pointer :: region_to_iface_map => null() !< (sparse) mapping from regional index to interface ixd
    integer(I4B), dimension(:), pointer :: regionalmodeloffset => null() !< the new offset to compactify the range of indices
    integer(I4B), pointer :: indexcount => null() !< counts the number of cells in the interface
 
    type(connectionstype), pointer :: connections => null() !< sparse matrix with the connections
    integer(I4B), dimension(:), pointer :: connectionmask => null() !< to mask out connections from the amat coefficient calculation
 
    type(interfacemaptype), pointer :: interfacemap => null() !< defining map for the interface
 
  contains
 
    ! public
    procedure, pass(this) :: construct
    procedure, private, pass(this) :: allocatescalars
    procedure, pass(this) :: destroy
    procedure, pass(this) :: addtoregionalmodels
    procedure, pass(this) :: connectprimaryexchange
    procedure, pass(this) :: extendconnection
    procedure, pass(this) :: getdiscretization
    procedure, pass(this) :: buildinterfacemap
 
    ! 'protected'
    procedure, pass(this) :: isperiodic
 
    ! private routines
    procedure, private, pass(this) :: connectcell
    procedure, private, pass(this) :: buildconnections
    procedure, private, pass(this) :: addneighbors
    procedure, private, pass(this) :: addneighborcell
    procedure, private, pass(this) :: addremoteneighbors
    procedure, private, pass(this) :: get_regional_offset
    generic, private :: getinterfaceindex => getinterfaceindexbycell, &
      getinterfaceindexbyindexmodel
    procedure, private, pass(this) :: getinterfaceindexbycell
    procedure, private, pass(this) :: getinterfaceindexbyindexmodel
    procedure, private, pass(this) :: registerinterfacecells
    procedure, private, pass(this) :: addtoglobalmap
    procedure, private, pass(this) :: compressglobalmap
    procedure, private, pass(this) :: sortinterfacegrid
    procedure, private, pass(this) :: makeprimaryconnections
    procedure, private, pass(this) :: connectneighborcells
    procedure, private, pass(this) :: fillconnectiondatainternal
    procedure, private, pass(this) :: fillconnectiondatafromexchanges
    procedure, private, pass(this) :: createconnectionmask
    procedure, private, pass(this) :: maskinternalconnections
    procedure, private, pass(this) :: setmaskonconnection
  end type
 
contains
 
  !> @brief Construct the GridConnection and allocate
  !! the data structures for the primary connections
  !<
  subroutine construct(this, model, nrOfPrimaries, connectionName)
    class(gridconnectiontype), intent(inout) :: this !> this instance
    class(numericalmodeltype), pointer, intent(in) :: model !> the model for which the interface is constructed
    integer(I4B) :: nrOfPrimaries !> the number of primary connections between the two models
    character(len=*) :: connectionName !> the name, for memory management mostly
    ! local
    class(virtualmodeltype), pointer :: v_model
 
    this%model => model
    this%memoryPath = create_mem_path(connectionname, 'GC')
 
    call this%allocateScalars()
 
    allocate (this%boundaryCells(nrofprimaries))
    allocate (this%connectedCells(nrofprimaries))
    allocate (this%idxToGlobal(2 * nrofprimaries))
 
    v_model => get_virtual_model(model%id)
    call this%addToRegionalModels(v_model)
 
    this%nrOfBoundaryCells = 0
 
    this%internalStencilDepth = 1
    this%exchangeStencilDepth = 1
    this%haloExchanges => null()
 
  end subroutine construct
 
  !> @brief Make connections for the primary exchange
  !<
  subroutine connectprimaryexchange(this, primEx)
    class(gridconnectiontype) :: this !< this grid connection
    class(disconnexchangetype), pointer :: primEx !< the primary exchange for this connection
    ! local
    integer(I4B) :: iconn
 
    ! store the primary exchange
    this%primaryExchange => primex
 
    ! connect the cells
    do iconn = 1, primex%nexg
      call this%connectCell(primex%nodem1(iconn), primex%v_model1, &
                            primex%nodem2(iconn), primex%v_model2)
    end do
 
  end subroutine connectprimaryexchange
 
  !> @brief Connect neighboring cells at the interface by
  !! storing them in the boundary cell and connected cell
  !! arrays
  !<
  subroutine connectcell(this, idx1, v_model1, idx2, v_model2)
    class(gridconnectiontype), intent(in) :: this !< this grid connection
    integer(I4B) :: idx1 !< local index cell 1
    class(virtualmodeltype), pointer :: v_model1 !< model of cell 1
    integer(I4B) :: idx2 !< local index cell 2
    class(virtualmodeltype), pointer :: v_model2 !< model of cell 2
    ! local
    type(globalcelltype), pointer :: bnd_cell, conn_cell
 
    this%nrOfBoundaryCells = this%nrOfBoundaryCells + 1
    if (this%nrOfBoundaryCells > size(this%boundaryCells)) then
      write (*, *) 'Error: nr of cell connections exceeds '// &
        'capacity in grid connection, terminating...'
      call ustop()
    end if
 
    bnd_cell => this%boundaryCells(this%nrOfBoundaryCells)%cell
    conn_cell => this%connectedCells(this%nrOfBoundaryCells)%cell
    if (v_model1 == this%model) then
      bnd_cell%index = idx1
      bnd_cell%v_model => v_model1
      conn_cell%index = idx2
      conn_cell%v_model => v_model2
    else if (v_model2 == this%model) then
      bnd_cell%index = idx2
      bnd_cell%v_model => v_model2
      conn_cell%index = idx1
      conn_cell%v_model => v_model1
    else
      write (*, *) 'Error: unable to connect cells outside the model'
      call ustop()
    end if
 
  end subroutine connectcell
 
  !> @brief Add a model to a list of all regional models
  !<
  subroutine addtoregionalmodels(this, v_model)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection
    class(virtualmodeltype), pointer :: v_model !< the model to add to the region
    ! local
    class(*), pointer :: vm_obj
 
    vm_obj => v_model
    if (.not. this%regionalModels%ContainsObject(vm_obj)) then
      call this%regionalModels%Add(vm_obj)
    end if
 
  end subroutine addtoregionalmodels
 
  !> @brief Extend the connection topology to deal with
  !! higher levels of connectivity (neighbors-of-neighbors, etc.)
  !!
  !! The following steps are taken:
  !! 1. Recursively add interior neighbors (own model) up to the specified depth
  !! 2. Recursively add exterior neighbors
  !! 3. Allocate a (sparse) mapping table for the region
  !! 4. Build connection object for the interface grid, and the mask
  !<
  subroutine extendconnection(this)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection
    ! local
    integer(I4B) :: remoteDepth, localDepth
    integer(I4B) :: icell
    integer(I4B) :: imod, regionSize, offset
    class(virtualmodeltype), pointer :: v_model
    !integer(I4B), pointer :: nr_nodes
 
    ! we need (stencildepth-1) extra cells for the interior
    remotedepth = this%exchangeStencilDepth
    localdepth = 2 * this%internalStencilDepth - 1
    if (localdepth < remotedepth) then
      localdepth = remotedepth
    end if
 
    ! first add the neighbors for the interior
    ! (possibly extending into other models)
    do icell = 1, this%nrOfBoundaryCells
      call this%addNeighbors(this%boundaryCells(icell), localdepth, &
                             this%connectedCells(icell)%cell, .true.)
    end do
    ! and for the exterior
    do icell = 1, this%nrOfBoundaryCells
      call this%addNeighbors(this%connectedCells(icell), remotedepth, &
                             this%boundaryCells(icell)%cell, .false.)
    end do
 
    ! set up mapping for the region (models participating in interface model grid)
    allocate (this%regionalModelOffset(this%regionalModels%Count()))
    regionsize = 0
    offset = 0
    do imod = 1, this%regionalModels%Count()
      v_model => get_virtual_model_from_list(this%regionalModels, imod)
      regionsize = regionsize + v_model%dis_nodes%get()
      this%regionalModelOffset(imod) = offset
      offset = offset + v_model%dis_nodes%get()
    end do
    ! init to -1, meaning 'interface index was not assigned yet'
    allocate (this%region_to_iface_map(regionsize))
    this%region_to_iface_map = -1
 
    call this%buildConnections()
 
  end subroutine extendconnection
 
  !> @brief Builds a sparse matrix holding all cell connections,
  !< with new indices, and stores the mapping to the global ids
  subroutine buildconnections(this)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    ! local
    integer(I4B) :: icell, iconn
    integer(I4B), dimension(:), allocatable :: nnz
    type(sparsematrix), pointer :: sparse
    integer(I4B) :: ierror
    type(connectionstype), pointer :: conn
 
    ! Recursively generate interface cell indices, fill map to global cells,
    ! and add to region lookup table
    this%indexCount = 0
    do icell = 1, this%nrOfBoundaryCells
      call this%registerInterfaceCells(this%boundaryCells(icell))
    end do
    do icell = 1, this%nrOfBoundaryCells
      call this%registerInterfaceCells(this%connectedCells(icell))
    end do
    this%nrOfCells = this%indexCount
 
    ! compress lookup table
    call this%compressGlobalMap()
 
    ! sort interface indexes such that 'n > m' means 'n below m'
    call this%sortInterfaceGrid()
 
    ! allocate a map from interface index to global coordinates
    call mem_allocate(this%idxToGlobalIdx, this%nrOfCells, &
                      'IDXTOGLOBALIDX', this%memoryPath)
 
    ! create sparse data structure, to temporarily hold connections
    allocate (sparse)
    allocate (nnz(this%nrOfCells))
    nnz = initnrneighbors + 1
    call sparse%init(this%nrOfCells, this%nrOfCells, nnz)
 
    ! now (recursively) add connections to sparse, start with
    ! the primary connections (n-m from the exchange files)
    call this%makePrimaryConnections(sparse)
    ! then into own domain
    do icell = 1, this%nrOfBoundaryCells
      call this%connectNeighborCells(this%boundaryCells(icell), sparse)
    end do
    ! and same for the neighbors of connected cells
    do icell = 1, this%nrOfBoundaryCells
      call this%connectNeighborCells(this%connectedCells(icell), sparse)
    end do
 
    ! create connections object
    allocate (this%connections)
    conn => this%connections
    call conn%allocate_scalars(this%memoryPath)
    conn%nodes = this%nrOfCells
    conn%nja = sparse%nnz
    conn%njas = (conn%nja - conn%nodes) / 2
    call conn%allocate_arrays()
    do iconn = 1, conn%njas
      conn%anglex(iconn) = -999.
    end do
 
    ! fill connection from sparse
    call sparse%filliaja(conn%ia, conn%ja, ierror)
    if (ierror /= 0) then
      write (*, *) 'Error filling ia/ja in GridConnection: terminating...'
      call ustop()
    end if
    call fillisym(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym)
    call filljas(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym, conn%jas)
    call sparse%destroy()
 
    ! fill connection data (ihc, cl1, cl2, etc.) using data
    ! from models and exchanges
    call this%fillConnectionDataInternal()
    call this%fillConnectionDataFromExchanges()
 
    ! set the masks on connections
    call this%createConnectionMask()
 
  end subroutine buildconnections
 
  !< @brief Routine for finding neighbors-of-neighbors, recursively
  !<
  recursive subroutine addneighbors(this, cellNbrs, depth, mask, interior)
    use simmodule, only: ustop
    class(gridconnectiontype), intent(inout) :: this !< this grid connection
    type(cellwithnbrstype), intent(inout), target :: cellnbrs !< cell to add to
    integer(I4B), intent(inout) :: depth !< current depth (typically decreases in recursion)
    type(globalcelltype), optional :: mask !< mask to excluded back-and-forth connection between cells
    logical(LGP) :: interior !< when true, we are adding from the exchange back into the model
    ! local
    type(globalcelltype), pointer :: cell
    integer(I4B) :: ipos, ipos_start, ipos_end
    integer(I4B) :: nbridx, inbr
    integer(I4B) :: newdepth
 
    ! readability
    cell => cellnbrs%cell
 
    ! if depth == 1, then we are not adding neighbors but use
    ! the boundary and connected cell only
    if (depth < 2) then
      return
    end if
    newdepth = depth - 1
 
    ! find neighbors local to this cell by looping through grid connections
    ipos_start = cell%v_model%con_ia%get(cell%index) + 1
    ipos_end = cell%v_model%con_ia%get(cell%index + 1) - 1
    do ipos = ipos_start, ipos_end
      nbridx = cell%v_model%con_ja%get(ipos)
      call this%addNeighborCell(cellnbrs, nbridx, cellnbrs%cell%v_model, mask)
    end do
 
    ! add remote nbr using the data from the exchanges
    call this%addRemoteNeighbors(cellnbrs, mask)
 
    ! now find nbr-of-nbr
    do inbr = 1, cellnbrs%nrOfNbrs
 
      ! are we leaving the model through another exchange?
      if (interior .and. cellnbrs%cell%v_model == this%model) then
        if (.not. cellnbrs%neighbors(inbr)%cell%v_model == this%model) then
          ! decrement by 1, because the connection we are crossing is not
          ! calculated by this interface
          newdepth = newdepth - 1
        end if
      end if
      ! and add neighbors with the new depth
      call this%addNeighbors(cellnbrs%neighbors(inbr), newdepth, &
                             cellnbrs%cell, interior)
    end do
 
  end subroutine addneighbors
 
  !> @brief Add cell neighbors across models using the stored exchange
  !! data structures
  subroutine addremoteneighbors(this, cellNbrs, mask)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(cellwithnbrstype), intent(inout) :: cellNbrs !< cell to add to
    type(globalcelltype), optional :: mask !< a mask to exclude back-and-forth connections
    ! local
    integer(I4B) :: ix, iexg
    class(virtualexchangetype), pointer :: v_exchange
    class(virtualmodeltype), pointer :: v_m1, v_m2
 
    ! loop over all exchanges
    do ix = 1, this%haloExchanges%size
 
      v_exchange => get_virtual_exchange(this%haloExchanges%at(ix))
      v_m1 => v_exchange%v_model1
      v_m2 => v_exchange%v_model2
 
      ! loop over n-m links in the exchange
      if (cellnbrs%cell%v_model == v_m1) then
        do iexg = 1, v_exchange%nexg%get()
          if (v_exchange%nodem1%get(iexg) == cellnbrs%cell%index) then
            ! we have a link, now add foreign neighbor
            call this%addNeighborCell( &
              cellnbrs, v_exchange%nodem2%get(iexg), v_m2, mask)
          end if
        end do
      end if
      ! and the reverse
      if (cellnbrs%cell%v_model == v_m2) then
        do iexg = 1, v_exchange%nexg%get()
          if (v_exchange%nodem2%get(iexg) == cellnbrs%cell%index) then
            ! we have a link, now add foreign neighbor
            call this%addNeighborCell( &
              cellnbrs, v_exchange%nodem1%get(iexg), v_m1, mask)
          end if
        end do
      end if
 
    end do
 
  end subroutine addremoteneighbors
 
  !> @brief Add neighboring cell to tree structure
  !<
  subroutine addneighborcell(this, cellNbrs, newNbrIdx, v_nbr_model, mask)
    class(gridconnectiontype), intent(in) :: this !< this grid connection instance
    type(cellwithnbrstype), intent(inout) :: cellNbrs !< the root cell which to add to
    integer(I4B), intent(in) :: newNbrIdx !< the neighboring cell's index
    class(virtualmodeltype), pointer :: v_nbr_model !< the model where the new neighbor lives
    type(globalcelltype), optional :: mask !< don't add connections to this cell (optional)
 
    if (present(mask)) then
      if (newnbridx == mask%index .and. mask%v_model == v_nbr_model) then
        return
      end if
    end if
 
    call cellnbrs%addNbrCell(newnbridx, v_nbr_model)
 
  end subroutine addneighborcell
 
  !> @brief Recursively set interface cell indexes and
  !< add to the region-to-interface lookup table
  recursive subroutine registerinterfacecells(this, cellWithNbrs)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(cellwithnbrstype) :: cellwithnbrs !< the cell from where to start registering neighbors
    ! local
    integer(I4B) :: offset, inbr
    integer(I4B) :: regionidx ! unique idx in the region (all connected models)
    integer(I4B) :: ifaceidx ! unique idx in the interface grid
 
    offset = this%get_regional_offset(cellwithnbrs%cell%v_model)
    regionidx = offset + cellwithnbrs%cell%index
    ifaceidx = this%getInterfaceIndex(cellwithnbrs%cell)
    if (ifaceidx == -1) then
      this%indexCount = this%indexCount + 1
      ifaceidx = this%indexCount
      call this%addToGlobalMap(ifaceidx, cellwithnbrs%cell)
      this%region_to_iface_map(regionidx) = ifaceidx
    end if
 
    ! and also for its neighbors
    do inbr = 1, cellwithnbrs%nrOfNbrs
      call this%registerInterfaceCells(cellwithnbrs%neighbors(inbr))
    end do
 
  end subroutine registerinterfacecells
 
  !> @brief Add entry to lookup table, inflating when necessary
  !<
  subroutine addtoglobalmap(this, ifaceIdx, cell)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    integer(I4B), intent(in) :: ifaceIdx !< unique idx in the interface grid
    type(globalcelltype), intent(in) :: cell !< the global cell
    ! local
    integer(I4B) :: i, currentSize, newSize
    type(globalcelltype), dimension(:), pointer :: tempMap
 
    ! inflate?
    currentsize = size(this%idxToGlobal)
    if (ifaceidx > currentsize) then
      newsize = nint(1.5 * currentsize)
      allocate (tempmap(newsize))
      do i = 1, currentsize
        tempmap(i) = this%idxToGlobal(i)
      end do
 
      deallocate (this%idxToGlobal)
      this%idxToGlobal => tempmap
    end if
 
    this%idxToGlobal(ifaceidx) = cell
 
  end subroutine addtoglobalmap
 
  !> @brief Compress lookup table to get rid of unused entries
  !<
  subroutine compressglobalmap(this)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    ! local
    type(globalcelltype), dimension(:), pointer :: tempMap
 
    if (size(this%idxToGlobal) > this%nrOfCells) then
      allocate (tempmap(this%nrOfCells))
      tempmap(1:this%nrOfCells) = this%idxToGlobal(1:this%nrOfCells)
      deallocate (this%idxToGlobal)
      allocate (this%idxToGlobal(this%nrOfCells))
      this%idxToGlobal(1:this%nrOfCells) = tempmap(1:this%nrOfCells)
      deallocate (tempmap)
    end if
 
  end subroutine compressglobalmap
 
  !> @brief Soft cell ids in the interface grid such that
  !< id_1 < id_2 means that cell 1 lies above cell 2
  subroutine sortinterfacegrid(this)
    use gridsorting, only: quicksortgrid
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    ! local
    integer(I4B), dimension(:), allocatable :: newToOldIdx
    integer(I4B), dimension(:), allocatable :: oldToNewIdx
    integer(I4B) :: idxOld
    integer(I4B) :: i
    type(globalcelltype), dimension(:), allocatable :: sortedGlobalMap
    integer(I4B), dimension(:), allocatable :: sortedRegionMap
 
    ! sort based on coordinates
    newtooldidx = (/(i, i=1, size(this%idxToGlobal))/)
    call quicksortgrid(newtooldidx, size(newtooldidx), this%idxToGlobal)
 
    ! and invert
    allocate (oldtonewidx(size(newtooldidx)))
    do i = 1, size(oldtonewidx)
      oldtonewidx(newtooldidx(i)) = i
    end do
 
    ! reorder global table
    allocate (sortedglobalmap(size(this%idxToGlobal)))
    do i = 1, size(newtooldidx)
      sortedglobalmap(i) = this%idxToGlobal(newtooldidx(i))
    end do
    do i = 1, size(newtooldidx)
      this%idxToGlobal(i) = sortedglobalmap(i)
    end do
    deallocate (sortedglobalmap)
 
    ! reorder regional lookup table
    allocate (sortedregionmap(size(this%region_to_iface_map)))
    do i = 1, size(sortedregionmap)
      if (this%region_to_iface_map(i) /= -1) then
        idxold = this%region_to_iface_map(i)
        sortedregionmap(i) = oldtonewidx(idxold)
      else
        sortedregionmap(i) = -1
      end if
    end do
    do i = 1, size(sortedregionmap)
      this%region_to_iface_map(i) = sortedregionmap(i)
    end do
    deallocate (sortedregionmap)
 
  end subroutine sortinterfacegrid
 
  !> @brief Add primary connections to the sparse data structure
  !<
  subroutine makeprimaryconnections(this, sparse)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(sparsematrix), pointer :: sparse !< the sparse data structure to hold the connections
    ! local
    integer(I4B) :: icell
    integer(I4B) :: ifaceIdx, ifaceIdxNbr
 
    do icell = 1, this%nrOfBoundaryCells
      ifaceidx = this%getInterfaceIndex(this%boundaryCells(icell)%cell)
      ifaceidxnbr = this%getInterfaceIndex(this%connectedCells(icell)%cell)
 
      ! add diagonals to sparse
      call sparse%addconnection(ifaceidx, ifaceidx, 1)
      call sparse%addconnection(ifaceidxnbr, ifaceidxnbr, 1)
 
      ! and cross terms
      call sparse%addconnection(ifaceidx, ifaceidxnbr, 1)
      call sparse%addconnection(ifaceidxnbr, ifaceidx, 1)
    end do
 
  end subroutine makeprimaryconnections
 
  !> @brief Recursively add higher order connections (from
  !! cells neighboring the primarily connected cells) to the
  !< sparse data structure
  recursive subroutine connectneighborcells(this, cell, sparse)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(cellwithnbrstype) :: cell !< the cell whose connections is to be added
    type(sparsematrix), pointer :: sparse !< the sparse data structure to hold the connections
    ! local
    integer(I4B) :: ifaceidx, ifaceidxnbr ! unique idx in the interface grid
    integer(I4B) :: inbr
 
    ifaceidx = this%getInterfaceIndex(cell%cell)
    do inbr = 1, cell%nrOfNbrs
      ifaceidxnbr = this%getInterfaceIndex(cell%neighbors(inbr)%cell)
 
      call sparse%addconnection(ifaceidxnbr, ifaceidxnbr, 1)
      call sparse%addconnection(ifaceidx, ifaceidxnbr, 1)
      call sparse%addconnection(ifaceidxnbr, ifaceidx, 1)
 
      ! recurse
      call this%connectNeighborCells(cell%neighbors(inbr), sparse)
    end do
 
  end subroutine connectneighborcells
 
  !> @brief Fill connection data (ihc, cl1, ...) for
  !< connections between cells within the same model.
  subroutine fillconnectiondatainternal(this)
    use constantsmodule, only: dpi, dtwopi
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    ! local
    type(connectionstype), pointer :: conn
    integer(I4B) :: n, m, ipos, isym, ipos_orig, isym_orig
    type(globalcelltype), pointer :: ncell, mcell
    class(virtualmodeltype), pointer :: v_m !< pointer to virtual model (readability)
 
    conn => this%connections
 
    do n = 1, conn%nodes
      do ipos = conn%ia(n) + 1, conn%ia(n + 1) - 1
        m = conn%ja(ipos)
        if (n > m) cycle
 
        isym = conn%jas(ipos)
        ncell => this%idxToGlobal(n)
        mcell => this%idxToGlobal(m)
        if (ncell%v_model == mcell%v_model) then
 
          ! for readability
          v_m => ncell%v_model
 
          ! within same model, straight copy
          ipos_orig = getcsrindex(ncell%index, mcell%index, &
                                  v_m%con_ia%get_array(), v_m%con_ja%get_array())
 
          if (ipos_orig == 0) then
            ! periodic boundary conditions can add connections between cells in
            ! the same model, but they are dealt with through the exchange data
            if (this%isPeriodic(ncell%index, mcell%index)) cycle
 
            ! this should not be possible
            write (*, *) 'Error: cannot find cell connection in model grid'
            call ustop()
          end if
 
          isym_orig = v_m%con_jas%get(ipos_orig)
          conn%hwva(isym) = v_m%con_hwva%get(isym_orig)
          conn%ihc(isym) = v_m%con_ihc%get(isym_orig)
          if (ncell%index < mcell%index) then
            conn%cl1(isym) = v_m%con_cl1%get(isym_orig)
            conn%cl2(isym) = v_m%con_cl2%get(isym_orig)
            conn%anglex(isym) = v_m%con_anglex%get(isym_orig)
          else
            conn%cl1(isym) = v_m%con_cl2%get(isym_orig)
            conn%cl2(isym) = v_m%con_cl1%get(isym_orig)
            conn%anglex(isym) = mod(v_m%con_anglex%get(isym_orig) + dpi, dtwopi)
          end if
        end if
      end do
    end do
 
  end subroutine fillconnectiondatainternal
 
  !> @brief Fill connection data (ihc, cl1, ...) for
  !< all exchanges
  subroutine fillconnectiondatafromexchanges(this)
    use constantsmodule, only: dpio180
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    ! local
    integer(I4B) :: inx, iexg
    integer(I4B) :: ipos, isym
    integer(I4B) :: nOffset, mOffset, nIfaceIdx, mIfaceIdx
    class(virtualexchangetype), pointer :: v_exg
    class(virtualmodeltype), pointer :: v_m1, v_m2
    type(connectionstype), pointer :: conn
 
    conn => this%connections
 
    do inx = 1, this%haloExchanges%size
      v_exg => get_virtual_exchange(this%haloExchanges%at(inx))
 
      v_m1 => v_exg%v_model1
      v_m2 => v_exg%v_model2
 
      if (v_exg%ianglex%get() > 0) then
        conn%ianglex = 1
      end if
 
      noffset = this%get_regional_offset(v_m1)
      moffset = this%get_regional_offset(v_m2)
      do iexg = 1, v_exg%nexg%get()
        nifaceidx = this%region_to_iface_map(noffset + v_exg%nodem1%get(iexg))
        mifaceidx = this%region_to_iface_map(moffset + v_exg%nodem2%get(iexg))
        ! not all nodes from the exchanges are part of the interface grid
        ! (think of exchanges between neighboring models, and their neighbors)
        if (nifaceidx == -1 .or. mifaceidx == -1) then
          cycle
        end if
 
        ipos = conn%getjaindex(nifaceidx, mifaceidx)
        ! (see prev. remark) sometimes the cells are in the interface grid,
        ! but the connection isn't. This can happen for leaf nodes of the grid.
        if (ipos == 0) then
          ! no match, safely cycle
          cycle
        end if
        isym = conn%jas(ipos)
 
        ! note: cl1 equals L_nm: the length from cell n to the shared
        ! face with cell m (and cl2 analogously for L_mn)
        if (nifaceidx < mifaceidx) then
          conn%cl1(isym) = v_exg%cl1%get(iexg)
          conn%cl2(isym) = v_exg%cl2%get(iexg)
          if (v_exg%ianglex%get() > 0) then
            conn%anglex(isym) = &
              v_exg%auxvar%get(v_exg%ianglex%get(), iexg) * dpio180
          end if
        else
          conn%cl1(isym) = v_exg%cl2%get(iexg)
          conn%cl2(isym) = v_exg%cl1%get(iexg)
          if (v_exg%ianglex%get() > 0) then
            conn%anglex(isym) = mod(v_exg%auxvar%get(v_exg%ianglex%get(), iexg) &
                                    + 180.0_dp, 360.0_dp) * dpio180
          end if
        end if
        conn%hwva(isym) = v_exg%hwva%get(iexg)
        conn%ihc(isym) = v_exg%ihc%get(iexg)
 
      end do
    end do
 
  end subroutine fillconnectiondatafromexchanges
 
  !> @brief Create the connection masks
  !!
  !! The level indicates the nr of connections away from
  !! the remote neighbor, the diagonal term holds the negated
  !! value of their nearest connection. We end with setting
  !< a normalized mask to the connections object: 0 or 1
  subroutine createconnectionmask(this)
    class(gridconnectiontype), intent(inout) :: this !< instance of this grid connection
    ! local
    integer(I4B) :: icell, inbr, n, ipos
    integer(I4B) :: level, newMask
    type(cellwithnbrstype), pointer :: cell, nbrCell
 
    ! set all masks to zero to begin with
    do ipos = 1, this%connections%nja
      call this%connections%set_mask(ipos, 0)
    end do
 
    ! remote connections remain masked
    ! now set mask for exchange connections (level == 1)
    level = 1
    do icell = 1, this%nrOfBoundaryCells
      call this%setMaskOnConnection(this%boundaryCells(icell), &
                                    this%connectedCells(icell), level)
      ! for cross-boundary connections, we need to apply the mask to both n-m and m-n,
      ! because if the upper triangular one is disabled, its transposed (lower triangular)
      ! counter part is skipped in the NPF calculation as well.
      call this%setMaskOnConnection(this%connectedCells(icell), &
                                    this%boundaryCells(icell), level)
    end do
 
    ! now extend mask recursively into the internal domain (level > 1)
    do icell = 1, this%nrOfBoundaryCells
      cell => this%boundaryCells(icell)
      do inbr = 1, cell%nrOfNbrs
        nbrcell => this%boundaryCells(icell)%neighbors(inbr)
        level = 2 ! this is incremented within the recursion
        call this%maskInternalConnections(this%boundaryCells(icell), &
                                          this%boundaryCells(icell)% &
                                          neighbors(inbr), level)
      end do
    end do
 
    ! set normalized mask:
    ! =1 for links with connectivity <= interior stencil depth
    ! =0 otherwise
    do n = 1, this%connections%nodes
      ! set diagonals to zero
      call this%connections%set_mask(this%connections%ia(n), 0)
 
      do ipos = this%connections%ia(n) + 1, this%connections%ia(n + 1) - 1
        newmask = 0
        if (this%connections%mask(ipos) > 0) then
          if (this%connections%mask(ipos) < this%internalStencilDepth + 1) then
            newmask = 1
          end if
        end if
        ! set mask on off-diag
        call this%connections%set_mask(ipos, newmask)
      end do
    end do
 
  end subroutine createconnectionmask
 
  !> @brief Recursively mask connections, increasing the level as we go
  !<
  recursive subroutine maskinternalconnections(this, cell, nbrCell, level)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(cellwithnbrstype), intent(inout) :: cell !< cell 1 in the connection to mask
    type(cellwithnbrstype), intent(inout) :: nbrcell !< cell 2 in the connection to mask
    integer(I4B), intent(in) :: level
    ! local
    integer(I4B) :: inbr, newlevel
 
    ! only set the mask for internal connections, leaving the
    ! others at 0
    if (cell%cell%v_model == this%model .and. &
        nbrcell%cell%v_model == this%model) then
      ! this will set a mask on both diagonal, and both cross terms
      call this%setMaskOnConnection(cell, nbrcell, level)
      call this%setMaskOnConnection(nbrcell, cell, level)
    end if
 
    ! recurse on nbrs-of-nbrs
    newlevel = level + 1
    do inbr = 1, nbrcell%nrOfNbrs
      call this%maskInternalConnections(nbrcell, &
                                        nbrcell%neighbors(inbr), &
                                        newlevel)
    end do
 
  end subroutine maskinternalconnections
 
  !> @brief Set a mask on the connection from a cell to its neighbor,
  !! (and not the transposed!) not overwriting the current level
  !< of a connection when it is smaller
  subroutine setmaskonconnection(this, cell, nbrCell, level)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(cellwithnbrstype), intent(inout) :: cell !< cell 1 in the connection
    type(cellwithnbrstype), intent(inout) :: nbrCell !< cell 2 in the connection
    integer(I4B), intent(in) :: level !< the level value to set the mask to
    ! local
    integer(I4B) :: ifaceIdx, ifaceIdxNbr
    integer(I4B) :: iposdiag, ipos
    integer(I4B) :: currentLevel
 
    ifaceidx = this%getInterfaceIndex(cell%cell)
    ifaceidxnbr = this%getInterfaceIndex(nbrcell%cell)
 
    ! diagonal
    iposdiag = this%connections%getjaindex(ifaceidx, ifaceidx)
    currentlevel = this%connections%mask(iposdiag)
    if (currentlevel == 0 .or. level < currentlevel) then
      call this%connections%set_mask(iposdiag, level)
    end if
    ! cross term
    ipos = this%connections%getjaindex(ifaceidx, ifaceidxnbr)
    currentlevel = this%connections%mask(ipos)
    if (currentlevel == 0 .or. level < currentlevel) then
      call this%connections%set_mask(ipos, level)
    end if
 
  end subroutine setmaskonconnection
 
  !> @brief Get interface index from global cell
  !<
  function getinterfaceindexbycell(this, cell) result(iface_idx)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    type(globalcelltype), intent(in) :: cell !< the global cell to get the interface index for
    integer(I4B) :: iface_idx !< the index in the interface model
    ! local
    integer(I4B) :: offset, region_idx
 
    offset = this%get_regional_offset(cell%v_model)
    region_idx = offset + cell%index
    iface_idx = this%region_to_iface_map(region_idx)
 
  end function getinterfaceindexbycell
 
  !> @brief Get interface index from a model pointer and the local index
  !<
  function getinterfaceindexbyindexmodel(this, index, v_model) result(iface_idx)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    integer(I4B) :: index !< the local cell index
    class(virtualmodeltype), pointer :: v_model !< the cell's model
    integer(I4B) :: iface_idx !< the index in the interface model
    ! local
    integer(I4B) :: offset, region_idx
 
    offset = this%get_regional_offset(v_model)
    region_idx = offset + index
    iface_idx = this%region_to_iface_map(region_idx)
 
  end function getinterfaceindexbyindexmodel
 
  !> @brief Get the offset for a regional model
  !<
  function get_regional_offset(this, v_model) result(offset)
    class(gridconnectiontype), intent(inout) :: this !< this grid connection instance
    class(virtualmodeltype), pointer :: v_model !< the model to get the offset for
    integer(I4B) :: offset !< the index offset in the regional domain
    ! local
    integer(I4B) :: im
    class(virtualmodeltype), pointer :: vm
 
    offset = 0
    do im = 1, this%regionalModels%Count()
      vm => get_virtual_model_from_list(this%regionalModels, im)
      if (vm == v_model) then
        offset = this%regionalModelOffset(im)
        return
      end if
    end do
 
  end function get_regional_offset
 
  !> @brief Allocate scalar data
  !<
  subroutine allocatescalars(this)
    use memorymanagermodule, only: mem_allocate
    class(gridconnectiontype) :: this !< this grid connection instance
 
    call mem_allocate(this%nrOfBoundaryCells, 'NRBNDCELLS', this%memoryPath)
    call mem_allocate(this%indexCount, 'IDXCOUNT', this%memoryPath)
    call mem_allocate(this%nrOfCells, 'NRCELLS', this%memoryPath)
 
  end subroutine allocatescalars
 
  !> @brief Sets the discretization (DISU) after all
  !! preprocessing by this grid connection has been done,
  !< this comes after disu_cr
  subroutine getdiscretization(this, disu)
    use connectionsmodule
    use sparsemodule, only: sparsematrix
    class(gridconnectiontype) :: this !< the grid connection
    class(disutype), pointer :: disu !< the target disu object
    ! local
    integer(I4B) :: icell, nrOfCells, idx
    class(virtualmodeltype), pointer :: v_model
    real(DP) :: xglo, yglo
 
    ! the following is similar to dis_df
    nrofcells = this%nrOfCells
    disu%nodes = nrofcells
    disu%nodesuser = nrofcells
    disu%nja = this%connections%nja
 
    call disu%allocate_arrays()
    ! these are otherwise allocated in dis%read_dimensions
    call disu%allocate_arrays_mem()
 
    ! fill data
    do icell = 1, nrofcells
      idx = this%idxToGlobal(icell)%index
      disu%top(icell) = -huge(1.0_dp)
      disu%bot(icell) = -huge(1.0_dp)
      disu%area(icell) = -huge(1.0_dp)
    end do
 
    ! grid connections follow from GridConnection:
    disu%con => this%connections
    disu%njas = disu%con%njas
 
    ! copy cell x,y
    do icell = 1, nrofcells
      idx = this%idxToGlobal(icell)%index
      v_model => this%idxToGlobal(icell)%v_model
 
      ! we are merging grids with possibly (likely) different origins,
      ! transform to global coordinates:
      call dis_transform_xy(v_model%dis_xc%get(idx), &
                            v_model%dis_yc%get(idx), &
                            v_model%dis_xorigin%get(), &
                            v_model%dis_yorigin%get(), &
                            v_model%dis_angrot%get(), &
                            xglo, yglo)
 
      ! NB: usernodes equals internal nodes for interface
      disu%cellxy(1, icell) = xglo
      disu%xc(icell) = xglo
      disu%cellxy(2, icell) = yglo
      disu%yc(icell) = yglo
    end do
 
    ! if vertices will be needed, it will look like this:
    !
    ! 1. determine total nr. of verts
    ! 2. allocate vertices list
    ! 3. create sparse
    ! 4. get vertex data per cell, add functions to base
    ! 5. add vertex (x,y) to list and connectivity to sparse
    ! 6. generate ia/ja from sparse
 
  end subroutine getdiscretization
 
  !> @brief Build interface map object for outside use
  subroutine buildinterfacemap(this)
    use basemodelmodule, only: basemodeltype, getbasemodelfromlist
    use stlvecintmodule
    class(gridconnectiontype) :: this !< this grid connection
    ! local
    integer(I4B) :: i, j, iloc, jloc
    integer(I4B) :: im, ix, mid, n
    integer(I4B) :: ipos, ipos_model
    type(stlvecint) :: model_ids
    type(stlvecint) :: src_idx_tmp, tgt_idx_tmp, sign_tmp
    class(virtualexchangetype), pointer :: v_exg
    class(virtualmodeltype), pointer :: vm
    class(virtualmodeltype), pointer :: v_model1, v_model2
    type(interfacemaptype), pointer :: imap
    integer(I4B), dimension(:), pointer, contiguous :: ia_ptr, ja_ptr
 
    allocate (this%interfaceMap)
    imap => this%interfaceMap
 
    ! first get the participating models
    call model_ids%init()
    do i = 1, this%nrOfCells
      call model_ids%push_back_unique(this%idxToGlobal(i)%v_model%id)
    end do
 
    ! initialize the map
    call imap%init(model_ids%size, this%haloExchanges%size)
 
    ! for each model part of this interface, ...
    do im = 1, model_ids%size
      mid = model_ids%at(im)
      imap%model_ids(im) = mid
      vm => get_virtual_model(mid)
      imap%model_names(im) = vm%name
      call src_idx_tmp%init()
      call tgt_idx_tmp%init()
 
      ! store the node map for this model
      do i = 1, this%nrOfCells
        if (mid == this%idxToGlobal(i)%v_model%id) then
          call src_idx_tmp%push_back(this%idxToGlobal(i)%index)
          call tgt_idx_tmp%push_back(i)
        end if
      end do
 
      ! and copy into interface map
      allocate (imap%node_maps(im)%src_idx(src_idx_tmp%size))
      allocate (imap%node_maps(im)%tgt_idx(tgt_idx_tmp%size))
      do i = 1, src_idx_tmp%size
        imap%node_maps(im)%src_idx(i) = src_idx_tmp%at(i)
        imap%node_maps(im)%tgt_idx(i) = tgt_idx_tmp%at(i)
      end do
 
      call src_idx_tmp%destroy()
      call tgt_idx_tmp%destroy()
 
      ! and for connections
      call src_idx_tmp%init()
      call tgt_idx_tmp%init()
 
      ! store the connection map for this model
      do i = 1, this%nrOfCells
        if (mid /= this%idxToGlobal(i)%v_model%id) cycle
        do ipos = this%connections%ia(i), this%connections%ia(i + 1) - 1
          j = this%connections%ja(ipos)
          if (mid /= this%idxToGlobal(j)%v_model%id) cycle
 
          ! i and j are now in same model (mid)
          iloc = this%idxToGlobal(i)%index
          jloc = this%idxToGlobal(j)%index
          ia_ptr => this%idxToGlobal(i)%v_model%con_ia%get_array()
          ja_ptr => this%idxToGlobal(i)%v_model%con_ja%get_array()
          ipos_model = getcsrindex(iloc, jloc, ia_ptr, ja_ptr)
          call src_idx_tmp%push_back(ipos_model)
          call tgt_idx_tmp%push_back(ipos)
        end do
      end do
 
      ! copy into interface map
      allocate (imap%conn_maps(im)%src_idx(src_idx_tmp%size))
      allocate (imap%conn_maps(im)%tgt_idx(tgt_idx_tmp%size))
      do i = 1, src_idx_tmp%size
        imap%conn_maps(im)%src_idx(i) = src_idx_tmp%at(i)
        imap%conn_maps(im)%tgt_idx(i) = tgt_idx_tmp%at(i)
      end do
 
      call src_idx_tmp%destroy()
      call tgt_idx_tmp%destroy()
 
    end do
 
    call model_ids%destroy()
 
    ! for each exchange that is part of this interface
    do ix = 1, this%haloExchanges%size
 
      ! all exchanges in this list should have at
      ! least one relevant connection for this map
      v_exg => get_virtual_exchange(this%haloExchanges%at(ix))
      v_model1 => v_exg%v_model1
      v_model2 => v_exg%v_model2
 
      imap%exchange_ids(ix) = v_exg%id
      imap%exchange_names(ix) = v_exg%name
 
      call src_idx_tmp%init()
      call tgt_idx_tmp%init()
      call sign_tmp%init()
 
      do n = 1, v_exg%nexg%get()
        i = this%getInterfaceIndex(v_exg%nodem1%get(n), v_model1)
        j = this%getInterfaceIndex(v_exg%nodem2%get(n), v_model2)
        if (i == -1 .or. j == -1) cycle ! not all exchange nodes are part of the interface
        ipos = this%connections%getjaindex(i, j)
        if (ipos == 0) then
          ! this can typically happen at corner points for a larger
          ! stencil (XT3D), when both i and j are included in the
          ! interface grid as leaf nodes, but their connection is not.
          ! (c.f. 'test_gwf_ifmod_mult_exg.py')
          cycle
        end if
        call src_idx_tmp%push_back(n)
        call tgt_idx_tmp%push_back(ipos)
        call sign_tmp%push_back(1)
 
        ! and the reverse connection:
        call src_idx_tmp%push_back(n)
        call tgt_idx_tmp%push_back(this%connections%isym(ipos))
        call sign_tmp%push_back(-1)
      end do
 
      allocate (imap%exchange_maps(ix)%src_idx(src_idx_tmp%size))
      allocate (imap%exchange_maps(ix)%tgt_idx(tgt_idx_tmp%size))
      allocate (imap%exchange_maps(ix)%sign(sign_tmp%size))
      do i = 1, src_idx_tmp%size
        imap%exchange_maps(ix)%src_idx(i) = src_idx_tmp%at(i)
        imap%exchange_maps(ix)%tgt_idx(i) = tgt_idx_tmp%at(i)
        imap%exchange_maps(ix)%sign(i) = sign_tmp%at(i)
      end do
 
      call src_idx_tmp%destroy()
      call tgt_idx_tmp%destroy()
      call sign_tmp%destroy()
 
    end do
 
    ! set the primary exchange idx
    ! findloc cannot be used until gfortran 9...
    imap%prim_exg_idx = -1
    do i = 1, imap%nr_exchanges
      if (imap%exchange_names(i) == this%primaryExchange%name) then
        imap%prim_exg_idx = i
        exit
      end if
    end do
 
    ! sanity check
    ! do i = 1, interfaceMap%nr_models
    !   if (size(interfaceMap%node_map(i)%src_idx) == 0) then
    !     write(*,*) 'Error: empty node map in interface for ', &
    !                this%primaryExchange%name
    !     call ustop()
    !   end if
    !   if (size(interfaceMap%connection_map(i)%src_idx) == 0) then
    !     write(*,*) 'Error: empty connection map in interface for ', &
    !                this%primaryExchange%name
    !     call ustop()
    !   end if
    ! end do
    ! do i = 1, interfaceMap%nr_exchanges
    !   if (size(interfaceMap%exchange_map(i)%src_idx) == 0) then
    !     write(*,*) 'Error: empty exchange map in interface for ', &
    !                this%primaryExchange%name
    !     call ustop()
    !   end if
    ! end do
 
  end subroutine buildinterfacemap
 
  !> @brief Deallocate grid connection resources
  !<
  subroutine destroy(this)
    use memorymanagermodule, only: mem_deallocate
    class(gridconnectiontype) :: this !< this grid connection instance
 
    call mem_deallocate(this%nrOfBoundaryCells)
    call mem_deallocate(this%indexCount)
    call mem_deallocate(this%nrOfCells)
 
    ! arrays
    deallocate (this%idxToGlobal)
    deallocate (this%boundaryCells)
    deallocate (this%connectedCells)
 
    call mem_deallocate(this%idxToGlobalIdx)
 
  end subroutine destroy
 
  !> @brief Test if the connection between nodes within
  !< the same model is periodic
  function isperiodic(this, n, m) result(periodic)
    class(gridconnectiontype), intent(in) :: this !< this grid connection instance
    integer(I4B), intent(in) :: n !< first node of the connection
    integer(I4B), intent(in) :: m !< second node of the connection
    logical :: periodic !< true when periodic
    ! local
    integer(I4B) :: icell
 
    periodic = .false.
    do icell = 1, this%nrOfBoundaryCells
      if (.not. this%boundaryCells(icell)%cell%v_model == &
          this%connectedCells(icell)%cell%v_model) then
        cycle
      end if
 
      ! one way
      if (this%boundaryCells(icell)%cell%index == n) then
        if (this%connectedCells(icell)%cell%index == m) then
          periodic = .true.
          return
        end if
      end if
      ! or the other
      if (this%boundaryCells(icell)%cell%index == m) then
        if (this%connectedCells(icell)%cell%index == n) then
          periodic = .true.
          return
        end if
      end if
 
    end do
 
  end function
 
end module gridconnectionmodule