subroutine check_mesh_from_sph_pb_EM_file(correct, rc)
type(ESMF_Mesh) :: mesh
logical :: correct
integer :: rc
integer, pointer :: nodeIds(:),nodeOwners(:),nodeMask(:)
real(ESMF_KIND_R8), pointer :: nodeCoords(:)
real(ESMF_KIND_R8), pointer :: ownedNodeCoords(:)
integer :: numNodes, numOwnedNodes, numOwnedNodesTst
integer :: numElems,numOwnedElemsTst
integer :: numElemConns
integer :: numTriElems, numQuadElems, numPolyBreakElems
real(ESMF_KIND_R8), pointer :: elemCoords(:)
integer, pointer :: elemIds(:),elemTypes(:),elemConn(:)
integer, allocatable :: elemMask(:)
real(ESMF_KIND_R8),allocatable :: elemArea(:)
integer :: petCount, localPet
type(ESMF_VM) :: vm
type(ESMF_DistGrid) :: elemDistgrid
type(ESMF_Field) :: nodeField, elemField
integer :: i,j,k
integer :: numNodesTst, numElemsTst,numElemConnsTst
integer,allocatable :: elemIdsTst(:)
integer,allocatable :: elemTypesTst(:)
integer,allocatable :: elemConnTst(:)
real(ESMF_KIND_R8),allocatable :: elemAreaTst(:)
integer,allocatable :: nodeIdsTst(:)
real(ESMF_KIND_R8),allocatable :: nodeCoordsTst(:)
integer,allocatable :: nodeOwnersTst(:)
integer,allocatable :: nodeMaskTst(:)
integer,allocatable :: elemMaskTst(:)
real(ESMF_KIND_R8), allocatable :: elemCoordsTst(:)
logical :: nodeMaskIsPresentTst
logical :: elemMaskIsPresentTst
logical :: elemAreaIsPresentTst
logical :: elemCoordsIsPresentTst
integer :: pdim, sdim
type(ESMF_CoordSys_Flag) :: coordSys
! get global VM
call ESMF_VMGetGlobal(vm, rc=rc)
if (rc /= ESMF_SUCCESS) return
call ESMF_VMGet(vm, localPet=localPet, petCount=petCount, rc=rc)
if (rc /= ESMF_SUCCESS) return
! return with an error if not 1 or 4 PETs
if ((petCount /= 1) .and. (petCount /=4)) then
rc=ESMF_FAILURE
return
endif
! Setup mesh info depending on the
! number of PETs
if (petCount .eq. 1) then
! Fill in node data
numNodes=16
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/1,2,3,4,5,6,7,8, &
9,10,11,12,13,14,&
15,16/)
!! node Mask
allocate(nodeMask(numNodes))
nodeMask=(/1,0,1,0,1,0,1,0, &
1,0,1,0,1,0,&
1,0/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/0.0,0.0, & ! 1
1.0,0.0, & ! 2
2.0,0.0, & ! 3
3.0,0.0, & ! 4
0.0,1.0, & ! 5
1.0,1.0, & ! 6
2.0,1.0, & ! 7
3.0,1.0, & ! 8
0.0,2.0, & ! 9
1.0,2.0, & ! 10
2.0,2.0, & ! 11
3.0,2.0, & ! 12
0.0,3.0, & ! 13
1.0,3.0, & ! 14
2.0,3.0, & ! 15
3.0,3.0 /) ! 16
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=0 ! everything on proc 0
! Fill in elem data
numTriElems=4
numQuadElems=6
numPolyBreakElems=1
numElems=numTriElems+numQuadElems+numPolyBreakElems
numElemConns=3*numTriElems+ &
4*numQuadElems+ &
7*numPolyBreakElems
!! elem ids
allocate(elemIds(numElems))
elemIds=(/1,2,3,4,5,6,7,8,9,10,11/)
!! elem types
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! 1
ESMF_MESHELEMTYPE_TRI, & ! 2
ESMF_MESHELEMTYPE_TRI, & ! 3
ESMF_MESHELEMTYPE_TRI, & ! 4
ESMF_MESHELEMTYPE_QUAD, & ! 5
ESMF_MESHELEMTYPE_QUAD, & ! 6
ESMF_MESHELEMTYPE_QUAD, & ! 7
ESMF_MESHELEMTYPE_QUAD, & ! 8
ESMF_MESHELEMTYPE_QUAD, & ! 9
ESMF_MESHELEMTYPE_QUAD, & ! 10
7/) ! 11
!! elem coords
allocate(elemCoords(2*numElems))
elemCoords=(/0.25,0.75, & ! 1
0.75,0.25, & ! 2
1.25,0.75, & ! 3
1.75,0.25, & ! 4
2.5,0.5, & ! 5
0.5,1.5, & ! 6
1.5,1.5, & ! 7
2.5,1.5, & ! 8
0.5,2.5, & ! 9
1.5,2.5, & ! 10
2.5,2.5/) ! 11
!! elem conn
allocate(elemConn(numElemConns))
elemConn=(/1,6,5, & ! 1
1,2,6, & ! 2
2,7,6, & ! 3
2,3,7, & ! 4
3,4,8,7, & ! 5
5,6,10,9, & ! 6
6,7,11,10, & ! 7
7,8,12,11, & ! 8
9,10,14,13, & ! 9
10,11,15,14, & ! 10
11,12,16, ESMF_MESH_POLYBREAK, 11,16,15/) ! 11
!! elem mask
allocate(elemMask(numElems))
elemMask=(/1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0/)
!! elem area
!! (need kinds because not R4 doesn't rep. close enough)
allocate(elemArea(numElems))
elemArea=(/0.51_ESMF_KIND_R8, 0.52_ESMF_KIND_R8, 0.53_ESMF_KIND_R8, 0.54_ESMF_KIND_R8, &
1.5_ESMF_KIND_R8, 1.6_ESMF_KIND_R8, 1.7_ESMF_KIND_R8, 1.8_ESMF_KIND_R8, &
1.9_ESMF_KIND_R8, 1.1_ESMF_KIND_R8, 1.11_ESMF_KIND_R8/)
else if (petCount .eq. 4) then
! Setup mesh data depending on PET
if (localPet .eq. 0) then
! Fill in node data
numNodes=4
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/1,2,5,6/)
!! node Mask
allocate(nodeMask(numNodes))
nodeMask=(/1,0,1,0/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/0.0,0.0, & ! 1
1.0,0.0, & ! 2
0.0,1.0, & ! 5
1.0,1.0 /) ! 6
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=0 ! everything on proc 0
! Fill in elem data
numTriElems=2
numQuadElems=0
numElems=numTriElems+numQuadElems
numElemConns=3*numTriElems+4*numQuadElems
!! elem ids
allocate(elemIds(numElems))
elemIds=(/1,2/)
!! elem types
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! 1
ESMF_MESHELEMTYPE_TRI/) ! 2
!! elem coords
allocate(elemCoords(2*numElems))
elemCoords=(/0.25,0.75, & ! 1
0.75,0.25/) ! 2
!! elem conn
allocate(elemConn(numElemConns))
elemConn=(/1,4,3, & ! 1
1,2,4 /) ! 2
!! elem mask
allocate(elemMask(numElems))
elemMask=(/1,0/)
!! elem area
!! (need kinds because not R4 doesn't rep. close enough)
allocate(elemArea(numElems))
elemArea=(/0.51_ESMF_KIND_R8, 0.52_ESMF_KIND_R8/)
else if (localPet .eq. 1) then
! Fill in node data
numNodes=6
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/2,3,4,6,7,8/)
!! node Mask
allocate(nodeMask(numNodes))
nodeMask=(/0,1,0,0,1,0/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/1.0,0.0, & ! 2
2.0,0.0, & ! 3
3.0,0.0, & ! 4
1.0,1.0, & ! 6
2.0,1.0, & ! 7
3.0,1.0 /) ! 8
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! 2
1, & ! 3
1, & ! 4
0, & ! 6
1, & ! 7
1/) ! 8
! Fill in elem data
numTriElems=2
numQuadElems=1
numElems=numTriElems+numQuadElems
numElemConns=3*numTriElems+4*numQuadElems
!! elem ids
allocate(elemIds(numElems))
elemIds=(/3,4,5/)
!! elem types
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! 3
ESMF_MESHELEMTYPE_TRI, & ! 4
ESMF_MESHELEMTYPE_QUAD/) ! 5
!! elem coords
allocate(elemCoords(2*numElems))
elemCoords=(/1.25,0.75, & ! 3
1.75,0.25, & ! 4
2.50,0.50/) ! 5
!! elem conn
allocate(elemConn(numElemConns))
elemConn=(/1,5,4, & ! 3
1,2,5, & ! 4
2,3,6,5/) ! 5
!! elem mask
allocate(elemMask(numElems))
elemMask=(/1,1,1/)
!! elem area
!! (need kinds because not R4 doesn't rep. close enough)
allocate(elemArea(numElems))
elemArea=(/0.53_ESMF_KIND_R8,0.54_ESMF_KIND_R8,1.5_ESMF_KIND_R8/)
else if (localPet .eq. 2) then
! Fill in node data
numNodes=9
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/5,6,7, &
9,10,11, &
13,14,15/)
!! node Mask
allocate(nodeMask(numNodes))
nodeMask=(/1,0,1, &
1,0,1, &
1,0,1/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/0.0,1.0, & ! 5
1.0,1.0, & ! 6
2.0,1.0, & ! 7
0.0,2.0, & ! 9
1.0,2.0, & ! 10
2.0,2.0, & ! 11
0.0,3.0, & ! 13
1.0,3.0, & ! 14
2.0,3.0/) ! 15
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! 5
0, & ! 6
1, & ! 7
2, & ! 9
2, & ! 10
2, & ! 11
2, & ! 13
2, & ! 14
2/) ! 15
! Fill in elem data
numTriElems=0
numQuadElems=4
numElems=numTriElems+numQuadElems
numElemConns=3*numTriElems+4*numQuadElems
!! elem ids
allocate(elemIds(numElems))
elemIds=(/6,7,9,10/)
!! elem types
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! 6
ESMF_MESHELEMTYPE_QUAD, & ! 7
ESMF_MESHELEMTYPE_QUAD, & ! 9
ESMF_MESHELEMTYPE_QUAD/) ! 10
!! elem coords
allocate(elemCoords(2*numElems))
elemCoords=(/0.5,1.5, & ! 6
1.5,1.5, & ! 7
0.5,2.5, & ! 9
1.5,2.5/) ! 10
!! elem conn
allocate(elemConn(numElemConns))
elemConn=(/1,2,5,4, & ! 6
2,3,6,5, & ! 7
4,5,8,7, & ! 9
5,6,9,8/) ! 10
!! elem mask
allocate(elemMask(numElems))
elemMask=(/1,1,1,1/)
!! elem area
!! (need kinds because not R4 doesn't rep. close enough)
allocate(elemArea(numElems))
elemArea=(/1.6_ESMF_KIND_R8,1.7_ESMF_KIND_R8,1.9_ESMF_KIND_R8,1.1_ESMF_KIND_R8/)
else if (localPet .eq. 3) then
! Fill in node data
numNodes=6
!! node ids
allocate(nodeIds(numNodes))
nodeIds=(/7,8,11,12,15,16/)
!! node Mask
allocate(nodeMask(numNodes))
nodeMask=(/1,0,1,0,1,0/)
!! node Coords
allocate(nodeCoords(numNodes*2))
nodeCoords=(/2.0,1.0, & ! 7
3.0,1.0, & ! 8
2.0,2.0, & ! 11
3.0,2.0, & ! 12
2.0,3.0, & ! 15
3.0,3.0 /) ! 16
!! node owners
allocate(nodeOwners(numNodes))
nodeOwners=(/1, & ! 7
1, & ! 8
2, & ! 11
3, & ! 12
2, & ! 15
3/) ! 16
! Fill in elem data
numTriElems=0
numQuadElems=1
numPolyBreakElems=1
numElems=numTriElems+numQuadElems+numPolyBreakElems
numElemConns=3*numTriElems+ &
4*numQuadElems+ &
7*numPolyBreakElems
!! elem ids
allocate(elemIds(numElems))
elemIds=(/8,11/)
!! elem types
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! 8
7/) ! 11
!! elem coords
allocate(elemCoords(2*numElems))
elemCoords=(/2.5,1.5, & ! 8
2.5,2.5/) ! 11
!! elem conn
allocate(elemConn(numElemConns))
elemConn=(/1,2,4,3, & ! 8
11,12,16, ESMF_MESH_POLYBREAK, 11,16,15/) ! 11
!! elem mask
allocate(elemMask(numElems))
elemMask=(/1,0/)
!! elem area
!! (need kinds because not R4 doesn't rep. close enough)
allocate(elemArea(numElems))
elemArea=(/1.8_ESMF_KIND_R8,1.11_ESMF_KIND_R8/)
endif
endif
! Create distgrid to ensure that elements are created on the
! same PET as described above
elemdistgrid=ESMF_DistGridCreate(elemIds, rc=rc)
if (rc /= ESMF_SUCCESS) return
! Read mesh from file that's the same as the one described by info set up above
mesh=ESMF_MeshCreate("data/test_sph_polybreak_esmf.nc", &
fileformat=ESMF_FILEFORMAT_ESMFMESH, &
addUserArea=.true., &
elementDistgrid=elemDistgrid, &
rc=rc)
if (rc /= ESMF_SUCCESS) return
! DEBUG OUTPUT
!call ESMF_MeshWrite(mesh,"mesh_pb_esmf",rc=rc)
!if (rc .ne. ESMF_SUCCESS) return
! Init correct to true before looking for problems
correct=.true.
! Get dim and coord info
call ESMF_MeshGet(mesh, &
parametricDim=pdim, &
spatialDim=sdim, &
coordSys=coordSys, &
rc=rc)
if (rc /= ESMF_SUCCESS) return
! Check Counts
if (pdim .ne. 2) correct=.false.
if (sdim .ne. 2) correct=.false.
if (coordSys .ne. ESMF_COORDSYS_SPH_DEG) correct=.false.
! Get counts
call ESMF_MeshGet(mesh, &
nodeCount=numNodesTst, &
elementCount=numElemsTst, &
! elementConnCount=numElemConnsTst, & ! DOESN'T WORK YET WITH NGONS
rc=rc)
if (rc /= ESMF_SUCCESS) return
! Check Counts
if (numNodes .ne. numNodesTst) correct=.false.
if (numElems .ne. numElemsTst) correct=.false.
! if (numElemConns .ne. numElemConnsTst) correct=.false. ! DOESN'T WORK YET WITH NGONS
! Debugging
!write(*,*) "numNodes=",numNodes,numNodesTst
!write(*,*) "numElems=",numElems,numElemsTst
!write(*,*) "numElemConns=",numElemConns,numElemConnsTst
! Get is present information
call ESMF_MeshGet(mesh, &
nodeMaskIsPresent=nodeMaskIsPresentTst, &
elementMaskIsPresent=elemMaskIsPresentTst, &
elementAreaIsPresent=elemAreaIsPresentTst, &
elementCoordsIsPresent=elemCoordsIsPresentTst, &
rc=rc)
if (rc /= ESMF_SUCCESS) return
! Debug
! write(*,*) "nodeMaskIsPresent=",nodeMaskIsPresentTst
! write(*,*) "elemMaskIsPresent=",elemMaskIsPresentTst
! write(*,*) "elemAreaIsPresent=",elemAreaIsPresentTst
! write(*,*) "elemCoordsIsPresent=",elemCoordsIsPresentTst
! Check is present info
if (.not. nodeMaskIsPresentTst) correct=.false.
if (.not. elemMaskIsPresentTst) correct=.false.
if (.not. elemAreaIsPresentTst) correct=.false.
if (.not. elemCoordsIsPresentTst) correct=.false.
! Allocate space for tst arrays
allocate(nodeIdsTst(numNodesTst))
allocate(nodeCoordsTst(2*numNodesTst))
allocate(nodeOwnersTst(numNodesTst))
allocate(nodeMaskTst(numNodesTst))
! allocate(elemIdsTst(numElemsTst)) ! DOESN'T WORK YET WITH NGONS
! allocate(elemTypesTst(numElemsTst)) ! DOESN'T WORK YET WITH NGONS
! allocate(elemConnTst(numElemConnsTst)) ! DOESN'T WORK YET WITH NGONS
! allocate(elemMaskTst(numElemsTst)) ! DOESN'T WORK YET WITH NGONS
! allocate(elemAreaTst(numElemsTst)) ! DOESN'T WORK YET WITH NGONS
! allocate(elemCoordsTst(2*numElemsTst)) ! DOESN'T WORK YET WITH NGONS
! XMRKX
! Get Information
call ESMF_MeshGet(mesh, &
nodeIds=nodeIdsTst, &
nodeCoords=nodeCoordsTst, &
nodeOwners=nodeOwnersTst, &
nodeMask=nodeMaskTst, &
! elementIds=elemIdsTst, & ! DOESN'T WORK YET WITH NGONS
! elementTypes=elemTypesTst, & ! DOESN'T WORK YET WITH NGONS
! elementConn=elemConnTst, & ! DOESN'T WORK YET WITH NGONS
! elementMask=elemMaskTst, & ! DOESN'T WORK YET WITH NGONS
! elementArea=elemAreaTst, & ! DOESN'T WORK YET WITH NGONS
! elementCoords=elemCoordsTst, & ! DOESN'T WORK YET WITH NGONS
rc=rc)
if (rc /= ESMF_SUCCESS) return
! do i=1,numElemsTst
! write(*,*) localPet,"# ",elemIdsTst(i),"ea=",elemAreaTst(i),"ec=",elemCoordsTst(2*i-1),elemCoordsTst(2*i),"em=",elemMaskTst(i)
! enddo
! Check node ids
do i=1,numNodesTst
if (nodeIds(i) .ne. nodeIdsTst(i)) correct=.false.
enddo
! Check node Coords
k=1
do i=1,numNodesTst ! Loop over nodes
do j=1,2 ! Loop over coord spatial dim
if (nodeCoords(k) .ne. nodeCoordsTst(k)) correct=.false.
k=k+1
enddo
enddo
! Check node Owners
do i=1,numNodesTst
if (nodeOwners(i) .ne. nodeOwnersTst(i)) correct=.false.
enddo
! Check node Mask
do i=1,numNodesTst
if (nodeMask(i) .ne. nodeMaskTst(i)) correct=.false.
enddo
! Debug output
! write(*,*) "nodeMask =",nodeMask
! write(*,*) "nodeMaskTst=",nodeMaskTst
! Check elem ids
! DOESN'T WORK YET WITH NGONS
!do i=1,numElemsTst
! if (elemIds(i) .ne. elemIdsTst(i)) correct=.false.
!enddo
! Check elem Types
! DOESN'T WORK YET WITH NGONS
!do i=1,numElemsTst
! if (elemTypes(i) .ne. elemTypesTst(i)) correct=.false.
!enddo
! Check elem Connections
! DOESN'T WORK YET WITH NGONS
!do i=1,numElemConnsTst
! if (elemConn(i) .ne. elemConnTst(i)) correct=.false.
!enddo
! Check elem mask
! DOESN'T WORK YET WITH NGONS
!do i=1,numElems
! if (elemMask(i) .ne. elemMaskTst(i)) correct=.false.
!enddo
! Check elem area
! DOESN'T WORK YET WITH NGONS
!do i=1,numElems
! if (elemArea(i) .ne. elemAreaTst(i)) correct=.false.
!enddo
! Check elem Coords
! DOESN'T WORK YET WITH NGONS
!k=1
!do i=1,numElemsTst ! Loop over nodes
! do j=1,2 ! Loop over coord spatial dim
! if (elemCoords(k) .ne. elemCoordsTst(k)) correct=.false.
! k=k+1
! enddo
!enddo
! Debug output
! write(*,*) "elemCoords =",elemCoords
! write(*,*) "elemCoordsTst=",elemCoordsTst
! deallocate node data
deallocate(nodeIds)
deallocate(nodeCoords)
deallocate(nodeOwners)
deallocate(nodeMask)
! deallocate elem data
deallocate(elemIds)
deallocate(elemTypes)
deallocate(elemCoords)
deallocate(elemConn)
deallocate(elemMask)
deallocate(elemArea)
! Deallocate tst Arrays
deallocate(nodeIdsTst)
deallocate(nodeCoordsTst)
deallocate(nodeOwnersTst)
deallocate(nodeMaskTst)
!deallocate(elemIdsTst) ! DOESN'T WORK YET WITH NGONS
!deallocate(elemTypesTst) ! DOESN'T WORK YET WITH NGONS
!deallocate(elemConnTst) ! DOESN'T WORK YET WITH NGONS
!deallocate(elemMaskTst) ! DOESN'T WORK YET WITH NGONS
!deallocate(elemAreaTst) ! DOESN'T WORK YET WITH NGONS
!deallocate(elemCoordsTst) ! DOESN'T WORK YET WITH NGONS
! Get rid of Mesh
call ESMF_MeshDestroy(mesh, rc=rc)
if (rc /= ESMF_SUCCESS) return
! Get rid of elemDistgrid
call ESMF_DistGridDestroy(elemDistgrid, rc=rc)
if (rc /= ESMF_SUCCESS) return
! Return success
rc=ESMF_SUCCESS
end subroutine check_mesh_from_sph_pb_EM_file