[cig-commits] r6699 - in short/3D/PyLith/trunk: . libsrc/faults
libsrc/feassemble modulesrc/faults
brad at geodynamics.org
brad at geodynamics.org
Thu Apr 26 15:03:46 PDT 2007
Author: brad
Date: 2007-04-26 15:03:45 -0700 (Thu, 26 Apr 2007)
New Revision: 6699
Modified:
short/3D/PyLith/trunk/TODO
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/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
Log:
Setup basic layout of fault implementation for kinematic source.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/TODO 2007-04-26 22:03:45 UTC (rev 6699)
@@ -4,10 +4,14 @@
Cleanup
- Add check in MeshIO 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.
+ 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.
Use pylith/utils/sievefwd.hh and pylith/utils/sievetypes.hh in unit
tests. Remove typedefs when possible (take advantage of typedefs in
@@ -16,22 +20,32 @@
Use MeshIOAscii to get mesh information into Python tests?
0. Start implementing faults
- a. Add creation of cohesive elements
- Add tests for interpolated meshes
+ a. Add creation of cohesive cells
+ i. Add tests for interpolated meshes and constrained cohesive cells.
b. Implement slip time function
- Similar to material. Set up parameters using spatial databases.
-
+ i. Add Python
+ ii. Add unit tests
c. Start implementing integrator for faults
- Initialize - calculate orientation
- How do we create kinematic source
+ i. Update Reference cell
+ (1) Compute basis and basisDeriv at vertices
+ (2) Add unit tests
+ ii. Update Quadrature
+ (1) Add basis and basisDeriv at vertices to Quadrature
+ (2) Add computeGeometryQuad/computeGeometryVert to Quadrature
+ (3) Add unit tests
+ iii. FaultCohesive
+ (1) orientXD
+ (2) Unit tests
+ iv. FaultCohesiveKin
+ (1) initialize(), setField(), integrateResidual(), integrateJacobian()
+ (2) Add Python unit tests
+ (3) Add C++ unit tests
- d. Create unit tests at C++ level
- e. Create unit tests at Python level
-
1. Finish implementing ExplicitElasticity and Explicit
- a. Double check loops for calcConstant() and calcJacobian()
- b. Create unit test (construction of constant term and Jacobian)
+ a. Replace integrateConstant() with integrateResidual()
+ 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
@@ -115,7 +129,7 @@
Outputs:
* PETSc Mesh on each processor (restricted to processor)
- Note: Called from topology.Mesh.py (need Python bindings for Sieve function).
+ [DONE]
4. Global refinement of mesh.
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2007-04-26 22:03:45 UTC (rev 6699)
@@ -84,7 +84,7 @@
* @param mesh PETSc mesh
*/
virtual
- void adjustTopology(ALE::Obj<ALE::Mesh>* mesh) = 0;
+ void adjustTopology(const ALE::Obj<ALE::Mesh>& mesh) = 0;
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
@@ -95,7 +95,7 @@
* be up-dip direction).
*/
virtual
- void initialize(ALE::Obj<ALE::Mesh>* mesh,
+ void initialize(const ALE::Obj<ALE::Mesh>& mesh,
const double_array& upDir) = 0;
// PROTECTED METHODS //////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-04-26 22:03:45 UTC (rev 6699)
@@ -17,11 +17,8 @@
#include "CohesiveTopology.hh" // USES CohesiveTopology::create()
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
-#include "pylith/utils/array.hh" // USES double_array
#include <assert.h> // USES assert()
-#include <sstream> // USES std::ostringstream
-#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
// Default constructor.
@@ -50,35 +47,17 @@
// ----------------------------------------------------------------------
// Adjust mesh topology for fault implementation.
void
-pylith::faults::FaultCohesive::adjustTopology(ALE::Obj<ALE::Mesh>* mesh)
+pylith::faults::FaultCohesive::adjustTopology(const ALE::Obj<ALE::Mesh>& mesh)
{ // adjustTopology
- assert(0 != mesh);
assert("" != label());
// Get group of vertices associated with fault
const ALE::Obj<int_section_type>& groupField =
- (*mesh)->getIntSection(label());
+ mesh->getIntSection(label());
assert(!groupField.isNull());
- CohesiveTopology::create(_faultMesh, *mesh, groupField, id());
+ CohesiveTopology::create(_faultMesh, mesh, groupField, id());
} // adjustTopology
-// ----------------------------------------------------------------------
-// Initialize fault. Determine orientation and setup boundary
-void
-pylith::faults::FaultCohesive::initialize(ALE::Obj<ALE::Mesh>* mesh,
- const double_array& upDir)
-{ // initialize
- assert(0 != mesh);
- assert(0 != _quadrature);
-
- if (3 != upDir.size())
- throw std::runtime_error("Up direction for fault orientation must be "
- "a vector with 3 components.");
-
-
-
-} // initialize
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-04-26 22:03:45 UTC (rev 6699)
@@ -23,6 +23,7 @@
#include "pylith/feassemble/Integrator.hh" // ISA Integrator
#include "pylith/utils/sievefwd.hh" // HOLDSA PETSc Mesh
+#include "pylith/utils/petscfwd.h" // USES PetscMat
/// Namespace for pylith package
namespace pylith {
@@ -53,7 +54,7 @@
*
* @param mesh PETSc mesh
*/
- void adjustTopology(ALE::Obj<ALE::Mesh>* mesh);
+ void adjustTopology(const ALE::Obj<ALE::Mesh>& mesh);
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
@@ -64,9 +65,40 @@
* be up-dip direction).
*/
virtual
- void initialize(ALE::Obj<ALE::Mesh>* mesh,
- const double_array& upDir);
+ void initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const double_array& upDir) = 0;
+ /** 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 disp Displacement field
+ * @param mesh Finite-element mesh
+ */
+ virtual
+ void setField(const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh) = 0;
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -76,6 +108,73 @@
*/
FaultCohesive(const FaultCohesive& m);
+ /** Compute weighted orientation of fault for cohesive cell between
+ * 1-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = numPts*spaceDim*spaceDim
+ * index = iPt*spaceDim*spaceDim + iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ void orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir);
+
+ /** Compute weighted orientation of fault for cohesive cell between
+ * 2-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = numPts*spaceDim*spaceDim
+ * index = iPt*spaceDim*spaceDim + iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ void orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir);
+
+ /** Compute weighted orientation of fault for cohesive cell between
+ * 3-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = numPts*spaceDim*spaceDim
+ * index = iPt*spaceDim*spaceDim + iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ void orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir);
+
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-04-26 22:03:45 UTC (rev 6699)
@@ -16,6 +16,12 @@
#include "EqKinSrc.hh" // USES EqKinSrc
+#include "pylith/utils/array.hh" // USES double_array
+
+#include <assert.h> // USES assert()
+#include <sstream> // USES std::ostringstream
+#include <stdexcept> // USES std::runtime_error
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void) :
@@ -48,5 +54,50 @@
delete _eqsrc; _eqsrc = (0 != src) ? src->clone() : 0;
} // eqsrc
+// ----------------------------------------------------------------------
+// Initialize fault. Determine orientation and setup boundary
+void
+pylith::faults::FaultCohesiveKin::initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const double_array& upDir)
+{ // initialize
+ assert(0 != _quadrature);
+
+ if (3 != upDir.size())
+ throw std::runtime_error("Up direction for fault orientation must be "
+ "a vector with 3 components.");
+
+
+
+} // initialize
+// ----------------------------------------------------------------------
+// Integrate contribution of cohesive cells to residual term.
+void
+pylith::faults::FaultCohesiveKin::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateResidual
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute Jacobian matrix (A) associated with operator.
+void
+pylith::faults::FaultCohesiveKin::integrateJacobian(
+ PetscMat* mat,
+ const ALE::Obj<real_section_type>& dispT,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateJacobian
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Set field.
+void
+pylith::faults::FaultCohesiveKin::setField(
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // setField
+} // setField
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-04-26 22:03:45 UTC (rev 6699)
@@ -59,6 +59,45 @@
*/
void eqsrc(EqKinSrc* src);
+ /** Initialize fault. Determine orientation and setup boundary
+ * condition parameters.
+ *
+ * @param mesh PETSc mesh
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ void initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const double_array& upDir);
+
+ /** 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
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Compute Jacobian matrix (A) associated with operator.
+ *
+ * @param mat Sparse matrix
+ * @param disp Displacement field
+ * @param mesh Finite-element mesh
+ */
+ void integrateJacobian(PetscMat* mat,
+ const ALE::Obj<real_section_type>& dispT,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Set field.
+ *
+ * @param disp Displacement field
+ * @param mesh Finite-element mesh
+ */
+ void setField(const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-04-26 22:03:45 UTC (rev 6699)
@@ -95,7 +95,7 @@
void material(const materials::ElasticMaterial* m);
/** Integrate constant term (b) for dynamic elasticity term
- * for 3-D finite elements.
+ * 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]
Modified: short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2007-04-26 06:20:17 UTC (rev 6698)
+++ short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2007-04-26 22:03:45 UTC (rev 6699)
@@ -87,7 +87,7 @@
try {
assert(0 != objVptr);
ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) meshVptr;
- ((pylith::faults::Fault*) objVptr)->adjustTopology(mesh);
+ ((pylith::faults::Fault*) objVptr)->adjustTopology(*mesh);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
More information about the cig-commits
mailing list