[cig-commits] r7034 - in short/3D/PyLith/trunk: libsrc
libsrc/feassemble libsrc/materials pylith
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Jun 1 13:13:00 PDT 2007
Author: willic3
Date: 2007-06-01 13:13:00 -0700 (Fri, 01 Jun 2007)
New Revision: 7034
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
short/3D/PyLith/trunk/libsrc/materials/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
short/3D/PyLith/trunk/libsrc/materials/Material.icc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
short/3D/PyLith/trunk/pylith/Makefile.am
Log:
Fixed Maxwell stuff so things compile now, but they haven't been tested.
Added stuff to use needNewJacobian flag.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-06-01 20:13:00 UTC (rev 7034)
@@ -51,6 +51,7 @@
materials/ElasticMaterial.cc \
materials/ElasticPlaneStrain.cc \
materials/ElasticPlaneStress.cc \
+ materials/MaxwellIsotropic3D.cc \
meshio/BinaryIO.cc \
meshio/GMVFile.cc \
meshio/GMVFileAscii.cc \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-06-01 20:13:00 UTC (rev 7034)
@@ -64,6 +64,15 @@
} // material
// ----------------------------------------------------------------------
+// Determine whether we need to recompute the Jacobian.
+bool
+pylith::feassemble::ElasticityImplicit::needNewJacobian(void)
+{ // needNewJacobian
+ _needNewJacobian = _material->needNewJacobian();
+ return _needNewJacobian;
+} // needNewJacobian
+
+// ----------------------------------------------------------------------
// Integrate constributions to residual term (r) for operator.
void
pylith::feassemble::ElasticityImplicit::integrateResidual(
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-06-01 20:13:00 UTC (rev 7034)
@@ -129,6 +129,12 @@
*/
void updateState(const ALE::Obj<real_section_type>& field,
const ALE::Obj<Mesh>& mesh);
+
+ /** Determine whether we need to recompute the Jacobian.
+ *
+ * @returns True if Jacobian needs to be recomputed, false otherwise.
+ */
+ bool needNewJacobian(void);
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
@@ -139,6 +145,11 @@
/// Not implemented
const ElasticityImplicit& operator=(const ElasticityImplicit&);
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ bool _needNewJacobian;
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2007-06-01 20:13:00 UTC (rev 7034)
@@ -26,6 +26,8 @@
ElasticPlaneStrain.icc \
ElasticPlaneStress.hh \
ElasticPlaneStress.icc \
+ MaxwellIsotropic3D.hh \
+ MaxwellIsotropic3D.icc \
Material.hh \
Material.icc
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2007-06-01 20:13:00 UTC (rev 7034)
@@ -30,6 +30,7 @@
// Default constructor.
pylith::materials::Material::Material(void) :
_dt(0.0),
+ _dtm(0.0),
_parameters(0),
_dimension(0),
_needNewJacobian(false),
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2007-06-01 20:13:00 UTC (rev 7034)
@@ -190,6 +190,7 @@
protected :
double _dt; ///< Current time step
+ double _dtm; ///< Previous time step
///< Manager of parameters for physical properties of material
topology::FieldsManager* _parameters;
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.icc 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.icc 2007-06-01 20:13:00 UTC (rev 7034)
@@ -60,6 +60,7 @@
inline
void
pylith::materials::Material::timestep(const double dt) {
+ _dtm = _dt;
_dt = dt;
}
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-06-01 20:13:00 UTC (rev 7034)
@@ -119,8 +119,8 @@
assert(_MaxwellIsotropic3D::numParameters == numParams);
const int numDBValues = dbValues.size();
assert(_MaxwellIsotropic3D::numDBValues == numDBValues);
- for (int i=0; i< numParams; ++i)
- assert(_MaxwellIsotropic3D;;numParamValues[i] == (*paramVals)[i].size());
+ for (int i=0; i < numParams; ++i)
+ assert(_MaxwellIsotropic3D::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_MaxwellIsotropic3D::didDensity];
const double vs = dbValues[_MaxwellIsotropic3D::didVs];
@@ -166,8 +166,9 @@
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::MaxwellIsotropic3D::_calcDensity(double_array* const density,
- const std::vector<double_array&> parameters)
+pylith::materials::MaxwellIsotropic3D::_calcDensity(
+ double_array* const density,
+ const std::vector<double_array>& parameters)
{ // _calcDensity
assert(0 != density);
assert(1 == density->size());
@@ -251,6 +252,7 @@
double devStrainTpdt = 0.0;
double devStrainT = 0.0;
double devStressTpdt = 0.0;
+ double visStrain = 0.0;
for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][iComp] -
@@ -259,7 +261,7 @@
dq*(devStrainTpdt - devStrainT);
devStressTpdt = elasFac*visStrain;
// Later I will want to put in initial stresses.
- (*stress)[iComp] =diag[iComp]*smeanTpdt+devStressTpdt;
+ (*stress)[iComp] =diag[iComp]*meanStressTpdt+devStressTpdt;
} // for
} //else
} // _calcStress
@@ -277,9 +279,9 @@
assert(_MaxwellIsotropic3D::numParameters == parameters.size());
assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
- const double density = parameters[_MaxwellIsotropic3D::pidDensity];
- const double mu = parameters[_MaxwellIsotropic3D::pidMu];
- const double lambda = parameters[_MaxwellIsotropic3D::pidLambda];
+ const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
+ const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
+ const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
const double lambda2mu = lambda + 2.0 * mu;
@@ -354,7 +356,7 @@
// Update state variables.
void
pylith::materials::MaxwellIsotropic3D::_updateState(
- const std::vector<double_array>& parameters,
+ std::vector<double_array>& parameters,
const double_array& totalStrain)
{ // _updateState
assert(_MaxwellIsotropic3D::numParameters == parameters.size());
@@ -384,7 +386,7 @@
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
if (useElasticBehavior()) {
- for (int iComp=0; iComp < tensorSize; ++iComp) {
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
parameters[_MaxwellIsotropic3D::pidStrainT][iComp] = totalStrain[iComp];
parameters[_MaxwellIsotropic3D::pidVisStrain][iComp] =
totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
@@ -417,6 +419,7 @@
const double expFac = exp(-_dt/maxwelltime);
double devStrainTpdt = 0.0;
double devStrainT = 0.0;
+ double visStrain = 0.0;
for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][iComp] -
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2007-06-01 20:13:00 UTC (rev 7034)
@@ -49,6 +49,14 @@
/// Destructor
~MaxwellIsotropic3D(void);
+ /** Get flag indicating whether Jacobian matrix must be reformed for
+ * current state.
+ *
+ * @returns True if Jacobian matrix must be reformed, false otherwise.
+ */
+
+ bool needNewJacobian(void);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -76,7 +84,7 @@
*
* @returns Number of parameters
*/
- int _numParamValues(int_array* numValues) const;
+ void _numParamValues(int_array* numValues) const;
/** Compute parameters from values in spatial database.
*
@@ -137,6 +145,14 @@
const std::vector<double_array>& parameters,
const double_array& totalStrain);
+ /** Update state variables after solve.
+ *
+ * @param parameters Parameters at locations.
+ * @param totalStrain Total strain at locations.
+ */
+ void _updateState(std::vector<double_array>& parameters,
+ const double_array& totalStrain);
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -146,6 +162,7 @@
/// Not implemented
const MaxwellIsotropic3D& operator=(const MaxwellIsotropic3D& m);
+
}; // class MaxwellIsotropic3D
#include "MaxwellIsotropic3D.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2007-06-01 20:13:00 UTC (rev 7034)
@@ -14,5 +14,21 @@
#error "MaxwellIsotropic3D.icc can only be included from MaxwellIsotropic3D.hh"
#endif
+/* Get flag indicating whether Jacobian matrix must be reformed for
+ * current state.
+ */
+inline
+bool
+pylith::materials::MaxwellIsotropic3D::needNewJacobian(void) {
+ // Jacobian should be reformed for the elastic solution, for the first
+ // time step, and any time the step size changes.
+ if (_dt == _dtm)
+ _needNewJacobian = false;
+ else
+ _needNewJacobian = true;
+ return _needNewJacobian;
+}
+
+
// End of file
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2007-06-01 19:08:06 UTC (rev 7033)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2007-06-01 20:13:00 UTC (rev 7034)
@@ -53,6 +53,7 @@
materials/Homogeneous.py \
materials/Material.py \
materials/MaterialsBin.py \
+ materials/MaxwellIsotropic3D.py \
meshio/__init__.py \
meshio/MeshIO.py \
meshio/MeshIOAscii.py \
More information about the cig-commits
mailing list