Write from Fortran
From XdmfWeb
Jump to navigationJump to search
Writing Xdmf from Fortran
Xdmf can be generated in many different manners. Using the low level HDF5 library and print statements is certainly one of them. Utilizing the XDMF API, however, provides some convenient advantages. Consider the following Fortran program which build a grid of Hexahedra with some node and cell centered values :
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Create a Grid of Hexahedron Centered at 0,0,0
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE CreateGrid( IDIM, JDIM, KDIM, XYZ, ICONN )
INTEGER IDIM, JDIM, KDIM
REAL*8 XYZ
DIMENSION XYZ( 3, IDIM, JDIM, KDIM )
INTEGER ICONN
DIMENSION ICONN ( 8, ( IDIM - 1 ) * ( JDIM - 1 ) * ( KDIM - 1 ))
INTEGER I, J, K, IDX
REAL*8 X, Y, Z, DX, DY, DZ
C Print *, 'Size = ', IDIM, JDIM, KDIM
PRINT *, 'Initialze Problem'
C XYZ Values of Nodes
C From -1 to 1
DX = 2.0 / ( IDIM - 1 )
DY = 2.0 / ( JDIM - 1 )
DZ = 2.0 / ( KDIM - 1 )
Z = -1.0
DO 112 K= 1, KDIM
Y = -1.0
DO 111 J= 1, JDIM
X = -1.0
DO 110 I= 1, IDIM
XYZ( 1, I, J, K ) = X
XYZ( 2, I, J, K ) = Y
XYZ( 3, I, J, K ) = Z
X = X + DX
110 CONTINUE
Y = Y + DY
111 CONTINUE
Z = Z + DZ
112 CONTINUE
C Connections
IDX = 1
DO 122 K= 0, KDIM - 2
DO 121 J= 0, JDIM - 2
DO 120 I= 1, IDIM - 1
ICONN( 1, IDX ) = ( K * JDIM * IDIM ) + ( J * IDIM ) + I
ICONN( 2, IDX ) = ( K * JDIM * IDIM ) + ( J * IDIM ) + I + 1
ICONN( 3, IDX ) = ( ( K + 1 ) * JDIM * IDIM ) + ( J * IDIM ) + I + 1
ICONN( 4, IDX ) = ( ( K + 1 ) * JDIM * IDIM ) + ( J * IDIM ) + I
ICONN( 5, IDX ) = ( K * JDIM * IDIM ) + ( ( J + 1 ) * IDIM ) + I
ICONN( 6, IDX ) = ( K * JDIM * IDIM ) + ( ( J + 1 ) * IDIM ) + I + 1
ICONN( 7, IDX ) = ( ( K + 1 ) * JDIM * IDIM ) +
C ( ( J + 1 ) * IDIM ) + I + 1
ICONN( 8, IDX ) = ( ( K + 1 ) * JDIM * IDIM ) +
C ( ( J + 1 ) * IDIM ) + I
IDX = IDX + 1
120 CONTINUE
121 CONTINUE
122 CONTINUE
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Create a Node Centered Solution Field
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE NodeData( IDIM, JDIM, KDIM, XYZ, NCVALUES)
INTEGER IDIM, JDIM, KDIM
REAL*8 XYZ
DIMENSION XYZ( 3, IDIM, JDIM, KDIM )
REAL*8 NCVALUES
DIMENSION NCVALUES( IDIM, JDIM, KDIM )
INTEGER I, J, K
REAL*8 X, Y, Z
PRINT *, 'Calculating Node Centered Data'
DO 212, K=1, KDIM
DO 211, J=1, JDIM
DO 210, I=1, IDIM
X = XYZ( 1, I, J, K )
Y = XYZ( 2, I, J, K )
Z = XYZ( 3, I, J, K )
NCVALUES( I, J, K ) = SQRT( ( X * X ) + ( Y * Y ) + ( Z * Z ))
210 CONTINUE
211 CONTINUE
212 CONTINUE
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Create a Cell Centered Solution Field
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE CellData( IDIM, JDIM, KDIM, ITER, KICKER, XYZ, CCVALUES)
INTEGER IDIM, JDIM, KDIM, ITER, KICKER
REAL*8 XYZ
DIMENSION XYZ( 3, IDIM, JDIM, KDIM )
REAL*8 CCVALUES
DIMENSION CCVALUES( IDIM - 1, JDIM - 1, KDIM - 1 )
INTEGER I, J, K
PRINT *, 'Calculating Cell Centered Data for Iteration ', ITER
DO 312, K=1, KDIM - 1
DO 311, J=1, JDIM - 1
DO 310, I=1, IDIM - 1
X = XYZ( 1, I, J, K )
CCVALUES( I, J, K ) =
C SIN( ( ( X + 1 ) * IDIM * KICKER ) / 3 * ITER ) /
C EXP( X / ( 1.0 * ITER ) )
310 CONTINUE
311 CONTINUE
312 CONTINUE
C Waste Time
DO 313 I=1, 1000000
X = 0.1 * ITER / I
Y = SQRT( X * X )
Z = EXP( Y )
313 CONTINUE
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Main Program :
C Initialize Grid
C Initialize Node Centered Data
C For Iteration = 1 to 10
C Generate Cell Centered Data
C Done
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
PROGRAM HexMesh
PARAMETER ( IDIM = 11 )
PARAMETER ( JDIM = 13 )
PARAMETER ( KDIM = 15 )
REAL*8 XYZ
DIMENSION XYZ( 3, IDIM, JDIM, KDIM )
REAL*8 NCVALUES
DIMENSION NCVALUES( IDIM, JDIM, KDIM )
C
REAL*8 CCVALUES
DIMENSION CCVALUES( IDIM - 1, JDIM - 1, KDIM - 1 )
C
INTEGER ICONN
DIMENSION ICONN ( 8, ( IDIM - 1 ) * ( JDIM - 1 ) * ( KDIM - 1 ))
C
INTEGER ITER, KICKER, NITER, NARG
INTEGER IUNIT
CHARACTER*80 ARGIN
C
NARG = IARGC()
IF( NARG .GE. 1 ) THEN
CALL GETARG( 1, ARGIN )
READ( ARGIN, '(I)') NITER
ELSE
NITER = 10
ENDIF
CALL CreateGrid ( IDIM, JDIM, KDIM, XYZ, ICONN )
CALL NodeData( IDIM, JDIM, KDIM, XYZ, NCVALUES)
C
IUNIT = 14
OPEN( IUNIT, FILE='XYZ.dat', STATUS='unknown' )
REWIND IUNIT
WRITE ( IUNIT, * ) IDIM * JDIM * KDIM
WRITE ( IUNIT, * ) XYZ
CLOSE ( IUNIT )
C
IUNIT = 14
OPEN( IUNIT, FILE='CONN.dat', STATUS='unknown' )
REWIND IUNIT
WRITE ( IUNIT, * ) 'Hex', ( IDIM - 1 ) * ( JDIM - 1 ) * ( KDIM - 1 )
WRITE ( IUNIT, * ) ICONN
CLOSE ( IUNIT )
C
IUNIT = 14
OPEN( IUNIT, FILE='NodeValues.dat', STATUS='unknown' )
REWIND IUNIT
WRITE ( IUNIT, * ) NCVALUES
CLOSE ( IUNIT )
C
IUNIT = 14
OPEN( IUNIT, FILE='CellValues.dat', STATUS='unknown' )
REWIND IUNIT
C
KICKER = NITER
DO 1000 ITER = 1, NITER
CALL CellData( IDIM, JDIM, KDIM, ITER, KICKER, XYZ, CCVALUES)
WRITE ( IUNIT, * ) CCVALUES
1000 CONTINUE
CLOSE ( IUNIT )
C
END
Running the program produces four files XYZ.dat, CONN.dat, NodeValues.dat and CellValues.dat. Now let's suppose we want to generate Xdmf. We could use the HDF5 fortran bindings and produce the XML via print statements (we could also shove bamboo shoots under our fingernails). A better approach would be to create a 'C' wrapper routine that we call from fortran which contains all of the necessary Xdmf calls to produce valid Xdmf XML and HDF5 files :
#include <Xdmf.h>
// We want the filenames to be based on the iteration
// and padded with zeros
using std::setw;
using std::setfill;
// This works with g77. Different compilers require different
// name mangling.
#define XdmfWrite xdmfwrite_
//
// C/C++ expect NULL terminated strings. Here is a conversion subroutine.
char *
DemoConvertFortranString( char *FtnName ) {
static char Name[80];
char *np;
memcpy(Name, FtnName, 79 );
Name[79] = '\0';
np = &Name[78];
while( ( np > Name ) && ( *np <= ' ') ) {
np--;
}
*np = '\0';
return( Name );
}
//
// C++ will mangle the name based on the argument list. This tells the
// compiler not to mangle the name so we can call it from 'C' (but
// really Fortran in this case)
//
extern "C" {
//
void
XdmfWrite( char *FtnName, int *Iteration,
int *NumberOfPoints, int *NumberOfHex, XdmfFloat64 *Points,
XdmfInt32 *Conns, XdmfFloat64 *NodeData,
XdmfFloat64 *CellData){
char *Name;
char FullName[80];
ostrstream DataName(FullName, 80);
XdmfDOM dom;
XdmfRoot root;
XdmfDomain domain;
XdmfGrid grid;
XdmfTime time;
XdmfTopology *topology;
XdmfGeometry *geometry;
XdmfAttribute nodedata;
XdmfAttribute celldata;
XdmfArray *array;
Name = DemoConvertFortranString( FtnName );
root.SetDOM(&dom);
root.SetVersion(2.0);
root.Build();
// Domain
root.Insert(&domain);
// Grid
grid.SetName("Demonstration Grid");
domain.Insert(&grid);
time.SetTimeType(XDMF_TIME_SINGLE);
time.SetValue(0.001 * *Iteration);
grid.Insert(&time);
// Topology
topology = grid.GetTopology();
topology->SetTopologyType(XDMF_HEX);
topology->SetNumberOfElements(*NumberOfHex);
// Fortran is 1 based while c++ is 0 based so
// Either subtract 1 from all connections or specify a BaseOffset
topology->SetBaseOffset(1);
// If you haven't assigned an XdmfArray, GetConnectivity() will create one.
array = topology->GetConnectivity();
array->SetNumberOfElements(*NumberOfHex * 8);
array->SetValues(0, Conns, *NumberOfHex * 8);
// C++ string hocus pocus.
// We're actually building the string in FullName[] but were using streams.
// the DatasetName will be Demo_00001.h5:/Conns.
DataName.seekp(0);
DataName << Name << "_" << setw(5) << setfill('0') << *Iteration << ".h5:/Conns" << ends;
// Where the data will actually be written
array->SetHeavyDataSetName(FullName);
// Geometry
geometry = grid.GetGeometry();
geometry->SetGeometryType(XDMF_GEOMETRY_XYZ);
geometry->SetNumberOfPoints(*NumberOfPoints);
array = geometry->GetPoints();
array->SetNumberType(XDMF_FLOAT64_TYPE);
array->SetValues(0, Points, *NumberOfPoints * 3);
DataName.seekp(0);
DataName << Name << "_" << setw(5) << setfill('0') << *Iteration << ".h5:/Points" << ends;
array->SetHeavyDataSetName(FullName);
// Node Data
nodedata.SetName("Node Scalar");
nodedata.SetAttributeCenter(XDMF_ATTRIBUTE_CENTER_NODE);
nodedata.SetAttributeType(XDMF_ATTRIBUTE_TYPE_SCALAR);
array = nodedata.GetValues();
array->SetNumberType(XDMF_FLOAT64_TYPE);
array->SetNumberOfElements(*NumberOfPoints);
array->SetValues(0, NodeData, *NumberOfPoints);
DataName.seekp(0);
DataName << Name << "_" << setw(5) << setfill('0') << *Iteration << ".h5:/NodeData" << ends;
array->SetHeavyDataSetName(FullName);
// Cell Data
celldata.SetName("Cell Scalar");
celldata.SetAttributeCenter(XDMF_ATTRIBUTE_CENTER_CELL);
celldata.SetAttributeType(XDMF_ATTRIBUTE_TYPE_SCALAR);
array = celldata.GetValues();
array->SetNumberType(XDMF_FLOAT64_TYPE);
array->SetNumberOfElements(*NumberOfHex);
array->SetValues(0, CellData, *NumberOfHex);
DataName.seekp(0);
DataName << Name << "_" << setw(5) << setfill('0') << *Iteration << ".h5:/CellData" << ends;
array->SetHeavyDataSetName(FullName);
// Attach and Write
grid.Insert(&nodedata);
grid.Insert(&celldata);
// Build is recursive ... it will be called on all of the child nodes.
// This updates the DOM and writes the HDF5
root.Build();
// Write the XML
DataName.seekp(0);
DataName << Name << "_" << setw(5) << setfill('0') << *Iteration << ".xmf" << ends;
dom.Write(FullName);
}
}
Now we modify the Original Fortran code to call our new 'C' wrapper subroutine XdmfWrite()
Original :
KICKER = NITER
DO 1000 ITER = 1, NITER
CALL CellData( IDIM, JDIM, KDIM, ITER, KICKER, XYZ, CCVALUES)
WRITE ( IUNIT, * ) CCVALUES
1000 CONTINUE
CLOSE ( IUNIT )
C
END
Modified :
INHEX = ( IDIM - 1 ) * ( JDIM - 1 ) * ( KDIM - 1 )
INPNT = IDIM * JDIM * KDIM
KICKER = NITER
DO 1000 ITER = 1, NITER
CALL CellData( IDIM, JDIM, KDIM, ITER, KICKER, XYZ, CCVALUES)
WRITE ( IUNIT, * ) CCVALUES
CALL XDMFWRITE( 'Demo', ITER, INPNT, INHEX, XYZ, ICONN, NCVALUES, CCVALUES)
1000 CONTINUE
CLOSE ( IUNIT )
C
END
The program now produces one HDF5 file and one XML file for each iteration. Here is one of the XML output files :
<?xml version="1.0" ?> <!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []> <Xdmf xmlns:xi="http://www.w3.org/2003/XInclude" Version="2"> <Domain> <Grid Name="Demonstration Grid" GridType="Uniform"> <Topology TopologyType="Hexahedron" NumberOfElements="1680 " BaseOffset="1"> <DataItem Dimensions="13440 " NumberType="Float" Precision="4" Format="HDF">Demo_00004.h5:/Conns</DataItem> </Topology> <Geometry GeometryType="XYZ"> <DataItem Dimensions="6435 " NumberType="Float" Precision="4" Format="HDF">Demo_00004.h5:/Points</DataItem> </Geometry>
The final step would be to group all of these together so we could animate them over time. Instead of duplicating the XML in another file for each grid, we can use a little XML magic to pull them together :
<?xml version="1.0" ?> <!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []> <Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.0"> <Domain> <Grid GridType="Collection" CollectionType="Temporal"> <xi:include href="Demo_00001.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00002.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00003.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00004.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00005.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00006.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00007.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00008.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00009.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> <xi:include href="Demo_00010.xmf" xpointer="xpointer(//Xdmf/Domain/Grid)" /> </Grid> </Domain> </Xdmf>
This uses XInclude and XPointer terminology to create a Temporal Collection Grid that consists of the first (in our case, the only) grid in each file. This XML file could be generated programmtically or by a post processing script.
