[cig-commits] r6970 - in short/3D/PyLith/trunk: . doc doc/uml
examples/twotri3 libsrc libsrc/faults libsrc/feassemble
libsrc/materials libsrc/meshio libsrc/topology
modulesrc/feassemble pylith pylith/faults pylith/feassemble
pylith/problems unittests/libtests/faults
unittests/libtests/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri May 25 13:57:28 PDT 2007
Author: brad
Date: 2007-05-25 13:57:26 -0700 (Fri, 25 May 2007)
New Revision: 6970
Added:
short/3D/PyLith/trunk/doc/uml/
short/3D/PyLith/trunk/doc/uml/PyLith.vpp
short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicit.py
short/3D/PyLith/trunk/pylith/feassemble/ElasticityImplicit.py
short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
Removed:
short/3D/PyLith/trunk/doc/PyLith.vpp
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.icc
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.icc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh
short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py
short/3D/PyLith/trunk/pylith/feassemble/ImplicitElasticity.py
short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py
short/3D/PyLith/trunk/pylith/feassemble/IntegratorImplicit.py
short/3D/PyLith/trunk/pylith/problems/EqDeformation.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestExplicitElasticity.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestImplicitElasticity.py
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/examples/twotri3/implicit.cfg
short/3D/PyLith/trunk/examples/twotri3/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/faults/Fault.cc
short/3D/PyLith/trunk/libsrc/faults/Fault.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/meshio/BinaryIO.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
short/3D/PyLith/trunk/libsrc/topology/FieldsManager.cc
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/pylith/PyLithApp.py
short/3D/PyLith/trunk/pylith/faults/FaultsBin.py
short/3D/PyLith/trunk/pylith/faults/SingleFault.py
short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
short/3D/PyLith/trunk/pylith/feassemble/__init__.py
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/pylith/problems/Problem.py
short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh
short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py
Log:
Propagated reorganization of Integrator and Constraint into Elasticity integrators and faults conditions. Removed superfluous C++ objects (IntegratorExplicit, IntegratorImplicit). Added updateState() to Integrators.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/TODO 2007-05-25 20:57:26 UTC (rev 6970)
@@ -2,44 +2,23 @@
MAIN PRIORITIES (Brad)
======================================================================
-Cleanup
+1. Unit tests for Integrator stuff.
- Add check at beginning of simulation to make sure cell dimension
- matches space dimension. We don't have enough information to
- determine the fault orientation information for a cell in a higher
- dimension coordinate space. We also haven't implemented the
- elasticity stuff for lower dimension cells in a higher
- dimensions. Since Sieve can handle lower order meshes in higher
- dimensions, we probably want to allow it in MeshIO, just not in the
- simulation where we don't support it.
+ b. Integrator
+ iii. C++ unit tests
+ iv. Python unit tests
-1. Switch to more uniform implementation of integrators.
+ e. ElasticityExplicit
+ ii. Python unit tests
- Formulation HASA FieldsManager (solnfields)
+ f. ElasticityImplicit
+ ii. Python unit tests
- Integrators have
- integrateResidual(residual, solnfields, mesh, timeId)
- integrateJacobian(jacobian, solnfields, mesh, timeId)
+ g. Add unit test for IntegratorElasticity::calcTotalStrain
- timeId = {'t+dt', 't', 't-dt'}
+ i. Add in use of useElasticBehavior()
-2. Finish implementing ExplicitElasticity and Explicit
- a. Replace integrateConstant() with integrateResidual()
-
- {f}-[A]{u}, {u} is "guess" (zero for explicit)
-
- a. Double check loops for integrateResidual() and integrateJacobian()
- b. Create unit test (construction of residual term and Jacobian)
-
- c. Add unit test for IntegratorElasticity::calcTotalStrain
-
-2a. Cleanup integrators and boundary conditions.
-
- a. Categorize into "integrators" and "constraints"
- b. Add feature to materials and integrators wherein they indicate
- if the Jacobian needs to be reformed at the current time step.
-
-3. Add dualBasis to Quadrature.
+2. Add dualBasis to Quadrature.
a. Python
ReferenceCell
FIATSimplex
@@ -49,7 +28,7 @@
c. C++ unit tests
d. Python unit tests
-4. Implement faults for kinematic source
+3. Implement faults for kinematic source
a. Creation of cohesive cells
i. Add tests for interpolated meshes.
Double check consistency in ordering of vertices (positive/negative).
@@ -94,9 +73,9 @@
constructor
initialize()
-5. Implement absorbing boundary conditions
+4. Implement absorbing boundary conditions
-6. Create suite of simple full test cases
+5. Create suite of simple full test cases
======================================================================
SECONDARY PRIORITIES
Deleted: short/3D/PyLith/trunk/doc/PyLith.vpp
===================================================================
(Binary files differ)
Copied: short/3D/PyLith/trunk/doc/uml/PyLith.vpp (from rev 6947, short/3D/PyLith/trunk/doc/PyLith.vpp)
===================================================================
(Binary files differ)
Modified: short/3D/PyLith/trunk/examples/twotri3/implicit.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twotri3/implicit.cfg 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/examples/twotri3/implicit.cfg 2007-05-25 20:57:26 UTC (rev 6970)
@@ -4,5 +4,5 @@
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
-[pylithapp.eqdeformation]
+[pylithapp.timedependent]
formulation = pylith.problems.Implicit
Modified: short/3D/PyLith/trunk/examples/twotri3/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twotri3/pylithapp.cfg 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/examples/twotri3/pylithapp.cfg 2007-05-25 20:57:26 UTC (rev 6970)
@@ -27,7 +27,7 @@
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
-[pylithapp.eqdeformation]
+[pylithapp.timedependent]
total_time = 1.0*s
default_dt = 1.0*s
dimension = 2
@@ -35,10 +35,10 @@
# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------
-[pylithapp.eqdeformation.materials]
+[pylithapp.timedependent.materials]
material = pylith.materials.ElasticPlaneStrain
-[pylithapp.eqdeformation.materials.material]
+[pylithapp.timedependent.materials.material]
label = elastic material
id = 1
db.iohandler.filename = matprops.spatialdb
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-05-25 20:57:26 UTC (rev 6970)
@@ -33,12 +33,10 @@
faults/FaultCohesiveDyn.cc \
faults/SlipTimeFn.cc \
feassemble/Constraint.cc \
- feassemble/ExplicitElasticity.cc \
- feassemble/ImplicitElasticity.cc \
+ feassemble/Elasticity.cc \
+ feassemble/ElasticityExplicit.cc \
+ feassemble/ElasticityImplicit.cc \
feassemble/Integrator.cc \
- feassemble/IntegratorElasticity.cc \
- feassemble/IntegratorExplicit.cc \
- feassemble/IntegratorImplicit.cc \
feassemble/Quadrature.cc \
feassemble/Quadrature1D.cc \
feassemble/Quadrature1Din2D.cc \
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -28,13 +28,5 @@
{ // destructor
} // destructor
-// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::faults::Fault::Fault(const Fault& f) :
- _id(f._id),
- _label(f._label)
-{ // copy constructor
-} // copy constructor
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -55,13 +55,6 @@
virtual
~Fault(void);
- /** Create copy of fault.
- *
- * @returns Copy of fault.
- */
- virtual
- Fault* clone(void) const = 0;
-
/** Set identifier of fault.
*
* @param value Fault identifier
@@ -110,16 +103,13 @@
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Copy constructor.
- *
- * @param m Fault to copy
- */
- Fault(const Fault& m);
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
/// Not implemented
+ Fault(const Fault& m);
+
+ /// Not implemented
const Fault& operator=(const Fault& m);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -38,16 +38,6 @@
} // destructor
// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::faults::FaultCohesive::FaultCohesive(const FaultCohesive& f) :
- Fault(f),
- Integrator(f),
- _faultMesh(new ALE::Obj<ALE::Mesh>)
-{ // copy constructor
- *_faultMesh = *f._faultMesh;
-} // copy constructor
-
-// ----------------------------------------------------------------------
// Adjust mesh topology for fault implementation.
void
pylith::faults::FaultCohesive::adjustTopology(const ALE::Obj<ALE::Mesh>& mesh)
@@ -172,16 +162,5 @@
} // for
} // _orient3D
-// ----------------------------------------------------------------------
-// Get size (fiber dimension) of orientation information.
-int
-pylith::faults::FaultCohesive::_orientationSize(void) const
-{ // _orientationSize
- assert(0 != _quadrature);
- const int size = _quadrature->cellDim()*_quadrature->spaceDim();
- return size;
-} // _orientationSize
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -20,7 +20,6 @@
#define pylith_faults_faultcohesive_hh
#include "Fault.hh" // ISA Fault
-#include "pylith/feassemble/Integrator.hh" // ISA Integrator
#include "pylith/utils/sievefwd.hh" // HOLDSA PETSc Mesh
#include "pylith/utils/petscfwd.h" // USES PetscMat
@@ -35,8 +34,7 @@
/// C++ abstract base class for a fault surface implemented with
/// cohesive elements.
-class pylith::faults::FaultCohesive : public Fault,
- public feassemble::Integrator
+class pylith::faults::FaultCohesive : public Fault
{ // class FaultCohesive
friend class TestFaultCohesive; // unit testing
@@ -56,39 +54,6 @@
*/
void adjustTopology(const ALE::Obj<ALE::Mesh>& mesh);
- /** Integrate contribution of cohesive cells to residual term.
- *
- * @param residual Residual field (output)
- * @param disp Displacement field at time t
- * @param mesh Finite-element mesh
- */
- virtual
- void integrateResidual(const ALE::Obj<real_section_type>& residual,
- const ALE::Obj<real_section_type>& disp,
- const ALE::Obj<Mesh>& mesh) = 0;
-
- /** Compute Jacobian matrix (A) associated with operator.
- *
- * @param mat Sparse matrix
- * @param disp Displacement field
- * @param mesh Finite-element mesh
- */
- virtual
- void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh) = 0;
-
- /** Set field.
- *
- * @param t Current time
- * @param disp Displacement field
- * @param mesh Finite-element mesh
- */
- virtual
- void setField(const double t,
- const ALE::Obj<real_section_type>& disp,
- const ALE::Obj<Mesh>& mesh) = 0;
-
// PROTECTED TYPEDEFS /////////////////////////////////////////////////
protected :
@@ -102,12 +67,6 @@
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Copy constructor.
- *
- * @param m Fault to copy
- */
- FaultCohesive(const FaultCohesive& m);
-
/** Cohesive cells use Lagrange multiplier constraints?
*
* @returns True if implementation using Lagrange multiplier
@@ -191,16 +150,13 @@
const double_array& upDir,
const int numLocs);
- /** Get size (fiber dimension) of orientation information.
- *
- * @returns Size of orientation information.
- */
- int _orientationSize(void) const;
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
/// Not implemented
+ FaultCohesive(const FaultCohesive& m);
+
+ /// Not implemented
const FaultCohesive& operator=(const FaultCohesive& m);
// PROTECTED MEMBERS ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -15,6 +15,7 @@
#include "FaultCohesiveDyn.hh" // implementation of object methods
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "pylith/utils/array.hh" // USES double_array
#include <assert.h> // USES assert()
@@ -34,13 +35,6 @@
} // destructor
// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(const FaultCohesiveDyn& f) :
- FaultCohesive(f)
-{ // copy constructor
-} // copy constructor
-
-// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
void
pylith::faults::FaultCohesiveDyn::initialize(const ALE::Obj<ALE::Mesh>& mesh,
@@ -50,6 +44,7 @@
assert(0 != _quadrature);
assert(0 != _faultMesh);
assert(!_faultMesh->isNull());
+ assert(0 != cs);
if (3 != upDir.size())
throw std::runtime_error("Up direction for fault orientation must be "
@@ -60,7 +55,9 @@
ALE::Obj<real_section_type> orientation =
new real_section_type((*_faultMesh)->comm(), (*_faultMesh)->debug());
assert(!orientation.isNull());
- const int orientationSize = _orientationSize();
+ const int cellDim = (*_faultMesh)->getDimension();
+ const int spaceDim = cs->spaceDim();
+ const int orientationSize = cellDim*spaceDim;
orientation->setFiberDimension((*_faultMesh)->depthStratum(0),
orientationSize);
(*_faultMesh)->allocate(orientation);
@@ -205,7 +202,7 @@
void
pylith::faults::FaultCohesiveDyn::integrateResidual(
const ALE::Obj<real_section_type>& residual,
- const ALE::Obj<real_section_type>& disp,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
#if 0
@@ -255,20 +252,10 @@
void
pylith::faults::FaultCohesiveDyn::integrateJacobian(
PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateJacobian
} // integrateJacobian
-// ----------------------------------------------------------------------
-// Set field.
-void
-pylith::faults::FaultCohesiveDyn::setField(
- const double t,
- const ALE::Obj<real_section_type>& disp,
- const ALE::Obj<Mesh>& mesh)
-{ // setField
-} // setField
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -23,7 +23,8 @@
#if !defined(pylith_faults_faultcohesivedyn_hh)
#define pylith_faults_faultcohesivedyn_hh
-#include "FaultCohesive.hh"
+#include "FaultCohesive.hh" // ISA FaultCohesive
+#include "pylith/feassemble/Integrator.hh" // ISA Constraint
/// Namespace for pylith package
namespace pylith {
@@ -35,7 +36,8 @@
/// @brief C++ implementation for a fault surface with spontaneous
/// (dynamic) slip implemented with cohesive elements.
-class pylith::faults::FaultCohesiveDyn : public FaultCohesive
+class pylith::faults::FaultCohesiveDyn : public FaultCohesive,
+ public feassemble::Integrator
{ // class FaultCohesiveDyn
friend class TestFaultCohesiveDyn; // unit testing
@@ -49,12 +51,6 @@
virtual
~FaultCohesiveDyn(void);
- /** Create copy of fault.
- *
- * @returns Copy of fault.
- */
- Fault* clone(void) const;
-
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
*
@@ -71,42 +67,26 @@
/** Integrate contribution of cohesive cells to residual term.
*
* @param residual Residual field (output)
- * @param disp Displacement field at time t
+ * @param fields Solution fields
* @param mesh Finite-element mesh
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
- const ALE::Obj<real_section_type>& disp,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
/** Compute Jacobian matrix (A) associated with operator.
*
* @param mat Sparse matrix
- * @param disp Displacement field
+ * @param fields Solution fields
* @param mesh Finite-element mesh
*/
void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
- /** Set field.
- *
- * @param t Current time
- * @param disp Displacement field
- * @param mesh Finite-element mesh
- */
- void setField(const double t,
- const ALE::Obj<real_section_type>& disp,
- const ALE::Obj<Mesh>& mesh);
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Copy constructor.
- *
- * @param m Fault to copy
- */
- FaultCohesiveDyn(const FaultCohesiveDyn& m);
-
/** Cohesive cells use Lagrange multiplier constraints?
*
* @returns True if implementation using Lagrange multiplier
@@ -118,6 +98,9 @@
private :
/// Not implemented
+ FaultCohesiveDyn(const FaultCohesiveDyn& m);
+
+ /// Not implemented
const FaultCohesiveDyn& operator=(const FaultCohesiveDyn& m);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -14,13 +14,6 @@
#error "FaultCohesiveDyn.icc can only be included from FaultCohesiveDyn.hh"
#endif
-// Create copy of fault.
-inline
-pylith::faults::Fault*
-pylith::faults::FaultCohesiveDyn::clone(void) const {
- return new FaultCohesiveDyn(*this);
-} // clone
-
// Cohesive cells use Lagrange multiplier constraints?
inline
bool
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -17,8 +17,11 @@
#include "EqKinSrc.hh" // USES EqKinSrc
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "pylith/utils/array.hh" // USES double_array
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
#include <assert.h> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -38,16 +41,6 @@
} // destructor
// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::faults::FaultCohesiveKin::FaultCohesiveKin(const FaultCohesiveKin& f) :
- FaultCohesive(f),
- _eqsrc(0)
-{ // copy constructor
- if (0 != f._eqsrc)
- _eqsrc = f._eqsrc->clone();
-} // copy constructor
-
-// ----------------------------------------------------------------------
// Set kinematic earthquake source.
void
pylith::faults::FaultCohesiveKin::eqsrc(EqKinSrc* src)
@@ -75,7 +68,9 @@
ALE::Obj<real_section_type> orientation =
new real_section_type((*_faultMesh)->comm(), (*_faultMesh)->debug());
assert(!orientation.isNull());
- const int orientationSize = _orientationSize();
+ const int cellDim = (*_faultMesh)->getDimension();
+ const int spaceDim = cs->spaceDim();
+ const int orientationSize = cellDim*spaceDim;
orientation->setFiberDimension((*_faultMesh)->depthStratum(0),
orientationSize);
(*_faultMesh)->allocate(orientation);
@@ -86,8 +81,8 @@
assert(!coordinates.isNull());
// Set orientation method
- const int cellDim = _quadrature->cellDim();
- const int spaceDim = _quadrature->spaceDim();
+ assert(cellDim == _quadrature->cellDim());
+ assert(spaceDim == _quadrature->spaceDim());
orient_fn_type orientFn;
switch (cellDim)
{ // switch
@@ -220,17 +215,29 @@
void
pylith::faults::FaultCohesiveKin::integrateResidual(
const ALE::Obj<real_section_type>& residual,
- const ALE::Obj<real_section_type>& disp,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
+ assert(!residual.isNull());
+ assert(0 != fields);
+ assert(!mesh.isNull());
+ assert(0 != _faultMesh);
+ assert(!_faultMesh->isNull());
+
// Subtract constraint forces (which are disp at the constraint
// DOF) to residual; contributions are at DOF of normal vertices (i and j)
+ // Get cell information
const ALE::Obj<Mesh::label_sequence>& cells =
(*_faultMesh)->heightStratum(0);
+ assert(!cells.isNull());
const Mesh::label_sequence::iterator cBegin = cells->begin();
const Mesh::label_sequence::iterator cEnd = cells->end();
+ // Get section information
+ const ALE::Obj<real_section_type>& disp = fields->getHistoryItem(1);
+ assert(!disp.isNull());
+
// Allocate vector for cell values (if necessary)
_initCellVector();
@@ -268,24 +275,37 @@
void
pylith::faults::FaultCohesiveKin::integrateJacobian(
PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateJacobian
+ assert(0 != mat);
+ assert(0 != fields);
+ assert(!mesh.isNull());
+ assert(0 != _faultMesh);
+ assert(!_faultMesh->isNull());
// Add constraint information to Jacobian matrix; these are the
// direction cosines. Entries are associated with vertices ik, jk,
// ki, and kj.
+ // Get cell informatino
const ALE::Obj<Mesh::label_sequence>& cells =
(*_faultMesh)->heightStratum(0);
const Mesh::label_sequence::iterator cBegin = cells->begin();
const Mesh::label_sequence::iterator cEnd = cells->end();
+ // Get section information
+ const ALE::Obj<real_section_type>& disp = fields->getHistoryItem(1);
+ assert(!disp.isNull());
+
+ const int cellDim = _quadrature->cellDim();
+ const int spaceDim = _quadrature->spaceDim();
+ const int orientationSize = cellDim*spaceDim;
+
// Allocate matrix for cell values (if necessary)
_initCellMatrix();
const int numVertices = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
for (Mesh::label_sequence::iterator c_iter=cBegin;
c_iter != cEnd;
++c_iter) {
@@ -297,7 +317,6 @@
const int numConstraintVert = numVertices / 3;
assert(numVertices == numConstraintVert * 3);
- const int orientationSize = _orientationSize();
for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
// Get orientations at constraint vertex
const real_section_type::value_type* constraintOrient =
@@ -338,8 +357,8 @@
// Assemble cell contribution into PETSc Matrix
const ALE::Obj<Mesh::order_type>& globalOrder =
- mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
- err = updateOperator(*mat, mesh, dispT, globalOrder,
+ mesh->getFactory()->getGlobalOrder(mesh, "default", disp);
+ err = updateOperator(*mat, mesh, disp, globalOrder,
*c_iter, _cellMatrix, ADD_VALUES);
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
@@ -347,6 +366,26 @@
} // integrateJacobian
// ----------------------------------------------------------------------
+// Set number of degrees of freedom that are constrained at points
+void
+pylith::faults::FaultCohesiveKin::setConstraintSizes(
+ const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh)
+{ // setConstraintSizes
+ throw std::logic_error("FaultCohesiveKin::setConstraintSizes() not implemented.");
+} // setConstraintSizes
+
+// ----------------------------------------------------------------------
+// Set which degrees of freedom are constrained at points in field.
+void
+pylith::faults::FaultCohesiveKin::setConstraints(
+ const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh)
+{ // setConstraints
+ throw std::logic_error("FaultCohesiveKin::setConstraints() not implemented.");
+} // setConstraints
+
+// ----------------------------------------------------------------------
// Set field.
void
pylith::faults::FaultCohesiveKin::setField(
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -29,7 +29,9 @@
#if !defined(pylith_faults_faultcohesivekin_hh)
#define pylith_faults_faultcohesivekin_hh
-#include "FaultCohesive.hh"
+#include "FaultCohesive.hh" // ISA FaultCohesive
+#include "pylith/feassemble/Constraint.hh" // ISA Constraint
+#include "pylith/feassemble/Integrator.hh" // ISA Integrator
/// Namespace for pylith package
namespace pylith {
@@ -43,7 +45,9 @@
/// C++ implementation for a fault surface with kinematic (prescribed)
/// slip implemented with cohesive elements.
-class pylith::faults::FaultCohesiveKin : public FaultCohesive
+class pylith::faults::FaultCohesiveKin : public FaultCohesive,
+ public feassemble::Integrator,
+ public feassemble::Constraint
{ // class FaultCohesiveKin
friend class TestFaultCohesiveKin; // unit testing
@@ -57,12 +61,6 @@
virtual
~FaultCohesiveKin(void);
- /** Create copy of fault.
- *
- * @returns Copy of fault.
- */
- Fault* clone(void) const;
-
/** Set kinematic earthquake source.
*
* @param src Kinematic earthquake source.
@@ -82,31 +80,49 @@
const spatialdata::geocoords::CoordSys* cs,
const double_array& upDir);
- /** Integrate contribution of cohesive cells to residual term.
+ /** Integrate contributions to residual term (r) for operator.
*
- * @param residual Residual field (output)
- * @param disp Displacement field at time t
+ * @param residual Field containing values for residual
+ * @param fields Solution fields
* @param mesh Finite-element mesh
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
- const ALE::Obj<real_section_type>& disp,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
- /** Compute Jacobian matrix (A) associated with operator.
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
*
* @param mat Sparse matrix
- * @param disp Displacement field
+ * @param fields Solution fields
* @param mesh Finite-element mesh
*/
void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
-
+
+ /** Set number of degrees of freedom that are constrained at points
+ * in field.
+ *
+ * @param field Solution field
+ * @param mesh PETSc mesh
+ */
+ void setConstraintSizes(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh);
+
+ /** Set which degrees of freedom are constrained at points in field.
+ *
+ * @param field Solution field
+ * @param mesh PETSc mesh
+ */
+ void setConstraints(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh);
+
/** Set field.
*
- * @param t Current time
- * @param disp Displacement field
- * @param mesh Finite-element mesh
+ * @param t Current time.
+ * @param disp Displacement field at time t.
+ * @param mesh Finite-element mesh.
*/
void setField(const double t,
const ALE::Obj<real_section_type>& disp,
@@ -115,12 +131,6 @@
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Copy constructor.
- *
- * @param m Fault to copy
- */
- FaultCohesiveKin(const FaultCohesiveKin& m);
-
/** Cohesive cells use Lagrange multiplier constraints?
*
* @returns True if implementation using Lagrange multiplier
@@ -132,6 +142,9 @@
private :
/// Not implemented
+ FaultCohesiveKin(const FaultCohesiveKin& m);
+
+ /// Not implemented
const FaultCohesiveKin& operator=(const FaultCohesiveKin& m);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -14,13 +14,6 @@
#error "FaultCohesiveKin.icc can only be included from FaultCohesiveKin.hh"
#endif
-// Create copy of fault.
-inline
-pylith::faults::Fault*
-pylith::faults::FaultCohesiveKin::clone(void) const {
- return new FaultCohesiveKin(*this);
-} // clone
-
// Cohesive cells use Lagrange multiplier constraints?
inline
bool
Copied: short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc (from rev 6964, short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,113 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "Elasticity.hh" // implementation of class methods
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::logic_error
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::Elasticity::calcTotalStrain1D(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int numBasis)
+{ // calcTotalStrain1D
+ assert(0 != strain);
+ assert(0 != disp);
+
+ const int dimension = 1;
+ const int numQuadPts = strain->size();
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(1 == (*strain)[iQuad].size());
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ (*strain)[iQuad][0] +=
+ basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
+ } // for
+} // calcTotalStrain1D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::Elasticity::calcTotalStrain2D(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int numBasis)
+{ // calcTotalStrain2D
+ assert(0 != strain);
+ assert(0 != disp);
+
+ const int dimension = 2;
+ const int numQuadPts = strain->size();
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(3 == (*strain)[iQuad].size());
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ (*strain)[iQuad][0] +=
+ basisDeriv[iQ+iBasis ] * disp[iBasis ];
+ (*strain)[iQuad][1] +=
+ basisDeriv[iQ+iBasis+1] * disp[iBasis+1];
+ (*strain)[iQuad][2] +=
+ 0.5 * (basisDeriv[iQ+iBasis+1] * disp[iBasis ] +
+ basisDeriv[iQ+iBasis ] * disp[iBasis+1]);
+ } // for
+ } // for
+} // calcTotalStrain2D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::Elasticity::calcTotalStrain3D(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int numBasis)
+{ // calcTotalStrain3D
+ assert(0 != strain);
+ assert(0 != disp);
+
+ const int dimension = 3;
+ const int numQuadPts = strain->size();
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(6 == (*strain)[iQuad].size());
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ (*strain)[iQuad][0] +=
+ basisDeriv[iQ+iBasis ] * disp[iBasis ];
+ (*strain)[iQuad][1] +=
+ basisDeriv[iQ+iBasis+1] * disp[iBasis+1];
+ (*strain)[iQuad][2] +=
+ basisDeriv[iQ+iBasis+2] * disp[iBasis+2];
+ (*strain)[iQuad][3] +=
+ 0.5 * (basisDeriv[iQ+iBasis+1] * disp[iBasis ] +
+ basisDeriv[iQ+iBasis ] * disp[iBasis+1]);
+ (*strain)[iQuad][4] +=
+ 0.5 * (basisDeriv[iQ+iBasis+2] * disp[iBasis+1] +
+ basisDeriv[iQ+iBasis+1] * disp[iBasis+2]);
+ (*strain)[iQuad][5] +=
+ 0.5 * (basisDeriv[iQ+iBasis+2] * disp[iBasis ] +
+ basisDeriv[iQ+iBasis ] * disp[iBasis+2]);
+ } // for
+ } // for
+} // calcTotalStrain3D
+
+
+// End of file
Copied: short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh (from rev 6964, short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,86 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/Elasticity.hh
+ *
+ * @brief C++ Utility class for general operations associated with
+ * integrating finite-element actions for solid elasticity.
+ */
+
+#if !defined(pylith_feassemble_elasticity_hh)
+#define pylith_feassemble_elasticity_hh
+
+#include "pylith/utils/array.hh" // USES double_array, std::vector
+
+namespace pylith {
+ namespace feassemble {
+ class Elasticity;
+ class TestElasticity;
+ } // feassemble
+} // pylith
+
+class pylith::feassemble::Elasticity
+{ // Elasticity
+ friend class TestElasticity; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param dimension Dimension of cell.
+ * @param numBasis Number of basis functions for cell.
+ */
+
+ static
+ void calcTotalStrain1D(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int numBasis);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ */
+
+ static
+ void calcTotalStrain2D(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int numBasis);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ */
+
+ static
+ void calcTotalStrain3D(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int numBasis);
+
+}; // Elasticity
+
+#endif // pylith_feassemble_elasticity_hh
+
+// End of file
Copied: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc (from rev 6964, short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,454 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ElasticityExplicit.hh" // implementation of class methods
+
+#include "Quadrature.hh" // USES Quadrature
+#include "Elasticity.hh" // USES Elasticity
+
+#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "petscmat.h" // USES PetscMat
+#include "spatialdata/spatialdb/SpatialDB.hh"
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::ElasticityExplicit::ElasticityExplicit(void) :
+ _dtm1(-1.0),
+ _material(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElasticityExplicit::~ElasticityExplicit(void)
+{ // destructor
+ delete _material; _material = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::ElasticityExplicit::timeStep(const double dt)
+{ // timeStep
+ if (_dt != -1.0)
+ _dtm1 = _dt;
+ else
+ _dtm1 = dt;
+ _dt = dt;
+ assert(_dt == _dtm1); // For now, don't allow variable time step
+} // timeStep
+
+// ----------------------------------------------------------------------
+// Set material.
+void
+pylith::feassemble::ElasticityExplicit::material(
+ const materials::ElasticMaterial* m)
+{ // material
+ delete _material; _material = (0 != m) ? m->clone() : 0;
+} // material
+
+// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityExplicit::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateResidual
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(!residual.isNull());
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ // Get cell information
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<real_section_type>& dispTpdt = fields->getHistoryItem(0);
+ const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
+ const ALE::Obj<real_section_type>& dispTmdt = fields->getHistoryItem(2);
+ assert(!dispTpdt.isNull());
+ assert(!dispT.isNull());
+ assert(!dispTmdt.isNull());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ /** :TODO:
+ *
+ * If cellDim and spaceDim are different, we need to transform
+ * displacements into cellDim, compute action, and transform result
+ * back into spaceDim. We get this information from the Jacobian and
+ * inverse of the Jacobian.
+ */
+ if (cellDim != spaceDim)
+ throw std::logic_error("Not implemented yet.");
+
+ // Allocate vector for cell values (if necessary)
+ _initCellVector();
+
+ // Allocate vector for total strain
+ int tensorSize = 0;
+ if (1 == cellDim)
+ tensorSize = 1;
+ else if (2 == cellDim)
+ tensorSize = 3;
+ else if (3 == cellDim)
+ tensorSize = 6;
+ else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'." << std::endl;
+ assert(0);
+ } // else
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+
+ // Set cell data in material
+ _material->initCellData(*c_iter, numQuadPts);
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input fields to cell
+ const real_section_type::value_type* dispTpdtCell =
+ mesh->restrict(dispTpdt, *c_iter);
+ const real_section_type::value_type* dispTCell =
+ mesh->restrict(dispT, *c_iter);
+ const real_section_type::value_type* dispTmdtCell =
+ mesh->restrict(dispTmdt, *c_iter);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Compute action for inertial terms
+ const std::vector<double_array>& density =
+ _material->calcDensity();
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt =
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basis[iQuad*numBasis+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQuad*numBasis+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis*spaceDim+iDim] +=
+ valIJ * (- dispTpdtCell[jBasis*spaceDim+iDim]
+ + 2.0 * dispTCell[jBasis*spaceDim+iDim]
+ - dispTmdtCell[jBasis*spaceDim+iDim]);
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Compute action for elastic terms
+ if (1 == cellDim) {
+ // Compute stresses
+ Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
+ // Compute elastic action
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis*cellDim ];
+ _cellVector[iBlock ] -= N1*s11;
+ } // for
+ } // for
+ PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*5));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ } else if (2 == cellDim) {
+ // Compute stresses
+ Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
+ // Compute elastic action
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s12 = stress[iQuad][2];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
+ const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ _cellVector[iBlock ] -= N1*s11 + N2*s12;
+ _cellVector[iBlock+1] -= N1*s12 + N2*s22;
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ } else if (3 == cellDim) {
+ // Compute stresses
+ Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
+ // Compute elastic action
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s33 = stress[iQuad][2];
+ const double s12 = stress[iQuad][3];
+ const double s23 = stress[iQuad][4];
+ const double s13 = stress[iQuad][5];
+
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBasis*cellDim+0];
+ const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ const double N3 = wt*basisDeriv[iQ+iBasis*cellDim+2];
+ _cellVector[iBlock ] -= N1*s11 + N2*s12 + N3*s13;
+ _cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
+ _cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'."
+ << std::endl;
+ assert(0);
+ } // if/else
+ // Assemble cell contribution into field
+ mesh->updateAdd(residual, *c_iter, _cellVector);
+ } // for
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicit::integrateJacobian(
+ PetscMat* mat,
+ topology::FieldsManager* fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != mat);
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ // Get cell information
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
+ assert(!dispT.isNull());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // Allocate vector for cell values (if necessary)
+ _initCellMatrix();
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+
+ // Set cell data in material
+ _material->initCellData(*c_iter, numQuadPts);
+
+ // Reset element vector to zero
+ _resetCellMatrix();
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Get material physical properties at quadrature points for this cell
+ const std::vector<double_array>& density = _material->calcDensity();
+
+ // Compute Jacobian for inertial terms
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt =
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const int iBlock = (iBasis * spaceDim + iDim) * (spaceDim * numBasis);
+ const int jBlock = (jBasis * spaceDim + iDim);
+ _cellMatrix[iBlock+jBlock] += valIJ;
+ } // for
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into PETSc Matrix
+ const ALE::Obj<Mesh::order_type>& globalOrder =
+ mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
+
+ err = updateOperator(*mat, mesh, dispT, globalOrder,
+ *c_iter, _cellMatrix, ADD_VALUES);
+ if (err)
+ throw std::runtime_error("Update to PETSc Mat failed.");
+ } // for
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Update state variables as needed.
+void
+pylith::feassemble::ElasticityExplicit::updateState(
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // updateState
+ assert(0 != _material);
+ assert(!disp.isNull());
+
+ // No need to update state if using elastic behavior
+ if (_material->useElasticBehavior())
+ return;
+
+ // Get cell information
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ // Allocate vector for total strain
+ int tensorSize = 0;
+ if (1 == cellDim)
+ tensorSize = 1;
+ else if (2 == cellDim)
+ tensorSize = 3;
+ else if (3 == cellDim)
+ tensorSize = 6;
+ else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'." << std::endl;
+ assert(0);
+ } // else
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+
+ // Set cell data in material
+ _material->initCellData(*c_iter, numQuadPts);
+
+ // Restrict input fields to cell
+ const real_section_type::value_type* dispCell =
+ mesh->restrict(disp, *c_iter);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ // Compute action for elastic terms
+ if (1 == cellDim) {
+ Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispCell, numBasis);
+ _material->updateState(totalStrain);
+ } else if (2 == cellDim) {
+ Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispCell, numBasis);
+ _material->updateState(totalStrain);
+ } else if (3 == cellDim) {
+ // Compute stresses
+ Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
+ dispCell, numBasis);
+ _material->updateState(totalStrain);
+ } else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'."
+ << std::endl;
+ assert(0);
+ } // if/else
+ } // for
+} // updateState
+
+
+// End of file
Copied: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh (from rev 6964, short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,156 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/ElasticityExplicit.hh
+ *
+ * @brief Explicit time integration of dynamic elasticity equation
+ * using finite-elements.
+ *
+ * Note: This object operates on a single finite-element family, which
+ * is defined by the quadrature and a database of material property
+ * parameters.
+ *
+ * Computes contributions to terms A and r in
+ *
+ * A(t) u(t+dt) = b(u(t), u(t-dt)),
+ *
+ * r = b - A u0(t+dt)
+ *
+ * where A(t) is a sparse matrix or vector, u(t+dt) is the field we
+ * want to compute at time t+dt, b is a vector that depends on the
+ * field at time t and t-dt, and u0 is zero at unknown DOF and set to
+ * the constrained values at known DOF.
+ *
+ * Contributions from elasticity include the intertial and stiffness
+ * terms, so this object computes the following portions of A and r:
+ *
+ * A = 1/(dt*dt) [M]
+ *
+ * r = (1/(dt*dt) [M])(- {u(t+dt)} + 2/(dt*dt){u(t)} - {u(t-dt)}) - [K]{u(t)}
+ *
+ * Translational inertia.
+ * - Residual action over cell
+ * \f[
+ * \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
+ * \f]
+ * - Jacobian action over cell
+ * \f[
+ * \int_{V^e} (\rho N^q N^q)_i \, dV
+ * \f]
+ *
+ * See governing equations section of user manual for more
+ * information.
+ */
+
+#if !defined(pylith_feassemble_elasticityexplicit_hh)
+#define pylith_feassemble_elasticityexplicit_hh
+
+#include "Integrator.hh" // ISA Integrator
+
+namespace pylith {
+ namespace feassemble {
+ class ElasticityExplicit;
+ class TestElasticityExplicit;
+ } // feassemble
+
+ namespace materials {
+ class ElasticMaterial;
+ } // feassemble
+} // pylith
+
+namespace spatialdata {
+ namespace spatialdb {
+ class SpatialDB; // USES SpatialDB
+ } // spatialdb
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+
+class pylith::feassemble::ElasticityExplicit : public Integrator
+{ // ElasticityExplicit
+ friend class TestElasticityExplicit; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElasticityExplicit(void);
+
+ /// Destructor
+ ~ElasticityExplicit(void);
+
+ /** Set material.
+ *
+ * @param m Elastic material.
+ */
+ void material(const materials::ElasticMaterial* m);
+
+ /** Set time step for advancing from time t to time t+dt.
+ *
+ * @param dt Time step
+ */
+ void timeStep(const double dt);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param mat Sparse matrix
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void integrateJacobian(PetscMat* mat,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Update state variables as needed.
+ *
+ * @param field Current solution field.
+ * @param mesh Finite-element mesh
+ */
+ void updateState(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<Mesh>& mesh);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ ElasticityExplicit(const ElasticityExplicit& i);
+
+ /// Not implemented
+ const ElasticityExplicit& operator=(const ElasticityExplicit&);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double _dtm1; ///< Time step for t-dt1 -> t
+
+ /// Elastic material associated with integrator
+ materials::ElasticMaterial* _material;
+
+}; // ElasticityExplicit
+
+#endif // pylith_feassemble_elasticityexplicit_hh
+
+
+// End of file
Copied: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc (from rev 6964, short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,602 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ElasticityImplicit.hh" // implementation of class methods
+
+#include "Quadrature.hh" // USES Quadrature
+#include "Elasticity.hh" // USES Elasticity
+
+#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "petscmat.h" // USES PetscMat
+#include "spatialdata/spatialdb/SpatialDB.hh"
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
+ _dtm1(-1.0),
+ _material(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElasticityImplicit::~ElasticityImplicit(void)
+{ // destructor
+ delete _material; _material = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::ElasticityImplicit::timeStep(const double dt)
+{ // timeStep
+ if (_dt != -1.0)
+ _dtm1 = _dt;
+ else
+ _dtm1 = dt;
+ _dt = dt;
+ assert(_dt == _dtm1); // For now, don't allow variable time step
+} // timeStep
+
+// ----------------------------------------------------------------------
+// Set material.
+void
+pylith::feassemble::ElasticityImplicit::material(
+ const materials::ElasticMaterial* m)
+{ // material
+ delete _material; _material = (0 != m) ? m->clone() : 0;
+} // material
+
+// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityImplicit::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateResidual
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(!residual.isNull());
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ // Get cell information
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
+ assert(!dispT.isNull());
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ // Allocate vector for cell values (if necessary)
+ _initCellVector();
+
+ // Allocate vector for total strain
+ int tensorSize = 0;
+ if (1 == cellDim)
+ tensorSize = 1;
+ else if (2 == cellDim)
+ tensorSize = 3;
+ else if (3 == cellDim)
+ tensorSize = 6;
+ else
+ throw std::logic_error("Tensor size not implemented for given cellDim.");
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+ // Loop over cells
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+
+ // Set cell data in material
+ _material->initCellData(*c_iter, numQuadPts);
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input fields to cell
+ const real_section_type::value_type* dispTCell =
+ mesh->restrict(dispT, *c_iter);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ if (cellDim != spaceDim)
+ throw std::logic_error("Not implemented yet.");
+
+ /* Comment out gravity section for now, until we figure out how to deal
+ with gravity vector.
+ // Get density at quadrature points for this cell
+ const std::vector<double_array>& density = _material->calcDensity();
+
+ // Compute action for element body forces
+ if (!grav.isNull()) {
+ const real_section_type::value_type* gravCell =
+ mesh->restrict(grav, cell);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double valI = wt*basis[iQ+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ _cellVector[iBlock+iDim] += valI*gravCell[iDim];
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(2+numBasis*(2+2*spaceDim)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } // if
+ */
+
+ // Compute B(transpose) * sigma, first computing strains
+ if (1 == cellDim) {
+ // Compute total strains and then use these to compute stresses
+ Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
+ // Compute elastic action
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
+ _cellVector[iBlock ] -= N1*s11;
+ } // for
+ } // for
+ PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*5));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ } else if (2 == cellDim) {
+ // Compute total strains and then use these to compute stresses
+ Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
+ // Compute elastic action
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s12 = stress[iQuad][2];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
+ const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ _cellVector[iBlock ] -= N1*s11 + N2*s12;
+ _cellVector[iBlock+1] -= N1*s12 + N2*s22;
+ } // for
+ } // for
+ PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ } else if (3 == cellDim) {
+ // Compute total strains and then use these to compute stresses
+ Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
+ // Compute elastic action
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s33 = stress[iQuad][2];
+ const double s12 = stress[iQuad][3];
+ const double s23 = stress[iQuad][4];
+ const double s13 = stress[iQuad][5];
+
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBasis*cellDim+0];
+ const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ const double N3 = wt*basisDeriv[iQ+iBasis*cellDim+2];
+ _cellVector[iBlock ] -= N1*s11 + N2*s12 + N3*s13;
+ _cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
+ _cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
+ } // for
+ } // for
+ PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'."
+ << std::endl;
+ assert(0);
+ } // if/else
+
+ // Assemble cell contribution into field
+ mesh->updateAdd(residual, *c_iter, _cellVector);
+ } // for
+} // integrateResidual
+
+
+// ----------------------------------------------------------------------
+// Compute stiffness matrix.
+void
+pylith::feassemble::ElasticityImplicit::integrateJacobian(
+ PetscMat* mat,
+ topology::FieldsManager* fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != mat);
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ PetscErrorCode err = 0;
+
+ // Get cell information
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
+ assert(!dispT.isNull());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ assert(dt > 0);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ if (cellDim != spaceDim)
+ throw std::logic_error("Not implemented yet.");
+
+ // Allocate vector for cell values (if necessary)
+ _initCellMatrix();
+
+ // Allocate vector for total strain
+ int tensorSize = 0;
+ if (1 == cellDim)
+ tensorSize = 1;
+ else if (2 == cellDim)
+ tensorSize = 3;
+ else if (3 == cellDim)
+ tensorSize = 6;
+ else
+ throw std::logic_error("Tensor size not implemented for given cellDim.");
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+ // Loop over cells
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+
+ // Set cell data in material
+ _material->initCellData(*c_iter, numQuadPts);
+
+ // Reset element vector to zero
+ _resetCellMatrix();
+
+ // Restrict input fields to cell
+ const real_section_type::value_type* dispTCell =
+ mesh->restrict(dispT, *c_iter);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ if (1 == cellDim) { // 1-D case
+ // Compute strains
+ Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+
+ // Get "elasticity" matrix at quadrature points for this cell
+ const std::vector<double_array>& elasticConsts =
+ _material->calcDerivElastic(totalStrain);
+
+ // Compute Jacobian for consistent tangent matrix
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double C1111 = elasticConsts[iQuad][0];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double valI = wt*basisDeriv[iQ+iBasis]*C1111;
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double valIJ = valI * basisDeriv[iQ+jBasis];
+ _cellMatrix[iBlock+jBlock] += valIJ;
+ } // for
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ } else if (2 == cellDim) { // 2-D case
+ // Compute strains
+ Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+
+ // Get "elasticity" matrix at quadrature points for this cell
+ const std::vector<double_array>& elasticConsts =
+ _material->calcDerivElastic(totalStrain);
+
+ // Compute Jacobian for consistent tangent matrix
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double C1111 = elasticConsts[iQuad][0];
+ const double C1122 = elasticConsts[iQuad][1];
+ const double C1112 = elasticConsts[iQuad][2];
+ const double C2222 = elasticConsts[iQuad][3];
+ const double C2212 = elasticConsts[iQuad][4];
+ const double C1212 = elasticConsts[iQuad][5];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double Nip = wt*basisDeriv[iQ+iBasis*cellDim ];
+ const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double Njp = basisDeriv[iQ+jBasis*cellDim ];
+ const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp;
+ _cellMatrix[iBlock +jBlock ] += ki0j0;
+ _cellMatrix[iBlock +jBlock+1] += ki0j1;
+ _cellMatrix[iBlock+1+jBlock ] += ki0j1;
+ _cellMatrix[iBlock+1+jBlock+1] += ki1j1;
+ } // for
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ } else if (3 == cellDim) { // 3-D case
+ // Compute strains
+ Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
+ // Get "elasticity" matrix at quadrature points for this cell
+ const std::vector<double_array>& elasticConsts =
+ _material->calcDerivElastic(totalStrain);
+
+ // Compute Jacobian for consistent tangent matrix
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double C1111 = elasticConsts[iQuad][ 0];
+ const double C1122 = elasticConsts[iQuad][ 1];
+ const double C1133 = elasticConsts[iQuad][ 2];
+ const double C1112 = elasticConsts[iQuad][ 3];
+ const double C1123 = elasticConsts[iQuad][ 4];
+ const double C1113 = elasticConsts[iQuad][ 5];
+ const double C2222 = elasticConsts[iQuad][ 6];
+ const double C2233 = elasticConsts[iQuad][ 7];
+ const double C2212 = elasticConsts[iQuad][ 8];
+ const double C2223 = elasticConsts[iQuad][ 9];
+ const double C2213 = elasticConsts[iQuad][10];
+ const double C3333 = elasticConsts[iQuad][11];
+ const double C3312 = elasticConsts[iQuad][12];
+ const double C3323 = elasticConsts[iQuad][13];
+ const double C3313 = elasticConsts[iQuad][14];
+ const double C1212 = elasticConsts[iQuad][15];
+ const double C1223 = elasticConsts[iQuad][16];
+ const double C1213 = elasticConsts[iQuad][17];
+ const double C2323 = elasticConsts[iQuad][18];
+ const double C2313 = elasticConsts[iQuad][19];
+ const double C1313 = elasticConsts[iQuad][20];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double Nip = wt*basisDeriv[iQ+iBasis*cellDim+0];
+ const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ const double Nir = wt*basisDeriv[iQ+iBasis*cellDim+2];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
+ const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
+ const double Njr = basisDeriv[iQ+jBasis*cellDim+2];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
+ C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
+ C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
+ const double ki0j2 =
+ C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
+ C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
+ C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
+ C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
+ const double ki1j2 =
+ C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
+ C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
+ C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
+ const double ki2j2 =
+ C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
+ C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
+ C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
+
+ _cellMatrix[iBlock +jBlock ] += ki0j0;
+ _cellMatrix[iBlock+1+jBlock ] += ki0j1;
+ _cellMatrix[iBlock+2+jBlock ] += ki0j2;
+ _cellMatrix[iBlock +jBlock+1] += ki0j1;
+ _cellMatrix[iBlock+1+jBlock+1] += ki1j1;
+ _cellMatrix[iBlock+2+jBlock+1] += ki1j2;
+ _cellMatrix[iBlock +jBlock+2] += ki0j2;
+ _cellMatrix[iBlock+1+jBlock+2] += ki1j2;
+ _cellMatrix[iBlock+2+jBlock+2] += ki2j2;
+ } // for
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'."
+ << std::endl;
+ assert(0);
+ } // if/else
+
+ // Assemble cell contribution into field. Not sure if this is correct for
+ // global stiffness matrix.
+ const ALE::Obj<Mesh::order_type>& globalOrder =
+ mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
+ err = updateOperator(*mat, mesh, dispT, globalOrder,
+ *c_iter, _cellMatrix, ADD_VALUES);
+ if (err)
+ throw std::runtime_error("Update to PETSc Mat failed.");
+ } // for
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Update state variables as needed.
+void
+pylith::feassemble::ElasticityImplicit::updateState(
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // updateState
+ assert(0 != _material);
+ assert(!disp.isNull());
+
+ // No need to update state if using elastic behavior
+ if (_material->useElasticBehavior())
+ return;
+
+ // Get cell information
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ // Allocate vector for total strain
+ int tensorSize = 0;
+ if (1 == cellDim)
+ tensorSize = 1;
+ else if (2 == cellDim)
+ tensorSize = 3;
+ else if (3 == cellDim)
+ tensorSize = 6;
+ else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'." << std::endl;
+ assert(0);
+ } // else
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+
+ // Set cell data in material
+ _material->initCellData(*c_iter, numQuadPts);
+
+ // Restrict input fields to cell
+ const real_section_type::value_type* dispCell =
+ mesh->restrict(disp, *c_iter);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ // Compute action for elastic terms
+ if (1 == cellDim) {
+ Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispCell, numBasis);
+ _material->updateState(totalStrain);
+ } else if (2 == cellDim) {
+ Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispCell, numBasis);
+ _material->updateState(totalStrain);
+ } else if (3 == cellDim) {
+ // Compute stresses
+ Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
+ dispCell, numBasis);
+ _material->updateState(totalStrain);
+ } else {
+ std::cerr << "Unknown case for cellDim '" << cellDim << "'."
+ << std::endl;
+ assert(0);
+ } // if/else
+ } // for
+} // updateState
+
+
+// End of file
Copied: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh (from rev 6964, short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.hh 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,155 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/ElasticityImplicit.hh
+ *
+ * @brief Implicit time integration of quasi-static elasticity equation
+ * using finite-elements.
+ *
+ * Note: This object operates on a single finite-element family, which
+ * is defined by the quadrature and a database of material property
+ * parameters.
+ *
+ * Computes contributions to terms A and r in
+ *
+ * A(t) u(t+dt) = b(u(t), u(t-dt)),
+ *
+ * r = b - A u0(t+dt)
+ *
+ * where A(t) is a sparse matrix or vector, u(t+dt) is the field we
+ * want to compute at time t+dt, b is a vector that depends on the
+ * field at time t and t-dt, and u0 is zero at unknown DOF and set to
+ * the constrained values at known DOF.
+ *
+ * Contributions to the RHS (b) include body forces, which are either
+ * independent of u (small strain case) or are computed based on the
+ * displacements at time t. The RHS also includes the internal force
+ * vector, which is either constant for a time step (small strain,
+ * elastic rheology) or changes with each iteration (large strain or
+ * non-elastic rheology). The internal force vector is subtracted from the
+ * existing force vector to get the residual. This object also computes
+ * the entire stiffness matrix (A).
+ *
+ * See governing equations section of user manual for more
+ * information.
+ */
+
+#if !defined(pylith_feassemble_elasticityimplicit_hh)
+#define pylith_feassemble_elasticityimplicit_hh
+
+#include "Integrator.hh" // ISA Integrator
+
+namespace pylith {
+ namespace feassemble {
+ class ElasticityImplicit;
+ class TestElasticityImplicit;
+ } // feassemble
+
+ namespace materials {
+ class ElasticMaterial;
+ } // feassemble
+} // pylith
+
+namespace spatialdata {
+ namespace spatialdb {
+ class SpatialDB; // USES SpatialDB
+ } // spatialdb
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+
+class pylith::feassemble::ElasticityImplicit : public Integrator
+{ // ElasticityImplicit
+ friend class TestElasticityImplicit; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElasticityImplicit(void);
+
+ /// Destructor
+ ~ElasticityImplicit(void);
+
+ /** Set material.
+ *
+ * @param m Elastic material.
+ */
+ void material(const materials::ElasticMaterial* m);
+
+ /** Set time step for advancing from time t to time t+dt.
+ *
+ * @param dt Time step
+ */
+ void timeStep(const double dt);
+
+ /** Integrate residual part of RHS for 3-D finite elements.
+ * Includes gravity and element internal force contribution.
+ *
+ * We assume that the effects of boundary conditions are already
+ * included in b (tractions, concentrated nodal forces, and
+ * contributions to internal force vector due to
+ * displacement/velocity BC). This routine computes the additional
+ * external loads due to body forces (not yet implemented) plus the
+ * element internal forces for the current stress state.
+ *
+ * @param b Residual field (output)
+ * @param fields Solution fields
+ * @param mesh Mesh object
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& b,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Compute Jacobian matrix (A) associated with operator.
+ *
+ * @param mat Sparse matrix
+ * @param fields Solution fields
+ * @param mesh Mesh object
+ */
+ void integrateJacobian(PetscMat* mat,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Update state variables as needed.
+ *
+ * @param field Current solution field.
+ * @param mesh Finite-element mesh
+ */
+ void updateState(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<Mesh>& mesh);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ ElasticityImplicit(const ElasticityImplicit& i);
+
+ /// Not implemented
+ const ElasticityImplicit& operator=(const ElasticityImplicit&);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double _dtm1; ///< Time step for t-dt1 -> t
+
+ /// Elastic material associated with integrator
+ materials::ElasticMaterial* _material;
+
+}; // ElasticityImplicit
+
+#endif // pylith_feassemble_elasticityimplicit_hh
+
+
+// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,351 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "ExplicitElasticity.hh" // implementation of class methods
-
-#include "Quadrature.hh" // USES Quadrature
-#include "IntegratorElasticity.hh" // USES IntegratorElasticity
-#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
-#include "pylith/utils/array.hh" // USES double_array
-
-#include "petscmat.h" // USES PetscMat
-#include "spatialdata/spatialdb/SpatialDB.hh"
-
-#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-
-// ----------------------------------------------------------------------
-// Constructor
-pylith::feassemble::ExplicitElasticity::ExplicitElasticity(void) :
- _material(0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor
-pylith::feassemble::ExplicitElasticity::~ExplicitElasticity(void)
-{ // destructor
- delete _material; _material = 0;
-} // destructor
-
-// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::feassemble::ExplicitElasticity::ExplicitElasticity(const ExplicitElasticity& i) :
- IntegratorExplicit(i),
- _material(0)
-{ // copy constructor
- if (0 != i._material)
- _material = i._material->clone();
-} // copy constructor
-
-// ----------------------------------------------------------------------
-// Set material.
-void
-pylith::feassemble::ExplicitElasticity::material(
- const materials::ElasticMaterial* m)
-{ // material
- delete _material; _material = (0 != m) ? m->clone() : 0;
-} // material
-
-// ----------------------------------------------------------------------
-// Integrate constant term (b) for dynamic elasticity term for 3-D
-// finite elements.
-void
-pylith::feassemble::ExplicitElasticity::integrateConstant(
- const ALE::Obj<real_section_type>& b,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& dispTmdt,
- const ALE::Obj<Mesh>& mesh)
-{ // integrateConstant
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(!b.isNull());
- assert(!dispT.isNull());
- assert(!dispTmdt.isNull());
- assert(!mesh.isNull());
-
- // Get information about section
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
-
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
-
- /** :TODO:
- *
- * If cellDim and spaceDim are different, we need to transform
- * displacements into cellDim, compute action, and transform result
- * back into spaceDim. We get this information from the Jacobian and
- * inverse of the Jacobian.
- */
- if (cellDim != spaceDim)
- throw std::logic_error("Not implemented yet.");
-
- // Allocate vector for cell values (if necessary)
- _initCellVector();
-
- // Allocate vector for total strain
- int tensorSize = 0;
- if (1 == cellDim)
- tensorSize = 1;
- else if (2 == cellDim)
- tensorSize = 3;
- else if (3 == cellDim)
- tensorSize = 6;
- else
- throw std::logic_error("Tensor size not implemented for given cellDim.");
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
-
- for (Mesh::label_sequence::iterator cellIter=cells->begin();
- cellIter != cellsEnd;
- ++cellIter) {
- // Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *cellIter);
-
- // Set cell data in material
- _material->initCellData(*cellIter, numQuadPts);
-
- // Reset element vector to zero
- _resetCellVector();
-
- // Restrict input fields to cell
- const real_section_type::value_type* dispTCell =
- mesh->restrict(dispT, *cellIter);
- const real_section_type::value_type* dispTmdtCell =
- mesh->restrict(dispTmdt, *cellIter);
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& basisDeriv = _quadrature->basisDeriv();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- // Compute action for cell
-
- // Compute action for inertial terms
- const std::vector<double_array>& density =
- _material->calcDensity();
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double valI = wt*basis[iQuad*numBasis+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const double valIJ = valI * basis[iQuad*numBasis+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis*spaceDim+iDim] +=
- valIJ * (2.0 * dispTCell[jBasis*spaceDim+iDim] -
- dispTmdtCell[jBasis*spaceDim+iDim]);
- } // for
- } // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- // Compute action for elastic terms
- if (1 == cellDim) {
- // Compute stresses
- IntegratorElasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis*cellDim ];
- _cellVector[iBlock ] -= N1*s11;
- } // for
- } // for
- PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*5));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- } else if (2 == cellDim) {
- // Compute stresses
- IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s12 = stress[iQuad][2];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
- const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
- _cellVector[iBlock ] -= N1*s11 + N2*s12;
- _cellVector[iBlock+1] -= N1*s12 + N2*s22;
- } // for
- } // for
- err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- } else if (3 == cellDim) {
- // Compute stresses
- IntegratorElasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s33 = stress[iQuad][2];
- const double s12 = stress[iQuad][3];
- const double s23 = stress[iQuad][4];
- const double s13 = stress[iQuad][5];
-
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double N1 = wt*basisDeriv[iQ+iBasis*cellDim+0];
- const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
- const double N3 = wt*basisDeriv[iQ+iBasis*cellDim+2];
- _cellVector[iBlock ] -= N1*s11 + N2*s12 + N3*s13;
- _cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
- _cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
- } // for
- } // for
- err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } // if/else
-
- // Assemble cell contribution into field
- mesh->updateAdd(b, *cellIter, _cellVector);
- } // for
-} // integrateConstant
-
-// ----------------------------------------------------------------------
-// Compute matrix associated with operator.
-void
-pylith::feassemble::ExplicitElasticity::integrateJacobian(
- PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh)
-{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != mat);
- assert(!dispT.isNull());
- assert(!mesh.isNull());
-
- // Get information about section
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
-
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
-
- // Allocate vector for cell values (if necessary)
- _initCellMatrix();
-
- for (Mesh::label_sequence::iterator cellIter=cells->begin();
- cellIter != cellsEnd;
- ++cellIter) {
- // Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *cellIter);
-
- // Set cell data in material
- _material->initCellData(*cellIter, numQuadPts);
-
- // Reset element vector to zero
- _resetCellMatrix();
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- // Get material physical properties at quadrature points for this cell
- const std::vector<double_array>& density = _material->calcDensity();
-
- // Compute Jacobian for cell
-
- // Compute Jacobian for inertial terms
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const double valI = wt*basis[iQ+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const double valIJ = valI * basis[iQ+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const int iBlock = (iBasis * spaceDim + iDim) * (spaceDim * numBasis);
- const int jBlock = (jBasis * spaceDim + iDim);
- _cellMatrix[iBlock+jBlock] += valIJ;
- } // for
- } // for
- } // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- // Assemble cell contribution into PETSc Matrix
- const ALE::Obj<Mesh::order_type>& globalOrder =
- mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
-
- err = updateOperator(*mat, mesh, dispT, globalOrder,
- *cellIter, _cellMatrix, ADD_VALUES);
- if (err)
- throw std::runtime_error("Update to PETSc Mat failed.");
- } // for
-} // integrateJacobian
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,152 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/ExplicitElasticity.hh
- *
- * @brief Explicit time integration of dynamic elasticity equation
- * using finite-elements.
- *
- * Note: This object operates on a single finite-element family, which
- * is defined by the quadrature and a database of material property
- * parameters.
- *
- * Computes contributions to terms A and b in A(t) u(t+dt) = b(u(t),
- * u(t-dt)), where A(t) is a sparse matrix or vector, u(t+dt) is the
- * field we want to compute at time t+dt and b is a vector that
- * depends on the field at time t and t-dt.
- *
- * Contributions from elasticity include the intertial and stiffness
- * terms, so this object computes the following portions of A and b:
- *
- * A = 1/(dt*dt) [M]
- *
- * b = 2/(dt*dt)[M]{u(t)} - 1/(dt*dt)[M]{u(t-dt)} - [K]{u(t)}
- *
- * Translational inertia.
- * - Residual action over cell
- * \f[
- * \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
- * \f]
- * - Jacobian action over cell
- * \f[
- * \int_{V^e} (\rho N^q N^q)_i \, dV
- * \f]
- * - Integrate and lump to form lumped matrix (field)
- *
- * See governing equations section of user manual for more
- * information.
- */
-
-#if !defined(pylith_feassemble_explicitelasticity_hh)
-#define pylith_feassemble_explicitelasticity_hh
-
-#include "IntegratorExplicit.hh" // ISA IntegratorExplicit
-
-namespace pylith {
- namespace feassemble {
- class ExplicitElasticity;
- class TestExplicitElasticity;
- } // feassemble
-
- namespace materials {
- class ElasticMaterial;
- } // feassemble
-} // pylith
-
-namespace spatialdata {
- namespace spatialdb {
- class SpatialDB; // USES SpatialDB
- } // spatialdb
- namespace geocoords {
- class CoordSys; // USES CoordSys
- } // geocoords
-} // spatialdata
-
-class pylith::feassemble::ExplicitElasticity : public IntegratorExplicit
-{ // ExplicitElasticity
- friend class TestExplicitElasticity; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /// Constructor
- ExplicitElasticity(void);
-
- /// Destructor
- ~ExplicitElasticity(void);
-
- /// Create a copy of this object.
- IntegratorExplicit* clone(void) const;
-
- /** Set material.
- *
- * @param m Elastic material.
- */
- void material(const materials::ElasticMaterial* m);
-
- /** Integrate constant term (b) for dynamic elasticity term
- * for finite elements.
- *
- * Compute b = 2/(dt*dt)[M]{u(t) - 1/(dt*dt)[M]{u(t-dt)} - [K]{u(t)}, where
- * [M] = density * [N]^T [N]
- *
- *
- * @param b Constant field (output)
- * @param dispT Displacement field at time t
- * @param dispTmdt Displacement field at time t-dt
- * @param mesh Finite-element mesh
- */
- void integrateConstant(const ALE::Obj<real_section_type>& b,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& dispTmdt,
- const ALE::Obj<Mesh>& mesh);
-
- /** Compute Jacobian matrix (A) associated with operator.
- *
- * @param mat Sparse matrix
- * @param dispT Displacement at time t
- * @param mesh Finite-element mesh
- */
- void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh);
-
-// PROTECTED METHODS ////////////////////////////////////////////////////
-protected :
-
- /** Copy constructor.
- *
- * @param i Integrator to copy
- */
- ExplicitElasticity(const ExplicitElasticity& i);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
- /// Not implemented
- const ExplicitElasticity& operator=(const ExplicitElasticity&);
-
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- /// Elastic material associated with integrator
- materials::ElasticMaterial* _material;
-
-}; // ExplicitElasticity
-
-#include "ExplicitElasticity.icc" // inline methods
-
-#endif // pylith_feassemble_explicitelasticity_hh
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.icc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.icc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,27 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#if !defined(pylith_feassemble_explicitelasticity_hh)
-#error "ExplicitElasticity.icc must be included only from ExplicitElasticity.hh"
-#else
-
-// Create a copy of this object.
-inline
-pylith::feassemble::IntegratorExplicit*
-pylith::feassemble::ExplicitElasticity::clone(void) const {
- return new ExplicitElasticity(*this);
-} // clone
-
-#endif
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,514 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "ImplicitElasticity.hh" // implementation of class methods
-
-#include "Quadrature.hh" // USES Quadrature
-#include "IntegratorElasticity.hh" // USES IntegratorElasticity
-#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
-#include "pylith/utils/array.hh" // USES double_array
-
-#include "petscmat.h" // USES PetscMat
-#include "spatialdata/spatialdb/SpatialDB.hh"
-
-#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-
-// ----------------------------------------------------------------------
-// Constructor
-pylith::feassemble::ImplicitElasticity::ImplicitElasticity(void) :
- _material(0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor
-pylith::feassemble::ImplicitElasticity::~ImplicitElasticity(void)
-{ // destructor
- delete _material; _material = 0;
-} // destructor
-
-// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::feassemble::ImplicitElasticity::ImplicitElasticity(const ImplicitElasticity& i) :
- IntegratorImplicit(i),
- _material(0)
-{ // copy constructor
- if (0 != i._material)
- _material = i._material->clone();
-} // copy constructor
-
-// ----------------------------------------------------------------------
-// Set material.
-void
-pylith::feassemble::ImplicitElasticity::material(
- const materials::ElasticMaterial* m)
-{ // material
- delete _material; _material = (0 != m) ? m->clone() : 0;
-} // material
-
-// ----------------------------------------------------------------------
-// Compute residual for quasi-static finite elements.
-// We assume that the effects of boundary conditions are already included
-// in b (tractions, concentrated nodal forces, and contributions to internal
-// force vector due to displacement/velocity BC).
-// This routine computes the additional external loads due to body forces
-// (not yet implemented) plus the element internal forces for the current
-// stress state.
-void
-pylith::feassemble::ImplicitElasticity::integrateResidual(
- const ALE::Obj<real_section_type>& b,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh)
-{ // integrateResidual
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(!b.isNull());
- assert(!dispT.isNull());
- assert(!mesh.isNull());
-
- // Get information about section
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
-
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
-
- // Allocate vector for cell values (if necessary)
- _initCellVector();
-
- // Allocate vector for total strain
- int tensorSize = 0;
- if (1 == cellDim)
- tensorSize = 1;
- else if (2 == cellDim)
- tensorSize = 3;
- else if (3 == cellDim)
- tensorSize = 6;
- else
- throw std::logic_error("Tensor size not implemented for given cellDim.");
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
-
- // Loop over cells
- for (Mesh::label_sequence::iterator cellIter=cells->begin();
- cellIter != cellsEnd;
- ++cellIter) {
- // Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *cellIter);
-
- // Set cell data in material
- _material->initCellData(*cellIter, numQuadPts);
-
- // Reset element vector to zero
- _resetCellVector();
-
- // Restrict input fields to cell
- const real_section_type::value_type* dispTCell =
- mesh->restrict(dispT, *cellIter);
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& basisDeriv = _quadrature->basisDeriv();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- if (cellDim != spaceDim)
- throw std::logic_error("Not implemented yet.");
-
- /* Comment out gravity section for now, until we figure out how to deal
- with gravity vector.
- // Get density at quadrature points for this cell
- const std::vector<double_array>& density = _material->calcDensity();
-
- // Compute action for cell
-
- // Compute action for element body forces
- if (!grav.isNull()) {
- const real_section_type::value_type* gravCell =
- mesh->restrict(grav, cell);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double valI = wt*basis[iQ+iBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- _cellVector[iBlock+iDim] += valI*gravCell[iDim];
- } // for
- } // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(2+numBasis*(2+2*spaceDim)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } // if
- */
-
- // Compute B(transpose) * sigma, first computing strains
- if (1 == cellDim) {
- // Compute total strains and then use these to compute stresses
- IntegratorElasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- // Need a way to tell calcStress whether to update state variables or not.
- // Otherwise, need a separate update routine.
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
- _cellVector[iBlock ] -= N1*s11;
- } // for
- } // for
- PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*5));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- } else if (2 == cellDim) {
- // Compute total strains and then use these to compute stresses
- IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s12 = stress[iQuad][2];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
- const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
- _cellVector[iBlock ] -= N1*s11 + N2*s12;
- _cellVector[iBlock+1] -= N1*s12 + N2*s22;
- } // for
- } // for
- PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- } else if (3 == cellDim) {
- // Compute total strains and then use these to compute stresses
- IntegratorElasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s33 = stress[iQuad][2];
- const double s12 = stress[iQuad][3];
- const double s23 = stress[iQuad][4];
- const double s13 = stress[iQuad][5];
-
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double N1 = wt*basisDeriv[iQ+iBasis*cellDim+0];
- const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
- const double N3 = wt*basisDeriv[iQ+iBasis*cellDim+2];
- _cellVector[iBlock ] -= N1*s11 + N2*s12 + N3*s13;
- _cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
- _cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
- } // for
- } // for
- PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } // if/else
-
- // Assemble cell contribution into field
- mesh->updateAdd(b, *cellIter, _cellVector);
- } // for
-} // integrateResidual
-
-
-// ----------------------------------------------------------------------
-// Compute stiffness matrix.
-void
-pylith::feassemble::ImplicitElasticity::integrateJacobian(
- PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh)
-{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != mat);
- assert(!dispT.isNull());
- assert(!mesh.isNull());
-
- // Clear stiffness matrix. Not sure if this is the correct way.
- PetscErrorCode err = MatZeroEntries(*mat);
-
- // Get information about section
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
-
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Get parameters used in integration.
- const double dt = _dt;
- assert(dt > 0);
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
-
- if (cellDim != spaceDim)
- throw std::logic_error("Not implemented yet.");
-
- // Allocate vector for cell values (if necessary)
- _initCellMatrix();
-
- // Allocate vector for total strain
- int tensorSize = 0;
- if (1 == cellDim)
- tensorSize = 1;
- else if (2 == cellDim)
- tensorSize = 3;
- else if (3 == cellDim)
- tensorSize = 6;
- else
- throw std::logic_error("Tensor size not implemented for given cellDim.");
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
-
- // Loop over cells
- for (Mesh::label_sequence::iterator cellIter=cells->begin();
- cellIter != cellsEnd;
- ++cellIter) {
- // Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *cellIter);
-
- // Set cell data in material
- _material->initCellData(*cellIter, numQuadPts);
-
- // Reset element vector to zero
- _resetCellMatrix();
-
- // Restrict input fields to cell
- const real_section_type::value_type* dispTCell =
- mesh->restrict(dispT, *cellIter);
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& basisDeriv = _quadrature->basisDeriv();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- // 1D Case
- if (1 == cellDim) {
- // Compute strains
- IntegratorElasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
-
- // Get "elasticity" matrix at quadrature points for this cell
- const std::vector<double_array>& elasticConsts =
- _material->calcDerivElastic(totalStrain);
-
- // Compute Jacobian for consistent tangent matrix
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double C1111 = elasticConsts[iQuad][0];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double valI = wt*basisDeriv[iQ+iBasis]*C1111;
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis * spaceDim;
- const double valIJ = valI * basisDeriv[iQ+jBasis];
- _cellMatrix[iBlock+jBlock] += valIJ;
- } // for
- } // for
- } // for
- err =
- PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- // 2D Case
- } else if (2 == cellDim) {
- // Compute strains
- IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
-
- // Get "elasticity" matrix at quadrature points for this cell
- const std::vector<double_array>& elasticConsts =
- _material->calcDerivElastic(totalStrain);
-
- // Compute Jacobian for consistent tangent matrix
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double C1111 = elasticConsts[iQuad][0];
- const double C1122 = elasticConsts[iQuad][1];
- const double C1112 = elasticConsts[iQuad][2];
- const double C2222 = elasticConsts[iQuad][3];
- const double C2212 = elasticConsts[iQuad][4];
- const double C1212 = elasticConsts[iQuad][5];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double Nip = wt*basisDeriv[iQ+iBasis*cellDim ];
- const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis * spaceDim;
- const double Njp = basisDeriv[iQ+jBasis*cellDim ];
- const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
- const double ki0j0 =
- C1111 * Nip * Njp + C1112 * Niq * Njp +
- C1112 * Nip * Njq + C1212 * Niq * Njq;
- const double ki0j1 =
- C1122 * Nip * Njq + C2212 * Niq * Njq +
- C1112 * Nip * Njp + C1212 * Niq * Njp;
- const double ki1j1 =
- C2222 * Niq * Njq + C2212 * Nip * Njq +
- C2212 * Niq * Njp + C1212 * Nip * Njp;
- _cellMatrix[iBlock +jBlock ] += ki0j0;
- _cellMatrix[iBlock +jBlock+1] += ki0j1;
- _cellMatrix[iBlock+1+jBlock ] += ki0j1;
- _cellMatrix[iBlock+1+jBlock+1] += ki1j1;
- } // for
- } // for
- } // for
- err =
- PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
- // 3D Case
- } else if (3 == cellDim) {
- // Compute strains
- IntegratorElasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
-
- // Get "elasticity" matrix at quadrature points for this cell
- const std::vector<double_array>& elasticConsts =
- _material->calcDerivElastic(totalStrain);
-
- // Compute Jacobian for consistent tangent matrix
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double C1111 = elasticConsts[iQuad][ 0];
- const double C1122 = elasticConsts[iQuad][ 1];
- const double C1133 = elasticConsts[iQuad][ 2];
- const double C1112 = elasticConsts[iQuad][ 3];
- const double C1123 = elasticConsts[iQuad][ 4];
- const double C1113 = elasticConsts[iQuad][ 5];
- const double C2222 = elasticConsts[iQuad][ 6];
- const double C2233 = elasticConsts[iQuad][ 7];
- const double C2212 = elasticConsts[iQuad][ 8];
- const double C2223 = elasticConsts[iQuad][ 9];
- const double C2213 = elasticConsts[iQuad][10];
- const double C3333 = elasticConsts[iQuad][11];
- const double C3312 = elasticConsts[iQuad][12];
- const double C3323 = elasticConsts[iQuad][13];
- const double C3313 = elasticConsts[iQuad][14];
- const double C1212 = elasticConsts[iQuad][15];
- const double C1223 = elasticConsts[iQuad][16];
- const double C1213 = elasticConsts[iQuad][17];
- const double C2323 = elasticConsts[iQuad][18];
- const double C2313 = elasticConsts[iQuad][19];
- const double C1313 = elasticConsts[iQuad][20];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * spaceDim;
- const double Nip = wt*basisDeriv[iQ+iBasis*cellDim+0];
- const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
- const double Nir = wt*basisDeriv[iQ+iBasis*cellDim+2];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis * spaceDim;
- const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
- const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
- const double Njr = basisDeriv[iQ+jBasis*cellDim+2];
- const double ki0j0 =
- C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
- C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
- C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
- const double ki0j1 =
- C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
- C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
- C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
- const double ki0j2 =
- C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
- C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
- C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
- const double ki1j1 =
- C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
- C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
- C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
- const double ki1j2 =
- C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
- C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
- C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
- const double ki2j2 =
- C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
- C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
- C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
-
- _cellMatrix[iBlock +jBlock ] += ki0j0;
- _cellMatrix[iBlock+1+jBlock ] += ki0j1;
- _cellMatrix[iBlock+2+jBlock ] += ki0j2;
- _cellMatrix[iBlock +jBlock+1] += ki0j1;
- _cellMatrix[iBlock+1+jBlock+1] += ki1j1;
- _cellMatrix[iBlock+2+jBlock+1] += ki1j2;
- _cellMatrix[iBlock +jBlock+2] += ki0j2;
- _cellMatrix[iBlock+1+jBlock+2] += ki1j2;
- _cellMatrix[iBlock+2+jBlock+2] += ki2j2;
- } // for
- } // for
- } // for
- err =
- PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } // if/else
-
- // Assemble cell contribution into field. Not sure if this is correct for
- // global stiffness matrix.
- const ALE::Obj<Mesh::order_type>& globalOrder =
- mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
- err = updateOperator(*mat, mesh, dispT, globalOrder,
- *cellIter, _cellMatrix, ADD_VALUES);
- if (err)
- throw std::runtime_error("Update to PETSc Mat failed.");
- } // for
-} // integrateJacobian
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,138 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/ImplicitElasticity.hh
- *
- * @brief Implicit time integration of quasi-static elasticity equation
- * using finite-elements.
- *
- * Note: This object operates on a single finite-element family, which
- * is defined by the quadrature and a database of material property
- * parameters.
- *
- * Computes contributions to terms A and b in A(t+dt) u(t+dt) = b(u(t),
- * u(t+dt)), where A(t+dt) is a sparse matrix, u(t+dt) is the
- * field we want to compute at time t+dt and b is a vector that
- * depends on the field at time t and t+dt.
- *
- * Contributions to the RHS (b) include body forces, which are either
- * independent of u (small strain case) or are computed based on the
- * displacements at time t. The RHS also includes the internal force
- * vector, which is either constant for a time step (small strain,
- * elastic rheology) or changes with each iteration (large strain or
- * non-elastic rheology). The internal force vector is subtracted from the
- * existing force vector to get the residual. This object also computes
- * the entire stiffness matrix (A).
- *
- * See governing equations section of user manual for more
- * information.
- */
-
-#if !defined(pylith_feassemble_implicitelasticity_hh)
-#define pylith_feassemble_implicitelasticity_hh
-
-#include "IntegratorImplicit.hh" // ISA IntegratorImplicit
-
-namespace pylith {
- namespace feassemble {
- class ImplicitElasticity;
- class TestImplicitElasticity;
- } // feassemble
-
- namespace materials {
- class ElasticMaterial;
- } // feassemble
-} // pylith
-
-namespace spatialdata {
- namespace spatialdb {
- class SpatialDB; // USES SpatialDB
- } // spatialdb
- namespace geocoords {
- class CoordSys; // USES CoordSys
- } // geocoords
-} // spatialdata
-
-class pylith::feassemble::ImplicitElasticity : public IntegratorImplicit
-{ // ImplicitElasticity
- friend class TestImplicitElasticity; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /// Constructor
- ImplicitElasticity(void);
-
- /// Destructor
- ~ImplicitElasticity(void);
-
- /// Create a copy of this object.
- IntegratorImplicit* clone(void) const;
-
- /** Set material.
- *
- * @param m Elastic material.
- */
- void material(const materials::ElasticMaterial* m);
-
- /** Integrate residual part of RHS for 3-D finite elements.
- * Includes gravity and element internal force contribution.
- *
- *
- * @param b Residual field (output)
- * @param dispT Displacement field at time t
- * @param mesh Mesh object
- */
- void integrateResidual(const ALE::Obj<real_section_type>& b,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh);
-
- /** Compute Jacobian matrix (A) associated with operator.
- *
- * @param mat Sparse matrix
- * @param dispT Displacement at time t
- * @param mesh Mesh object
- */
- void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<Mesh>& mesh);
-
-// PROTECTED METHODS ////////////////////////////////////////////////////
-protected :
-
- /** Copy constructor.
- *
- * @param i Integrator to copy
- */
- ImplicitElasticity(const ImplicitElasticity& i);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
- /// Not implemented
- const ImplicitElasticity& operator=(const ImplicitElasticity&);
-
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- /// Elastic material associated with integrator
- materials::ElasticMaterial* _material;
-
-}; // ImplicitElasticity
-
-#include "ImplicitElasticity.icc" // inline methods
-
-#endif // pylith_feassemble_implicitelasticity_hh
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.icc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.icc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,27 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#if !defined(pylith_feassemble_implicitelasticity_hh)
-#error "ImplicitElasticity.icc must be included only from ImplicitElasticity.hh"
-#else
-
-// Create a copy of this object.
-inline
-pylith::feassemble::IntegratorImplicit*
-pylith::feassemble::ImplicitElasticity::clone(void) const {
- return new ImplicitElasticity(*this);
-} // clone
-
-#endif
-
-
-// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -21,9 +21,11 @@
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::Integrator::Integrator(void) :
+ _dt(-1.0),
_quadrature(0),
_cellVector(0),
- _cellMatrix(0)
+ _cellMatrix(0),
+ _needNewJacobian(false)
{ // constructor
} // constructor
@@ -37,17 +39,6 @@
} // destructor
// ----------------------------------------------------------------------
-// Copy constructor
-pylith::feassemble::Integrator::Integrator(const Integrator& i) :
- _quadrature(0),
- _cellVector(0),
- _cellMatrix(0)
-{ // copy constructor
- if (0 != i._quadrature)
- _quadrature = i._quadrature->clone();
-} // copy constructor
-
-// ----------------------------------------------------------------------
// Set quadrature for integrating finite-element quantities.
void
pylith::feassemble::Integrator::quadrature(const Quadrature* q)
@@ -61,6 +52,40 @@
} // quadrature
// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::Integrator::timeStep(const double dt)
+{ // timeStep
+ _dt = dt;
+} // timeStep
+
+// ----------------------------------------------------------------------
+// Get stable time step for advancing from time t to time t+dt.
+double
+pylith::feassemble::Integrator::stableTimeStep(void) const
+{ // stableTimeStep
+ // Default is current time step
+ return _dt;
+} // stableTimeStep
+
+// ----------------------------------------------------------------------
+// Check whether Jacobian needs to be recomputed.
+bool
+pylith::feassemble::Integrator::needNewJacobian(void) const
+{ // needNewJacobian
+ return _needNewJacobian;
+} // needNewJacobian
+
+// ----------------------------------------------------------------------
+// Update state variables as needed.
+void
+pylith::feassemble::Integrator::updateState(
+ const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<Mesh>& mesh)
+{ // updateState
+} // updateState
+
+// ----------------------------------------------------------------------
// Initialize vector containing result of integration action for cell.
void
pylith::feassemble::Integrator::_initCellVector(void)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -24,7 +24,8 @@
#if !defined(pylith_feassemble_integrator_hh)
#define pylith_feassemble_integrator_hh
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "pylith/utils/sievetypes.hh" // USES real_section_type
+#include "pylith/utils/petscfwd.h" // USES PetscMat
namespace pylith {
namespace feassemble {
@@ -33,6 +34,10 @@
class Quadrature; // HOLDSA Quadrature
} // feassemble
+
+ namespace topology {
+ class FieldsManager;
+ } // topology
} // pylith
class pylith::feassemble::Integrator
@@ -56,15 +61,63 @@
*/
void quadrature(const Quadrature* q);
-// PROTECTED METHODS ////////////////////////////////////////////////////
-protected :
+ /** Set time step for advancing from time t to time t+dt.
+ *
+ * @param dt Time step
+ */
+ virtual
+ void timeStep(const double dt);
- /** Copy constructor.
+ /** Get stable time step for advancing from time t to time t+dt.
*
- * @param i Integrator to copy
+ * Default is current time step.
+ *
+ * @returns Time step
*/
- Integrator(const Integrator& i);
+ virtual
+ double stableTimeStep(void) const;
+ /** Check whether Jacobian needs to be recomputed.
+ *
+ * @returns True if Jacobian needs to be recomputed, false otherwise.
+ */
+ bool needNewJacobian(void) const;
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ virtual
+ void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh) = 0;
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param mat Sparse matrix
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ virtual
+ void integrateJacobian(PetscMat* mat,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh) = 0;
+
+ /** Update state variables as needed.
+ *
+ * @param field Current solution field.
+ * @param mesh Finite-element mesh
+ */
+ virtual
+ void updateState(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<Mesh>& mesh);
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
/// Initialize vector containing result of integration action for cell.
void _initCellVector(void);
@@ -80,12 +133,17 @@
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
+ // Not implemented.
+ Integrator(const Integrator& i);
+
/// Not implemented
const Integrator& operator=(const Integrator&);
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
+ double _dt; ///< Time step for t -> t+dt
+
Quadrature* _quadrature; ///< Quadrature for integrating finite-element
/// Vector local to cell containing result of integration action
@@ -94,6 +152,10 @@
/// Matrix local to cell containing result of integration
real_section_type::value_type* _cellMatrix;
+ /// True if we need to recompute Jacobian for operator, false otherwise.
+ /// Default is false;
+ bool _needNewJacobian;
+
}; // Integrator
#endif // pylith_feassemble_integrator_hh
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,113 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "IntegratorElasticity.hh" // implementation of class methods
-
-#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::logic_error
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::IntegratorElasticity::calcTotalStrain1D(
- std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double* disp,
- const int numBasis)
-{ // calcTotalStrain1D
- assert(0 != strain);
- assert(0 != disp);
-
- const int dimension = 1;
- const int numQuadPts = strain->size();
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(1 == (*strain)[iQuad].size());
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- (*strain)[iQuad][0] +=
- basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
- } // for
-} // calcTotalStrain1D
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::IntegratorElasticity::calcTotalStrain2D(
- std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double* disp,
- const int numBasis)
-{ // calcTotalStrain2D
- assert(0 != strain);
- assert(0 != disp);
-
- const int dimension = 2;
- const int numQuadPts = strain->size();
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(3 == (*strain)[iQuad].size());
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad][0] +=
- basisDeriv[iQ+iBasis ] * disp[iBasis ];
- (*strain)[iQuad][1] +=
- basisDeriv[iQ+iBasis+1] * disp[iBasis+1];
- (*strain)[iQuad][2] +=
- 0.5 * (basisDeriv[iQ+iBasis+1] * disp[iBasis ] +
- basisDeriv[iQ+iBasis ] * disp[iBasis+1]);
- } // for
- } // for
-} // calcTotalStrain2D
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::IntegratorElasticity::calcTotalStrain3D(
- std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double* disp,
- const int numBasis)
-{ // calcTotalStrain3D
- assert(0 != strain);
- assert(0 != disp);
-
- const int dimension = 3;
- const int numQuadPts = strain->size();
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(6 == (*strain)[iQuad].size());
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad][0] +=
- basisDeriv[iQ+iBasis ] * disp[iBasis ];
- (*strain)[iQuad][1] +=
- basisDeriv[iQ+iBasis+1] * disp[iBasis+1];
- (*strain)[iQuad][2] +=
- basisDeriv[iQ+iBasis+2] * disp[iBasis+2];
- (*strain)[iQuad][3] +=
- 0.5 * (basisDeriv[iQ+iBasis+1] * disp[iBasis ] +
- basisDeriv[iQ+iBasis ] * disp[iBasis+1]);
- (*strain)[iQuad][4] +=
- 0.5 * (basisDeriv[iQ+iBasis+2] * disp[iBasis+1] +
- basisDeriv[iQ+iBasis+1] * disp[iBasis+2]);
- (*strain)[iQuad][5] +=
- 0.5 * (basisDeriv[iQ+iBasis+2] * disp[iBasis ] +
- basisDeriv[iQ+iBasis ] * disp[iBasis+2]);
- } // for
- } // for
-} // calcTotalStrain3D
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,86 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/IntegratorElasticity.hh
- *
- * @brief C++ Utility class for general operations associated with
- * integrating finite-element actions for solid elasticity.
- */
-
-#if !defined(pylith_feassemble_integratorelasticity_hh)
-#define pylith_feassemble_integratorelasticity_hh
-
-#include "pylith/utils/array.hh" // USES double_array, std::vector
-
-namespace pylith {
- namespace feassemble {
- class IntegratorElasticity;
- class TestIntegratorElasticity;
- } // feassemble
-} // pylith
-
-class pylith::feassemble::IntegratorElasticity
-{ // IntegratorElasticity
- friend class TestIntegratorElasticity; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /** Compute total strain in at quadrature points of a cell.
- *
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param dimension Dimension of cell.
- * @param numBasis Number of basis functions for cell.
- */
-
- static
- void calcTotalStrain1D(std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double* disp,
- const int numBasis);
-
- /** Compute total strain in at quadrature points of a cell.
- *
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
- */
-
- static
- void calcTotalStrain2D(std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double* disp,
- const int numBasis);
-
- /** Compute total strain in at quadrature points of a cell.
- *
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
- */
-
- static
- void calcTotalStrain3D(std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double* disp,
- const int numBasis);
-
-}; // IntegratorElasticity
-
-#endif // pylith_feassemble_integratorelasticity_hh
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,66 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "IntegratorExplicit.hh" // implementation of class methods
-
-#include <assert.h> // USES assert()
-
-// ----------------------------------------------------------------------
-// Constructor
-pylith::feassemble::IntegratorExplicit::IntegratorExplicit(void) :
- Integrator(),
- _dt(-1.0),
- _dtm1(-1.0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor
-pylith::feassemble::IntegratorExplicit::~IntegratorExplicit(void)
-{ // destructor
-} // destructor
-
-// ----------------------------------------------------------------------
-// Copy constructor
-pylith::feassemble::IntegratorExplicit::IntegratorExplicit(const IntegratorExplicit& i) :
- Integrator(i),
- _dt(i._dt),
- _dtm1(i._dtm1)
-{ // copy constructor
-} // copy constructor
-
-// ----------------------------------------------------------------------
-// Set time step for advancing from time t to time t+dt.
-void
-pylith::feassemble::IntegratorExplicit::timeStep(const double dt)
-{ // timeStep
- if (_dt != -1.0)
- _dtm1 = _dt;
- else
- _dtm1 = dt;
- _dt = dt;
- assert(_dt == _dtm1); // For now, don't allow variable time step
-} // timeStep
-
-// ----------------------------------------------------------------------
-// Get stable time step for advancing from time t to time t+dt.
-double
-pylith::feassemble::IntegratorExplicit::stableTimeStep(void) const
-{ // stableTimeStep
- // Default is current time step
- return _dt;
-} // stableTimeStep
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,127 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/IntegratorExplicit.hh
- *
- * @brief Abstract base class for explicit time integration of
- * finite-element actions.
- *
- * Note: This object operates on a single finite-element family, which
- * is defined by the quadrature and a database of material property
- * parameters.
- *
- * Computes terms A and b in A(t) u(t+dt) = b(u(t), u(t-dt)), where
- * A(t) is a sparse matrix or vector, u(t+dt) is the field we want to
- * compute at time t+dt and b is a vector that depends on the field at
- * time t and t-dt.
- */
-
-#if !defined(pylith_feassemble_integratorexplicit_hh)
-#define pylith_feassemble_integratorexplicit_hh
-
-#include "pylith/utils/petscfwd.h" // USES PetscMat
-
-#include "Integrator.hh" // ISA Integrator
-
-namespace pylith {
- namespace feassemble {
- class IntegratorExplicit;
- class TestIntegratorExplicit;
- } // feassemble
-} // pylith
-
-class pylith::feassemble::IntegratorExplicit : public Integrator
-{ // Integrator
- friend class TestIntegratorExplicit; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /// Constructor
- IntegratorExplicit(void);
-
- /// Destructor
- virtual
- ~IntegratorExplicit(void);
-
- /// Create a copy of this object.
- virtual
- IntegratorExplicit* clone(void) const = 0;
-
- /** Set time step for advancing from time t to time t+dt.
- *
- * @param dt Time step
- */
- void timeStep(const double dt);
-
- /** Get stable time step for advancing from time t to time t+dt.
- *
- * Default is current time step.
- *
- * @returns Time step
- */
- virtual
- double stableTimeStep(void) const;
-
- /** Integrate constant term (b) for dynamic elasticity term
- * for finite elements.
- *
- * @param fieldOut Constant field (output)
- * @param fieldInT Input field at time t
- * @param fieldInTmdt Input field at time t-dt
- * @param mesh Finite-element mesh
- */
- virtual
- void integrateConstant(const ALE::Obj<real_section_type>& fieldOut,
- const ALE::Obj<real_section_type>& fieldInT,
- const ALE::Obj<real_section_type>& fieldInTmdt,
- const ALE::Obj<Mesh>& mesh) = 0;
-
- /** Compute Jacobian matrix (A) associated with operator.
- *
- * @param mat Sparse matrix
- * @param fieldIn Input field at time t
- * @param mesh Finite-element mesh
- */
- virtual
- void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& fieldIn,
- const ALE::Obj<Mesh>& mesh) = 0;
-
-// PROTECTED METHODS ////////////////////////////////////////////////////
-protected :
-
- /** Copy constructor.
- *
- * @param i Integrator to copy
- */
- IntegratorExplicit(const IntegratorExplicit& i);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
- /// Not implemented
- const IntegratorExplicit& operator=(const IntegratorExplicit&);
-
-// PROTECTED MEMBERS ////////////////////////////////////////////////////
-protected :
-
- double _dt; ///< Time step for t -> t+dt
- double _dtm1; ///< Time step for t-dt1 -> t
-
-}; // IntegratorExplicit
-
-#endif // pylith_feassemble_integratorexplicit_hh
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,66 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "IntegratorImplicit.hh" // implementation of class methods
-
-#include <assert.h> // USES assert()
-
-// ----------------------------------------------------------------------
-// Constructor
-pylith::feassemble::IntegratorImplicit::IntegratorImplicit(void) :
- Integrator(),
- _dt(-1.0),
- _dtm1(-1.0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor
-pylith::feassemble::IntegratorImplicit::~IntegratorImplicit(void)
-{ // destructor
-} // destructor
-
-// ----------------------------------------------------------------------
-// Copy constructor
-pylith::feassemble::IntegratorImplicit::IntegratorImplicit(const IntegratorImplicit& i) :
- Integrator(i),
- _dt(i._dt),
- _dtm1(i._dtm1)
-{ // copy constructor
-} // copy constructor
-
-// ----------------------------------------------------------------------
-// Set time step for advancing from time t to time t+dt.
-void
-pylith::feassemble::IntegratorImplicit::timeStep(const double dt)
-{ // timeStep
- if (_dt != -1.0)
- _dtm1 = _dt;
- else
- _dtm1 = dt;
- _dt = dt;
- assert(_dt == _dtm1); // For now, don't allow variable time step
-} // timeStep
-
-// ----------------------------------------------------------------------
-// Get stable time step for advancing from time t to time t+dt.
-double
-pylith::feassemble::IntegratorImplicit::stableTimeStep(void) const
-{ // stableTimeStep
- // Default is current time step
- return _dt;
-} // stableTimeStep
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,124 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/IntegratorImplicit.hh
- *
- * @brief Abstract base class for implicit time integration of
- * finite-element actions.
- *
- * Note: This object operates on a single finite-element family, which
- * is defined by the quadrature and a database of material property
- * parameters.
- *
- * Computes terms A and b in A(t) u(t+dt) = b(u(t), u(t-dt)), where
- * A(t) is a sparse matrix or vector, u(t+dt) is the field we want to
- * compute at time t+dt and b is a vector that depends on the field at
- * time t and t-dt.
- */
-
-#if !defined(pylith_feassemble_integratorimplicit_hh)
-#define pylith_feassemble_integratorimplicit_hh
-
-#include "pylith/utils/petscfwd.h" // USES PetscMat
-
-#include "Integrator.hh" // ISA Integrator
-
-namespace pylith {
- namespace feassemble {
- class IntegratorImplicit;
- class TestIntegratorImplicit;
- } // feassemble
-} // pylith
-
-class pylith::feassemble::IntegratorImplicit : public Integrator
-{ // Integrator
- friend class TestIntegratorImplicit; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /// Constructor
- IntegratorImplicit(void);
-
- /// Destructor
- virtual
- ~IntegratorImplicit(void);
-
- /// Create a copy of this object.
- virtual
- IntegratorImplicit* clone(void) const = 0;
-
- /** Set time step for advancing from time t to time t+dt.
- *
- * @param dt Time step
- */
- void timeStep(const double dt);
-
- /** Get stable time step for advancing from time t to time t+dt.
- *
- * Default is current time step.
- *
- * @returns Time step
- */
- virtual
- double stableTimeStep(void) const;
-
- /** Integrate residual for quasi-static finite elements.
- *
- * @param fieldOut Residual field (output)
- * @param fieldInT Input field at time t
- * @param mesh Mesh object
- */
- virtual
- void integrateResidual(const ALE::Obj<real_section_type>& fieldOut,
- const ALE::Obj<real_section_type>& fieldInT,
- const ALE::Obj<Mesh>& mesh) = 0;
-
- /** Compute Jacobian matrix (A) associated with operator.
- *
- * @param mat Sparse matrix
- * @param fieldIn Input field at time t
- * @param coordinates Field of cell vertex coordinates
- */
- virtual
- void integrateJacobian(PetscMat* mat,
- const ALE::Obj<real_section_type>& fieldIn,
- const ALE::Obj<Mesh>& mesh) = 0;
-
-// PROTECTED METHODS ////////////////////////////////////////////////////
-protected :
-
- /** Copy constructor.
- *
- * @param i Integrator to copy
- */
- IntegratorImplicit(const IntegratorImplicit& i);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
- /// Not implemented
- const IntegratorImplicit& operator=(const IntegratorImplicit&);
-
-// PROTECTED MEMBERS ////////////////////////////////////////////////////
-protected :
-
- double _dt; ///< Time step for t -> t+dt
- double _dtm1; ///< Time step for t-dt1 -> t
-
-}; // IntegratorImplicit
-
-#endif // pylith_feassemble_integratorimplicit_hh
-
-
-// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-05-25 20:57:26 UTC (rev 6970)
@@ -15,14 +15,10 @@
subpkginclude_HEADERS = \
Constraint.hh \
- ExplicitElasticity.hh \
- ExplicitElasticity.icc \
- ImplicitElasticity.hh \
- ImplicitElasticity.icc \
+ ElasticityExplicit.hh \
+ ElasticityImplicit.hh \
Integrator.hh \
- IntegratorElasticity.hh \
- IntegratorExplicit.hh \
- IntegratorImplicit.hh \
+ Elasticity.hh \
Quadrature.hh \
Quadrature.icc \
Quadrature1D.hh \
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -24,7 +24,8 @@
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(void) :
- _numQuadPts(0)
+ _numQuadPts(0),
+ _useElasticBehavior(true)
{ // constructor
} // constructor
@@ -38,7 +39,8 @@
// Copy constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(const ElasticMaterial& m) :
Material(m),
- _numQuadPts(0)
+ _numQuadPts(0),
+ _useElasticBehavior(m._useElasticBehavior)
{ // copy constructor
} // copy constructor
Modified: short/3D/PyLith/trunk/libsrc/meshio/BinaryIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/BinaryIO.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/meshio/BinaryIO.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -25,7 +25,8 @@
{ // readString
char* buffer = (numChars > 0) ? new char[numChars+1] : 0;
fin.read(buffer, sizeof(char)*numChars);
- buffer[numChars] = '\0';
+ if (numChars > 0)
+ buffer[numChars] = '\0';
// get string from buffer
std::string bufstring(buffer);
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -87,8 +87,10 @@
assert(coordinates.size() == numVertices*spaceDim);
assert(cells.size() == numCells*numCorners);
+ // :BUG: This causes a memory leak.
*_mesh = new Mesh(PETSC_COMM_WORLD, meshDim);
_mesh->addRef();
+
assert(!_mesh->isNull());
(*_mesh)->setDebug(_debug);
Modified: short/3D/PyLith/trunk/libsrc/topology/FieldsManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/FieldsManager.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/libsrc/topology/FieldsManager.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -46,9 +46,7 @@
throw std::runtime_error(msg.str());
} // if
- ALE::Obj<real_section_type> field =
- new real_section_type(_mesh->comm(), _mesh->debug());
- _real[name] = field;
+ _real[name] = new real_section_type(_mesh->comm(), _mesh->debug());
} // addReal
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2007-05-25 20:57:26 UTC (rev 6970)
@@ -20,10 +20,8 @@
#include "pylith/feassemble/Quadrature3D.hh"
#include "pylith/feassemble/Integrator.hh"
-#include "pylith/feassemble/IntegratorExplicit.hh"
-#include "pylith/feassemble/IntegratorImplicit.hh"
-#include "pylith/feassemble/ExplicitElasticity.hh"
-#include "pylith/feassemble/ImplicitElasticity.hh"
+#include "pylith/feassemble/ElasticityExplicit.hh"
+#include "pylith/feassemble/ElasticityImplicit.hh"
#include "pylith/feassemble/Constraint.hh"
@@ -430,76 +428,25 @@
return
- def _createHandle(self):
+ def integrateResidual(self, residual, fields, mesh):
"""
- Wrap pointer to C++ object in PyCObject.
+ Integrate contributions to residual term (r) for operator.
"""
- return PyCObject_FromVoidPtr(self.thisptr, Integrator_destructor)
-
-
- property quadrature:
- def __set__(self, q):
- """
- Set quadrature.
- """
- # create shim for method 'quadrature'
- #embed{ void Integrator_quadrature_set(void* objVptr, void* qVptr)
- try {
- assert(0 != objVptr);
- pylith::feassemble::Quadrature* quadrature =
- (pylith::feassemble::Quadrature*) qVptr;
- ((pylith::feassemble::Integrator*) objVptr)->quadrature(quadrature);
- } catch (const std::exception& err) {
- PyErr_SetString(PyExc_RuntimeError,
- const_cast<char*>(err.what()));
- } catch (...) {
- PyErr_SetString(PyExc_RuntimeError,
- "Caught unknown C++ exception.");
- } // try/catch
- #}embed
- if not q.name == "pylith_feassemble_Quadrature":
- raise TypeError, \
- "Argument must be extension module type 'Quadrature'."
- Integrator_quadrature_set(self.thisptr, ptrFromHandle(q))
-
-
-# ----------------------------------------------------------------------
-cdef class IntegratorExplicit(Integrator):
-
- def __init__(self):
- """
- Constructor.
- """
- Integrator.__init__(self)
- return
-
-
- def integrateConstant(self, fieldOut, fieldInT, fieldInTmdt, mesh):
- """
- Integrate constant term (b) for dynamic elasticity term for 3-D
- finite elements.
- """
- # create shim for method 'integrateConstant'
- #embed{ void IntegratorExplicit_integrateConstant(void* objVptr, void* fieldOutVptr, void* fieldInTVptr, void* fieldInTmdtVptr, void* meshVptr)
- typedef ALE::Mesh Mesh;
- typedef ALE::Mesh::real_section_type real_section_type;
-
+ # create shim for method 'integrateResidual'
+ #embed{ void Integrator_integrateResidual(void* objVptr, void* residualVptr, void* fieldsVptr, void* meshVptr)
try {
assert(0 != objVptr);
- assert(0 != fieldOutVptr);
- assert(0 != fieldInTVptr);
- assert(0 != fieldInTmdtVptr);
+ assert(0 != residualVptr);
+ assert(0 != fieldsVptr);
assert(0 != meshVptr);
- ALE::Obj<Mesh>* mesh =
- (ALE::Obj<Mesh>*) meshVptr;
- ALE::Obj<real_section_type>* fieldOut =
- (ALE::Obj<real_section_type>*) fieldOutVptr;
- ALE::Obj<real_section_type>* fieldInT =
- (ALE::Obj<real_section_type>*) fieldInTVptr;
- ALE::Obj<real_section_type>* fieldInTmdt =
- (ALE::Obj<real_section_type>*) fieldInTmdtVptr;
- ((pylith::feassemble::IntegratorExplicit*) objVptr)->integrateConstant(
- *fieldOut, *fieldInT, *fieldInTmdt, *mesh);
+ ALE::Obj<ALE::Mesh>* mesh =
+ (ALE::Obj<ALE::Mesh>*) meshVptr;
+ ALE::Obj<pylith::real_section_type>* residual =
+ (ALE::Obj<pylith::real_section_type>*) residualVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ ((pylith::feassemble::Integrator*) objVptr)->integrateResidual(*residual,
+ fields, *mesh);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -511,40 +458,34 @@
"Caught unknown C++ exception.");
} // try/catch
#}embed
- cdef void* fieldOutVptr
- cdef void* fieldInTVptr
- cdef void* fieldInTmdtVptr
- fieldOutVptr = PyCObject_AsVoidPtr(fieldOut)
- fieldInTVptr = PyCObject_AsVoidPtr(fieldInT)
- fieldInTmdtVptr = PyCObject_AsVoidPtr(fieldInTmdt)
- IntegratorExplicit_integrateConstant(self.thisptr,
- fieldOutVptr, fieldInTVptr,
- fieldInTmdtVptr,
- ptrFromHandle(mesh))
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument must be extension module type 'Mesh'."
+ Integrator_integrateResidual(self.thisptr,
+ PyCObject_AsVoidPtr(residual),
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh))
return
- def integrateJacobian(self, mat, fieldIn, mesh):
+ def integrateJacobian(self, mat, fields, mesh):
"""
- Compute matrix (A) associated with operator.
+ Compute contributions to Jacobian matrix (A) associated with operator.
"""
# create shim for method 'integrateJacobian'
- #embed{ void IntegratorExplicit_integrateJacobian(void* objVptr, void* matVptr, void* fieldInVptr, void* meshVptr)
- typedef ALE::Mesh Mesh;
- typedef ALE::Mesh::real_section_type real_section_type;
-
+ #embed{ void Integrator_integrateJacobian(void* objVptr, void* matVptr, void* fieldsVptr, void* meshVptr)
try {
assert(0 != objVptr);
assert(0 != matVptr);
- assert(0 != fieldInVptr);
+ assert(0 != fieldsVptr);
assert(0 != meshVptr);
- ALE::Obj<Mesh>* mesh =
- (ALE::Obj<Mesh>*) meshVptr;
+ ALE::Obj<ALE::Mesh>* mesh =
+ (ALE::Obj<ALE::Mesh>*) meshVptr;
PetscMat* mat = (PetscMat*) matVptr;
- ALE::Obj<real_section_type>* fieldIn =
- (ALE::Obj<real_section_type>*) fieldInVptr;
- ((pylith::feassemble::IntegratorExplicit*) objVptr)->integrateJacobian(
- mat, *fieldIn, *mesh);
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ ((pylith::feassemble::Integrator*) objVptr)->integrateJacobian(
+ mat, fields, *mesh);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -556,70 +497,31 @@
"Caught unknown C++ exception.");
} // try/catch
#}embed
- cdef void* matVptr
- cdef void* fieldInVptr
- matVptr = PyCObject_AsVoidPtr(mat)
- fieldInVptr = PyCObject_AsVoidPtr(fieldIn)
- IntegratorExplicit_integrateJacobian(self.thisptr, matVptr, fieldInVptr,
- ptrFromHandle(mesh))
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type 'Mesh'."
+ Integrator_integrateJacobian(self.thisptr,
+ PyCObject_AsVoidPtr(mat),
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh))
return
- property timeStep:
- def __set__(self, dt):
- """
- Set timeStep.
- """
- # create shim for method 'timeStep'
- #embed{ void IntegratorExplicit_timeStep_set(void* objVptr, double dt)
- try {
- assert(0 != objVptr);
- ((pylith::feassemble::IntegratorExplicit*) objVptr)->timeStep(dt);
- } catch (const std::exception& err) {
- PyErr_SetString(PyExc_RuntimeError,
- const_cast<char*>(err.what()));
- } catch (...) {
- PyErr_SetString(PyExc_RuntimeError,
- "Caught unknown C++ exception.");
- } // try/catch
- #}embed
- IntegratorExplicit_timeStep_set(self.thisptr, dt)
-
-
-# ----------------------------------------------------------------------
-cdef class IntegratorImplicit(Integrator):
-
- def __init__(self):
+ def updateState(self, field, mesh):
"""
- Constructor.
+ Update state variables as needed.
"""
- Integrator.__init__(self)
- return
-
-
- def integrateResidual(self, fieldOut, fieldInT, mesh):
- """
- Integrate constant term (b) for quasi-static elasticity term for
- finite elements.
- """
- # create shim for method 'integrateResidual'
- #embed{ void IntegratorImplicit_integrateResidual(void* objVptr, void* fieldOutVptr, void* fieldInTVptr, void* meshVptr)
- typedef ALE::Mesh Mesh;
- typedef ALE::Mesh::real_section_type real_section_type;
-
+ # create shim for method 'updateState'
+ #embed{ void Integrator_updateState(void* objVptr, void* fieldVptr, void* meshVptr)
try {
assert(0 != objVptr);
- assert(0 != fieldOutVptr);
- assert(0 != fieldInTVptr);
+ assert(0 != fieldVptr);
assert(0 != meshVptr);
- ALE::Obj<Mesh>* mesh =
- (ALE::Obj<Mesh>*) meshVptr;
- ALE::Obj<real_section_type>* fieldOut =
- (ALE::Obj<real_section_type>*) fieldOutVptr;
- ALE::Obj<real_section_type>* fieldInT =
- (ALE::Obj<real_section_type>*) fieldInTVptr;
- ((pylith::feassemble::IntegratorImplicit*) objVptr)->integrateResidual(
- *fieldOut, *fieldInT, *mesh);
+ ALE::Obj<pylith::real_section_type>* field =
+ (ALE::Obj<pylith::real_section_type>*) fieldVptr;
+ ALE::Obj<ALE::Mesh>* mesh =
+ (ALE::Obj<ALE::Mesh>*) meshVptr;
+ ((pylith::feassemble::Integrator*) objVptr)->updateState(*field, *mesh);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -631,67 +533,57 @@
"Caught unknown C++ exception.");
} // try/catch
#}embed
- cdef void* fieldOutVptr
- cdef void* fieldInTVptr
- fieldOutVptr = PyCObject_AsVoidPtr(fieldOut)
- fieldInTVptr = PyCObject_AsVoidPtr(fieldInT)
- IntegratorImplicit_integrateResidual(self.thisptr,
- fieldOutVptr, fieldInTVptr,
- ptrFromHandle(mesh))
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type 'Mesh'."
+ Integrator_updateState(self.thisptr, PyCObject_AsVoidPtr(field),
+ ptrFromHandle(mesh))
return
- def integrateJacobian(self, mat, fieldIn, mesh):
+ def _createHandle(self):
"""
- Compute matrix (A) associated with operator.
+ Wrap pointer to C++ object in PyCObject.
"""
- # create shim for method 'integrateJacobian'
- #embed{ void IntegratorImplicit_integrateJacobian(void* objVptr, void* matVptr, void* fieldInVptr, void* meshVptr)
- typedef ALE::Mesh Mesh;
- typedef ALE::Mesh::real_section_type real_section_type;
+ return PyCObject_FromVoidPtr(self.thisptr, Integrator_destructor)
- try {
- assert(0 != objVptr);
- assert(0 != matVptr);
- assert(0 != fieldInVptr);
- assert(0 != meshVptr);
- ALE::Obj<Mesh>* mesh =
- (ALE::Obj<Mesh>*) meshVptr;
- PetscMat* mat = (PetscMat*) matVptr;
- ALE::Obj<real_section_type>* fieldIn =
- (ALE::Obj<real_section_type>*) fieldInVptr;
- ((pylith::feassemble::IntegratorImplicit*) objVptr)->integrateJacobian(
- mat, *fieldIn, *mesh);
- } catch (const std::exception& err) {
- PyErr_SetString(PyExc_RuntimeError,
- const_cast<char*>(err.what()));
- } catch (const ALE::Exception& err) {
- PyErr_SetString(PyExc_RuntimeError,
- const_cast<char*>(err.msg().c_str()));
- } catch (...) {
- PyErr_SetString(PyExc_RuntimeError,
- "Caught unknown C++ exception.");
- } // try/catch
- #}embed
- cdef void* matVptr
- cdef void* fieldInVptr
- matVptr = PyCObject_AsVoidPtr(mat)
- fieldInVptr = PyCObject_AsVoidPtr(fieldIn)
- IntegratorImplicit_integrateJacobian(self.thisptr, matVptr, fieldInVptr,
- ptrFromHandle(mesh))
- return
+ property quadrature:
+ def __set__(self, q):
+ """
+ Set quadrature.
+ """
+ # create shim for method 'quadrature'
+ #embed{ void Integrator_quadrature_set(void* objVptr, void* qVptr)
+ try {
+ assert(0 != objVptr);
+ pylith::feassemble::Quadrature* quadrature =
+ (pylith::feassemble::Quadrature*) qVptr;
+ ((pylith::feassemble::Integrator*) objVptr)->quadrature(quadrature);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if not q.name == "pylith_feassemble_Quadrature":
+ raise TypeError, \
+ "Argument must be extension module type 'Quadrature'."
+ Integrator_quadrature_set(self.thisptr, ptrFromHandle(q))
+
property timeStep:
def __set__(self, dt):
"""
Set timeStep.
"""
# create shim for method 'timeStep'
- #embed{ void IntegratorImplicit_timeStep_set(void* objVptr, double dt)
+ #embed{ void Integrator_timeStep_set(void* objVptr, double dt)
try {
assert(0 != objVptr);
- ((pylith::feassemble::IntegratorImplicit*) objVptr)->timeStep(dt);
+ ((pylith::feassemble::Integrator*) objVptr)->timeStep(dt);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -700,21 +592,44 @@
"Caught unknown C++ exception.");
} // try/catch
#}embed
- IntegratorImplicit_timeStep_set(self.thisptr, dt)
+ Integrator_timeStep_set(self.thisptr, dt)
+ property needNewJacobian:
+ def __get__(self):
+ """
+ Set timeStep.
+ """
+ # create shim for method 'needNewJacobian'
+ #embed{ int Integrator_needNewJacobian_get(void* objVptr)
+ int result = 0;
+ try {
+ assert(0 != objVptr);
+ ((pylith::feassemble::Integrator*) objVptr)->needNewJacobian();
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ return result;
+ #}embed
+ return Integrator_needNewJacobian_get(self.thisptr)
+
+
# ----------------------------------------------------------------------
-cdef class ExplicitElasticity(IntegratorExplicit):
+cdef class ElasticityExplicit(Integrator):
def __init__(self):
"""
Constructor.
"""
# create shim for constructor
- #embed{ void* ExplicitElasticity_constructor()
+ #embed{ void* ElasticityExplicit_constructor()
void* result = 0;
try {
- result = (void*)(new pylith::feassemble::ExplicitElasticity);
+ result = (void*)(new pylith::feassemble::ElasticityExplicit);
assert(0 != result);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
@@ -726,8 +641,8 @@
return result;
#}embed
- IntegratorExplicit.__init__(self)
- self.thisptr = ExplicitElasticity_constructor()
+ Integrator.__init__(self)
+ self.thisptr = ElasticityExplicit_constructor()
self.handle = self._createHandle()
return
@@ -738,12 +653,12 @@
Set material.
"""
# create shim for method 'material'
- #embed{ void ExplicitElasticity_material_set(void* objVptr, void* mVptr)
+ #embed{ void ElasticityExplicit_material_set(void* objVptr, void* mVptr)
try {
assert(0 != objVptr);
pylith::materials::ElasticMaterial* material =
(pylith::materials::ElasticMaterial*) mVptr;
- ((pylith::feassemble::ExplicitElasticity*) objVptr)->material(material);
+ ((pylith::feassemble::ElasticityExplicit*) objVptr)->material(material);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -755,21 +670,21 @@
if not m.name == "pylith_materials_ElasticMaterial":
raise TypeError, \
"Argument must be extension module type 'ElasticMaterial'."
- ExplicitElasticity_material_set(self.thisptr, ptrFromHandle(m))
+ ElasticityExplicit_material_set(self.thisptr, ptrFromHandle(m))
# ----------------------------------------------------------------------
-cdef class ImplicitElasticity(IntegratorImplicit):
+cdef class ElasticityImplicit(Integrator):
def __init__(self):
"""
Constructor.
"""
# create shim for constructor
- #embed{ void* ImplicitElasticity_constructor()
+ #embed{ void* ElasticityImplicit_constructor()
void* result = 0;
try {
- result = (void*)(new pylith::feassemble::ImplicitElasticity);
+ result = (void*)(new pylith::feassemble::ElasticityImplicit);
assert(0 != result);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
@@ -781,8 +696,8 @@
return result;
#}embed
- IntegratorImplicit.__init__(self)
- self.thisptr = ImplicitElasticity_constructor()
+ Integrator.__init__(self)
+ self.thisptr = ElasticityImplicit_constructor()
self.handle = self._createHandle()
return
@@ -793,12 +708,12 @@
Set material.
"""
# create shim for method 'material'
- #embed{ void ImplicitElasticity_material_set(void* objVptr, void* mVptr)
+ #embed{ void ElasticityImplicit_material_set(void* objVptr, void* mVptr)
try {
assert(0 != objVptr);
pylith::materials::ElasticMaterial* material =
(pylith::materials::ElasticMaterial*) mVptr;
- ((pylith::feassemble::ImplicitElasticity*) objVptr)->material(material);
+ ((pylith::feassemble::ElasticityImplicit*) objVptr)->material(material);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -810,7 +725,7 @@
if not m.name == "pylith_materials_ElasticMaterial":
raise TypeError, \
"Argument must be extension module type 'ElasticMaterial'."
- ImplicitElasticity_material_set(self.thisptr, ptrFromHandle(m))
+ ElasticityImplicit_material_set(self.thisptr, ptrFromHandle(m))
# ----------------------------------------------------------------------
@@ -854,7 +769,7 @@
if not mesh.name == "pylith_topology_Mesh":
raise TypeError, \
- "Argument must be extension module type " \
+ "Argument 'mesh' must be extension module type " \
"'pylith::topology::Mesh'."
Constraint_setConstraintSizes(self.thisptr, PyCObject_AsVoidPtr(field),
ptrFromHandle(mesh))
@@ -885,7 +800,7 @@
if not mesh.name == "pylith_topology_Mesh":
raise TypeError, \
- "Argument must be extension module type " \
+ "Argument 'mesh' must be extension module type " \
"'pylith::topology::Mesh'."
Constraint_setConstraints(self.thisptr, PyCObject_AsVoidPtr(field), ptrFromHandle(mesh))
return
@@ -915,7 +830,7 @@
if not mesh.name == "pylith_topology_Mesh":
raise TypeError, \
- "Argument must be extension module type " \
+ "Argument 'mesh' must be extension module type " \
"'pylith::topology::Mesh'."
Constraint_setField(self.thisptr, t, PyCObject_AsVoidPtr(field),
ptrFromHandle(mesh))
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2007-05-25 20:57:26 UTC (rev 6970)
@@ -28,13 +28,12 @@
faults/SlipTimeFn.py \
feassemble/__init__.py \
feassemble/Constraint.py \
- feassemble/ExplicitElasticity.py \
- feassemble/ImplicitElasticity.py \
+ feassemble/ElasticityExplicit.py \
+ feassemble/ElasticityImplicit.py \
feassemble/FIATCell.py \
feassemble/FIATSimplex.py \
feassemble/Integrator.py \
- feassemble/IntegratorExplicit.py \
- feassemble/IntegratorImplicit.py \
+ feassemble/IntegratorElasticity.py \
feassemble/ReferenceCell.py \
feassemble/quadrature/Quadrature.py \
feassemble/quadrature/Quadrature1D.py \
@@ -64,7 +63,6 @@
problems/BCPrism.py \
problems/BCQuadrilateral.py \
problems/BoundaryConditions.py \
- problems/EqDeformation.py \
problems/Explicit.py \
problems/Formulation.py \
problems/Implicit.py \
Modified: short/3D/PyLith/trunk/pylith/PyLithApp.py
===================================================================
--- short/3D/PyLith/trunk/pylith/PyLithApp.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/PyLithApp.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -42,19 +42,14 @@
import pyre.inventory
- from pyre.units.time import second
- totalTime = pyre.inventory.dimensional("total_time", default=0.0*second,
- validator=pyre.inventory.greaterEqual(0.0*second))
- totalTime.meta['tip'] = "Time duration for simulation."
-
from pylith.topology.MeshImporter import MeshImporter
mesher = pyre.inventory.facility("mesh_generator", family="mesh_generator",
factory=MeshImporter)
mesher.meta['tip'] = "Generates or imports the computational mesh."
- from pylith.problems.EqDeformation import EqDeformation
+ from pylith.problems.TimeDependent import TimeDependent
problem = pyre.inventory.facility("problem", family="problem",
- factory=EqDeformation)
+ factory=TimeDependent)
problem.meta['tip'] = "Computational problem to solve."
# Dummy facility for passing options to PETSc
Modified: short/3D/PyLith/trunk/pylith/faults/FaultsBin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultsBin.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/faults/FaultsBin.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -14,7 +14,7 @@
##
## @brief Python container for faults.
##
-## Factory: faults_bin
+## Factory: interfaces_bin
from pyre.components.Component import Component
@@ -23,7 +23,7 @@
"""
Python container for faults.
- Factory: faults_bin
+ Factory: interfaces_bin
"""
# PUBLIC METHODS /////////////////////////////////////////////////////
@@ -33,7 +33,7 @@
Constructor.
"""
Component.__init__(self, name, facility="faults_bin")
- self.faults = []
+ self.ic = []
return
Modified: short/3D/PyLith/trunk/pylith/faults/SingleFault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/SingleFault.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/faults/SingleFault.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -14,7 +14,7 @@
##
## @brief Python faults container with one fault.
##
-## Factory: faults_bin
+## Factory: interfaces_bin
from FaultsBin import FaultsBin
@@ -23,7 +23,7 @@
"""
Python faults container with one material.
- Factory: faults_bin
+ Factory: interfaces_bin
"""
# INVENTORY //////////////////////////////////////////////////////////
@@ -66,13 +66,13 @@
Set attributes from inventory.
"""
FaultsBin._configure(self)
- self.faults = [self.inventory.fault]
+ self.ic = [self.inventory.fault]
return
# FACTORIES ////////////////////////////////////////////////////////////
-def faults_bin():
+def interfaces_bin():
"""
Factory associated with SingleFault.
"""
Copied: short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicit.py (from rev 6964, short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py)
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/ElasticityExplicit.py
+##
+## @brief Python object for explicit time integration of dynamic
+## elasticity equation using finite-elements.
+##
+## Factory: integrator
+
+from IntegratorElasticity import IntegratorElasticity
+
+# ElasticityExplicit class
+class ElasticityExplicit(IntegratorElasticity):
+ """
+ Python object for explicit time integration of dynamic elasticity
+ equation using finite-elements.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="elasticityexplicit"):
+ """
+ Constructor.
+ """
+ IntegratorElasticity.__init__(self, name)
+
+ import pylith.feassemble.feassemble as bindings
+ self.cppHandle = bindings.ElasticityExplicit()
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def integrator():
+ """
+ Factory associated with ElasticityExplicit.
+ """
+ return ElasticityExplicit()
+
+
+# End of file
Copied: short/3D/PyLith/trunk/pylith/feassemble/ElasticityImplicit.py (from rev 6964, short/3D/PyLith/trunk/pylith/feassemble/ImplicitElasticity.py)
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ImplicitElasticity.py 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/pylith/feassemble/ElasticityImplicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/ElasticityImplicit.py
+##
+## @brief Python object for implicit time integration of dynamic
+## elasticity equation using finite-elements.
+##
+## Factory: integrator
+
+from IntegratorElasticity import IntegratorElasticity
+
+# ElasticityImplicit class
+class ElasticityImplicit(IntegratorElasticity):
+ """
+ Python object for implicit time integration of dynamic elasticity
+ equation using finite-elements.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="elasticityimplicit"):
+ """
+ Constructor.
+ """
+ IntegratorElasticity.__init__(self, name)
+
+ import pylith.feassemble.feassemble as bindings
+ self.cppHandle = bindings.ElasticityImplicit()
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def integrator():
+ """
+ Factory associated with ElasticityImplicit.
+ """
+ return ElasticityImplicit()
+
+
+# End of file
Deleted: short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,62 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file pylith/feassemble/ExplicitElasticity.py
-##
-## @brief Python object for explicit time integration of dynamic
-## elasticity equation using finite-elements.
-##
-## Factory: integrator
-
-from IntegratorExplicit import IntegratorExplicit
-
-# ExplicitElasticity class
-class ExplicitElasticity(IntegratorExplicit):
- """
- Python object for explicit time integration of dynamic elasticity
- equation using finite-elements.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="explicitelasticity"):
- """
- Constructor.
- """
- IntegratorExplicit.__init__(self, name)
-
- import pylith.feassemble.feassemble as bindings
- self.cppHandle = bindings.ExplicitElasticity()
- return
-
-
- def initMaterial(self, mesh, material):
- """
- Initialize material properties.
- """
- self._info.log("Initializing integrator material '%s'." % material.label)
- material.initialize(mesh)
- self.material = material
- self.cppHandle.material = self.material.cppHandle
- return
-
-
-# FACTORIES ////////////////////////////////////////////////////////////
-
-def integrator():
- """
- Factory associated with ExplicitElasticity.
- """
- return ExplicitElasticity()
-
-
-# End of file
Deleted: short/3D/PyLith/trunk/pylith/feassemble/ImplicitElasticity.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ImplicitElasticity.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/ImplicitElasticity.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,62 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file pylith/feassemble/ImplicitElasticity.py
-##
-## @brief Python object for implicit time integration of dynamic
-## elasticity equation using finite-elements.
-##
-## Factory: integrator
-
-from IntegratorImplicit import IntegratorImplicit
-
-# ImplicitElasticity class
-class ImplicitElasticity(IntegratorImplicit):
- """
- Python object for implicit time integration of dynamic elasticity
- equation using finite-elements.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="implicitelasticity"):
- """
- Constructor.
- """
- IntegratorImplicit.__init__(self, name)
-
- import pylith.feassemble.feassemble as bindings
- self.cppHandle = bindings.ImplicitElasticity()
- return
-
-
- def initMaterial(self, mesh, material):
- """
- Initialize material properties.
- """
- self._info.log("Initializing integrator material '%s'." % material.label)
- material.initialize(mesh)
- self.material = material
- self.cppHandle.material = self.material.cppHandle
- return
-
-
-# FACTORIES ////////////////////////////////////////////////////////////
-
-def integrator():
- """
- Factory associated with ImplicitElasticity.
- """
- return ImplicitElasticity()
-
-
-# End of file
Modified: short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -17,10 +17,23 @@
##
## Factory: fe_integrator.
-from pyre.components.Component import Component
+def implementsIntegrator(obj):
+ """
+ Check whether object implements an integrator.
+ """
+ result = True
+ attrs = dir(obj)
+ if not "initQuadrature" in attrs or \
+ not "timeStep" in attrs or \
+ not "stableTimeStep" in attrs or \
+ not "integrateResidual" in attrs or \
+ not "integrateJacobian" in attrs:
+ result = False
+ return result
+
# Integrator class
-class Integrator(Component):
+class Integrator(object):
"""
Python abstract base class for integration of actions with
finite-elements.
@@ -30,12 +43,10 @@
# PUBLIC METHODS /////////////////////////////////////////////////////
- def __init__(self, name="integrator"):
+ def __init__(self):
"""
Constructor.
"""
- Component.__init__(self, name, facility="integrator")
- self.cppHandle = None
self.quadrature = None
self.mesh = None
return
@@ -53,14 +64,64 @@
"""
Initialize quadrature.
"""
- if self.cppHandle is None:
- raise ValueError("C++ handle not set.")
-
- self._info.log("Initialize quadrature")
+ assert(None != self.cppHandle)
quadrature.initialize()
self.quadrature = quadrature
self.cppHandle.quadrature = self.quadrature.cppHandle
return
+ def timeStep(self, dt):
+ """
+ Set time step for advancing from time t to time t+dt.
+ """
+ assert(None != self.cppHandle)
+ self.cppHandle.timeStep = dt.value
+ return
+
+
+ def stableTimeStep(self):
+ """
+ Get stable time step for advancing from time t to time t+dt.
+ """
+ assert(None != self.cppHandle)
+ return self.cppHandle.getStableTimeStep()
+
+
+ def integrateResidual(self, residual, fields):
+ """
+ Integrate contributions to residual term at time t.
+ """
+ assert(None != self.cppHandle)
+ self.cppHandle.integrateResidual(residual, fields, self.mesh.cppHandle)
+ return
+
+
+ def needNewJacobian(self):
+ """
+ Returns true if we need to recompute Jacobian matrix for operator,
+ false otherwise.
+ """
+ assert(None != self.cppHandle)
+ return self.cppHandle.needNewJacobian
+
+
+ def integrateJacobian(self, jacobian, fields):
+ """
+ Integrate contributions to Jacobian term at time t.
+ """
+ assert(None != self.cppHandle)
+ self.cppHandle.integrateJacobian(jacobian, fields, self.mesh.cppHandle)
+ return
+
+
+ def updateState(self, field):
+ """
+ Update state variables as needed.
+ """
+ assert(None != self.cppHandle)
+ self.cppHandle.updateState(field, self.mesh.cppHandle)
+ return
+
+
# End of file
Added: short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/IntegratorElasticity.py
+##
+## @brief Python object implementing sgeneral methods for time
+## integration of the elasticity equation using finite-elements.
+##
+## Factory: integrator
+
+from Integrator import Integrator
+
+# IntegratorElasticity class
+class IntegratorElasticity(Integrator):
+ """
+ Python object implementing sgeneral methods for time integration of
+ the elasticity equation using finite-elements.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="integratorelasticity"):
+ """
+ Constructor.
+ """
+ Integrator.__init__(self)
+ import journal
+ self._info = journal.info(name)
+ return
+
+
+ def initMaterial(self, mesh, material):
+ """
+ Initialize material properties.
+ """
+ self._info.log("Initializing integrator material '%s'." % material.label)
+ material.initialize(mesh)
+ self.material = material
+ self.cppHandle.material = self.material.cppHandle
+ return
+
+
+# End of file
Deleted: short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,73 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file pylith/feassemble/IntegratorExplicit.py
-##
-## @brief Python object for explicit time integration of actions with
-## finite-elements.
-##
-## Factory: integrator
-
-from Integrator import Integrator
-
-# IntegratorInertia class
-class IntegratorExplicit(Integrator):
- """
- Python object for explicit integration of operator actions with
- finite-elements.
-
- Factory: integrator.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="integratorexplicit"):
- """
- Constructor.
- """
- Integrator.__init__(self, name)
- return
-
-
- def timeStep(self, t):
- """
- Set time step for advancing from time t to time t+dt.
- """
- self.cppHandle.timeStep = t.value
- return
-
-
- def stableTimeStep(self):
- """
- Get stable time step for advancing from time t to time t+dt.
- """
- return self.cppHandle.getStableTimeStep()
-
-
- def integrateConstant(self, fieldOut, fieldInT, fieldInTmdt):
- """
- Integrate constant term for dynamic terms for finite-elements.
- """
- self.cppHandle.integrateConstant(fieldOut, fieldInT, fieldInTmdt,
- self.mesh.cppHandle)
- return
-
-
- def integrateJacobian(self, jacobian, fieldInT):
- """
- Integrate Jacobian term for dynamic terms for finite-elements.
- """
- self.cppHandle.integrateJacobian(jacobian, fieldInT, self.mesh.cppHandle)
- return
-
-
-# End of file
Deleted: short/3D/PyLith/trunk/pylith/feassemble/IntegratorImplicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/IntegratorImplicit.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/IntegratorImplicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,72 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file pylith/feassemble/IntegratorImplicit.py
-##
-## @brief Python object for implicit time integration of actions with
-## finite-elements.
-##
-## Factory: integrator
-
-from Integrator import Integrator
-
-# IntegratorInertia class
-class IntegratorImplicit(Integrator):
- """
- Python object for implicit integration of operator actions with
- finite-elements.
-
- Factory: integrator.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="integratorimplicit"):
- """
- Constructor.
- """
- Integrator.__init__(self, name)
- return
-
-
- def timeStep(self, t):
- """
- Set time step for advancing from time t to time t+dt.
- """
- self.cppHandle.timeStep = t.value
- return
-
-
- def stableTimeStep(self):
- """
- Get stable time step for advancing from time t to time t+dt.
- """
- return self.cppHandle.getStableTimeStep()
-
-
- def integrateResidual(self, fieldOut, fieldInT):
- """
- Integrate residual term for quasi-static terms for finite-elements.
- """
- self.cppHandle.integrateResidual(fieldOut, fieldInT, self.mesh.cppHandle)
- return
-
-
- def integrateJacobian(self, jacobian, fieldInT):
- """
- Integrate Jacobian term for quasi-static terms for finite-elements.
- """
- self.cppHandle.integrateJacobian(jacobian, fieldInT, self.mesh.cppHandle)
- return
-
-
-# End of file
Modified: short/3D/PyLith/trunk/pylith/feassemble/__init__.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/__init__.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/feassemble/__init__.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -14,14 +14,12 @@
##
## @brief Python PyLith finite-element assembler module initialization
-__all__ = ['ExplicitElasticity',
+__all__ = ['ElasticityExplicit',
+ 'ElasticityImplicit',
'FIATCell',
'FIATLagrange',
'FIATSimplex',
- 'ImplicitElasticity',
'Integrator',
- 'IntegratorExplicit',
- 'IntegratorImplicit',
'ReferenceCell',
'quadrature']
Deleted: short/3D/PyLith/trunk/pylith/problems/EqDeformation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/EqDeformation.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/problems/EqDeformation.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,100 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file pylith/problems/EqDeformation.py
-##
-## @brief Python EqDeformation for computing deformation associated
-## with earthquakes.
-##
-## Factory: problem.
-
-from TimeDependent import TimeDependent
-
-# EqDeformation class
-class EqDeformation(TimeDependent):
- """
- Python EqDeformation for computing deformation associated with
- earthquakes.
-
- Factory: problem.
- """
-
- # INVENTORY //////////////////////////////////////////////////////////
-
- class Inventory(TimeDependent.Inventory):
- """
- Python object for managing EqDeformation facilities and properties.
- """
-
- ## @class Inventory
- ## Python object for managing EqDeformation facilities and properties.
- ##
- ## \b Properties
- ## @li None
- ##
- ## \b Facilities
- ## @li \b faults Faults or interior slip surfaces.
-
- import pyre.inventory
-
- from pylith.faults.FaultsBin import FaultsBin
- faults = pyre.inventory.facility("faults", family="faults",
- factory=FaultsBin)
- faults.meta['tip'] = "Faults or interior slip surfaces."
-
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="eqdeformation"):
- """
- Constructor.
- """
- TimeDependent.__init__(self, name)
- return
-
-
- def checkpoint(self):
- """
- Save problem state for restart.
- """
- TimeDependent.checkpoint() # Save state of parent
-
- # Save state of this object
- raise NotImplementedError, \
- "EqDeformation::checkpoint() not implemented."
-
- # Save state of children
- self.faults.checkpoint()
- return
-
-
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _configure(self):
- """
- Set members based using inventory.
- """
- TimeDependent._configure(self)
- self.faults = self.inventory.faults
- return
-
-
-# FACTORIES ////////////////////////////////////////////////////////////
-
-def problem():
- """
- Factory associated with EqDeformation.
- """
- return EqDeformation()
-
-
-# End of file
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -38,25 +38,6 @@
Factory: pde_formulation.
"""
- # INVENTORY //////////////////////////////////////////////////////////
-
- class Inventory(Formulation.Inventory):
- """
- Python object for managing Explicit facilities and properties.
- """
-
- ## @class Inventory
- ## Python object for managing Explicit facilities and properties.
- ##
- ## \b Properties
- ## @li None
- ##
- ## \b Facilities
- ## @li None
-
- import pyre.inventory
-
-
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="explicit"):
@@ -71,44 +52,42 @@
"""
Get integrator for elastic material.
"""
- from pylith.feassemble.ExplicitElasticity import ExplicitElasticity
- return ExplicitElasticity()
+ from pylith.feassemble.ElasticityExplicit import ElasticityExplicit
+ return ElasticityExplicit()
- def initialize(self, mesh, materials, boundaryConditions, dimension, dt):
+ def initialize(self, mesh, materials, boundaryConditions,
+ interfaceConditions, dimension, dt):
"""
Initialize problem for explicit time integration.
"""
- self.integrators = []
Formulation.initialize(self, mesh, materials, boundaryConditions,
- dimension, dt)
+ interfaceConditions, dimension, dt)
- self._info.log("Initializing integrators.")
-
self._info.log("Creating fields and matrices.")
- self.dispT = mesh.createRealSection("dispT", dimension)
- self.dispTmdt = mesh.createRealSection("dispTmdt", dimension)
- self.dispTpdt = mesh.createRealSection("dispTpdt", dimension)
- self.constant = mesh.createRealSection("constant", dimension)
+ self.fields.addReal("dispT")
+ self.fields.addReal("dispTmdt")
+ self.fields.addReal("dispTpdt")
+ self.fields.addReal("residual")
+ self.fields.createHistory(["dispTpdt", "dispT", "dispTmdt"])
+ self.fields.setFiberDimension("dispT", dimension)
+ for constraint in self.constraints:
+ constraint.setConstraintSizes(self.fields.getReal("dispT"), mesh)
+ self.fields.allocate("dispT")
+ for constraint in self.constraints:
+ constraint.setConstraints(self.fields.getReal("dispT"), mesh)
+ self.fields.copyLayout("dispT")
+ self.jacobian = mesh.createMatrix(self.fields.getReal("residual"))
- # Setup constraints
- # STUFF GOES HERE
-
- mesh.allocateRealSection(self.dispT)
- mesh.allocateRealSection(self.dispTmdt)
- mesh.allocateRealSection(self.dispTpdt)
- mesh.allocateRealSection(self.constant)
-
- self.jacobian = mesh.createMatrix(self.constant)
-
- self._info.log("Integrating Jacobian of operator.")
+ self._info.log("Forming Jacobian of operator.")
+ import pylith.utils.petsc as petsc
+ #petsc.zeroMatrix(self.jacobian)
for integrator in self.integrators:
integrator.timeStep(dt)
- integrator.integrateJacobian(self.jacobian, self.dispT)
- import pylith.utils.petsc as petsc
+ integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
petsc.mat_assemble(self.jacobian)
- self.solver.initialize(mesh, self.dispTpdt)
+ self.solver.initialize(mesh, self.fields.getReal("dispTpdt"))
return
@@ -126,7 +105,22 @@
"""
Hook for doing stuff before advancing time step.
"""
- self._info.log("WARNING: Explicit::prestep() not implemented.")
+ dispTpdt = self.fields.getReal("dispTpdt")
+ for constraint in self.constraints:
+ constraint.setField(dispTpdt, t+dt, dt)
+
+ needNewJacobian = False
+ for integrator in self.integrators:
+ if integrator.needNewJacobian():
+ needNewJacobian = True
+ if needNewJacobian:
+ self._info.log("Reforming Jacobian of operator.")
+ import pylith.utils.petsc as petsc
+ petsc.zeroMatrix(self.jacobian)
+ for integrator in self.integrators:
+ integrator.timeStep(dt)
+ integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+ petsc.mat_assemble(self.jacobian)
return
@@ -135,14 +129,15 @@
Advance to next time step.
"""
self._info.log("Integrating constant term in operator.")
+ residual = self.fields.getReal("residual")
import pylith.topology.topology as bindings
- bindings.zeroRealSection(self.constant)
+ bindings.zeroRealSection(residual)
for integrator in self.integrators:
integrator.timeStep(dt)
- integrator.integrateConstant(self.constant, self.dispT, self.dispTmdt)
+ integrator.integrateResidual(residual, self.fields.cppHandle)
self._info.log("Solving equations.")
- self.solver.solve(self.dispTpdt, self.jacobian, self.constant)
+ self.solver.solve(self.fields.getReal("dispTpdt"), self.jacobian, residual)
return
@@ -150,10 +145,11 @@
"""
Hook for doing stuff after advancing time step.
"""
- tmp = self.dispTmdt
- self.dispTmdt = self.dispT
- self.dispT = self.dispTpdt
- self.dispTpdt = tmp
+ self.fields.shiftHistory()
+
+ self._info.log("Updating integrators states.")
+ for integrator in self.integrators:
+ integrator.updateState(self.fields.getReal("dispT"))
return
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -61,16 +61,25 @@
"""
Component.__init__(self, name, facility="pde_formulation")
self.integrators = None
+ self.constraints = None
+ self.fields = None
return
- def initialize(self, mesh, materials, boundaryConditions, dimension, dt):
+ def initialize(self, mesh, materials, boundaryConditions,
+ interfaceConditions, dimension, dt):
"""
Create integrators for each element family.
"""
- self._info.log("Initializing integrators.")
+ from pylith.feassemble.Integrator import implementsIntegrator
+ from pylith.feassemble.Constraint import implementsConstraint
+
+ from pylith.topology.FieldsManager import FieldsManager
+ self.fields = FieldsManager(mesh)
self.integrators = []
+ self.constraints = []
+ self._info.log("Initializing materials.")
for material in materials.materials:
if material.quadrature.spaceDim != dimension:
raise ValueError, \
@@ -78,15 +87,38 @@
"for material '%s' is for spatial dimension '%d'." % \
(dimension, material.label, material.quadrature.spaceDim)
integrator = self.elasticityIntegrator()
+ if not implementsIntegrator(integrator):
+ raise TypeError, \
+ "Could not use '%s' as an integrator for material '%s'. " \
+ "Functionality missing." % (integrator.name, material.label)
integrator.setMesh(mesh)
integrator.initQuadrature(material.quadrature)
integrator.initMaterial(mesh, material)
self.integrators.append(integrator)
+ self._info.log("Initializing boundary conditions.")
for bc in boundaryConditions.bc:
- integrator = bc
- integrator.initialize(mesh)
- self.integrators.append(integrator)
+ bc.initialize(mesh)
+ if implementsIntegrator(bc):
+ self.integrators.append(bc)
+ elif implementsConstraint(bc):
+ self.constraints.append(bc)
+ else:
+ raise TypeError, \
+ "Could not determine whether boundary condition '%s' is an " \
+ "integrator or a constraint." % bc.name
+
+ self._info.log("Initializing interior interfaces.")
+ for ic in interfaceConditions.ic:
+ ic.initialize(mesh)
+ if implementsIntegrator(ic):
+ self.integrators.append(ic)
+ elif implementsConstraint(ic):
+ self.constraints.append(ic)
+ else:
+ raise TypeError, \
+ "Could not determine whether interface condition '%s' is an " \
+ "integrator or a constraint." % ic.name
return
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -75,45 +75,36 @@
"""
Get integrator for elastic material.
"""
- from pylith.feassemble.ImplicitElasticity import ImplicitElasticity
- return ImplicitElasticity()
+ from pylith.feassemble.ElasticityImplicit import ElasticityImplicit
+ return ElasticityImplicit()
- def initialize(self, mesh, materials, boundaryConditions, dimension, dt):
+ def initialize(self, mesh, materials, boundaryConditions,
+ interfaceConditions, dimension, dt):
"""
Initialize problem for implicit time integration.
"""
- self.integrators = []
Formulation.initialize(self, mesh, materials, boundaryConditions,
- dimension, dt)
+ interfaceConditions, dimension, dt)
- self._info.log("Initializing integrators.")
-
self._info.log("Creating fields and matrices.")
- self.dispT = mesh.createRealSection("dispT", dimension)
- self.dispTBctpdt = mesh.createRealSection("dispTBctpdt", dimension)
- self.dispIncrement = mesh.createRealSection("dispIncrement", dimension)
- self.residual = mesh.createRealSection("residual", dimension)
+ self.fields.addReal("dispT")
+ self.fields.addReal("dispTBctpdt")
+ self.fields.addReal("dispIncr")
+ self.fields.addReal("residual")
+ self.fields.createHistory(["dispTBctpdt", "dispT"])
+ self.fields.setFiberDimension("dispT", dimension)
+ for constraint in self.constraints:
+ constraint.setConstraintSizes(self.fields.getReal("dispT"), mesh)
+ self.fields.allocate("dispT")
+ for constraint in self.constraints:
+ constraint.setConstraints(self.fields.getReal("dispT"), mesh)
+ self.fields.copyLayout("dispT")
+ self.jacobian = mesh.createMatrix(self.fields.getReal("dispT"))
- # Setup constraints
- # STUFF GOES HERE
+ self._solveElastic(mesh, t=0.0, dt=dt)
- mesh.allocateRealSection(self.dispT)
- mesh.allocateRealSection(self.dispTBctpdt)
- mesh.allocateRealSection(self.dispIncrement)
- mesh.allocateRealSection(self.residual)
-
- self.jacobian = mesh.createMatrix(self.residual)
-
- self._info.log("Integrating Jacobian of operator.")
- for integrator in self.integrators:
- integrator.timeStep(dt)
- integrator.integrateJacobian(self.jacobian, self.dispTBctpdt)
- import pylith.utils.petsc as petsc
- petsc.mat_assemble(self.jacobian)
-
- self.solver.initialize(mesh, self.dispIncrement)
- # self.mesh = mesh
+ self.solver.initialize(mesh, self.fields.getReal("dispIncr"))
return
@@ -131,17 +122,24 @@
"""
Hook for doing stuff before advancing time step.
"""
- # This will need to set dispTBctpdt to the BC at time step t+dt.
- # Non-constrained DOF are unaffected and will be equal to their
- # values from time step t.
- # In this routine I also need to integrate the tractions for step
- # t+dt, but I don't think the function for this is available yet.
- # from pylith.bc.Dirichlet import Dirichlet
+ # Set dispTBctpdt to the BC at time t+dt. Unconstrained DOF are
+ # unaffected and will be equal to their values at time t.
+ dispTBctpdt = self.fields.getReal("dispTBctpdt")
+ for constraint in self.constraints:
+ constraint.setField(dispTBctpdt, t+dt, dt)
- # dispbc = Dirichlet
-
- # dispbc.setField(t+dt, self.dispTBctpdt, self.mesh)
- self._info.log("WARNING: Implicit::prestep() not implemented.")
+ needNewJacobian = False
+ for integrator in self.integrators:
+ if integrator.needNewJacobian():
+ needNewJacobian = True
+ if needNewJacobian:
+ self._info.log("Reforming Jacobian of operator.")
+ import pylith.utils.petsc as petsc
+ petsc.zeroMatrix(self.jacobian)
+ for integrator in self.integrators:
+ integrator.timeStep(dt)
+ integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+ petsc.mat_assemble(self.jacobian)
return
@@ -150,14 +148,15 @@
Advance to next time step.
"""
self._info.log("Integrating residual term in operator.")
+ residual = self.fields.getReal("residual")
import pylith.topology.topology as bindings
- bindings.zeroRealSection(self.residual)
+ bindings.zeroRealSection(residual)
for integrator in self.integrators:
integrator.timeStep(dt)
- integrator.integrateResidual(self.residual, self.dispTBctpdt)
+ integrator.integrateResidual(residual, self.fields.cppHandle)
self._info.log("Solving equations.")
- self.solver.solve(self.dispIncrement, self.jacobian, self.residual)
+ self.solver.solve(self.fields.getReal("dispIncr"), self.jacobian, residual)
return
@@ -165,16 +164,21 @@
"""
Hook for doing stuff after advancing time step.
"""
- # This should give us the total displacements for time step t+dt, which
- # is renamed as time step t following the solve.
- # The vector dispTBctpdt contains the displacements from time step t
- # along with the displacement BC from time step t+dt. The displacement
- # increments computed from the residual are then added to this to give us
- # the total displacement field at time t+dt.
+ # This should give us the total displacements after stepping
+ # forward from t to t+dt, which is renamed as time step t for the
+ # next solve. The field dispTBctpdt contains the displacements
+ # from time step t along with the displacement BC from time step
+ # t+dt. The displacement increments computed from the residual
+ # are then added to this to give us the total displacement field
+ # at time t+dt.
+
# Need a real way to do the operation below.
- # It is commented out for now, and should be replaced with a call to PETSc.
- # self.dispT = self.dispTBctpdt+self.dispIncrement
- self.dispTBctpdt = self.dispT
+ # self.dispT = self.dispTBctpdt + self.dispIncr
+ self.fields.shiftHistory()
+
+ self._info.log("Updating integrators states.")
+ for integrator in self.integrators:
+ integrator.updateState(self.fields.getReal("dispT"))
return
@@ -188,6 +192,34 @@
return
+ def _solveElastic(self, mesh, t, dt):
+ """
+ Solve for elastic solution.
+ """
+ self._info.log("Computing elastic solution.")
+
+ self._info.log("Setting constraints.")
+ dispT = self.fields.getReal("dispT")
+ for constraint in self.constraints:
+ constraint.setField(dispT, t, dt)
+
+ self._info.log("Integrating Jacobian and residual of operator.")
+ for integrator in self.integrators:
+ integrator.timeStep(dt)
+ integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+ integrator.integrateResidual(self.fields.getReal("dispT"),
+ self.fields.cppHandle)
+ import pylith.utils.petsc as petsc
+ petsc.mat_assemble(self.jacobian)
+
+ self.solver.initialize(mesh, dispT)
+
+ self._info.log("Solving equations.")
+ self.solver.solve(dispT, self.jacobian, self.fields.getReal("residual"))
+
+ return
+
+
# FACTORIES ////////////////////////////////////////////////////////////
def pde_formulation():
Modified: short/3D/PyLith/trunk/pylith/problems/Problem.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Problem.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/problems/Problem.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -42,6 +42,8 @@
## \b Facilities
## @li \b materials Materials in problem.
## @li \b bc Boundary conditions.
+ ## @li \b interfaces Interior surfaces with constraints or
+ ## constitutive models.
import pyre.inventory
@@ -51,11 +53,17 @@
materials.meta['tip'] = "Materials in problem."
from BoundaryConditions import BoundaryConditions
- bc = pyre.inventory.facility("bc", family="bc",
+ bc = pyre.inventory.facility("bc", family="boundary_conditions",
factory=BoundaryConditions)
bc.meta['tip'] = "Boundary conditions."
-
+
+ from pylith.faults.FaultsBin import FaultsBin
+ interfaces = pyre.inventory.facility("interfaces", family="interfaces_bin",
+ factory=FaultsBin)
+ interfaces.meta['tip'] = "Interior surfaces with constraints or " \
+ "constitutive models."
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="problem"):
@@ -101,6 +109,7 @@
Component._configure(self)
self.materials = self.inventory.materials
self.bc = self.inventory.bc
+ self.interfaces = self.inventory.interfaces
return
Modified: short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/TimeDependent.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/pylith/problems/TimeDependent.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -90,9 +90,15 @@
bc/quadrature, etc.).
"""
self._info.log("Initializing problem.")
+
+ if self.dimension != mesh.dimension():
+ raise ValueError, \
+ "Spatial dimension of problem is '%d' but mesh contains cells " \
+ "for spatial dimension '%d'." % \
+ (self.dimension, mesh.dimension)
self.mesh = mesh
self.formulation.initialize(mesh, self.materials, self.bc,
- self.dimension, self.dt)
+ self.interfaces, self.dimension, self.dt)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -145,30 +145,6 @@
} // testAdjustTopologyHex8Lagrange
// ----------------------------------------------------------------------
-// Test _orientationSize().
-void
-pylith::faults::TestFaultCohesive::testOrientationSize(void)
-{ // testOrientationSize
- const int cellDim = 2;
- const int spaceDim = 3;
- const int numBasis = 1;
- const int numQuadPts = 1;
- const double basis[] = { 0.5, 0.5, 0.4 };
- const double basisDeriv[] = { 0.5, 0.3, -0.4 };
- const double quadPtsRef[] = { 0.0, 3.0 };
- const double quadWts[] = { 2.0 };
- const double minJacobian = 1.0;
-
- feassemble::Quadrature2Din3D q;
- q.initialize(basis, basisDeriv, quadPtsRef, quadWts,
- cellDim, numBasis, numQuadPts, spaceDim);
-
- FaultCohesiveKin fault;
- fault.quadrature(&q);
- CPPUNIT_ASSERT_EQUAL(cellDim*spaceDim, fault._orientationSize());
-} // testOrientationSize
-
-// ----------------------------------------------------------------------
// Test _orient1D().
void
pylith::faults::TestFaultCohesive::testOrient1D(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -53,8 +53,6 @@
CPPUNIT_TEST( testAdjustTopologyTet4Lagrange );
CPPUNIT_TEST( testAdjustTopologyHex8Lagrange );
- CPPUNIT_TEST( testOrientationSize );
-
CPPUNIT_TEST( testOrient1D );
CPPUNIT_TEST( testOrient2D );
CPPUNIT_TEST( testOrient3D );
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -19,25 +19,5 @@
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesiveKin );
-// ----------------------------------------------------------------------
-// Test clone()
-void
-pylith::faults::TestFaultCohesiveKin::testClone(void)
-{ // testClone
- const int id = 65;
- const std::string label = "fault ABC";
- FaultCohesiveKin fault;
- fault.id(id);
- fault.label(label.c_str());
-
- Fault* faultCopy = fault.clone();
- CPPUNIT_ASSERT(0 != faultCopy);
-
- CPPUNIT_ASSERT_EQUAL(id, faultCopy->id());
- CPPUNIT_ASSERT_EQUAL(label, faultCopy->label());
- delete faultCopy; faultCopy = 0;
-} // testClone
-
-
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -36,15 +36,11 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestFaultCohesiveKin );
- CPPUNIT_TEST( testClone );
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
- /// Test clone()
- void testClone(void);
-
}; // class TestFaultCohesiveKin
#endif // pylith_faults_testfaultcohesivekin_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc 2007-05-25 20:57:26 UTC (rev 6970)
@@ -14,7 +14,7 @@
#include "TestIntegrator.hh" // Implementation of class methods
-#include "pylith/feassemble/Integrator.hh" // USES Integrator
+#include "pylith/feassemble/ElasticityExplicit.hh" // USES ElasticityExplicit
#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
#include <petscmat.h>
@@ -23,24 +23,6 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestIntegrator );
// ----------------------------------------------------------------------
-// Test clone().
-void
-pylith::feassemble::TestIntegrator::testCopy(void)
-{ // testCopy
- // Test copy constructor by testing value of minJacobian value in quadrature
-
- Quadrature1D quadrature;
- const double minJacobian = 4.0;
- quadrature.minJacobian(minJacobian);
-
- Integrator iOrig;
- iOrig.quadrature(&quadrature);
-
- Integrator iCopy = Integrator(iOrig);
- CPPUNIT_ASSERT_EQUAL(minJacobian, iCopy._quadrature->minJacobian());
-} // testCopy
-
-// ----------------------------------------------------------------------
// Test quadrature().
void
pylith::feassemble::TestIntegrator::testQuadrature(void)
@@ -52,7 +34,7 @@
const double minJacobian = 4.0;
quadrature.minJacobian(minJacobian);
- Integrator integrator;
+ ElasticityExplicit integrator;
integrator.quadrature(&quadrature);
CPPUNIT_ASSERT_EQUAL(minJacobian, integrator._quadrature->minJacobian());
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh 2007-05-25 20:57:26 UTC (rev 6970)
@@ -38,7 +38,6 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestIntegrator );
- CPPUNIT_TEST( testCopy );
CPPUNIT_TEST( testQuadrature );
CPPUNIT_TEST_SUITE_END();
@@ -46,9 +45,6 @@
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
- /// Test copy constructor.
- void testCopy(void);
-
/// Test quadrature()
void testQuadrature(void);
Copied: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py (from rev 6964, short/3D/PyLith/trunk/unittests/pytests/feassemble/TestExplicitElasticity.py)
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestExplicitElasticity.py 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/feassemble/TestElasticityExplicit.py
+
+## @brief Unit testing of Python ElasticityExplicit object.
+
+import unittest
+from pylith.feassemble.ElasticityExplicit import ElasticityExplicit
+
+# ----------------------------------------------------------------------
+class TestElasticityExplicit(unittest.TestCase):
+ """
+ Unit testing of Python ElasticityExplicit object.
+ """
+
+ def test_initQuadrature(self):
+ """
+ Test initQuadrature().
+ """
+ from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
+ q = Quadrature2D()
+ minJacobian = 4.0e-02;
+ q.minJacobian = minJacobian
+ from pylith.feassemble.FIATSimplex import FIATSimplex
+ q.cell = FIATSimplex()
+ q.cell.shape = "triangle"
+ q.cell.order = 1
+ q.cell.degree = 1
+
+ integrator = ElasticityExplicit()
+ integrator.initQuadrature(q)
+ self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
+ return
+
+
+# End of file
Copied: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py (from rev 6964, short/3D/PyLith/trunk/unittests/pytests/feassemble/TestImplicitElasticity.py)
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestImplicitElasticity.py 2007-05-25 03:32:59 UTC (rev 6964)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/feassemble/TestElasticityImplicit.py
+
+## @brief Unit testing of Python ElasticityImplicit object.
+
+import unittest
+from pylith.feassemble.ElasticityImplicit import ElasticityImplicit
+
+# ----------------------------------------------------------------------
+class TestElasticityImplicit(unittest.TestCase):
+ """
+ Unit testing of Python ElasticityImplicit object.
+ """
+
+ def test_initQuadrature(self):
+ """
+ Test initQuadrature().
+ """
+ from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
+ q = Quadrature2D()
+ minJacobian = 4.0e-02;
+ q.minJacobian = minJacobian
+ from pylith.feassemble.FIATSimplex import FIATSimplex
+ q.cell = FIATSimplex()
+ q.cell.shape = "triangle"
+ q.cell.order = 1
+ q.cell.degree = 1
+
+ integrator = ElasticityImplicit()
+ integrator.initQuadrature(q)
+ self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
+ return
+
+
+# End of file
Deleted: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestExplicitElasticity.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestExplicitElasticity.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestExplicitElasticity.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,46 +0,0 @@
-#!/usr/bin/env python
-#
-# ======================================================================
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# {LicenseText}
-#
-# ======================================================================
-#
-
-## @file unittests/pytests/feassemble/TestExplicitElasticity.py
-
-## @brief Unit testing of Python ExplicitElasticity object.
-
-import unittest
-from pylith.feassemble.ExplicitElasticity import ExplicitElasticity
-
-# ----------------------------------------------------------------------
-class TestExplicitElasticity(unittest.TestCase):
- """
- Unit testing of Python ExplicitElasticity object.
- """
-
- def test_initQuadrature(self):
- """
- Test initQuadrature().
- """
- from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
- q = Quadrature2D()
- minJacobian = 4.0e-02;
- q.minJacobian = minJacobian
- from pylith.feassemble.FIATSimplex import FIATSimplex
- q.cell = FIATSimplex()
- q.cell.shape = "triangle"
- q.cell.order = 1
- q.cell.degree = 1
-
- integrator = ExplicitElasticity()
- integrator.initQuadrature(q)
- self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
- return
-
-
-# End of file
Deleted: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestImplicitElasticity.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestImplicitElasticity.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestImplicitElasticity.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -1,46 +0,0 @@
-#!/usr/bin/env python
-#
-# ======================================================================
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# {LicenseText}
-#
-# ======================================================================
-#
-
-## @file unittests/pytests/feassemble/TestImplicitElasticity.py
-
-## @brief Unit testing of Python ImplicitElasticity object.
-
-import unittest
-from pylith.feassemble.ImplicitElasticity import ImplicitElasticity
-
-# ----------------------------------------------------------------------
-class TestImplicitElasticity(unittest.TestCase):
- """
- Unit testing of Python ImplicitElasticity object.
- """
-
- def test_initQuadrature(self):
- """
- Test initQuadrature().
- """
- from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
- q = Quadrature2D()
- minJacobian = 4.0e-02;
- q.minJacobian = minJacobian
- from pylith.feassemble.FIATSimplex import FIATSimplex
- q.cell = FIATSimplex()
- q.cell.shape = "triangle"
- q.cell.order = 1
- q.cell.degree = 1
-
- integrator = ImplicitElasticity()
- integrator.initQuadrature(q)
- self.assertEqual(minJacobian, integrator.quadrature.minJacobian)
- return
-
-
-# End of file
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py 2007-05-25 20:51:55 UTC (rev 6969)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py 2007-05-25 20:57:26 UTC (rev 6970)
@@ -59,11 +59,11 @@
from TestIntegrator import TestIntegrator
suite.addTest(unittest.makeSuite(TestIntegrator))
- from TestExplicitElasticity import TestExplicitElasticity
- suite.addTest(unittest.makeSuite(TestExplicitElasticity))
+ from TestElasticityExplicit import TestElasticityExplicit
+ suite.addTest(unittest.makeSuite(TestElasticityExplicit))
- from TestImplicitElasticity import TestImplicitElasticity
- suite.addTest(unittest.makeSuite(TestImplicitElasticity))
+ from TestElasticityImplicit import TestElasticityImplicit
+ suite.addTest(unittest.makeSuite(TestElasticityImplicit))
return suite
More information about the cig-commits
mailing list