[cig-commits] r16473 - in long/3D/SNAC/trunk/Snac/plugins: . remesherSPR viscoplasticSPR
echoi at geodynamics.org
echoi at geodynamics.org
Wed Mar 31 08:24:04 PDT 2010
Author: echoi
Date: 2010-03-31 08:24:03 -0700 (Wed, 31 Mar 2010)
New Revision: 16473
Added:
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Make.mm
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Makefile.def
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ViscoPlastic.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/makefile
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/types.h
long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/units.h
Modified:
long/3D/SNAC/trunk/Snac/plugins/remesherSPR/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshCoords.c
long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshElements.c
long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshNodes.c
Log:
remesherSPR plugin now works.
viscoplasticSPR has been newly added, which is compatible with SPR-based remesehr.
Modified: long/3D/SNAC/trunk/Snac/plugins/remesherSPR/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesherSPR/Remesh.c 2010-03-31 13:44:34 UTC (rev 16472)
+++ long/3D/SNAC/trunk/Snac/plugins/remesherSPR/Remesh.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -49,14 +49,14 @@
void _SnacRemesher_Remesh( void* _context, void* data ) {
Snac_Context* context = (Snac_Context*)_context;
SnacRemesher_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
+ context->extensionMgr,
+ context,
+ SnacRemesher_ContextHandle );
Mesh* mesh = context->mesh;
SnacRemesher_Mesh* meshExt = ExtensionManager_Get(
- context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
+ context->meshExtensionMgr,
+ mesh,
+ SnacRemesher_MeshHandle );
Bool remesh;
Journal_DPrintf( context->debug, "In: %s\n", __func__ );
@@ -123,7 +123,7 @@
/* Remesh the coordinates. */
_SnacRemesher_NewCoords( context );
-
+
/* Interpolate current nodal values onto new coordinates. */
meshExt->newNodes = (Snac_Node*)ExtensionManager_Malloc( mesh->nodeExtensionMgr, mesh->nodeLocalCount );
_SnacRemesher_InterpolateNodes( context );
@@ -132,28 +132,36 @@
This simple copy works because bottoms nodes are always bottom
and remeshing doesn't change the node number. */
for( newNode_i = 0; newNode_i < mesh->nodeLocalCount; newNode_i++ ) {
- Snac_Node* dstNode =
- (Snac_Node*)ExtensionManager_At( context->mesh->nodeExtensionMgr,
- meshExt->newNodes,
- newNode_i );
+ Snac_Node* dstNode = (Snac_Node*)ExtensionManager_At( context->mesh->nodeExtensionMgr,
+ meshExt->newNodes,
+ newNode_i );
- Snac_Node* srcNode =
- Snac_Node_At( context, newNode_i );
+ Snac_Node* srcNode = Snac_Node_At( context, newNode_i );
dstNode->residualFr = srcNode->residualFr;
dstNode->residualFt = srcNode->residualFt;
}
- /* Copy accross the new coord, node & element information to the current arrays. */
- memcpy( mesh->nodeCoord, meshExt->newNodeCoords, mesh->nodeLocalCount * sizeof(Coord) );
- memcpy( mesh->node, meshExt->newNodes, mesh->nodeExtensionMgr->finalSize * mesh->nodeLocalCount );
+ /* Copy accross the new coord */
+ memcpy( mesh->nodeCoord, meshExt->newNodeCoords, mesh->nodeLocalCount * sizeof(Coord) );
+ /* Copy node structures to the current array */
+ /* The next single line can possibly replace the for loop.
+ memcpy( mesh->node, meshExt->newNodes, mesh->nodeExtensionMgr->finalSize * mesh->nodeLocalCount ); */
+ for( newNode_i = 0; newNode_i < mesh->nodeLocalCount; newNode_i++ ) {
+ Snac_Node* dstNode = Snac_Node_At( context, newNode_i );
+ Snac_Node* srcNode = (Snac_Node*)ExtensionManager_At( context->mesh->nodeExtensionMgr,
+ meshExt->newNodes,
+ newNode_i );
+ memcpy( dstNode, srcNode, mesh->nodeExtensionMgr->finalSize );
+ }
+
/* Simply average the recovered fields at the barycenter of each tet and run updateElements. */
- _SnacRemesher_InterpolateElements( context );
- _SnacRemesher_UpdateElements( context );
+ _SnacRemesher_InterpolateElements( context );
+ _SnacRemesher_UpdateElements( context );
/* Free some space, as it won't be needed until the next remesh. */
- ExtensionManager_Free( mesh->nodeExtensionMgr, meshExt->newNodes );
+ ExtensionManager_Free( mesh->nodeExtensionMgr, meshExt->newNodes );
meshExt->newNodes = NULL;
/*
Modified: long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshCoords.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshCoords.c 2010-03-31 13:44:34 UTC (rev 16472)
+++ long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshCoords.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -542,7 +542,6 @@
}
//fprintf(stderr,"dett=%e det0=%e det1=%e det2=%e shape=%e %e %e\n",dett,det[0],det[1],det[2],shape[0],shape[1],shape[2]);
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
meshExt->newNodeCoord[node_lI][0] = initialNodeCoord[0];
meshExt->newNodeCoord[node_lI][1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
//meshExt->newNodeCoord[node_lI][1] -= TOL*fabs(meshExt->newNodeCoord[node_lI][1]);
@@ -726,7 +725,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
meshExt->newNodeCoord[node_lI][0] = initialNodeCoord[0];
meshExt->newNodeCoord[node_lI][1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
meshExt->newNodeCoord[node_lI][2] = initialNodeCoord[2];
@@ -918,7 +916,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
meshExt->newNodeCoord[node_lI][0] = initialNodeCoord[0];
meshExt->newNodeCoord[node_lI][1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
meshExt->newNodeCoord[node_lI][2] = initialNodeCoord[2];
@@ -1106,7 +1103,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
meshExt->newNodeCoord[node_lI][0] = initialNodeCoord[0];
meshExt->newNodeCoord[node_lI][1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
meshExt->newNodeCoord[node_lI][2] = initialNodeCoord[2];
@@ -1380,7 +1376,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
x1[0] = initialNodeCoordS[0];
x1[1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
x1[2] = initialNodeCoordS[2];
@@ -1578,7 +1573,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
x1[0] = initialNodeCoordS[0];
x1[1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
x1[2] = initialNodeCoordS[2];
@@ -1782,7 +1776,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
x1[0] = initialNodeCoordS[0];
x1[1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
x1[2] = initialNodeCoordS[2];
@@ -1979,7 +1972,6 @@
shape[tNode_I] = det[tNode_I] / dett;
}
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
x1[0] = initialNodeCoordS[0];
x1[1] = y1*shape[0] + y2*shape[1] + y3*shape[2];
x1[2] = initialNodeCoordS[2];
Modified: long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshElements.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshElements.c 2010-03-31 13:44:34 UTC (rev 16472)
+++ long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshElements.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -84,8 +84,8 @@
/* if( element->tetra[tetra_I].material_I==2 ) */
/* fprintf(stderr,"\t int matId=%e (%d)\n",materialInd, (Material_Index)((materialInd/(100*Tetrahedra_Count))+0.5)); */
}
- if( element->material_I ==2 ) //element_lI==0 )
- fprintf(stderr,"int matId=%d\n",(Material_Index)(materialInd/10+0.5));
+/* if( element->material_I ==2 ) //element_lI==0 ) */
+/* fprintf(stderr,"int matId=%d\n",(Material_Index)(materialInd/10+0.5)); */
//element->material_I = ((Material_Index)(((double)materialInd/(Tetrahedra_Count))+0.5))/100;
element->material_I = (Material_Index)(materialInd/10+0.5);
@@ -113,21 +113,37 @@
Index i,j;
double strain[6];
double stress[6];
- double materialInd;
+ float materialInd;
+ double density;
/* Compute averages of recovered fields for this tet.*/
for(j=0;j<6;j++) {
strain[j] = 0.0;
stress[j] = 0.0;
for(i=0;i<4;i++) {
- strain[j] += (0.25f*Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][i])->strainSPR[j]);
- stress[j] += (0.25f*Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][i])->stressSPR[j]);
+ Index dstNodeNum = Element_Node_I(context, dstEltInd, TetraToNode[dstTetInd][i]);
+ Snac_Node* dstNode = Snac_Node_At( context, dstNodeNum );
+ strain[j] += 0.25f * dstNode->strainSPR[j];
+ stress[j] += 0.25f * dstNode->stressSPR[j];
+ /* The next two lines can possibly replace the above codes, but need testing. */
+ /* strain[j] += (0.25f*Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][i])->strainSPR[j]); */
+ /* stress[j] += (0.25f*Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][i])->stressSPR[j]); */
}
}
materialInd = 0.0;
- for(i=0;i<4;i++)
- materialInd += 0.25f*Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][i])->material_ISPR;
+ for(i=0;i<4;i++) {
+ Index dstNodeNum = Element_Node_I(context, dstEltInd, TetraToNode[dstTetInd][i]);
+ Snac_Node* dstNode = Snac_Node_At( context, dstNodeNum );
+ materialInd += 0.25f*(dstNode->material_ISPR);
+ }
+ density = 0.0;
+ for(i=0;i<4;i++) {
+ Index dstNodeNum = Element_Node_I(context, dstEltInd, TetraToNode[dstTetInd][i]);
+ Snac_Node* dstNode = Snac_Node_At( context, dstNodeNum );
+ density += 0.25f*(dstNode->densitySPR);
+ }
+
/* Assign the computed averages to the tets. */
dstElt->tetra[dstTetInd].strain[0][0] = strain[0];
dstElt->tetra[dstTetInd].strain[1][1] = strain[1];
@@ -149,7 +165,14 @@
dstElt->tetra[dstTetInd].stress[2][0] = stress[4];
dstElt->tetra[dstTetInd].stress[2][1] = stress[5];
- dstElt->tetra[dstTetInd].material_I = (Material_Index)(materialInd);
+ dstElt->tetra[dstTetInd].material_I = 0;
+/* if( dstElt->tetra[dstTetInd].material_I != 0 ) */
+/* fprintf(stderr,"%d %d %e %d\n",dstEltInd,dstTetInd,dstElt->tetra[dstTetInd].material_I, */
+/* (unsigned int)(dstElt->tetra[dstTetInd].material_I)); */
+
+ dstElt->tetra[dstTetInd].density = density;
+
+#if 0
if( dstElt->material_I == 2 )
fprintf(stderr,"el=%d tet=%d elmatI=%d matI=%e (%d) (%e %e %e %e)\n",dstEltInd,dstTetInd,dstElt->material_I,materialInd,
dstElt->tetra[dstTetInd].material_I,
@@ -158,4 +181,5 @@
Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][2])->material_ISPR,
Snac_Element_Node_P(context, dstEltInd, TetraToNode[dstTetInd][3])->material_ISPR
);
+#endif
}
Modified: long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshNodes.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshNodes.c 2010-03-31 13:44:34 UTC (rev 16472)
+++ long/3D/SNAC/trunk/Snac/plugins/remesherSPR/RemeshNodes.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -76,7 +76,7 @@
/*
** Populate arrays for recovered fields using the SPR method.
*/
- for( newNode_i = 0; newNode_i < mesh->nodeLocalCount; newNode_i++ )
+ for( newNode_i = 0; newNode_i < mesh->nodeLocalCount; newNode_i++ )
SnacRemesher_RecoverNode( context, contextExt, newNode_i );
/*
@@ -185,20 +185,20 @@
void interpolateNode( void* _context, Node_LocalIndex newNodeInd, Element_DomainIndex dEltInd ) {
Snac_Context* context = (Snac_Context*)_context;
SnacRemesher_Context* contextExt = ExtensionManager_Get( context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
- Mesh* mesh = context->mesh;
+ context,
+ SnacRemesher_ContextHandle );
+ Mesh* mesh = context->mesh;
SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
- NodeLayout* nLayout = mesh->layout->nodeLayout;
+ mesh,
+ SnacRemesher_MeshHandle );
+ NodeLayout* nLayout = mesh->layout->nodeLayout;
unsigned nEltNodes;
Node_DomainIndex* eltNodes;
- Coord newNodeCoord;
- Coord crds[8];
- double weights[4];
- unsigned tetNodeInds[4];
- unsigned eltNode_i;
+ Coord newNodeCoord;
+ Coord crds[8];
+ double weights[4];
+ unsigned tetNodeInds[4];
+ unsigned eltNode_i;
/* Extract the element's node indices. Note that there should always be eight of these. */
{
@@ -297,6 +297,8 @@
dstNode->material_ISPR = 0.0;
+ dstNode->densitySPR = 0.0;
+
/* Loop over each contributing node. */
for( tetNode_i = 0; tetNode_i < 4; tetNode_i++ ) {
Snac_Node* srcNode;
@@ -325,6 +327,8 @@
dstNode->stressSPR[5] += srcNode->stressSPR[5] * weights[tetNode_i];
dstNode->material_ISPR += srcNode->material_ISPR * weights[tetNode_i];
+
+ dstNode->densitySPR += srcNode->densitySPR * weights[tetNode_i];
}
}
@@ -334,7 +338,9 @@
{
Snac_Context* context = (Snac_Context*)_context;
Mesh* mesh = context->mesh;
- NodeLayout* nLayout = mesh->layout->nodeLayout;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+ NodeLayout* nLayout = layout->nodeLayout;
Snac_Node* node = Snac_Node_At( context, node_lI );
Coord* coord = Snac_NodeCoord_P( context, node_lI );
@@ -345,105 +351,163 @@
gsl_vector* vecaStrain[6];
gsl_vector* vecaStress[6];
gsl_vector* vecaMaterial_I;
+ gsl_vector* vecaDensity;
gsl_vector* vecbStrain[6];
gsl_vector* vecbStress[6];
gsl_vector* vecbMaterial_I;
+ gsl_vector* vecbDensity;
Index i,j;
-
+ IJK ijk;
+ Node_GlobalIndex node_gI = _MeshDecomp_Node_LocalToGlobal1D( decomp, node_lI );
+ Node_GlobalIndex gNodeI = decomp->nodeGlobal3DCounts[0];
+ Node_GlobalIndex gNodeJ = decomp->nodeGlobal3DCounts[1];
+ Node_GlobalIndex gNodeK = decomp->nodeGlobal3DCounts[2];
+ Node_GlobalIndex intNode_gI;
+ Node_DomainIndex intNode_lI;
+ unsigned int isBoundaryNode=0, patchCenterNum=1, patchCenterI, patchCenterList[2];
+
+ RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
+
+ /* If boundary node, find a (topologically?) nearest interior node.
+ Loosely following
+ Khoei and Gharehbaghi, 2007, The superconvergent patch recovery techinique and data transfer operators in 3D plasticity problems,
+ Finite Elements in Analysis and Design 43 (2007) 630-- 648. */
+ if( (gNodeI>2) && (ijk[0]==0) ) {
+ ijk[0] += 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeI>2) && (ijk[0]==decomp->nodeGlobal3DCounts[0]-1) ) {
+ ijk[0] -= 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeJ>2) && (ijk[1]==0) ) {
+ ijk[1] += 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeJ>2) && (ijk[1]==decomp->nodeGlobal3DCounts[1]-1) ) {
+ ijk[1] -= 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeK>2) && (ijk[2]==0) ) {
+ ijk[2] += 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeK>2) && (ijk[2]==decomp->nodeGlobal3DCounts[2]-1) ) {
+ ijk[2] -= 1;
+ isBoundaryNode = 1;
+ }
+
+ /* node_lI itself always becomes a patch center,
+ and if the current node is a boundary node, find an interior node and add it to the patch center list. */
+ patchCenterList[0] = node_lI;
+ if( isBoundaryNode ) {
+ patchCenterNum=2;
+ intNode_gI = ijk[0]+gNodeI*ijk[1]+gNodeI*gNodeJ*ijk[2];
+ patchCenterList[1] = Mesh_NodeMapGlobalToLocal( mesh, intNode_gI );
+ }
+
/* initialize gsl vectors and matrix. */
matA = gsl_matrix_alloc(4,4); gsl_matrix_set_zero( matA );
vecaMaterial_I = gsl_vector_alloc(4); gsl_vector_set_zero( vecaMaterial_I );
vecbMaterial_I = gsl_vector_alloc(4); gsl_vector_set_zero( vecbMaterial_I );
+ vecaDensity = gsl_vector_alloc(4); gsl_vector_set_zero( vecaDensity );
+ vecbDensity = gsl_vector_alloc(4); gsl_vector_set_zero( vecbDensity );
for(i=0;i<6;i++) {
vecaStrain[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecaStrain[i] );
vecaStress[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecaStress[i] );
vecbStrain[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecbStrain[i] );
vecbStress[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecbStress[i] );
}
-
- /* For each incident element, find inicident tets. */
- for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
- Element_DomainIndex element_dI = context->mesh->nodeElementTbl[node_lI][nodeElement_I];
+
+ /* For each patch center */
+ for( patchCenterI=0; patchCenterI < patchCenterNum; patchCenterI++ ) {
+ /* For each incident element, find inicident tets. */
+ for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
+ Element_DomainIndex element_dI = context->mesh->nodeElementTbl[patchCenterList[patchCenterI]][nodeElement_I];
- if( element_dI < mesh->elementDomainCount ) {
- Index elementTetra_I;
- Snac_Element* element = Snac_Element_At( context, element_dI );
+ if( element_dI < mesh->elementDomainCount ) {
+ Index elementTetra_I;
+ Snac_Element* element = Snac_Element_At( context, element_dI );
- /* Extract the element's node indices. Note that there should always be eight of these. */
- {
- Element_GlobalIndex element_gI;
+ /* Extract the element's node indices. Note that there should always be eight of these. */
+ {
+ Element_GlobalIndex element_gI;
- elements = Memory_Alloc_Array( Node_DomainIndex, nodeElementCount, "SnacRemesher" );
- element_gI = Mesh_ElementMapDomainToGlobal( mesh, element_dI );
- nLayout->buildElementNodes( nLayout, element_gI, elements );
- }
+ elements = Memory_Alloc_Array( Node_DomainIndex, nodeElementCount, "SnacRemesher" );
+ element_gI = Mesh_ElementMapDomainToGlobal( mesh, element_dI );
+ nLayout->buildElementNodes( nLayout, element_gI, elements );
+ }
- /* Convert global node indices to domain. */
- {
- unsigned eltNode_i;
+ /* Convert global node indices to domain. */
+ {
+ unsigned eltNode_i;
- for( eltNode_i = 0; eltNode_i < nodeElementCount; eltNode_i++ ) {
- elements[eltNode_i] = Mesh_NodeMapGlobalToDomain( mesh, elements[eltNode_i] );
+ for( eltNode_i = 0; eltNode_i < nodeElementCount; eltNode_i++ ) {
+ elements[eltNode_i] = Mesh_NodeMapGlobalToDomain( mesh, elements[eltNode_i] );
+ }
}
- }
- /* For each incident tetrahedron in the incident element,
- add up contributions to P, A, and b as in Zienkiewicz and Zhu (1992), p. 1336 */
- for( elementTetra_I = 0; elementTetra_I < Node_Element_Tetrahedra_Count;elementTetra_I++ ) {
- Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][elementTetra_I];
- Coord tCrds[4];
- double positionP[4] = {1.0,0.0,0.0,0.0};
- Index ii,jj;
+ /* For each incident tetrahedron in the incident element,
+ add up contributions to P, A, and b as in Zienkiewicz and Zhu (1992), p. 1336 */
+ for( elementTetra_I = 0; elementTetra_I < Node_Element_Tetrahedra_Count;elementTetra_I++ ) {
+ Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][elementTetra_I];
+ Coord tCrds[4];
+ double positionP[4] = {1.0,0.0,0.0,0.0};
+ Index ii,jj;
- /* Extract the tetrahedron's coordinates. */
- Vector_Set( tCrds[0], mesh->nodeCoord[elements[TetraToNode[tetra_I][0]]] );
- Vector_Set( tCrds[1], mesh->nodeCoord[elements[TetraToNode[tetra_I][1]]] );
- Vector_Set( tCrds[2], mesh->nodeCoord[elements[TetraToNode[tetra_I][2]]] );
- Vector_Set( tCrds[3], mesh->nodeCoord[elements[TetraToNode[tetra_I][3]]] );
+ /* Extract the tetrahedron's coordinates. */
+ Vector_Set( tCrds[0], mesh->nodeCoord[elements[TetraToNode[tetra_I][0]]] );
+ Vector_Set( tCrds[1], mesh->nodeCoord[elements[TetraToNode[tetra_I][1]]] );
+ Vector_Set( tCrds[2], mesh->nodeCoord[elements[TetraToNode[tetra_I][2]]] );
+ Vector_Set( tCrds[3], mesh->nodeCoord[elements[TetraToNode[tetra_I][3]]] );
- for(ii=1;ii<4;ii++)
- for(jj=0;jj<4;jj++)
- positionP[ii] += (0.25f * tCrds[jj][ii-1]);
+ for(ii=1;ii<4;ii++)
+ for(jj=0;jj<4;jj++)
+ positionP[ii] += (0.25f * tCrds[jj][ii-1]);
- for(ii=0;ii<4;ii++) {
- double tmp;
- tmp = gsl_vector_get(vecbStrain[0],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][0];
- gsl_vector_set(vecbStrain[0],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[1],ii) + positionP[ii]*element->tetra[tetra_I].strain[1][1];
- gsl_vector_set(vecbStrain[1],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[2],ii) + positionP[ii]*element->tetra[tetra_I].strain[2][2];
- gsl_vector_set(vecbStrain[2],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[3],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][1];
- gsl_vector_set(vecbStrain[3],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[4],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][2];
- gsl_vector_set(vecbStrain[4],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[5],ii) + positionP[ii]*element->tetra[tetra_I].strain[1][2];
- gsl_vector_set(vecbStrain[5],ii,tmp);
+ for(ii=0;ii<4;ii++) {
+ double tmp;
+ tmp = gsl_vector_get(vecbStrain[0],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][0];
+ gsl_vector_set(vecbStrain[0],ii,tmp);
+ tmp = gsl_vector_get(vecbStrain[1],ii) + positionP[ii]*element->tetra[tetra_I].strain[1][1];
+ gsl_vector_set(vecbStrain[1],ii,tmp);
+ tmp = gsl_vector_get(vecbStrain[2],ii) + positionP[ii]*element->tetra[tetra_I].strain[2][2];
+ gsl_vector_set(vecbStrain[2],ii,tmp);
+ tmp = gsl_vector_get(vecbStrain[3],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][1];
+ gsl_vector_set(vecbStrain[3],ii,tmp);
+ tmp = gsl_vector_get(vecbStrain[4],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][2];
+ gsl_vector_set(vecbStrain[4],ii,tmp);
+ tmp = gsl_vector_get(vecbStrain[5],ii) + positionP[ii]*element->tetra[tetra_I].strain[1][2];
+ gsl_vector_set(vecbStrain[5],ii,tmp);
- tmp = gsl_vector_get(vecbStress[0],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][0];
- gsl_vector_set(vecbStress[0],ii,tmp);
- tmp = gsl_vector_get(vecbStress[1],ii) + positionP[ii]*element->tetra[tetra_I].stress[1][1];
- gsl_vector_set(vecbStress[1],ii,tmp);
- tmp = gsl_vector_get(vecbStress[2],ii) + positionP[ii]*element->tetra[tetra_I].stress[2][2];
- gsl_vector_set(vecbStress[2],ii,tmp);
- tmp = gsl_vector_get(vecbStress[3],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][1];
- gsl_vector_set(vecbStress[3],ii,tmp);
- tmp = gsl_vector_get(vecbStress[4],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][2];
- gsl_vector_set(vecbStress[4],ii,tmp);
- tmp = gsl_vector_get(vecbStress[5],ii) + positionP[ii]*element->tetra[tetra_I].stress[1][2];
- gsl_vector_set(vecbStress[5],ii,tmp);
+ tmp = gsl_vector_get(vecbStress[0],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][0];
+ gsl_vector_set(vecbStress[0],ii,tmp);
+ tmp = gsl_vector_get(vecbStress[1],ii) + positionP[ii]*element->tetra[tetra_I].stress[1][1];
+ gsl_vector_set(vecbStress[1],ii,tmp);
+ tmp = gsl_vector_get(vecbStress[2],ii) + positionP[ii]*element->tetra[tetra_I].stress[2][2];
+ gsl_vector_set(vecbStress[2],ii,tmp);
+ tmp = gsl_vector_get(vecbStress[3],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][1];
+ gsl_vector_set(vecbStress[3],ii,tmp);
+ tmp = gsl_vector_get(vecbStress[4],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][2];
+ gsl_vector_set(vecbStress[4],ii,tmp);
+ tmp = gsl_vector_get(vecbStress[5],ii) + positionP[ii]*element->tetra[tetra_I].stress[1][2];
+ gsl_vector_set(vecbStress[5],ii,tmp);
- tmp = gsl_vector_get(vecbMaterial_I,ii) + positionP[ii]*((double)(element->tetra[tetra_I].material_I)+0.5);
- gsl_vector_set(vecbMaterial_I,ii,tmp);
+ tmp = gsl_vector_get(vecbMaterial_I,ii) + positionP[ii]*((double)(element->tetra[tetra_I].material_I)+0.5);
+ gsl_vector_set(vecbMaterial_I,ii,tmp);
+
+ tmp = gsl_vector_get(vecbDensity,ii) + positionP[ii]*element->tetra[tetra_I].density;
+ gsl_vector_set(vecbDensity,ii,tmp);
- for(jj=0;jj<4;jj++) {
- tmp = gsl_matrix_get(matA,ii,jj) + positionP[ii]*positionP[jj];
- gsl_matrix_set(matA,ii,jj,tmp);
- }
- } /* end of verteces of a tet. */
- } /* end of incident tets. */
- } /* if within my domain */
- } /* end of incident elements. */
+ for(jj=0;jj<4;jj++) {
+ tmp = gsl_matrix_get(matA,ii,jj) + positionP[ii]*positionP[jj];
+ gsl_matrix_set(matA,ii,jj,tmp);
+ }
+ } /* end of verteces of a tet. */
+ } /* end of incident tets. */
+ } /* if within my domain */
+ } /* end of incident elements. */
+ } /* end of patchCenterI */
/* compute parameter vectors. */
{
@@ -457,6 +521,7 @@
gsl_linalg_LU_solve (matA, p, vecbStress[i], vecaStress[i]);
}
gsl_linalg_LU_solve (matA, p, vecbMaterial_I, vecaMaterial_I);
+ gsl_linalg_LU_solve (matA, p, vecbDensity, vecaDensity);
/* printf ("x = \n"); */
/* gsl_vector_fprintf (stdout, x, "%g"); */
gsl_permutation_free (p);
@@ -478,14 +543,21 @@
node->stressSPR[j] += gsl_vector_get(vecaStress[j],i+1)*(*coord)[i];
}
}
+
node->material_ISPR = gsl_vector_get(vecaMaterial_I,0);
for(i=0;i<3;i++)
node->material_ISPR += gsl_vector_get(vecaMaterial_I,i+1)*(*coord)[i];
+ node->densitySPR = gsl_vector_get(vecaDensity,0);
+ for(i=0;i<3;i++)
+ node->densitySPR += gsl_vector_get(vecaDensity,i+1)*(*coord)[i];
+
/* free gsl vectors and matrix. */
gsl_matrix_free( matA );
gsl_vector_free( vecaMaterial_I );
gsl_vector_free( vecbMaterial_I );
+ gsl_vector_free( vecaDensity );
+ gsl_vector_free( vecbDensity );
for(i=0;i<6;i++) {
gsl_vector_free( vecaStrain[i] );
gsl_vector_free( vecaStress[i] );
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,601 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2004,
+** Pururav Thoutireddy, Center for Advanced Computing Research, Caltech.
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Pururav Thoutireddy, Center for Advanced Computing Research, Caltech. ( puru at cacr.caltech.edu)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Constitutive.c 3274 2007-03-27 20:25:29Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#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"
+#include "Register.h"
+#include <math.h>
+#include "Snac/Temperature/Temperature.h"
+#include <string.h>
+#include <assert.h>
+
+#ifndef PI
+#ifndef M_PIl
+#ifndef M_PI
+#define PI 4.0*atan(1.0)
+#else
+#define PI M_PI
+#endif
+#else
+#define PI M_PIl
+#endif
+#endif
+
+
+void SnacViscoPlastic_Constitutive( void* _context, Element_LocalIndex element_lI ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacViscoPlastic_Element* viscoplasticElement = ExtensionManager_Get( context->mesh->elementExtensionMgr, element, SnacViscoPlastic_ElementHandle );
+ SnacTemperature_Element* temperatureElement = ExtensionManager_Get( context->mesh->elementExtensionMgr, element, SnacTemperature_ElementHandle );
+ const Snac_Material* material = &context->materialProperty[element->material_I];
+
+ /*ccccc*/
+ MeshLayout* meshLayout = (MeshLayout*)context->meshLayout;
+ HexaMD* decomp = (HexaMD*)meshLayout->decomp;
+ IJK ijk;
+ Element_GlobalIndex element_gI = _MeshDecomp_Element_LocalToGlobal1D( decomp, element_lI );
+ EntryPoint* temperatureEP;
+
+ RegularMeshUtils_Element_1DTo3D( decomp, element_gI, &ijk[0], &ijk[1], &ijk[2] );
+ /*ccccc*/
+
+ temperatureEP = Context_GetEntryPoint( context, "Snac_EP_LoopElementsEnergy" );
+
+ /* If this is a ViscoPlastic material, calculate its stress. */
+ if ( material->rheology & Snac_Material_ViscoPlastic ) {
+ Tetrahedra_Index tetra_I;
+ 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;
+ double straind0,straind1,straind2,stressd0,stressd1,stressd2;
+ double Stress0[3][3];
+ double trace_strain;
+ double trace_stress;
+ double temp;
+ double vic1;
+ double vic2;
+ double VolumicStress;
+ double rviscosity=material->refvisc;
+ double rmu= material->mu;
+ double srJ2;
+ double avgTemp;
+ 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 */
+ /* values are from paper by Hall et. al., EPSL, 2003 */
+ double rstrainrate = material->refsrate;
+ double rTemp = material->reftemp;
+ double H = material->activationE; // kJ/mol
+ double srexponent = material->srexponent;
+ double srexponent1 = material->srexponent1;
+ double srexponent2 = material->srexponent2;
+ const double R=8.31448; // J/mol/K
+ /* elasto-plastic material properties */
+ double cohesion = 0.0f;
+ double frictionAngle = 0.0f;
+ double dilationAngle = 0.0f;
+ double hardening = 0.0f;
+ 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;
+
+ /* printf("Entered ViscoPlastic update \n"); */
+
+ /* Work out the plastic material properties of this element */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ double cn[3][3] = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
+ double s[3] = {0.0,0.0,0.0};
+ double alam, dep1, dep2, dep3, depm;
+ /*ccccc*/
+
+ stress = &element->tetra[tetra_I].stress;
+ strain = &element->tetra[tetra_I].strain;
+ plasticStrain = viscoplasticElement->plasticStrain[tetra_I];
+ viscosity = &viscoplasticElement->viscosity[tetra_I];
+
+ if( context->computeThermalStress ) {
+ (*stress)[0][0] += temperatureElement->thermalStress[tetra_I];
+ (*stress)[1][1] += temperatureElement->thermalStress[tetra_I];
+ (*stress)[2][2] += temperatureElement->thermalStress[tetra_I];
+ }
+
+ /* storing original stress in local array */
+ Stress0[0][0] = (*stress)[0][0];
+ Stress0[1][1] = (*stress)[1][1];
+ Stress0[2][2] = (*stress)[2][2];
+ Stress0[0][1] = (*stress)[0][1];
+ Stress0[0][2] = (*stress)[0][2];
+ Stress0[1][2] = (*stress)[1][2];
+ trace_stress = (*stress)[0][0] + (*stress)[1][1] + (*stress)[2][2];
+ /* trace_strain = element->tetra[tetra_I].volume/element->tetra[tetra_I].old_volume-1.0f; */
+ trace_strain = (*strain)[0][0] + (*strain)[1][1] + (*strain)[2][2];
+
+ /* printf(" Trace strain=%g\t Trace Stress=%g\n",trace_strain,trace_stress); */
+
+ /* Deviatoric Stresses and Strains */
+ straind0 = (*strain)[0][0] - (trace_strain) / 3.0f;
+ straind1 = (*strain)[1][1] - (trace_strain) / 3.0f;
+ straind2 = (*strain)[2][2] - (trace_strain) / 3.0f;
+
+ stressd0 = (*stress)[0][0] - (trace_stress) / 3.0f;
+ stressd1 = (*stress)[1][1] - (trace_stress) / 3.0f;
+ stressd2 = (*stress)[2][2] - (trace_stress) / 3.0f;
+
+ /* compute viscosity and add thermal stress */
+ if( temperatureEP ) {
+
+ srJ2 = sqrt(fabs(straind1*straind2+straind2*straind0+straind0*straind1 -(*strain[0][1])*(*strain[0][1])-(*strain[0][2])*(*strain[0][2])-(*strain[1][2])*(*strain[1][2])))/context->dt;
+ if(srJ2 == 0.0f) srJ2 = rstrainrate; // temporary. should be vmax/length_scale
+
+ avgTemp=0.0;
+ for(node_lI=0; node_lI<4; node_lI++) {
+ Snac_Node* contributingNode = Snac_Element_Node_P(
+ context,
+ element_lI,
+ TetraToNode[tetra_I][node_lI] );
+ SnacTemperature_Node* temperatureNodeExt = ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,contributingNode,
+ SnacTemperature_NodeHandle );
+
+ avgTemp += 0.25 * temperatureNodeExt->temperature;
+ assert( !isnan(avgTemp) && !isinf(avgTemp) );
+ }
+ // Hall et. al., 2004, G3
+ (*viscosity)= rviscosity*pow((srJ2/rstrainrate),(1./srexponent-1.))
+ *exp(H/R*(1./(avgTemp+273.15)-1./(rTemp+273.15)));
+ if((*viscosity) < material->vis_min) (*viscosity) = material->vis_min;
+ if((*viscosity) > material->vis_max) (*viscosity) = material->vis_max;
+ Journal_Firewall(
+ !isnan((*viscosity)) && !isinf((*viscosity)),
+ context->snacError,
+ "rvisc=%e Erattio=%e pow(E)=%e, dT=%e exp=%e\n",
+ rviscosity,
+ (srJ2/rstrainrate),
+ pow((srJ2/rstrainrate),
+ (1./srexponent-1.)),
+ exp(H/R*(1./(avgTemp+273.15)-1./(rTemp+273.15))),
+ (1./(avgTemp+273.15)-1./(rTemp+273.15)) );
+#if 0
+ // Lavier and Buck, JGR, 2002
+ (*viscosity) = pow(rviscosity,-1.0/srexponent1)*pow(srJ2,1.0/srexponent2-1)*exp(H/R/(avgTemp+273.15));
+#endif
+ }
+ else {
+ (*viscosity) = rviscosity;
+ Journal_Firewall(
+ !isnan((*viscosity)) && !isinf((*viscosity)),
+ context->snacError,
+ "(*viscosity) is nan or inf\n" );
+ }
+
+ /* Non dimensional parameters elastic/viscous */
+ temp = rmu / (2.0f* (*viscosity)) * context->dt;
+
+ vic1 = 1.0f - temp;
+ vic2 = 1.0f / (1.0f + temp);
+ /* printf("temp=%g\t rmu=%g\t viscosity=%g\t\n",temp,rmu,*viscosity); */
+ /* printf("trace_stress=%g\t trace_strain=%g\t vic1=%g\t vic2=%g\n",trace_stress,trace_strain,vic1,vic2); */
+ /* Deviatoric Stress Update */
+
+ stressd0 = (stressd0 * vic1 + 2.0f * rmu * straind0) * vic2 ;
+ stressd1 = (stressd1 * vic1 + 2.0f * rmu * straind1) * vic2 ;
+ stressd2 = (stressd2 * vic1 + 2.0f * rmu * straind2) * vic2 ;
+
+ (*stress)[0][1] =((*stress)[0][1] * vic1 + 2.0f * rmu * (*strain)[0][1]) * vic2;
+ (*stress)[0][2] =((*stress)[0][2] * vic1 + 2.0f * rmu * (*strain)[0][2]) * vic2;
+ (*stress)[1][2] =((*stress)[1][2] * vic1 + 2.0f * rmu * (*strain)[1][2]) * vic2;
+
+ /* Isotropic stress is elastic,
+ WARNING:volumic Strain may be better defined as
+ volumique change in the mesh */
+ VolumicStress = trace_stress / 3.0f + bulkm * trace_strain;
+
+ (*stress)[0][0] = stressd0 + VolumicStress;
+ (*stress)[1][1] = stressd1 + VolumicStress;
+ (*stress)[2][2] = stressd2 + VolumicStress;
+
+ principal_stresses(stress,s,cn);
+
+ /* compute friction and dilation angles based on accumulated plastic strain in tetrahedra */
+ /* Piece-wise linear softening */
+ /* Find current properties from linear interpolation */
+#if 0
+ /*ccccc*/
+ if(material->putSeeds && context->loop <= 1) {
+ if(ijk[1] >= decomp->elementGlobal3DCounts[1]-2)
+ if(ijk[0] == decomp->elementGlobal3DCounts[0]/2 ) {
+ //if(ijk[0] >= 10 && ijk[0] <= 14 ) {
+ viscoplasticElement->plasticStrain[tetra_I] = 1.1*material->plstrain[1];
+ plasticStrain = viscoplasticElement->plasticStrain[tetra_I];
+ fprintf(stderr,"loop=%d ijk=%d %d %d aps=%e plasticE=%e\n",
+ context->loop,ijk[0],ijk[1],ijk[2],plasticStrain,
+ viscoplasticElement->plasticStrain[tetra_I]);
+ }
+ }
+ /*ccccc*/
+#endif
+ for( i = 0; i < material->nsegments; i++ ) {
+ const double pl1 = material->plstrain[i];
+ const double pl2 = material->plstrain[i+1];
+
+ if( plasticStrain >= pl1 && plasticStrain <= pl2 ) {
+ const double tgf = (material->frictionAngle[i+1] - material->frictionAngle[i]) / (pl2 - pl1);
+ const double tgd = (material->dilationAngle[i+1] - material->dilationAngle[i]) / (pl2 - pl1);
+ const double tgc = (material->cohesion[i+1] - material->cohesion[i]) / (pl2 - pl1);
+
+ frictionAngle = material->frictionAngle[i] + tgf * (plasticStrain - pl1);
+ dilationAngle = material->dilationAngle[i] + tgd * (plasticStrain - pl1);
+ cohesion = material->cohesion[i] + tgc * (plasticStrain - pl1);
+ hardening = tgc;
+ }
+ }
+
+ if( frictionAngle > 0.0f ) {
+ tension_cutoff = material->ten_off;
+ if( frictionAngle > 0.0) {
+ tmp = cohesion / tan( frictionAngle * degrad);
+ if(tmp < tension_cutoff) tension_cutoff=tmp;
+ }
+ }
+ else {
+ /* frictionAngle < 0.0 violates second law of thermodynamics */
+ fprintf(stderr,"rank:%d elem:%d tet:%d plasticStrain=%e\n",element_lI,tetra_I,plasticStrain);
+ assert(0);
+ }
+
+ if( yieldcriterion == mohrcoulomb )
+ {
+
+ double sphi = sin( frictionAngle * degrad );
+ double spsi = sin( dilationAngle * degrad );
+ double anphi = (1.0f + sphi) / (1.0f - sphi);
+ double anpsi = (1.0f + spsi) / (1.0f - spsi);
+ double fs = s[0] - s[2] * anphi + 2 * cohesion * sqrt( anphi );
+ double ft = s[2] - tension_cutoff;
+ /* CHECK FOR COMPOSITE YIELD CRITERION */
+ ind=0;
+ if( fs < 0.0f || ft > 0.0f ) {
+ /*! Failure: shear or tensile */
+ double aP = sqrt( 1.0f + anphi * anphi ) + anphi;
+ double sP = tension_cutoff * anphi - 2 * cohesion * sqrt( anphi );
+ double h = s[2] - tension_cutoff + aP * ( s[0] - sP );
+
+ ind=1;
+
+ if( h < 0.0f ) {
+ /* !shear failure */
+ alam = fs / ( a1 - a2 * anpsi + a1 * anphi * anpsi - a2 * anphi + 2.0*sqrt(anphi)*hardening );
+ s[0] -= alam * ( a1 - a2 * anpsi );
+ s[1] -= alam * a2 * ( 1.0f - anpsi );
+ s[2] -= alam * ( a2 - a1 * anpsi );
+ dep1 = alam;
+ dep2 = 0.0f;
+ dep3 = -alam * anpsi;
+ }
+ else {
+ /* tensile failure */
+ alam = ft / a1;
+ s[0] -= alam * a2;
+ s[1] -= alam * a2;
+ s[2] -= alam * a1;
+ dep1 = 0.0f;
+ dep2 = 0.0f;
+ dep3 = alam;
+ }
+ }
+ else {
+ /* ! no failure - just elastic increment */
+
+ dep1 = 0.0f;
+ dep2 = 0.0f;
+ dep3 = 0.0f;
+ }
+ if(ind) {
+ unsigned int k,m,n;
+ /* Second invariant of accumulated plastic increment */
+ depm = ( dep1 + dep2 + dep3 ) / 3.0f;
+ viscoplasticElement->plasticStrain[tetra_I] += sqrt( 0.5f * ((dep1-depm) * (dep1-depm) + (dep2-depm) * (dep2-depm) + (dep3-depm) * (dep3-depm) + depm*depm) );
+
+ /* Stress projection back to euclidean coordinates */
+ memset( stress, 0, sizeof((*stress)) );
+ /* Resolve back to global axes */
+ 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];
+ }
+ }
+ }
+ }
+ }
+ else
+ Journal_Firewall( (0>1), "In %s: \"mohrcoulomb\" is the only available yield criterion.\n", __func__ );
+
+ depls += viscoplasticElement->plasticStrain[tetra_I]*element->tetra[tetra_I].volume;
+ totalVolume += element->tetra[tetra_I].volume;
+ }
+ /* volume-averaged accumulated plastic strain, aps */
+ viscoplasticElement->aps = depls/totalVolume;
+ }
+}
+
+
+int principal_stresses(StressTensor* stress, double sp[3], double cn[3][3])
+{
+
+ double **a,**v,*d;
+ int i,j;
+
+ a = dmatrix(1,3,1,3);
+ v = dmatrix(1,3,1,3);
+ d = dvector(1,3);
+
+ a[1][1] = (*stress)[0][0];
+ a[2][2] = (*stress)[1][1];
+ a[3][3] = (*stress)[2][2];
+ a[1][2] = (*stress)[0][1];
+ a[1][3] = (*stress)[0][2];
+ a[2][3] = (*stress)[1][2];
+ a[2][1] = a[1][2];
+ a[3][1] = a[1][3];
+ a[3][2] = a[2][3];
+
+ jacobi(a,d,v);
+
+ d[1] *= -1.0f;
+ d[2] *= -1.0f;
+ d[3] *= -1.0f;
+
+ eigsrt(d,v);
+
+ d[1] *= -1.0f;
+ d[2] *= -1.0f;
+ d[3] *= -1.0f;
+
+ for(i=0;i<3;i++) {
+ sp[i] = d[i+1];
+ for(j=0;j<3;j++) {
+ cn[i][j] = v[j+1][i+1];
+ }
+ }
+
+ free_dmatrix(a,1,3,1,3);
+ free_dmatrix(v,1,3,1,3);
+ free_dvector(d,1,3);
+
+ return(1);
+}
+
+int jacobi(double** a, double* d, double** v)
+{
+
+ int nrot = 0;
+ const unsigned int nmax = 100, n = 3;
+ double b[nmax], z[nmax], tresh,sm,g,h,t,theta,c,s,tau;
+
+ int i,j,ip,iq;
+
+ for(ip=1;ip<=n;ip++) {
+ for(iq=1;iq<=n;iq++) v[ip][iq] = 0.0f;
+ v[ip][ip] = 1.0f;
+ }
+
+ for(ip=1;ip<=n;ip++) {
+ b[ip] = d[ip] = a[ip][ip];
+ z[ip] = 0.0f;
+ }
+
+ for(i=1;i<=50;i++) {
+ sm = 0.0f;
+ for(ip=1;ip<=n-1;ip++) {
+ for(iq=ip+1;iq<=n;iq++)
+ sm += fabs(a[ip][iq]);
+ }
+ if(sm == 0.0f)
+ return(0);
+
+ if(i < 4) {
+ tresh = 0.2f * sm / ( n*n );
+ }
+ else {
+ tresh = 0.0f;
+ }
+
+ for(ip=1;ip<=n-1;ip++) {
+ for(iq=ip+1;iq<=n;iq++) {
+ g = 100.0f*fabs(a[ip][iq]);
+ if( (i > 4) && (double)(fabs(d[ip])+g) == (double)fabs(d[ip]) && (double)(fabs(d[iq])+g) == (double)fabs(d[iq]))
+ a[ip][iq] = 0.0f;
+ else if( fabs(a[ip][iq]) > tresh ) {
+ h = d[iq] - d[ip];
+ if( (fabs(h)+g) == fabs(h) )
+ t = a[ip][iq] / h;
+ else {
+ theta = 0.5f * h / (a[ip][iq]);
+ t = 1.0f / (fabs(theta) + sqrt(1.0f + theta*theta));
+ if(theta < 0.0f) t *= -1.0f;
+ }
+ c = 1.0f / sqrt(1.0f + t*t);
+ s = t * c;
+ tau = s / (1.0f + c);
+ h = t * a[ip][iq];
+ z[ip] -= h;
+ z[iq] += h;
+ d[ip] -= h;
+ d[iq] += h;
+ a[ip][iq] = 0.0f;
+ for(j=1;j<=ip-1;j++) {
+ ROTATE(a,j,ip,j,iq);
+ }
+ for(j=ip+1;j<=iq-1;j++) {
+ ROTATE(a,ip,j,j,iq);
+ }
+ for(j=iq+1;j<=n;j++) {
+ ROTATE(a,ip,j,iq,j);
+ }
+ for(j=1;j<=n;j++) {
+ ROTATE(v,j,ip,j,iq);
+ }
+
+ ++nrot;
+ }
+ }
+ }
+ for(ip=1;ip<=n;ip++) {
+ b[ip] += z[ip];
+ d[ip] = b[ip];
+ z[ip] = 0.0f;
+ }
+ }
+ assert(i<50);
+
+ return 1;
+}
+
+
+int eigsrt(double* d, double** v)
+{
+
+ const unsigned int n = 3;
+ int i,j,k;
+ double p;
+
+ for(i=1;i<n;i++) {
+ k = i;
+ p = d[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=1;j<=n;j++) {
+ p = v[j][i];
+ v[j][i] = v[j][k];
+ v[j][k] = p;
+ }
+ }
+ }
+ return(0);
+}
+
+double **dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+ long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+ double **m;
+
+ /* allocate pointers to rows */
+ m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
+ if (!m) nrerror("allocation failure 1 in matrix()");
+ m += NR_END;
+ m -= nrl;
+ /* allocate rows and set pointers to them */
+ m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
+ if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+ m[nrl] += NR_END;
+ m[nrl] -= ncl;
+ for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+ /* return pointer to array of pointers to rows */
+ return m;
+}
+
+
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+ free((FREE_ARG) (m[nrl]+ncl-NR_END));
+ free((FREE_ARG) (m+nrl-NR_END));
+}
+
+
+void nrerror(char error_text[])
+{
+ fprintf(stderr,"Numerical Recipes run-time error...\n");
+ fprintf(stderr,"%s\n",error_text);
+ fprintf(stderr,"...now exiting to system...\n");
+ exit(1);
+}
+
+
+double *dvector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+ double *v;
+ v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
+ if (!v) nrerror("allocation failure in dvector()");
+ return v-nl+NR_END;
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+ int *v;
+ v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+ if (!v) nrerror("allocation failure in ivector()");
+ return v-nl+NR_END;
+}
+
+void free_dvector(double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,70 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2004,
+** Pururav Thoutireddy, Center for Advanced Computing Research, Caltech.
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Pururav Thoutireddy, Center for Advanced Computing Research, Caltech. ( puru at cacr.caltech.edu)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+** Calculates the plastic properties for a given element.
+**
+** Assumptions:
+** None as yet.
+**
+** Comments:
+** None as yet.
+**
+** $Id: Constitutive.h 3274 2007-03-27 20:25:29Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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);
+#define SWAP(g,h) {y=(g);(g)=(h);(h)=y;}
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+
+
+#define TINY 1.0e-20
+#define RADIX 2.0
+#define NR_END 1
+#define FREE_ARG char*
+
+#ifndef __SnacViscoPlastic_h__
+#define __SnacViscoPlastic_h__
+
+ void SnacViscoPlastic_Constitutive( void* context, Element_LocalIndex element_lI );
+
+ double** dmatrix(long nrl, long nrh, long ncl, long nch);
+ double *dvector(long nl, long nh);
+ int *ivector(long nl, long nh);
+ void free_ivector(int *v, long nl, long nh);
+ void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
+ void free_dvector(double *v, long nl, long nh);
+ void nrerror(char error_text[]);
+ int jacobi(double** a, double* d, double** v);
+ int eigsrt(double* d, double** v);
+#endif /* __SnacViscoPlastic_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,112 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: ConstructExtensions.c 3234 2006-06-23 04:12:40Z LaetitiaLePourhiet $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Context.h"
+#include "Register.h"
+#include "Element.h"
+#include "ConstructExtensions.h"
+#include <assert.h>
+#include <limits.h>
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+void _SnacViscoPlastic_ConstructExtensions( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Snac_Element tmpElement;
+ SnacViscoPlastic_Element* tmpElementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ &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( 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
+
+
+ /* Prepare the dump file */
+ sprintf( tmpBuf, "%s/plStrain.%u", context->outputPath, context->rank );
+ if( (contextExt->plStrainOut = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->plStrainOut /* failed to open file for writing */ );
+ abort();
+ }
+ sprintf( tmpBuf, "%s/plStrainCP.%u", context->outputPath, context->rank );
+ if( (contextExt->plStrainCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->plStrainCheckpoint /* 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 */ );
+ abort();
+ }
+ sprintf( tmpBuf, "%s/viscosityCP.%u", context->outputPath, context->rank );
+ if( (contextExt->viscCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->viscCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ConstructExtensions.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,45 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: ConstructExtensions.h 1958 2004-08-26 19:08:58Z puru $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_ConstructExtensions_h__
+#define __SnacViscoPlastic_ConstructExtensions_h__
+
+ void _SnacViscoPlastic_ConstructExtensions( void* _context, void* data );
+
+#endif /* __SnacViscoPlastic_ConstructExtensions_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,45 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Context.c 1958 2004-08-26 19:08:58Z puru $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Context.h"
+#include <stdio.h>
+
+void SnacViscoPlastic_Context_Print( void* _context ) {
+ SnacViscoPlastic_Context* self = (SnacViscoPlastic_Context*)_context;
+
+ printf( "SnacViscoPlastic_Context:\n" );
+ printf( "\tplstrain (ptr): %p", self->plStrainOut );
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Context.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,56 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: Context.h 3104 2005-07-14 22:16:41Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_Context_h__
+#define __SnacViscoPlastic_Context_h__
+
+ /* Context Information */
+ struct _SnacViscoPlastic_Context {
+ /* Dumping */
+ FILE* plStrainOut;
+ FILE* viscOut;
+ FILE* plStrainCheckpoint;
+ FILE* viscCheckpoint;
+ FILE* plstrainTensorOut;
+ };
+
+ /* Print the contents of the context extension */
+ void SnacViscoPlastic_Context_Print( void* _context );
+
+#endif /* __SnacViscoPlastic_Context_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,49 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Element.c 1958 2004-08-26 19:08:58Z puru $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Element.h"
+#include <stdio.h>
+
+void SnacViscoPlastic_Element_Print( void* element ) {
+ SnacViscoPlastic_Element* self = (SnacViscoPlastic_Element*)element;
+ Tetrahedra_Index tetra_I;
+
+ printf( "SnacViscoPlastic_Element:\n" );
+ printf( "\tplasticStrain: " );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ printf( "%g, ", self->plasticStrain[tetra_I] );
+ }
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Element.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,54 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+** Snac plastic properties for an element.
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: Element.h 3274 2007-03-27 20:25:29Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_Element_h__
+#define __SnacViscoPlastic_Element_h__
+
+ /* Element Information */
+ struct _SnacViscoPlastic_Element {
+ Strain plasticStrain[Tetrahedra_Count];
+ double viscosity[Tetrahedra_Count];
+ Strain aps;
+ };
+
+ /* Print the contents of an Element */
+ void SnacViscoPlastic_Element_Print( void* element );
+
+#endif /* __SnacViscoPlastic_Element_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,138 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: InitialConditions.c 3261 2006-11-14 21:18:05Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Element.h"
+#include "InitialConditions.h"
+#include "Register.h"
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+#define min(a,b) ((a)<(b) ? (a):(b))
+#ifndef PATH_MAX
+#define PATH_MAX 1024
+#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;
+
+ // loop over the phase using the dictionary
+
+ Dictionary_Entry_Value* materialList = Dictionary_Get( context->dictionary, "materials" );
+ 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);
+ PhaseI++;
+ 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;
+ }
+ 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;
+ 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 );
+
+ /* read in restart file to reconstruct the previous plastic strain.*/
+ 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 );
+ const Snac_Material* material = &context->materialProperty[element->material_I];
+
+ vis_min = context->materialProperty[element->material_I].vis_min;
+ if( material->yieldcriterion == mohrcoulomb ) {
+ Tetrahedra_Index tetra_I;
+ double depls = 0.0f;
+ double totalVolume = 0.0f;
+
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ double tetraPlStrain;
+
+ /* not the actual viscosity at restartTimestep, but doesn't affect the calculation afterwards */
+ viscoplasticElement->viscosity[tetra_I] = vis_min;
+
+ fscanf( plStrainIn, "%le", &tetraPlStrain );
+ viscoplasticElement->plasticStrain[tetra_I] = tetraPlStrain;
+ depls += viscoplasticElement->plasticStrain[tetra_I]*element->tetra[tetra_I].volume;
+ totalVolume += element->tetra[tetra_I].volume;
+ }//for tets
+ /* volume-averaged accumulated plastic strain, aps */
+ viscoplasticElement->aps = depls/totalVolume;
+ }//if(mohrcoulomb)
+ }//for elements
+ if( plStrainIn )
+ fclose( plStrainIn );
+ }// if restarting.
+ else {
+ /* Set the plastic element initial conditions */
+ 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.0;
+ vis_min = context->materialProperty[element->material_I].vis_min;
+ Tetrahedra_Index tetra_I;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ viscoplasticElement->plasticStrain[tetra_I] = 0.0;
+ viscoplasticElement->viscosity[tetra_I] = vis_min;
+ }//for tets
+ }//for elements
+ }//else
+}//function
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/InitialConditions.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,45 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: InitialConditions.h 3131 2005-08-13 01:48:13Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_InitialConditions_h__
+#define __SnacViscoPlastic_InitialConditions_h__
+
+ void SnacViscoPlastic_InitialConditions( void* _context, void* data );
+
+#endif /* __SnacViscoPlastic_InitialConditions_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Make.mm
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Make.mm (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Make.mm 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,67 @@
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+##
+## Copyright (C), 2003,
+## Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+## Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+## University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+##
+## Authors:
+## Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+## Stevan M. Quenette, Visitor in Geophysics, Caltech.
+## Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+## Luc Lavier, Research Scientist, Caltech.
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+##
+## $Id: Make.mm 955 2004-03-04 18:43:37Z SteveQuenette $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+include Makefile.def
+
+PROJECT = Snac
+PACKAGE = ${def_mod}module
+
+PROJ_LIB = $(BLD_LIBDIR)/$(PACKAGE).a
+PROJ_DLL = $(BLD_LIBDIR)/$(PACKAGE).$(EXT_SO)
+PROJ_TMPDIR = $(BLD_TMPDIR)/$(PROJECT)/$(PACKAGE)
+PROJ_CLEAN += $(PROJ_LIB) $(PROJ_DLL)
+PROJ_INCDIR = $(BLD_INCDIR)/${def_inc}
+
+PROJ_SRCS = ${def_srcs}
+PROJ_CC_FLAGS += -I$(BLD_INCDIR)/$(PROJECT) -I$(BLD_INCDIR)/Snac -I$(BLD_INCDIR)/StGermain -I$(STGERMAIN_INCDIR)/ -I$(STGERMAIN_INCDIR)/StGermain `xml2-config --cflags` -DCURR_MODULE_NAME=\"${def_mod}\"
+PROJ_LIBRARIES = -L$(BLD_LIBDIR) -L$(STGERMAIN_LIBDIR)/ -lSnac -lStGermain `xml2-config --libs` $(MPI_LIBPATH) $(MPI_LIBS)
+LCCFLAGS =
+
+# I keep file lists to build a monolith .so from a set of .a's
+PROJ_OBJS_IN_TMP = ${addprefix $(PROJECT)/$(PACKAGE)/, ${addsuffix .o, ${basename $(PROJ_SRCS)}}}
+PROJ_OBJLIST = $(BLD_TMPDIR)/$(PROJECT).$(PACKAGE).objlist
+
+all: $(PROJ_LIB) DLL createObjList export
+
+DLL: product_dirs $(PROJ_OBJS)
+ $(CC) -o $(PROJ_DLL) $(PROJ_OBJS) $(COMPILER_LCC_SOFLAGS) $(LCCFLAGS) $(PROJ_LIBRARIES) $(EXTERNAL_LIBPATH) $(EXTERNAL_LIBS)
+
+
+
+createObjList::
+ @echo ${PROJ_OBJS_IN_TMP} | cat > ${PROJ_OBJLIST}
+
+#export:: export-headers
+export:: export-headers export-libraries
+EXPORT_HEADERS = ${def_hdrs}
+EXPORT_LIBS = $(PROJ_LIB) $(PROJ_DLL)
+
+check::
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Makefile.def
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Makefile.def (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Makefile.def 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,61 @@
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+##
+## Copyright (C), 2003,
+## Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+## Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+## University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+##
+## Authors:
+## Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+## Stevan M. Quenette, Visitor in Geophysics, Caltech.
+## Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+## Luc Lavier, Research Scientist, Caltech.
+## Luc Lavier, Caltech.
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+##
+## $Id: Makefile.def 1677 2004-07-20 10:30:34Z SteveQuenette $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def_mod = SnacViscoPlastic
+def_inc = Snac/ViscoPlastic
+
+def_srcs = \
+ Node.c \
+ Element.c \
+ Mesh.c \
+ Context.c \
+ ConstructExtensions.c \
+ Constitutive.c \
+ Output.c \
+ Register.c \
+ InitialConditions.c \
+ Remesh.c \
+
+def_hdrs = \
+ types.h \
+ Node.h \
+ Element.h \
+ Mesh.h \
+ Context.h \
+ ConstructExtensions.h \
+ Constitutive.h \
+ Output.h \
+ Register.h \
+ InitialConditions.h \
+ Remesh.h \
+ ViscoPlastic.h
+
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,52 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Mesh.c 3095 2005-07-13 09:50:46Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Mesh.h"
+
+#include <stdio.h>
+
+
+/*
+** Due to the rewrite this func is severely useless. Need to fix it.
+*/
+
+void SnacViscoPlastic_Mesh_Print( void* mesh, Stream* stream ) {
+ SnacViscoPlastic_Mesh* self = (SnacViscoPlastic_Mesh*)mesh;
+
+ Journal_Printf( stream, "SnacViscoPlastic_Mesh:\n" );
+
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Mesh.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,55 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+** Snac exchanger properties extensions to the mesh construct.
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: Mesh.h 3173 2005-11-21 23:47:09Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_Mesh_h__
+#define __SnacViscoPlastic_Mesh_h__
+
+
+ /* Mesh Information */
+ struct _SnacViscoPlastic_Mesh {
+ Snac_Node* newNodes;
+ };
+
+ /* Print the contents of a mesh */
+ void SnacViscoPlastic_Mesh_Print( void* mesh, Stream* stream );
+
+
+#endif /* __SnacViscoPlastic_Mesh_h__ */
+
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,47 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Node.c 1792 2004-07-30 05:42:39Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Node.h"
+#include <stdio.h>
+
+void SnacViscoPlastic_Node_Print( void* node ) {
+ SnacViscoPlastic_Node* self = (SnacViscoPlastic_Node*)node;
+ Tetrahedra_Index tetra_I;
+
+ printf( "SnacViscoPlastic_Node:\n" );
+ printf( "\tplasticStrainSPR: %g,\n", self->plStrainSPR );
+ printf( "\tviscositySPR: %g", self->viscositySPR );
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Node.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,53 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+** Snac plastic properties for an element.
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: Node.h 1720 2004-07-22 18:33:41Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_Node_h__
+#define __SnacViscoPlastic_Node_h__
+
+ /* Node Information */
+ struct _SnacViscoPlastic_Node {
+ Strain plStrainSPR;
+ double viscositySPR;
+ };
+
+ /* Print the contents of an Node */
+ void SnacViscoPlastic_Node_Print( void* element );
+
+#endif /* __SnacViscoPlastic_Node_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,235 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Output.c 3125 2005-07-25 21:32:45Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Output.h"
+#include "Context.h"
+#include "Element.h"
+#include "Register.h"
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+
+
+void _SnacViscoPlastic_WritePlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+
+ if( isTimeToDump( context ) )
+ _SnacViscoPlastic_DumpPlasticStrain( context );
+
+ if( isTimeToCheckpoint( context ) )
+ _SnacViscoPlastic_CheckpointPlasticStrain( context );
+}
+
+
+void _SnacViscoPlastic_DumpPlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ 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 );
+ float plasticStrain = elementExt->aps;
+ fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainOut );
+ }
+ fflush( contextExt->plStrainOut );
+}
+
+
+void _SnacViscoPlastic_CheckpointPlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ 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 );
+ Tetrahedra_Index tetra_I;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float plasticStrain = elementExt->plasticStrain[tetra_I];
+ fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainCheckpoint );
+ }
+ fflush( contextExt->plStrainCheckpoint );
+ }
+}
+
+
+void _SnacViscoPlastic_WriteViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+
+ if( isTimeToDump( context ) )
+ _SnacViscoPlastic_DumpViscosity( context );
+
+ if( isTimeToCheckpoint( context ) )
+ _SnacViscoPlastic_CheckpointViscosity( context );
+}
+
+
+void _SnacViscoPlastic_DumpViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ 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 );
+ Tetrahedra_Index tetra_I;
+ double viscosity = 0.0f;
+ float logviscosity = 0.0f;
+
+ /* Take average of tetra viscosity for the element */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
+ viscosity += elementExt->viscosity[tetra_I]/Tetrahedra_Count;
+ assert(viscosity>0.0);
+ logviscosity = log10(viscosity);
+ fwrite( &logviscosity, sizeof(float), 1, contextExt->viscOut );
+ }
+ fflush( contextExt->viscOut );
+}
+
+
+void _SnacViscoPlastic_CheckpointViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ 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 );
+ Tetrahedra_Index tetra_I;
+ double viscosity = 0.0f;
+ float logviscosity = 0.0f;
+
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ viscosity = elementExt->viscosity[tetra_I];
+ assert(viscosity>0.0);
+ logviscosity = log10(viscosity);
+ fwrite( &logviscosity, sizeof(float), 1, contextExt->viscCheckpoint );
+ }
+ }
+ 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
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Output.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,43 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Output.h 3104 2005-07-14 22:16:41Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef _SnacViscoPlastic_Output_
+#define _SnacViscoPlastic_Output_
+
+ void _SnacViscoPlastic_WritePlasticStrain( void* _context );
+ void _SnacViscoPlastic_WriteViscosity( void* _context );
+ void _SnacViscoPlastic_DumpPlasticStrain( void* _context );
+ void _SnacViscoPlastic_CheckpointPlasticStrain( void* _context );
+ void _SnacViscoPlastic_DumpViscosity( void* _context );
+ void _SnacViscoPlastic_CheckpointViscosity( void* _context );
+ void _SnacViscoPlastic_DumpPlasticStrainTensor( void* _context );
+
+#endif
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,168 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Register.c 3243 2006-10-12 09:04:00Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Node.h"
+#include "Element.h"
+#include "Mesh.h"
+#include "Context.h"
+#include "Constitutive.h"
+#include "ConstructExtensions.h"
+#include "Output.h"
+#include "InitialConditions.h"
+#include "Remesh.h"
+#include "Register.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 ) {
+ return PluginsManager_Submit( pluginsMgr,
+ SnacViscoPlastic_Type,
+ "0",
+ _SnacViscoPlastic_DefaultNew );
+}
+
+
+void* _SnacViscoPlastic_DefaultNew( Name name ) {
+ return _Codelet_New( sizeof(Codelet),
+ SnacViscoPlastic_Type,
+ _Codelet_Delete,
+ _Codelet_Print,
+ _Codelet_Copy,
+ _SnacViscoPlastic_DefaultNew,
+ _SnacViscoPlastic_Construct,
+ _Codelet_Build,
+ _Codelet_Initialise,
+ _Codelet_Execute,
+ _Codelet_Destroy,
+ name );
+}
+
+
+void _SnacViscoPlastic_Construct( void* component, Stg_ComponentFactory* cf, void* data ) {
+ Snac_Context* context;
+ EntryPoint* recoverNodeEP;
+ EntryPoint* interpolateNodeEP;
+ EntryPoint* interpolateElementEP;
+
+ /* Retrieve context. */
+ context = (Snac_Context*)Stg_ComponentFactory_ConstructByName( cf, "context", Snac_Context, True, data );
+
+ #ifdef DEBUG
+ 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 ),
+ SnacViscoPlastic_Type,
+ SnacViscoPlastic_Constitutive,
+ SnacViscoPlastic_Type );
+ EntryPoint_InsertBefore(
+ Context_GetEntryPoint( context, AbstractContext_EP_Initialise ),
+ "SnacTimeStepZero",
+ SnacViscoPlastic_Type,
+ SnacViscoPlastic_InitialConditions,
+ SnacViscoPlastic_Type );
+ EntryPoint_Prepend( /* Dump the initial plastic strain */
+ Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
+ "SnacViscoPlastic_WritePlasticStrain",
+ _SnacViscoPlastic_WritePlasticStrain,
+ SnacViscoPlastic_Type );
+ EntryPoint_Append( /* and dump each loop */
+ Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
+ "SnacViscoPlastic_WritePlasticStrain",
+ _SnacViscoPlastic_WritePlasticStrain,
+ SnacViscoPlastic_Type );
+ EntryPoint_Prepend( /* Dump the initial viscosity */
+ Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
+ "SnacViscoPlastic_WriteViscosity",
+ _SnacViscoPlastic_WriteViscosity,
+ SnacViscoPlastic_Type );
+ EntryPoint_Append( /* and dump visocisty each loop */
+ Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
+ "SnacViscoPlastic_WriteViscosity",
+ _SnacViscoPlastic_WriteViscosity,
+ SnacViscoPlastic_Type );
+
+ /* Add extensions to the interpolate element entry point, but it will only exist if the remesher is loaded. */
+ recoverNodeEP = Context_GetEntryPoint( context, "SnacRemesher_EP_RecoverNode" );
+ if( recoverNodeEP ) {
+ EntryPoint_Append(
+ recoverNodeEP,
+ SnacViscoPlastic_Type,
+ _SnacViscoPlastic_RecoverNode,
+ SnacViscoPlastic_Type );
+ }
+ interpolateNodeEP = Context_GetEntryPoint( context, "SnacRemesher_EP_InterpolateNode" );
+ if( interpolateNodeEP ) {
+ EntryPoint_Append(
+ interpolateNodeEP,
+ SnacViscoPlastic_Type,
+ _SnacViscoPlastic_InterpolateNode,
+ SnacViscoPlastic_Type );
+ }
+ interpolateElementEP = Context_GetEntryPoint( context, "SnacRemesher_EP_InterpolateElement" );
+ if( interpolateElementEP ) {
+ EntryPoint_Append(
+ interpolateElementEP,
+ SnacViscoPlastic_Type,
+ _SnacViscoPlastic_InterpolateElement,
+ SnacViscoPlastic_Type );
+ }
+
+ /* Construct. */
+ _SnacViscoPlastic_ConstructExtensions( context, data );
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Register.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,61 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+** Only one can be registered at the time (limitation).
+**
+** Comments:
+**
+** $Id: Register.h 3243 2006-10-12 09:04:00Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_Register_h__
+#define __SnacViscoPlastic_Register_h__
+
+ /* 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 );
+
+ void* _SnacViscoPlastic_DefaultNew( Name name );
+
+ void _SnacViscoPlastic_Construct( void* component, Stg_ComponentFactory* cf, void* data );
+
+#endif /* __SnacViscoPlastic_Register_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.c 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,324 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Pururav Thoutireddy,
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Pururav Thoutireddy, Staff Scientist, Caltech
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: Remesh.c 3250 2006-10-23 06:15:18Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Node.h"
+#include "Element.h"
+#include "Mesh.h"
+#include "Remesh.h"
+#include "Register.h"
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+
+
+void _SnacViscoPlastic_RecoverNode( void* _context, unsigned nodeInd )
+{
+ Snac_Context* context = (Snac_Context*)_context;
+ Mesh* mesh = context->mesh;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+ NodeLayout* nLayout = layout->nodeLayout;
+
+ Snac_Node* node = Snac_Node_At( context, nodeInd );
+ SnacViscoPlastic_Node* nodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr,
+ node,
+ SnacViscoPlastic_NodeHandle );
+ Coord* coord = Snac_NodeCoord_P( context, nodeInd );
+
+ Index nodeElementCount = context->mesh->nodeElementCountTbl[nodeInd];
+ Index nodeElement_I;
+ Element_DomainIndex* elements;
+ gsl_matrix* matA;
+ gsl_vector* vecaPlStrain;
+ gsl_vector* vecbPlStrain;
+ gsl_vector* vecaViscosity;
+ gsl_vector* vecbViscosity;
+ Index i;
+ IJK ijk;
+ Node_GlobalIndex node_gI = _MeshDecomp_Node_LocalToGlobal1D( decomp, nodeInd );
+ Node_GlobalIndex gNodeI = decomp->nodeGlobal3DCounts[0];
+ Node_GlobalIndex gNodeJ = decomp->nodeGlobal3DCounts[1];
+ Node_GlobalIndex gNodeK = decomp->nodeGlobal3DCounts[2];
+ Node_GlobalIndex intNode_gI;
+ Node_DomainIndex intNode_lI;
+ unsigned int isBoundaryNode=0, patchCenterNum=1, patchCenterI, patchCenterList[2];
+
+ RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
+
+ /* If boundary node, find a (topologically?) nearest interior node.
+ Loosely following
+ Khoei and Gharehbaghi, 2007, The superconvergent patch recovery techinique and data transfer operators in 3D plasticity problems,
+ Finite Elements in Analysis and Design 43 (2007) 630-- 648. */
+ if( (gNodeI>2) && (ijk[0]==0) ) {
+ ijk[0] += 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeI>2) && (ijk[0]==decomp->nodeGlobal3DCounts[0]-1) ) {
+ ijk[0] -= 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeJ>2) && (ijk[1]==0) ) {
+ ijk[1] += 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeJ>2) && (ijk[1]==decomp->nodeGlobal3DCounts[1]-1) ) {
+ ijk[1] -= 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeK>2) && (ijk[2]==0) ) {
+ ijk[2] += 1;
+ isBoundaryNode = 1;
+ }
+ if( (gNodeK>2) && (ijk[2]==decomp->nodeGlobal3DCounts[2]-1) ) {
+ ijk[2] -= 1;
+ isBoundaryNode = 1;
+ }
+ /* node_lI itself always becomes a patch center,
+ and if the current node is a boundary node, find an interior node and add it to the patch center list. */
+ patchCenterList[0] = nodeInd;
+ if( isBoundaryNode ) {
+ patchCenterNum=2;
+ intNode_gI = ijk[0]+gNodeI*ijk[1]+gNodeI*gNodeJ*ijk[2];
+ patchCenterList[1] = Mesh_NodeMapGlobalToLocal( mesh, intNode_gI );
+ }
+
+ /* initialize gsl vectors and matrix. */
+ matA = gsl_matrix_alloc(4,4); gsl_matrix_set_zero( matA );
+ vecaPlStrain = gsl_vector_alloc(4); gsl_vector_set_zero( vecaPlStrain );
+ vecbPlStrain = gsl_vector_alloc(4); gsl_vector_set_zero( vecbPlStrain );
+ vecaViscosity = gsl_vector_alloc(4); gsl_vector_set_zero( vecaViscosity );
+ vecbViscosity = gsl_vector_alloc(4); gsl_vector_set_zero( vecbViscosity );
+
+ /* For each patch center */
+ for( patchCenterI=0; patchCenterI < patchCenterNum; patchCenterI++ ) {
+ /* For each incident element, find inicident tets. */
+ for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
+ Element_DomainIndex element_dI = context->mesh->nodeElementTbl[patchCenterList[patchCenterI]][nodeElement_I];
+
+ if( element_dI < mesh->elementDomainCount ) {
+ Index elementTetra_I;
+ Snac_Element* element = Snac_Element_At( context, element_dI );
+ SnacViscoPlastic_Element* elementExt = ExtensionManager_Get( context->mesh->elementExtensionMgr,
+ element,
+ SnacViscoPlastic_ElementHandle );
+
+ /* Extract the element's node indices. Note that there should always be eight of these. */
+ {
+ Element_GlobalIndex element_gI;
+
+ elements = Memory_Alloc_Array( Node_DomainIndex, nodeElementCount, "SnacRemesher" );
+ element_gI = Mesh_ElementMapDomainToGlobal( mesh, element_dI );
+ nLayout->buildElementNodes( nLayout, element_gI, elements );
+ }
+
+ /* Convert global node indices to domain. */
+ {
+ unsigned eltNode_i;
+
+ for( eltNode_i = 0; eltNode_i < nodeElementCount; eltNode_i++ ) {
+ elements[eltNode_i] = Mesh_NodeMapGlobalToDomain( mesh, elements[eltNode_i] );
+ }
+ }
+
+ /* For each incident tetrahedron in the incident element,
+ add up contributions to P, A, and b as in Zienkiewicz and Zhu (1992), p. 1336 */
+ for( elementTetra_I = 0; elementTetra_I < Node_Element_Tetrahedra_Count;elementTetra_I++ ) {
+ Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][elementTetra_I];
+ Coord tCrds[4];
+ double positionP[4] = {1.0,0.0,0.0,0.0};
+ Index ii,jj;
+
+ /* Extract the tetrahedron's coordinates. */
+ Vector_Set( tCrds[0], mesh->nodeCoord[elements[TetraToNode[tetra_I][0]]] );
+ Vector_Set( tCrds[1], mesh->nodeCoord[elements[TetraToNode[tetra_I][1]]] );
+ Vector_Set( tCrds[2], mesh->nodeCoord[elements[TetraToNode[tetra_I][2]]] );
+ Vector_Set( tCrds[3], mesh->nodeCoord[elements[TetraToNode[tetra_I][3]]] );
+
+ for(ii=1;ii<4;ii++)
+ for(jj=0;jj<4;jj++)
+ positionP[ii] += (0.25f * tCrds[jj][ii-1]);
+
+ for(ii=0;ii<4;ii++) {
+ double tmp = gsl_vector_get(vecbPlStrain,ii) + positionP[ii]*elementExt->plasticStrain[tetra_I];
+ gsl_vector_set(vecbPlStrain,ii,tmp);
+ /* tmp = gsl_vector_get(vecbViscosity,ii) + positionP[ii]*elementExt->viscosity[tetra_I]; */
+ /* gsl_vector_set(vecbViscosity,ii,tmp); */
+
+ for(jj=0;jj<4;jj++) {
+ tmp = gsl_matrix_get(matA,ii,jj) + positionP[ii]*positionP[jj];
+ gsl_matrix_set(matA,ii,jj,tmp);
+ }
+ } /* end of verteces of a tet. */
+/* if(context->rank==0 && nodeInd==5375) */
+/* fprintf(stderr,"rank:%d node:%d elem:%d tet:%d plStrain: %e\n", */
+/* context->rank,nodeInd,element_dI,tetra_I,elementExt->plasticStrain[tetra_I]); */
+ } /* end of incident tets. */
+ } /* if within my domain */
+ } /* end of incident elements. */
+ } /* end of patchCenterI */
+
+
+ /* compute parameter vectors. */
+ {
+ int s;
+ gsl_permutation * p = gsl_permutation_alloc (4);
+
+ gsl_linalg_LU_decomp (matA, p, &s);
+ gsl_linalg_LU_solve (matA, p, vecbPlStrain, vecaPlStrain);
+/* gsl_linalg_LU_solve (matA, p, vecbViscosity, vecaViscosity); */
+ gsl_permutation_free (p);
+ }
+
+ /* Recover using the parameter vectors. */
+ nodeExt->plStrainSPR = gsl_vector_get(vecaPlStrain,0);
+/* nodeExt->viscositySPR = gsl_vector_get(vecaViscosity,0); */
+ for(i=0;i<3;i++) {
+ nodeExt->plStrainSPR += gsl_vector_get(vecaPlStrain,i+1)*(*coord)[i];
+/* nodeExt->viscositySPR += gsl_vector_get(vecaViscosity,i+1)*(*coord)[i]; */
+ }
+
+/* if(context->rank==0) */
+/* fprintf(stderr,"rank:%d node:%d plStrainSPR:%e\n",context->rank,nodeInd, */
+/* nodeExt->plStrainSPR); */
+
+ /* free gsl vectors and matrix. */
+ gsl_matrix_free( matA );
+ gsl_vector_free( vecaPlStrain );
+ gsl_vector_free( vecbPlStrain );
+ gsl_vector_free( vecaViscosity );
+ gsl_vector_free( vecbViscosity );
+
+ /* Free the element node array. */
+ FreeArray( elements );
+
+ /* end of recovery. */
+}
+
+
+void _SnacViscoPlastic_InterpolateNode( void* _context,
+ unsigned nodeInd, unsigned elementInd, unsigned tetInd,
+ unsigned* tetNodeInds, double* weights,
+ Snac_Node* dstNodes )
+{
+ Snac_Context* context = (Snac_Context*)_context;
+ Mesh* mesh = context->mesh;
+ NodeLayout* nLayout = mesh->layout->nodeLayout;
+ Snac_Node* dstNode = (Snac_Node*)ExtensionManager_At( context->mesh->nodeExtensionMgr,
+ dstNodes,
+ nodeInd );
+ SnacViscoPlastic_Node* dstNodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr,
+ dstNode,
+ SnacViscoPlastic_NodeHandle );
+ Index tetNode_i;
+
+/* SnacViscoPlastic_Node_Print( dstNodeExt ); */
+
+ /* Extract the element's node indices. Note that there should always be eight of these. */
+ Node_DomainIndex* eltNodes;
+ unsigned int nEltNodes = 8;
+ {
+ Element_GlobalIndex gEltInd;
+
+ eltNodes = Memory_Alloc_Array( Node_DomainIndex, nEltNodes, "SnacViscoPlastic" );
+ gEltInd = Mesh_ElementMapDomainToGlobal( mesh, elementInd );
+ nLayout->buildElementNodes( nLayout, gEltInd, eltNodes );
+ }
+
+ /* Convert global node indices to local. */
+ {
+ unsigned eltNode_i;
+
+ for( eltNode_i = 0; eltNode_i < nEltNodes; eltNode_i++ ) {
+ eltNodes[eltNode_i] = Mesh_NodeMapGlobalToDomain( mesh, eltNodes[eltNode_i] );
+ }
+ }
+
+ /* Clear the plastic strain. */
+ dstNodeExt->plStrainSPR = 0.0;
+/* dstNodeExt->viscositySPR = 0.0; */
+
+ /* Loop over each contributing node. */
+/* fprintf(stderr,"node=%d ",nodeInd ); */
+ for( tetNode_i = 0; tetNode_i < 4; tetNode_i++ ) {
+
+ /* Where is this contibution coming from? */
+ Snac_Node* srcNode = Snac_Node_At( context, eltNodes[tetNodeInds[tetNode_i]] );
+ SnacViscoPlastic_Node* srcNodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr,
+ srcNode,
+ SnacViscoPlastic_NodeHandle );
+
+ dstNodeExt->plStrainSPR += srcNodeExt->plStrainSPR * weights[tetNode_i];
+/* dstNodeExt->viscositySPR += srcNodeExt->viscositySPR * weights[tetNode_i]; */
+
+/* fprintf(stderr,"%d:(%d %e %e) ",tetNode_i, eltNodes[tetNodeInds[tetNode_i]], */
+/* srcNodeExt->plStrainSPR, weights[tetNode_i]); */
+ }
+ if( dstNodeExt->plStrainSPR < 0 ) dstNodeExt->plStrainSPR = 0.0;
+/* fprintf(stderr," %e\n", dstNodeExt->plStrainSPR); */
+
+}
+
+
+void _SnacViscoPlastic_InterpolateElement( void* _context,
+ Element_LocalIndex dstEltInd,
+ Tetrahedra_Index dstTetInd )
+{
+ Snac_Context* context = (Snac_Context*)_context;
+ Snac_Element* dstElt = Snac_Element_At( context, dstEltInd );
+ SnacViscoPlastic_Element* dstEltExt = ExtensionManager_Get( context->mesh->elementExtensionMgr,
+ dstElt,
+ SnacViscoPlastic_ElementHandle );
+ double plasticStrain = 0.0;
+/* double viscosity = 0.0; */
+ Index i;
+
+/* fprintf(stderr,"element:%d tet=%d\n",dstEltInd,dstTetInd); */
+ for(i=0;i<4;i++) {
+ Snac_Node* node = Snac_Element_Node_P( context, dstEltInd, TetraToNode[dstTetInd][i] );
+ SnacViscoPlastic_Node* nodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr,
+ node,
+ SnacViscoPlastic_NodeHandle );
+ plasticStrain += 0.25f*nodeExt->plStrainSPR;
+/* viscosity += 0.25f*nodeExt->viscositySPR; */
+
+/* if( nodeExt->plStrainSPR != 0.0 ) */
+/* fprintf(stderr,"\t%d:(%d %e)\n", i, TetraToNode[dstTetInd][i], nodeExt->plStrainSPR); */
+ }
+ dstEltExt->plasticStrain[dstTetInd] = plasticStrain;
+/* dstEltExt->viscosity[dstTetInd] = viscosity; */
+}
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Remesh.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,50 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: Remesh.h 3250 2006-10-23 06:15:18Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_Remesh_h__
+#define __SnacViscoPlastic_Remesh_h__
+
+void _SnacViscoPlastic_RecoverNode( void* context, unsigned nodeInd );
+void _SnacViscoPlastic_InterpolateNode( void* context,
+ unsigned nodeInd, unsigned elementInd, unsigned tetInd,
+ unsigned* tetNodeInds, double* weights,
+ Snac_Node* dstNodes );
+void _SnacViscoPlastic_InterpolateElement( void* _context, unsigned dstElementInd, unsigned dstTetInd );
+
+#endif /* __SnacViscoPlastic_Remesh_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ViscoPlastic.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ViscoPlastic.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/ViscoPlastic.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,57 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+** Only one can be registered at the time (limitation).
+**
+** Comments:
+**
+** $Id: ViscoPlastic.h 2202 2004-10-19 08:43:09Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_h__
+#define __SnacViscoPlastic_h__
+
+#include "types.h"
+#include "Node.h"
+#include "Element.h"
+#include "Mesh.h"
+#include "Context.h"
+#include "InitialConditions.h"
+#include "Constitutive.h"
+#include "Output.h"
+#include "ConstructExtensions.h"
+#include "ViscoPlastic.h"
+#include "Remesh.h"
+#include "Register.h"
+
+#endif /* __SnacViscoPlastic_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/makefile
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/makefile (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/makefile 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,53 @@
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+##
+## Copyright (C), 2003,
+## Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+## Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+## University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+##
+## Authors:
+## Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+## Stevan M. Quenette, Visitor in Geophysics, Caltech.
+## Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+## Luc Lavier, Research Scientist, Caltech.
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+##
+## $Id: Makefile.rules 1792 2004-07-30 05:42:39Z SteveQuenette $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+# obtain defaults for required variables according to system and project location, and then run the build.
+ifndef PROJ_ROOT
+ PROJ_ROOT=../..
+endif
+include ${PROJ_ROOT}/Makefile.system
+
+include Makefile.def
+
+mod = ${def_mod}
+includes = ${def_inc}
+
+SRCS = ${def_srcs}
+
+HDRS = ${def_hdrs}
+
+PROJ_LIBS = ${def_libs}
+EXTERNAL_LIBS = -L${STGERMAIN_LIBDIR} -lSnac -lStGermain
+EXTERNAL_INCLUDES = -I${STGERMAIN_INCDIR}/StGermain -I${STGERMAIN_INCDIR}
+
+packages = MPI XML MATH
+
+include ${PROJ_ROOT}/Makefile.vmake
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/types.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/types.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/types.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,52 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+** Plastic types.
+**
+** Assumptions:
+** None as yet.
+**
+** Comments:
+** None as yet.
+**
+** $Id: types.h 1958 2004-08-26 19:08:58Z puru $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_types_h__
+#define __SnacViscoPlastic_types_h__
+
+ /* Plastic */
+ typedef struct _SnacViscoPlastic_Node SnacViscoPlastic_Node;
+ typedef struct _SnacViscoPlastic_Element SnacViscoPlastic_Element;
+ typedef struct _SnacViscoPlastic_Mesh SnacViscoPlastic_Mesh;
+ typedef struct _SnacViscoPlastic_Context SnacViscoPlastic_Context;
+
+#endif /* __SnacViscoPlastic_types_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/units.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/units.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/units.h 2010-03-31 15:24:03 UTC (rev 16473)
@@ -0,0 +1,47 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+** None as yet.
+**
+** Comments:
+** None as yet.
+**
+** $Id: units.h 1968 2004-08-27 21:55:29Z puru $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacViscoPlastic_units_h__
+#define __SnacViscoPlastic_units_h__
+
+ typedef double ViscoPlastic;
+
+#endif /* __SnacViscoPlastic_units_h__ */
More information about the CIG-COMMITS
mailing list