Write from Fortran

Writing Xdmf from Fortran

Note code snippets on this page are from version 2 or 1 and need updating for xdmf 3.

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>
     
     <Attribute Name="Node Scalar" AttributeType="Scalar" Center="Node">
       <DataItem Dimensions="2145 " NumberType="Float" Precision="4" Format="HDF">Demo_00004.h5:/NodeData</DataItem>
     </Attribute>
     <Attribute Name="Cell Scalar" AttributeType="Scalar" Center="Cell">
       <DataItem Dimensions="1680 " NumberType="Float" Precision="4" Format="HDF">Demo_00004.h5:/CellData</DataItem>
     </Attribute>
   </Grid>
 </Domain>
</Xdmf>

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.