[cig-commits] r16110 - in short/3D/PyLith/trunk: . libsrc/faults unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Fri Dec 18 13:44:28 PST 2009
Author: brad
Date: 2009-12-18 13:44:27 -0800 (Fri, 18 Dec 2009)
New Revision: 16110
Added:
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_initialtract.spatialdb
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc
short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am
Log:
Worked on setting up unit tests for friction using tri3 mesh as initial test case.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/TODO 2009-12-18 21:44:27 UTC (rev 16110)
@@ -2,53 +2,33 @@
CURRENT ISSUES/PRIORITIES
======================================================================
-Status?
+* Uniform global refinement for tets with faults
- Field split
- uniform global refinement for tets with faults
- customized SNES line search
- Laplacian based preconditioner
+* Better preconditioning
+ Better PETSc settings or Laplacian based preconditioner
-Data structure to hold label of fault Lagrange vertices and
+* Data structure to hold label of fault Lagrange vertices and
conventional vertices
array of structure (lagrange, positive, negative, fault)
no mesh queries? Just loops?
-Drucker-Prager elastoplastic (or viscoelastoplastic?)
+* Drucker-Prager elastoplastic (or viscoelastoplastic?)
-Lumped solver
+* Friction
-Setup fault constitutive models
+* Lumped solver
+* Updates to manual
+ + tutorial
+ + governing equations
+* Cleanup
+ + memory model
+ + full-scale testing
+ + Code cleanup
----------------------------------------------------------------------
-TODO WELL BEFORE WORKSHOP 2010
-----------------------------------------------------------------------
-
-Tutorial
- 3d/hex8
- 1. Dirichlet BC (static)
- 2. Neumann BC (static)
- 3. Earthquake rupture (static)
- 4. Dirichlet BC (quasi-static)
- 5. Dirichlet + Neumann (quasi-static)
- 6. Multiple earthquake rupture + creep (quasi-static)
- 7. Earthquake rupture + creep + Dirichlet BC (quasi-static)
- 8. Add power-law rheology (quasi-static)
- 9. Static friction + Dirichlet BC (static)
- 10. Static friction + Dirchlet BC (quasi-static)
- 11. Rate- and state-friction + Dirichlet BC (quasi-static)??
-
-Write up description of Savage and Prescott (1978) benchmark and
-distribute to Greg Lyzenga and Jay Parker
-
-Rewrite governing equations to include fault implementation
-
-Rewrite bulk constitutive models
-
-----------------------------------------------------------------------
FRICTION
----------------------------------------------------------------------
@@ -77,6 +57,7 @@
Rate- and state-friction with aging law
Rate- and state-friction with slip law
+
----------------------------------------------------------------------
LUMPED SOLVER
----------------------------------------------------------------------
@@ -100,14 +81,34 @@
adjust displacements
----------------------------------------------------------------------
->>>>>>> .merge-right.r16093
-POST RELEASE 1.4
+TODO WELL BEFORE WORKSHOP 2010
----------------------------------------------------------------------
-Matt
+Tutorial
+ 3d/hex8
+ 1. Dirichlet BC (static)
+ 2. Neumann BC (static)
+ 3. Earthquake rupture (static)
+ 4. Dirichlet BC (quasi-static)
+ 5. Dirichlet + Neumann (quasi-static)
+ 6. Multiple earthquake rupture + creep (quasi-static)
+ 7. Earthquake rupture + creep + Dirichlet BC (quasi-static)
+ 8. Add power-law rheology (quasi-static)
+ 9. Static friction + Dirichlet BC (static)
+ 10. Static friction + Dirchlet BC (quasi-static)
+ 11. Rate- and state-friction + Dirichlet BC (quasi-static)??
- field split
+Write up description of Savage and Prescott (1978) benchmark and
+distribute to Greg Lyzenga and Jay Parker
+Rewrite governing equations to include fault implementation
+
+Rewrite bulk constitutive models
+
+----------------------------------------------------------------------
+POST RELEASE 1.4
+----------------------------------------------------------------------
+
Brad
Memory model
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc 2009-12-18 21:44:27 UTC (rev 16110)
@@ -49,7 +49,8 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesiveDynL::FaultCohesiveDynL(void)
+pylith::faults::FaultCohesiveDynL::FaultCohesiveDynL(void) :
+ _dbInitialTract(0)
{ // constructor
} // constructor
@@ -68,6 +69,7 @@
FaultCohesive::deallocate();
// :TODO: Use shared pointers for initial database
+ _dbInitialTract = 0;
} // deallocate
// ----------------------------------------------------------------------
@@ -683,7 +685,8 @@
slipSection->view("SLIP");
areaSection->view("AREA");
- tractionInitialSection->view("INITIAL TRACTION");
+ if (0 != _dbInitialTract)
+ tractionInitialSection->view("INITIAL TRACTION");
dispTSection->view("DISP (t)");
dispTIncrSection->view("DISP INCR (t->t+dt)");
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2009-12-18 21:44:27 UTC (rev 16110)
@@ -46,6 +46,7 @@
TestFaultCohesiveKinSrcsTet4.cc \
TestFaultCohesiveKinSrcsHex8.cc \
TestFaultCohesiveDynL.cc \
+ TestFaultCohesiveDynLTri3.cc \
test_faults.cc
@@ -73,7 +74,8 @@
TestFaultCohesiveKinSrcsQuad4.hh \
TestFaultCohesiveKinSrcsTet4.hh \
TestFaultCohesiveKinSrcsHex8.hh \
- TestFaultCohesiveDynL.hh
+ TestFaultCohesiveDynL.hh \
+ TestFaultCohesiveDynLTri3.hh
# Source files associated with testing data
testfaults_SOURCES += \
@@ -131,7 +133,8 @@
data/CohesiveKinSrcsDataQuad4.cc \
data/CohesiveKinSrcsDataTet4.cc \
data/CohesiveKinSrcsDataHex8.cc \
- data/CohesiveDynLData.cc
+ data/CohesiveDynLData.cc \
+ data/CohesiveDynLDataTri3.cc
noinst_HEADERS += \
data/CohesiveData.hh \
@@ -188,7 +191,8 @@
data/CohesiveKinSrcsDataQuad4.hh \
data/CohesiveKinSrcsDataTet4.hh \
data/CohesiveKinSrcsDataHex8.hh \
- data/CohesiveDynLData.hh
+ data/CohesiveDynLData.hh \
+ data/CohesiveDynLDataTri3.hh
AM_CPPFLAGS = \
$(PYTHON_EGG_CPPFLAGS) -I$(PYTHON_INCDIR) \
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc 2009-12-18 21:44:27 UTC (rev 16110)
@@ -97,7 +97,7 @@
topology::Mesh mesh;
FaultCohesiveDynL fault;
topology::SolutionFields fields(mesh);
- _initialize(&mesh, &fault, &fields, _data->fieldIncrStick);
+ _initialize(&mesh, &fault, &fields);
const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
@@ -203,7 +203,97 @@
void
pylith::faults::TestFaultCohesiveDynL::testConstrainSolnSpaceSlip(void)
{ // testConstrainSolnSpaceSlip
- // STUFF GOES HERE (Brad)
+ assert(0 != _data);
+
+ topology::Mesh mesh;
+ FaultCohesiveDynL fault;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
+ topology::Jacobian jacobian(fields, "seqdense");
+ _setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrSlip);
+
+ const int spaceDim = _data->spaceDim;
+
+ const double t = 2.134;
+ const double dt = 0.01;
+ fault.timeStep(dt);
+ fault.constrainSolnSpace(&fields, t, jacobian);
+
+ //residual.view("RESIDUAL"); // DEBUGGING
+
+ { // Check solution values
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields.get("dispIncr(t->t+dt)").section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+
+ const double* valsE = _data->fieldIncrOpenE;
+ int iVertex = 0;
+ const int fiberDimE = spaceDim;
+ const double tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex) {
+ const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const double* vals = dispIncrSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ for (int i = 0; i < fiberDimE; ++i) {
+ const int index = iVertex * spaceDim + i;
+ const double valE = valsE[index];
+ if (fabs(valE) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ } // for
+ } // for
+ } // Check solution values
+
+ { // Check slip values
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ SieveSubMesh::renumbering_type& renumbering =
+ faultSieveMesh->getRenumbering();
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ const ALE::Obj<RealSection>& slipSection =
+ fault._fields->get("slip").section();
+ CPPUNIT_ASSERT(!slipSection.isNull());
+
+ const double* valsE = _data->slipSlipE;
+ int iVertex = 0;
+ const int fiberDimE = spaceDim;
+ const double tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter, ++iVertex) {
+ const int fiberDim = slipSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const double* vals = slipSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ for (int i = 0; i < fiberDimE; ++i) {
+ const int index = iVertex * spaceDim + i;
+ const double valE = valsE[index];
+ if (fabs(valE) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ } // for
+ } // for
+ } // Check slip values
+
} // testConstrainSolnSpaceSlip
// ----------------------------------------------------------------------
@@ -232,7 +322,9 @@
topology::Mesh mesh;
FaultCohesiveDynL fault;
topology::SolutionFields fields(mesh);
- _initialize(&mesh, &fault, &fields, _data->fieldIncrStick);
+ _initialize(&mesh, &fault, &fields);
+ topology::Jacobian jacobian(fields, "seqdense");
+ _setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrStick);
// Test with dbInitialTract
// :TODO: STUFF GOES HERE (Brad)
@@ -247,8 +339,7 @@
pylith::faults::TestFaultCohesiveDynL::_initialize(
topology::Mesh* const mesh,
FaultCohesiveDynL* const fault,
- topology::SolutionFields* const fields,
- const double* const fieldIncr)
+ topology::SolutionFields* const fields)
{ // _initialize
CPPUNIT_ASSERT(0 != mesh);
CPPUNIT_ASSERT(0 != fault);
@@ -302,19 +393,98 @@
fault->initialize(*mesh, upDir, normalDir);
// Setup fields
- fields->add("residual", "residual");
fields->add("disp(t)", "displacement");
fields->add("dispIncr(t->t+dt)", "displacement_increment");
fields->solutionName("dispIncr(t->t+dt)");
const int spaceDim = _data->spaceDim;
- topology::Field<topology::Mesh>& residual = fields->get("residual");
- residual.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
- residual.allocate();
- fields->copyLayout("residual");
+ topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
+ disp.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
+ disp.allocate();
+ fields->copyLayout("disp(t)");
} // _initialize
// ----------------------------------------------------------------------
+// Set values for fields and Jacobian.
+void
+pylith::faults::TestFaultCohesiveDynL::_setFieldsJacobian(
+ topology::Mesh* const mesh,
+ FaultCohesiveDynL* const fault,
+ topology::SolutionFields* const fields,
+ topology::Jacobian* const jacobian,
+ const double* const fieldIncr)
+{ // _initialize
+ CPPUNIT_ASSERT(0 != mesh);
+ CPPUNIT_ASSERT(0 != fault);
+ CPPUNIT_ASSERT(0 != fields);
+ CPPUNIT_ASSERT(0 != jacobian);
+ CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(0 != fieldIncr);
+
+ const int spaceDim = _data->spaceDim;
+
+ // Get vertices in mesh
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Set displacement values
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
+ int iVertex = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex)
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+
+ // Set increment values
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ iVertex = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex)
+ dispIncrSection->updatePoint(*v_iter, &fieldIncr[iVertex*spaceDim]);
+
+ // Setup Jacobian matrix
+ const int nrows = dispIncrSection->sizeWithBC();
+ const int ncols = nrows;
+ int nrowsM = 0;
+ int ncolsM = 0;
+ PetscMat jacobianMat = jacobian->matrix();
+ MatGetSize(jacobianMat, &nrowsM, &ncolsM);
+ CPPUNIT_ASSERT_EQUAL(nrows, nrowsM);
+ CPPUNIT_ASSERT_EQUAL(ncols, ncolsM);
+
+ int_array rows(nrows);
+ int_array cols(ncols);
+ for (int iRow=0; iRow < nrows; ++iRow)
+ rows[iRow] = iRow;
+ for (int iCol=0; iCol < ncols; ++iCol)
+ cols[iCol] = iCol;
+ MatSetValues(jacobianMat, nrows, &rows[0], ncols, &cols[0], _data->jacobian, INSERT_VALUES);
+ jacobian->assemble("final_assembly");
+
+ // Set Jacobian diagonal
+ fields->add("Jacobian diagonal", "jacobian_diagonal");
+ topology::Field<topology::Mesh>& jacobianDiag =
+ fields->get("Jacobian diagonal");
+ const topology::Field<topology::Mesh>& disp =
+ fields->get("disp(t)");
+ jacobianDiag.cloneSection(disp);
+ jacobianDiag.createVector();
+ jacobianDiag.createScatter();
+ MatGetDiagonal(jacobian->matrix(), jacobianDiag.vector());
+ jacobianDiag.scatterVectorToSection();
+} // _setFieldsJacobian
+
+// ----------------------------------------------------------------------
// Determine if vertex is a Lagrange multiplier constraint vertex.
bool
pylith::faults::TestFaultCohesiveDynL::_isConstraintVertex(const int vertex) const
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.hh 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.hh 2009-12-18 21:44:27 UTC (rev 16110)
@@ -27,9 +27,7 @@
#include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
#include "pylith/feassemble/feassemblefwd.hh" // HOLDSA Quadrature
#include "spatialdata/spatialdb/spatialdbfwd.hh" // HOLDSA SpatialDB
-
#include <vector> // HASA std::vector
-
/// Namespace for pylith package
namespace pylith {
namespace faults {
@@ -40,27 +38,26 @@
} // pylith
/// C++ unit testing for FaultCohesiveDynL
-class pylith::faults::TestFaultCohesiveDynL : public CppUnit::TestFixture
-{ // class TestFaultCohesiveDynL
+class pylith::faults::TestFaultCohesiveDynL: public CppUnit::TestFixture { // class TestFaultCohesiveDynL
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUITE( TestFaultCohesiveDynL );
+CPPUNIT_TEST_SUITE( TestFaultCohesiveDynL );
- CPPUNIT_TEST( testConstructor );
- CPPUNIT_TEST( testDBInitialTract );
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testDBInitialTract );
- // Tests in derived classes:
- // testInitialize()
- // testConstrainSolnSpaceStick()
- // testConstrainSolnSpaceSlip()
- // testConstrainSolnSpaceOpen()
- // testUpdateStateVars()
- // testCalcTractions()
+ // Tests in derived classes:
+ // testInitialize()
+ // testConstrainSolnSpaceStick()
+ // testConstrainSolnSpaceSlip()
+ // testConstrainSolnSpaceOpen()
+ // testUpdateStateVars()
+ // testCalcTractions()
CPPUNIT_TEST_SUITE_END();
// PROTECTED MEMBERS //////////////////////////////////////////////////
-protected :
+protected:
CohesiveDynLData* _data; ///< Data for testing
feassemble::Quadrature<topology::SubMesh>* _quadrature; ///< Fault quad.
@@ -68,7 +65,7 @@
bool _flipFault; ///< If true, flip fault orientation.
// PUBLIC METHODS /////////////////////////////////////////////////////
-public :
+public:
/// Setup testing data.
void setUp(void);
@@ -101,20 +98,36 @@
void testCalcTractions(void);
// PRIVATE METHODS ////////////////////////////////////////////////////
-private :
+private:
/** Initialize FaultCohesiveDynL interface condition.
*
* @param mesh PETSc mesh to initialize
* @param fault Cohesive fault interface condition to initialize.
* @param fields Solution fields.
- * @param fieldIncrVals Values for solution increment field.
*/
void _initialize(topology::Mesh* const mesh,
- FaultCohesiveDynL* const fault,
- topology::SolutionFields* const fields,
- const double* const fieldIncrVals);
+ FaultCohesiveDynL* const fault,
+ topology::SolutionFields* const fields);
+ /** Set values for fields and Jacobian.
+ *
+ * @pre Must call _initialize() before _setFieldsJacobian to set solution
+ * field. Note: Call to Jacobian constructor should be after call to
+ * _initialize() so that the solution field is setup.
+ *
+ * @param mesh PETSc mesh to initialize
+ * @param fault Cohesive fault interface condition to initialize.
+ * @param fields Solution fields.
+ * @param jacobian Jacobian sparse matrix.
+ * @param fieldIncrVals Values for solution increment field.
+ */
+ void _setFieldsJacobian(topology::Mesh* const mesh,
+ FaultCohesiveDynL* const fault,
+ topology::SolutionFields* const fields,
+ topology::Jacobian* const jacobian,
+ const double* const fieldIncrVals);
+
/** Determine if vertex is a Lagrange multiplier constraint vertex.
*
* @param vertex Label of vertex.
@@ -122,11 +135,9 @@
* @returns True if vertex is a constraint vertex, false otherwise.
*/
bool _isConstraintVertex(const int vertex) const;
-
}; // class TestFaultCohesiveDynL
#endif // pylith_faults_testfaultcohesivedynl_hh
-
// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.cc 2009-12-18 21:44:27 UTC (rev 16110)
@@ -0,0 +1,42 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestFaultCohesiveDynLTri3.hh" // Implementation of class methods
+
+#include "data/CohesiveDynLDataTri3.hh" // USES CohesiveDynLDataTri3
+
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature<SubMesh>
+#include "pylith/feassemble/GeometryLine2D.hh" // USES GeometryLine2D
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesiveDynLTri3 );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::faults::TestFaultCohesiveDynLTri3::setUp(void)
+{ // setUp
+ TestFaultCohesiveDynL::setUp();
+ _data = new CohesiveDynLDataTri3();
+
+ CPPUNIT_ASSERT(0 != _quadrature);
+ feassemble::GeometryLine2D geometry;
+ _quadrature->refGeometry(&geometry);
+
+ _flipFault = true;
+} // setUp
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynLTri3.hh 2009-12-18 21:44:27 UTC (rev 16110)
@@ -0,0 +1,60 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/faults/TestFaultCohesiveDynLTri3.hh
+ *
+ * @brief C++ TestFaultCohesiveDynLTri3 object.
+ *
+ * C++ unit testing for FaultCohesiveDynL for mesh with 2-D triangular cells.
+ */
+
+#if !defined(pylith_faults_testfaultcohesivedynltri3_hh)
+#define pylith_faults_testfaultcohesivedynltri3_hh
+
+#include "TestFaultCohesiveDynL.hh" // ISA TestFaultCohesiveDynL
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace faults {
+ class TestFaultCohesiveDynLTri3;
+ } // bc
+} // pylith
+
+/// C++ unit testing for FaultCohesiveDynL for mesh with 2-D triangular cells.
+class pylith::faults::TestFaultCohesiveDynLTri3 : public TestFaultCohesiveDynL
+{ // class TestFaultCohesiveDynLTri3
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestFaultCohesiveDynLTri3 );
+
+ CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testConstrainSolnSpaceStick );
+ CPPUNIT_TEST( testConstrainSolnSpaceSlip );
+ CPPUNIT_TEST( testConstrainSolnSpaceOpen );
+ CPPUNIT_TEST( testUpdateStateVars );
+ CPPUNIT_TEST( testCalcTractions );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+}; // class TestFaultCohesiveDynLTri3
+
+#endif // pylith_faults_testfaultcohesivedynltri3_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.cc 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.cc 2009-12-18 21:44:27 UTC (rev 16110)
@@ -37,9 +37,9 @@
initialTractions(0),
area(0),
fieldIncrSlipE(0),
- slipSlip(0),
+ slipSlipE(0),
fieldIncrOpenE(0),
- slipOpen(0),
+ slipOpenE(0),
constraintVertices(0),
constraintCells(0),
numConstraintVert(0)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.hh 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLData.hh 2009-12-18 21:44:27 UTC (rev 16110)
@@ -53,7 +53,7 @@
//@{
int id; ///< Fault material identifier
char* label; ///< Label for fault
- char* initialTractFilename; ///< Name of db for final slip
+ char* initialTractFilename; ///< Name of db for initial tractions.
//@}
/// @name Input fields
@@ -71,9 +71,9 @@
double* area; ///< Expected values for fault area.
double* initialTractions; ///< Expected values for initial tractions.
double* fieldIncrSlipE; ///< Expected values for solution increment for slipping case.
- double* slipSlip; ///< Expected values for slip for slipping case.
+ double* slipSlipE; ///< Expected values for slip for slipping case.
double* fieldIncrOpenE; ///< Expected values for solution increment for opening case.
- double* slipOpen; ///< Expected values for slip for opening case.
+ double* slipOpenE; ///< Expected values for slip for opening case.
int* constraintVertices; ///< Expected points for constraint vertices
int* constraintCells; ///< Expected cells for constraint vertices
Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc 2009-12-18 21:44:27 UTC (rev 16110)
@@ -0,0 +1,387 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/* Original mesh
+ *
+ * Cells are 0-1, vertices are 2-5.
+ *
+ * 3
+ * /|\
+ * / | \
+ * / | \
+ * / | \
+ * 2 | 5
+ * \ | /
+ * \ | /
+ * \ | /
+ * \|/
+ * 4
+ *
+ *
+ * After adding cohesive elements
+ *
+ * Cells are 0-1, 8, vertices are 2-7.
+ *
+ * 6 -8- 3
+ * /| |\
+ * / | | \
+ * / | | \
+ * / | | \
+ * 2 | | 5
+ * \ | | /
+ * \ | | /
+ * \ | | /
+ * \| |/
+ * 7 -9- 4
+ */
+
+#include "CohesiveDynLDataTri3.hh"
+
+const char* pylith::faults::CohesiveDynLDataTri3::_meshFilename =
+ "data/tri3.mesh";
+
+const int pylith::faults::CohesiveDynLDataTri3::_spaceDim = 2;
+
+const int pylith::faults::CohesiveDynLDataTri3::_cellDim = 1;
+
+const int pylith::faults::CohesiveDynLDataTri3::_numBasis = 2;
+
+const int pylith::faults::CohesiveDynLDataTri3::_numQuadPts = 1;
+
+const double pylith::faults::CohesiveDynLDataTri3::_quadPts[] = {
+ 0.0,
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_quadWts[] = {
+ 2.0,
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_basis[] = {
+ 0.5,
+ 0.5
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_basisDeriv[] = {
+ -0.5,
+ 0.5
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_verticesRef[] = {
+ -1.0, 1.0
+};
+
+const int pylith::faults::CohesiveDynLDataTri3::_id = 10;
+
+const char* pylith::faults::CohesiveDynLDataTri3::_label = "fault";
+
+const char* pylith::faults::CohesiveDynLDataTri3::_initialTractFilename =
+ "data/tri3_initialtract.spatialdb";
+
+const double pylith::faults::CohesiveDynLDataTri3::_fieldT[] = {
+ 8.1, 9.1,
+ 8.2, 9.2, // 3
+ 8.3, 9.3, // 4
+ 8.4, 9.4,
+ 8.5, 9.5, // 6
+ 8.7, 9.7, // 7
+ 8.6, 9.6, // 8
+ 8.8, 9.8, // 9
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_jacobian[] = {
+ 1.0, 0.0, // 2x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 1.0, // 2y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 3x
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0,+1.0, // 8
+ 0.0, 0.0,
+ 0.0, 0.0, // 3y
+ 0.0, 1.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ +1.0, 0.0, // 8
+ 0.0, 0.0,
+ 0.0, 0.0, // 4x
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0,+1.0, // 9
+ 0.0, 0.0, // 4y
+ 0.0, 0.0,
+ 0.0, 1.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ +1.0, 0.0, // 9
+ 0.0, 0.0, // 5x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 5y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 1.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 6x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0,-1.0, // 8
+ 0.0, 0.0,
+ 0.0, 0.0, // 6y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 1.0,
+ 0.0, 0.0,
+ -1.0, 0.0, // 8
+ 0.0, 0.0,
+ 0.0, 0.0, // 7x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0,-1.0, // 9
+ 0.0, 0.0, // 7y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 1.0,
+ 0.0, 0.0,
+ -1.0, 0.0, // 9
+
+ 0.0, 0.0, // 8x
+ 0.0,+1.0, // 3
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0,-1.0, // 6
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 8y
+ +1.0, 0.0, // 3
+ 0.0, 0.0,
+ 0.0, 0.0,
+ -1.0, 0.0, // 6
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+
+ 0.0, 0.0, // 9x
+ 0.0, 0.0,
+ 0.0,+1.0, // 4
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0,-1.0, // 7
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 9y
+ 0.0, 0.0,
+ +1.0, 0.0, // 4
+ 0.0, 0.0,
+ 0.0, 0.0,
+ -1.0, 0.0, // 7
+ 0.0, 0.0,
+ 0.0, 0.0,
+};
+
+// ----------------------------------------------------------------------
+// Computed values
+// ----------------------------------------------------------------------
+
+const double pylith::faults::CohesiveDynLDataTri3::_orientation[] = {
+ 0.0, -1.0, -1.0, 0.0,
+ 0.0, -1.0, -1.0, 0.0
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_area[] = {
+ 1.0,
+ 1.0,
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_initialTractions[] = {
+ 1.0, -2.0,
+ 1.1, -2.1,
+};
+
+
+const int pylith::faults::CohesiveDynLDataTri3::_numConstraintVert = 2;
+const int pylith::faults::CohesiveDynLDataTri3::_constraintVertices[] = {
+ 8, 9
+};
+const int pylith::faults::CohesiveDynLDataTri3::_constraintCells[] = {
+ 11, 11
+};
+
+// ----------------------------------------------------------------------
+// Stick case
+// ----------------------------------------------------------------------
+// Input
+const double pylith::faults::CohesiveDynLDataTri3::_fieldIncrStick[] = {
+ 8.1, 9.1,
+ 8.2, 9.2, // 3
+ 8.3, 9.3, // 4
+ 8.4, 9.4,
+ 8.5, 9.5, // 6
+ 8.7, 9.7, // 7
+ 8.6, 9.6, // 8
+ 8.8, 9.8, // 9
+};
+
+// No change in fieldIncr
+// Zero slip
+
+// ----------------------------------------------------------------------
+// Slip case
+// ----------------------------------------------------------------------
+// Input
+const double pylith::faults::CohesiveDynLDataTri3::_fieldIncrSlip[] = {
+ 8.1, 9.1,
+ 8.2, 9.2, // 3
+ 8.3, 9.3, // 4
+ 8.4, 9.4,
+ 8.5, 9.5, // 6
+ 8.7, 9.7, // 7
+ 8.6, 9.6, // 8
+ 8.8, 9.8, // 9
+};
+
+// Output
+const double pylith::faults::CohesiveDynLDataTri3::_fieldIncrSlipE[] = {
+ 8.1, 9.1,
+ 8.2, 9.2, // 3
+ 8.3, 9.3, // 4
+ 8.4, 9.4,
+ 8.5, 9.5, // 6
+ 8.7, 9.7, // 7
+ 8.6, 9.6, // 8
+ 8.8, 9.8, // 9
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_slipSlipE[] = {
+ 0.0, 0.0,
+ 0.0, 0.0,
+};
+
+// ----------------------------------------------------------------------
+// Open case
+// ----------------------------------------------------------------------
+// Input
+const double pylith::faults::CohesiveDynLDataTri3::_fieldIncrOpen[] = {
+ 8.1, 9.1,
+ 8.2, 9.2, // 3
+ 8.3, 9.3, // 4
+ 8.4, 9.4,
+ 8.5, 9.5, // 6
+ 8.7, 9.7, // 7
+ 8.6, 9.6, // 8
+ 8.8, 9.8, // 9
+};
+
+// Output
+const double pylith::faults::CohesiveDynLDataTri3::_fieldIncrOpenE[] = {
+ 8.1, 9.1,
+ 8.2, 9.2, // 3
+ 8.3, 9.3, // 4
+ 8.4, 9.4,
+ 8.5, 9.5, // 6
+ 8.7, 9.7, // 7
+ 8.6, 9.6, // 8
+ 8.8, 9.8, // 9
+};
+
+const double pylith::faults::CohesiveDynLDataTri3::_slipOpenE[] = {
+ 0.0, 0.0,
+ 0.0, 0.0,
+};
+
+// ----------------------------------------------------------------------
+pylith::faults::CohesiveDynLDataTri3::CohesiveDynLDataTri3(void)
+{ // constructor
+ meshFilename = const_cast<char*>(_meshFilename);
+ spaceDim = _spaceDim;
+ cellDim = _cellDim;
+ numBasis = _numBasis;
+ numQuadPts = _numQuadPts;
+ quadPts = const_cast<double*>(_quadPts);
+ quadWts = const_cast<double*>(_quadWts);
+ basis = const_cast<double*>(_basis);
+ basisDeriv = const_cast<double*>(_basisDeriv);
+ verticesRef = const_cast<double*>(_verticesRef);
+ id = _id;
+ label = const_cast<char*>(_label);
+ initialTractFilename = const_cast<char*>(_initialTractFilename);
+
+ fieldT = const_cast<double*>(_fieldT);
+ jacobian = const_cast<double*>(_jacobian);
+ orientation = const_cast<double*>(_orientation);
+ area = const_cast<double*>(_area);
+
+ constraintVertices = const_cast<int*>(_constraintVertices);
+ constraintCells = const_cast<int*>(_constraintCells);
+ numConstraintVert = _numConstraintVert;
+
+ // Stick
+ fieldIncrStick = const_cast<double*>(_fieldIncrStick);
+
+ // Slip
+ fieldIncrSlip = const_cast<double*>(_fieldIncrSlip);
+ fieldIncrSlipE = const_cast<double*>(_fieldIncrSlipE);
+ slipSlipE = const_cast<double*>(_slipSlipE);
+
+ // Open
+ fieldIncrOpen = const_cast<double*>(_fieldIncrOpen);
+ fieldIncrOpenE = const_cast<double*>(_fieldIncrOpenE);
+ slipOpenE = const_cast<double*>(_slipOpenE);
+} // constructor
+
+pylith::faults::CohesiveDynLDataTri3::~CohesiveDynLDataTri3(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.hh 2009-12-18 21:44:27 UTC (rev 16110)
@@ -0,0 +1,78 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_faults_cohesivedynldatatri3_hh)
+#define pylith_faults_cohesivedynldatatri3_hh
+
+#include "CohesiveDynLData.hh"
+
+namespace pylith {
+ namespace faults {
+ class CohesiveDynLDataTri3;
+ } // pylith
+} // faults
+
+class pylith::faults::CohesiveDynLDataTri3 : public CohesiveDynLData
+{
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public:
+
+ /// Constructor
+ CohesiveDynLDataTri3(void);
+
+ /// Destructor
+ ~CohesiveDynLDataTri3(void);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private:
+
+ static const char* _meshFilename; ///< Filename of input mesh
+
+ static const int _spaceDim; ///< Number of dimensions in vertex coordinates
+ static const int _cellDim; ///< Number of dimensions associated with cell
+
+ static const int _numBasis; ///< Number of vertices in cell
+ static const int _numQuadPts; ///< Number of quadrature points
+ static const double _quadPts[]; ///< Coordinates of quad pts in ref cell
+ static const double _quadWts[]; ///< Weights of quadrature points
+ static const double _basis[]; ///< Basis fns at quadrature points
+ static const double _basisDeriv[]; ///< Derivatives of basis fns at quad pts
+ static const double _verticesRef[]; ///< Coordinates of vertices in ref cell (dual basis)
+
+ static const int _id; ///< Fault material identifier
+ static const char* _label; ///< Label for fault
+ static const char* _initialTractFilename; ///< Name of db for initial tractions.
+
+ static const double _fieldT[]; ///< Solution field at time t.
+ static const double _fieldIncrStick[]; ///< Solution increment at time t for stick case.
+ static const double _fieldIncrSlip[]; ///< Solution increment at time t for slip case.
+ static const double _fieldIncrOpen[]; ///< Solution increment at time t for opening case.
+ static const double _jacobian[]; ///< Jacobian sparse matrix.
+
+ static const double _orientation[]; ///< Expected values for fault orientation.
+ static const double _area[]; ///< Expected values for fault area.
+ static const double _initialTractions[]; ///< Expected values for initial tractions.
+ static const double _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
+ static const double _slipSlipE[]; ///< Expected values for slip for slip case.
+ static const double _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
+ static const double _slipOpenE[]; ///< Expected values for slip for opening case.
+ static const int _constraintVertices[]; ///< Expected points for constraint vertices
+ static const int _constraintCells[]; ///< Expected cells for constraint vertices
+ static const int _numConstraintVert; ///< Number of constraint vertices
+
+};
+
+#endif // pylith_faults_cohesivedynldatatri3_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am 2009-12-18 17:47:14 UTC (rev 16109)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am 2009-12-18 21:44:27 UTC (rev 16110)
@@ -34,6 +34,7 @@
tri3d_finalslip.spatialdb \
tri3d_sliptime.spatialdb \
tri3d_risetime.spatialdb \
+ tri3_initialtract.spatialdb \
quad4.mesh \
quad4b.mesh \
quad4c.mesh \
Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_initialtract.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_initialtract.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_initialtract.spatialdb 2009-12-18 21:44:27 UTC (rev 16110)
@@ -0,0 +1,15 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 2
+ value-names = traction-shear traction-normal
+ value-units = MPa MPa
+ num-locs = 2
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+0.0 +1.0 1.0 -2.0
+0.0 -1.0 1.1 -2.1
More information about the CIG-COMMITS
mailing list