[cig-commits] r19193 - long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI
echoi at geodynamics.org
echoi at geodynamics.org
Mon Nov 14 07:57:18 PST 2011
Author: echoi
Date: 2011-11-14 07:57:18 -0800 (Mon, 14 Nov 2011)
New Revision: 19193
Added:
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Node.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Node.h
Modified:
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Element.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/types.h
Log:
Made the whole plugin consistent with viscoplasticSPR, which has improved output routines.
-Some unnecessary structures such as mesh and node extensions will be removed after testing.
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Constitutive.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -1,35 +1,35 @@
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- **
- ** 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 $
- **
- **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+**
+** 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>
@@ -83,15 +83,13 @@
Node_LocalIndex node_lI;
/* viscoplastic material properties */
- const double bulkm = material->lambda + 2.0f * material->mu/3.0f;
- StressTensor* stress;
- StrainTensor* strain;
- StrainTensor* plasticstrainTensor;
- StrainTensor* creepstrainTensor;
- Strain plasticStrain;
- ViscoPlastic* viscosity;
+ 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 Stress0[3][3];
double trace_strain;
double trace_stress;
double temp;
@@ -99,10 +97,10 @@
double vic2;
double VolumicStress;
double rviscosity=material->refvisc;
- double rmu= material->mu;
+ double rmu= material->mu;
double srJ2;
double avgTemp;
- plModel yieldcriterion=material->yieldcriterion;
+ plModel yieldcriterion=material->yieldcriterion;
/* For now reference values of viscosity, second invariant of deviatoric */
/* strain rate and reference temperature are being hard wired ( these specific */
@@ -113,47 +111,35 @@
double srexponent = material->srexponent;
double srexponent1 = material->srexponent1;
double srexponent2 = material->srexponent2;
- const double R=8.31448; // J/mol/K
+ 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 sphi=0.0;
- double spsi=0.0;
- double st=0.0;
- double tolerance;
+ 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;
+ 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;
- double anphi=0.0,anpsi=0.0;/* for mohr-coulomb model */
- double qphi=0.0,qpsi=0.0,kphi=0.0; /* for drucker-prager model */
+ unsigned int i;
+ double tmp=0.0;
const double a1 = material->lambda + 2.0f * material->mu ;
const double a2 = material->lambda ;
- double J2,I1; /* Stress Invariants for drucker-prager model */
- int ind=0,ite,mite;
- int shearfailure,tensilefailure;
- double tolerf,previousfs,previousft,palam;
+ 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, dep0,dep1, dep2, dep3, dep01, dep02, dep12, depm;
- double fs,ft,fs0,ft0,aP,sP,h;
- unsigned int k, m, n;
+ 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];
- plasticstrainTensor = &viscoplasticElement->plasticstrainTensor[tetra_I];
- creepstrainTensor = &viscoplasticElement->creepstrainTensor[tetra_I];
if( context->computeThermalStress ) {
(*stress)[0][0] += temperatureElement->thermalStress[tetra_I];
@@ -192,12 +178,12 @@
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] );
+ context,
+ element_lI,
+ TetraToNode[tetra_I][node_lI] );
SnacTemperature_Node* temperatureNodeExt = ExtensionManager_Get(
- context->mesh->nodeExtensionMgr,contributingNode,
- SnacTemperature_NodeHandle );
+ context->mesh->nodeExtensionMgr,contributingNode,
+ SnacTemperature_NodeHandle );
avgTemp += 0.25 * temperatureNodeExt->temperature;
assert( !isnan(avgTemp) && !isinf(avgTemp) );
@@ -208,15 +194,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)) );
+ !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));
@@ -225,9 +211,9 @@
else {
(*viscosity) = rviscosity;
Journal_Firewall(
- !isnan((*viscosity)) && !isinf((*viscosity)),
- context->snacError,
- "(*viscosity) is nan or inf\n" );
+ !isnan((*viscosity)) && !isinf((*viscosity)),
+ context->snacError,
+ "(*viscosity) is nan or inf\n" );
}
/* Non dimensional parameters elastic/viscous */
@@ -248,8 +234,8 @@
(*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 */
+ 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;
@@ -273,54 +259,66 @@
context->loop,ijk[0],ijk[1],ijk[2],plasticStrain,
viscoplasticElement->plasticStrain[tetra_I]);
}
- }
- /*ccccc*/
+ }
+ /*ccccc*/
#endif
- for( i = 0; i < material->nsegments; i++ ) {
- const double pl1 = material->plstrain[i];
- const double pl2 = material->plstrain[i+1];
+ 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);
+ 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;
- }
+ 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 ) {
- st = material->ten_off;
- if( frictionAngle > 0.0) {
- tmp = cohesion / tan( frictionAngle * degrad);
- if(tmp < st)st=tmp;
- }
+ }
+
+ 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 {
+ if( plasticStrain < 0.0 ) {
+ viscoplasticElement->plasticStrain[tetra_I] = 0.0;
+ plasticStrain = viscoplasticElement->plasticStrain[tetra_I];
+ fprintf(stderr,"Warning: negative plastic strain. Setting to zero, but check if remesher is on and this happended for an external tet. rank:%d elem:%d tet:%d plasticStrain=%e frictionAngle=%e\n",context->rank,element_lI,tetra_I,plasticStrain,frictionAngle);
+ frictionAngle = material->frictionAngle[0];
+ dilationAngle = material->dilationAngle[0];
+ cohesion = material->cohesion[0];
+ }
else {
/* frictionAngle < 0.0 violates second law of thermodynamics */
+ fprintf(stderr,"Error due to an unknown reason: rank:%d elem:%d tet:%d plasticStrain=%e frictionAngle=%e\n",context->rank,element_lI,tetra_I,plasticStrain,frictionAngle);
assert(0);
}
+ }
- if( yieldcriterion == mohrcoulomb )
+ if( yieldcriterion == mohrcoulomb )
{
- sphi = sin( frictionAngle * degrad );
- spsi = sin( dilationAngle * degrad );
- anphi = (1.0f + sphi) / (1.0f - sphi);
- anpsi = (1.0f + spsi) / (1.0f - spsi);
- fs = s[0] - s[2] * anphi + 2 * cohesion * sqrt( anphi );
- ft = s[2] - st;
+
+ 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 ) {
- ind=1;
-
/*! 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 );
- aP = sqrt( 1.0f + anphi * anphi ) + anphi;
- sP = st * anphi - 2 * cohesion * sqrt( anphi );
- h = s[2] - st + aP * ( s[0] - sP );
+ ind=1;
if( h < 0.0f ) {
/* !shear failure */
@@ -350,8 +348,8 @@
dep2 = 0.0f;
dep3 = 0.0f;
}
- if(ind)
- {
+ 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) );
@@ -368,391 +366,251 @@
}
}
}
+ else
+ Journal_Firewall( (0>1), "In %s: \"mohrcoulomb\" is the only available yield criterion.\n", __func__ );
- if( yieldcriterion == druckerprager )
- {
- qphi = tan( frictionAngle * degrad);
- qpsi = tan( dilationAngle * degrad);
- kphi = cohesion;
- /* qphi = 0.5*(6*sin(frictionAngle * degrad)/(sqrt(3.)*(3.-sin(frictionAngle * degrad)))+6*sin(frictionAngle * degrad)/(sqrt(3.)*(3.+sin(frictionAngle * degrad)))); */
- /* kphi = 0.5*(6*cohesion*cos(frictionAngle * degrad)/(sqrt(3.)*(3.-sin(frictionAngle * degrad)))+6*cohesion*cos(frictionAngle * degrad)/(sqrt(3.)*(3.+sin(frictionAngle * degrad)))); */
- /* qpsi = 0.5*(6*sin(dilationAngle * degrad)/(sqrt(3.)*(3.-sin(dilationAngle * degrad)))+6*sin(dilationAngle * degrad)/(sqrt(3.)*(3.+sin(dilationAngle * degrad)))); */
+ /* linear healing: applied whether this tet has yielded or not.
+ Parameters are hardwired for now, but should be given through an input file. */
+ /* viscoplasticElement->plasticStrain[tetra_I] *= (1.0/(1.0+context->dt/1.0e+12)); */
+ viscoplasticElement->plasticStrain[tetra_I] *= (1.0/(1.0+context->dt/(ind?1.0e+13:5.0e+11)));
- I1 = ((*stress)[0][0]+(*stress)[1][1]+(*stress)[2][2])/3.;
- J2 = sqrt(0.5*(stressd0*stressd0+stressd1*stressd1+stressd2*stressd2)
- +(*stress)[0][1]*(*stress)[0][1]+(*stress)[0][2]*(*stress)[0][2]+(*stress)[1][2]*(*stress)[1][2]);
- fs = J2 + I1 * qphi - kphi;
- ft = I1 - st;
+ 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;
+ }
+}
- /* CHECK FOR COMPOSITE YIELD CRITERION */
- ind=0;
- if( fs > 0.0f || ft > 0.0f )
- {
- fs0=fs;
- ft0=ft;
- ind=1;
- dep0=dep1=dep2=dep01=dep02=dep12=0.0;
- tolerance = material->constitutivetolerance;
- if(fs > ft)
- {
- tolerf=tolerance*fs;
- }
- else
- {
- tolerf=tolerance*ft;
- }
+int principal_stresses(StressTensor* stress, double sp[3], double cn[3][3])
+{
- printf(" ENTERED PLASTIC UPDATE \n");
- mite = material->maxiterations;
- previousfs=fs;
- previousft=ft;
- for(ite=0;ite<mite;ite++)
- {
- /*! Failure: shear or tensile */
+ double **a,**v,*d;
+ int i,j;
- aP = sqrt( 1.0f + qphi * qphi ) - qphi;
- sP = kphi - qphi*st;
- h = J2 - sP - aP * ( I1 - st );
- shearfailure=tensilefailure=0;
- if( h > 0.0f )
- {
- shearfailure=1;
- /* !shear failure */
- alam = fs / (bulkm*qphi*qpsi + rmu);
- if(alam < 0.0)alam=palam*(previousfs/(previousfs-fs));
- if(isnan(alam) || isinf(alam))printf("alam = %g \n",alam);
- if(isnan(J2) || isinf(J2))printf("alam = %g \n",alam);
- dep0 = alam*(0.5*stressd0/J2+qpsi/3.);
- dep1 = alam*(0.5*stressd1/J2+qpsi/3.);
- dep2 = alam*(0.5*stressd2/J2+qpsi/3.);
- dep01 = alam*0.5*((*stress)[0][1])/J2;
- dep02 = alam*0.5*((*stress)[0][2])/J2;
- dep12 = alam*0.5*((*stress)[1][2])/J2;
- }
- else
- {
- tensilefailure=1;
- /* tensile failure */
- alam = ft / bulkm;
- if(alam < 0.0)alam=palam*(previousft/(previousfs-ft));
- dep0 =alam/3.;
- dep1 =alam/3.;
- dep2 =alam/3.;
- }
- if(mite==1) break;
- /* newly added stuff */
+ a = dmatrix(1,3,1,3);
+ v = dmatrix(1,3,1,3);
+ d = dvector(1,3);
- depm = (dep0 + dep1 + dep2) / 3.0f;
+ 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];
- trace_strain = (*strain)[0][0] + (*strain)[1][1] + (*strain)[2][2];
- /* Deviatoric Stresses and Strains */
- straind0 = (*strain)[0][0] - dep0 - (trace_strain) / 3.0f+depm;
- straind1 = (*strain)[1][1] - dep1 - (trace_strain) / 3.0f+depm;
- straind2 = (*strain)[2][2] - dep2 - (trace_strain) / 3.0f+depm;
+ jacobi(a,d,v);
- trace_stress= Stress0[0][0]+Stress0[1][1]+Stress0[2][2];
+ d[1] *= -1.0f;
+ d[2] *= -1.0f;
+ d[3] *= -1.0f;
- stressd0 = Stress0[0][0] - trace_stress/3.;
- stressd1 = Stress0[1][1] - trace_stress/3.;
- stressd2 = Stress0[2][2] - trace_stress/3.;
+ eigsrt(d,v);
- /* viscous relaxation with reduced plastic strain to account for plasticity */
- stressd0 = (stressd0 * vic1 + 2.0f * rmu * straind0) * vic2 ;
- stressd1 = (stressd1 * vic1 + 2.0f * rmu * straind1) * vic2 ;
- stressd2 = (stressd2 * vic1 + 2.0f * rmu * straind2) * vic2 ;
+ d[1] *= -1.0f;
+ d[2] *= -1.0f;
+ d[3] *= -1.0f;
- (*stress)[0][1] =(Stress0[0][1] * vic1 + 2.0f * rmu * ((*strain)[0][1]-dep01)) * vic2;
- (*stress)[0][2] =(Stress0[0][2] * vic1 + 2.0f * rmu * ((*strain)[0][2]-dep02)) * vic2;
- (*stress)[1][2] =(Stress0[1][2] * vic1 + 2.0f * rmu * ((*strain)[1][2]-dep12)) * vic2;
-
- VolumicStress = trace_stress / 3.0f + bulkm * (trace_strain - 3.*depm);
-
- (*stress)[0][0] = stressd0 + VolumicStress;
- (*stress)[1][1] = stressd1 + VolumicStress;
- (*stress)[2][2] = stressd2 + VolumicStress;
- I1 = ((*stress)[0][0]+(*stress)[1][1]+(*stress)[2][2])/3.;
- stressd0=(*stress)[0][0]-I1;
- stressd1=(*stress)[1][1]-I1;
- stressd2=(*stress)[2][2]-I1;
- J2 = sqrt(0.5*(stressd0*stressd0+stressd1*stressd1+stressd2*stressd2)
- +(*stress)[0][1]*(*stress)[0][1]+(*stress)[0][2]*(*stress)[0][2]+(*stress)[1][2]*(*stress)[1][2]);
- previousfs=fs;
- previousft=ft;
- palam=alam;
- fs = J2 + I1 * qphi - kphi;
- ft = I1 - st;
- if(fs < tolerf && ft < tolerf) break;
-
- }
- plasticStrain=viscoplasticElement->plasticStrain[tetra_I]+sqrt( 0.5f * ((dep0-depm) * (dep0-depm) + (dep1-depm) * (dep1-depm) + (dep2-depm) * (dep2-depm)) + dep01*dep01 + dep02*dep02 + dep12*dep12);
-
-
- if(isnan((*stress)[0][0])||isinf((*stress)[0][0]))printf("stress[0][0]=%g\n",(*stress)[0][0]);
- if(isnan((*stress)[1][1])||isinf((*stress)[1][1]))printf("stress[1][1]=%g\n",(*stress)[1][1]);
- if(isnan((*stress)[2][2])||isinf((*stress)[2][2]))printf("stress[2][2]=%g\n",(*stress)[2][2]);
- if(isnan((*stress)[0][1])||isinf((*stress)[0][1]))printf("stress[0][1]=%g\n",(*stress)[0][1]);
- if(isnan((*stress)[0][2])||isinf((*stress)[0][2]))printf("stress[0][2]=%g\n",(*stress)[0][2]);
- if(isnan((*stress)[1][2])||isinf((*stress)[1][2]))printf("stress[1][2]=%g\n",(*stress)[1][2]);
- }
- else
- {
- /* ! no failure - just elastic increment */
-
- dep0 = dep1 = dep2 = dep01 = dep02 = 0.0f;
- }
-
- if(ind)
- {
- (*plasticstrainTensor)[0][0]+=dep0;
- (*plasticstrainTensor)[1][1]+=dep1;
- (*plasticstrainTensor)[2][2]+=dep2;
- (*plasticstrainTensor)[0][1]+=dep01;
- (*plasticstrainTensor)[0][2]+=dep02;
- (*plasticstrainTensor)[1][2]+=dep12;
- /* Second invariant of plastic strain */
- depm = ( (*plasticstrainTensor)[0][0]+(*plasticstrainTensor)[1][1]+(*plasticstrainTensor)[2][2] ) / 3.0f;
- viscoplasticElement->plasticStrain[tetra_I]
- = sqrt(0.5f*(((*plasticstrainTensor)[0][0]-depm) * ((*plasticstrainTensor)[0][0]-depm) + ((*plasticstrainTensor)[1][1]-depm) * ((*plasticstrainTensor)[1][1]-depm) + ((*plasticstrainTensor)[2][2]-depm) * ((*plasticstrainTensor)[2][2]-depm)) + (*plasticstrainTensor)[0][1]*(*plasticstrainTensor)[0][1] + (*plasticstrainTensor)[0][2]*(*plasticstrainTensor)[0][2] + (*plasticstrainTensor)[1][2]*(*plasticstrainTensor)[1][2]);
- }
- }
- 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;
+ 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);
- int principal_stresses(StressTensor* stress, double sp[3], double cn[3][3])
- {
+ return(1);
+}
- double **a,**v,*d;
- int i,j;
+int jacobi(double** a, double* d, double** v)
+{
- a = dmatrix(1,3,1,3);
- v = dmatrix(1,3,1,3);
- d = dvector(1,3);
+ int nrot = 0;
+ const unsigned int nmax = 100, n = 3;
+ double b[nmax], z[nmax], tresh,sm,g,h,t,theta,c,s,tau;
- 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];
+ int i,j,ip,iq;
- jacobi(a,d,v);
+ for(ip=1;ip<=n;ip++) {
+ for(iq=1;iq<=n;iq++) v[ip][iq] = 0.0f;
+ v[ip][ip] = 1.0f;
+ }
- 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);
+ for(ip=1;ip<=n;ip++) {
+ b[ip] = d[ip] = a[ip][ip];
+ z[ip] = 0.0f;
}
- 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(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);
- for(ip=1;ip<=n;ip++) {
- b[ip] = d[ip] = a[ip][ip];
- z[ip] = 0.0f;
+ if(i < 4) {
+ tresh = 0.2f * sm / ( n*n );
}
+ else {
+ tresh = 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-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;
+ 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;
+int eigsrt(double* d, double** v)
+{
- for(i=1;i<n;i++) {
- k = i;
- p = d[i];
+ const unsigned int n = 3;
+ int i,j,k;
+ double p;
- for(j=i+1;j<=n;j++) {
- if(d[j] >= p) {
- k = j;
- p = d[j];
- }
+ 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;
- }
+ }
+ 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);
}
+ 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;
+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;
- }
+ /* 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 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);
- }
+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;
- }
+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;
- }
+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_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));
- }
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ConstructExtensions.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -33,77 +33,57 @@
#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 "ViscoPlastic.h"
#include <assert.h>
#include <limits.h>
#ifndef PATH_MAX
#define PATH_MAX 1024
#endif
-void _SnacVelocity_VariableCondition( Index index, Variable_Index var_I, void* _context, void* result ){
- Snac_Context* context = (Snac_Context*)_context;
-
- double* vy = (double*)result;
- if(context->timeStep > 0.02 * context->maxTimeSteps)
- (*vy)=0.0;
- else
- (*vy) = 0.005*erfc(context->currentTime*4.0);
-
-}
-
-
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 );
+ 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 };
+ /* 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
+ // #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
- ConditionFunction_Register_Add(
- context->condFunc_Register,
- ConditionFunction_New( _SnacVelocity_VariableCondition, "variablevelBC" ) );
-
/* Prepare the dump file */
sprintf( tmpBuf, "%s/plStrain.%u", context->outputPath, context->rank );
if( (contextExt->plStrainOut = fopen( tmpBuf, "w+" )) == NULL ) {
@@ -115,6 +95,11 @@
assert( contextExt->plStrainCheckpoint /* failed to open file for writing */ );
abort();
}
+ sprintf( tmpBuf, "%s/avgPlStrainCP.%u", context->outputPath, context->rank );
+ if( (contextExt->avgPlStrainCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->avgPlStrainCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
sprintf( tmpBuf, "%s/viscosity.%u", context->outputPath, context->rank );
if( (contextExt->viscOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->viscOut /* failed to open file for writing */ );
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Context.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -46,6 +46,7 @@
FILE* plStrainOut;
FILE* viscOut;
FILE* plStrainCheckpoint;
+ FILE* avgPlStrainCheckpoint;
FILE* viscCheckpoint;
FILE* plstrainTensorOut;
};
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -0,0 +1,62 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: Destory.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 "ViscoPlastic.h"
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+
+
+void _SnacViscoPlastic_Destroy( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+
+ /* Close any of the following if present. */
+ if( contextExt->plStrainOut )
+ fclose( contextExt->plStrainOut );
+ if( contextExt->plStrainCheckpoint )
+ fclose( contextExt->plStrainCheckpoint );
+ if( contextExt->avgPlStrainCheckpoint )
+ fclose( contextExt->avgPlStrainCheckpoint );
+ if( contextExt->viscOut )
+ fclose( contextExt->viscOut );
+ if( contextExt->viscCheckpoint )
+ fclose( contextExt->viscCheckpoint );
+
+}
+
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Destroy.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -0,0 +1,37 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: Destroy.h 3104 2005-07-14 22:16:41Z EunseoChoi $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef _SnacViscoPlastic_Destroy_
+#define _SnacViscoPlastic_Destroy_
+
+ void _SnacViscoPlastic_Destroy( void* _context );
+
+#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Element.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Element.h 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Element.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -46,8 +46,6 @@
Strain plasticStrain[Tetrahedra_Count];
double viscosity[Tetrahedra_Count];
Strain aps;
- StrainTensor plasticstrainTensor[Tetrahedra_Count];
- StrainTensor creepstrainTensor[Tetrahedra_Count];
};
/* Print the contents of an Element */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/InitialConditions.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -33,7 +33,6 @@
#include <StGermain/StGermain.h>
#include <StGermain/FD/FD.h>
#include "Snac/Snac.h"
-#include "units.h"
#include "types.h"
#include "Element.h"
#include "InitialConditions.h"
@@ -85,10 +84,13 @@
}
if( context->restartTimestep > 0 ) {
FILE* plStrainIn;
+ FILE* avgPlStrainIn;
char path[PATH_MAX];
sprintf(path, "%s/snac.plStrain.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
Journal_Firewall( (plStrainIn = fopen(path,"r")) != NULL, "Can't find %s", path );
+ sprintf(path, "%s/snac.avgPlStrain.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
+ Journal_Firewall( (avgPlStrainIn = fopen(path,"r")) != NULL, "Can't find %s", path );
/* read in restart file to reconstruct the previous plastic strain.*/
for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
@@ -98,29 +100,30 @@
const Snac_Material* material = &context->materialProperty[element->material_I];
vis_min = context->materialProperty[element->material_I].vis_min;
+ /* To-Do: Read in viscosity and assign it to tets and elements. */
if( material->yieldcriterion == mohrcoulomb ) {
Tetrahedra_Index tetra_I;
- double depls = 0.0f;
- double totalVolume = 0.0f;
+ double avgPlStrain = 0.0f;
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- double tetraPlStrain;
+ double tetraPlStrain = 0.0f;
/* 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
+ }/* for tets */
/* volume-averaged accumulated plastic strain, aps */
- viscoplasticElement->aps = depls/totalVolume;
- }//if(mohrcoulomb)
- }//for elements
+ fscanf( avgPlStrainIn, "%le", &avgPlStrain );
+ viscoplasticElement->aps = avgPlStrain;
+ }/* if(mohrcoulomb) */
+ }/* for elements */
if( plStrainIn )
fclose( plStrainIn );
- }// if restarting.
+ if( avgPlStrainIn )
+ fclose( avgPlStrainIn );
+ }/* if restarting.*/
else {
/* Set the plastic element initial conditions */
for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Makefile.def 2011-11-14 15:57:18 UTC (rev 19193)
@@ -34,19 +34,24 @@
def_inc = Snac/ViscoPlastic
def_srcs = \
+ Node.c \
Element.c \
+ Mesh.c \
Context.c \
ConstructExtensions.c \
Constitutive.c \
Output.c \
Register.c \
InitialConditions.c \
+ Destroy.c \
Remesh.c \
def_hdrs = \
types.h \
- units.h \
+ units.h \
+ Node.h \
Element.h \
+ Mesh.h \
Context.h \
ConstructExtensions.h \
Constitutive.h \
@@ -54,5 +59,6 @@
Register.h \
InitialConditions.h \
Remesh.h \
+ Destroy.h \
ViscoPlastic.h
Added: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -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/viscoplastic_BI/Mesh.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Mesh.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -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/viscoplastic_BI/Node.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Node.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Node.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -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/viscoplastic_BI/Node.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Node.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Node.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -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__ */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Output.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -97,13 +97,19 @@
context->mesh->elementExtensionMgr,
element,
SnacViscoPlastic_ElementHandle );
+ float avgPlasticStrain = elementExt->aps;
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 );
+
+
+ fwrite( &avgPlasticStrain, sizeof(float), 1, contextExt->avgPlStrainCheckpoint );
}
+ fflush( contextExt->avgPlStrainCheckpoint );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -33,24 +33,28 @@
#include <StGermain/StGermain.h>
#include <StGermain/FD/FD.h>
#include "Snac/Snac.h"
-#include "types.h"
+#include "ConstructExtensions.h"
+#include "Constitutive.h"
+#include "Context.h"
+#include "Destroy.h"
+#include "InitialConditions.h"
#include "Element.h"
-#include "Context.h"
-#include "Constitutive.h"
-#include "ConstructExtensions.h"
+#include "Mesh.h"
+#include "Node.h"
#include "Output.h"
-#include "InitialConditions.h"
+#include "Register.h"
#include "Remesh.h"
-#include "Register.h"
+#include "types.h"
#include <stdio.h>
/* Textual name of this class */
const Type SnacViscoPlastic_Type = "SnacViscoPlastic";
+ExtensionInfo_Index SnacViscoPlastic_NodeHandle;
ExtensionInfo_Index SnacViscoPlastic_ElementHandle;
+ExtensionInfo_Index SnacViscoPlastic_MeshHandle;
ExtensionInfo_Index SnacViscoPlastic_ContextHandle;
-
Index _SnacViscoPlastic_Register( PluginsManager* pluginsMgr ) {
return PluginsManager_Submit( pluginsMgr,
SnacViscoPlastic_Type,
@@ -87,12 +91,16 @@
#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( "\telement extension handle: %u\n", SnacViscoPlastic_ElementHandle );
+ 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 */
@@ -127,15 +135,20 @@
"SnacViscoPlastic_WriteViscosity",
_SnacViscoPlastic_WriteViscosity,
SnacViscoPlastic_Type );
+ EntryPoint_Append( /* perform clean-ups like closing output files. */
+ Context_GetEntryPoint( context, AbstractContext_EP_DestroyExtensions ),
+ "SnacViscoPlastic_Destory",
+ _SnacViscoPlastic_Destroy,
+ SnacViscoPlastic_Type );
/* Add extensions to the interpolate element entry point, but it will only exist if the remesher is loaded. */
interpolateElementEP = Context_GetEntryPoint( context, "SnacRemesher_EP_InterpolateElement" );
if( interpolateElementEP ) {
EntryPoint_Append(
- interpolateElementEP,
- SnacViscoPlastic_Type,
- _SnacViscoPlastic_InterpolateElement,
- SnacViscoPlastic_Type );
+ interpolateElementEP,
+ SnacViscoPlastic_Type,
+ _SnacViscoPlastic_InterpolateElement,
+ SnacViscoPlastic_Type );
}
/* Construct. */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Register.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -44,8 +44,12 @@
/* Textual name of this class */
extern const Type SnacViscoPlastic_Type;
+ extern ExtensionInfo_Index SnacViscoPlastic_NodeHandle;
+
extern ExtensionInfo_Index SnacViscoPlastic_ElementHandle;
-
+
+ extern ExtensionInfo_Index SnacViscoPlastic_MeshHandle;
+
extern ExtensionInfo_Index SnacViscoPlastic_ContextHandle;
Index _SnacViscoPlastic_Register( PluginsManager* pluginsMgr );
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/Remesh.c 2011-11-14 15:57:18 UTC (rev 19193)
@@ -87,13 +87,14 @@
elementExt->plasticStrain[dstTetInd] = 0.0;
for(coef_I=0;coef_I<4;coef_I++) {
- /* The actual src elements are the apexes of the tet in the old barycenter grid. */
- Snac_Element* srcElt = (Snac_Element*)ExtensionManager_At( context->mesh->elementExtensionMgr,
- dstElements,
- meshExt->orderedToDomain[eltdI[TetraToNode[srcTetInd][coef_I]]] );
- SnacViscoPlastic_Element* srcEltExt = ExtensionManager_Get( context->mesh->elementExtensionMgr,
- srcElt,
- SnacViscoPlastic_ElementHandle );
+ /* The actual src elements are the four apexes of a tet (srcTetInd) in the old barycenter grid. */
+ Snac_Element* srcElt = Snac_Element_At( context, eltdI[TetraToNode[srcTetInd][coef_I]]);
+ SnacViscoPlastic_Element* srcEltExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ srcElt,
+ SnacViscoPlastic_ElementHandle );
+ /* Weights are associated only with destination element but not on the tet level.
+ So, "dstTetInd" is used in both source and destination terms. */
elementExt->plasticStrain[dstTetInd] += meshExt->barcord[dstEltInd].L[coef_I]*srcEltExt->plasticStrain[dstTetInd];
}
}
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/ViscoPlastic.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -41,16 +41,19 @@
#ifndef __SnacViscoPlastic_h__
#define __SnacViscoPlastic_h__
-#include "units.h"
-#include "types.h"
+#include "ViscoPlastic.h"
+#include "Destroy.h"
+#include "Context.h"
+#include "Constitutive.h"
+#include "ConstructExtensions.h"
#include "Element.h"
-#include "Context.h"
#include "InitialConditions.h"
-#include "Constitutive.h"
+#include "Mesh.h"
+#include "Node.h"
#include "Output.h"
-#include "ConstructExtensions.h"
-#include "ViscoPlastic.h"
+#include "Register.h"
#include "Remesh.h"
-#include "Register.h"
+#include "types.h"
+#include "units.h"
#endif /* __SnacViscoPlastic_h__ */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/types.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/types.h 2011-11-14 05:06:30 UTC (rev 19192)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic_BI/types.h 2011-11-14 15:57:18 UTC (rev 19193)
@@ -44,7 +44,9 @@
#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__ */
More information about the CIG-COMMITS
mailing list