[cig-commits] r16141 - in short/3D/PyLith/trunk: libsrc/faults unittests/libtests/faults unittests/libtests/faults/data
surendra at geodynamics.org
surendra at geodynamics.org
Sun Jan 17 22:48:21 PST 2010
Author: surendra
Date: 2010-01-17 22:48:20 -0800 (Sun, 17 Jan 2010)
New Revision: 16141
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc
Log:
Unittests for frictions (stick & open cases) are updated and working
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc 2010-01-17 23:13:23 UTC (rev 16140)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc 2010-01-18 06:48:20 UTC (rev 16141)
@@ -19,7 +19,7 @@
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
-#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/Mesh.hh" // USES MeshslipTpdtVertex
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc 2010-01-17 23:13:23 UTC (rev 16140)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDynL.cc 2010-01-18 06:48:20 UTC (rev 16141)
@@ -193,9 +193,115 @@
void
pylith::faults::TestFaultCohesiveDynL::testConstrainSolnSpaceStick(void)
{ // testConstrainSolnSpaceStick
- // STUFF GOES HERE (Surendra)
- // No change to dispIncr field
- // Slip field should be zero
+ 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->fieldIncrStick);
+
+ 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
+ // Lagrange multipliers should be adjusted according to friction
+ // (No change to dispIncr field).
+ 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();
+
+ // Get section containing solution (disp + Lagrange multipliers)
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields.get("dispIncr(t->t+dt)").section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+
+ // Get expected values
+ const double* valsE = _data->fieldIncrStick; // No change in dispIncr
+ int iVertex = 0; // variable to use as index into valsE array
+ const int fiberDimE = spaceDim; // number of values per point
+ const double tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex) { // loop over all vertices in mesh
+ // Check fiber dimension (number of values at point)
+ const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const double* vals = dispIncrSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ // Check values at point
+ 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
+ // Slip values should be adjusted based on the change in the
+ // Lagrange multipliers (slip should be zero).
+
+ // Get fault vertex info
+ 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();
+
+ // Get section containing slip
+ const ALE::Obj<RealSection>& slipSection =
+ fault._fields->get("slip").section();
+ CPPUNIT_ASSERT(!slipSection.isNull());
+
+ // Get expected values
+ // const double* valsE = _data->slipStickE;
+ double_array slipStickE(4);
+ slipStickE = 0.0;
+ const double* valsE = &slipStickE[0];
+ int iVertex = 0; // variable to use as index into valsE array
+ const int fiberDimE = spaceDim; // number of values per point
+ const double tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd;
+ ++v_iter, ++iVertex) { // loop over fault vertices
+ // Check fiber dimension (number of values at point)
+ const int fiberDim = slipSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const double* vals = slipSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ // Check values at point
+ 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
+
} // testConstrainSolnSpaceStick
// ----------------------------------------------------------------------
@@ -319,7 +425,112 @@
void
pylith::faults::TestFaultCohesiveDynL::testConstrainSolnSpaceOpen(void)
{ // testConstrainSolnSpaceOpen
- // STUFF GOES HERE (Surendra)
+ 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->fieldIncrOpen);
+
+ 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
+ // Lagrange multipliers should be adjusted according to friction
+ // as reflected in the fieldIncrOpenE data member.
+ 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();
+
+ // Get section containing solution (disp + Lagrange multipliers)
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields.get("dispIncr(t->t+dt)").section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+
+ // Get expected values
+ const double* valsE = _data->fieldIncrOpenE; // Expected values for dispIncr
+ int iVertex = 0; // variable to use as index into valsE array
+ const int fiberDimE = spaceDim; // number of values per point
+ const double tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex) { // loop over all vertices in mesh
+ // Check fiber dimension (number of values at point)
+ const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const double* vals = dispIncrSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ // Check values at point
+ 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
+ // Slip values should be adjusted based on the change in the
+ // Lagrange multipliers as reflected in the slipOpenE data member.
+
+ // Get fault vertex info
+ 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();
+
+ // Get section containing slip
+ const ALE::Obj<RealSection>& slipSection =
+ fault._fields->get("slip").section();
+ CPPUNIT_ASSERT(!slipSection.isNull());
+
+ // Get expected values
+ const double* valsE = _data->slipOpenE;
+ int iVertex = 0; // variable to use as index into valsE array
+ const int fiberDimE = spaceDim; // number of values per point
+ const double tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd;
+ ++v_iter, ++iVertex) { // loop over fault vertices
+ // Check fiber dimension (number of values at point)
+ const int fiberDim = slipSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const double* vals = slipSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ // Check values at point
+ 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
+
} // testConstrainSolnSpaceOpen
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc 2010-01-17 23:13:23 UTC (rev 16140)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynLDataTri3.cc 2010-01-18 06:48:20 UTC (rev 16141)
@@ -265,14 +265,14 @@
// ----------------------------------------------------------------------
// 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
+ 1.1, 29.1,
+ 1.2, 29.2, // 3
+ 1.3, 29.3, // 4
+ 1.4, 29.4,
+ 1.5, 29.5, // 6
+ 1.7, 29.7, // 7
+ 1.6, -29.6, // 8
+ 1.8, -29.8, // 9
};
// No change in fieldIncr
@@ -317,31 +317,31 @@
// ----------------------------------------------------------------------
// 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
+ 9.1, 10.1,
+ 9.2, 10.2, // 3
+ 9.3, 10.3, // 4
+ 9.4, 10.4,
+ 9.5, 10.5, // 6
+ 9.7, 10.7, // 7
+ 9.6, 10.6, // 8
+ 9.8, 10.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
+ 9.1, 10.1,
+ 9.2, 10.2, // 3
+ 9.3, 10.3, // 4
+ 9.4, 10.4,
+ 9.5, 10.5, // 6
+ 9.7, 10.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,
+ 36.4, 40.4,
+ 37.2, 41.2,
};
// ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list