[cig-commits] r20540 - long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI
echoi at geodynamics.org
echoi at geodynamics.org
Tue Jul 24 14:47:35 PDT 2012
Author: echoi
Date: 2012-07-24 14:47:34 -0700 (Tue, 24 Jul 2012)
New Revision: 20540
Modified:
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h
Log:
Made compatible with the new two-step mapping of element variables.
Cleaned up legacies from the SPR version.
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -35,7 +35,6 @@
#include <StGermain/StGermain.h>
#include <StGermain/FD/FD.h>
#include "Snac/Snac.h"
-#include "units.h"
#include "types.h"
#include "Element.h"
#include "Constitutive.h"
@@ -80,16 +79,16 @@
/* If this is a ViscoPlastic material, calculate its stress. */
if ( material->rheology & Snac_Material_ViscoPlastic ) {
Tetrahedra_Index tetra_I;
- Node_LocalIndex node_lI;
+ Node_LocalIndex node_lI;
/* viscoplastic material properties */
const double bulkm = material->lambda + 2.0f * material->mu/3.0f;
StressTensor* stress;
StrainTensor* strain;
- Strain plasticStrain;
- ViscoPlastic* viscosity;
+ Strain plasticStrain;
+ double* viscosity;
double straind0,straind1,straind2,stressd0,stressd1,stressd2;
- double Stress0[3][3];
+ double Stress0[3][3];
double trace_strain;
double trace_stress;
double temp;
@@ -97,10 +96,10 @@
double vic2;
double VolumicStress;
double rviscosity=material->refvisc;
- double rmu= material->mu;
+ double rmu= material->mu;
double srJ2;
double avgTemp;
- plModel yieldcriterion=material->yieldcriterion;
+ plModel yieldcriterion=material->yieldcriterion;
/* For now reference values of viscosity, second invariant of deviatoric */
/* strain rate and reference temperature are being hard wired ( these specific */
@@ -120,13 +119,15 @@
double tension_cutoff=0.0;
const double degrad = PI / 180.0f;
double totalVolume=0.0f,depls=0.0f;
- int principal_stresses(StressTensor* stress,double sp[],double cn[3][3]);
unsigned int i;
double tmp=0.0;
const double a1 = material->lambda + 2.0f * material->mu ;
const double a2 = material->lambda ;
int ind=0;
+ int principal_stresses(StressTensor* stress,double sp[],double cn[3][3]);
+ int principal_stresses_orig(StressTensor* stress, double sp[3], double cn[3][3]);
+
/* printf("Entered ViscoPlastic update \n"); */
/* Work out the plastic material properties of this element */
@@ -360,7 +361,8 @@
for( m = 0; m < 3; m++ ) {
for( n = m; n < 3; n++ ) {
for( k = 0; k < 3; k++ ) {
- (*stress)[m][n] += cn[k][m] * cn[k][n] * s[k];
+ /* (*stress)[m][n] += cn[k][m] * cn[k][n] * s[k]; */
+ (*stress)[m][n] += cn[m][k] * cn[n][k] * s[k];
}
}
}
@@ -382,8 +384,279 @@
}
}
+int principal_stresses(StressTensor* stress, double d[3], double V[3][3])
+{
-int principal_stresses(StressTensor* stress, double sp[3], double cn[3][3])
+ int i,j;
+ double e[3];
+ void tred2(int n, double d[3], double e[3], double V[3][3]);
+ void tql2 (int n, double d[3], double e[3], double V[3][3]);
+
+ for (i = 0; i < 3; i++) {
+ for (j = 0; j < 3; j++) {
+ V[i][j] = (*stress)[i][j];
+ }
+ }
+
+ // Tridiagonalize.
+ tred2(3, d, e, V);
+
+ // Diagonalize.
+ tql2(3, d, e, V);
+
+ /* fprintf(stderr,"eig 1: %e eigV: %e %e %e\n",d[0],V[0][0],V[1][0],V[2][0]); */
+ /* fprintf(stderr,"eig 2: %e eigV: %e %e %e\n",d[1],V[0][1],V[1][1],V[2][1]); */
+ /* fprintf(stderr,"eig 3: %e eigV: %e %e %e\n",d[2],V[0][2],V[1][2],V[2][2]); */
+ return(1);
+}
+
+
+// Symmetric Householder reduction to tridiagonal form.
+
+void tred2(int n, double d[3], double e[3], double V[3][3]) {
+
+ int i,j,k;
+ // This is derived from the Algol procedures tred2 by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ for (j = 0; j < n; j++) {
+ d[j] = V[n-1][j];
+ }
+
+ // Householder reduction to tridiagonal form.
+
+ for (i = n-1; i > 0; i--) {
+
+ // Scale to avoid under/overflow.
+
+ double scale = 0.0;
+ double h = 0.0;
+ for (k = 0; k < i; k++) {
+ scale = scale + abs(d[k]);
+ }
+ if (scale == 0.0) {
+ e[i] = d[i-1];
+ for (j = 0; j < i; j++) {
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ V[j][i] = 0.0;
+ }
+ } else {
+
+ double f, g, hh;
+ // Generate Householder vector.
+ for (k = 0; k < i; k++) {
+ d[k] /= scale;
+ h += d[k] * d[k];
+ }
+ f = d[i-1];
+ g = sqrt(h);
+ if (f > 0) {
+ g = -g;
+ }
+ e[i] = scale * g;
+ h = h - f * g;
+ d[i-1] = f - g;
+ for (j = 0; j < i; j++) {
+ e[j] = 0.0;
+ }
+
+ // Apply similarity transformation to remaining columns.
+
+ for (j = 0; j < i; j++) {
+ f = d[j];
+ V[j][i] = f;
+ g = e[j] + V[j][j] * f;
+ for (k = j+1; k <= i-1; k++) {
+ g += V[k][j] * d[k];
+ e[k] += V[k][j] * f;
+ }
+ e[j] = g;
+ }
+ f = 0.0;
+ for (j = 0; j < i; j++) {
+ e[j] /= h;
+ f += e[j] * d[j];
+ }
+ hh = f / (h + h);
+ for (j = 0; j < i; j++) {
+ e[j] -= hh * d[j];
+ }
+ for (j = 0; j < i; j++) {
+ f = d[j];
+ g = e[j];
+ for (k = j; k <= i-1; k++) {
+ V[k][j] -= (f * e[k] + g * d[k]);
+ }
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ }
+ }
+ d[i] = h;
+ }
+
+ // Accumulate transformations.
+
+ for (i = 0; i < n-1; i++) {
+ double h, g;
+ V[n-1][i] = V[i][i];
+ V[i][i] = 1.0;
+ h = d[i+1];
+ if (h != 0.0) {
+ for (k = 0; k <= i; k++) {
+ d[k] = V[k][i+1] / h;
+ }
+ for (j = 0; j <= i; j++) {
+ g = 0.0;
+ for (k = 0; k <= i; k++) {
+ g += V[k][i+1] * V[k][j];
+ }
+ for (k = 0; k <= i; k++) {
+ V[k][j] -= g * d[k];
+ }
+ }
+ }
+ for (k = 0; k <= i; k++) {
+ V[k][i+1] = 0.0;
+ }
+ }
+ for (j = 0; j < n; j++) {
+ d[j] = V[n-1][j];
+ V[n-1][j] = 0.0;
+ }
+ V[n-1][n-1] = 1.0;
+ e[0] = 0.0;
+}
+
+// Symmetric tridiagonal QL algorithm.
+
+void tql2 (int n, double d[3], double e[3], double V[3][3]) {
+
+ // This is derived from the Algol procedures tql2, by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+ int i,j,k,l;
+ double f = 0.0;
+ double tst1 = 0.0;
+ double eps = pow(2.0,-52.0);
+
+ for (i = 1; i < n; i++) {
+ e[i-1] = e[i];
+ }
+ e[n-1] = 0.0;
+
+ for (l = 0; l < n; l++) {
+
+ int m = l;
+
+ // Find small subdiagonal element
+ tst1 = MAX(tst1,abs(d[l]) + abs(e[l]));
+
+ // Original while-loop from Java code
+ while (m < n) {
+ if (abs(e[m]) <= eps*tst1) {
+ break;
+ }
+ m++;
+ }
+
+
+ // If m == l, d[l] is an eigenvalue,
+ // otherwise, iterate.
+
+ if (m > l) {
+ int iter = 0;
+ do {
+ // Compute implicit shift
+
+ double g = d[l];
+ double p = (d[l+1] - g) / (2.0 * e[l]);
+ double r = hypot(p,1.0);
+ double dl1, h, c, c2, c3, el1, s, s2;
+
+ iter = iter + 1; // (Could check iteration count here.)
+
+ if (p < 0) {
+ r = -r;
+ }
+ d[l] = e[l] / (p + r);
+ d[l+1] = e[l] * (p + r);
+ dl1 = d[l+1];
+ h = g - d[l];
+ for (i = l+2; i < n; i++) {
+ d[i] -= h;
+ }
+ f = f + h;
+
+ // Implicit QL transformation.
+
+ p = d[m];
+ c = 1.0;
+ c2 = c;
+ c3 = c;
+ el1 = e[l+1];
+ s = 0.0;
+ s2 = 0.0;
+ for (i = m-1; i >= l; i--) {
+ c3 = c2;
+ c2 = c;
+ s2 = s;
+ g = c * e[i];
+ h = c * p;
+ r = hypot(p,e[i]);
+ e[i+1] = s * r;
+ s = e[i] / r;
+ c = p / r;
+ p = c * d[i] - s * g;
+ d[i+1] = h + s * (c * g + s * d[i]);
+
+ // Accumulate transformation.
+
+ for (k = 0; k < n; k++) {
+ h = V[k][i+1];
+ V[k][i+1] = s * V[k][i] + c * h;
+ V[k][i] = c * V[k][i] - s * h;
+ }
+ }
+ p = -s * s2 * c3 * el1 * e[l] / dl1;
+ e[l] = s * p;
+ d[l] = c * p;
+
+ // Check for convergence.
+
+ } while (abs(e[l]) > eps*tst1);
+ }
+ d[l] = d[l] + f;
+ e[l] = 0.0;
+ }
+
+ // Sort eigenvalues and corresponding vectors.
+
+ for (i = 0; i < n-1; i++) {
+ double p = d[i];
+ k = i;
+ for (j = i+1; j < n; j++) {
+ if (d[j] < p) {
+ k = j;
+ p = d[j];
+ }
+ }
+ if (k != i) {
+ d[k] = d[i];
+ d[i] = p;
+ for (j = 0; j < n; j++) {
+ p = V[j][i];
+ V[j][i] = V[j][k];
+ V[j][k] = p;
+ }
+ }
+ }
+}
+
+int principal_stresses_orig(StressTensor* stress, double sp[3], double cn[3][3])
{
double **a,**v,*d;
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -52,38 +52,11 @@
&tmpElement,
SnacViscoPlastic_ElementHandle );
char tmpBuf[PATH_MAX];
- /* Because plastic strain is not an array by itself, we must the "complex" constructor for Variable... the info needs to be
- * wrapped this generic way... */
- Index viscoplasticOffsetCount = 1;
- SizeT viscoplasticOffsets[] = {
- (SizeT)((char*)&tmpElementExt->aps - (char*)&tmpElement) };
- Variable_DataType viscoplasticDataTypes[] = { Variable_DataType_Double };
- Index viscoplasticDataTypeCounts[] = { 1 };
- // #if DEBUG
+#if DEBUG
if( context->rank == 0 ) printf( "In %s()\n", __func__ );
- // #endif
- /* Create the StGermain variable aps, which is stored on an element extension */
- Variable_New(
- "plStrain",
- viscoplasticOffsetCount,
- viscoplasticOffsets,
- viscoplasticDataTypes,
- viscoplasticDataTypeCounts,
- 0,
- &ExtensionManager_GetFinalSize( context->mesh->elementExtensionMgr ),
- &context->mesh->layout->decomp->elementDomainCount,
- (void**)&context->mesh->element,
- context->variable_Register );
- int element_lI;
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( context, element_lI );
- SnacViscoPlastic_Element* viscoplasticElement = ExtensionManager_Get( context->mesh->elementExtensionMgr, element,
- SnacViscoPlastic_ElementHandle );
- viscoplasticElement->aps = 0.0f;
- }//for
+#endif
-
/* Prepare the dump file */
sprintf( tmpBuf, "%s/plStrain.%u", context->outputPath, context->rank );
if( (contextExt->plStrainOut = fopen( tmpBuf, "w+" )) == NULL ) {
@@ -95,11 +68,6 @@
assert( contextExt->plStrainCheckpoint /* failed to open file for writing */ );
abort();
}
- sprintf( tmpBuf, "%s/avgPlStrainCP.%u", context->outputPath, context->rank );
- if( (contextExt->avgPlStrainCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
- assert( contextExt->avgPlStrainCheckpoint /* failed to open file for writing */ );
- abort();
- }
sprintf( tmpBuf, "%s/viscosity.%u", context->outputPath, context->rank );
if( (contextExt->viscOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->viscOut /* failed to open file for writing */ );
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h 2012-07-24 21:47:34 UTC (rev 20540)
@@ -46,9 +46,7 @@
FILE* plStrainOut;
FILE* viscOut;
FILE* plStrainCheckpoint;
- FILE* avgPlStrainCheckpoint;
FILE* viscCheckpoint;
- FILE* plstrainTensorOut;
};
/* Print the contents of the context extension */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -51,8 +51,6 @@
fclose( contextExt->plStrainOut );
if( contextExt->plStrainCheckpoint )
fclose( contextExt->plStrainCheckpoint );
- if( contextExt->avgPlStrainCheckpoint )
- fclose( contextExt->avgPlStrainCheckpoint );
if( contextExt->viscOut )
fclose( contextExt->viscOut );
if( contextExt->viscCheckpoint )
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -46,51 +46,50 @@
#endif
void SnacViscoPlastic_InitialConditions( void* _context, void* data ) {
- Snac_Context* context = (Snac_Context*)_context;
- Element_LocalIndex element_lI;
- double mu,vis_min,dt_maxwellI;
- double dt_maxwell = 0;
- const double mfrac = 9.0e-01;
+ Snac_Context* context = (Snac_Context*)_context;
+ Element_LocalIndex element_lI;
+ double mu,vis_min,dt_maxwellI;
+ double dt_maxwell = 0;
+ const double mfrac = 9.0e-01;
// loop over the phase using the dictionary
Dictionary_Entry_Value* materialList = Dictionary_Get( context->dictionary, "materials" );
- int PhaseI = 0;
+ int PhaseI = 0;
if( materialList ) {
Dictionary_Entry_Value* materialEntry = Dictionary_Entry_Value_GetFirstElement( materialList );
/* loop around the phases to initialize rheology */
while( materialEntry ) {
context->materialProperty[PhaseI].rheology |= Snac_Material_ViscoPlastic;
- mu = context->materialProperty[PhaseI].mu;
- vis_min = context->materialProperty[PhaseI].vis_min;
- dt_maxwellI = mfrac*vis_min/mu;
- dt_maxwell = (PhaseI==0)?dt_maxwellI:min(dt_maxwellI,dt_maxwell);
+ mu = context->materialProperty[PhaseI].mu;
+ vis_min = context->materialProperty[PhaseI].vis_min;
+ dt_maxwellI = mfrac*vis_min/mu;
+ dt_maxwell = (PhaseI==0)?dt_maxwellI:min(dt_maxwellI,dt_maxwell);
PhaseI++;
- materialEntry = materialEntry->next;
+ materialEntry = materialEntry->next;
}
}
else {
context->materialProperty[PhaseI].rheology |= Snac_Material_ViscoPlastic;
- mu = context->materialProperty[PhaseI].mu;
- vis_min = context->materialProperty[PhaseI].vis_min;
- dt_maxwell = mfrac*vis_min/mu;
+ mu = context->materialProperty[PhaseI].mu;
+ vis_min = context->materialProperty[PhaseI].vis_min;
+ dt_maxwell = mfrac*vis_min/mu;
}
+
if( context->dt > dt_maxwell) {
if(context->dtType == Snac_DtType_Constant) {
fprintf(stderr,"dt(%e) should be smaller than dt_maxwell(%e) (mu=%e vis_min=%e mfrac=%e)\n",context->dt,dt_maxwell,mu,vis_min,mfrac);
}
else context->dt = dt_maxwell;
}
+
if( context->restartTimestep > 0 ) {
FILE* plStrainIn;
- FILE* avgPlStrainIn;
char path[PATH_MAX];
sprintf(path, "%s/snac.plStrain.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
Journal_Firewall( (plStrainIn = fopen(path,"r")) != NULL, "Can't find %s", path );
- sprintf(path, "%s/snac.avgPlStrain.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
- Journal_Firewall( (avgPlStrainIn = fopen(path,"r")) != NULL, "Can't find %s", path );
/* read in restart file to reconstruct the previous plastic strain.*/
for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
@@ -103,7 +102,7 @@
/* To-Do: Read in viscosity and assign it to tets and elements. */
if( material->yieldcriterion == mohrcoulomb ) {
Tetrahedra_Index tetra_I;
- double avgPlStrain = 0.0f;
+ double depls = 0.0f;
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
double tetraPlStrain = 0.0f;
@@ -113,16 +112,13 @@
fscanf( plStrainIn, "%le", &tetraPlStrain );
viscoplasticElement->plasticStrain[tetra_I] = tetraPlStrain;
+ depls += viscoplasticElement->plasticStrain[tetra_I];
}/* for tets */
- /* volume-averaged accumulated plastic strain, aps */
- fscanf( avgPlStrainIn, "%le", &avgPlStrain );
- viscoplasticElement->aps = avgPlStrain;
+ viscoplasticElement->aps = depls/Tetrahedra_Count;
}/* if(mohrcoulomb) */
}/* for elements */
if( plStrainIn )
fclose( plStrainIn );
- if( avgPlStrainIn )
- fclose( avgPlStrainIn );
}/* if restarting.*/
else {
/* Set the plastic element initial conditions */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def 2012-07-24 21:47:34 UTC (rev 20540)
@@ -34,9 +34,7 @@
def_inc = Snac/ViscoPlastic
def_srcs = \
- Node.c \
Element.c \
- Mesh.c \
Context.c \
ConstructExtensions.c \
Constitutive.c \
@@ -48,10 +46,7 @@
def_hdrs = \
types.h \
- units.h \
- Node.h \
Element.h \
- Mesh.h \
Context.h \
ConstructExtensions.h \
Constitutive.h \
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -105,11 +105,7 @@
fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainCheckpoint );
}
fflush( contextExt->plStrainCheckpoint );
-
-
- fwrite( &avgPlasticStrain, sizeof(float), 1, contextExt->avgPlStrainCheckpoint );
}
- fflush( contextExt->avgPlStrainCheckpoint );
}
@@ -188,54 +184,3 @@
}
fflush( contextExt->viscCheckpoint );
}
-
-
-#if 0
-void _SnacViscoPlastic_DumpPlasticStrainTensor( void* _context ) {
- Snac_Context* context = (Snac_Context*) _context;
- SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacViscoPlastic_ContextHandle );
-
-#if DEBUG
- printf( "In %s()\n", __func__ );
-#endif
-
-
- if( context->timeStep ==0 || (context->timeStep-1) % context->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
-
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( context, element_lI );
- SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
- context->mesh->elementExtensionMgr,
- element,
- SnacViscoPlastic_ElementHandle );
- const Snac_Material* material = &context->materialProperty[element->material_I];
- Tetrahedra_Index tetra_I;
-
-
- if( material->yieldcriterion == druckerprager ) {
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- Index i,j;
- float tensor[3][3];
- for(i=0;i<3;i++)
- for(j=0;j<3;j++) {
- tensor[i][j] = elementExt->plasticstrainTensor[tetra_I][i][j];
- fwrite( &tensor[i][j], sizeof(float), 1, contextExt->plstrainTensorOut );
- }
- }
- }
- else if( material->yieldcriterion == mohrcoulomb ) {
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- float tetraStrain;
- tetraStrain = elementExt->plasticStrain[tetra_I];
- fwrite( &tetraStrain, sizeof(float), 1, contextExt->plstrainTensorOut );
- }
- }
- }
- fflush( contextExt->plstrainTensorOut );
- }
-}
-#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.h 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.h 2012-07-24 21:47:34 UTC (rev 20540)
@@ -38,6 +38,5 @@
void _SnacViscoPlastic_CheckpointPlasticStrain( void* _context );
void _SnacViscoPlastic_DumpViscosity( void* _context );
void _SnacViscoPlastic_CheckpointViscosity( void* _context );
- void _SnacViscoPlastic_DumpPlasticStrainTensor( void* _context );
#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -33,26 +33,23 @@
#include <StGermain/StGermain.h>
#include <StGermain/FD/FD.h>
#include "Snac/Snac.h"
+#include "Snac/Remesher/Remesher.h"
+#include "types.h"
#include "ConstructExtensions.h"
#include "Constitutive.h"
#include "Context.h"
#include "Destroy.h"
#include "InitialConditions.h"
#include "Element.h"
-#include "Mesh.h"
-#include "Node.h"
#include "Output.h"
#include "Register.h"
#include "Remesh.h"
-#include "types.h"
#include <stdio.h>
/* Textual name of this class */
const Type SnacViscoPlastic_Type = "SnacViscoPlastic";
-ExtensionInfo_Index SnacViscoPlastic_NodeHandle;
ExtensionInfo_Index SnacViscoPlastic_ElementHandle;
-ExtensionInfo_Index SnacViscoPlastic_MeshHandle;
ExtensionInfo_Index SnacViscoPlastic_ContextHandle;
Index _SnacViscoPlastic_Register( PluginsManager* pluginsMgr ) {
@@ -82,6 +79,7 @@
void _SnacViscoPlastic_Construct( void* component, Stg_ComponentFactory* cf, void* data ) {
Snac_Context* context;
EntryPoint* interpolateElementEP;
+ EntryPoint* copyElementEP;
/* Retrieve context. */
context = (Snac_Context*)Stg_ComponentFactory_ConstructByName( cf, "context", Snac_Context, True, data );
@@ -90,19 +88,9 @@
printf( "In: _SnacViscoPlastic_Register( void*, void* )\n" );
#endif
- /* Add extensions to nodes, elements and the context */
- SnacViscoPlastic_NodeHandle = ExtensionManager_Add( context->mesh->nodeExtensionMgr, SnacViscoPlastic_Type, sizeof(SnacViscoPlastic_Node) );
SnacViscoPlastic_ElementHandle = ExtensionManager_Add( context->mesh->elementExtensionMgr, SnacViscoPlastic_Type, sizeof(SnacViscoPlastic_Element) );
- SnacViscoPlastic_MeshHandle = ExtensionManager_Add( context->meshExtensionMgr, SnacViscoPlastic_Type, sizeof(SnacViscoPlastic_Mesh) );
SnacViscoPlastic_ContextHandle = ExtensionManager_Add( context->extensionMgr, SnacViscoPlastic_Type, sizeof(SnacViscoPlastic_Context) );
- #ifdef DEBUG
- printf( "\tcontext extension handle: %u\n", SnacViscoPlastic_ContextHandle );
- printf( "\tmesh extension handle: %u\n", SnacViscoPlastic_MeshHandle );
- printf( "\tnode extension handle: %u\n", SnacViscoPlastic_ElementHandle );
- printf( "\telement extension handle: %u\n", SnacViscoPlastic_NodeHandle );
- #endif
-
/* Add extensions to the entry points */
EntryPoint_Append(
Context_GetEntryPoint( context, Snac_EP_Constitutive ),
@@ -151,6 +139,16 @@
SnacViscoPlastic_Type );
}
+ /* Add extensions to the copyElement entry point, but it will only exist if the remesher is loaded. */
+ copyElementEP = Context_GetEntryPoint( context, "SnacRemesher_EP_CopyElement" );
+ if( copyElementEP ) {
+ EntryPoint_Append(
+ copyElementEP,
+ SnacViscoPlastic_Type,
+ _SnacViscoPlastic_CopyElement,
+ SnacViscoPlastic_Type );
+ }
+
/* Construct. */
_SnacViscoPlastic_ConstructExtensions( context, data );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h 2012-07-24 21:47:34 UTC (rev 20540)
@@ -44,12 +44,8 @@
/* Textual name of this class */
extern const Type SnacViscoPlastic_Type;
- extern ExtensionInfo_Index SnacViscoPlastic_NodeHandle;
-
extern ExtensionInfo_Index SnacViscoPlastic_ElementHandle;
- extern ExtensionInfo_Index SnacViscoPlastic_MeshHandle;
-
extern ExtensionInfo_Index SnacViscoPlastic_ContextHandle;
Index _SnacViscoPlastic_Register( PluginsManager* pluginsMgr );
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c 2012-07-24 21:47:34 UTC (rev 20540)
@@ -45,28 +45,23 @@
void _SnacViscoPlastic_InterpolateElement( void* _context,
Element_LocalIndex dstEltInd,
Tetrahedra_Index dstTetInd,
- Snac_Element* dstElements,
+ SnacRemesher_Element* dstElements,
Element_DomainIndex srcEltInd,
Tetrahedra_Index srcTetInd )
{
- Snac_Context* context = (Snac_Context*)_context;
- Mesh* mesh = context->mesh;
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
+ Snac_Context* context = (Snac_Context*)_context;
+ Mesh* mesh = context->mesh;
+ SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
mesh,
SnacRemesher_MeshHandle );
- HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
- Snac_Element* element = (Snac_Element*)ExtensionManager_At( context->mesh->elementExtensionMgr,
- dstElements,
- dstEltInd );
- SnacViscoPlastic_Element* elementExt = ExtensionManager_Get( context->mesh->elementExtensionMgr,
- element,
- SnacViscoPlastic_ElementHandle );
- Element_DomainIndex eltdI[8],eldI,eldJ,eldK;
- Index coef_I;
- Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
- Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
- Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
- enum { threeD, xy, undefined } geomType;
+ HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
+ SnacRemesher_Element* dstElt = &dstElements[dstEltInd];
+ Element_DomainIndex eltdI[8],eldI,eldJ,eldK;
+ Index coef_I;
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
+ enum { threeD, xy, undefined } geomType;
#ifdef DEBUG
printf( "element_lI: %u, fromElement_lI: %u\n", dstElementInd, srcElementInd );
@@ -95,7 +90,7 @@
eltdI[6] = (eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ;
eltdI[7] = eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ;
- elementExt->plasticStrain[dstTetInd] = 0.0;
+ dstElt->plasticStrain[dstTetInd] = 0.0;
for(coef_I=0;coef_I<4;coef_I++) {
/* The actual src elements are the four apexes of a tet (srcTetInd) in the old barycenter grid. */
Snac_Element* srcElt = Snac_Element_At( context,
@@ -106,7 +101,7 @@
SnacViscoPlastic_ElementHandle );
/* Weights are associated only with destination element but not on the tet level.
So, "dstTetInd" is used in both source and destination terms. */
- elementExt->plasticStrain[dstTetInd] += meshExt->barcord[dstEltInd].L[coef_I]*srcEltExt->plasticStrain[dstTetInd];
+ dstElt->plasticStrain[dstTetInd] += meshExt->barcord[dstEltInd].L[coef_I]*srcEltExt->plasticStrain[dstTetInd];
}
break;
case xy:
@@ -120,7 +115,7 @@
eltdI[2] = (eldI+1) + (eldJ+1)*neldI;
eltdI[3] = eldI + (eldJ+1)*neldI;
- elementExt->plasticStrain[dstTetInd] = 0.0;
+ dstElt->plasticStrain[dstTetInd] = 0.0;
for(coef_I=0;coef_I<3;coef_I++) {
/* The actual src elements are the four apexes of a tet (srcTetInd) in the old barycenter grid. */
Snac_Element* srcElt = Snac_Element_At( context,
@@ -131,8 +126,30 @@
SnacViscoPlastic_ElementHandle );
/* Weights are associated only with destination element but not on the tet level.
So, "dstTetInd" is used in both source and destination terms. */
- elementExt->plasticStrain[dstTetInd] += meshExt->barcord[dstEltInd].L[coef_I]*srcEltExt->plasticStrain[dstTetInd];
+ dstElt->plasticStrain[dstTetInd] += meshExt->barcord[dstEltInd].L[coef_I]*srcEltExt->plasticStrain[dstTetInd];
}
break;
}
}
+
+/*
+ Copy interpolated values stored in newElement array
+ back to the original element array.
+*/
+void _SnacViscoPlastic_CopyElement( void* _context,
+ Element_LocalIndex eltInd,
+ Tetrahedra_Index tetInd,
+ SnacRemesher_Element* srcEltArray )
+{
+
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacRemesher_Element* srcElt = &srcEltArray[eltInd];
+ Snac_Element* dstElt = Snac_Element_At( context, eltInd );
+ SnacViscoPlastic_Element* dstEltExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ dstElt,
+ SnacViscoPlastic_ElementHandle );
+
+ dstEltExt->plasticStrain[tetInd] = srcElt->plasticStrain[tetInd];
+
+}
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.h 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.h 2012-07-24 21:47:34 UTC (rev 20540)
@@ -43,8 +43,13 @@
void _SnacViscoPlastic_InterpolateElement( void* _context,
Element_LocalIndex dstEltInd,
Tetrahedra_Index dstTetInd,
- Snac_Element* dstElements,
+ SnacRemesher_Element* dstElements,
Element_DomainIndex srcEltInd,
Tetrahedra_Index srcTetInd );
-
+
+void _SnacViscoPlastic_CopyElement( void* _context,
+ Element_LocalIndex EltInd,
+ Tetrahedra_Index tetInd,
+ SnacRemesher_Element* srcEltArray );
+
#endif /* __SnacViscoPlastic_Remesh_h__ */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h 2012-07-24 21:46:28 UTC (rev 20539)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h 2012-07-24 21:47:34 UTC (rev 20540)
@@ -48,12 +48,9 @@
#include "ConstructExtensions.h"
#include "Element.h"
#include "InitialConditions.h"
-#include "Mesh.h"
-#include "Node.h"
#include "Output.h"
#include "Register.h"
-#include "Remesh.h"
+/* #include "Remesh.h" */
#include "types.h"
-#include "units.h"
#endif /* __SnacViscoPlastic_h__ */
More information about the CIG-COMMITS
mailing list