[cig-commits] commit: Make Q2 XDMF output work
Mercurial
hg at geodynamics.org
Sun Oct 2 12:22:52 PDT 2011
changeset: 873:f4fe195bb604
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sun Oct 02 12:21:15 2011 -0700
files: Utils/src/XDMFGenerator.cxx
description:
Make Q2 XDMF output work
diff -r f5bd1cb7d3f6 -r f4fe195bb604 Utils/src/XDMFGenerator.cxx
--- a/Utils/src/XDMFGenerator.cxx Thu Sep 29 15:09:34 2011 -0700
+++ b/Utils/src/XDMFGenerator.cxx Sun Oct 02 12:21:15 2011 -0700
@@ -111,185 +111,247 @@ void XDMFGenerator_GenerateAll( void* _c
}
void _XDMFGenerator_WriteFieldSchema( UnderworldContext* context, Stream* stream ) {
- FieldVariable* fieldVar = NULL;
- FeVariable* feVar = NULL;
- Mesh* mesh = NULL;
- FeMesh* feMesh = NULL;
- unsigned nDims;
- unsigned totalVerts;
- Index maxNodes;
- Index elementGlobalSize;
- Element_GlobalIndex gElement_I;
- Index var_I = 0;
- Bool saveCoords = Dictionary_GetBool_WithDefault( context->dictionary, (Dictionary_Entry_Key)"saveCoordsWithFields", False );
- Stream* errorStream = Journal_Register( Error_Type, (Name)CURR_MODULE_NAME );
- char* variableType = NULL;
- char* topologyType = NULL;
- unsigned componentCount = LiveComponentRegister_GetCount(stgLiveComponentRegister);
- unsigned compI;
- Stg_Component* stgComp;
+ FieldVariable* fieldVar = NULL;
+ FeVariable* feVar = NULL;
+ Mesh* mesh = NULL;
+ FeMesh* feMesh = NULL;
+ unsigned nDims;
+ unsigned totalVerts;
+ Index maxNodes;
+ Index elementGlobalSize;
+ Element_GlobalIndex gElement_I;
+ Index var_I = 0;
+ Bool saveCoords = Dictionary_GetBool_WithDefault( context->dictionary, (Dictionary_Entry_Key)"saveCoordsWithFields", False );
+ Stream* errorStream = Journal_Register( Error_Type, (Name)CURR_MODULE_NAME );
+ char* variableType = NULL;
+ char* topologyType = NULL;
+ unsigned componentCount = LiveComponentRegister_GetCount(stgLiveComponentRegister);
+ unsigned compI;
+ Stg_Component* stgComp;
- /** search for entire live component register for feMesh types **/
- for( compI = 0 ; compI < componentCount ; compI++ ){
- stgComp = LiveComponentRegister_At( stgLiveComponentRegister, compI );
- /* check that component is of type FeMesh, and that its element family is linear */
- if ( Stg_Class_IsInstance( stgComp, FeMesh_Type ) && !strcmp( ((FeMesh*)stgComp)->feElFamily, "linear" ) ) {
- mesh = ( Mesh*)stgComp;
- feMesh = (FeMesh*)stgComp;
+ /** search for entire live component register for feMesh types **/
+ for( compI = 0 ; compI < componentCount ; compI++ ){
+ stgComp = LiveComponentRegister_At( stgLiveComponentRegister, compI );
+ /* check that component is of type FeMesh, and that its element family is linear */
+ if(Stg_Class_IsInstance(stgComp,FeMesh_Type)
+ && (!strcmp(((FeMesh*)stgComp)->feElFamily,"linear")
+ || !strcmp(((FeMesh*)stgComp)->feElFamily,"quadratic"))) {
+ mesh = ( Mesh*)stgComp;
+ feMesh = (FeMesh*)stgComp;
- nDims = Mesh_GetDimSize( mesh );
- totalVerts = Mesh_GetGlobalSize( mesh, (MeshTopology_Dim)0 );
- elementGlobalSize = FeMesh_GetElementGlobalSize(mesh);
+ nDims = Mesh_GetDimSize( mesh );
+ totalVerts = Mesh_GetGlobalSize( mesh, (MeshTopology_Dim)0 );
+ elementGlobalSize = FeMesh_GetElementGlobalSize(mesh);
- /* get connectivity array size */
- if (mesh->nElTypes == 1)
- maxNodes = FeMesh_GetElementNodeSize( mesh, 0);
- else {
- /* determine the maximum number of nodes each element has */
- maxNodes = 0;
- for ( gElement_I = 0 ; gElement_I < FeMesh_GetElementGlobalSize(mesh); gElement_I++ ) {
- unsigned numNodes;
- numNodes = FeMesh_GetElementNodeSize( mesh, gElement_I);
- if( maxNodes < numNodes ) maxNodes = numNodes;
+ /* get connectivity array size */
+ if (mesh->nElTypes == 1)
+ maxNodes = FeMesh_GetElementNodeSize( mesh, 0);
+ else {
+ /* determine the maximum number of nodes each element has */
+ maxNodes = 0;
+ for ( gElement_I = 0 ; gElement_I < FeMesh_GetElementGlobalSize(mesh); gElement_I++ ) {
+ unsigned numNodes;
+ numNodes = FeMesh_GetElementNodeSize( mesh, gElement_I);
+ if( maxNodes < numNodes ) maxNodes = numNodes;
+ }
+ }
+ /** now write all the xdmf geometry info **/
+ /**----------------------- START GEOMETRY ------------------------------------------------------------------------------------------------------------------- **/
+ if(nDims==2 && (maxNodes==4 || maxNodes==9)){
+ Stg_asprintf( &topologyType, "Quadrilateral" );
+ } else if(nDims==3 && (maxNodes==8 || maxNodes==27)) {
+ Stg_asprintf( &topologyType, "Hexahedron" );
+ } else {
+ Journal_DPrintf( errorStream, "\n\n Error: number of element nodes %u not supported by XDMF generator...\n should be 4 or 9 (2D quadrilateral) "
+ "or should be 8 or 27 (3D hexahedron). \n\n", maxNodes );
+ Stg_asprintf( &topologyType, "UNKNOWN_POSSIBLY_ERROR" );
+ }
+ /** first create the grid which is applicable to the checkpointed fevariables **/
+ Journal_Printf( stream, " <Grid Name=\"FEM_Grid_%s\">\n\n", feMesh->name);
+ Journal_Printf( stream, " <Time Value=\"%f\" />\n\n", context->currentTime );
+ /** now print out topology info, only quadrilateral elements are supported at the moment **/
+
+ if(maxNodes==9 || maxNodes==27)
+ {
+ for(uint i=0;i<maxNodes;++i)
+ {
+ Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"C%u\">\n",
+ elementGlobalSize, i);
+ Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n",
+ i, elementGlobalSize);
+ Journal_Printf(stream," <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u %u\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n",
+ elementGlobalSize, maxNodes, feMesh->name, context->timeStep );
+ Journal_Printf(stream," </DataItem>\n");
}
- }
- /** now write all the xdmf geometry info **/
- /**----------------------- START GEOMETRY ------------------------------------------------------------------------------------------------------------------- **/
- if( maxNodes == 4 ){
- Stg_asprintf( &topologyType, "Quadrilateral" );
- } else if ( maxNodes == 8 ) {
- Stg_asprintf( &topologyType, "Hexahedron" );
- } else {
- Journal_DPrintf( errorStream, "\n\n Error: number of element nodes %u not supported by XDMF generator...\n should be 4 (2D quadrilateral) "
- "or should be 8 (3D hexahedron). \n\n", maxNodes );
- Stg_asprintf( &topologyType, "UNKNOWN_POSSIBLY_ERROR" );
- }
- /** first create the grid which is applicable to the checkpointed fevariables **/
- Journal_Printf( stream, " <Grid Name=\"FEM_Grid_%s\">\n\n", feMesh->name);
- Journal_Printf( stream, " <Time Value=\"%f\" />\n\n", context->currentTime );
- /** now print out topology info, only quadrilateral elements are supported at the moment **/
- Journal_Printf( stream, " <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n", topologyType, elementGlobalSize );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u %u\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n", elementGlobalSize, maxNodes, feMesh->name, context->timeStep );
- Journal_Printf( stream, " </Topology>\n\n" );
- Journal_Printf( stream, " <Geometry Type=\"XYZ\">\n" );
+ int num_nodes=(maxNodes==9) ? 4 : 8;
+ Journal_Printf(stream," <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n",
+ topologyType,elementGlobalSize*num_nodes);
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u %u\" Function=\"JOIN(",
+ elementGlobalSize*num_nodes, num_nodes);
+ for(int j=0;j<num_nodes-1;++j)
+ {
+ Journal_Printf(stream,"$%u; ",j);
+ }
+ Journal_Printf(stream,"$%u)\">\n",num_nodes-1);
- Stg_asprintf( &variableType, "NumberType=\"Float\" Precision=\"8\"" );
- if( nDims == 2 ){
- /** note that for 2d, we feed back a quasi 3d array, with the 3rd Dof zeroed. so in effect we always work in 3d.
- this is done because paraview in particular seems to do everything in 3d, and if you try and give it a 2d vector
- or array, it complains.... and hence the verbosity of the following 2d definitions**/
- Journal_Printf( stream, " <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0, $1, 0*$1)\">\n", totalVerts );
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"XCoords\">\n", totalVerts );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 0 1 1 %u 1 </DataItem>\n", totalVerts );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u 2\">Mesh.%s.%05d.h5:/vertices</DataItem>\n", variableType, totalVerts, feMesh->name, context->timeStep );
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"YCoords\">\n", totalVerts );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 1 1 1 %u 1 </DataItem>\n", totalVerts );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u 2\">Mesh.%s.%05d.h5:/vertices</DataItem>\n", variableType, totalVerts, feMesh->name, context->timeStep );
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " </DataItem>\n" );
- } else if ( nDims == 3 ) {
- /** in 3d we simply feed back the 3d hdf5 array, nice and easy **/
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u 3\">Mesh.%s.%05d.h5:/vertices</DataItem>\n", variableType, totalVerts, feMesh->name, context->timeStep );
- } else {
- Journal_DPrintf( errorStream, "\n\n Error: Mesh vertex location is not of dofCount 2 or 3.\n\n" );
- }
- Journal_Printf( stream, " </Geometry>\n\n" );
- /**----------------------- FINISH GEOMETRY ------------------------------------------------------------------------------------------------------------------- **/
+ int first_mapping[8]={0,1,4,3,9,10,13,12};
+ int mapping[8][8];
+ for(int i=0;i<8;++i)
+ {
+ mapping[0][i]=first_mapping[i];
+ mapping[1][i]=first_mapping[i]+1;
+ mapping[2][i]=first_mapping[i]+3;
+ mapping[3][i]=first_mapping[i]+4;
+ mapping[4][i]=first_mapping[i]+9;
+ mapping[5][i]=first_mapping[i]+10;
+ mapping[6][i]=first_mapping[i]+12;
+ mapping[7][i]=first_mapping[i]+13;
+ }
+
+ for(int j=0;j<num_nodes;++j)
+ {
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u %u\" Function=\"JOIN(",
+ elementGlobalSize, num_nodes);
+ for(int k=0;k<num_nodes-1;++k)
+ {
+ Journal_Printf(stream,"$%u, ",k);
+ }
+ Journal_Printf(stream,"$%u)\">\n",num_nodes-1);
+ for(int k=0;k<num_nodes;++k)
+ {
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C%u\"] </DataItem>\n",
+ feMesh->name, mapping[j][k]);
+ }
+ Journal_Printf(stream," </DataItem>\n");
+ }
+ Journal_Printf(stream," </DataItem>\n");
+ Journal_Printf(stream," </Topology>\n\n" );
+ }
+ else
+ {
+ Journal_Printf( stream, " <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n", topologyType, elementGlobalSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u %u\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n", elementGlobalSize, maxNodes, feMesh->name, context->timeStep );
+ Journal_Printf( stream, " </Topology>\n\n" );
+ }
+ Journal_Printf( stream, " <Geometry Type=\"XYZ\">\n" );
+
+ Stg_asprintf( &variableType, "NumberType=\"Float\" Precision=\"8\"" );
+ if( nDims == 2 ){
+ /** note that for 2d, we feed back a quasi 3d array, with the 3rd Dof zeroed. so in effect we always work in 3d.
+ this is done because paraview in particular seems to do everything in 3d, and if you try and give it a 2d vector
+ or array, it complains.... and hence the verbosity of the following 2d definitions**/
+ Journal_Printf( stream, " <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0, $1, 0*$1)\">\n", totalVerts );
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"XCoords\">\n", totalVerts );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 0 1 1 %u 1 </DataItem>\n", totalVerts );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u 2\">Mesh.%s.%05d.h5:/vertices</DataItem>\n", variableType, totalVerts, feMesh->name, context->timeStep );
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"YCoords\">\n", totalVerts );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 1 1 1 %u 1 </DataItem>\n", totalVerts );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u 2\">Mesh.%s.%05d.h5:/vertices</DataItem>\n", variableType, totalVerts, feMesh->name, context->timeStep );
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </DataItem>\n" );
+ } else if ( nDims == 3 ) {
+ /** in 3d we simply feed back the 3d hdf5 array, nice and easy **/
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u 3\">Mesh.%s.%05d.h5:/vertices</DataItem>\n", variableType, totalVerts, feMesh->name, context->timeStep );
+ } else {
+ Journal_DPrintf( errorStream, "\n\n Error: Mesh vertex location is not of dofCount 2 or 3.\n\n" );
+ }
+ Journal_Printf( stream, " </Geometry>\n\n" );
+ /**----------------------- FINISH GEOMETRY ------------------------------------------------------------------------------------------------------------------- **/
/** now write FeVariable data **/
for ( var_I = 0; var_I < context->fieldVariable_Register->objects->count; var_I++ ) {
- fieldVar = FieldVariable_Register_GetByIndex( context->fieldVariable_Register, var_I );
+ fieldVar = FieldVariable_Register_GetByIndex( context->fieldVariable_Register, var_I );
- if ( Stg_Class_IsInstance( fieldVar, FeVariable_Type ) ) {
- feVar = (FeVariable*)fieldVar;
- if ( (feVar->isCheckpointedAndReloaded && (context->timeStep % context->checkpointEvery == 0)) ||
- (feVar->isCheckpointedAndReloaded && (context->checkpointAtTimeInc && (context->currentTime >= context->nextCheckpointTime))) ||
- (feVar->isSavedData && (context->timeStep % context->saveDataEvery == 0)) ){
- FeMesh* feVarMesh = NULL;
- /** check what type of generator was used to know where elementMesh is **/
- if( Stg_Class_IsInstance( feVar->feMesh->generator, C0Generator_Type)) feVarMesh = (FeMesh*)(((C0Generator*)feVar->feMesh->generator)->elMesh);
- if( Stg_Class_IsInstance( feVar->feMesh->generator, CartesianGenerator_Type)) feVarMesh = feVar->feMesh;
- if( Stg_Class_IsInstance( feVar->feMesh->generator, MeshAdaptor_Type)) feVarMesh = feVar->feMesh;
- /** make sure that the fevariable femesh is the same as that used above for the geometry definition, if so proceed **/
- if( feVarMesh == feMesh ){
- char* centering = NULL;
- Index offset = 0;
- Index meshSize = Mesh_GetGlobalSize( feVar->feMesh, (MeshTopology_Dim)0 );
- Index dofCountIndex;
- Index dofAtEachNodeCount;
+ if ( Stg_Class_IsInstance( fieldVar, FeVariable_Type ) ) {
+ feVar = (FeVariable*)fieldVar;
+ if ( (feVar->isCheckpointedAndReloaded && (context->timeStep % context->checkpointEvery == 0)) ||
+ (feVar->isCheckpointedAndReloaded && (context->checkpointAtTimeInc && (context->currentTime >= context->nextCheckpointTime))) ||
+ (feVar->isSavedData && (context->timeStep % context->saveDataEvery == 0)) ){
+ FeMesh* feVarMesh = NULL;
+ /** check what type of generator was used to know where elementMesh is **/
+ if( Stg_Class_IsInstance( feVar->feMesh->generator, C0Generator_Type)) feVarMesh = (FeMesh*)(((C0Generator*)feVar->feMesh->generator)->elMesh);
+ if( Stg_Class_IsInstance( feVar->feMesh->generator, CartesianGenerator_Type)) feVarMesh = feVar->feMesh;
+ if( Stg_Class_IsInstance( feVar->feMesh->generator, MeshAdaptor_Type)) feVarMesh = feVar->feMesh;
+ /** make sure that the fevariable femesh is the same as that used above for the geometry definition, if so proceed **/
+ if( feVarMesh == feMesh ){
+ char* centering = NULL;
+ Index offset = 0;
+ Index meshSize = Mesh_GetGlobalSize( feVar->feMesh, (MeshTopology_Dim)0 );
+ Index dofCountIndex;
+ Index dofAtEachNodeCount;
- /**----------------------- START ATTRIBUTES ------------------------------------------------------------------------------------------------------------------- **/
- /** if coordinates are being stored with feVariable, account for this **/
- if( saveCoords) offset = nDims;
- /** all feVariables are currently stored as doubles **/
- Stg_asprintf( &variableType, "NumberType=\"Float\" Precision=\"8\"" );
- /** determine whether feVariable data is cell centered (like Pressure), or on the nodes (like Velocity) **/
- if( meshSize == elementGlobalSize ){
- Stg_asprintf( ¢ering, "Cell" );
- } else if( meshSize == totalVerts ) {
- Stg_asprintf( ¢ering, "Node" );
- } else {
- /* unknown/unsupported type */
- Stg_asprintf( ¢ering, "UNKNOWN_POSSIBLY_ERROR" );
- }
- /** how many degrees of freedom does the fevariable have? **/
- dofAtEachNodeCount = feVar->fieldComponentCount;
- if ( dofAtEachNodeCount == 1 ) {
- Journal_Printf( stream, " <Attribute Type=\"Scalar\" Center=\"%s\" Name=\"%s\">\n", centering, feVar->name);
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" >\n", meshSize );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", offset, meshSize );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " </Attribute>\n\n" );
- } else if ( dofAtEachNodeCount == 2 ){
- /** note that for 2d, we feed back a quasi 3d array, with the 3rd Dof zeroed. so in effect we always work in 3d.
- this is done because paraview in particular seems to do everything in 3d, and if you try and give it a 2d vector
- or array, it complains.... and hence the verbosity of the following 2d definitions **/
- Journal_Printf( stream, " <Attribute Type=\"Vector\" Center=\"%s\" Name=\"%s\">\n", centering, feVar->name);
- Journal_Printf( stream, " <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0, $1, 0*$1)\">\n", meshSize );
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"XValue\">\n", meshSize );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", offset, meshSize );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"YValue\">\n", meshSize );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", (offset+1), meshSize );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " </Attribute>\n\n" );
- } else if ( dofAtEachNodeCount == 3 ) {
- /** in 3d we simply feed back the 3d hdf5 array, nice and easy **/
- Journal_Printf( stream, " <Attribute Type=\"Vector\" Center=\"%s\" Name=\"%s\">\n", centering, feVar->name);
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 3\" >\n", meshSize );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 3 </DataItem>\n", offset, meshSize );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " </Attribute>\n\n" );
- } else {
- /** where there are more than 3 components, we write each one out as a scalar **/
- for(dofCountIndex = 0 ; dofCountIndex < dofAtEachNodeCount ; ++dofCountIndex){
- Journal_Printf( stream, " <Attribute Type=\"Scalar\" Center=\"%s\" Name=\"%s-Component-%u\">\n", centering, feVar->name, dofCountIndex);
- Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" >\n", meshSize );
- Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", (offset+dofCountIndex), meshSize );
- Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
- Journal_Printf( stream, " </DataItem>\n" );
- Journal_Printf( stream, " </Attribute>\n" );
- }
- Journal_Printf( stream, "\n" );
- }
- /**----------------------- END ATTRIBUTES ------------------------------------------------------------------------------------------------------------------- **/
- if(centering) Memory_Free( centering );
- }
- }
- }
+ /**----------------------- START ATTRIBUTES ------------------------------------------------------------------------------------------------------------------- **/
+ /** if coordinates are being stored with feVariable, account for this **/
+ if( saveCoords) offset = nDims;
+ /** all feVariables are currently stored as doubles **/
+ Stg_asprintf( &variableType, "NumberType=\"Float\" Precision=\"8\"" );
+ /** determine whether feVariable data is cell centered (like Pressure), or on the nodes (like Velocity) **/
+ if( meshSize == elementGlobalSize ){
+ Stg_asprintf( ¢ering, "Cell" );
+ } else if( meshSize == totalVerts ) {
+ Stg_asprintf( ¢ering, "Node" );
+ } else {
+ /* unknown/unsupported type */
+ Stg_asprintf( ¢ering, "UNKNOWN_POSSIBLY_ERROR" );
+ }
+ /** how many degrees of freedom does the fevariable have? **/
+ dofAtEachNodeCount = feVar->fieldComponentCount;
+ if ( dofAtEachNodeCount == 1 ) {
+ Journal_Printf( stream, " <Attribute Type=\"Scalar\" Center=\"%s\" Name=\"%s\">\n", centering, feVar->name);
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" >\n", meshSize );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", offset, meshSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </Attribute>\n\n" );
+ } else if ( dofAtEachNodeCount == 2 ){
+ /** note that for 2d, we feed back a quasi 3d array, with the 3rd Dof zeroed. so in effect we always work in 3d.
+ this is done because paraview in particular seems to do everything in 3d, and if you try and give it a 2d vector
+ or array, it complains.... and hence the verbosity of the following 2d definitions **/
+ Journal_Printf( stream, " <Attribute Type=\"Vector\" Center=\"%s\" Name=\"%s\">\n", centering, feVar->name);
+ Journal_Printf( stream, " <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0, $1, 0*$1)\">\n", meshSize );
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"XValue\">\n", meshSize );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", offset, meshSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"YValue\">\n", meshSize );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", (offset+1), meshSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </Attribute>\n\n" );
+ } else if ( dofAtEachNodeCount == 3 ) {
+ /** in 3d we simply feed back the 3d hdf5 array, nice and easy **/
+ Journal_Printf( stream, " <Attribute Type=\"Vector\" Center=\"%s\" Name=\"%s\">\n", centering, feVar->name);
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 3\" >\n", meshSize );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 3 </DataItem>\n", offset, meshSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </Attribute>\n\n" );
+ } else {
+ /** where there are more than 3 components, we write each one out as a scalar **/
+ for(dofCountIndex = 0 ; dofCountIndex < dofAtEachNodeCount ; ++dofCountIndex){
+ Journal_Printf( stream, " <Attribute Type=\"Scalar\" Center=\"%s\" Name=\"%s-Component-%u\">\n", centering, feVar->name, dofCountIndex);
+ Journal_Printf( stream, " <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" >\n", meshSize );
+ Journal_Printf( stream, " <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 %u 1 1 %u 1 </DataItem>\n", (offset+dofCountIndex), meshSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" %s Dimensions=\"%u %u\">%s.%05d.h5:/data</DataItem>\n", variableType, meshSize, (offset + dofAtEachNodeCount), feVar->name, context->timeStep);
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </Attribute>\n" );
+ }
+ Journal_Printf( stream, "\n" );
+ }
+ /**----------------------- END ATTRIBUTES ------------------------------------------------------------------------------------------------------------------- **/
+ if(centering) Memory_Free( centering );
+ }
+ }
+ }
}
- Journal_Printf( stream, " </Grid>\n\n" );
+ Journal_Printf( stream, " </Grid>\n\n" );
if(variableType) Memory_Free( variableType );
if(topologyType) Memory_Free( topologyType );
- }
- }
+ }
+ }
}
void _XDMFGenerator_WriteSwarmSchema( UnderworldContext* context, Stream* stream ) {
More information about the CIG-COMMITS
mailing list