[cig-commits] r14159 - in long/3D/SNAC/trunk/Snac/plugins: . hillSlope hillSlope/archive hillSlope/archive/tmp hillSlope/archive/unused
cstark at geodynamics.org
cstark at geodynamics.org
Thu Feb 26 09:49:08 PST 2009
Author: cstark
Date: 2009-02-26 09:48:54 -0800 (Thu, 26 Feb 2009)
New Revision: 14159
Added:
long/3D/SNAC/trunk/Snac/plugins/hillSlope/
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/HillSlope.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Make.mm
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Makefile.def
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Context.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c.preshellsort
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Register.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts2.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/cuts3.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/hillslope_demo1.xml
long/3D/SNAC/trunk/Snac/plugins/hillSlope/make.log
long/3D/SNAC/trunk/Snac/plugins/hillSlope/makefile
long/3D/SNAC/trunk/Snac/plugins/hillSlope/types.h
Log:
Test commit of hillSlope plugin.
HillSlope does/provides the following:
- generates a deformed mesh in the shape of a flat-ramp-flat hillslope cross-section (slopes in x direction)
- tracks establishment of elastic equilibrium
- creates weak points ("plastic seeds") within subdomain of mesh
- creates a trigger weak point at chosen location
- turns on these weak points once elastic eqm is achieved
- provides flags to control/report above behavior
Since this is a test commit, will wait to provide more details.
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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: Build.c $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Register.h"
+#include "Build.h"
+
+void _SnacHillSlope_Build( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+/* SnacHillSlope_Context* contextExt = ExtensionManager_Get( */
+/* context->extensionMgr, */
+/* context, */
+/* SnacHillSlope_ContextHandle ); */
+
+ Journal_Printf( context->debug, "In: %s\n", __func__ );
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Build.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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.
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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: Build.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_Build_h__
+#define __SnacHillSlope_Build_h__
+
+ void _SnacHillSlope_Build( void* _context, void* data );
+
+#endif /* __SnacHillSlope_Build_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,184 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "InitialConditions.h"
+#include "ConstructExtensions.h"
+#include <assert.h>
+#include <limits.h>
+
+//#define DEBUG
+
+void _SnacHillSlope_ConstructExtensions( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacHillSlope_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+/* Dictionary* hillSlopeBCsDict; */
+/* char tmpBuf[PATH_MAX]; */
+
+ /* Because hillSlope is not an array by itself, we must the "complex" constructor for Variable... the info needs to be
+ * wrapped this generic way... */
+/* Index hillSlopeOffsetCount = 1; */
+/* SizeT hillSlopeOffsets[] = { (SizeT)((char*)&tmpNodeExt->hillSlope - (char*)&tmpNode) }; */
+/* Variable_DataType hillSlopeDataTypes[] = { Variable_DataType_Double }; */
+/* Index hillSlopeDataTypeCounts[] = { 1 }; */
+ #ifdef DEBUG
+ printf( "In %s()\n", __func__ );
+ #endif
+
+
+ /* Create the StGermain variable hillSlope, which is stored on a node extension */
+/* Variable_New( */
+/* "hillSlope", */
+/* hillSlopeOffsetCount, */
+/* hillSlopeOffsets, */
+/* hillSlopeDataTypes, */
+/* hillSlopeDataTypeCounts, */
+/* 0, */
+/* &ExtensionManager_GetFinalSize( context->mesh->nodeExtensionMgr ), */
+/* &context->mesh->layout->decomp->nodeDomainCount, */
+/* (void**)&context->mesh->node, */
+/* context->variable_Register ); */
+
+
+
+ /* HillSlope variables */
+ contextExt->slopeAngle = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "slopeAngle",
+ Dictionary_Entry_Value_FromDouble( 0.0f ) ) );
+
+ contextExt->rngSeed = Dictionary_Entry_Value_AsUnsignedInt(
+ Dictionary_GetDefault( context->dictionary, "rngSeed",
+ Dictionary_Entry_Value_FromUnsignedInt( 1 ) ) );
+
+ contextExt->fractionWeakPoints = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "fractionWeakPoints",
+ Dictionary_Entry_Value_FromDouble( 0.02f ) ) );
+
+ contextExt->resolveDepth = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "resolveDepth",
+ Dictionary_Entry_Value_FromDouble( -1.0f ) ) );
+
+ contextExt->leftFlatFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "leftFlatFraction",
+ Dictionary_Entry_Value_FromDouble( 0.05f ) ) );
+ contextExt->rightFlatFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "rightFlatFraction",
+ Dictionary_Entry_Value_FromDouble( -1.0f ) ) );
+ contextExt->rampFlatSmoothFactor = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "rampFlatSmoothFactor",
+ Dictionary_Entry_Value_FromDouble( 3.0f ) ) );
+
+ contextExt->xSubDomainFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "xSubDomainFraction",
+ Dictionary_Entry_Value_FromDouble( 0.5f ) ) );
+ contextExt->ySubDomainFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "ySubDomainFraction",
+ Dictionary_Entry_Value_FromDouble( 0.4f ) ) );
+ contextExt->zSubDomainFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "zSubDomainFraction",
+ Dictionary_Entry_Value_FromDouble( 0.75f ) ) );
+
+ contextExt->xTriggerPointFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "xTriggerPointFraction",
+ Dictionary_Entry_Value_FromDouble( 0.5f ) ) );
+ contextExt->yTriggerPointFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "yTriggerPointFraction",
+ Dictionary_Entry_Value_FromDouble( 0.4f ) ) );
+ contextExt->zTriggerPointFraction = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "zTriggerPointFraction",
+ Dictionary_Entry_Value_FromDouble( 0.5f ) ) );
+
+ contextExt->weakPointCohesion = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "weakPointCohesion",
+ Dictionary_Entry_Value_FromDouble( 4.0e+05f ) ) );
+ contextExt->triggerPointCohesion = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "triggerPointCohesion",
+ Dictionary_Entry_Value_FromDouble( 1.0e+04f ) ) );
+
+ contextExt->trackLevel = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "trackLevel",
+ Dictionary_Entry_Value_FromDouble( 0.0f ) ) );
+
+ contextExt->startThreshold = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "startThreshold",
+ Dictionary_Entry_Value_FromDouble( 1.0f ) ) );
+ contextExt->stopThreshold = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "stopThreshold",
+ Dictionary_Entry_Value_FromDouble( 1.0f ) ) );
+
+ contextExt->startedTrackingFlag = Dictionary_Entry_Value_AsBool(
+ Dictionary_GetDefault( context->dictionary, "startedTrackingFlag",
+ Dictionary_Entry_Value_FromBool( 0 ) ) );
+ contextExt->elasticStabilizedFlag = Dictionary_Entry_Value_AsBool(
+ Dictionary_GetDefault( context->dictionary, "elasticStabilizedFlag",
+ Dictionary_Entry_Value_FromBool( 0 ) ) );
+ contextExt->solveElasticEqmOnlyFlag = Dictionary_Entry_Value_AsBool(
+ Dictionary_GetDefault( context->dictionary, "solveElasticEqmOnlyFlag",
+ Dictionary_Entry_Value_FromBool( 0 ) ) );
+ contextExt->seedingCompletedFlag = Dictionary_Entry_Value_AsBool(
+ Dictionary_GetDefault( context->dictionary, "seedingCompletedFlag",
+ Dictionary_Entry_Value_FromBool( 0 ) ) );
+
+ /*
+ * Open the output streams for mesh data
+ */
+/* sprintf( tmpBuf, "%s/shearStress.%u", context->outputPath, context->rank ); */
+/* if( (contextExt->shearStress = fopen( tmpBuf, "w+" )) == NULL ) { */
+/* assert( contextExt->shearStress ); */
+/* } */
+/* sprintf( tmpBuf, "%s/failurePotential.%u", context->outputPath, context->rank ); */
+/* if( (contextExt->failurePotential = fopen( tmpBuf, "w+" )) == NULL ) { */
+/* assert( contextExt->failurePotential); */
+/* } */
+
+
+ /* Build the hillSlope IC and BC managers */
+/* hillSlopeBCsDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( context->dictionary, "hillSlopeBCs" ) ); */
+/* contextExt->hillSlopeBCs = CompositeVC_New("tempBC", */
+/* context->variable_Register, */
+/* context->condFunc_Register, */
+/* hillSlopeBCsDict, */
+/* context->mesh ); */
+
+
+/* #ifdef DEBUG */
+/* fprintf( stderr, "In %s()\n", __func__ ); */
+/* #endif */
+
+/* Journal_Printf( context->debug, "slopeAngle: %g\n", contextExt->slopeAngle ); */
+
+}
Property changes on: long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
___________________________________________________________________
Name: svn:executable
+ *
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,48 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_ConstructExtensions_h__
+#define __SnacHillSlope_ConstructExtensions_h__
+
+ void _SnacHillSlope_ConstructExtensions( void* _context, void* data );
+
+#endif /* __SnacHillSlope_ConstructExtensions_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 SnacHillSlope_Context_Print( void* _context, Stream* stream ) {
+ // SnacHillSlope_Context* self = (SnacHillSlope_Context*)_context;
+
+ Journal_Printf( stream, "SnacHillSlope_Context:\n" );
+
+/* Print( self->slopeAngle, stream ); */
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,81 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_Context_h__
+#define __SnacHillSlope_Context_h__
+
+ /* Context Information */
+ struct _SnacHillSlope_Context {
+ /* For _SnacHillSlope_Top2BottomSweep condition function */
+ double slopeAngle;
+ unsigned int rngSeed;
+ double fractionWeakPoints;
+ double resolveDepth;
+ double leftFlatFraction;
+ double rightFlatFraction;
+ double rampFlatSmoothFactor;
+ double xSubDomainFraction;
+ double ySubDomainFraction;
+ double zSubDomainFraction;
+ double xTriggerPointFraction;
+ double yTriggerPointFraction;
+ double zTriggerPointFraction;
+ double weakPointCohesion;
+ double triggerPointCohesion;
+ double trackLevel;
+ double startThreshold;
+ double stopThreshold;
+ int startedTrackingFlag;
+ int elasticStabilizedFlag;
+ int solveElasticEqmOnlyFlag;
+ int seedingCompletedFlag;
+ /* Data output to file */
+/* FILE* shearStress; */
+/* FILE* failurePotential; */
+ };
+
+ /* Print the contents of the context extension */
+ void SnacHillSlope_Context_Print( void* _context, Stream* stream );
+
+#endif /* __SnacHillSlope_Context_h__ */
+
+
Property changes on: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
___________________________________________________________________
Name: svn:executable
+ *
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,406 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ **
+ ** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+ ** 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: CreateWeakPoints.c $
+ **
+ **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "CreateWeakPoints.h"
+#include "Snac/Plastic/Plastic.h"
+
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#ifndef TRUE
+#define TRUE 1
+#endif
+#ifndef FALSE
+#define FALSE 0
+#endif
+
+#ifndef PI
+ #ifndef M_PIl
+ #ifndef M_PI
+ #define PI 3.14159265358979323846
+ #else
+ #define PI M_PI
+ #endif
+ #else
+ #define PI M_PIl
+ #endif
+#endif
+#define min(a,b) ((a)<(b) ? (a):(b))
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+#define DEBUG
+
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * SnacHillSlope_CreateWeakPoints --
+ *
+ * tbd
+ *
+ * Returns:
+ * Nothing directly
+ *
+ * Side effects:
+ * Altered mesh with weak points induced by forced plastic deformation
+ *
+ *----------------------------------------------------------------------
+ */
+void SnacHillSlope_CreateWeakPoints( void* _context ) {
+ Snac_Context *context = (Snac_Context*)_context;
+ SnacHillSlope_Context *contextExt = ExtensionManager_Get(context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+ Mesh *mesh = context->mesh;
+ MeshLayout *layout = (MeshLayout*)mesh->layout;
+ HexaMD *decomp = (HexaMD*)layout->decomp;
+
+ double fractionWeakPoints = contextExt->fractionWeakPoints;
+
+ const int full_I_element_range=decomp->elementGlobal3DCounts[0];
+ const int full_J_element_range=decomp->elementGlobal3DCounts[1];
+ const int full_K_element_range=decomp->elementGlobal3DCounts[2];
+
+ Element_LocalIndex element_lI;
+ Element_GlobalIndex element_gI;
+ int index, index_I,index_J,index_K;
+ IJK ijk;
+
+ /* int restart = 0; */
+ /* Dictionary_Entry_Value *pluginsList, *plugin; */
+ /* #ifdef DEBUG */
+ /* int imax=0; */
+ /* #endif */
+ /* pluginsList = PluginsManager_GetPluginsList( context->dictionary ); */
+ /* if (pluginsList) { */
+ /* plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList); */
+ /* while ( plugin ) { */
+ /* if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ), */
+ /* "SnacRestart" ) ) { */
+ /* restart = 1; */
+ /* break; */
+ /* } */
+ /* plugin = plugin->next; */
+ /* } */
+ /* } */
+ /* if( restart ) */
+ /* return; */
+
+
+
+ /*
+ * Bail now if initial elastic equilibrium has not been reached
+ */
+ if(!contextExt->elasticStabilizedFlag || contextExt->seedingCompletedFlag) return;
+
+ // fprintf(stderr, "CWP\n");
+
+ /*
+ * Insert weak (low cohesion) points at N random elements by forcing plastic strain
+ */
+ if(fractionWeakPoints>=0.0 && fractionWeakPoints<=1.0) {
+ unsigned int rngSeed = contextExt->rngSeed;
+ double x_subdomainFraction = contextExt->xSubDomainFraction;
+ double y_subdomainFraction = contextExt->ySubDomainFraction;
+ double z_subdomainFraction = contextExt->zSubDomainFraction;
+ unsigned int numberSubDomainPoints,numberWeakPoints;
+ unsigned int *shuffleIndexListPtr, *randomNumberListPtr;
+ int subdomain_I_element_range,subdomain_J_element_range,subdomain_K_element_range;
+ /* Report HillSlope plugin variables picked up (?) from xml parameter file */
+ Journal_Printf( context->snacInfo, "\tRNG seed = %u\n", rngSeed );
+ Journal_Printf( context->snacInfo, "\tFraction of weak points = %g\n", fractionWeakPoints );
+ Journal_Printf( context->snacInfo, "\tx subdomain fraction = %g\n", x_subdomainFraction );
+ Journal_Printf( context->snacInfo, "\ty subdomain fraction = %g\n", y_subdomainFraction );
+ Journal_Printf( context->snacInfo, "\tz subdomain fraction = %g\n", z_subdomainFraction );
+
+ /*
+ * Define portion of mesh (cells) that need weak point "seeding"
+ * - lots of sloppy float/int casting back and forth here
+ */
+ subdomain_I_element_range = (int)(full_I_element_range*x_subdomainFraction);
+ subdomain_J_element_range = (int)(full_J_element_range*y_subdomainFraction);
+ subdomain_K_element_range = (int)(full_K_element_range*z_subdomainFraction);
+
+ Journal_Printf( context->snacInfo, "\tx subdomain element range = %d/%d\n", subdomain_I_element_range,full_I_element_range );
+ Journal_Printf( context->snacInfo, "\ty subdomain element range = %d/%d\n", subdomain_J_element_range,full_J_element_range );
+ Journal_Printf( context->snacInfo, "\tz subdomain element range = %d/%d\n\n", subdomain_K_element_range,full_K_element_range );
+
+ numberSubDomainPoints = subdomain_I_element_range*subdomain_J_element_range*subdomain_K_element_range;
+ numberWeakPoints = (unsigned int)((float)numberSubDomainPoints*fractionWeakPoints);
+ Journal_Printf(context->snacInfo, "Number of weak points = %u/%u\n", numberWeakPoints,numberSubDomainPoints);
+
+ /*
+ * Seed RNG (random number generator) using seed value from input xml file
+ */
+ srand( rngSeed );
+ shuffleIndexListPtr=(unsigned int *)malloc((size_t)(numberSubDomainPoints*sizeof(unsigned int)));
+ randomNumberListPtr=(unsigned int *)malloc((size_t)(numberSubDomainPoints*sizeof(unsigned int)));
+ /*
+ * Create (1) ordered sequence of element indices, (2) parallel list of random variates
+ */
+ index=0;
+ for(index_I = 0; index_I < subdomain_I_element_range; index_I++) {
+ for(index_J = 0; index_J < subdomain_J_element_range; index_J++) {
+ for(index_K = 0; index_K < subdomain_K_element_range; index_K++) {
+ /*
+ * Work out the element index from the i,j,k and the required location of the subdomain
+ */
+ element_gI = index_I+(full_I_element_range-subdomain_I_element_range)/2
+ + full_I_element_range*(full_J_element_range-1-index_J)
+ + full_I_element_range*full_J_element_range*(index_K+(full_K_element_range-subdomain_K_element_range)/2);
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+ shuffleIndexListPtr[index]=element_lI;
+ randomNumberListPtr[index]=rand();
+ index++;
+ if(index-1>numberSubDomainPoints) break;
+ }
+ if(index>0 && index-1>numberSubDomainPoints) break;
+ }
+ if(index>0 && index-1>numberSubDomainPoints) {fprintf(stderr, "We miscalculated the subdomain index\n"); break;}
+ }
+ /*
+ * Sort the list of sub-domain local (ie global set pared down for this CPU)
+ * indices according to the random numbers
+ * - i.e., randomize the index list
+ */
+ ShellSort( randomNumberListPtr, shuffleIndexListPtr, numberSubDomainPoints );
+#ifdef DEBUG
+ for(index = 0; index < numberSubDomainPoints; index++) {
+ Journal_Printf( context->snacInfo, "\t Subdomain point #%d = %d, RNG var=%d\n",
+ index, shuffleIndexListPtr[index],randomNumberListPtr[index] );
+ }
+#endif
+ /*
+ * Now the subdomain index list is shuffled and the random number list is ordered
+ */
+
+
+ /*
+ * Loop over the shuffled subdomain index list BUT only up to the number of weak points within it
+ * - i.e., subset the subdomain and allocate weak points to only a random fraction of them
+ */
+ for(index = 0; index < numberWeakPoints; index++) {
+ element_lI = shuffleIndexListPtr[index];
+
+#ifdef DEBUG
+ Journal_Printf(context->snacInfo, "Shuffled weak point # %d -> %d \n", index, element_lI);
+#endif
+
+ /*
+ * At each point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+ Snac_Element *element;
+ Snac_Material *material;
+ SnacPlastic_Element *plasticElement;
+ Tetrahedra_Index tetra_I;
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get( mesh->elementExtensionMgr, element, SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = PlasticStrainFromCohesion(material,(double)contextExt->weakPointCohesion);
+#ifdef DEBUG
+ if(tetra_I==0) Journal_Printf( context->snacInfo,"timeStep=%d et=%d %d plasticE: %e -> %e\n",context->timeStep, element_lI, tetra_I, material->plstrain[1], plasticElement->plasticStrain[tetra_I] );
+#endif
+
+ }
+ } // End if
+
+ } // End for
+ /*
+ * Free up allocated memory for shuffling arrays
+ */
+ free(shuffleIndexListPtr);
+ free(randomNumberListPtr);
+ } // End if (check fraction of weak points is in valid range [0,1]
+
+
+ fprintf(stderr, "CWP trigger pt\n");
+
+
+ /*
+ * Trigger point assignment
+ */
+ if(contextExt->triggerPointCohesion>=0.0) {
+ ijk[0]=(int)(full_I_element_range*contextExt->xTriggerPointFraction);
+ ijk[1]=(int)(full_J_element_range*contextExt->yTriggerPointFraction);
+ ijk[1]= full_J_element_range-1-ijk[1];
+ ijk[2]=(int)(full_K_element_range*contextExt->zTriggerPointFraction);
+ /*
+ * Calculate "global" (across all processors) element index and then find local equivalent (this CPU)
+ */
+ element_gI = ijk[0] + full_I_element_range*ijk[1] + full_I_element_range*full_J_element_range*ijk[2];
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+#ifdef DEBUG
+ Journal_Printf(context->snacInfo, "Trigger point: ijk = %d,%d,%d -> %d -> %d <? %d\n", ijk[0], ijk[1], ijk[2], element_gI, element_lI, mesh->elementLocalCount);
+#endif
+ /*
+ * At trigger point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+ Snac_Element *element;
+ Snac_Material *material;
+ SnacPlastic_Element *plasticElement;
+ Tetrahedra_Index tetra_I;
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get( mesh->elementExtensionMgr,element,SnacPlastic_ElementHandle );
+ /* Force each of 5*2 tetrahedra to have extra plastic strain to impose lower cohesion indirectly */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = PlasticStrainFromCohesion(material,(double)contextExt->triggerPointCohesion);
+#ifdef DEBUG
+ if(tetra_I==0) Journal_Printf(context->snacInfo, "timeStep=%d ijkt=%d %d %d %d plasticE=%e\n", context->timeStep,ijk[0],ijk[1],ijk[2], tetra_I, plasticElement->plasticStrain[tetra_I] );
+#endif
+ }
+ }
+ } // End if
+
+ /*
+ * Flag that seeing has been done and don't do any more
+ */
+ contextExt->seedingCompletedFlag = TRUE;
+}
+
+
+
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * PlasticStrainFromCohesion --
+ *
+ * Inverts plastic strain from known cohesion
+ * using material def from Plastic plugin
+ *
+ * Takes:
+ * material (ptr to struct)
+ * cohesion (double)
+ *
+ * Returns:
+ * plastic strain (double)
+ *
+ * Side effects:
+ * Modifies input vectors
+ *
+ *----------------------------------------------------------------------
+ */
+
+double
+PlasticStrainFromCohesion(Snac_Material *material, double cohesion)
+{
+ int i;
+ double plasticStrain=material->plstrain[0];
+
+ // fprintf(stderr, "In pl\n");
+ for( i = 0; i < material->nsegments; i++ ) {
+ const double pl1 = material->plstrain[i];
+ const double pl2 = material->plstrain[i+1];
+ const double coh1 = material->cohesion[i];
+ const double coh2 = material->cohesion[i+1];
+ plasticStrain = pl1+(pl2-pl1)*( (cohesion-coh1)/(coh2-coh1) );
+ if( plasticStrain >= pl1 && plasticStrain <= pl2 ) {
+#ifdef DEBUG
+ fprintf(stderr, "Returning plastic strain = %g from cohesion = %g (%g, %g)\n", plasticStrain, cohesion, coh1,coh2);
+#endif
+ return plasticStrain;
+ }
+ }
+#ifdef DEBUG
+ fprintf(stderr, "Bailing... default plastic strain = %g from cohesion = %g \n", plasticStrain, cohesion);
+#endif
+ return material->plstrain[material->nsegments];
+}
+
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * ShellSort --
+ *
+ * Shell sort vector (unsigned int) and coevally shuffle input
+ * (index, unsigned int) vector
+ *
+ * Returns:
+ * Nothing
+ *
+ * Side effects:
+ * Modifies input vectors
+ *
+ *----------------------------------------------------------------------
+ */
+
+void
+ShellSort(vecPtr,idxPtr,n)
+ unsigned int vecPtr[]; /* Input unsorted data vector (e.g. random numbers) */ /* was double */
+ unsigned int idxPtr[]; /* Input index vector, to be sorted according to vecPtr */
+ unsigned long n; /* Data length */
+{
+ unsigned long incr;
+ incr=1;
+ do {
+ incr *= 3;
+ incr++;
+ } while (incr <= n);
+ do {
+ unsigned long i,j;
+ incr /= 3;
+ for (i=incr;i<n;i++) {
+ unsigned int v; /* was double */
+ unsigned int w;
+ v=vecPtr[i];
+ w=idxPtr[i];
+ j=i;
+ while (vecPtr[j-incr] > v) {
+ vecPtr[j]=vecPtr[j-incr];
+ idxPtr[j]=idxPtr[j-incr];
+ j -= incr;
+ if (j < incr) {
+ break;
+ }
+ }
+ vecPtr[j]=v;
+ idxPtr[j]=w;
+ }
+ } while (incr > 1);
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,44 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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
+**
+** Role:
+**
+** Assumptions:
+** None as yet.
+**
+** Comments:
+** None as yet.
+**
+** $Id: CreateWeakPoints.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Snac_CreateWeakPoints_h__
+#define __Snac_CreateWeakPoints_h__
+
+void SnacHillSlope_CreateWeakPoints( void* context );
+double PlasticStrainFromCohesion(Snac_Material *material, double cohesion);
+void ShellSort( unsigned int* vecPtr, unsigned int* idxPtr, unsigned long n );
+
+#endif /* __Snac_CreateWeakPoints_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,48 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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: DeleteExtensions.c $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Register.h"
+#include "DeleteExtensions.h"
+
+void _SnacHillSlope_DeleteExtensions( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacHillSlope_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+
+ Journal_Printf( context->debug, "In: %s\n", __func__ );
+/* fclose( contextExt->shearStress ); */
+/* fclose( contextExt->failurePotential ); */
+
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/DeleteExtensions.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,48 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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: DeleteExtensions.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_DeleteExtensions_h__
+#define __SnacHillSlope_DeleteExtensions_h__
+
+ void _SnacHillSlope_DeleteExtensions( void* _context, void* data );
+
+#endif /* __SnacHillSlope_DeleteExtensions_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/HillSlope.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/HillSlope.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/HillSlope.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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: HillSlope.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_h__
+#define __SnacHillSlope_h__
+
+ #include "types.h"
+ #include "Context.h"
+ #include "Register.h"
+ #include "ConstructExtensions.h"
+ #include "Build.h"
+ #include "InitialConditions.h"
+ #include "DeleteExtensions.h"
+
+#endif /* __SnacHillSlope_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,241 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "InitialConditions.h"
+#include "Snac/Plastic/Plastic.h"
+
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#ifndef PI
+ #ifndef M_PIl
+ #ifndef M_PI
+ #define PI 3.14159265358979323846
+ #else
+ #define PI M_PI
+ #endif
+ #else
+ #define PI M_PIl
+ #endif
+#endif
+#define min(a,b) ((a)<(b) ? (a):(b))
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+
+//#define DEBUG
+
+void _SnacHillSlope_InitialConditions( void* _context, void* data ) {
+ Snac_Context *context = (Snac_Context*)_context;
+ SnacHillSlope_Context *contextExt = ExtensionManager_Get(context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+ Mesh *mesh = context->mesh;
+ MeshLayout *layout = (MeshLayout*)mesh->layout;
+ HexaMD *decomp = (HexaMD*)layout->decomp;
+ BlockGeometry *geometry = (BlockGeometry*)layout->elementLayout->geometry;
+ Node_GlobalIndex node_g, node_gI, node_gJ;
+
+ const int full_I_node_range=decomp->nodeGlobal3DCounts[0];
+ const int full_J_node_range=decomp->nodeGlobal3DCounts[1];
+
+ double new_x[full_I_node_range], new_y[full_J_node_range];
+ double reg_dx, reg_dy, new_xMax, new_xMin, new_xHalf, midOffset, old_yMax, new_yMax, new_yMin;
+
+ double slopeAngle = contextExt->slopeAngle * M_PI/180.0;
+ double resolveHeight = geometry->max[1]-contextExt->resolveDepth;
+
+ double xLeftFlatRamp = (contextExt->leftFlatFraction<0.0 ? 0.0 :
+ (contextExt->leftFlatFraction>1.0 ? 1.0 : contextExt->leftFlatFraction));
+ double xRightRampFlat = (contextExt->rightFlatFraction<0.0 ?
+ 1.0-contextExt->leftFlatFraction : (contextExt->rightFlatFraction>1.0 ?
+ 1.0 : contextExt->rightFlatFraction));
+ const double smoothFactor = contextExt->rampFlatSmoothFactor;
+
+
+ /* int restart = 0; */
+ /* Dictionary_Entry_Value *pluginsList, *plugin; */
+ /* #ifdef DEBUG */
+ /* int imax=0; */
+ /* #endif */
+ /* pluginsList = PluginsManager_GetPluginsList( context->dictionary ); */
+ /* if (pluginsList) { */
+ /* plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList); */
+ /* while ( plugin ) { */
+ /* if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ), */
+ /* "SnacRestart" ) ) { */
+ /* restart = 1; */
+ /* break; */
+ /* } */
+ /* plugin = plugin->next; */
+ /* } */
+ /* } */
+ /* if( restart ) */
+ /* return; */
+
+
+#ifdef DEBUG
+ printf( "In: %s\n", __func__ );
+#endif
+ fprintf(stderr, "Slope angle = %g degrees\n", slopeAngle/(M_PI/180.0));
+
+ /* Report HillSlope plugin variables picked up (?) from xml parameter file */
+ Journal_Printf( context->snacInfo, "\n\tSlope angle = %g degrees\n", slopeAngle/(M_PI/180.0) );
+
+ reg_dx = (geometry->max[0]-geometry->min[0])/(double)(full_I_node_range-1);
+ reg_dy = (geometry->max[1]-geometry->min[1])/(double)(full_J_node_range-1);
+ /*
+ * Decide at what height (prior to shearing) to switch from rombhoid to parallelopiped geometry,
+ * i.e., at what depth to resolve using "regular" elements parallel to the surface
+ * versus "irregular" elements deformed to fit to a horizontal bottom
+ * If contextExt->resolveDepth is negative, resolveHeight is set to zero and the whole domain is "regular"
+ */
+ if(contextExt->resolveDepth<0.0 || resolveHeight<0.0) {
+ resolveHeight = 0.0;
+ } else if(resolveHeight>geometry->max[1]) {
+ resolveHeight = geometry->max[1];
+ }
+ for( node_gI = 0; node_gI < full_I_node_range; node_gI++ )
+ new_x[node_gI] = geometry->min[0] + node_gI*reg_dx;
+
+ for( node_gJ = 0; node_gJ <= (full_J_node_range-1); node_gJ++ )
+ new_y[node_gJ] = geometry->min[1] + node_gJ*reg_dy;
+/* if((full_J_node_range-1)<4) { */
+/* for( node_gJ = 0; node_gJ <= (full_J_node_range-1); node_gJ++ ) */
+/* new_y[node_gJ] = geometry->min[1] + node_gJ*reg_dy; */
+/* } else { */
+/* for( node_gJ = 0; node_gJ <= (full_J_node_range-1)/4; node_gJ++ ) */
+/* new_y[node_gJ] = geometry->min[1] + node_gJ*reg_dy*(stretchFactor/(1.0+stretchFactor)); */
+/* for( ; node_gJ <= (full_J_node_range-1); node_gJ++ ) */
+/* new_y[node_gJ] = geometry->min[1] + ((full_J_node_range-1)/4)*reg_dy*(stretchFactor/(1.0+stretchFactor)) */
+/* + (node_gJ-((full_J_node_range-1)/4))*reg_dy*(1.0/(1.0+stretchFactor)); */
+/* } */
+ /*
+ * Assume highest point in x,y is at the last mesh node
+ */
+ new_xMin = new_x[0];
+ new_xMax = new_x[full_I_node_range-1];
+ new_xHalf = new_xMin+(new_xMax-new_xMin)/2.0;
+ new_yMin = new_y[0];
+ old_yMax = new_yMax = new_y[full_J_node_range-1]; /*+ tan(mesh_slope)*geometry->max[0]; */
+
+ xLeftFlatRamp=new_xMin+(new_xMax-new_xMin)*xLeftFlatRamp;
+ xRightRampFlat=new_xMin+(new_xMax-new_xMin)*xRightRampFlat;
+
+ fprintf(stderr, "Left flat <= %g -> %g\n", contextExt->leftFlatFraction, xLeftFlatRamp);
+ fprintf(stderr, "Right flat >= %g -> %g\n", contextExt->rightFlatFraction, xRightRampFlat);
+ fprintf(stderr, "Ramp-flat smoothness = %g\n", smoothFactor);
+ fprintf(stderr, "Flat/ramp/flat ranges: %g <-> %g <-> %g <-> %g\n", new_xMin, xLeftFlatRamp, xRightRampFlat, new_xMax);
+ /*
+ * Loop across all nodes in mesh
+ * - total of x*y*z nodes is given by context->mesh->nodeGlobalCount+1
+ */
+ for( node_g = 0; node_g < context->mesh->nodeGlobalCount; node_g++ ) {
+ Node_LocalIndex node_l = Mesh_NodeMapGlobalToLocal( mesh, node_g );
+ Index i_g;
+ Index j_g;
+ Index k_g;
+ Coord* coord = 0;
+ Coord tmpCoord;
+
+ /*
+ * If a local node, directly change the node coordinates else use a temporary location
+ */
+ if( node_l < context->mesh->nodeLocalCount ) { /* a local node */
+ coord = Snac_NodeCoord_P( context, node_l );
+ }
+ else {
+ /*
+ * If the node is running on another CPU, use dummy ptr
+ * - wouldn't it be better to just skip (continue) to next node?
+ */
+ coord = &tmpCoord;
+ }
+ /*
+ * Convert 1d mesh node index to 3d i,j,k mesh node position
+ */
+ RegularMeshUtils_Node_1DTo3D( decomp, node_g, &i_g, &j_g, &k_g );
+
+ /*
+ * Do nothing to x coordinate
+ */
+ (*coord)[0] = new_x[i_g];
+
+ /*exp(-smoothFactor*(xRightRampFlat-new_x[i_g])/old_yMax)
+ * Transform y coordinates by shearing mesh parallel to vertical axis
+ */
+ (*coord)[1] = new_y[j_g];
+ midOffset=( (xRightRampFlat-(xRightRampFlat-new_xHalf)*( 1.0-exp(-smoothFactor*(xRightRampFlat-new_xHalf)/old_yMax) ))
+ - (new_xHalf-xLeftFlatRamp)*( 1.0-exp(-smoothFactor*(new_xHalf-xLeftFlatRamp)/old_yMax) ) )*tan(slopeAngle);
+ if(new_x[i_g]<=new_xHalf) {
+ /* LHS */
+ (*coord)[1] +=
+ ( (new_x[i_g]<xLeftFlatRamp ? 0 :
+ (new_x[i_g]-xLeftFlatRamp)*( 1.0-exp(-smoothFactor*(new_x[i_g]-xLeftFlatRamp)/old_yMax) )*tan(slopeAngle) )
+ *(new_y[j_g]>=resolveHeight ? 1.0 : new_y[j_g]/resolveHeight) );
+ } else {
+ /* RHS */
+ (*coord)[1] +=
+ ( (new_x[i_g]>=xRightRampFlat ? xRightRampFlat*tan(slopeAngle)-midOffset :
+ (xRightRampFlat-(xRightRampFlat-new_x[i_g])*( 1.0-exp(-smoothFactor*(xRightRampFlat-new_x[i_g])/old_yMax) ))*tan(slopeAngle) - midOffset)
+ *(new_y[j_g]>=resolveHeight ? 1.0 : new_y[j_g]/resolveHeight) );
+
+ }
+ /*
+ * Track maximum y as mesh deforms
+ */
+ if((*coord)[1]>new_yMax) new_yMax=(*coord)[1];
+
+#ifdef DEBUG
+ fprintf( stderr, "\tnode_l: %2u, node_g: %2u, i: %2u, j: %2u, k: %2u, ",
+ node_l,
+ node_g,
+ i_g,
+ j_g,
+ k_g );
+ fprintf( stderr, "x: %12g, y: %12g, z: %12g\n", (*coord)[0], (*coord)[1], (*coord)[2] );
+#endif
+ }
+ /*
+ * Reset vertical max
+ */
+ geometry->max[1] = new_yMax;
+
+}
+
+
Property changes on: long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.c
___________________________________________________________________
Name: svn:executable
+ *
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/InitialConditions.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_InitialConditions_h__
+#define __SnacHillSlope_InitialConditions_h__
+
+ void _SnacHillSlope_InitialConditions( void* _context, void* data );
+/* double PlasticStrainFromCohesion(Snac_Material *material, double cohesion); */
+/* void ShellSort( unsigned int* vecPtr, unsigned int* idxPtr, unsigned long n ); */
+
+#endif /* __SnacHillSlope_InitialConditions_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Make.mm
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Make.mm (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Make.mm 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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 2921 2005-05-10 06:44:06Z EunseoChoi $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+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/hillSlope/Makefile.def
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Makefile.def (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Makefile.def 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,60 @@
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+##
+## 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.
+##
+## Mods by:
+## Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+##
+## 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 $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def_mod = SnacHillSlope
+def_inc = Snac/HillSlope
+
+def_srcs = \
+ Register.c \
+ ConstructExtensions.c \
+ Build.c \
+ Context.c \
+ InitialConditions.c \
+ Track.c \
+ CreateWeakPoints.c \
+ DeleteExtensions.c
+
+def_hdrs = \
+ types.h \
+ Register.h \
+ ConstructExtensions.h \
+ Build.h \
+ Context.h \
+ InitialConditions.h \
+ Track.h \
+ CreateWeakPoints.h \
+ DeleteExtensions.h \
+ HillSlope.h
+
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,142 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Context.h"
+#include "Track.h"
+#include "CreateWeakPoints.h"
+#include "ConstructExtensions.h"
+#include "Build.h"
+#include "InitialConditions.h"
+#include "DeleteExtensions.h"
+#include "Register.h"
+#include "Snac/Plastic/Plastic.h"
+
+#include <stdio.h>
+
+//#define DEBUG
+
+/* Textual name of this class */
+const Type SnacHillSlope_Type = "SnacHillSlope";
+
+ExtensionInfo_Index SnacHillSlope_ContextHandle = -1;
+
+
+Index _SnacHillSlope_Register( PluginsManager* pluginsMgr ) {
+ return PluginsManager_Submit( pluginsMgr,
+ SnacHillSlope_Type,
+ "0",
+ _SnacHillSlope_DefaultNew );
+}
+
+
+void* _SnacHillSlope_DefaultNew( Name name ) {
+ return _Codelet_New( sizeof(Codelet),
+ SnacHillSlope_Type,
+ _Codelet_Delete,
+ _Codelet_Print,
+ _Codelet_Copy,
+ _SnacHillSlope_DefaultNew,
+ _SnacHillSlope_Construct,
+ _Codelet_Build,
+ _Codelet_Initialise,
+ _Codelet_Execute,
+ _Codelet_Destroy,
+ name );
+}
+
+
+void _SnacHillSlope_Construct( void* component, Stg_ComponentFactory* cf, void* data ) {
+ Snac_Context* context;
+
+ #ifdef DEBUG
+ fprintf(stderr, "Entering Register.c...\n");
+ #endif
+
+ /* Retrieve context. */
+ context = (Snac_Context*)Stg_ComponentFactory_ConstructByName( cf, "context", Snac_Context, True, data );
+
+ Journal_Printf( context->debug, "In: %s\n", __func__ );
+
+
+ /* Add extensions to nodes, elements and the context */
+ SnacHillSlope_ContextHandle = ExtensionManager_Add(
+ context->extensionMgr,
+ SnacHillSlope_Type,
+ sizeof(SnacHillSlope_Context) );
+
+ /* Add extensions to the entry points */
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, AbstractContext_EP_Build ),
+ SnacHillSlope_Type,
+ _SnacHillSlope_Build,
+ SnacHillSlope_Type );
+
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, AbstractContext_EP_Initialise ),
+ SnacHillSlope_Type,
+ _SnacHillSlope_InitialConditions,
+ SnacHillSlope_Type );
+
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, Snac_EP_LoopNodesMomentum ),
+ SnacHillSlope_Type,
+ SnacHillSlope_Track,
+ SnacHillSlope_Type );
+
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, Snac_EP_LoopNodesMomentum ),
+ SnacHillSlope_Type,
+ SnacHillSlope_CreateWeakPoints,
+ SnacHillSlope_Type );
+
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, AbstractContext_EP_DestroyExtensions ),
+ SnacHillSlope_Type,
+ _SnacHillSlope_DeleteExtensions,
+ SnacHillSlope_Type );
+
+
+
+ /* Construct. */
+ _SnacHillSlope_ConstructExtensions( context, data );
+
+ #ifdef DEBUG
+ fprintf(stderr, "In Register.c\n");
+ #endif
+
+ #ifdef DEBUG
+ fprintf(stderr, "...leaving Register.c\n");
+ #endif
+
+
+}
Property changes on: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c
___________________________________________________________________
Name: svn:executable
+ *
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,44 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_Register_h__
+#define __SnacHillSlope_Register_h__
+
+ /* Textual name of this class */
+ extern const Type SnacHillSlope_Type;
+
+ /* Handles to extensions for quicker access */
+ extern ExtensionInfo_Index SnacHillSlope_ContextHandle;
+
+ Index _SnacHillSlope_Register( PluginsManager* pluginsMgr );
+
+ void* _SnacHillSlope_DefaultNew( Name name );
+
+ void _SnacHillSlope_Construct( void* component, Stg_ComponentFactory* cf, void* data );
+
+#endif /* __SnacHillSlope_Register_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,268 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ **
+ ** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+ ** 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: Track.c $
+ **
+ **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "Track.h"
+#include <string.h>
+#include <math.h>
+
+#ifndef TRUE
+#define TRUE 1
+#endif
+#ifndef FALSE
+#define FALSE 0
+#endif
+
+void SnacHillSlope_Track( void* _context ) {
+ Snac_Context *context = (Snac_Context*)_context;
+ SnacHillSlope_Context *contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+ Mesh *mesh = context->mesh;
+ MeshLayout *layout = (MeshLayout*)mesh->layout;
+ HexaMD *decomp = (HexaMD*)layout->decomp;
+ BlockGeometry *geometry = (BlockGeometry*)layout->elementLayout->geometry;
+ const int full_I_node_range=decomp->nodeGlobal3DCounts[0];
+ const int full_J_node_range=decomp->nodeGlobal3DCounts[1];
+ const int full_K_node_range=decomp->nodeGlobal3DCounts[2];
+ int index_I,index_J,index_K;
+
+ static double *yGridOldPtr, *yGridOlderPtr;
+ double max_yVelocity, max_yAcceln;
+ static double unit_yVelocity=0.0, unit_yAcceln=0.0;
+ static double old_max_yVelocity=0.0, old_max_yAcceln=0.0;
+/* static double min_yVelocity=0.0, min_yAcceln=0.0; */
+ static char fallingFlag=FALSE;
+
+ const double trackLevel=(double)contextExt->trackLevel;
+ const double startThreshold=(double)(contextExt->startThreshold<=1?
+ (contextExt->startThreshold>=0 ? contextExt->startThreshold : 0) : 1);
+ const double stopThreshold=(double)(contextExt->stopThreshold<=1?
+ (contextExt->stopThreshold>=0 ? contextExt->stopThreshold : 0) : 1);
+
+
+/* if (context->timeStep % context->dumpEvery == 0) { */
+/* Journal_Printf( context->snacInfo,"timeStep=%d (in track)\n", context->timeStep ); */
+/* } */
+ /*
+ * Set up a tracking grids (slices of mesh) to allow t instance to be compared with t-1, t-2 instances
+ */
+ if(context->timeStep==1) {
+ yGridOldPtr=(double *)malloc((size_t)(full_I_node_range*full_K_node_range*sizeof(double)));
+ yGridOlderPtr=(double *)malloc((size_t)(full_I_node_range*full_K_node_range*sizeof(double)));
+ for(index_I = 0; index_I < full_I_node_range; index_I++) {
+ for(index_K = 0; index_K < full_K_node_range; index_K++) {
+ *(yGridOlderPtr+index_I + full_I_node_range*index_K)=-1.0E19;
+ *(yGridOldPtr+index_I + full_I_node_range*index_K)=-1.0E19;
+ }
+ }
+ }
+
+ /*
+ * Compare t instance with t-1, t-2 instances, calculating slice vel, acceln and identifying maxima
+ */
+ max_yVelocity=0.0;
+ max_yAcceln=0.0;
+ for(index_I = 0; index_I < full_I_node_range; index_I++) {
+ for(index_K = 0; index_K < full_K_node_range; index_K++) {
+ Node_GlobalIndex node_gI;
+ Node_LocalIndex node_lI;
+ Coord *coordPtr = 0;
+ double node_yElevation, node_yVelocity, node_yAcceln;
+ double *tmp_yGridOldPtr, *tmp_yGridOlderPtr;
+
+ index_J=(int)(((double)full_J_node_range-1.0)*(1.0-trackLevel));
+ node_gI = index_I + full_I_node_range*index_J + full_I_node_range*full_J_node_range*index_K;
+ node_lI = Mesh_NodeMapGlobalToLocal( mesh, node_gI );
+/* Journal_Printf( context->snacInfo,"\t\t%d,%d,%d -> %d ->%d \n", index_I,index_J,index_K,node_gI,node_lI ); */
+/* node_lI = _MeshDecomp_Node_GlobalToLocal1D( decomp, node_gI ); */
+/* Journal_Printf( context->snacInfo,"\t\t%d,%d,%d -> %d ->%d \n", index_I,index_J,index_K,node_gI,node_lI ); */
+
+ /* If a local node, read its elevation, if not, give a dummy value */
+ if( node_lI < context->mesh->nodeLocalCount ) { /* a local node */
+ coordPtr = Snac_NodeCoord_P( context, node_lI );
+ node_yElevation = (*coordPtr)[1]+(*coordPtr)[1]+(*coordPtr)[1];
+ }
+ else {
+ node_yElevation=-1.0E29;
+ }
+ tmp_yGridOldPtr = yGridOldPtr+index_I + full_I_node_range*index_K;
+ tmp_yGridOlderPtr = yGridOlderPtr+index_I + full_I_node_range*index_K;
+ /*
+ * Start assessing the vertical surface motions after two time steps
+ */
+ node_yVelocity = node_yElevation-*tmp_yGridOldPtr;
+ node_yAcceln = node_yVelocity-(*tmp_yGridOldPtr-*tmp_yGridOlderPtr);
+ if(context->timeStep>=3){
+ if(fabs(node_yVelocity)>max_yVelocity)
+ max_yVelocity = fabs(node_yVelocity);
+ if(fabs(node_yAcceln)>max_yAcceln)
+ max_yAcceln = fabs(node_yAcceln);
+/* if(fabs(node_yVelocity)<min_yVelocity && fabs(node_yVelocity)>0.0) */
+/* min_yVelocity = fabs(node_yVelocity); */
+/* if(fabs(node_yAcceln)<min_yAcceln && fabs(node_yAcceln)>0.0) */
+/* min_yAcceln = fabs(node_yAcceln); */
+ }
+/* Journal_Printf( context->snacInfo, */
+/* "%d,%d,%d -> %d ->%d : y=%0.8f / oy=%0.8f / ooy=%0.8f -> v%g (%g) , a=%g (%g)\n", */
+/* index_I,index_J,index_K,node_gI,node_lI, */
+/* node_yElevation, *tmp_yGridOldPtr, *tmp_yGridOlderPtr, */
+/* node_yVelocity,max_yVelocity, node_yAcceln,max_yAcceln */
+/* ); */
+ /*
+ * Record this elevation field and push previous back to "OlderPtr" array
+ */
+ *tmp_yGridOlderPtr = *tmp_yGridOldPtr;
+ *tmp_yGridOldPtr = node_yElevation;
+ }
+ }
+/* if(!contextExt->startedTrackingFlag && unit_yVelocity>min_yVelocity && min_yVelocity<1.0) */
+/* unit_yVelocity = min_yVelocity; */
+/* if(!contextExt->startedTrackingFlag && unit_yAcceln>min_yAcceln && min_yAcceln<1.0) */
+/* unit_yAcceln = min_yAcceln; */
+ if(unit_yVelocity==0.0 && max_yVelocity>0.0)
+ unit_yVelocity = max_yVelocity;
+ if(unit_yAcceln==0.0 && max_yAcceln>0.0)
+ unit_yAcceln = max_yAcceln;
+ if(!contextExt->startedTrackingFlag && max_yVelocity>=unit_yVelocity && context->timeStep>=4)
+ contextExt->startedTrackingFlag=TRUE;
+
+ /*
+ * Decide whether to stop or to continue simulation
+ */
+ if(contextExt->startedTrackingFlag && !contextExt->elasticStabilizedFlag){
+
+ fallingFlag = CheckFallingFn(max_yVelocity,max_yAcceln,old_max_yVelocity,old_max_yAcceln);
+ if(CheckStabilizingFn(max_yVelocity/unit_yVelocity,
+ max_yAcceln/unit_yAcceln, stopThreshold, fallingFlag)==TRUE
+ && context->maxTimeSteps!=context->timeStep) {
+ /*
+ * Stabilizing! Therefore terminate processing at the next time step
+ */
+ contextExt->elasticStabilizedFlag=TRUE;
+ if(contextExt->solveElasticEqmOnlyFlag)
+ context->maxTimeSteps=context->timeStep+1;
+
+/* Journal_Printf( context->snacInfo,"Stabilizing (falling?=%d) on level %d (%g)\n", */
+/* fallingFlag, index_J, trackLevel ); */
+ } else {
+ /*
+ * Still working...
+ */
+/* if (context->timeStep % context->dumpEvery == 0) { */
+/* Journal_Printf( context->snacInfo,"Changing (rising?=%d) on level %d (%g) (t steps=%d/%d)\n", */
+/* fallingFlag, index_J, trackLevel, */
+/* context->timeStep, context->maxTimeSteps ); */
+/* } */
+ }
+ }
+
+
+/* if (context->timeStep % context->dumpEvery == 0) { */
+/* Journal_Printf( context->snacInfo,"\tTracking?=%d unit vel=%g, unit accel=%g\n", */
+/* contextExt->startedTrackingFlag, */
+/* unit_yVelocity/geometry->max[1], unit_yAcceln/geometry->max[1] */
+/* ); */
+/* Journal_Printf( context->snacInfo,"\tmax vel=%g (was %g) , max accel=%g (was %g)\n", */
+/* max_yVelocity/(unit_yVelocity>0.0?unit_yVelocity:1), */
+/* old_max_yVelocity/(unit_yVelocity>0.0?unit_yVelocity:1), */
+/* max_yAcceln/(unit_yAcceln>0.0?unit_yAcceln:1), */
+/* old_max_yAcceln/(unit_yAcceln>0.0?unit_yAcceln:1) */
+/* ); */
+/* Journal_Printf( context->snacInfo,"\tfalling?=%d tracking?=%d stabilized?=%d\n", */
+/* fallingFlag, contextExt->startedTrackingFlag, contextExt->elasticStabilizedFlag); */
+/* } */
+
+
+ /*
+ * Record the current mesh slice velocity and acceln for use next iteration
+ */
+ old_max_yVelocity = max_yVelocity;
+ old_max_yAcceln = max_yAcceln;
+
+
+#ifdef DEBUG
+#endif
+}
+
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * CheckFallingFn --
+ *
+ * Decide whether solution is heading away from or towards equilibrium
+ *
+ * Returns:
+ * True=falling/False=rising
+ *
+ * Side effects:
+ * None
+ *
+ *----------------------------------------------------------------------
+ */
+
+char
+CheckFallingFn(double velocity, double acceln, double old_velocity, double old_acceln)
+{
+ /*
+ * Return TRUE if either velocity or acceln is dropping or constant
+ */
+/* fprintf(stderr, "%g->%g %g->%g\n", old_velocity,velocity,old_acceln,acceln); */
+ return ( ((velocity<old_velocity) && (acceln<old_acceln)) ? TRUE : FALSE );
+}
+
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * CheckStabilizingFn --
+ *
+ * Decide whether solution is stabilizing
+ *
+ * Returns:
+ * True=stabilizing/False=still changing
+ *
+ * Side effects:
+ * None
+ *
+ *----------------------------------------------------------------------
+ */
+
+char
+CheckStabilizingFn(double velocity, double acceln, double stopThreshold, char fallingFlag)
+{
+ return ( fallingFlag & (velocity<=stopThreshold /* && acceln<stopThreshold */ ? TRUE : FALSE) );
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,44 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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
+**
+** Role:
+**
+** Assumptions:
+** None as yet.
+**
+** Comments:
+** None as yet.
+**
+** $Id: Track.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Snac_Track_h__
+#define __Snac_Track_h__
+
+void SnacHillSlope_Track( void* context );
+char CheckStabilizingFn(double velocity, double acceln, double stopThreshold, char risingOrFallingFlag);
+char CheckFallingFn(double velocity, double acceln, double old_velocity, double old_acceln);
+
+#endif /* __Snac_Track_h__ */
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/ConstructExtensions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/ConstructExtensions.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,111 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "InitialConditions.h"
+#include "ConstructExtensions.h"
+#include <assert.h>
+#include <limits.h>
+
+//#define DEBUG
+
+void _SnacHillSlope_ConstructExtensions( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacHillSlope_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+/* Dictionary* hillSlopeBCsDict; */
+/* char tmpBuf[PATH_MAX]; */
+
+ /* Because hillSlope is not an array by itself, we must the "complex" constructor for Variable... the info needs to be
+ * wrapped this generic way... */
+/* Index hillSlopeOffsetCount = 1; */
+/* SizeT hillSlopeOffsets[] = { (SizeT)((char*)&tmpNodeExt->hillSlope - (char*)&tmpNode) }; */
+/* Variable_DataType hillSlopeDataTypes[] = { Variable_DataType_Double }; */
+/* Index hillSlopeDataTypeCounts[] = { 1 }; */
+ #ifdef DEBUG
+ printf( "In %s()\n", __func__ );
+ #endif
+
+
+ /* Create the StGermain variable hillSlope, which is stored on a node extension */
+/* Variable_New( */
+/* "hillSlope", */
+/* hillSlopeOffsetCount, */
+/* hillSlopeOffsets, */
+/* hillSlopeDataTypes, */
+/* hillSlopeDataTypeCounts, */
+/* 0, */
+/* &ExtensionManager_GetFinalSize( context->mesh->nodeExtensionMgr ), */
+/* &context->mesh->layout->decomp->nodeDomainCount, */
+/* (void**)&context->mesh->node, */
+/* context->variable_Register ); */
+
+
+
+ /* HillSlope variables */
+ contextExt->slopeAngle = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "slopeAngle", Dictionary_Entry_Value_FromDouble( 30.0f ) ) );
+
+/* contextExt->rngSeed = Dictionary_Entry_Value_AsUnsignedInt( */
+/* Dictionary_GetDefault( context->dictionary, "rngSeed", Dictionary_Entry_Value_FromUnsignedInt( 1 ) ) ); */
+
+/* contextExt->fractionPlasticSeeds = Dictionary_Entry_Value_AsDouble( */
+/* Dictionary_GetDefault( context->dictionary, "fractionPlasticSeeds", Dictionary_Entry_Value_FromDouble( 0.02f ) ) ); */
+
+
+
+
+ /* Build the hillSlope IC and BC managers */
+/* hillSlopeBCsDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( context->dictionary, "hillSlopeBCs" ) ); */
+/* contextExt->hillSlopeBCs = CompositeVC_New("tempBC", */
+/* context->variable_Register, */
+/* context->condFunc_Register, */
+/* hillSlopeBCsDict, */
+/* context->mesh ); */
+
+
+/* #ifdef DEBUG */
+/* fprintf( stderr, "In %s()\n", __func__ ); */
+/* #endif */
+
+/* Journal_Printf( context->debug, "slopeAngle: %g\n", contextExt->slopeAngle ); */
+
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Context.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Context.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,60 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_Context_h__
+#define __SnacHillSlope_Context_h__
+
+ /* Context Information */
+ struct _SnacHillSlope_Context {
+ /* For _SnacHillSlope_Top2BottomSweep condition function */
+ double slopeAngle;
+/* unsigned int rngSeed; */
+/* double fractionWeakPoints; */
+
+ };
+
+ /* Print the contents of the context extension */
+ void SnacHillSlope_Context_Print( void* _context, Stream* stream );
+
+#endif /* __SnacHillSlope_Context_h__ */
+
+
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,261 @@
+s/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <time.h>
+#include <math.h>
+
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Context.h"
+#include "Register.h"
+#include "InitialConditions.h"
+#include "Snac/Plastic/Plastic.h"
+
+
+#ifndef PI
+ #ifndef M_PIl
+ #ifndef M_PI
+ #define PI 3.1415927
+ #else
+ #define PI M_PI
+ #endif
+ #else
+ #define PI M_PIl
+ #endif
+#endif
+
+#define min(a,b) ((a)<(b) ? (a):(b))
+
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+
+//#define DEBUG
+
+void _SnacHillSlope_InitialConditions( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacHillSlope_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+ Mesh* mesh = context->mesh;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+ BlockGeometry* geometry = (BlockGeometry*)layout->elementLayout->geometry;
+ Node_GlobalIndex node_g;
+ Node_GlobalIndex node_gI;
+ Node_GlobalIndex node_gJ;
+ double new_x[decomp->nodeGlobal3DCounts[0]];
+ const double reg_dx = (geometry->max[0]-geometry->min[0])/(decomp->nodeGlobal3DCounts[0] - 1);
+ double new_y[decomp->nodeGlobal3DCounts[1]];
+ const double reg_dy = (geometry->max[1]-geometry->min[1])/(decomp->nodeGlobal3DCounts[1] - 1);
+ double new_xMax;
+ double new_xMin;
+ double new_yMax;
+ double new_yMin;
+ double slope = contextExt->slopeAngle * M_PI/180.0;
+/* unsigned int rngSeed = contextExt->rngSeed; */
+/* double fractionWeakPoints = contextExt->fractionWeakPoints; */
+
+ int restart = 0;
+ Dictionary_Entry_Value* pluginsList;
+ Dictionary_Entry_Value* plugin;
+
+ Element_LocalIndex element_lI;
+ Element_GlobalIndex element_gI;
+ int index_I,index_J,global_I_range,global_J_range,global_K_range;
+ IJK ijk;
+ Snac_Element* element;
+ const Snac_Material* material;
+ Tetrahedra_Index tetra_I;
+ SnacPlastic_Element* plasticElement;
+ float rand_num;
+ int NUM_SEED;
+ int incZ = 0;
+
+ pluginsList = PluginsManager_GetPluginsList( context->dictionary );
+ if (pluginsList) {
+ plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList);
+ while ( plugin ) {
+ if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ),
+ "SnacRestart" ) ) {
+ restart = 1;
+ break;
+ }
+ plugin = plugin->next;
+ }
+ }
+ if( restart )
+ return;
+
+
+ #ifdef DEBUG
+ printf( "In: %s\n", __func__ );
+ fprintf(stderr, "Slope angle = %g degrees\n", slope/(M_PI/180.0));
+ #endif
+
+ /* Report HillSlope plugin variables picked up (?) from xml parameter file */
+ Journal_Printf( context->snacInfo, "\n\tSlope angle = %g degrees\n", contextExt->slopeAngle );
+/* Journal_Printf( context->snacInfo, "\tRNG seed = %u\n", contextExt->rngSeed ); */
+/* Journal_Printf( context->snacInfo, "\tFraction of weak points = %g\n\n", contextExt->fractionWeakPoints ); */
+
+
+ for( node_gI = 0; node_gI < decomp->nodeGlobal3DCounts[0]; node_gI++ )
+ new_x[node_gI] = geometry->min[0] + node_gI*reg_dx;
+
+ for( node_gJ = 0; node_gJ < decomp->nodeGlobal3DCounts[1]; node_gJ++ )
+ new_y[node_gJ] = geometry->min[1] + node_gJ*reg_dy;
+
+ /*
+ for( node_gJ = 0; node_gJ < decomp->nodeGlobal3DCounts[1]; node_gJ++ ) {
+ double reg_y;
+ if(node_gJ==0)
+ new_y[node_gJ] = geometry->min[1];
+ else if(node_gJ>0) {
+ reg_y = geometry->min[1] + node_gJ*reg_dy;
+ new_y[node_gJ] = new_y[node_gJ-1] + reg_dy/(5.0/PI*atan((reg_y-ym)/Ly)+3.0);
+ }
+ }
+ */
+
+
+ /*
+ * Assume highest point in x,y is at the last mesh node
+ */
+ new_xMin = new_x[0];
+ new_xMax = new_x[(decomp->nodeGlobal3DCounts[0] - 1)];
+ new_yMin = new_y[0];
+ new_yMax = new_y[(decomp->nodeGlobal3DCounts[1] - 1)]; /*+ tan(mesh_slope)*geometry->max[0]; */
+ /*
+ * Loop across all nodes in mesh
+ * - total of x*y*z nodes is given by context->mesh->nodeGlobalCount+1
+ */
+ for( node_g = 0; node_g < context->mesh->nodeGlobalCount; node_g++ ) {
+ Node_LocalIndex node_l = _MeshDecomp_Node_GlobalToLocal1D( decomp, node_g );
+ Index i_g;
+ Index j_g;
+ Index k_g;
+ Coord* coord = 0;
+ Coord tmpCoord;
+
+ /* If a local node, directly change the node coordinates and initial tpr, else use a temporary location */
+ if( node_l < context->mesh->nodeLocalCount ) { /* a local node */
+ coord = Snac_NodeCoord_P( context, node_l );
+ }
+ else {
+ coord = &tmpCoord;
+ }
+
+ RegularMeshUtils_Node_1DTo3D( decomp, node_g, &i_g, &j_g, &k_g );
+
+ (*coord)[0] = new_x[i_g];
+ (*coord)[0] *= geometry->max[0]/new_xMax;
+
+ (*coord)[1] = new_y[j_g] + tan(slope)*(*coord)[0];
+ (*coord)[1] *= geometry->max[1]/new_yMax;
+
+ //fprintf(stderr,"i_g=%d coord=%e newX=%e min=%e max=%e\n",
+ // i_g,(*coord)[0],new_x[i_g],geometry->min[0],geometry->max[0]);
+
+ #ifdef DEBUG
+ printf( "\tnode_l: %2u, node_g: %2u, i: %2u, j: %2u, k: %2u, ",
+ node_l,
+ node_g,
+ i_g,
+ j_g,
+ k_g );
+ printf( "x: %12g, y: %12g, z: %12g\n", (*coord)[0], (*coord)[1], (*coord)[2] );
+ #endif
+ }
+
+/* srand( rngSeed ); */
+ srand( 1 );
+ global_I_range = decomp->elementGlobal3DCounts[0]/4;
+ global_J_range = 1; //decomp->elementGlobal3DCounts[1];
+ global_K_range = 1; //decomp->elementGlobal3DCounts[2]; //10;
+ NUM_SEED = decomp->elementGlobal3DCounts[2] * global_I_range * 0.02;
+
+ for(index_I = 0; index_I < NUM_SEED; index_I++) {
+ rand_num = rand()/(RAND_MAX+1.0);
+ ijk[2] = decomp->elementGlobal3DCounts[2]*rand_num;
+ rand_num = rand()/(RAND_MAX+1.0);
+ ijk[0] = global_I_range * rand_num;
+ ijk[0] = decomp->elementGlobal3DCounts[0]/2 - (ijk[0]-global_I_range/2);
+
+ incZ = 0;
+ //ijk[2] -= global_K_range/2;
+ while( incZ < global_K_range && ijk[2] >=0 && ijk[2] < decomp->elementGlobal3DCounts[2] ) {
+ ijk[1] = decomp->elementGlobal3DCounts[1]-1;
+ for(index_J = 0; index_J < global_J_range; index_J++) {
+ element_gI = ijk[0] + decomp->elementGlobal3DCounts[0]*ijk[1] + decomp->elementGlobal3DCounts[0]*decomp->elementGlobal3DCounts[1]*ijk[2];
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+ if(element_lI < mesh->elementLocalCount) {
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get(
+ mesh->elementExtensionMgr,
+ element,
+ SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = 1.1 * material->plstrain[1];
+
+ Journal_Printf(
+ context->snacInfo,
+ "timeStep=%d ijkt=%d %d %d %d plasticE=%e\n",
+ context->timeStep,
+ ijk[0],
+ ijk[1],
+ ijk[2],
+ tetra_I,
+ plasticElement->plasticStrain[tetra_I] );
+ }
+ //fprintf(stderr,"me=%d seeds(%d/%d) ijk=%d %d %d\n",context->rank,index_I,NUM_SEED,ijk[0],ijk[1],ijk[2]);
+ }
+ ijk[1]--;
+ }
+ ijk[2]++;
+ incZ++;
+ }
+ }
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c.preshellsort
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c.preshellsort (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c.preshellsort 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,433 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "InitialConditions.h"
+#include "Snac/Plastic/Plastic.h"
+
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#ifndef PI
+ #ifndef M_PIl
+ #ifndef M_PI
+ #define PI 3.1415927
+ #else
+ #define PI M_PI
+ #endif
+ #else
+ #define PI M_PIl
+ #endif
+#endif
+
+
+
+#define min(a,b) ((a)<(b) ? (a):(b))
+
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+
+//void _SnacHillSlope_ShellSort( unsigned int* vecPtr, unsigned int* idxPtr, long n );
+
+//#define DEBUG
+
+void _SnacHillSlope_InitialConditions( void* _context, void* data ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacHillSlope_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacHillSlope_ContextHandle );
+ Mesh* mesh = context->mesh;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+ BlockGeometry* geometry = (BlockGeometry*)layout->elementLayout->geometry;
+ Node_GlobalIndex node_g;
+ Node_GlobalIndex node_gI;
+ Node_GlobalIndex node_gJ;
+
+ const int full_I_node_range=decomp->nodeGlobal3DCounts[0];
+ const int full_J_node_range=decomp->nodeGlobal3DCounts[1];
+ // const int full_K_node_range=decomp->nodeGlobal3DCounts[2];
+
+ const int full_I_element_range=decomp->elementGlobal3DCounts[0];
+ const int full_J_element_range=decomp->elementGlobal3DCounts[1];
+ const int full_K_element_range=decomp->elementGlobal3DCounts[2];
+
+ double new_x[full_I_node_range];
+ double reg_dx;
+ double new_y[full_J_node_range];
+ double reg_dy;
+ double new_xMax;
+ double new_xMin;
+ double new_yMax;
+ double new_yMin;
+
+ double slopeAngle = contextExt->slopeAngle * M_PI/180.0;
+ unsigned int rngSeed = contextExt->rngSeed;
+ double fractionWeakPoints = contextExt->fractionWeakPoints;
+ double x_subdomain_fraction = contextExt->xSubDomainFraction;
+ double y_subdomain_fraction = contextExt->ySubDomainFraction;
+ double z_subdomain_fraction = contextExt->zSubDomainFraction;
+
+ int restart = 0;
+ Dictionary_Entry_Value* pluginsList;
+ Dictionary_Entry_Value* plugin;
+
+ Element_LocalIndex element_lI;
+ Element_GlobalIndex element_gI;
+ int index_I,index_J,subdomain_I_element_range,
+ subdomain_J_element_range,subdomain_K_element_range;
+ IJK ijk;
+ Snac_Element* element;
+ const Snac_Material* material;
+ Tetrahedra_Index tetra_I;
+ SnacPlastic_Element* plasticElement;
+ double randomNumber;
+ unsigned int numberSubDomainPoints;
+ unsigned int numberWeakPoints;
+
+ unsigned int* shuffleIndexListPtr;
+ unsigned int* randomNumberListPtr;
+
+ int imax=0;
+
+ pluginsList = PluginsManager_GetPluginsList( context->dictionary );
+ if (pluginsList) {
+ plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList);
+ while ( plugin ) {
+ if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ),
+ "SnacRestart" ) ) {
+ restart = 1;
+ break;
+ }
+ plugin = plugin->next;
+ }
+ }
+ if( restart )
+ return;
+
+
+ #ifdef DEBUG
+ printf( "In: %s\n", __func__ );
+ fprintf(stderr, "Slope angle = %g degrees\n", slopeAngle/(M_PI/180.0));
+ #endif
+
+ /* Report HillSlope plugin variables picked up (?) from xml parameter file */
+ Journal_Printf( context->snacInfo, "\n\tSlope angle = %g degrees\n", slopeAngle/(M_PI/180.0) );
+ Journal_Printf( context->snacInfo, "\tRNG seed = %u\n", rngSeed );
+ Journal_Printf( context->snacInfo, "\tFraction of weak points = %g\n", fractionWeakPoints );
+ Journal_Printf( context->snacInfo, "\tx subdomain fraction = %g\n", x_subdomain_fraction );
+ Journal_Printf( context->snacInfo, "\ty subdomain fraction = %g\n", y_subdomain_fraction );
+ Journal_Printf( context->snacInfo, "\tz subdomain fraction = %g\n", z_subdomain_fraction );
+
+ reg_dx = (geometry->max[0]-geometry->min[0])/(full_I_node_range-1);
+ reg_dy = (geometry->max[1]-geometry->min[1])/(full_J_node_range-1);
+/* Journal_Printf( context->snacInfo, "\t new_x[1] = %g\n\n", geometry->min[0] + 1*reg_dx ); */
+/* Journal_Printf( context->snacInfo, "\t new_y[1] = %g\n\n", geometry->min[1] + 1*reg_dy ); */
+
+ for( node_gI = 0; node_gI < full_I_node_range; node_gI++ )
+ new_x[node_gI] = geometry->min[0] + node_gI*reg_dx;
+
+ for( node_gJ = 0; node_gJ < full_J_node_range; node_gJ++ )
+ new_y[node_gJ] = geometry->min[1] + node_gJ*reg_dy;
+
+ /*
+ for( node_gJ = 0; node_gJ < full_J_node_range; node_gJ++ ) {
+ double reg_y;
+ if(node_gJ==0)
+ new_y[node_gJ] = geometry->min[1];
+ else if(node_gJ>0) {
+ reg_y = geometry->min[1] + node_gJ*reg_dy;
+ new_y[node_gJ] = new_y[node_gJ-1] + reg_dy/(5.0/PI*atan((reg_y-ym)/Ly)+3.0);
+ }
+ }
+ */
+
+
+ /*
+ * Assume highest point in x,y is at the last mesh node
+ */
+ new_xMin = new_x[0];
+ new_xMax = new_x[full_I_node_range-1];
+ new_yMin = new_y[0];
+ new_yMax = new_y[full_J_node_range-1]; /*+ tan(mesh_slope)*geometry->max[0]; */
+ /*
+ * Loop across all nodes in mesh
+ * - total of x*y*z nodes is given by context->mesh->nodeGlobalCount+1
+ */
+ for( node_g = 0; node_g < context->mesh->nodeGlobalCount; node_g++ ) {
+ Node_LocalIndex node_l = _MeshDecomp_Node_GlobalToLocal1D( decomp, node_g );
+ Index i_g;
+ Index j_g;
+ Index k_g;
+ Coord* coord = 0;
+ Coord tmpCoord;
+
+ /* If a local node, directly change the node coordinates and initial tpr, else use a temporary location */
+ if( node_l < context->mesh->nodeLocalCount ) { /* a local node */
+ coord = Snac_NodeCoord_P( context, node_l );
+ }
+ else {
+ coord = &tmpCoord;
+ }
+
+ RegularMeshUtils_Node_1DTo3D( decomp, node_g, &i_g, &j_g, &k_g );
+
+ (*coord)[0] = new_x[i_g];
+ (*coord)[0] *= geometry->max[0]/new_xMax;
+
+ (*coord)[1] = new_y[j_g] + tan(slopeAngle)*(*coord)[0];
+ (*coord)[1] *= geometry->max[1]/new_yMax;
+
+ //fprintf(stderr,"i_g=%d coord=%e newX=%e min=%e max=%e\n",
+ // i_g,(*coord)[0],new_x[i_g],geometry->min[0],geometry->max[0]);
+
+ #ifdef DEBUG
+ fprintf( stderr, "\tnode_l: %2u, node_g: %2u, i: %2u, j: %2u, k: %2u, ",
+ node_l,
+ node_g,
+ i_g,
+ j_g,
+ k_g );
+ fprintf( stderr, "x: %12g, y: %12g, z: %12g\n", (*coord)[0], (*coord)[1], (*coord)[2] );
+ #endif
+ }
+
+
+ // return;
+
+ /*
+ * Define portion of mesh (cells) that need weak point "seeding"
+ * - lots of sloppy float/int casting back and forth here
+ */
+ subdomain_I_element_range = (int)(full_I_element_range*x_subdomain_fraction);
+ subdomain_J_element_range = (int)(full_J_element_range*y_subdomain_fraction);
+ subdomain_K_element_range = (int)(full_K_element_range*z_subdomain_fraction);
+
+ Journal_Printf( context->snacInfo, "\tx subdomain element range = %d/%d\n", subdomain_I_element_range,full_I_element_range );
+ Journal_Printf( context->snacInfo, "\ty subdomain element range = %d/%d\n", subdomain_J_element_range,full_J_element_range );
+ Journal_Printf( context->snacInfo, "\tz subdomain element range = %d/%d\n\n", subdomain_K_element_range,full_K_element_range );
+
+ numberSubDomainPoints = subdomain_I_element_range*subdomain_J_element_range*subdomain_K_element_range;
+ numberWeakPoints = (unsigned int)((float)numberSubDomainPoints*fractionWeakPoints);
+ Journal_Printf(context->snacInfo, "Number of weak points = %u/%u\n", numberWeakPoints,numberSubDomainPoints);
+
+ /*
+ * Seed RNG (random number generator) using seed value from input xml file
+ */
+ srand( rngSeed );
+ shuffleIndexListPtr=(unsigned int *)malloc((size_t)(numberSubDomainPoints*sizeof(unsigned int)));
+ randomNumberListPtr=(unsigned int *)malloc((size_t)(numberSubDomainPoints*sizeof(unsigned int)));
+ /*
+ * Create (1) ordered sequence of element indices, (2) parallel list of random variates
+ */
+ for(index_I = 0; index_I < numberSubDomainPoints; index_I++) {
+ shuffleIndexListPtr[index_I]=index_I;
+ randomNumberListPtr[index_I]=rand();
+ }
+ /*
+ *
+ */
+ ShellSort( randomNumberListPtr, shuffleIndexListPtr, numberSubDomainPoints );
+ /*
+ * Loop over total number of weak points and locate each one
+ */
+ for(index_I = 0; index_I < numberWeakPoints; index_I++) {
+ /*
+ * Select random i,j,k indices
+ * - converts int random number into double, scales, and casts back to int. Ugh.
+ */
+ randomNumber = (double)rand()/((double)RAND_MAX+1.0);
+ ijk[0] = (int)((subdomain_I_element_range) * randomNumber);
+ ijk[0] = ijk[0]+(full_I_element_range-subdomain_I_element_range)/2;
+
+ randomNumber = (double)rand()/((double)RAND_MAX+1.0);
+ ijk[1] = (int)((subdomain_J_element_range) * randomNumber);
+ ijk[1] = full_J_element_range-1-ijk[1];
+
+ randomNumber = (double)rand()/((double)RAND_MAX+1.0);
+ ijk[2] = (int)((subdomain_K_element_range) * randomNumber);
+ ijk[2] = ijk[2]+(full_K_element_range-subdomain_K_element_range)/2;
+
+ #ifdef DEBUG
+ if (imax<ijk[1]) imax=ijk[1];
+ #endif
+ /*
+ * Calculate "global" (across all processors) element index and then find local equivalent (this CPU)
+ */
+ element_gI = ijk[0] + full_I_element_range*ijk[1]
+ + full_I_element_range*full_J_element_range*ijk[2];
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+ #ifdef DEBUG
+ Journal_Printf(context->snacInfo, "Weak point: ijk = %d,%d,%d -> %d -> %d <? %d\n",
+ ijk[0], ijk[1], ijk[2], element_gI, element_lI, mesh->elementLocalCount);
+ #endif
+
+ /*
+ * At each point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get(
+ mesh->elementExtensionMgr,
+ element,
+ SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = contextExt->weakPointCohesion * material->plstrain[1];
+#ifdef DEBUG
+ if(tetra_I==0) Journal_Printf(
+ context->snacInfo,
+ "timeStep=%d ijkt=%d %d %d %d plasticE=%e\n",
+ context->timeStep,
+ ijk[0],
+ ijk[1],
+ ijk[2],
+ tetra_I,
+ plasticElement->plasticStrain[tetra_I] );
+#endif
+
+ }
+ } // End if
+
+ } // End for
+
+
+ /*
+ * Trigger point assignment
+ */
+ ijk[0]=(int)(full_I_element_range*contextExt->xTriggerPointFraction);
+ ijk[1]=(int)(full_J_element_range*contextExt->yTriggerPointFraction);
+ ijk[1]= full_J_element_range-ijk[1];
+ ijk[2]=(int)(full_K_element_range*contextExt->zTriggerPointFraction);
+ /*
+ * Calculate "global" (across all processors) element index and then find local equivalent (this CPU)
+ */
+ element_gI = ijk[0] + full_I_element_range*ijk[1]
+ + full_I_element_range*full_J_element_range*ijk[2];
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+ Journal_Printf(context->snacInfo, "Trigger point: ijk = %d,%d,%d -> %d -> %d <? %d\n", ijk[0], ijk[1], ijk[2], element_gI, element_lI, mesh->elementLocalCount);
+ #ifdef DEBUG
+ #endif
+ /*
+ * At trigger point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get(
+ mesh->elementExtensionMgr,
+ element,
+ SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = contextExt->triggerPointCohesion * material->plstrain[1];
+
+ if(tetra_I==0) Journal_Printf(
+ context->snacVerbose,
+ "timeStep=%d ijkt=%d %d %d %d plasticE=%e\n",
+ context->timeStep,
+ ijk[0],
+ ijk[1],
+ ijk[2],
+ tetra_I,
+ plasticElement->plasticStrain[tetra_I] );
+ }
+ }
+}
+
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * ShellSort --
+ *
+ * Shell sort vector (double) and coevally rearrange input
+ * (index, int) vector
+ *
+ * Results:
+ * None
+ *
+ * Side effects:
+ * Modifies input vectors
+ *
+ *----------------------------------------------------------------------
+ */
+
+void
+ShellSort(vecPtr,idxPtr,n)
+ unsigned int vecPtr[]; /* Input unsorted data vector (e.g. random numbers) */ /* was double */
+ unsigned int idxPtr[]; /* Input index vector, to be sorted according to vecPtr */
+ unsigned long n; /* Data length */
+{
+ unsigned long i,j,incr;
+ unsigned int v; /* was double */
+ unsigned int w;
+ incr=1;
+ do {
+ incr *= 3;
+ incr++;
+ } while (incr <= n);
+ do {
+ incr /= 3;
+ for (i=incr;i<n;i++) {
+ v=vecPtr[i];
+ w=idxPtr[i];
+ j=i;
+ while (vecPtr[j-incr] > v) {
+ vecPtr[j]=vecPtr[j-incr];
+ idxPtr[j]=idxPtr[j-incr];
+ j -= incr;
+ if (j < incr) {
+ break;
+ }
+ }
+ vecPtr[j]=v;
+ idxPtr[j]=w;
+ }
+ } while (incr > 1);
+}
Property changes on: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/InitialConditions.c.preshellsort
___________________________________________________________________
Name: svn:executable
+ *
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Register.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/Register.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,137 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "types.h"
+#include "Context.h"
+#include "ConstructExtensions.h"
+#include "Build.h"
+#include "InitialConditions.h"
+#include "DeleteExtensions.h"
+#include "Register.h"
+#include "Snac/Plastic/Plastic.h"
+
+#include <stdio.h>
+
+//#define DEBUG
+
+/* Textual name of this class */
+const Type SnacHillSlope_Type = "SnacHillSlope";
+
+ExtensionInfo_Index SnacHillSlope_ContextHandle = -1;
+
+
+Index _SnacHillSlope_Register( PluginsManager* pluginsMgr ) {
+ return PluginsManager_Submit( pluginsMgr,
+ SnacHillSlope_Type,
+ "0",
+ _SnacHillSlope_DefaultNew );
+}
+
+
+void* _SnacHillSlope_DefaultNew( Name name ) {
+ return _Codelet_New( sizeof(Codelet),
+ SnacHillSlope_Type,
+ _Codelet_Delete,
+ _Codelet_Print,
+ _Codelet_Copy,
+ _SnacHillSlope_DefaultNew,
+ _SnacHillSlope_Construct,
+ _Codelet_Build,
+ _Codelet_Initialise,
+ _Codelet_Execute,
+ _Codelet_Destroy,
+ name );
+}
+
+
+void _SnacHillSlope_Construct( void* component, Stg_ComponentFactory* cf, void* data ) {
+ Snac_Context* context;
+
+ #ifdef DEBUG
+ fprintf(stderr, "Entering Register.c...\n");
+ #endif
+
+ /* Retrieve context. */
+ context = (Snac_Context*)Stg_ComponentFactory_ConstructByName( cf, "context", Snac_Context, True, data );
+
+ Journal_Printf( context->debug, "In: %s\n", __func__ );
+
+
+ /* Add extensions to nodes, elements and the context */
+ SnacHillSlope_ContextHandle = ExtensionManager_Add(
+ context->extensionMgr,
+ SnacHillSlope_Type,
+ sizeof(SnacHillSlope_Context) );
+
+ /* Add extensions to the entry points */
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, AbstractContext_EP_Build ),
+ SnacHillSlope_Type,
+ _SnacHillSlope_Build,
+ SnacHillSlope_Type );
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, AbstractContext_EP_Initialise ),
+ SnacHillSlope_Type,
+ _SnacHillSlope_InitialConditions,
+ SnacHillSlope_Type );
+
+/* EntryPoint_Append( */
+/* Context_GetEntryPoint( context, Snac_EP_Constitutive ), */
+/* SnacHillSlope_Type, */
+/* _SnacHillSlope_Constitutive, */
+/* SnacHillSlope_Type ); */
+ EntryPoint_Append(
+ Context_GetEntryPoint( context, AbstractContext_EP_DestroyExtensions ),
+ SnacHillSlope_Type,
+ _SnacHillSlope_DeleteExtensions,
+ SnacHillSlope_Type );
+
+
+ /* Construct. */
+ _SnacHillSlope_ConstructExtensions( context, data );
+
+ #ifdef DEBUG
+ fprintf(stderr, "In Register.c\n");
+ #endif
+
+ #ifdef DEBUG
+ fprintf(stderr, "...leaving Register.c\n");
+ #endif
+
+
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,66 @@
+ /*
+ * Loop over total number of weak points and locate each one
+ */
+ for(index_I = 0; index_I < numberWeakPoints; index_I++) {
+ /*
+ * Select random i,j,k indices
+ * - converts int random number into double, scales, and casts back to int. Ugh.
+ */
+ randomNumber = (double)rand()/((double)RAND_MAX+1.0);
+ ijk[0] = (int)((subdomain_I_element_range) * randomNumber);
+ ijk[0] = ijk[0]+(full_I_element_range-subdomain_I_element_range)/2;
+
+ randomNumber = (double)rand()/((double)RAND_MAX+1.0);
+ ijk[1] = (int)((subdomain_J_element_range) * randomNumber);
+ ijk[1] = full_J_element_range-1-ijk[1];
+
+ randomNumber = (double)rand()/((double)RAND_MAX+1.0);
+ ijk[2] = (int)((subdomain_K_element_range) * randomNumber);
+ ijk[2] = ijk[2]+(full_K_element_range-subdomain_K_element_range)/2;
+
+ #ifdef DEBUG
+ if (imax<ijk[1]) imax=ijk[1];
+ #endif
+ /*
+ * Calculate "global" (across all processors) element index and then find local equivalent (this CPU)
+ */
+ element_gI = ijk[0] + full_I_element_range*ijk[1]
+ + full_I_element_range*full_J_element_range*ijk[2];
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+ #ifdef DEBUG
+ Journal_Printf(context->snacInfo, "Weak point: ijk = %d,%d,%d -> %d -> %d <? %d\n",
+ ijk[0], ijk[1], ijk[2], element_gI, element_lI, mesh->elementLocalCount);
+ #endif
+
+ /*
+ * At each point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get(
+ mesh->elementExtensionMgr,
+ element,
+ SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = contextExt->weakPointCohesion * material->plstrain[1];
+#ifdef DEBUG
+ if(tetra_I==0) Journal_Printf(
+ context->snacInfo,
+ "timeStep=%d ijkt=%d %d %d %d plasticE=%e\n",
+ context->timeStep,
+ ijk[0],
+ ijk[1],
+ ijk[2],
+ tetra_I,
+ plasticElement->plasticStrain[tetra_I] );
+#endif
+
+ }
+ } // End if
+
+ } // End for
+
+
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts2.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts2.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/cuts2.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,164 @@
+ if(fractionWeakPoints>=0.0 && fractionWeakPoints<=1.0) {
+ unsigned int rngSeed = contextExt->rngSeed;
+ double x_subdomainFraction = contextExt->xSubDomainFraction;
+ double y_subdomainFraction = contextExt->ySubDomainFraction;
+ double z_subdomainFraction = contextExt->zSubDomainFraction;
+ unsigned int numberSubDomainPoints,numberWeakPoints;
+ unsigned int *shuffleIndexListPtr, *randomNumberListPtr;
+ int subdomain_I_element_range,subdomain_J_element_range,subdomain_K_element_range;
+
+
+ /* Report HillSlope plugin variables picked up (?) from xml parameter file */
+ Journal_Printf( context->snacInfo, "\n\tSlope angle = %g degrees\n", slopeAngle/(M_PI/180.0) );
+ Journal_Printf( context->snacInfo, "\tRNG seed = %u\n", rngSeed );
+ Journal_Printf( context->snacInfo, "\tFraction of weak points = %g\n", fractionWeakPoints );
+ Journal_Printf( context->snacInfo, "\tx subdomain fraction = %g\n", x_subdomainFraction );
+ Journal_Printf( context->snacInfo, "\ty subdomain fraction = %g\n", y_subdomainFraction );
+ Journal_Printf( context->snacInfo, "\tz subdomain fraction = %g\n", z_subdomainFraction );
+
+ /*
+ * Define portion of mesh (cells) that need weak point "seeding"
+ * - lots of sloppy float/int casting back and forth here
+ */
+ subdomain_I_element_range = (int)(full_I_element_range*x_subdomainFraction);
+ subdomain_J_element_range = (int)(full_J_element_range*y_subdomainFraction);
+ subdomain_K_element_range = (int)(full_K_element_range*z_subdomainFraction);
+
+ Journal_Printf( context->snacInfo, "\tx subdomain element range = %d/%d\n", subdomain_I_element_range,full_I_element_range );
+ Journal_Printf( context->snacInfo, "\ty subdomain element range = %d/%d\n", subdomain_J_element_range,full_J_element_range );
+ Journal_Printf( context->snacInfo, "\tz subdomain element range = %d/%d\n\n", subdomain_K_element_range,full_K_element_range );
+
+ numberSubDomainPoints = subdomain_I_element_range*subdomain_J_element_range*subdomain_K_element_range;
+ numberWeakPoints = (unsigned int)((float)numberSubDomainPoints*fractionWeakPoints);
+ Journal_Printf(context->snacInfo, "Number of weak points = %u/%u\n", numberWeakPoints,numberSubDomainPoints);
+
+ /*
+ * Seed RNG (random number generator) using seed value from input xml file
+ */
+ srand( rngSeed );
+ shuffleIndexListPtr=(unsigned int *)malloc((size_t)(numberSubDomainPoints*sizeof(unsigned int)));
+ randomNumberListPtr=(unsigned int *)malloc((size_t)(numberSubDomainPoints*sizeof(unsigned int)));
+ /*
+ * Create (1) ordered sequence of element indices, (2) parallel list of random variates
+ */
+ index=0;
+ for(index_I = 0; index_I < subdomain_I_element_range; index_I++) {
+ for(index_J = 0; index_J < subdomain_J_element_range; index_J++) {
+ for(index_K = 0; index_K < subdomain_K_element_range; index_K++) {
+ /*
+ * Work out the element index from the i,j,k and the required location of the subdomain
+ */
+ element_gI = index_I+(full_I_element_range-subdomain_I_element_range)/2
+ + full_I_element_range*(full_J_element_range-1-index_J)
+ + full_I_element_range*full_J_element_range*(index_K+(full_K_element_range-subdomain_K_element_range)/2);
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+ shuffleIndexListPtr[index]=element_lI;
+ randomNumberListPtr[index]=rand();
+ index++;
+ if(index-1>numberSubDomainPoints) break;
+ }
+ if(index>0 && index-1>numberSubDomainPoints) break;
+ }
+ if(index>0 && index-1>numberSubDomainPoints) {fprintf(stderr, "We miscalculated the subdomain index\n"); break;}
+ }
+ /*
+ * Sort the list of sub-domain local (ie global set pared down for this CPU) indices according to the random numbers
+ * - i.e., randomize the index list
+ */
+ ShellSort( randomNumberListPtr, shuffleIndexListPtr, numberSubDomainPoints );
+#ifdef DEBUG
+ for(index = 0; index < numberSubDomainPoints; index++) {
+ Journal_Printf( context->snacInfo, "\t Subdomain point #%d = %d, RNG var=%d\n",
+ index, shuffleIndexListPtr[index],randomNumberListPtr[index] );
+ }
+#endif
+ /*
+ * Now the subdomain index list is shuffled and the random number list is ordered
+ */
+
+
+ /*
+ * Loop over the shuffled subdomain index list BUT only up to the number of weak points within it
+ * - i.e., subset the subdomain and allocate weak points to only a random fraction of them
+ */
+ for(index = 0; index < numberWeakPoints; index++) {
+ element_lI = shuffleIndexListPtr[index];
+
+#ifdef DEBUG
+ Journal_Printf(context->snacInfo, "Shuffled weak point # %d -> %d \n", index, element_lI);
+#endif
+
+ /*
+ * At each point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+ Snac_Element *element;
+ Snac_Material *material;
+ SnacPlastic_Element *plasticElement;
+ Tetrahedra_Index tetra_I;
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get( mesh->elementExtensionMgr, element, SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = PlasticStrainFromCohesion(material,(double)contextExt->weakPointCohesion);
+#ifdef DEBUG
+ if(tetra_I==0) Journal_Printf( context->snacInfo,"timeStep=%d et=%d %d plasticE: %e -> %e\n",context->timeStep, element_lI, tetra_I, material->plstrain[1], plasticElement->plasticStrain[tetra_I] );
+#endif
+
+ }
+ } // End if
+
+ } // End for
+ /*
+ * Free up allocated memory for shuffling arrays
+ */
+ free(shuffleIndexListPtr);
+ free(randomNumberListPtr);
+ } // End if (check fraction of weak points is in valid range [0,1]
+
+ /*
+ * Trigger point assignment
+ */
+ if(contextExt->triggerPointCohesion>=0.0) {
+ ijk[0]=(int)(full_I_element_range*contextExt->xTriggerPointFraction);
+ ijk[1]=(int)(full_J_element_range*contextExt->yTriggerPointFraction);
+ ijk[1]= full_J_element_range-1-ijk[1];
+ ijk[2]=(int)(full_K_element_range*contextExt->zTriggerPointFraction);
+ /*
+ * Calculate "global" (across all processors) element index and then find local equivalent (this CPU)
+ */
+ element_gI = ijk[0] + full_I_element_range*ijk[1] + full_I_element_range*full_J_element_range*ijk[2];
+ element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+#ifdef DEBUG
+ Journal_Printf(context->snacInfo, "Trigger point: ijk = %d,%d,%d -> %d -> %d <? %d\n", ijk[0], ijk[1], ijk[2], element_gI, element_lI, mesh->elementLocalCount);
+#endif
+ /*
+ * At trigger point, force low cohesion by imposing a degree of plastic strain
+ */
+ if(element_lI < mesh->elementLocalCount) {
+ Snac_Element *element;
+ Snac_Material *material;
+ SnacPlastic_Element *plasticElement;
+ Tetrahedra_Index tetra_I;
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get( mesh->elementExtensionMgr,element,SnacPlastic_ElementHandle );
+ /* Force each of 5*2 tetrahedra to have extra plastic strain to impose lower cohesion indirectly */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = PlasticStrainFromCohesion(material,(double)contextExt->triggerPointCohesion);
+#ifdef DEBUG
+ if(tetra_I==0) Journal_Printf(context->snacInfo, "timeStep=%d ijkt=%d %d %d %d plasticE=%e\n", context->timeStep,ijk[0],ijk[1],ijk[2], tetra_I, plasticElement->plasticStrain[tetra_I] );
+#endif
+ }
+ }
+ } // End if
+
+
+
+
+
+
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,422 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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 1095 2004-03-28 00:51:42Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "units.h"
+#include "types.h"
+#include "Register.h"
+#include "Context.h"
+#include "Node.h"
+#include "Element.h"
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+extern void effectiveDensity( void* _context );
+
+void _SnacTemperature_Top2BottomSweep( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* result ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
+ double* temperature = (double*)result;
+ MeshLayout* meshLayout = (MeshLayout*)context->meshLayout;
+ HexaMD* decomp = (HexaMD*)meshLayout->decomp;
+ Node_GlobalIndex node_gI;
+ IJK ijk;
+ const Index jCount = decomp->nodeGlobal3DCounts[1];
+
+ node_gI = context->mesh->nodeL2G[node_lI];
+ RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
+
+ *temperature = contextExt->bottomTemp + (contextExt->topTemp - contextExt->bottomTemp) / (jCount - 1) * ijk[1];
+
+ #if 0
+ printf( "Top: %8g, Bottom: %8g, node_I: %3u, node_gI: %3u, iijk: { %3u, %3u, %3u, }, *temperature: %g\n", contextExt->topTemp, contextExt->bottomTemp, node_lI, node_gI, ijk[0], ijk[1], ijk[2], *temperature );
+ #endif
+}
+
+
+void _SnacTemperature_Top2BottomSweep_Spherical(
+ Node_LocalIndex node_lI,
+ Variable_Index var_I,
+ void* _context,
+ void* result )
+{
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
+ Coord* coord = Snac_NodeCoord_P( context, node_lI );
+
+ double* temperature = (double*)result;
+ Dictionary* meshStruct= Dictionary_Entry_Value_AsDictionary(
+ Dictionary_Get( context->dictionary, "mesh" ) );
+ double rMin = Dictionary_Entry_Value_AsDouble( Dictionary_Get( meshStruct, "rMin" ) );
+ double rMax = Dictionary_Entry_Value_AsDouble( Dictionary_Get( meshStruct, "rMax" ) );
+ double r = sqrt((*coord)[0]*(*coord)[0] + (*coord)[1]*(*coord)[1] + (*coord)[2]*(*coord)[2]);
+
+ assert( (rMin != 0.0 && rMax != 0.0) );
+ r /= 6371000.0;
+ rMin /= 6371000.0;
+ rMax /= 6371000.0;
+
+ (*temperature) = contextExt->bottomTemp + (contextExt->topTemp - contextExt->bottomTemp) / (rMax - rMin) * (r - rMin);
+
+ #if 0
+ printf( "Top: %8g, Bottom: %8g, node_lI: %3u, node_gI: %3u, iijk: { %3u, %3u, %3u, }, *temperature: %g\n",
+ contextExt->topTemp,
+ contextExt->bottomTemp,
+ node_lI,
+ node_gI,
+ ijk[0],
+ ijk[1],
+ ijk[2],
+ *temperature );
+ #endif
+}
+
+void _SnacTemperature_Citcom_Compatible(
+ Node_LocalIndex node_lI,
+ Variable_Index var_I,
+ void* _context,
+ void* result )
+{
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
+ Mesh* mesh = context->mesh;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+ Node_GlobalIndex node_gI = _MeshDecomp_Node_LocalToGlobal1D( decomp, node_lI );
+ Coord* coord = Snac_NodeCoord_P( context, node_lI );
+
+ const Snac_Material* material = &context->materialProperty[0];
+ const double rTemp = contextExt->bottomTemp;
+ const double R = 6371000.0;
+ const double kappa = 1.0e-06;
+ double scalet = R*R/kappa/(1.0e+06*365.25*24.0*3600.0);
+ double age = 5.01f, temp=0.0f;
+
+ double* temperature = (double*)result;
+ Dictionary* meshStruct= Dictionary_Entry_Value_AsDictionary(
+ Dictionary_Get( context->dictionary, "mesh" ) );
+ double rMin = Dictionary_Entry_Value_AsDouble( Dictionary_Get( meshStruct, "rMin" ) );
+ double rMax = Dictionary_Entry_Value_AsDouble( Dictionary_Get( meshStruct, "rMax" ) );
+ double r = sqrt((*coord)[0]*(*coord)[0] + (*coord)[1]*(*coord)[1] + (*coord)[2]*(*coord)[2]);
+ IJK ijk;
+ const Node_GlobalIndex midI = (decomp->nodeGlobal3DCounts[0] + 1) / 2 - 1;
+ Node_GlobalIndex lmidI;
+
+ RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
+
+ /* for Cartesian case */
+ rMin = R - 4.0e+04f;
+ rMax = R;
+ r = rMax + (*coord)[1];
+ /*ccccc*/
+
+ assert( (rMin != 0.0 && rMax != 0.0) );
+ r /= R;
+ rMin /= R;
+ rMax /= R;
+
+ if( ijk[0] >= midI-5 && ijk[0] <= midI+5 ) {
+ age = 0.01f + 1.0f*abs(ijk[0]-midI)/5.0*5.01;
+ }
+
+ temp = (rMax-r) * 0.5f / sqrt(age/scalet);
+ *temperature = rTemp * erf(temp);
+ if( (*temperature) < 0.0) (*temperature) = 0.0f;
+}
+
+void _SnacTemperature_InitialConditions( void* _context ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
+ Element_LocalIndex element_lI;
+
+ int restart = 0;
+ Dictionary_Entry_Value* pluginsList;
+ Dictionary_Entry_Value* plugin;
+
+ pluginsList = PluginsManager_GetPluginsList( context->dictionary );
+ if (pluginsList) {
+ plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList);
+ while ( plugin ) {
+ if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ),
+ "SnacRestart" ) ) {
+ restart = 1;
+ break;
+ }
+ plugin = plugin->next;
+ }
+ }
+
+ Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
+
+ /* Temperature ICs are applied onto the "nodeICs" of Snac, and hence do not need to be repeated, but we must apply BCs */
+
+ /* In case of restarting, Temperature IC is still applied here */
+ if( restart ) {
+ FILE* fp;
+ Node_LocalIndex node_lI;
+ char path[PATH_MAX];
+
+ sprintf(path, "%s/snac.temp.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ Journal_Firewall( ( (fp = fopen(path,"r")) != NULL ), context->snacError, "Cannot find %s", path );
+
+ /* read in restart file to construct the initial temperature */
+ for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+ SnacTemperature_Node* nodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr, node, SnacTemperature_NodeHandle );
+ double T;
+ fscanf( fp, "%le", &T);
+
+ nodeExt->temperature = T;
+ nodeExt->temperature0 = nodeExt->temperature;
+ }
+ if( fp )
+ fclose( fp );
+ }
+ else {
+ Node_LocalIndex node_lI;
+
+ /* if not restarting, just copy initialized T to T0. */
+ for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+ SnacTemperature_Node* nodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr, node, SnacTemperature_NodeHandle );
+ nodeExt->temperature0 = nodeExt->temperature;
+ }
+ }
+ VariableCondition_Apply( contextExt->temperatureBCs, context );
+
+#if 0
+ {
+ Node_LocalIndex node_lI;
+
+ for( node_lI=0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+ SnacTemperature_Node* nodeExt = ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,
+ node,
+ SnacTemperature_NodeHandle );
+
+ printf( "node_lI: %3u, nodeExt->temperature: %g\n", node_lI, nodeExt->temperature );
+ }
+ }
+#endif
+
+ /* Apply element heat flux ICs */
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacTemperature_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacTemperature_ElementHandle );
+
+ memset( elementExt->heatFlux, 0, sizeof(elementExt->heatFlux) );
+ }
+
+ effectiveDensity( context );
+
+}
+
+void SnacTemperature_BoundaryConditions( void* _context ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
+
+ #if 0
+ Node_GlobalIndex node_gI, node_gI2;
+ Node_LocalIndex node_lI;
+ IJK ijk;
+ MeshLayout* meshLayout = (MeshLayout*)context->meshLayout;
+ HexaMeshTopology* topology = (HexaMeshTopology*)meshLayout->topology;
+ RegularMeshDecomp* decomp = (RegularMeshDecomp*)meshLayout->decomp;
+ #endif
+
+
+
+ Journal_Printf( context->debug, "In %s(): re-applying temp BCs\n", __func__ );
+
+ VariableCondition_Apply( contextExt->temperatureBCs, context );
+}
+
+void SnacTemperature_HeatFluxBC( void* _context ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ Element_LocalIndex element_lI;
+
+ Journal_Printf( context->debug, "In %s(): re-applying temp BCs\n", __func__ );
+
+ /* no heat flux bc for the vertical walls */
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ IJK ijk;
+ Tetrahedra_Index tetra_I;
+ MeshLayout* meshLayout = (MeshLayout*)context->meshLayout;
+ HexaMD* decomp = (HexaMD*)meshLayout->decomp;
+ Element_GlobalIndex element_gI = decomp->elementMapLocalToGlobal( decomp, element_lI );
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacTemperature_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacTemperature_ElementHandle );
+ double g[4],Fn,normal[3];
+ int i;
+
+ RegularMeshUtils_Element_1DTo3D( decomp, element_gI, &ijk[0], &ijk[1], &ijk[2] );
+
+ if( ijk[0] == 0 || ijk[0] == decomp->elementGlobal3DCounts[0]-1 || ijk[2] == 0 || ijk[2] == decomp->elementGlobal3DCounts[2]-1) {
+ if( ijk[0] == 0 ) {
+ normal[0] = 0.25f * (element->tetra[0].surface[2].normal[0] +
+ element->tetra[1].surface[2].normal[0] +
+ element->tetra[5].surface[0].normal[0] +
+ element->tetra[6].surface[0].normal[0]);
+ normal[1] = 0.25f * (element->tetra[0].surface[2].normal[1] +
+ element->tetra[1].surface[2].normal[1] +
+ element->tetra[5].surface[0].normal[1] +
+ element->tetra[6].surface[0].normal[1]);
+ normal[2] = 0.25f * (element->tetra[0].surface[2].normal[2] +
+ element->tetra[1].surface[2].normal[2] +
+ element->tetra[5].surface[0].normal[2] +
+ element->tetra[6].surface[0].normal[2]);
+ }
+ if(ijk[0] == decomp->elementGlobal3DCounts[0]-1) {
+ normal[0] = 0.25f * (element->tetra[2].surface[0].normal[0] +
+ element->tetra[3].surface[0].normal[0] +
+ element->tetra[7].surface[0].normal[0] +
+ element->tetra[8].surface[0].normal[0]);
+ normal[1] = 0.25f * (element->tetra[2].surface[0].normal[1] +
+ element->tetra[3].surface[0].normal[1] +
+ element->tetra[7].surface[0].normal[1] +
+ element->tetra[8].surface[0].normal[1]);
+ normal[2] = 0.25f * (element->tetra[2].surface[0].normal[2] +
+ element->tetra[3].surface[0].normal[2] +
+ element->tetra[7].surface[0].normal[2] +
+ element->tetra[8].surface[0].normal[2]);
+ }
+ if(ijk[2] == 0 ) {
+ normal[0] = 0.25f * (element->tetra[0].surface[1].normal[0] +
+ element->tetra[2].surface[2].normal[0] +
+ element->tetra[6].surface[1].normal[0] +
+ element->tetra[7].surface[1].normal[0]);
+ normal[1] = 0.25f * (element->tetra[0].surface[1].normal[1] +
+ element->tetra[2].surface[2].normal[1] +
+ element->tetra[6].surface[1].normal[1] +
+ element->tetra[7].surface[1].normal[1]);
+ normal[2] = 0.25f * (element->tetra[0].surface[1].normal[2] +
+ element->tetra[2].surface[2].normal[2] +
+ element->tetra[6].surface[1].normal[2] +
+ element->tetra[7].surface[1].normal[2]);
+ }
+ if(ijk[2] == decomp->elementGlobal3DCounts[2]-1) {
+ normal[0] = 0.25f * (element->tetra[1].surface[1].normal[0] +
+ element->tetra[3].surface[2].normal[0] +
+ element->tetra[5].surface[1].normal[0] +
+ element->tetra[8].surface[1].normal[0]);
+ normal[1] = 0.25f * (element->tetra[1].surface[1].normal[1] +
+ element->tetra[3].surface[2].normal[1] +
+ element->tetra[5].surface[1].normal[1] +
+ element->tetra[8].surface[1].normal[1]);
+ normal[2] = 0.25f * (element->tetra[1].surface[1].normal[2] +
+ element->tetra[3].surface[2].normal[2] +
+ element->tetra[5].surface[1].normal[2] +
+ element->tetra[8].surface[1].normal[2]);
+ }
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ Snac_Element_Tetrahedra* tetra = &element->tetra[tetra_I];
+ Snac_Node* contributingNode[4] = {
+ Snac_Element_Node_P( context, element_lI, TetraToNode[tetra_I][0] ),
+ Snac_Element_Node_P( context, element_lI, TetraToNode[tetra_I][1] ),
+ Snac_Element_Node_P( context, element_lI, TetraToNode[tetra_I][2] ),
+ Snac_Element_Node_P( context, element_lI, TetraToNode[tetra_I][3] ) };
+ SnacTemperature_Node* nodeExt[4] = {
+ ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,
+ contributingNode[0],
+ SnacTemperature_NodeHandle ),
+ ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,
+ contributingNode[1],
+ SnacTemperature_NodeHandle ),
+ ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,
+ contributingNode[2],
+ SnacTemperature_NodeHandle ),
+ ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,
+ contributingNode[3],
+ SnacTemperature_NodeHandle ),
+ };
+
+ g[0] = -1.0f * ( nodeExt[0]->temperature * tetra->surface[1].normal[0] * tetra->surface[1].area +
+ nodeExt[1]->temperature * tetra->surface[2].normal[0] * tetra->surface[2].area +
+ nodeExt[2]->temperature * tetra->surface[3].normal[0] * tetra->surface[3].area +
+ nodeExt[3]->temperature * tetra->surface[0].normal[0] * tetra->surface[0].area )
+ / 3.0f / element->tetra[tetra_I].volume;
+ g[1] = -1.0f * ( nodeExt[0]->temperature * tetra->surface[1].normal[1] * tetra->surface[1].area +
+ nodeExt[1]->temperature * tetra->surface[2].normal[1] * tetra->surface[2].area +
+ nodeExt[2]->temperature * tetra->surface[3].normal[1] * tetra->surface[3].area +
+ nodeExt[3]->temperature * tetra->surface[0].normal[1] * tetra->surface[0].area )
+ / 3.0f / element->tetra[tetra_I].volume;
+ g[2] = -1.0f * ( nodeExt[0]->temperature * tetra->surface[1].normal[2] * tetra->surface[1].area +
+ nodeExt[1]->temperature * tetra->surface[2].normal[2] * tetra->surface[2].area +
+ nodeExt[2]->temperature * tetra->surface[3].normal[2] * tetra->surface[3].area +
+ nodeExt[3]->temperature * tetra->surface[0].normal[2] * tetra->surface[0].area )
+ / 3.0f / element->tetra[tetra_I].volume;
+
+ Fn = 0.0f;
+ for( i=0; i<3;i++)
+ Fn += g[i] * normal[i];
+
+ elementExt->heatFlux[tetra_I][0] -= Fn*normal[0];
+ elementExt->heatFlux[tetra_I][1] -= Fn*normal[1];
+ elementExt->heatFlux[tetra_I][2] -= Fn*normal[2];
+ }
+ }
+ }
+}
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/VariableConditions.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,59 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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
+**
+** Role:
+** Handles the temperature initial and boundary conditions
+**
+** Assumptions:
+** None as yet.
+**
+** Comments:
+** None as yet.
+**
+** $Id: Heat.h 1095 2004-03-28 00:51:42Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Snac_Temperature_VariableConditions_h__
+#define __Snac_Temperature_VariableConditions_h__
+
+ void _SnacTemperature_InitialConditions( void* _context );
+ void SnacTemperature_BoundaryConditions( void* _context );
+ void SnacTemperature_HeatFluxBC( void* _context );
+ void _SnacTemperature_Top2BottomSweep( Node_LocalIndex node_lI, Variable_Index var_I, void* context, void* result );
+ void _SnacTemperature_Top2BottomSweep_Spherical(
+ Node_LocalIndex node_lI,
+ Variable_Index var_I,
+ void* context,
+ void* result );
+ void _SnacTemperature_Citcom_Compatible(
+ Node_LocalIndex node_lI,
+ Variable_Index var_I,
+ void* context,
+ void* result );
+
+#endif
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/cuts3.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/cuts3.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/tmp/cuts3.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,54 @@
+
+ /* Calculate tetrahedra heat flux (1st pass)*/
+ elementExt->heatFlux[tetra_I][0] = 1.0f * material->thermal_conduct * (
+ contributingNodeExt[0]->temperature *
+ element->tetra[tetra_I].surface[1].normal[0] *
+ element->tetra[tetra_I].surface[1].area +
+ contributingNodeExt[1]->temperature *
+ element->tetra[tetra_I].surface[2].normal[0] *
+ element->tetra[tetra_I].surface[2].area +
+ contributingNodeExt[2]->temperature *
+ element->tetra[tetra_I].surface[3].normal[0] *
+ element->tetra[tetra_I].surface[3].area +
+ contributingNodeExt[3]->temperature *
+ element->tetra[tetra_I].surface[0].normal[0] *
+ element->tetra[tetra_I].surface[0].area ) /
+ 3.0f / element->tetra[tetra_I].volume;
+ elementExt->heatFlux[tetra_I][1] = 1.0f * material->thermal_conduct * (
+ contributingNodeExt[0]->temperature *
+ element->tetra[tetra_I].surface[1].normal[1] *
+ element->tetra[tetra_I].surface[1].area +
+ contributingNodeExt[1]->temperature *
+ element->tetra[tetra_I].surface[2].normal[1] *
+ element->tetra[tetra_I].surface[2].area +
+ contributingNodeExt[2]->temperature *
+ element->tetra[tetra_I].surface[3].normal[1] *
+ element->tetra[tetra_I].surface[3].area +
+ contributingNodeExt[3]->temperature *
+ element->tetra[tetra_I].surface[0].normal[1] *
+ element->tetra[tetra_I].surface[0].area ) /
+ 3.0f / element->tetra[tetra_I].volume;
+ elementExt->heatFlux[tetra_I][2] = 1.0f * material->thermal_conduct * (
+ contributingNodeExt[0]->temperature *
+ element->tetra[tetra_I].surface[1].normal[2] *
+ element->tetra[tetra_I].surface[1].area +
+ contributingNodeExt[1]->temperature *
+ element->tetra[tetra_I].surface[2].normal[2] *
+ element->tetra[tetra_I].surface[2].area +
+ contributingNodeExt[2]->temperature *
+ element->tetra[tetra_I].surface[3].normal[2] *
+ element->tetra[tetra_I].surface[3].area +
+ contributingNodeExt[3]->temperature *
+ element->tetra[tetra_I].surface[0].normal[2] *
+ element->tetra[tetra_I].surface[0].area ) /
+ 3.0f / element->tetra[tetra_I].volume;
+
+ /* thermal stress */
+ if( context->computeThermalStress ) {
+ temp = 0.25f * ( contributingNodeExt[0]->temperature + contributingNodeExt[1]->temperature + contributingNodeExt[2]->temperature + contributingNodeExt[3]->temperature );
+ temp0 = 0.25f * ( contributingNodeExt[0]->temperature0 + contributingNodeExt[1]->temperature0 + contributingNodeExt[2]->temperature0 + contributingNodeExt[3]->temperature0 );
+ /* Thermal stresses (alfa_v = 3.e-5 1/K) */
+ /* elementExt->thermalStress[tetra_I] = -1.0f*alpha*bulkm*(context->timeStep==1?(1.0*(temp-1300.0)):(temp-temp0)); */
+ elementExt->thermalStress[tetra_I] = -1.0f*alpha*bulkm*(temp-temp0);
+ Journal_Firewall( (!isinf(elementExt->thermalStress[tetra_I]) && !isnan(elementExt->thermalStress[tetra_I])), " timeStep %d thermalStress for element %d tetra %d (T: %e %e) is inf or nan!!\n",context->timeStep,element_lI,tetra_I,temp,temp0);
+ }
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.c 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,79 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "Node.h"
+#include "Register.h"
+#include <stdio.h>
+
+void _SnacHillSlope_DumpShear( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ Node_LocalIndex node_lI;
+ StressTensor* stress = &element->tetra[tetra_I].stress;
+ StrainTensor* strain = &element->tetra[tetra_I].strain;
+ Strain plasticStrain = 0.0f;
+
+
+ /* Compute elastic stress first */
+ trace_strain = (*strain)[0][0] + (*strain)[1][1] + (*strain)[2][2];
+
+ (*stress)[0][0] += (2.0f * material->mu) * (*strain)[0][0] + material->lambda * (trace_strain );
+ (*stress)[1][1] += (2.0f * material->mu) * (*strain)[1][1] + material->lambda * (trace_strain );
+ (*stress)[2][2] += (2.0f * material->mu) * (*strain)[2][2] + material->lambda * (trace_strain );
+ (*stress)[0][1] += (2.0f * material->mu) * (*strain)[0][1];
+ (*stress)[0][2] += (2.0f * material->mu) * (*strain)[0][2];
+ (*stress)[1][2] += (2.0f * material->mu) * (*strain)[1][2];
+
+ principal_stresses(stress,s,cn);
+
+
+
+ if( context->timeStep == 0 || (context->timeStep - 1) % context->dumpEvery == 0 ) {
+ #if DEBUG
+ printf( "In %s()\n", __func__ );
+ #endif
+
+ for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+/* SnacHillSlope_Node* nodeExt = ExtensionManager_Get( */
+/* context->mesh->nodeExtensionMgr, */
+/* node, */
+/* SnacHillSlope_NodeHandle ); */
+ float shearStress = nodeExt->hillSlope;
+ fwrite( &hillSlope, sizeof(float), 1, contextExt->hillSlopeOut );
+ }
+ fflush( contextExt->hillSlopeOut );
+ }
+}
+
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/archive/unused/Output.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,34 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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 $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef _SnacHillSlope_Output_
+#define _SnacHillSlope_Output_
+
+ void _SnacHillSlope_DumpShear( void* _context );
+
+#endif
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/hillslope_demo1.xml
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/hillslope_demo1.xml (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/hillslope_demo1.xml 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,264 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+
+<!-- StGermain-Snac input file -->
+<!-- $Id: hillslope.xml $ -->
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+ <!-- StGermain simulation parameters -->
+ <param name="start"> 0 </param>
+ <param name="outputPath">results_demo1</param>
+ <param name="dumpEvery"> 10 </param>
+ <param name="maxTimeSteps"> 1000 </param> <!-- seconds -->
+
+ <!-- Snac variables -->
+ <param name="density"> 2500 </param> <!-- kg/m^3 -->
+ <param name="gravity"> 9.8 </param> <!-- m/s^2 -->
+ <param name="demf"> 0.8 </param> <!-- dimensionless force damping factor -->
+ <param name="dtType"> constant </param> <!-- time marching method: dynamic/constant -->
+ <!--<param name="timeStep"> 1e+02 </param> seconds/step (used if dtType==constant) -->
+ <param name="timeStep"> 1e+00 </param> <!-- seconds/step (used if dtType==constant) -->
+ <param name="forceCalcType"> complete </param>
+ <param name="decomposedAxis"> 0 </param> <!-- hack: 0=X, 1=Y, 2=Z. Should and will eventually be automatically discovered-->
+ <param name="storeForces"> no </param> <!-- store residual forces at each node? -->
+ <param name="forceCheckSum"> no </param> <!-- check sum of forces is consistent -->
+ <param name="topo_kappa"> 0.0 </param> <!-- parameter for topographic smoothing -->
+ <param name="alpha"> 0 </param> <!-- 1/K, volumetric thermal expansion coefficient -->
+
+ <!-- Extension modules -->
+ <!--
+ <param> SnacRemesher </param>
+ <param> SnacWinklerForce </param>
+ <param> SnacSoilPlSeeds </param>
+ <param> SnacTemperature </param>
+ -->
+ <list name="extensions">
+ <param> SnacPlastic </param>
+ <param> SnacHydroStaticIC </param>
+ <param> SnacHillSlope </param> <!-- Requires SnacPlastic -->
+ </list>
+
+ <struct name="mesh">
+ <param name="shadowDepth"> 1 </param>
+ <param name="decompDims"> 1 </param>
+
+ <!-- Mesh size --> <!-- subtract 1 to get # of elements -->
+ <param name="meshSizeI"> 20 </param> <param name="meshSizeJ"> 5 </param> <param name="meshSizeK"> 2 </param>
+ <param name="maxX"> 380 </param> <param name="maxY"> 100 </param> <param name="maxZ"> 20 </param>
+ <param name="minX"> 0 </param> <param name="minY"> 0 </param> <param name="minZ"> 0 </param>
+ <!-- Initial geometry - meters NB: maxY is presheared-mesh thickness -->
+
+<!--
+ <param name="meshSizeI"> 201 </param> <param name="meshSizeJ"> 51 </param> <param name="meshSizeK"> 2 </param>
+ <param name="meshSizeI"> 12 </param> <param name="meshSizeJ"> 7 </param> <param name="meshSizeK"> 2 </param>
+ <param name="meshSizeI"> 61 </param> <param name="meshSizeJ"> 41 </param> <param name="meshSizeK"> 2 </param>
+ <param name="meshSizeI"> 44 </param> <param name="meshSizeJ"> 36 </param> <param name="meshSizeK"> 2 </param>
+ <param name="meshSizeI"> 22 </param> <param name="meshSizeJ"> 12 </param> <param name="meshSizeK"> 18 </param>
+ <param name="meshSizeI"> 4 </param> <param name="meshSizeJ"> 5 </param> <param name="meshSizeK"> 3 </param>
+-->
+
+ <!-- Remeshing -->
+ <param name="meshType"> cartesian </param>
+ <param name="buildNodeNeighbourTbl"> True </param>
+ </struct>
+
+
+
+ <!-- HillSlope plugin parameters -->
+ <param name="slopeAngle"> 30.0 </param> <!-- degrees -->
+ <param name="rngSeed"> 1 </param> <!-- dimensionless -->
+ <param name="leftFlatFraction"> 0.05 </param> <!-- Fractional x-distance to bottom ramp-flat transition -->
+ <param name="rightFlatFraction"> 0.95 </param> <!-- Fractional x-distance to top flat-ramp transition -->
+ <param name="rampFlatSmoothFactor"> 3.0 </param> <!-- Smoothness of ramp-flat/flat-ramp transitions -->
+ <param name="resolveDepth"> 75 </param> <!-- meters, measured from top surface -->
+ <param name="elasticStabilizedFlag"> false </param> <!-- if true, weak seeding applied immediately -->
+
+ <param name="fractionWeakPoints"> 0.5 </param> <!-- dimensionless -->
+ <param name="xSubDomainFraction"> 1 </param> <!-- Horizontal (left-right) around center -->
+ <param name="ySubDomainFraction"> 0.9 </param> <!-- Sheared vertical (measured from top surface) -->
+ <param name="zSubDomainFraction"> 1.0 </param> <!-- Horizontal (back-front) around center -->
+ <param name="xTriggerPointFraction"> 0.5 </param> <!-- dimensionless -->
+ <param name="yTriggerPointFraction"> 0.45 </param> <!-- dimensionless, measured from top surface -->
+ <param name="zTriggerPointFraction"> 0.5 </param> <!-- dimensionless -->
+
+ <param name="weakPointCohesion"> 1.0e+05 </param> <!-- Pa -->
+ <param name="triggerPointCohesion"> 1.0e+04 </param> <!-- Pa -->
+ <param name="trackLevelFraction"> 0.0 </param> <!-- y mesh watch level fraction (measured from surface) -->
+ <param name="startThreshold"> 1.0 </param> <!-- internal units, elastic eqm reached if vel & acc drop below this threshold -->
+ <param name="stopThreshold"> 1.0 </param> <!-- internal units, stop run if vel & acc drop below this threshold -->
+ <param name="startedTrackingFlag"> false </param> <!-- must be false -->
+ <param name="solveElasticEqmOnlyFlag"> true </param> <!-- if true, solution stops at approx elastic eqm -->
+
+
+ <!-- Elastic material parameters -->
+ <param name="lambda"> 1.0e+7 </param>
+ <param name="mu"> 1.0e+7 </param>
+
+ <!-- Plastic material parameters -->
+ <param name="yieldcriterion"> mohrcoulomb </param>
+ <!-- define a piecewise coh-pl strain "curve" -->
+ <param name="nsegments"> 3 </param> <!-- curve has nsegments pieces, nsegments+1 nodes -->
+ <param name="plstrain0"> 0.0 </param> <!-- dimensionless -->
+ <param name="plstrain1"> 0.001 </param> <!-- dimensionless -->
+ <param name="plstrain2"> 0.002 </param> <!-- dimensionless -->
+ <param name="plstrain3"> 1.0e+05 </param> <!-- dimensionless -->
+ <param name="frictionAngle0"> 25 </param> <!-- degrees -->
+ <param name="frictionAngle1"> 25 </param> <!-- degrees -->
+ <param name="frictionAngle2"> 25 </param> <!-- degrees -->
+ <param name="frictionAngle3"> 25 </param> <!-- degrees -->
+ <param name="dilationAngle0"> 5.0 </param> <!-- degrees -->
+ <param name="dilationAngle1"> 5.0 </param> <!-- degrees -->
+ <param name="dilationAngle2"> 5.0 </param> <!-- degrees -->
+ <param name="dilationAngle3"> 5.0 </param> <!-- degrees -->
+ <param name="cohesion0"> 1.0e+06 </param> <!-- Pa -->
+ <param name="cohesion1"> 1.0e+05 </param> <!-- Pa -->
+ <param name="cohesion2"> 1.0e+04 </param> <!-- Pa -->
+ <param name="cohesion3"> 1.0e+00 </param> <!-- Pa -->
+ <param name="ten_off"> 1.0e+13 </param> <!-- limit tensional stress, yield forced above this -->
+
+ <!-- Remesher info -->
+ <!--
+ <param name="remeshCondition"> onBothTimeStepLength </param>
+ <param name="remeshCondition"> onTimeStep </param>
+ <param name="remeshCondition"> onMinLengthScale </param>
+ -->
+ <param name="remeshCondition"> onMinLengthScale </param>
+ <param name="remeshTimeStepCriterion"> 15000 </param>
+ <param name="remeshLengthCriterion"> 0.7 </param>
+ <param name="bottomRestore"> on </param>
+
+ <!--
+ <param name="fractionWeakPoints"> 1.0 </param>
+ <param name="xSubDomainFraction"> 0.5 </param>
+ <param name="ySubDomainFraction"> 0.5 </param>
+ <param name="zSubDomainFraction"> 0.5 </param>
+ <param name="xTriggerPointFraction"> 0.5 </param>
+ <param name="yTriggerPointFraction"> 0.4 </param>
+ <param name="zTriggerPointFraction"> 0.5 </param>
+ -->
+
+
+
+
+
+
+ <!-- node ICs -->
+ <struct name="nodeICs">
+ <list name="vcList">
+ <struct>
+ <param name="type"> AllNodesVC </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> vx </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ <struct>
+ <param name="name"> vy </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ <struct>
+ <param name="name"> vz </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
+ <!-- element ICs -->
+ <struct name="elementICs">
+ <list name="vcList">
+ <struct>
+ <param name="type"> AllElementsVC </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> elementMaterial </param>
+ <param name="type"> int </param>
+ <param name="value"> 0 </param> <!-- no dimensions -->
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
+ <!-- Velocity BCs -->
+ <struct name="velocityBCs">
+ <list name="vcList">
+ <struct>
+ <!-- "Free" left (x=xmin) boundary = unspecified vy, vz, zero vx -->
+ <param name="type"> WallVC </param>
+ <param name="wall"> left </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> vx </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ </list>
+ </struct>
+ <struct>
+ <!-- "Free" right (x=xmax) boundary = unspecified vy, vz, zero vx -->
+ <param name="type"> WallVC </param>
+ <param name="wall"> right </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> vx </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ </list>
+ </struct>
+ <struct>
+ <!-- "Free" top (y=ymax) boundary = unspecified vx, vy, vz -->
+ <!-- "Pinned" bottom (y=ymin) boundary = unspecified zero vx, vy, vz -->
+ <param name="type"> WallVC </param>
+ <param name="wall"> bottom </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> vx </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ <struct>
+ <param name="name"> vy </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ <struct>
+ <param name="name"> vz </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ </list>
+ </struct>
+ <struct>
+ <!-- "Free" back (z=zmin) boundary = unspecified vx, vy, zero vz -->
+ <param name="type"> WallVC </param>
+ <param name="wall"> back </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> vz </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ </list>
+ </struct>
+ <struct>
+ <!-- "Free" front (z=zmax) boundary = unspecified vx, vy, zero vz -->
+ <param name="type"> WallVC </param>
+ <param name="wall"> front </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> vz </param>
+ <param name="type"> double </param>
+ <param name="value"> 0.0 </param> <!-- m/s -->
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
+</StGermainData>
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/make.log
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/make.log (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/make.log 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,3 @@
+make[1]: Entering directory `/usr/local/SNAC/Snac/plugins/hillSlope'
+make[1]: Nothing to be done for `all_build'.
+make[1]: Leaving directory `/usr/local/SNAC/Snac/plugins/hillSlope'
Added: long/3D/SNAC/trunk/Snac/plugins/hillSlope/makefile
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/makefile (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/makefile 2009-02-26 17:48:54 UTC (rev 14159)
@@ -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 3049 2005-06-22 07:46:15Z AlanLo $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+# 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/hillSlope/types.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/types.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/types.h 2009-02-26 17:48:54 UTC (rev 14159)
@@ -0,0 +1,51 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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.
+**
+** Mods by:
+** Colin Stark, Doherty Research Scientist, Lamont-Doherty Earth Observatory (cstark at ldeo.columbia.edu)
+**
+** 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: types.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __SnacHillSlope_types_h__
+#define __SnacHillSlope_types_h__
+
+ /* Hillslope params */
+ typedef struct _SnacHillSlope_Context SnacHillSlope_Context;
+
+#endif /* __SnacHillSlope_types_h__ */
More information about the CIG-COMMITS
mailing list