[cig-commits] r7759 - in long/3D/Gale/trunk: .
src/Underworld/plugins/Output/VTKOutput
walter at geodynamics.org
walter at geodynamics.org
Fri Jul 27 22:26:04 PDT 2007
Author: walter
Date: 2007-07-27 22:26:03 -0700 (Fri, 27 Jul 2007)
New Revision: 7759
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c
Log:
r1906 at earth: boo | 2007-07-27 22:21:55 -0700
Make the vtk xml output work in parallel with fields. May not work if visualizing with multiple processors
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1900
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1906
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c 2007-07-27 23:33:47 UTC (rev 7758)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c 2007-07-28 05:26:03 UTC (rev 7759)
@@ -367,20 +367,26 @@
/* Print out the coordinates of the mesh. */
-void VTKOutput_print_coords(FILE *fp, FeMesh *feMesh, int nDims, int n) {
- Node_LocalIndex lNode_I;
+void VTKOutput_print_coords(FILE *fp, FeMesh *feMesh, Grid *grid, int nDims,
+ int lower[3], int upper[3]) {
+ IJK ijk;
fprintf(fp," <Points>\n");
fprintf(fp," <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n");
- for ( lNode_I = 0; lNode_I < FeMesh_GetNodeLocalSize(feMesh);
- lNode_I++ ) {
- double *coord = Mesh_GetVertex(feMesh, lNode_I);
- if(nDims==2)
- fprintf(fp, "%.15g %.15g 0\n", coord[0], coord[1]);
- else
- fprintf(fp, "%.15g %.15g %.15g\n", coord[0],
- coord[1], coord[2]);
- }
+
+ for(ijk[2]=lower[2];ijk[2]<upper[2];++ijk[2])
+ for(ijk[1]=lower[1];ijk[1]<upper[1];++ijk[1])
+ for(ijk[0]=lower[0];ijk[0]<upper[0];++ijk[0])
+ {
+ double *coord;
+ unsigned local;
+ Mesh_GlobalToDomain(feMesh,MT_VERTEX,Grid_Project(grid,ijk),&local);
+ coord=Mesh_GetVertex(feMesh,local);
+ if(nDims==2)
+ fprintf(fp, "%.15g %.15g 0\n", coord[0], coord[1]);
+ else
+ fprintf(fp, "%.15g %.15g %.15g\n", coord[0], coord[1], coord[2]);
+ }
fprintf(fp," </DataArray>\n </Points>\n");
}
@@ -394,11 +400,16 @@
Index var_I;
int header_printed=0;
int nDims, pn, i;
- int nx, ny, nz, n, pnx, pny, pnz;
+ int lower[3], upper[3], p_lower[3], p_upper[3];
Name field_filename, pressure_filename;
FILE *fp, *field_fp, *pressure_fp, *pfp, *pfield_fp, *ppressure_fp;
+ /* We need to save the grids to be used later when printing out the
+ pressure coordinates and to map between 1D and 3D indices for the
+ values. */
+ Grid *elGrid, *vertGrid;
+
/* Open the file */
Stg_asprintf( &pressure_filename, "%s/pressure.%d.%05d.vts",
@@ -439,6 +450,9 @@
FeVariable* feVar;
Node_LocalIndex lNode_I;
Dof_Index dofAtEachNodeCount;
+ int *low, *up;
+ Grid *grid;
+ IJK ijk;
feVar=(FeVariable*)fieldVar;
if(!header_printed)
@@ -453,64 +467,87 @@
if(!strcmp(gen->type,"SurfaceAdaptor"))
gen=(CartesianGenerator *)((SurfaceAdaptor *)(gen))->generator;
+ elGrid=gen->elGrid;
+ vertGrid=gen->vertGrid;
nDims=gen->nDims;
- nx=pnx=gen->range[0];
- if(gen->origin[0]==0)
- ++nx;
+ p_lower[2]=lower[2]=0;
+ p_upper[2]=upper[2]=1;
+ for(i=0;i<nDims;++i)
+ {
+ p_lower[i]=gen->origin[i];
+ p_upper[i]=p_lower[i]+gen->range[i];
- ny=pny=gen->range[1];
- if(gen->origin[1]==0)
- ++ny;
+ /* The grid is split up by elements, so for the pressure
+ grid, we simply add the ghost zones. Ghost zones are
+ only added at the top, not the bottom. */
+ if(p_upper[i]!=gen->elGrid->sizes[i])
+ ++p_upper[i];
- if(nDims==3)
- {
- nz=pnz=gen->range[2];
- if(gen->origin[2]==0)
- ++nz;
+ /* Then we get the vertices by adding appropriate
+ offsets to the element grid. */
+ upper[i]=p_upper[i]+1;
+ lower[i]=p_lower[i];
+ if(p_lower[i]!=0)
+ --lower[i];
}
- else
- nz=pnz=1;
- n=nx*ny*nz;
- pn=pnx*pny*pnz;
-
/* Now that we have the extents, write the header. */
fprintf(field_fp,"<?xml version=\"1.0\"?>\n\
<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n\
- <StructuredGrid WholeExtent=\"0 %d 0 %d 0 %d\">\n\
- <Piece Extent=\"0 %d 0 %d 0 %d\">\n\
- <CellData></CellData>\n",nx-1,ny-1,nz-1,nx-1,ny-1,nz-1);
+ <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n\
+ <Piece Extent=\"%d %d %d %d %d %d\">\n\
+ <CellData></CellData>\n",
+ lower[0],upper[0]-1,lower[1],upper[1]-1,lower[2],upper[2]-1,
+ lower[0],upper[0]-1,lower[1],upper[1]-1,lower[2],upper[2]-1);
fprintf(pressure_fp,"<?xml version=\"1.0\"?>\n\
<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n\
- <StructuredGrid WholeExtent=\"0 %d 0 %d 0 %d\">\n\
- <Piece Extent=\"0 %d 0 %d 0 %d\">\n\
- <CellData></CellData>\n",pnx-1,pny-1,pnz-1,pnx-1,pny-1,pnz-1);
+ <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n\
+ <Piece Extent=\"%d %d %d %d %d %d\">\n\
+ <CellData></CellData>\n",
+ p_lower[0],p_upper[0]-1,p_lower[1],p_upper[1]-1,
+ p_lower[2],p_upper[2]-1,
+ p_lower[0],p_upper[0]-1,p_lower[1],p_upper[1]-1,
+ p_lower[2],p_upper[2]-1);
/* Write the coordinates for the fields, but not the pressure */
- VTKOutput_print_coords(field_fp, feVar->feMesh, nDims, n);
+ VTKOutput_print_coords(field_fp, feVar->feMesh, gen->vertGrid, nDims,
+ lower, upper);
fprintf(field_fp," <PointData Scalars=\"StrainRateInvariantField\" Vectors=\"VelocityField\" Tensors=\"Stress\">\n");
/* Write out parallel control file headers. */
if(myRank==0)
{
- fprintf(pfield_fp,"<?xml version=\"1.0\"?>\n\
+ int global[3], p_global[3];
+ global[0]=gen->vertGrid->sizes[0];
+ global[1]=gen->vertGrid->sizes[1];
+ if(nDims==3)
+ {
+ global[2]=gen->elGrid->sizes[2];
+ p_global[2]=global[2]-1;
+ }
+ else
+ {
+ global[2]=p_global[2]=1;
+ }
+ fprintf(pfield_fp,"<?xml version=\"1.0\"?>\n\
<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n\
<PStructuredGrid GhostLevel=\"0\" WholeExtent=\"0 %d 0 %d 0 %d\">\n\
<PCellData></PCellData>\n\
<PPoints>\n\
<PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n\
</PPoints>\n\
- <PPointData Scalars=\"StrainRateInvariantField\" Vectors=\"VelocityField\" Tensors=\"Stress\">\n",nx-1,ny-1,nz-1);
- fprintf(ppressure_fp,"<?xml version=\"1.0\"?>\n\
-<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n\
+ <PPointData Scalars=\"StrainRateInvariantField\" Vectors=\"VelocityField\" Tensors=\"Stress\">\n",global[0]-1,global[1]-1,global[2]-1);
+ fprintf(ppressure_fp,"<?xml version=\"1.0\"?>\n\
+<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n \
<PStructuredGrid GhostLevel=\"0\" WholeExtent=\"0 %d 0 %d 0 %d\">\n\
<PCellData></PCellData>\n\
<PPoints>\n\
<PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n\
</PPoints>\n\
- <PPointData Scalars=\"PressureField\">\n",pnx-1,pny-1,pnz-1);
+ <PPointData Scalars=\"PressureField\">\n",
+ global[0]-2,global[1]-2,p_global[2]-1);
}
header_printed=1;
}
@@ -518,13 +555,20 @@
/* Write the coordinates for the pressure, and set the
appropriate file pointer to output for this variable. */
if(!strcmp(feVar->name,"PressureField")){
- VTKOutput_print_coords(pressure_fp, feVar->feMesh, nDims, pn);
+ VTKOutput_print_coords(pressure_fp, feVar->feMesh, elGrid, nDims,
+ p_lower, p_upper);
fprintf(pressure_fp," <PointData Scalars=\"PressureField\">\n");
fp=pressure_fp;
pfp=ppressure_fp;
+ low=p_lower;
+ up=p_upper;
+ grid=elGrid;
} else {
fp=field_fp;
pfp=pfield_fp;
+ low=lower;
+ up=upper;
+ grid=vertGrid;
}
/* Finally, output the fields. For now, just output every field */
@@ -558,55 +602,63 @@
break;
/* Unknown */
default:
- printf("Bad number of degrees of freedom for variable %s: %d\n",
+ fprintf(stderr,
+ "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*nDims)
- {
- /* If writing scalars or 3D objects, then just write the
- values. Otherwise, we have to fill in the 3D
- components with zeros. */
- case 2:
- case 3:
- case 9:
- case 27:
- for ( dof_I = 0; dof_I < dofAtEachNodeCount; dof_I++ ) {
- fprintf(fp, "%.15g ", variableValues[dof_I] );
+
+ for(ijk[2]=low[2];ijk[2]<up[2];++ijk[2])
+ for(ijk[1]=low[1];ijk[1]<up[1];++ijk[1])
+ for(ijk[0]=low[0];ijk[0]<up[0];++ijk[0])
+ {
+ double variableValues[MAX_FIELD_COMPONENTS];
+ Dof_Index dof_I;
+ unsigned local;
+ Mesh_GlobalToDomain(feVar->feMesh,MT_VERTEX,
+ Grid_Project(grid,ijk),&local);
+ FeVariable_GetValueAtNode(feVar,local,variableValues);
+
+ switch(dofAtEachNodeCount*nDims)
+ {
+ /* If writing scalars or 3D objects, then just write
+ the values. Otherwise, we have to fill in the 3D
+ components with zeros. */
+ case 2:
+ case 3:
+ case 9:
+ case 27:
+ for ( dof_I = 0; dof_I < dofAtEachNodeCount; dof_I++ ) {
+ fprintf(fp, "%.15g ", variableValues[dof_I] );
+ }
+ break;
+ case 4:
+ fprintf(fp, "%.15g %.15g 0", variableValues[0],
+ variableValues[1] );
+ break;
+ case 6:
+ /* Ordering is xx, yy, xy */
+ fprintf(fp, "%.15g %.15g 0 %.15g %.15g 0 0 0 0 ",
+ variableValues[0], variableValues[2],
+ variableValues[2], variableValues[1]);
+ break;
+ case 8:
+ fprintf(fp, "%.15g %.15g 0 %.15g %.15g 0 0 0 0 ",
+ variableValues[0], variableValues[1],
+ variableValues[2], variableValues[3]);
+ break;
+ case 18:
+ /* Ordering is xx, yy, zz, xy, xz, yz */
+ fprintf(fp, "%.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g ",
+ variableValues[0], variableValues[3],
+ variableValues[4], variableValues[3],
+ variableValues[1], variableValues[5],
+ variableValues[4], variableValues[5],
+ variableValues[2]);
+ break;
+ }
+ fprintf(fp, "\n" );
}
- break;
- case 4:
- fprintf(fp, "%.15g %.15g 0", variableValues[0],
- variableValues[1] );
- break;
- case 6:
- /* Ordering is xx, yy, xy */
- fprintf(fp, "%.15g %.15g 0 %.15g %.15g 0 0 0 0 ",
- variableValues[0], variableValues[2],
- variableValues[2], variableValues[1]);
- break;
- case 8:
- fprintf(fp, "%.15g %.15g 0 %.15g %.15g 0 0 0 0 ",
- variableValues[0], variableValues[1],
- variableValues[2], variableValues[3]);
- break;
- case 18:
- /* Ordering is xx, yy, zz, xy, xz, yz */
- fprintf(fp, "%.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g ",
- variableValues[0], variableValues[3],
- variableValues[4], variableValues[3],
- variableValues[1], variableValues[5],
- variableValues[4], variableValues[5],
- variableValues[2]);
- break;
- }
- fprintf(fp, "\n" );
- }
fprintf(fp," </DataArray>\n");
}
}
@@ -625,12 +677,14 @@
fprintf(pfield_fp," </PPointData>\n");
for(i=0;i<nprocs;++i)
{
+ fprintf(pfield_fp," <Piece Extent=\"%d %d %d %d %d %d\"\n\
+ Source=\"fields.%d.%05d.vts\"/>\n",lower[0],upper[0]-1,
+ lower[1],upper[1]-1,lower[2],upper[2]-1,
+ i,self->timeStep);
fprintf(ppressure_fp," <Piece Extent=\"%d %d %d %d %d %d\"\n\
- Source=\"pressure.%d.%05d.vts\"/>\n",0,pnx-1,0,pny-1,0,pnz-1,
+ Source=\"pressure.%d.%05d.vts\"/>\n",p_lower[0],p_upper[0]-1,
+ p_lower[1],p_upper[1]-1,p_lower[2],p_upper[2]-1,
i,self->timeStep);
- fprintf(pfield_fp," <Piece Extent=\"%d %d %d %d %d %d\"\n\
- Source=\"fields.%d.%05d.vts\"/>\n",0,nx-1,0,ny-1,0,nz-1,
- i,self->timeStep);
}
fprintf(ppressure_fp," </PStructuredGrid>\n\
</VTKFile>\n");
More information about the cig-commits
mailing list