test_pointlist_from_mesh_elems_wmask Subroutine

subroutine test_pointlist_from_mesh_elems_wmask(rc)

Arguments

Type IntentOptional Attributes Name
integer, intent(out) :: rc

Source Code

  subroutine test_pointlist_from_mesh_elems_wmask(rc)
    integer, intent(out) :: rc

    integer :: localrc

    !LOCAL VARIABLES:
    type(ESMF_PointList) :: pointlist
    type(ESMF_VM) :: vm
    integer :: maxpts, mydims, mypts, myid
    type(ESMF_Mesh) :: myMesh
    integer :: numNodes
    integer :: numTriElems, numQuadElems, numTotElems
    integer, allocatable :: elemIds(:),elemTypes(:),elemConn(:)
    integer, allocatable :: nodeIds(:), nodeOwners(:)
    real(ESMF_KIND_R8), allocatable :: nodeCoords(:)
    real(ESMF_KIND_R8), allocatable :: elemCoords(:)
    integer :: petCount,localPet
    integer, allocatable :: elemMask(:)
    integer(ESMF_KIND_I4) :: maskValues(2)
    integer :: local_pts
    real(ESMF_KIND_R8), dimension(2) :: test_coords
    real(ESMF_KIND_R8) test_coordx,test_coordy
    real(ESMF_KIND_R8) my_err1,my_err2,my_err3

    ! get global VM
    call ESMF_VMGetGlobal(vm, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    call ESMF_VMGet(vm, localPet=localPet, petCount=petCount, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! If we don't have 1 or 4 PETS then exit successfully
    if ((petCount .ne. 1) .and. (petCount .ne. 4)) then
      print*,'ERROR:  test must be run using exactly 1 or 4 PETS - detected ',petCount
      rc=ESMF_FAILURE
      return
    endif

    ! setup source Mesh
    if (petCount .eq. 1) then
      ! Set number of nodes
      numNodes=9

      ! Allocate and fill the node id array.
      allocate(nodeIds(numNodes))

      nodeIds=(/1,2,3,4,5,6,7,8,9/)

      ! Allocate and fill node coordinate array.
      ! Since this is a 2D Mesh the size is 2x the
      ! number of nodes.
      allocate(nodeCoords(2*numNodes))
      nodeCoords=(/0.0,0.0, & ! node id 1
                   1.0,0.0, & ! node id 2
                   2.0,0.0, & ! node id 3
                   0.0,1.0, & ! node id 4
                   1.0,1.0, & ! node id 5
                   2.0,1.0, & ! node id 6
                   0.0,2.0, & ! node id 7
                   1.0,2.0, & ! node id 8
                   2.0,2.0 /) ! node id 9

      allocate(nodeOwners(numNodes))
      nodeOwners=0 ! everything on PET 0

      ! Set the number of each type of element, plus the total number.
      numQuadElems=3
      numTriElems=2
      numTotElems=numQuadElems+numTriElems


      ! Allocate and fill the element id array.
      allocate(elemIds(numTotElems))
      elemIds=(/1,2,3,4,5/)

      ! Allocate and fill the element mask array.
      ! set masks on elements 1,3,5
      allocate(elemMask(numTotElems))
      elemMask=(/1,0,2,0,3/)

      ! Allocate and fill the element topology type array.
      allocate(elemTypes(numTotElems))
      elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! elem id 1
                  ESMF_MESHELEMTYPE_TRI,  & ! elem id 2
                  ESMF_MESHELEMTYPE_TRI,  & ! elem id 3
                  ESMF_MESHELEMTYPE_QUAD, & ! elem id 4
                  ESMF_MESHELEMTYPE_QUAD/)  ! elem id 5


      ! Allocate and fill the element connection type array.
      ! Note that entries in this array refer to the
      ! positions in the nodeIds, etc. arrays and that
      ! the order and number of entries for each element
      ! reflects that given in the Mesh options
      ! section for the corresponding entry
      ! in the elemTypes array. The number of
      ! entries in this elemConn array is the
      ! number of nodes in a quad. (4) times the
      ! number of quad. elements plus the number
      ! of nodes in a triangle (3) times the number
      ! of triangle elements.
      allocate(elemConn(4*numQuadElems+3*numTriElems))
      elemConn=(/1,2,5,4, &  ! elem id 1
                 2,3,5,   &  ! elem id 2
                 3,6,5,   &  ! elem id 3
                 4,5,8,7, &  ! elem id 4
                 5,6,9,8/)   ! elem id 5

      !! elem coords
      allocate(elemCoords(2*numTotElems))
      elemCoords=(/0.5,0.5, & ! 1
                   1.1,0.1, & ! 2
                   1.9,0.9, & ! 3
                   0.5,1.5, & ! 4
                   1.5,1.5/)  ! 5

      local_pts=3
      test_coordx=1.1
      test_coordy=0.1

    else if (petCount .eq. 4) then
      ! Setup mesh data depending on PET
      if (localPET .eq. 0) then !!! This part only for PET 0
        ! Set number of nodes
        numNodes=4

        ! Allocate and fill the node id array.
        allocate(nodeIds(numNodes))
        nodeIds=(/1,2,4,5/)

        ! Allocate and fill node coordinate array.
        ! Since this is a 2D Mesh the size is 2x the
        ! number of nodes.
        allocate(nodeCoords(2*numNodes))
        nodeCoords=(/0.0,0.0, & ! node id 1
                     1.0,0.0, & ! node id 2
                     0.0,1.0, & ! node id 4
                     1.0,1.0 /) ! node id 5

        ! Allocate and fill the node owner array.
        allocate(nodeOwners(numNodes))
        nodeOwners=(/0, & ! node id 1
                     0, & ! node id 2
                     0, & ! node id 4
                     0/)  ! node id 5

        ! Set the number of each type of element, plus the total number.
        numQuadElems=1
        numTriElems=0
        numTotElems=numQuadElems+numTriElems

        ! Allocate and fill the element id array.
        allocate(elemIds(numTotElems))
        elemIds=(/1/)

        ! Allocate and fill the element mask array.
        ! set masks on elements 1,3,5
        allocate(elemMask(numTotElems))
        elemMask=(/1/)

        ! Allocate and fill the element topology type array.
        allocate(elemTypes(numTotElems))
        elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 1

        ! Allocate and fill the element connection type array.
        ! Note that entry are local indices
        allocate(elemConn(4*numQuadElems+3*numTriElems))
        elemConn=(/1,2,4,3/) ! elem id 1

        !! elem coords
        allocate(elemCoords(2*numTotElems))
        elemCoords=(/0.5,0.5/)  ! 1

        local_pts=0
        test_coordx=0.5
        test_coordy=0.5

      else if (localPET .eq. 1) then !!! This part only for PET 1
        ! Set number of nodes
        numNodes=4

        ! Allocate and fill the node id array.
        allocate(nodeIds(numNodes))
        nodeIds=(/2,3,5,6/)

        ! Allocate and fill node coordinate array.
        ! Since this is a 2D Mesh the size is 2x the
        ! number of nodes.
        allocate(nodeCoords(2*numNodes))
        nodeCoords=(/1.0,0.0, & ! node id 2
                     2.0,0.0, & ! node id 3
                     1.0,1.0, & ! node id 5
                     2.0,1.0 /) ! node id 6

        ! Allocate and fill the node owner array.
        allocate(nodeOwners(numNodes))
        nodeOwners=(/0, & ! node id 2
                     1, & ! node id 3
                     0, & ! node id 5
                     1/)  ! node id 6

        ! Set the number of each type of element, plus the total number.
        numQuadElems=0
        numTriElems=2
        numTotElems=numQuadElems+numTriElems

        ! Allocate and fill the element id array.
        allocate(elemIds(numTotElems))
        elemIds=(/2,3/)

        ! Allocate and fill the element mask array.
        ! set masks on elements 1,3,5
        allocate(elemMask(numTotElems))
        elemMask=(/0,2/)

        ! Allocate and fill the element topology type array.
        allocate(elemTypes(numTotElems))
        elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! elem id 2
                    ESMF_MESHELEMTYPE_TRI/)  ! elem id 3

        ! Allocate and fill the element connection type array.
        allocate(elemConn(4*numQuadElems+3*numTriElems))
        elemConn=(/1,2,3, & ! elem id 2
                   2,4,3/)  ! elem id 3

        !! elem coords
        allocate(elemCoords(2*numTotElems))
        elemCoords=(/1.1,0.1, & ! 2
                     1.9,0.9/)  ! 3

        local_pts=1
        test_coordx=1.1
        test_coordy=0.1

     else if (localPET .eq. 2) then !!! This part only for PET 2
        ! Set number of nodes
        numNodes=4

        ! Allocate and fill the node id array.
        allocate(nodeIds(numNodes))
        nodeIds=(/4,5,7,8/)

        ! Allocate and fill node coordinate array.
        ! Since this is a 2D Mesh the size is 2x the
        ! number of nodes.
        allocate(nodeCoords(2*numNodes))
        nodeCoords=(/0.0,1.0, & ! node id 4
                     1.0,1.0, & ! node id 5
                     0.0,2.0, & ! node id 7
                     1.0,2.0 /) ! node id 8

        ! Allocate and fill the node owner array.
        ! Since this Mesh is all on PET 0, it's just set to all 0.
        allocate(nodeOwners(numNodes))
        nodeOwners=(/0, & ! node id 4
                     0, & ! node id 5
                     2, & ! node id 7
                     2/)  ! node id 8

                     ! Set the number of each type of element, plus the total number.
        numQuadElems=1
        numTriElems=0
        numTotElems=numQuadElems+numTriElems

        ! Allocate and fill the element id array.
        allocate(elemIds(numTotElems))
        elemIds=(/4/)

        ! Allocate and fill the element mask array.
        ! set masks on elements 1,3,5
        allocate(elemMask(numTotElems))
        elemMask=(/0/)

        ! Allocate and fill the element topology type array.
        allocate(elemTypes(numTotElems))
        elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 4

        ! Allocate and fill the element connection type array.
        allocate(elemConn(4*numQuadElems+3*numTriElems))
        elemConn=(/1,2,4,3/) ! elem id 4

        !! elem coords
        allocate(elemCoords(2*numTotElems))
        elemCoords=(/0.5,1.5/)  ! 4

        local_pts=1
        test_coordx=0.5
        test_coordy=1.5

      else if (localPET .eq. 3) then !!! This part only for PET 3
        ! Set number of nodes
        numNodes=4

        ! Allocate and fill the node id array.
        allocate(nodeIds(numNodes))
        nodeIds=(/5,6,8,9/)

        ! Allocate and fill node coordinate array.
        ! Since this is a 2D Mesh the size is 2x the
        ! number of nodes.
        allocate(nodeCoords(2*numNodes))
        nodeCoords=(/1.0,1.0, &  ! node id 5
                     2.0,1.0, &  ! node id 6
                     1.0,2.0, &  ! node id 8
                     2.0,2.0 /)  ! node id 9

        ! Allocate and fill the node owner array.
        allocate(nodeOwners(numNodes))
        nodeOwners=(/0, & ! node id 5
                     1, & ! node id 6
                     2, & ! node id 8
                     3/)  ! node id 9

                     ! Set the number of each type of element, plus the total number.
        numQuadElems=1
        numTriElems=0
        numTotElems=numQuadElems+numTriElems

        ! Allocate and fill the element id array.
        allocate(elemIds(numTotElems))
        elemIds=(/5/)

        ! Allocate and fill the element mask array.
        ! set masks on elements 1,3,5
        allocate(elemMask(numTotElems))
        elemMask=(/3/)

        ! Allocate and fill the element topology type array.
        allocate(elemTypes(numTotElems))
        elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 5

        ! Allocate and fill the element connection type array.
        allocate(elemConn(4*numQuadElems+3*numTriElems))
        elemConn=(/1,2,4,3/) ! elem id 5

        !! elem coords
        allocate(elemCoords(2*numTotElems))
        elemCoords=(/1.5,1.5/)  ! 5

        local_pts=1
        test_coordx=1.5
        test_coordy=1.5

      endif
    endif



    ! Create Mesh structure in 1 step
    myMesh=ESMF_MeshCreate(parametricDim=2,spatialDim=2, &
         coordSys=ESMF_COORDSYS_CART, &
         nodeIds=nodeIds, nodeCoords=nodeCoords, &
         nodeOwners=nodeOwners, elementIds=elemIds,&
         elementTypes=elemTypes, elementConn=elemConn, &
         elementMask=elemMask, elementCoords=elemCoords, &
         rc=localrc)

    if (localrc /=ESMF_SUCCESS) then
      print*,'ERROR:  trouble creating mesh'
      rc=ESMF_FAILURE
      return
    endif

    ! After the creation we are through with the arrays, so they may be
    ! deallocated.
    deallocate(nodeIds)
    deallocate(nodeCoords)
    deallocate(nodeOwners)
    deallocate(elemMask)
    deallocate(elemCoords)
    deallocate(elemIds)
    deallocate(elemTypes)
    deallocate(elemConn)

    ! convert mask values
    maskValues=(/1,2/)

    maxpts=-99
    mypts=-99
    mydims=-99
    myid=-99

    pointlist=ESMF_PointListCreate(myMesh, ESMF_MESHLOC_ELEMENT, &
                                   maskValues=maskValues, rc=localrc)
    if (localrc /= ESMF_SUCCESS) then
      print*,'ERROR:  trouble creating pointlist'
      rc=ESMF_FAILURE
      return
    endif       

    call ESMF_PointListGet(pointlist, dims=mydims, numpts=mypts, maxpts=maxpts, rc=localrc)
    if (localrc /= ESMF_SUCCESS) then
       print*,'ERROR: trouble accessing pointlist data with get routine'
       rc=ESMF_FAILURE
       return
    endif

    if (maxpts .ne. local_pts .or. mypts .ne. local_pts .or. mydims .ne. 2) then
       print*,'ERROR:  unexpected values for newly created pointlist:'
       print*,'maxpts should be: ',local_pts,' got: ',maxpts
       print*,'numpts should be: ',local_pts,' got: ',mypts
       print*,'dims should be: 2  got: ',mydims
       rc=ESMF_FAILURE
       return
    endif


