[cig-commits] r14571 - in long/3D/SNAC/trunk/Snac: libSnac/src plugins/restarter snac2restart snac2vtk
echoi at geodynamics.org
echoi at geodynamics.org
Thu Apr 2 12:38:46 PDT 2009
Author: echoi
Date: 2009-04-02 12:38:45 -0700 (Thu, 02 Apr 2009)
New Revision: 14571
Modified:
long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c
long/3D/SNAC/trunk/Snac/snac2restart/snac2restart.c
long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c
Log:
1. Modified the way of writing stressTensor.
* stressTensor outputs now have 6 components due to symmetry.
* The ordering is 00,11,22,01,02,12.
- The reading and writing routines in
libSnac/src/Context.c,
snac2vtk/snac2vtk.c,
snac2restart/snac2restart.c
plugins/restarter/InitialConditions.c
have been modified accordingly.
2. phaseIndex is now written by a journal stream.
3. snac2restart.c improved
* Fixed a bug in snac2restart.c, which caused a failure in assigning a correct writePath in the case of argNum==3.
* The "usage" message has been made clearer.
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2009-04-02 19:18:25 UTC (rev 14570)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2009-04-02 19:38:45 UTC (rev 14571)
@@ -262,7 +262,7 @@
GetOffsetOfMember( tmpElement, material_I ),
GetOffsetOfMember( tmpElement, strainRate ),
GetOffsetOfMember( tmpElement, stress ),
- GetOffsetOfMember( tmpElement, hydroPressure )};
+ GetOffsetOfMember( tmpElement, hydroPressure )};
Variable_DataType elementDataTypes[] = {
Variable_DataType_Int,
Variable_DataType_Double,
@@ -601,7 +601,6 @@
fclose( self->timeStepInfo );
fclose( self->simInfo );
fclose( self->stressTensorOut );
- fclose( self->phaseIndexOut );
/* Parallelisation information */
if( self->parallel ) {
@@ -985,7 +984,7 @@
Journal_Dump( self->forceOut, &(self->timeStep) );
- _Snac_Context_DumpPhaseIndex( self );
+ Journal_Dump( self->phaseIndexOut, &(self->timeStep) );
KeyCall( self, self->syncK, EntryPoint_Class_VoidPtr_CallCast* )( KeyHandle(self,self->syncK), self );
/* _Snac_Context_Sync( self ); */
@@ -1064,7 +1063,7 @@
Mesh_Sync( self->mesh );
}
- _Snac_Context_DumpPhaseIndex( self );
+ Journal_Dump( self->phaseIndexOut, &(self->timeStep) );
}
@@ -1213,11 +1212,11 @@
Stream_Flush( self->strainRateOut );
Stream_Flush( self->stressOut );
Stream_Flush( self->hydroPressureOut );
+ Stream_Flush( self->phaseIndexOut );
Stream_Flush( self->coordOut );
Stream_Flush( self->velOut );
Stream_Flush( self->forceOut );
fflush( self->stressTensorOut );
- fflush( self->phaseIndexOut );
fflush( NULL );
}
@@ -1263,6 +1262,16 @@
self->dumpEvery,
tmpBuf );
+ /* Create the phase index dumping stream */
+ self->phaseIndexOut = Journal_Register( VariableDumpStream_Type, "PhaseIndex" );
+ sprintf( tmpBuf, "%s/phaseIndex.%u", self->outputPath, self->rank );
+ VariableDumpStream_SetVariable(
+ self->phaseIndexOut,
+ Variable_Register_GetByName( self->variable_Register, "elementMaterial" ),
+ self->mesh->elementLocalCount,
+ self->dumpEvery,
+ tmpBuf );
+
/* Create the coords dumping stream */
self->coordOut = Journal_Register( VariableDumpStream_Type, "Coord" );
sprintf( tmpBuf, "%s/coord.%u", self->outputPath, self->rank );
@@ -1298,11 +1307,6 @@
if( (self->stressTensorOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( self->stressTensorOut /* failed to open file for writing */ );
}
- /* Create the phaseIndex dump file */
- sprintf( tmpBuf, "%s/phaseIndex.%u", self->outputPath, self->rank );
- if( (self->phaseIndexOut = fopen( tmpBuf, "w+" )) == NULL ) {
- assert( self->phaseIndexOut /* failed to open file for writing */ );
- }
}
void _Snac_Context_AdjustDump( Snac_Context* self ) {
@@ -1366,34 +1370,20 @@
void _Snac_Context_DumpStressTensor( Snac_Context* self ) {
-
+
/* fprintf(stderr, "Dumping ? stress tensor: ts=%d df=%d\n",self->timeStep,self->dumpEvery); */
if( self->timeStep ==0 || (self->timeStep-1) % self->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
- fprintf(stderr, "r=%d, ts=%d/%d: Dumping stress tensor given dump freq=%d\n", self->rank, self->timeStep, self->maxTimeSteps,self->dumpEvery);
-
- for( element_lI = 0; element_lI < self->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( self, element_lI );
- /* Take average of tetra viscosity for the element */
- Tetrahedra_Index tetra_I;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- float tensor[3][3] = { {element->tetra[tetra_I].stress[0][0], element->tetra[tetra_I].stress[0][1], element->tetra[tetra_I].stress[0][2]},
- {element->tetra[tetra_I].stress[0][1], element->tetra[tetra_I].stress[1][1], element->tetra[tetra_I].stress[1][2]},
- {element->tetra[tetra_I].stress[0][2], element->tetra[tetra_I].stress[1][2], element->tetra[tetra_I].stress[2][2]} };
- fwrite( &tensor, sizeof(float), 9, self->stressTensorOut );
- }
- }
- }
-}
-
-void _Snac_Context_DumpPhaseIndex( Snac_Context* self ) {
-
- if( self->timeStep ==0 || (self->timeStep-1) % self->dumpEvery == 0 ) {
Element_LocalIndex element_lI;
-
+
for( element_lI = 0; element_lI < self->mesh->elementLocalCount; element_lI++ ) {
Snac_Element* element = Snac_Element_At( self, element_lI );
- fwrite( &element->material_I, sizeof(unsigned int), 1, self->phaseIndexOut );
+ /* Take average of tetra viscosity for the element */
+ Tetrahedra_Index tetra_I;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float stressVector[6] = { element->tetra[tetra_I].stress[0][0], element->tetra[tetra_I].stress[1][1], element->tetra[tetra_I].stress[2][2],
+ element->tetra[tetra_I].stress[0][1], element->tetra[tetra_I].stress[0][2], element->tetra[tetra_I].stress[1][2]};
+ fwrite( &stressVector, sizeof(float), 6, self->stressTensorOut );
+ }
}
- }
+ }
}
Modified: long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c 2009-04-02 19:18:25 UTC (rev 14570)
+++ long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c 2009-04-02 19:38:45 UTC (rev 14571)
@@ -139,7 +139,10 @@
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
double S[3][3];
Index i,j;
- fscanf( fp, "%le %le %le %le %le %le %le %le %le", &S[0][0],&S[0][1],&S[0][2],&S[1][0],&S[1][1],&S[1][2],&S[2][0],&S[2][1],&S[2][2]);
+ fscanf( fp, "%le %le %le %le %le %le", &S[0][0],&S[1][1],&S[2][2],&S[0][1],&S[0][2],&S[1][2]);
+ S[1][0] = S[0][1];
+ S[2][0] = S[0][2];
+ S[2][1] = S[1][2];
for(i=0;i<3;i++)
for(j=0;j<3;j++)
element->tetra[tetra_I].stress[i][j] = S[i][j];
Modified: long/3D/SNAC/trunk/Snac/snac2restart/snac2restart.c
===================================================================
--- long/3D/SNAC/trunk/Snac/snac2restart/snac2restart.c 2009-04-02 19:18:25 UTC (rev 14570)
+++ long/3D/SNAC/trunk/Snac/snac2restart/snac2restart.c 2009-04-02 19:38:45 UTC (rev 14571)
@@ -37,8 +37,9 @@
#include <assert.h>
#include <stdlib.h>
#include <string.h>
+#include <math.h>
+#include <limits.h>
-#include <limits.h>
#ifndef PATH_MAX
#define PATH_MAX 1024
#endif
@@ -46,7 +47,6 @@
#define Tetrahedra_Count 10
#endif
-
void ConvertTimeStep( int rank, unsigned int dumpIteration, unsigned int simTimeStep, double time );
char readPath[PATH_MAX];
@@ -67,56 +67,73 @@
const double velocityDampingFactor = 1.0;
+const unsigned numStressVectorComponent = 6; /* 6 components in stress vector */
+const unsigned numStressComponentsPerElement = 60; /* 6 components times 10 tets per element */
void checkArgumentType( int argNum, char* arglist[] ) {
if( argNum == 1 ) {
- sprintf( readPath, "." );
- sprintf( writePath, "." );
- fprintf(stderr,"Reading Snac state files for the last time step from %s/\n",readPath);
- fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
- return;
- }
- else if( argNum == 2 ) {
- if( atoi(arglist[1]) ) {
- restartStep = atoi(arglist[1]);
- sprintf( readPath, "." );
- sprintf( writePath, "." );
- fprintf(stderr,"Reading Snac state files for time step ts=%d from %s/\n",restartStep,readPath);
- fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
- return;
- }
- else {
- char zeroChar[PATH_MAX];
- sprintf(zeroChar,"0");
- if( !strcmp( arglist[1], zeroChar) ) {
- fprintf(stderr,"\nDon't try to restart from 0th time step. Just rerun the case!!\n\n");
- exit(0);
- }
- else {
- sprintf( readPath, "%s", arglist[1]);
+ sprintf( readPath, "." );
sprintf( writePath, "." );
fprintf(stderr,"Reading Snac state files for the last time step from %s/\n",readPath);
fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
return;
- }
- }
}
+ else if( argNum == 2 ) {
+ if( atoi(arglist[1]) ) {
+ restartStep = atoi(arglist[1]);
+ sprintf( readPath, "." );
+ sprintf( writePath, "." );
+ fprintf(stderr,"Reading Snac state files for time step ts=%d from %s/\n",restartStep,readPath);
+ fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
+ return;
+ }
+ else {
+ char zeroChar[PATH_MAX];
+ sprintf(zeroChar,"0");
+ if( !strcmp( arglist[1], zeroChar) ) {
+ fprintf(stderr,"\nDon't try to restart from 0th time step. Just rerun the case!!\n\n");
+ exit(0);
+ }
+ else {
+ sprintf( readPath, "%s", arglist[1]);
+ sprintf( writePath, "." );
+ fprintf(stderr,"Reading Snac state files for the last time step from %s/\n",readPath);
+ fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
+ return;
+ }
+ }
+ }
else if( argNum == 3 ) {
- if( atoi(arglist[1]) && !atoi(arglist[2]) ) {
- restartStep = atoi(arglist[1]);
- sprintf( readPath, "%s", arglist[2]);
- fprintf(stderr,"Reading Snac state files for time step ts=%d from %s/\n",restartStep,readPath);
- fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
- return;
- }
- else {
- fprintf(stderr,"Wrong argument type\n\tUsage: %s [integer timeStep] [input file path]\n",arglist[0]);
- exit(0);
- }
+ if( atoi(arglist[1]) && !atoi(arglist[2]) ) {
+ restartStep = atoi(arglist[1]);
+ sprintf( readPath, "%s", arglist[2]);
+ sprintf( writePath, "%s", arglist[2]);
+ fprintf(stderr,"Reading Snac state files for time step ts=%d from %s/\n",restartStep,readPath);
+ fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
+ return;
+ }
+ else {
+ fprintf(stderr,"Wrong argument type\n\tUsage: %s [integer timeStep] [your \"outputPath\"] [path to write restart files]\n",arglist[0]);
+ exit(0);
+ }
}
+ else if( argNum == 4 ) {
+ if( atoi(arglist[1]) && !atoi(arglist[2]) && !atoi(arglist[3])) {
+ restartStep = atoi(arglist[1]);
+ sprintf( readPath, "%s", arglist[2]);
+ sprintf( writePath, "%s", arglist[3]);
+ fprintf(stderr,"Reading Snac state files for time step ts=%d from %s/\n",restartStep,readPath);
+ fprintf(stderr,"Writing Snac restart files to %s/\n",writePath);
+ return;
+ }
+ else {
+ fprintf(stderr,"Wrong argument type\n\tUsage: %s [integer timeStep] [your \"outputPath\"] [path to write restart files]\n",arglist[0]);
+ exit(0);
+ }
+ }
else if( argNum > 3 ) {
- fprintf(stderr,"Wrong number of arguments\n\tUsage: %s [timeStep] [input file path]\n",arglist[0]);
- exit(0);
+ fprintf(stderr,"Wrong number of arguments\n\tUsage: %s [timeStep] [your \"outputPath\"] [path to write restart files]\n",arglist[0]);
+ exit(0);
}
}
@@ -142,13 +159,13 @@
fprintf(stderr,"Attempting to read from %s ...",tmpBuf);
if( (simIn = fopen( tmpBuf, "r" )) == NULL ) {
if( rank == 0 ) {
- /* failed to open file for reading */
- fprintf(stderr, " failed - no such file\n");
- exit(0);
+ /* failed to open file for reading */
+ fprintf(stderr, " failed - no such file\n");
+ exit(0);
}
else {
- fprintf(stderr," no such file\n");
- break;
+ fprintf(stderr,"no files to process for rank %d.\n(Make sure the max rank is indeed %d).\n", rank, rank-1);
+ break;
}
} else {
fprintf(stderr," with success\n");
@@ -269,14 +286,13 @@
unsigned int tetra_I;
float minLength;
-
/* Write out position array */
sprintf( tmpBuf, "%s/snac.coord.%i.%06u.restart", writePath, rank, simTimeStep );
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
- /* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
- exit(0);
+ /* failed to open file for writing */
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
+ exit(0);
}
fseek( coordIn, dumpIteration * nodeLocalCount * sizeof(float) * 3, SEEK_SET );
for( node_gI = 0; node_gI < nodeLocalCount; node_gI++ ) {
@@ -291,9 +307,9 @@
sprintf( tmpBuf, "%s/snac.initCoord.%i.%06u.restart", writePath, rank, simTimeStep );
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
- /* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
- exit(0);
+ /* failed to open file for writing */
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
+ exit(0);
}
fseek( coordIn, 0 * nodeLocalCount * sizeof(float) * 3, SEEK_SET );
for( node_gI = 0; node_gI < nodeLocalCount; node_gI++ ) {
@@ -308,16 +324,16 @@
sprintf( tmpBuf, "%s/snac.vel.%i.%06u.restart", writePath, rank, simTimeStep );
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
- /* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
- exit(0);
+ /* failed to open file for writing */
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
+ exit(0);
}
/*
* CPS hack: need to try damping/zeroing of node velocities in expts where we are
* restarting from assumed static elastic equilibrium
*
* The variable velocityDampingFactor should be passed from the command line
- * BEWARE that it is fixed at zero for now.
+ * BEWARE that it is fixed at 1.0 for now.
*/
fseek( velIn, dumpIteration * nodeLocalCount * sizeof(float) * 3, SEEK_SET );
for( node_gI = 0; node_gI < nodeLocalCount; node_gI++ ) {
@@ -333,9 +349,9 @@
sprintf( tmpBuf, "%s/snac.force.%i.%06u.restart", writePath, rank, simTimeStep );
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
- /* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
- exit(0);
+ /* failed to open file for writing */
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
+ exit(0);
}
for( node_gI = 0; node_gI < nodeLocalCount; node_gI++ ) {
float force;
@@ -350,24 +366,38 @@
sprintf( tmpBuf, "%s/snac.stressTensor.%i.%06u.restart", writePath, rank, simTimeStep );
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
- /* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
- exit(0);
+ /* failed to open file for writing */
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
+ exit(0);
}
- fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float) * 9 * Tetrahedra_Count, SEEK_SET );
+ fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float) * numStressComponentsPerElement, SEEK_SET );
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- float tensor[3][3];
- fread( &tensor, sizeof(float), 9, stressTensorIn );
- fprintf( restartOut, "%.9e %.9e %.9e %.9e %.9e %.9e %.9e %.9e %.9e\n",
- tensor[0][0], tensor[0][1], tensor[0][2],
- tensor[1][0], tensor[1][1], tensor[1][2],
- tensor[2][0], tensor[2][1], tensor[2][2] );
- fflush( restartOut );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float stressTensorArray[numStressVectorComponent];
+
+ if ( fread( stressTensorArray, sizeof(float), numStressVectorComponent, stressTensorIn )==0 ) {
+ if (feof(stressTensorIn)) {
+ fprintf(stderr, "Error (reached EOF prematurely) while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
+ exit(1);
+ } else {
+ fprintf(stderr, "Error while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
+ exit(1);
+ }
+ }
+ if(isnan(stressTensorArray[0]) || isnan(stressTensorArray[1])
+ || isnan(stressTensorArray[2]) || isnan(stressTensorArray[3])
+ || isnan(stressTensorArray[4]) || isnan(stressTensorArray[5])) {
+ fprintf(stderr,"NaN in stress tensor file\n");
+ abort();
+ }
+ fprintf( restartOut, "%.9e %.9e %.9e %.9e %.9e %.9e\n",
+ stressTensorArray[0],stressTensorArray[1],stressTensorArray[2],
+ stressTensorArray[3],stressTensorArray[4],stressTensorArray[5] );
+ fflush( restartOut );
+ }
}
- }
if( restartOut )
- fclose( restartOut );
+ fclose( restartOut );
/* Write out temperature array */
if(doTemp) {
@@ -375,7 +405,7 @@
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
/* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
exit(0);
}
@@ -425,7 +455,7 @@
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
/* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
exit(0);
}
@@ -445,7 +475,7 @@
fprintf(stderr,"\tWriting out %s\n",tmpBuf);
if( (restartOut = fopen( tmpBuf, "w" )) == NULL ) {
/* failed to open file for writing */
- fprintf(stderr, "Failed to open for writing\n");
+ fprintf(stderr, "Failed to open %s for writing\n", tmpBuf);
exit(0);
}
fscanf( minLengthIn, "%e", &minLength );
Modified: long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c
===================================================================
--- long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c 2009-04-02 19:18:25 UTC (rev 14570)
+++ long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c 2009-04-02 19:38:45 UTC (rev 14571)
@@ -53,8 +53,10 @@
#ifndef PATH_MAX
#define PATH_MAX 1024
#endif
+#ifndef Tetrahedra_Count
+ #define Tetrahedra_Count 10
+#endif
-
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
a[k][l]=h+s*(g-h*tau);
@@ -121,6 +123,8 @@
int doHPr = 1;
int doVisc = 1;
double failureAngle = 30.0;
+const unsigned numStressVectorComponent = 6; /* 6 components in stress vector */
+const unsigned numStressComponentsPerElement = 60; /* 6 components times 10 tets per element */
int main( int argc, char* argv[])
{
@@ -551,25 +555,25 @@
* Write out the pressure information
*/
fprintf( vtkOut, " <DataArray type=\"Float32\" Name=\"Pressure\" format=\"ascii\">\n");
- if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*9*10, SEEK_SET )!=0) {
- fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
- rank, rankI, rankJ, rankK,
- dumpIteration, nodeLocalCount );
- exit(1);
+ if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float) * numStressComponentsPerElement, SEEK_SET )!=0) {
+ fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
+ rank, rankI, rankJ, rankK,
+ dumpIteration, nodeLocalCount );
+ exit(1);
}
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
- struct stressMeasures elementStressMeasures;
- /*
- * Build average stress tensor for element and derive useful stress measures
- */
- DeriveStressMeasures(stressTensorIn, elementStressTensor, &elementStressMeasures);
- /*
- * Write the chosen derived stress value to Paraview file
- */
- fprintf( vtkOut, "%g ", -elementStressMeasures.pressure );
+ double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
+ struct stressMeasures elementStressMeasures;
+ /*
+ * Build average stress tensor for element and derive useful stress measures
+ */
+ DeriveStressMeasures(stressTensorIn, elementStressTensor, &elementStressMeasures);
+ /*
+ * Write the chosen derived stress value to Paraview file
+ */
+ fprintf( vtkOut, "%g ", -elementStressMeasures.pressure );
#ifdef DEBUG
- if (element_gI<10) fprintf( stderr, "Element pressure %d: %g at angle %g\n", element_gI, elementStressMeasures.pressure, elementStressMeasures.failureAngle);
+ if (element_gI<10) fprintf( stderr, "Element pressure %d: %g at angle %g\n", element_gI, elementStressMeasures.pressure, elementStressMeasures.failureAngle);
#endif
}
fprintf( vtkOut, " </DataArray>\n");
@@ -622,7 +626,7 @@
fprintf(stderr, "Error (reached EOF prematurely) while reading Snac shear stress output file: rank=%d: %d, %d, %d, dump iteration=%d, node=%d/%d\n",
rank, rankI, rankJ, rankK,
dumpIteration, node_gI, nodeLocalCount );
- exit(1);
+ exit(1);
} else {
fprintf(stderr, "Error while reading Snac shear stress output file: rank=%d: %d, %d, %d, dump iteration=%d, node=%d/%d\n",
rank, rankI, rankJ, rankK,
@@ -639,14 +643,14 @@
* Write out the maxShearStress information
*/
fprintf( vtkOut, " <DataArray type=\"Float32\" Name=\"Max. shear stress\" format=\"ascii\">\n");
- if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*9*10, SEEK_SET )!=0) {
+ if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*numStressComponentsPerElement, SEEK_SET )!=0) {
fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
rank, rankI, rankJ, rankK,
dumpIteration, nodeLocalCount );
exit(1);
}
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
+ double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
struct stressMeasures elementStressMeasures;
/*
* Build average stress tensor for element and derive useful stress measures
@@ -667,14 +671,14 @@
* Write out the vonMisesStress information
*/
fprintf( vtkOut, " <DataArray type=\"Float32\" Name=\"Von Mises stress\" format=\"ascii\">\n");
- if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*9*10, SEEK_SET )!=0) {
+ if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*numStressComponentsPerElement, SEEK_SET )!=0) {
fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
rank, rankI, rankJ, rankK,
dumpIteration, nodeLocalCount );
exit(1);
}
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
+ double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
struct stressMeasures elementStressMeasures;
/*
* Build average stress tensor for element and derive useful stress measures
@@ -694,14 +698,14 @@
* Write out the slopeShearStress information
*/
fprintf( vtkOut, " <DataArray type=\"Float32\" Name=\"Shear stress @%gd\" format=\"ascii\">\n",failureAngle);
- if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*9*10, SEEK_SET )!=0) {
+ if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*numStressComponentsPerElement, SEEK_SET )!=0) {
fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
rank, rankI, rankJ, rankK,
dumpIteration, nodeLocalCount );
exit(1);
}
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
+ double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
struct stressMeasures elementStressMeasures;
/*
* Build average stress tensor for element and derive useful stress measures
@@ -722,14 +726,14 @@
* Write out the slopeNormalStress information
*/
fprintf( vtkOut, " <DataArray type=\"Float32\" Name=\"Normal stress @%gd\" format=\"ascii\">\n", failureAngle);
- if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*9*10, SEEK_SET )!=0) {
+ if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*numStressComponentsPerElement, SEEK_SET )!=0) {
fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
rank, rankI, rankJ, rankK,
dumpIteration, nodeLocalCount );
exit(1);
}
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
+ double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
struct stressMeasures elementStressMeasures;
/*
* Build average stress tensor for element and derive useful stress measures
@@ -748,14 +752,14 @@
/* Write out the failurePotential information */
fprintf( vtkOut, " <DataArray type=\"Float32\" Name=\"Failure potential\" format=\"ascii\">\n");
- if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*9*10, SEEK_SET )!=0) {
+ if (fseek( stressTensorIn, dumpIteration * elementLocalCount * sizeof(float)*numStressComponentsPerElement, SEEK_SET )!=0) {
fprintf(stderr, "Cannot find read required portion of Snac stress tensor output file: rank=%d: %d, %d, %d, dump iteration=%d, node count=%d\n",
rank, rankI, rankJ, rankK,
dumpIteration, nodeLocalCount );
exit(1);
}
for( element_gI = 0; element_gI < elementLocalCount; element_gI++ ) {
- double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
+ double elementStressTensor[3][3]={{0,0,0},{0,0,0},{0,0,0}};
struct stressMeasures elementStressMeasures;
/*
* Build average stress tensor for element and derive useful stress measures
@@ -776,7 +780,6 @@
-
/*
* Write out the phase information
*/
@@ -936,6 +939,7 @@
fprintf( vtkOut1, " <PDataArray type=\"Float32\" Name=\"Shear stress @%gd\"/>\n",failureAngle);
fprintf( vtkOut1, " <PDataArray type=\"Float32\" Name=\"Normal stress @%gd\"/>\n",failureAngle);
fprintf( vtkOut1, " <PDataArray type=\"Float32\" Name=\"Failure potential @%gd\"/>\n",failureAngle);
+ fprintf( vtkOut1, " <PDataArray type=\"Float32\" Name=\"Phase\"/>\n");
if( doVisc )
fprintf( vtkOut1, " <PDataArray type=\"Float32\" Name=\"Viscosity\"/>\n");
@@ -1229,45 +1233,47 @@
DeriveStressMeasures(FILE *stressTensorIn, double elementStressTensor[3][3], struct stressMeasures *elementStressMeasures)
{
int tetra_I;
- const int numberTetrahedra=10;
double failureAngleRadians=elementStressMeasures->failureAngle*M_PI/180.0;
double normalVector[3],slopeParallelVector[3],tractionVector[3];
double tmp;
- /*
- * Read all tetrahedral stress tensors for this element gI
- */
- for( tetra_I = 0; tetra_I < 10; tetra_I++ ) {
- float stressTensorArray[3][3];
- if ( fread( stressTensorArray, sizeof(float), 9, stressTensorIn )==0 ) {
- if (feof(stressTensorIn)) {
- fprintf(stderr, "Error (reached EOF prematurely) while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
- exit(1);
- } else {
- fprintf(stderr, "Error while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
- exit(1);
- }
- }
- /*
- * Report error and bail if we pick up NaNs in any of the stress components
- */
- if(isnan(stressTensorArray[0][0]) || isnan(stressTensorArray[1][1])
- || isnan(stressTensorArray[2][2]) || isnan(stressTensorArray[0][1])
- || isnan(stressTensorArray[0][2]) || isnan(stressTensorArray[1][2]))
- fprintf(stderr,"NaN in stress tensor file\n");
- /*
- * Build average stress tensor for element by summing tetrahedral tensor components
- * - even though it's symmetric, do for all 9 components in case we pick the wrong ones before diagonalization
- */
- elementStressTensor[0][0]+=stressTensorArray[0][0]/(double)numberTetrahedra;
- elementStressTensor[1][1]+=stressTensorArray[1][1]/(double)numberTetrahedra;
- elementStressTensor[2][2]+=stressTensorArray[2][2]/(double)numberTetrahedra;
- elementStressTensor[0][1]+=stressTensorArray[0][1]/(double)numberTetrahedra;
- elementStressTensor[0][2]+=stressTensorArray[0][2]/(double)numberTetrahedra;
- elementStressTensor[1][2]+=stressTensorArray[1][2]/(double)numberTetrahedra;
- elementStressTensor[1][0]+=stressTensorArray[1][0]/(double)numberTetrahedra;
- elementStressTensor[2][0]+=stressTensorArray[2][0]/(double)numberTetrahedra;
- elementStressTensor[2][1]+=stressTensorArray[2][1]/(double)numberTetrahedra;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float stressTensorArray[numStressVectorComponent];
+
+ /*
+ * Read all tetrahedral stress tensors for this element gI
+ */
+ if ( fread( stressTensorArray, sizeof(float), numStressVectorComponent, stressTensorIn )==0 ) {
+ if (feof(stressTensorIn)) {
+ fprintf(stderr, "Error (reached EOF prematurely) while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
+ exit(1);
+ } else {
+ fprintf(stderr, "Error while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
+ exit(1);
+ }
+ }
+ /*
+ * Report error and bail if we pick up NaNs in any of the stress components.
+ */
+ if(isnan(stressTensorArray[0]) || isnan(stressTensorArray[1])
+ || isnan(stressTensorArray[2]) || isnan(stressTensorArray[3])
+ || isnan(stressTensorArray[4]) || isnan(stressTensorArray[5])) {
+ fprintf(stderr,"NaN in stress tensor file\n");
+ abort();
+ }
+ /*
+ * Build average stress tensor for element by summing tetrahedral tensor components
+ * - even though it's symmetric, do for all 9 components in case we pick the wrong ones before diagonalization
+ */
+ elementStressTensor[0][0]+=stressTensorArray[0]/(double)Tetrahedra_Count;
+ elementStressTensor[1][1]+=stressTensorArray[1]/(double)Tetrahedra_Count;
+ elementStressTensor[2][2]+=stressTensorArray[2]/(double)Tetrahedra_Count;
+ elementStressTensor[0][1]+=stressTensorArray[3]/(double)Tetrahedra_Count;
+ elementStressTensor[0][2]+=stressTensorArray[4]/(double)Tetrahedra_Count;
+ elementStressTensor[1][2]+=stressTensorArray[5]/(double)Tetrahedra_Count;
+ elementStressTensor[1][0]+=stressTensorArray[3]/(double)Tetrahedra_Count;
+ elementStressTensor[2][0]+=stressTensorArray[4]/(double)Tetrahedra_Count;
+ elementStressTensor[2][1]+=stressTensorArray[5]/(double)Tetrahedra_Count;
}
/*
* Diagonalize and find principal stresses from mean stress tensor for element
More information about the CIG-COMMITS
mailing list