[cig-commits] r6006 - in long/3D/Gale/trunk: .
src/Underworld/plugins/Output/VTKOutput
walter at geodynamics.org
walter at geodynamics.org
Sun Feb 11 00:51:37 PST 2007
Author: walter
Date: 2007-02-11 00:51:37 -0800 (Sun, 11 Feb 2007)
New Revision: 6006
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c
Log:
r1547 at earth: boo | 2007-02-11 00:48:52 -0800
Make VTK output work (but not in parallel)
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1546
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1547
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c 2007-02-11 08:51:32 UTC (rev 6005)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c 2007-02-11 08:51:37 UTC (rev 6006)
@@ -80,6 +80,7 @@
void VTKOutput_particles(IntegrationPointsSwarm* picswarm, int stepping,
char *outputPath, int timeStep, int dim);
+void VTKOutput_fields(void *context);
void VTKOutput( void* _context ) {
UnderworldContext* context = (UnderworldContext*)_context;
@@ -96,7 +97,7 @@
(dictionary,"particleStepping",1),
context->outputPath, context->timeStep,
context->dim);
-/* VTKOutput_fields(); */
+ VTKOutput_fields(context);
}
void VTKOutput_particles(IntegrationPointsSwarm* picswarm, int stepping,
@@ -192,7 +193,7 @@
switch(iteration)
{
case 0:
- Journal_Printf(stream,"POINT_DATA %d\nSCALARS Viscosity float 1\nLOOKUP_TABLE default\n",
+ Journal_Printf(stream,"\nPOINT_DATA %d\nSCALARS Viscosity float 1\nLOOKUP_TABLE default\n",
(num_particles-1)/stepping+1);
break;
case 1:
@@ -205,3 +206,179 @@
Stream_CloseFile( stream );
Memory_Free( filename );
}
+
+/* Print out the coordinates of the mesh. */
+
+void VTKOutput_print_coords(Stream *stream, FeMesh *feMesh, int nDims, int n) {
+ Node_LocalIndex lNode_I;
+ for ( lNode_I = 0; lNode_I < FeMesh_GetNodeLocalSize(feMesh);
+ lNode_I++ ) {
+ double *coord = Mesh_GetVertex(feMesh, lNode_I);
+ if(nDims==2)
+ Journal_Printf( stream, "%.15g %.15g 0\n", coord[0], coord[1]);
+ else
+ Journal_Printf( stream, "%.15g %.15g %.15g\n", coord[0],
+ coord[1], coord[2]);
+ }
+ Journal_Printf(stream,"POINT_DATA %d\n",n);
+}
+
+
+/* Pressure is stored on cell centers, while everything else is stored
+ on cell vertices. So we have to make two different files, one for
+ pressure, and one for everything else. */
+void VTKOutput_fields(void *context) {
+
+ FiniteElementContext* self = (FiniteElementContext*) context;
+ Index var_I;
+ int header_printed=0;
+ int nDims, pn;
+
+ Name field_filename, pressure_filename;
+ Stream* field_stream;
+ Stream* pressure_stream;
+ Stream* stream;
+
+ /* Open the file */
+
+ field_stream = Journal_Register( MPIStream_Type, "fields" );
+ pressure_stream = Journal_Register( MPIStream_Type, "pressure" );
+ Stg_asprintf( &pressure_filename, "%s/pressure.%05d.vtk", self->outputPath,
+ self->timeStep);
+ Stg_asprintf( &field_filename, "%s/fields.%05d.vtk", self->outputPath,
+ self->timeStep);
+
+ Stream_RedirectFile( field_stream, field_filename );
+ Stream_RedirectFile( pressure_stream, pressure_filename );
+
+ /* First, output the coordinates. We have to do a huge song and
+ dance just to get the extents of the mesh. */
+
+ for ( var_I = 0; var_I < self->fieldVariable_Register->objects->count;
+ var_I++ ) {
+ FieldVariable* fieldVar;
+ fieldVar = FieldVariable_Register_GetByIndex( self->fieldVariable_Register,
+ var_I );
+ if (Stg_Class_IsInstance( fieldVar, FeVariable_Type )) {
+ FeVariable* feVar;
+ Node_LocalIndex lNode_I;
+ Dof_Index dofAtEachNodeCount;
+
+ feVar=(FeVariable*)fieldVar;
+ if(!header_printed)
+ {
+ int nx, ny, nz, n, pnx, pny, pnz;
+ Grid *grid;
+ Mesh *mesh;
+
+ mesh=(Mesh*)(feVar->feMesh);
+ grid=*(Grid**)ExtensionManager_Get
+ ( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
+ nx=grid->sizes[0];
+ ny=grid->sizes[1];
+ if(grid->nDims==3)
+ nz=grid->sizes[2];
+ else
+ nz=1;
+ n=grid->nPoints;
+
+ pnx=(nx==1 ? 1 : nx-1);
+ pny=(ny==1 ? 1 : ny-1);
+ pnz=(nz==1 ? 1 : nz-1);
+ pn=pnx*pny*pnz;
+
+ /* Now that we have the extents, write the header. */
+
+ Journal_Printf(field_stream,"# vtk DataFile Version 2.0\nGale\nASCII\nDATASET STRUCTURED_GRID\nDIMENSIONS %d %d %d\nPOINTS %d float\n",nx,ny,nz,n);
+ Journal_Printf(pressure_stream,"# vtk DataFile Version 2.0\nGale\nASCII\nDATASET STRUCTURED_GRID\nDIMENSIONS %d %d %d\nPOINTS %d float\n",pnx,pny,pnz,pn);
+
+ /* Write the coordinates */
+ nDims=grid->nDims;
+ VTKOutput_print_coords(field_stream, feVar->feMesh, nDims, n);
+ header_printed=1;
+ }
+
+ /* Write the coordinates for the pressure, and set the
+ appropriate stream to output for this variable. */
+ if(!strcmp(feVar->name,"PressureField")){
+ VTKOutput_print_coords(pressure_stream, feVar->feMesh, nDims, pn);
+ stream=pressure_stream;
+ } else {
+ stream=field_stream;
+ }
+
+ /* Finally, output the fields. For now, just output every field */
+
+ dofAtEachNodeCount = feVar->fieldComponentCount;
+
+ switch(dofAtEachNodeCount)
+ {
+ /* Scalars */
+ case 1:
+ Journal_Printf(stream,"SCALARS %s float 1\nLOOKUP_TABLE default\n",
+ feVar->name);
+ break;
+ /* Vectors */
+ case 2:
+ case 3:
+ Journal_Printf(stream,"VECTORS %s float\n",feVar->name);
+ break;
+ /* Rank 2 Tensors */
+ case 4:
+ case 6:
+ case 9:
+ Journal_Printf(stream,"TENSORS %s float\n",feVar->name);
+ break;
+ /* Unknown */
+ default:
+ printf("Bad number of degrees of freedom for variable %s: %d\n",
+ feVar->name, dofAtEachNodeCount);
+ abort();
+ }
+ for ( lNode_I = 0; lNode_I < FeMesh_GetNodeLocalSize(feVar->feMesh);
+ lNode_I++ ) {
+ double variableValues[MAX_FIELD_COMPONENTS];
+ Dof_Index dof_I;
+ FeVariable_GetValueAtNode( feVar, lNode_I, variableValues );
+ switch(dofAtEachNodeCount)
+ {
+ /* If writing scalars or 3D objects, then just write the
+ values. Otherwise, we have to fill in the 3D
+ components with zeros. */
+ case 1:
+ case 3:
+ case 9:
+ for ( dof_I = 0; dof_I < dofAtEachNodeCount; dof_I++ ) {
+ Journal_Printf(stream, "%.15g ", variableValues[dof_I] );
+ }
+ break;
+ case 2:
+ Journal_Printf(stream, "%.15g %.15g 0", variableValues[0],
+ variableValues[1] );
+ break;
+ case 4:
+ Journal_Printf(stream, "%.15g %.15g 0 %.15g %.15g 0 0 0 0 ",
+ variableValues[0], variableValues[1],
+ variableValues[2], variableValues[3]);
+ break;
+ case 6:
+ Journal_Printf(stream, "%.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g ",
+ variableValues[0], variableValues[1],
+ variableValues[2], variableValues[1],
+ variableValues[3], variableValues[4],
+ variableValues[2], variableValues[4],
+ variableValues[5]);
+ break;
+ }
+ Journal_Printf(stream, "\n" );
+ }
+ }
+ }
+ Stream_CloseFile( pressure_stream );
+ Memory_Free( pressure_filename );
+ Stream_CloseFile( field_stream );
+ Memory_Free( field_filename );
+}
+
More information about the cig-commits
mailing list