!    call ESMF_PointListPrint(pointlist)
!    if (localrc /= ESMF_SUCCESS) then
!       rc=ESMF_FAILURE
!       return
!    endif

    if (local_pts .gt. 0) then

      !locations values are zero based
      call ESMF_PointListGetForLoc(pointlist,0,loc_coords=test_coords,rc=localrc)
      if (localrc /= ESMF_SUCCESS) then
        print*,'ERROR:  trouble accessing pointlist data with get for location routine'
        rc=ESMF_FAILURE
        return
      endif

      my_err1 = abs(test_coordx - test_coords(1))
      my_err2 = abs(test_coordy - test_coords(2))
      if (my_err1 .gt. .0001 .or. my_err2 .gt. .0001) then
        print*,'ERROR:  unexpected coordinates for queried pointlist location:'
        print*,'expected ( ',test_coordx,' , ',test_coordy,' )  got  (',test_coords(1),',',test_coords(2),')'
        rc=ESMF_FAILURE
        return
      endif


    endif

    call ESMF_PointListDestroy(pointlist,rc=localrc)
    if (localrc /= ESMF_SUCCESS) then
      print*,'ERROR:  trouble destroying pointlist'
      rc=ESMF_FAILURE
      return
    endif       

    call ESMF_MeshDestroy(myMesh, rc=localrc)
    if (localrc /=ESMF_SUCCESS) then
      print*,'ERROR:  trouble destroying mesh'
      rc=ESMF_FAILURE
      return
    endif

    rc=ESMF_SUCCESS
  end subroutine test_pointlist_from_mesh_elems_wmask