[cig-commits] r16469 - in short/3D/PyLith/trunk: . libsrc/faults libsrc/feassemble libsrc/friction libsrc/problems modulesrc/problems unittests/libtests/faults unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Mar 30 16:04:49 PDT 2010
Author: brad
Date: 2010-03-30 16:04:47 -0700 (Tue, 30 Mar 2010)
New Revision: 16469
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc
short/3D/PyLith/trunk/libsrc/problems/Explicit.cc
short/3D/PyLith/trunk/libsrc/problems/Explicit.hh
short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
short/3D/PyLith/trunk/libsrc/problems/Implicit.cc
short/3D/PyLith/trunk/libsrc/problems/Implicit.hh
short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh
Log:
Use _calcRateFields to compute rate fields. Removed obsolete needVelocity().
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/TODO 2010-03-30 23:04:47 UTC (rev 16469)
@@ -29,6 +29,9 @@
+ Large deformation formulation
+ Lumped solver
+* valgrind segfaults
+ cases with 1-D mesh - start with meshio
+
SECONDARY PRIORITIES
* Uniform global refinement for tets with faults
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -58,7 +58,6 @@
_jacobian(0),
_ksp(0)
{ // constructor
- _needVelocity = true;
} // constructor
// ----------------------------------------------------------------------
@@ -119,6 +118,8 @@
// Setup fault constitutive model.
assert(0 != _friction);
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
_friction->initialize(*_faultMesh, _quadrature, _fields->get("area"));
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -31,7 +31,7 @@
_gravityField(0),
_logger(0),
_needNewJacobian(true),
- _needVelocity(false),
+ _isJacobianSymmetric(true),
_useSolnIncr(false)
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-03-30 23:04:47 UTC (rev 16469)
@@ -112,7 +112,7 @@
*
* @returns True if integrator needs velocity for computation.
*/
- bool needVelocity(void) const;
+ bool isJacobianSymmetric(void) const;
/** Set flag for setting constraints for total field solution or
* incremental field solution.
@@ -313,7 +313,7 @@
/// True if we need to compute velocity field, false otherwise.
/// Default is false;
- bool _needVelocity;
+ bool _isJacobianSymmetric;
/// Flag indicating whether to set constraints for a total field
/// solution or an incremental field solution
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -34,8 +34,8 @@
template<typename quadrature_type>
inline
bool
-pylith::feassemble::Integrator<quadrature_type>::needVelocity(void) const {
- return _needVelocity;
+pylith::feassemble::Integrator<quadrature_type>::isJacobianSymmetric(void) const {
+ return _isJacobianSymmetric;
} // needsVelocity
// Set flag for setting constraints for total field solution or
Modified: short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -107,7 +107,7 @@
assert(0 != quadrature);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Friction");
+ //logger.stagePush("Friction");
// Get quadrature information
const int numQuadPts = quadrature->numQuadPts();
@@ -362,7 +362,7 @@
_propertiesVertex.resize(_numPropsVertex);
_stateVarsVertex.resize(_numVarsVertex);
- logger.stagePop();
+ //logger.stagePop();
} // initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/problems/Explicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Explicit.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/problems/Explicit.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -34,14 +34,18 @@
} // destructor
// ----------------------------------------------------------------------
-// Compute velocity at time t.
+// Compute velocity and acceleration at time t.
void
-pylith::problems::Explicit::_calcVelocity(void)
-{ // _calcVelocity
+pylith::problems::Explicit::_calcRateFields(void)
+{ // _calcRateFields
assert(0 != _fields);
// vel(t) = (disp(t+dt) - disp(t-dt)) / (2*dt)
// = (dispIncr(t+dt) + disp(t) - disp(t-dt)) / (2*dt)
+ //
+ // acc(t) = (disp(t+dt) - 2*disp(t) + disp(t-dt)) / (dt*dt)
+ // = (dispIncr(t+dt) - disp(t) + disp(t-dt)) / (dt*dt)
+
const double dt = _dt;
topology::Field<topology::Mesh>& dispIncr = _fields->get("dispIncr(t->t+dt)");
@@ -53,6 +57,11 @@
_fields->add("velocity(t)", "velocity");
topology::Field<topology::Mesh>& velocity = _fields->get("velocity(t)");
velocity.cloneSection(dispIncr);
+
+ _fields->add("acceleration(t)", "acceleration");
+ topology::Field<topology::Mesh>& acceleration =
+ _fields->get("acceleration(t)");
+ acceleration.cloneSection(dispIncr);
} // if
// Get sections.
@@ -69,11 +78,16 @@
_fields->get("disp(t-dt)").section();
assert(!dispTmdtSection.isNull());
- double_array velocityVertex(spaceDim);
- topology::Field<topology::Mesh>& velocity = _fields->get("velocity(t)");
- const ALE::Obj<RealSection>& velocitySection = velocity.section();
- assert(!velocitySection.isNull());
+ double_array velVertex(spaceDim);
+ const ALE::Obj<RealSection>& velSection =
+ _fields->get("velocity(t)").section();
+ assert(!velSection.isNull());
+ double_array accVertex(spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ _fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+
// Get mesh vertices.
const ALE::Obj<SieveMesh>& sieveMesh = dispIncr.mesh().sieveMesh();
assert(!sieveMesh.isNull());
@@ -92,14 +106,20 @@
dispTVertex.size());
dispTmdtSection->restrictPoint(*v_iter, &dispTmdtVertex[0],
dispTmdtVertex.size());
- velocityVertex = (dispIncrVertex + dispTVertex - dispTmdtVertex) / (2.0 * dt);
+
+ velVertex = (dispIncrVertex + dispTVertex - dispTmdtVertex) / (2.0 * dt);
+ accVertex = (dispIncrVertex - dispTVertex + dispTmdtVertex) / (dt * dt);
- assert(velocitySection->getFiberDimension(*v_iter) == spaceDim);
- velocitySection->updatePoint(*v_iter, &velocityVertex[0]);
+ assert(velSection->getFiberDimension(*v_iter) == spaceDim);
+ velSection->updatePoint(*v_iter, &velVertex[0]);
+
+ assert(accSection->getFiberDimension(*v_iter) == spaceDim);
+ accSection->updatePoint(*v_iter, &accVertex[0]);
} // for
- PetscLogFlops(vertices->size() * spaceDim);
-} // _calcVelocity
+ PetscLogFlops(vertices->size() * 6*spaceDim);
+} // _calcRateFields
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/problems/Explicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Explicit.hh 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/problems/Explicit.hh 2010-03-30 23:04:47 UTC (rev 16469)
@@ -44,8 +44,8 @@
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
- /// Compute velocity at time t.
- void _calcVelocity(void);
+ /// Compute velocity and acceleration at time t.
+ void _calcRateFields(void);
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -32,7 +32,7 @@
_jacobian(0),
_jacobianLumped(0),
_fields(0),
- _needVelocity(false)
+ _isJacobianSymmetric(false)
{ // constructor
} // constructor
@@ -64,10 +64,10 @@
// ----------------------------------------------------------------------
// Get flag indicating whether we need to compute velocity at time t.
bool
-pylith::problems::Formulation::needVelocity(void) const
-{ // needVelocity
- return _needVelocity;
-} // needVelocity
+pylith::problems::Formulation::isJacobianSymmetric(void) const
+{ // isJacobianSymmetric
+ return _isJacobianSymmetric;
+} // isJacobianSymmetric
// ----------------------------------------------------------------------
// Set integrators over the mesh.
@@ -102,18 +102,18 @@
{ // initialize
// Determine whether we need to compute velocity
- _needVelocity = false;
+ _isJacobianSymmetric = false;
int numIntegrators = _meshIntegrators.size();
for (int i=0; i < numIntegrators; ++i) {
assert(0 != _meshIntegrators[i]);
- if (_meshIntegrators[i]->needVelocity())
- _needVelocity = true;
+ if (_meshIntegrators[i]->isJacobianSymmetric())
+ _isJacobianSymmetric = true;
} // fpr
numIntegrators = _submeshIntegrators.size();
for (int i=0; i < numIntegrators; ++i) {
assert(0 != _submeshIntegrators[i]);
- if (_submeshIntegrators[i]->needVelocity())
- _needVelocity = true;
+ if (_submeshIntegrators[i]->isJacobianSymmetric())
+ _isJacobianSymmetric = true;
} // for
} // initialize
@@ -169,9 +169,8 @@
solution.scatterVectorToSection(*tmpSolutionVec);
} // if
- // Compute velocity if necessary.
- if (_needVelocity)
- _calcVelocity();
+ // Compute rate fields.
+ _calcRateFields();
// Set residual to zero.
topology::Field<topology::Mesh>& residual = _fields->get("residual");
@@ -226,9 +225,8 @@
solution.scatterVectorToSection(*tmpSolutionVec);
} // if
- // Compute velocity if necessary.
- if (_needVelocity)
- _calcVelocity();
+ // Compute rate fields.
+ _calcRateFields();
// Set residual to zero.
topology::Field<topology::Mesh>& residual = _fields->get("residual");
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-03-30 23:04:47 UTC (rev 16469)
@@ -61,11 +61,11 @@
*/
const topology::SolutionFields& fields(void) const;
- /** Get flag indicating whether we need to compute velocity at time t.
+ /** Get flag indicating whether Jacobian is symmetric.
*
- * @returns True if velocity is needed, otherwise false.
+ * @returns True if Jacobian is symmetric, otherwise false.
*/
- bool needVelocity(void) const;
+ bool isJacobianSymmetric(void) const;
/** Set handles to integrators over the mesh.
*
@@ -155,7 +155,7 @@
/// Compute velocity at time t.
virtual
- void _calcVelocity(void) = 0;
+ void _calcRateFields(void) = 0;
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
@@ -172,7 +172,7 @@
///< Integrators over lower-dimensional subdomains of the mesh.
std::vector<IntegratorSubMesh*> _submeshIntegrators;
- bool _needVelocity; ///< Integrator(s) need velocity.
+ bool _isJacobianSymmetric; ///< Is system Jacobian symmetric?
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/problems/Implicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Implicit.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/problems/Implicit.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -36,8 +36,8 @@
// ----------------------------------------------------------------------
// Compute velocity at time t.
void
-pylith::problems::Implicit::_calcVelocity(void)
-{ // _calcVelocity
+pylith::problems::Implicit::_calcRateFields(void)
+{ // _calcRateFields
assert(0 != _fields);
// vel(t) = (disp(t+dt) - disp(t)) / dt
@@ -60,10 +60,10 @@
const ALE::Obj<RealSection>& dispIncrSection = dispIncr.section();
assert(!dispIncrSection.isNull());
- double_array velocityVertex(spaceDim);
- topology::Field<topology::Mesh>& velocity = _fields->get("velocity(t)");
- const ALE::Obj<RealSection>& velocitySection = velocity.section();
- assert(!velocitySection.isNull());
+ double_array velVertex(spaceDim);
+ const ALE::Obj<RealSection>& velSection =
+ _fields->get("velocity(t)").section();
+ assert(!velSection.isNull());
// Get mesh vertices.
const ALE::Obj<SieveMesh>& sieveMesh = dispIncr.mesh().sieveMesh();
@@ -79,14 +79,14 @@
++v_iter) {
dispIncrSection->restrictPoint(*v_iter, &dispIncrVertex[0],
dispIncrVertex.size());
- velocityVertex = dispIncrVertex / dt;
+ velVertex = dispIncrVertex / dt;
- assert(velocitySection->getFiberDimension(*v_iter) == spaceDim);
- velocitySection->updatePoint(*v_iter, &velocityVertex[0]);
+ assert(velSection->getFiberDimension(*v_iter) == spaceDim);
+ velSection->updatePoint(*v_iter, &velVertex[0]);
} // for
PetscLogFlops(vertices->size() * spaceDim);
-} // _calcVelocity
+} // _calcRateFields
// End of file
Modified: short/3D/PyLith/trunk/libsrc/problems/Implicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Implicit.hh 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/libsrc/problems/Implicit.hh 2010-03-30 23:04:47 UTC (rev 16469)
@@ -45,7 +45,7 @@
protected :
/// Compute velocity at time t.
- void _calcVelocity(void);
+ void _calcRateFields(void);
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-03-30 23:04:47 UTC (rev 16469)
@@ -40,11 +40,11 @@
*/
const pylith::topology::SolutionFields& fields(void) const;
- /** Get flag indicating whether we need to compute velocity at time t.
+ /** Get flag indicating whether Jacobian is symmetric.
*
- * @returns True if velocity is needed, otherwise false.
+ * @returns True if Jacobian is symmetric, otherwise false.
*/
- bool needVelocity(void) const;
+ bool isJacobianSymmetric(void) const;
/** Set handles to integrators over the mesh.
*
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -77,16 +77,6 @@
} // testConstructor
// ----------------------------------------------------------------------
-// Test needVelocity().
-void
-pylith::faults::TestFaultCohesiveDyn::testNeedVelocity(void)
-{ // testNeedVelocity
- FaultCohesiveDyn fault;
-
- CPPUNIT_ASSERT_EQUAL(true, fault.needVelocity());
-} // testNeedVelocity
-
-// ----------------------------------------------------------------------
// Test dbInitialTract().
void
pylith::faults::TestFaultCohesiveDyn::testDBInitialTract(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2010-03-30 23:04:47 UTC (rev 16469)
@@ -46,7 +46,6 @@
CPPUNIT_TEST_SUITE( TestFaultCohesiveDyn );
CPPUNIT_TEST( testConstructor );
- CPPUNIT_TEST( testNeedVelocity );
CPPUNIT_TEST( testDBInitialTract );
// Tests in derived classes:
@@ -81,9 +80,6 @@
/// Test constructor.
void testConstructor(void);
- /// Test needVelocity().
- void testNeedVelocity(void);
-
/// Test dbInitialTract().
void testDBInitialTract(void);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc 2010-03-30 23:04:47 UTC (rev 16469)
@@ -49,17 +49,17 @@
} // testStableTimeStep
// ----------------------------------------------------------------------
-// Test needVelocity().
+// Test isJacobianSymmetric().
void
-pylith::feassemble::TestIntegrator::testNeedVelocity(void)
-{ // testNeedVelocity
+pylith::feassemble::TestIntegrator::testIsJacobianSymmetric(void)
+{ // testIsJacobianSymmetric
ElasticityExplicit integrator;
- CPPUNIT_ASSERT_EQUAL(false, integrator.needVelocity());
+ CPPUNIT_ASSERT_EQUAL(true, integrator.isJacobianSymmetric());
- integrator._needVelocity = true;
- CPPUNIT_ASSERT_EQUAL(true, integrator.needVelocity());
-} // testNeedVelocity
+ integrator._isJacobianSymmetric = false;
+ CPPUNIT_ASSERT_EQUAL(false, integrator.isJacobianSymmetric());
+} // testIsJacobianSymmetric
// ----------------------------------------------------------------------
// Test quadrature().
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh 2010-03-30 23:03:29 UTC (rev 16468)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh 2010-03-30 23:04:47 UTC (rev 16469)
@@ -42,7 +42,7 @@
CPPUNIT_TEST( testTimeStep );
CPPUNIT_TEST( testStableTimeStep );
- CPPUNIT_TEST( testNeedVelocity );
+ CPPUNIT_TEST( testIsJacobianSymmetric );
CPPUNIT_TEST( testQuadrature );
CPPUNIT_TEST( testNormalizer );
@@ -65,8 +65,8 @@
/// Test stableTimeStep().
void testStableTimeStep(void);
- /// Test needVelocity().
- void testNeedVelocity(void);
+ /// Test isJacobianSymmetric().
+ void testIsJacobianSymmetric(void);
/// Test quadrature().
void testQuadrature(void);
More information about the CIG-COMMITS
mailing list