[cig-commits] commit: Properly implement XDMF support for discontinuous elements in 2D. 3D still just makes tetrahedra.
Mercurial
hg at geodynamics.org
Fri Oct 7 23:43:21 PDT 2011
changeset: 884:60cca4257758
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Oct 07 23:41:41 2011 -0700
files: Utils/src/XDMFGenerator.cxx
description:
Properly implement XDMF support for discontinuous elements in 2D. 3D still just makes tetrahedra.
diff -r 64ca21dfcaed -r 60cca4257758 Utils/src/XDMFGenerator.cxx
--- a/Utils/src/XDMFGenerator.cxx Fri Oct 07 23:40:49 2011 -0700
+++ b/Utils/src/XDMFGenerator.cxx Fri Oct 07 23:41:41 2011 -0700
@@ -50,6 +50,8 @@
#include <stdlib.h>
#include <string.h>
+#include <string>
+#include <sstream>
#ifndef MASTER
#define MASTER 0
@@ -135,20 +137,45 @@ namespace {
}
}
+namespace {
+ void write_continuous_geometry(Stream *stream, FeMesh *feMesh,
+ const std::string &element_family,
+ const int ×tep,
+ const int &nDims,
+ const std::string &topologyType);
+ void write_continuous_FeVariable(Stream *stream,
+ FeMesh *feMesh,
+ FeVariable *feVar,
+ const int ×tep,
+ const int &nDims,
+ Bool saveCoords);
+ void write_discontinuous_geometry(Stream *stream, Mesh *mesh,
+ const int ×tep,
+ const int &nDims);
+ void write_discontinuous_FeVariable(Stream *stream,
+ const std::string &mesh_name,
+ const std::string &var_name,
+ const int ×tep,
+ const int &num_elements);
+ void pack_together_rows(Stream *stream, const std::string &mesh_name,
+ const std::string &coord_name,
+ const int &number_in_element,
+ unsigned *sizes);
+ int pack_together_rows(Stream *stream, const std::string &mesh_name,
+ const std::string &coord_name,
+ const int &number_in_element,
+ const std::string &prefix,
+ const int &start,
+ const int &end,
+ unsigned *sizes);
+}
+
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;
std::string topologyType;
unsigned componentCount = LiveComponentRegister_GetCount(stgLiveComponentRegister);
unsigned compI;
@@ -167,24 +194,9 @@ void _XDMFGenerator_WriteFieldSchema( Un
|| element_family=="quadratic") {
mesh = ( Mesh*)stgComp;
feMesh = (FeMesh*)stgComp;
+ const unsigned nDims=Mesh_GetDimSize( mesh );
+ int elementGlobalSize = FeMesh_GetElementGlobalSize(feMesh);
- 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;
- }
- }
/** now write all the xdmf geometry info **/
/**----------------------- START GEOMETRY --------------------------- **/
if(element_family=="linear" || element_family=="quadratic")
@@ -201,250 +213,16 @@ void _XDMFGenerator_WriteFieldSchema( Un
Journal_Printf( stream, " <Time Value=\"%f\" />\n\n", context->currentTime );
/** now print out topology info **/
- /* Quadratic elements */
- if(element_family=="quadratic")
+ /* Discontinuous elements are written separately. */
+ if(element_family=="linear-inner" && nDims==2)
{
- /* First separate each node of each element into its own
- variable C0-C8 or C0-C26 */
- 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");
- }
- /* Print out the topology */
- int num_nodes=(maxNodes==9) ? 4 : 8;
- Journal_Printf(stream," <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n",
- topologyType.c_str(),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);
-
- 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;
- }
-
- /* Print out the quadrilaterals made up of only some of the
- nodes */
- 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" );
+ write_discontinuous_geometry(stream,mesh,context->timeStep,nDims);
}
-
- /* Discontinuous linear elements */
- else if(element_family=="linear-inner")
- {
- if(nDims==2)
- {
- Grid** grid=
- (Grid** )ExtensionManager_Get
- (mesh->info, mesh,
- ExtensionManager_GetHandle(mesh->info, (Name)"elementGrid"));
-
- unsigned *sizes=Grid_GetSizes(*grid);
-
- /* Create a bunch of subsets of the connectivity points.
- The range operator [a:b] does not seem to work with
- Paraview, so we have to manually create some
- subsets. */
- /* Create C0, a set of the points of the first node in
- each element */
- Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"C0\">\n",
- elementGlobalSize);
- Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 0 1 1 %u 1 </DataItem>\n",
- elementGlobalSize);
- Journal_Printf(stream," <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u 3\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n",
- elementGlobalSize, feMesh->name, context->timeStep);
- Journal_Printf(stream," </DataItem>\n");
-
- /* Create C0x, a set of the points of the first node in each
- element except for the elements on the right. */
- /* First create a strips that go from the left to the
- right, stopping one element short of the edge. */
- for(uint i=0;i<sizes[1];++i)
- {
- Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"C0x%u\">\n",
- sizes[0]-1,i);
- Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> %u 0 1 1 %u 1 </DataItem>\n",
- i*sizes[0],sizes[0]-1);
- Journal_Printf(stream," <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u 3\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n",
- elementGlobalSize, feMesh->name, context->timeStep);
- Journal_Printf(stream," </DataItem>\n");
- }
-
- /* Now join the strips together to create C0x */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u\" Function=\"JOIN(",
- (sizes[0]-1)*sizes[1]);
- for(uint i=0;i<sizes[1]-1;++i)
- Journal_Printf(stream,"$%u; ",i);
- Journal_Printf(stream,"$%u)\" Name=\"C0x\">\n",
- sizes[1]-1);
-
- for(uint i=0;i<sizes[1];++i)
- {
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0x%u\"]</DataItem>\n",
- feMesh->name,i);
- }
- Journal_Printf(stream," </DataItem>\n");
-
- /* Use C0 and C0x to create all of the triangles */
- int total_triangles=elementGlobalSize + 2*(elementGlobalSize-sizes[0])
- + 2*(elementGlobalSize-sizes[1]-sizes[0]+1)
- + elementGlobalSize-sizes[0]
- + 2*(sizes[0]-1);
- Journal_Printf(stream," <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n",
- topologyType.c_str(),
- total_triangles);
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0; $1; $2; $3; $4; $5; $6; $7)\">\n",
- total_triangles);
- /* Triangles for the three points in the element */
- Journal_Printf(stream," <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u 3\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n",
- elementGlobalSize, feMesh->name,
- context->timeStep );
- /* The numbering for xdmf is really wacked. In JOIN($0
- +1, $0+2), the second value will actually be $0+3,
- because it will add 1 and use that as the current $0.
- Very annoying. */
- /* Triangle for nodes 1,2 in the bottom element and 1 in
- the top element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN(1 + $0, 1 + $0, %u + $0)\">\n",
- elementGlobalSize-sizes[0],3*sizes[0]-1);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0\"] </DataItem>\n",
- feMesh->name);
- Journal_Printf(stream," </DataItem>\n" );
-
- /* Triangle for nodes 0,2 in the bottom element and 0 in
- the top element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0, 2 + $0, %u + $0)\">\n",
- elementGlobalSize-sizes[0],3*sizes[0]-2);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0\"] </DataItem>\n",
- feMesh->name);
- Journal_Printf(stream," </DataItem>\n" );
-
- /* Triangle for nodes 1 in the left element, 0 in the
- right element, and 1 in the top element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN(1 + $0, 2 + $0, %u + $0)\">\n",
- elementGlobalSize-sizes[1]-sizes[0]+1,3*sizes[0]-2);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0x\"] </DataItem>\n",
- feMesh->name);
- Journal_Printf(stream," </DataItem>\n" );
-
- /* Triangle for nodes 0 in the right element, 1 in the
- top element, and 0 in the top right element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN(3 + $0, %u + $0, 2 + $0)\">\n",
- elementGlobalSize-sizes[1]-sizes[0]+1,3*sizes[0]-3+1);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0x\"] </DataItem>\n",
- feMesh->name);
- Journal_Printf(stream," </DataItem>\n" );
-
- /* Triangle for nodes 2 in the left element and nodes
- 0,1 in the top element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN(2 + $0, %u + $0, 1 + $0)\">\n",
- elementGlobalSize-sizes[0],3*sizes[0]-2);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0\"] </DataItem>\n",
- feMesh->name);
- Journal_Printf(stream," </DataItem>\n" );
-
- /* Triangle for the strip on the top side. It is made
- up of nodes 1,2 in the left element and node 0 in
- the right element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN(%u + $0, 1 + $0, 1 + $0)\">\n",
- sizes[0]-1,(sizes[1]-1)*3*sizes[0]+1);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0x\"] </DataItem>\n",
- feMesh->name);
- Journal_Printf(stream," </DataItem>\n" );
-
- /* Triangle for the strip on the right side. It is made
- up of node 2 in the left element and nodes 0,2 in
- the right element. */
- Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN(%u + $0, 1 + $0, 2 + $0)\">\n",
- sizes[0]-1,(sizes[1]-1)*3*sizes[0]+2);
- Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"C0x\"] </DataItem>\n",
- feMesh->name);
- 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.c_str(), 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" );
- }
- }
- /* Continuous linear elements */
else
{
- Journal_Printf( stream, " <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n", topologyType.c_str(), 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" );
+ write_continuous_geometry(stream,feMesh,element_family,
+ context->timeStep,nDims,topologyType);
}
- Journal_Printf( stream, " <Geometry Type=\"XYZ\">\n" );
-
-
- /* Print out coordinates of the vertices */
- 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;
@@ -479,81 +257,26 @@ void _XDMFGenerator_WriteFieldSchema( Un
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" );
+ if(element_family=="linear-inner" && nDims==2)
+ {
+ write_discontinuous_FeVariable(stream,feMesh->name,feVar->name,
+ context->timeStep,
+ elementGlobalSize);
}
- Journal_Printf( stream, "\n" );
- }
- /**----------------------- END ATTRIBUTES ------------------------------------------------------------------------------------------------------------------- **/
- if(centering) Memory_Free( centering );
+ else
+ {
+ write_continuous_FeVariable
+ (stream,feMesh,feVar,context->timeStep,nDims,
+ Dictionary_GetBool_WithDefault
+ (context->dictionary,
+ (Dictionary_Entry_Key)"saveCoordsWithFields",False));
+
+ }
}
}
}
}
- Journal_Printf( stream, " </Grid>\n\n" );
- if(variableType) Memory_Free( variableType );
+ Journal_Printf( stream, " </Grid>\n\n" );
}
}
}
@@ -712,7 +435,7 @@ void _XDMFGenerator_WriteSwarmSchema( Un
}
- Journal_Printf( stream, " </Grid>\n\n" );
+ Journal_Printf( stream, " </Grid>\n\n" );
if(variableType) Memory_Free( variableType );
if(filename_part) Memory_Free( filename_part );
@@ -769,3 +492,511 @@ void _XDMFGenerator_SendInfo( Underworld
}
+namespace {
+
+ void write_continuous_geometry(Stream *stream, FeMesh *feMesh,
+ const std::string &element_family,
+ const int ×tep,
+ const int &nDims,
+ const std::string &topologyType)
+ {
+ int totalVerts = Mesh_GetGlobalSize( feMesh, (MeshTopology_Dim)0 );
+ int elementGlobalSize = FeMesh_GetElementGlobalSize(feMesh);
+
+ uint maxNodes;
+ /* get connectivity array size */
+ if (feMesh->nElTypes == 1)
+ maxNodes = FeMesh_GetElementNodeSize( feMesh, 0);
+ else {
+ /* determine the maximum number of nodes each element has */
+ maxNodes = 0;
+ for (uint gElement_I = 0; gElement_I < FeMesh_GetElementGlobalSize(feMesh);
+ gElement_I++ ) {
+ uint numNodes;
+ numNodes = FeMesh_GetElementNodeSize( feMesh, gElement_I);
+ if( maxNodes < numNodes ) maxNodes = numNodes;
+ }
+ }
+ /* Quadratic elements */
+ if(element_family=="quadratic")
+ {
+ /* First separate each node of each element into its own
+ variable C0-C8 or C0-C26 */
+ 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, timestep );
+ Journal_Printf(stream," </DataItem>\n");
+ }
+ /* Print out the topology */
+ int num_nodes=(maxNodes==9) ? 4 : 8;
+ Journal_Printf(stream," <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n",
+ topologyType.c_str(),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);
+
+ 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;
+ }
+
+ /* Print out the quadrilaterals made up of only some of the
+ nodes */
+ 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" );
+ }
+
+ /* Continuous linear elements */
+ else
+ {
+ Journal_Printf( stream, " <Topology Type=\"%s\" NumberOfElements=\"%u\"> \n",
+ topologyType.c_str(), elementGlobalSize );
+ Journal_Printf( stream, " <DataItem Format=\"HDF\" DataType=\"Int\" Dimensions=\"%u %u\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n",
+ elementGlobalSize, maxNodes, feMesh->name,
+ timestep );
+ Journal_Printf( stream, " </Topology>\n\n" );
+ }
+
+ /* Print out coordinates of the vertices */
+ Journal_Printf( stream, " <Geometry Type=\"XYZ\">\n" );
+ std::string 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.c_str(), totalVerts, feMesh->name, 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.c_str(), totalVerts, feMesh->name, 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.c_str(), totalVerts, feMesh->name, timestep );
+ } else {
+ Stream* errorStream=Journal_Register(Error_Type,(Name)CURR_MODULE_NAME);
+ Journal_DPrintf( errorStream, "\n\n Error: Mesh vertex location is not of dofCount 2 or 3.\n\n" );
+ }
+ Journal_Printf( stream, " </Geometry>\n\n" );
+ }
+
+ void write_continuous_FeVariable(Stream *stream, FeMesh *feMesh,
+ FeVariable *feVar,
+ const int ×tep,
+ const int &nDims,
+ Bool saveCoords)
+ {
+ std::string centering;
+ Index offset = 0;
+ Index meshSize = Mesh_GetGlobalSize( feVar->feMesh, (MeshTopology_Dim)0 );
+
+ /**----------------------- START ATTRIBUTES ----------------- **/
+ /** if coordinates are being stored with feVariable,
+ account for this **/
+ if( saveCoords) offset = nDims;
+ /** all feVariables are currently stored as doubles **/
+ std::string variableType="NumberType=\"Float\" Precision=\"8\"";
+ /** determine whether feVariable data is cell centered (like
+ Pressure with P0 elements), or on the nodes (like
+ Velocity) **/
+ uint totalVerts(Mesh_GetGlobalSize(feMesh,(MeshTopology_Dim)0));
+ uint elementGlobalSize(FeMesh_GetElementGlobalSize(feMesh));
+ if(meshSize == elementGlobalSize){
+ centering="Cell";
+ } else if(meshSize == totalVerts) {
+ centering="Node";
+ } else {
+ /* unknown/unsupported type */
+ centering="UNKNOWN_POSSIBLY_ERROR";
+ }
+
+ /** how many degrees of freedom does the fevariable have? **/
+ uint dofAtEachNodeCount = feVar->fieldComponentCount;
+ if (dofAtEachNodeCount == 1) {
+ Journal_Printf( stream, " <Attribute Type=\"Scalar\" Center=\"%s\" Name=\"%s\">\n", centering.c_str(), 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.c_str(), meshSize, (offset + dofAtEachNodeCount), feVar->name, 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.c_str(), 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.c_str(), meshSize, (offset + dofAtEachNodeCount), feVar->name, 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.c_str(), meshSize, (offset + dofAtEachNodeCount), feVar->name, 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.c_str(), 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.c_str(), meshSize, (offset + dofAtEachNodeCount), feVar->name, 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(uint dofCountIndex = 0 ; dofCountIndex < dofAtEachNodeCount ; ++dofCountIndex){
+ Journal_Printf( stream, " <Attribute Type=\"Scalar\" Center=\"%s\" Name=\"%s-Component-%u\">\n", centering.c_str(), 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.c_str(), meshSize, (offset + dofAtEachNodeCount), feVar->name, timestep);
+ Journal_Printf( stream, " </DataItem>\n" );
+ Journal_Printf( stream, " </Attribute>\n" );
+ }
+ Journal_Printf( stream, "\n" );
+ }
+ }
+
+ void write_discontinuous_geometry(Stream *stream, Mesh *mesh,
+ const int ×tep,
+ const int &nDims)
+ {
+ if(nDims==2)
+ {
+ Grid** grid=
+ (Grid** )ExtensionManager_Get
+ (mesh->info, mesh,
+ ExtensionManager_GetHandle(mesh->info, (Name)"elementGrid"));
+
+ unsigned *sizes=Grid_GetSizes(*grid);
+ unsigned num_elements(sizes[0]*sizes[1]);
+ unsigned vert_sizes[]={2*sizes[0]+1,2*sizes[1]+1};
+ unsigned num_vertices=vert_sizes[0] * vert_sizes[1];
+
+ char
+ *vert_mesh_name(((InnerGenerator*)(mesh->generator))->elMesh->name);
+
+ /* We need a number that goes from 0 to num_elements. We hack
+ it with a WHERE clause that is always true */
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 1\" Function=\"WHERE( $0 != -1)\" Name=\"N\">\n",
+ num_elements);
+ Journal_Printf(stream," <DataItem Format=\"HDF\" NumberType=\"Float\" Precision=\"8\" Dimensions=\"%u 3\">Mesh.%s.%05d.h5:/connectivity</DataItem>\n",
+ num_elements,mesh->name,timestep);
+ Journal_Printf(stream," </DataItem>\n");
+
+ /* Write out the Quadrilaterals connectivity. There are 4 quads per element. */
+ Journal_Printf(stream," <Topology Type=\"Quadrilateral\" NumberOfElements=\"%u\">\n",
+ num_elements*4);
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 4\" Function=\"($0; $1; $2; $3)\">\n",
+ num_elements*4);
+ /* Quad 0 */
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 4\" Function=\"(9*$0), (1 + 9*$0), (4 + 9*$0), (3 + 9*$0)\">\n",
+ num_elements);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"N\"] </DataItem>\n",
+ mesh->name);
+ Journal_Printf(stream," </DataItem>\n");
+ /* Quad 1 */
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 4\" Function=\"(1 + 9*$0), (2 + 9*$0), (5 + 9*$0), (4 + 9*$0)\">\n",
+ num_elements);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"N\"] </DataItem>\n",
+ mesh->name);
+ Journal_Printf(stream," </DataItem>\n");
+ /* Quad 2 */
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 4\" Function=\"(3 + 9*$0), (4 + 9*$0), (7 + 9*$0), (6 + 9*$0)\">\n",
+ num_elements);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"N\"] </DataItem>\n",
+ mesh->name);
+ Journal_Printf(stream," </DataItem>\n");
+ /* Quad 3 */
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 4\" Function=\"(4 + 9*$0), (5 + 9*$0), (8 + 9*$0), (7 + 9*$0)\">\n",
+ num_elements);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"N\"] </DataItem>\n",
+ mesh->name);
+ Journal_Printf(stream," </DataItem>\n");
+ Journal_Printf(stream," </DataItem>\n");
+ Journal_Printf(stream," </Topology>\n");
+
+
+ /* Write out the coordinates. This is a bit complicated,
+ because we want separate vertices for each element, so the
+ sides of elements can not share. */
+
+ /* First get the coordinates from the vertex mesh */
+ Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"XCoords\">\n",
+ num_vertices);
+ Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 0 1 1 %u 1 </DataItem>\n",
+ num_vertices);
+ Journal_Printf(stream," <DataItem Format=\"HDF\" NumberType=\"Float\" Precision=\"8\" Dimensions=\"%u 2\">Mesh.%s.%05d.h5:/vertices</DataItem>\n",
+ num_vertices,vert_mesh_name,timestep);
+ Journal_Printf(stream," </DataItem>\n");
+
+ Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"YCoords\">\n",
+ num_vertices);
+ Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> 0 1 1 1 %u 1 </DataItem>\n",
+ num_vertices);
+ Journal_Printf(stream," <DataItem Format=\"HDF\" NumberType=\"Float\" Precision=\"8\" Dimensions=\"%u 2\">Mesh.%s.%05d.h5:/vertices</DataItem>\n",
+ num_vertices,vert_mesh_name,timestep);
+ Journal_Printf(stream," </DataItem>\n");
+
+ std::string direction[]={"X", "Y"};
+ for(int d=0;d<2;++d)
+ {
+ for(int number_in_element=0;number_in_element<9;++number_in_element)
+ {
+ pack_together_rows(stream,mesh->name,direction[d],
+ number_in_element,sizes);
+ }
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 1\" Function=\"$0, $1, $2, $3, $4, $5, $6, $7, $8\" Name=\"%s\">\n",
+ 9*num_elements,direction[d].c_str());
+ for(int i=0;i<9;++i)
+ {
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"%s%d\"] </DataItem>\n",
+ mesh->name,direction[d].c_str(),i);
+ }
+ Journal_Printf(stream," </DataItem>\n");
+ }
+ Journal_Printf(stream," <Geometry Type=\"XYZ\">\n");
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 3\" Function=\"JOIN($0, $1, 0*$1)\">\n",
+ 9*num_elements);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"X\"] </DataItem>\n",
+ mesh->name);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"Y\"] </DataItem>\n",
+ mesh->name);
+ Journal_Printf(stream," </DataItem>\n");
+ Journal_Printf(stream," </Geometry>\n");
+ }
+ /* nDim==3 */
+ else
+ {
+ abort();
+ }
+ }
+
+ void pack_together_rows(Stream *stream, const std::string &mesh_name,
+ const std::string &coord_name,
+ const int &number_in_element,
+ unsigned *sizes)
+ {
+ std::stringstream ss;
+ ss << coord_name
+ << number_in_element;
+ pack_together_rows(stream,mesh_name,coord_name,number_in_element,
+ ss.str(),0,sizes[1],sizes);
+ }
+
+ int pack_together_rows(Stream *stream, const std::string &mesh_name,
+ const std::string &coord_name,
+ const int &number_in_element,
+ const std::string &prefix,
+ const int &start,
+ const int &end,
+ unsigned *sizes)
+ {
+ int total_size(0);
+ int n_elements;
+ int row_start[]={0,1,2,2*sizes[0]+1,2*sizes[0]+2,2*sizes[0]+3,
+ 2*(2*sizes[0]+1),2*(2*sizes[0]+1)+1,2*(2*sizes[0]+1)+2};
+ if(end-start<10)
+ {
+ for(int i=start;i<end;++i)
+ {
+ Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"%s%u\">\n",
+ sizes[0],prefix.c_str(),i-start);
+ Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> %u 0 2 1 %u 1 </DataItem>\n",
+ i*(2*sizes[0]+1)*2 + row_start[number_in_element],
+ sizes[0]);
+ Journal_Printf(stream," <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_%s\"]/DataItem[@Name=\"%sCoords\"] </DataItem>\n",
+ mesh_name.c_str(),coord_name.c_str());
+ Journal_Printf(stream," </DataItem>\n");
+ }
+ n_elements=end-start;
+ total_size=sizes[0]*n_elements;
+ }
+ else
+ {
+ for(int i=0;i<10;++i)
+ {
+ int new_start=((end-start)*i)/10;
+ int new_end=((end-start)*(i+1))/10;
+ std::stringstream ss;
+ ss << prefix
+ << i;
+ total_size+=pack_together_rows(stream,mesh_name,coord_name,
+ number_in_element,
+ ss.str(),new_start,
+ new_end,sizes);
+ }
+ n_elements=10;
+ }
+ std::stringstream ss;
+ ss << " <DataItem ItemType=\"Function\" Dimensions=\""
+ << total_size
+ << " 1\" Function=\"";
+ for(int i=0;i<n_elements-1;++i)
+ ss << "$" << i << "; ";
+ ss << "$" << n_elements-1
+ << "\" Name=\""
+ << prefix
+ << "\">\n";
+ Journal_Printf(stream,ss.str().c_str());
+ for(int i=0;i<n_elements;++i)
+ {
+ std::stringstream s;
+ s << " <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_"
+ << mesh_name
+ << "\"]/DataItem[@Name=\""
+ << prefix
+ << i
+ << "\"] </DataItem>\n";
+ Journal_Printf(stream,s.str().c_str());
+ }
+ Journal_Printf(stream," </DataItem>\n");
+
+ return total_size;
+ }
+
+ void write_three_component_function(Stream *stream,
+ const std::string &start,
+ const std::string &end,
+ const int &num_elements,
+ const std::string &equation,
+ const std::string &prefix,
+ const std::string &result,
+ const std::string &first,
+ const std::string &second,
+ const std::string &third)
+ {
+ std::stringstream ss;
+ ss << " <DataItem ItemType=\"Function\" Dimensions=\""
+ << num_elements
+ << " 1\" Function=\""
+ << equation
+ << "\" Name=\""
+ << prefix
+ << result
+ << "\">\n";
+ Journal_Printf(stream, ss.str().c_str());
+
+ Journal_Printf(stream,(start + first + end).c_str());
+ Journal_Printf(stream,(start + second + end).c_str());
+ Journal_Printf(stream,(start + third + end).c_str());
+
+ Journal_Printf(stream," </DataItem>\n");
+ }
+
+ void write_discontinuous_FeVariable(Stream *stream,
+ const std::string &mesh_name,
+ const std::string &var_name,
+ const int ×tep,
+ const int &num_elements)
+ {
+ for(int i=0;i<3;++i)
+ {
+ Journal_Printf(stream," <DataItem ItemType=\"HyperSlab\" Dimensions=\"%u 1\" Name=\"%s%d\">\n",
+ num_elements,var_name.c_str(),i);
+ Journal_Printf(stream," <DataItem Dimensions=\"3 2\" Format=\"XML\"> %u 0 3 1 %u 1 </DataItem>\n",
+ i,num_elements);
+ Journal_Printf(stream," <DataItem Format=\"HDF\" NumberType=\"Float\" Precision=\"8\" Dimensions=\"%u 1\">%s.%05d.h5:/data</DataItem>\n",
+ 3*num_elements,var_name.c_str(),timestep);
+ Journal_Printf(stream," </DataItem>\n");
+ }
+
+ /* Write P_avg, dP/dx, dP/dy */
+ std::string start(" <DataItem Reference=\"XML\">/Xdmf/Domain/Grid[@Name=\"FEM_Grid_"
+ + mesh_name
+ + "\"]/DataItem[@Name=\""
+ + var_name);
+ std::string end("\"] </DataItem>\n");
+ write_three_component_function(stream,start,end,num_elements,
+ "($2 + ($1 + $0) / 2) / 2",
+ var_name,"_avg","0","1","2");
+ write_three_component_function(stream,start,end,num_elements,
+ "($1 - $0)",var_name,"_dx","0","1","2");
+ write_three_component_function(stream,start,end,num_elements,
+ "($2 - ($1 + $0) / 2)",
+ var_name,"_dy","0","1","2");
+
+ /* Write P at each vertex */
+ write_three_component_function(stream,start,end,num_elements,
+ "($0 - $1 - $2)",var_name,
+ "_0","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 - $2",var_name,"_1","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 + $1 - $2",var_name,
+ "_2","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 - $1",var_name,"_3","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0",var_name,"_4","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 + $1",var_name,"_5","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 - $1 + $2",var_name,
+ "_6","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 + $2",var_name,"_7","_avg","_dx","_dy");
+ write_three_component_function(stream,start,end,num_elements,
+ "$0 + $1 + $2",var_name,
+ "_8","_avg","_dx","_dy");
+
+ /* Bundle them all together for the field */
+ Journal_Printf(stream," <Attribute Type=\"Scalar\" Center=\"Node\" Name=\"%s\">\n",
+ var_name.c_str());
+ Journal_Printf(stream," <DataItem ItemType=\"Function\" Dimensions=\"%u 1\" Function=\"$0, $1, $2, $3, $4, $5, $6, $7, $8\">\n",
+ num_elements*9);
+ for(int i=0;i<9;++i)
+ {
+ std::stringstream ss;
+ ss << " "
+ << start
+ << "_"
+ << i
+ << end;
+ Journal_Printf(stream,ss.str().c_str());
+ }
+ Journal_Printf(stream," </DataItem>\n");
+ Journal_Printf(stream," </Attribute>\n");
+ }
+}
More information about the CIG-COMMITS
mailing list