[cig-commits] r8998 - in short/3D/PyLith/trunk/libsrc: . materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Jan 11 08:47:12 PST 2008
Author: willic3
Date: 2008-01-11 08:47:12 -0800 (Fri, 11 Jan 2008)
New Revision: 8998
Added:
short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.cc
short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.hh
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
Log:
Factored out function that is used by all linear Maxwell viscoelastic
models (including generalized Maxwell model) and put it in static
function in ViscoelasticMaxwell.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2008-01-10 18:06:09 UTC (rev 8997)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2008-01-11 16:47:12 UTC (rev 8998)
@@ -68,6 +68,7 @@
materials/ElasticPlaneStrain.cc \
materials/ElasticPlaneStress.cc \
materials/MaxwellIsotropic3D.cc \
+ materials/ViscoelasticMaxwell.cc \
meshio/BinaryIO.cc \
meshio/DataWriter.cc \
meshio/DataWriterVTK.cc \
Modified: short/3D/PyLith/trunk/libsrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2008-01-10 18:06:09 UTC (rev 8997)
+++ short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2008-01-11 16:47:12 UTC (rev 8998)
@@ -24,7 +24,8 @@
MaxwellIsotropic3D.hh \
MaxwellIsotropic3D.icc \
Material.hh \
- Material.icc
+ Material.icc \
+ ViscoelasticMaxwell.hh
noinst_HEADERS =
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-01-10 18:06:09 UTC (rev 8997)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-01-11 16:47:12 UTC (rev 8998)
@@ -14,6 +14,8 @@
#include "MaxwellIsotropic3D.hh" // implementation of object methods
+#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
+
#include "pylith/utils/array.hh" // USES double_array
#include "petsc.h" // USES PetscLogFlopsNoCheck
@@ -268,29 +270,9 @@
parameters[_MaxwellIsotropic3D::pidStrainT+2])/3.0;
PetscLogFlopsNoCheck(11);
- // The code below should probably be in a separate function since it
- // is used more than once. I should also probably cover the possibility
- // that Maxwell time is zero (although this should never happen).
- const double timeFrac = 1.0e-5;
- const int numTerms = 5;
- double dq = 0.0;
- if(maxwelltime < timeFrac*_dt) {
- double fSign = 1.0;
- double factorial = 1.0;
- double fraction = 1.0;
- dq = 1.0;
- for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
- factorial *= iTerm;
- fSign *= -1.0;
- fraction *= _dt/maxwelltime;
- dq += fSign*fraction/factorial;
- } // for
- PetscLogFlopsNoCheck(1+7*numTerms);
- } else {
- dq = maxwelltime*(1.0-exp(-_dt/maxwelltime))/_dt;
- PetscLogFlopsNoCheck(7);
- } // else
+ // Time integration.
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
const double expFac = exp(-_dt/maxwelltime);
const double elasFac = 2.0*mu;
double devStrainTpdt = 0.0;
@@ -396,28 +378,8 @@
const double mu2 = 2.0 * mu;
const double bulkmodulus = lambda + mu2/3.0;
- const double timeFrac = 1.0e-5;
- const int numTerms = 5;
- double dq = 0.0;
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
- PetscLogFlopsNoCheck(3);
- if(maxwelltime < timeFrac*_dt) {
- double fSign = 1.0;
- double factorial = 1.0;
- double fraction = 1.0;
- dq = 1.0;
- for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
- factorial *= iTerm;
- fSign *= -1.0;
- fraction *= _dt/maxwelltime;
- dq += fSign*fraction/factorial;
- } // for
- PetscLogFlopsNoCheck(1+7*numTerms);
- } else {
- dq = maxwelltime*(1.0-exp(-_dt/maxwelltime))/_dt;
- PetscLogFlopsNoCheck(7);
- } // else
-
const double visFac = mu*dq/3.0;
elasticConsts[ 0] = bulkmodulus + 4.0*visFac; // C1111
elasticConsts[ 1] = bulkmodulus - 2.0*visFac; // C1122
@@ -529,28 +491,7 @@
PetscLogFlopsNoCheck(6);
- // The code below should probably be in a separate function since it
- // is used more than once. I should also probably cover the possibility
- // that Maxwell time is zero (although this should never happen).
- const double timeFrac = 1.0e-5;
- const int numTerms = 5;
- double dq = 0.0;
- if (maxwelltime < timeFrac*_dt) {
- double fSign = 1.0;
- double factorial = 1.0;
- double fraction = 1.0;
- dq = 1.0;
- for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
- factorial *= iTerm;
- fSign *= -1.0;
- fraction *= _dt/maxwelltime;
- dq += fSign*fraction/factorial;
- } // for
- PetscLogFlopsNoCheck(1 + 7 * numTerms);
- } else {
- dq = maxwelltime*(1.0-exp(-_dt/maxwelltime))/_dt;
- PetscLogFlopsNoCheck(7);
- } // else
+ double dq = ViscoelasticMaxwell::computeVisStrain(_dt, maxwelltime);
const double expFac = exp(-_dt/maxwelltime);
double devStrainTpdt = 0.0;
Added: short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.cc 2008-01-10 18:06:09 UTC (rev 8997)
+++ short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.cc 2008-01-11 16:47:12 UTC (rev 8998)
@@ -0,0 +1,65 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "ViscoelasticMaxwell.hh" // implementation of object methods
+
+#include "petsc.h" // USES PetscLogFlopsNoCheck
+
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Compute viscous strain parameter for a linear Maxwell model.
+double
+pylith::materials::ViscoelasticMaxwell::computeVisStrain(const double dt,
+ const double maxwelltime)
+{ // check parameters and define cutoff values
+ if (maxwelltime <= 0.0)
+ throw std::runtime_error("Maxwell time must be > 0.");
+ const double timeFrac = 1.0e-10;
+ const int numTerms = 5;
+
+ // Compute viscous strain parameter.
+ // The ratio of dt and maxwelltime should never approach timeFrac for any
+ // reasonable computation, but I have put in alternative solutions just in
+ // case.
+ double dq = 0.0;
+ // Use series expansion if dt is very small, since default solution blows
+ // up otherwise.
+ if(dt < timeFrac*maxwelltime) {
+ double fSign = 1.0;
+ double factorial = 1.0;
+ double fraction = 1.0;
+ dq = 1.0;
+ for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
+ factorial *= iTerm;
+ fSign *= -1.0;
+ fraction *= dt/maxwelltime;
+ dq += fSign*fraction/factorial;
+ } // for
+ PetscLogFlopsNoCheck(8*(numTerms-1));
+ // Throw away exponential term if maxwelltime is very small.
+ } else if (maxwelltime < timeFrac*dt) {
+ dq = maxwelltime/dt;
+ PetscLogFlopsNoCheck(1);
+ // Default solution.
+ } else{
+ dq = maxwelltime*(1.0-exp(-dt/maxwelltime))/dt;
+ PetscLogFlopsNoCheck(6);
+ } // else
+
+ return dq;
+} // computeVisStrain
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.hh 2008-01-10 18:06:09 UTC (rev 8997)
+++ short/3D/PyLith/trunk/libsrc/materials/ViscoelasticMaxwell.hh 2008-01-11 16:47:12 UTC (rev 8998)
@@ -0,0 +1,51 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/ViscoelasticMaxwell.hh
+ *
+ * @brief C++ ViscoelasticMaxwell object.
+ *
+ * This class contains a single function that can be used by any
+ * linear Maxwell viscoelastic class.
+ */
+
+#if !defined(pylith_materials_viscoelasticmaxwell_hh)
+#define pylith_materials_viscoelasticmaxwell_hh
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class ViscoelasticMaxwell;
+ } // materials
+
+} // pylith
+
+/// C++ abstract base class for Material object.
+class pylith::materials::ViscoelasticMaxwell
+{ // class ViscoelasticMaxwell
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /** Compute viscous strain parameter.
+ *
+ * @returns Viscous strain parameter.
+ */
+ static double computeVisStrain(const double dt,
+ const double maxwelltime);
+
+}; // class ViscoelasticMaxwell
+
+#endif // pylith_materials_viscoelasticmaxwell_hh
+
+
+// End of file
More information about the cig-commits
mailing list