Dear Petsc Team, I am playing around with DMPlex, using it to generate the Mesh for the ICON weather model(http://doi.org/10.1002/2015MS000431), which employs a triangle mesh horizontally and columns, vertically.
This results in a grid, looking like prisms, where top and bottom faces
are triangles and side faces are rectangles.
I was delighted to see that I could export the triangle DMPlex (2d Mesh)
to hdf5 and use petsc_gen_xdmf.py to then visualize the mesh in
visit/paraview.
This is especially nice when exporting petscsections/vectors directly to
VTK.
I then tried the same approach for the prism grid in 3D.
I attached the code for one single cell, as well as the output in hdf5.
However, trying to convert the hdf5 output, it fails with:
make prism.xmf
$PETSC_DIR/bin/petsc_gen_xdmf.py prism.h5
Traceback (most recent call last):
File
"/software/meteo/xenial/x86_64/petsc/master/debug_gcc/..//bin/petsc_gen_xdmf.py",
line 241, in <module>
generateXdmf(f)
File
"/software/meteo/xenial/x86_64/petsc/master/debug_gcc/..//bin/petsc_gen_xdmf.py",
line 235, in generateXdmf
Xdmf(xdmfFilename).write(hdfFilename, topoPath, numCells,
numCorners, cellDim, geomPath, numVertices, spaceDim, time, vfields,
cfields)
File
"/software/meteo/xenial/x86_64/petsc/master/debug_gcc/..//bin/petsc_gen_xdmf.py",
line 193, in write
self.writeSpaceGridHeader(fp, numCells, numCorners, cellDim, spaceDim)
File
"/software/meteo/xenial/x86_64/petsc/master/debug_gcc/..//bin/petsc_gen_xdmf.py",
line 75, in writeSpaceGridHeader
''' % (self.cellMap[cellDim][numCorners], numCells, "XYZ" if
spaceDim > 2 else "XY"))
KeyError: 6
Also, if I try to export a vector directly to vtk, visit and paraview
fail to open it.
My question is:
Is this a general limitation of these output formats, that I can not mix
faces with 3 and 4 vertices or is it a limitation of the
petsc_gen_xdmf.py or the VTK Viewer.
I'd also welcome any thoughts on the prism mesh in general.
Is it that uncommon to use and do you foresee other complications with it?
I fear I cannot change the discretization of the host model but maybe it
makes sense to use a different grid for my radiative transfer code?
Many thanks,
Fabian
include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules
prism.xmf:: prism.h5
${PETSC_DIR}/bin/petsc_gen_xdmf.py prism.h5
prism.h5:: plex_prism
./plex_prism -show_plex ::ascii_info_detail
./plex_prism -show_plex hdf5:prism.h5
plex_prism:: plex_prism.F90
${PETSC_FCOMPILE} -c plex_prism.F90
${FLINKER} plex_prism.o -o plex_prism ${PETSC_LIB}
clean::
rm -rf *.o prism.h5 prism.xmf plex_prism
program main
#include "petsc/finclude/petsc.h"
use petsc
implicit none
PetscErrorCode :: ierr
PetscInt :: i
DM :: dm
call PetscInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)
call create_plex()
call PetscObjectViewFromOptions(dm, PETSC_NULL_VEC, "-show_plex", ierr); CHKERRQ(ierr)
call DMDestroy(dm, ierr);CHKERRQ(ierr)
call PetscFinalize(ierr)
contains
subroutine create_plex()
call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRQ(ierr)
call PetscObjectSetName(dm, 'testplex', ierr);CHKERRQ(ierr)
call DMSetDimension(dm, 3, ierr);CHKERRQ(ierr)
call DMPlexSetChart(dm, 0, 21, ierr); CHKERRQ(ierr)
! Preallocation
! cell has 5 faces
call DMPlexSetConeSize(dm, 0, 5, ierr); CHKERRQ(ierr)
! Faces have 3 or 4 edges
do i=1,2
call DMPlexSetConeSize(dm, i, 3, ierr); CHKERRQ(ierr)
enddo
do i=3,5
call DMPlexSetConeSize(dm, i, 4, ierr); CHKERRQ(ierr)
enddo
! Edges have 2 vertices
do i=6,14
call DMPlexSetConeSize(dm, i, 2, ierr); CHKERRQ(ierr)
enddo
call DMSetUp(dm, ierr); CHKERRQ(ierr) ! Allocate space for cones
! Setup Connections
! 0 cell has faces 1(bottom) and 2(top) as well as three faces on sides
call DMPlexSetCone(dm, 0, [1,2,3,4,5], ierr); CHKERRQ(ierr)
! bottom face has 3 edges
call DMPlexSetCone(dm, 1, [ 6, 7, 8], ierr); CHKERRQ(ierr)
! top face has 3 edges
call DMPlexSetCone(dm, 2, [12,13,14], ierr); CHKERRQ(ierr)
! side faces have 4 edges
call DMPlexSetCone(dm, 3, [ 7, 9,11,13], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 4, [ 6, 9,10,12], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 5, [ 8,10,11,14], ierr); CHKERRQ(ierr)
! and the corresponding vertices
call DMPlexSetCone(dm, 6, [15,16], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 7, [15,17], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 8, [16,17], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 9, [15,18], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 10, [16,19], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 11, [17,20], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 12, [18,19], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 13, [18,20], ierr); CHKERRQ(ierr)
call DMPlexSetCone(dm, 14, [19,20], ierr); CHKERRQ(ierr)
call DMPlexSymmetrize(dm, ierr); CHKERRQ(ierr)
call DMPlexStratify(dm, ierr); CHKERRQ(ierr)
call set_coords()
end subroutine
subroutine set_coords()
PetscScalar, pointer :: coords(:)
Vec :: coordinates
PetscInt :: dimEmbed, coordSize, vStart, vEnd, pStart, pEnd, voff
PetscSection :: coordSection
call DMGetCoordinateDim(dm, dimEmbed, ierr); CHKERRQ(ierr)
call DMGetCoordinateSection(dm, coordSection, ierr); CHKERRQ(ierr)
call PetscSectionSetNumFields(coordSection, 1, ierr); CHKERRQ(ierr)
call PetscSectionSetUp(coordSection, ierr); CHKERRQ(ierr)
call PetscSectionSetFieldComponents(coordSection, 0, dimEmbed, ierr); CHKERRQ(ierr)
call DMPlexGetChart(dm, pStart, pEnd, ierr); CHKERRQ(ierr)
call DMPlexGetDepthStratum (dm, 0, vStart, vEnd, ierr); CHKERRQ(ierr) ! vertices
call PetscSectionSetChart(coordSection, pStart, pEnd, ierr);CHKERRQ(ierr)
do i = vStart, vEnd-1
call PetscSectionSetDof(coordSection, i, dimEmbed, ierr); CHKERRQ(ierr)
call PetscSectionSetFieldDof(coordSection, i, 0, dimEmbed, ierr); CHKERRQ(ierr)
enddo
call PetscSectionSetUp(coordSection, ierr); CHKERRQ(ierr)
call PetscSectionGetStorageSize(coordSection, coordSize, ierr); CHKERRQ(ierr)
call VecCreate(PETSC_COMM_SELF, coordinates, ierr); CHKERRQ(ierr)
call VecSetSizes(coordinates, coordSize, PETSC_DETERMINE, ierr);CHKERRQ(ierr)
call VecSetBlockSize(coordinates, dimEmbed, ierr);CHKERRQ(ierr)
call VecSetType(coordinates, VECSTANDARD, ierr);CHKERRQ(ierr)
call PetscObjectSetName(coordinates, "coordinates", ierr); CHKERRQ(ierr)
call VecGetArrayF90(coordinates, coords, ierr); CHKERRQ(ierr)
i = 15; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,0,0]
i = 16; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,2,0]
i = 17; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [2,1,0]
i = 18; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,0,1]
i = 19; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,2,1]
i = 20; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [2,1,1]
call VecRestoreArrayF90(coordinates, coords, ierr); CHKERRQ(ierr)
call DMSetCoordinatesLocal(dm, coordinates, ierr);CHKERRQ(ierr)
call PetscObjectViewFromOptions(coordinates, PETSC_NULL_VEC, "-show_plex_coordinates", ierr); CHKERRQ(ierr)
call VecDestroy(coordinates, ierr);CHKERRQ(ierr)
end subroutine
end program
prism.h5
Description: HDF file
signature.asc
Description: OpenPGP digital signature
