[cig-commits] r14902 - in short/3D/PyLith/trunk: . libsrc/faults pylith/tests tests/1d/line2
brad at geodynamics.org
brad at geodynamics.org
Wed May 6 13:37:51 PDT 2009
Author: brad
Date: 2009-05-06 13:37:51 -0700 (Wed, 06 May 2009)
New Revision: 14902
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/pylith/tests/Fault.py
short/3D/PyLith/trunk/pylith/tests/StateVariables.py
short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py
Log:
Finished 1-D tests.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-05-06 20:35:55 UTC (rev 14901)
+++ short/3D/PyLith/trunk/TODO 2009-05-06 20:37:51 UTC (rev 14902)
@@ -10,7 +10,7 @@
1-D
1. axial disp [DONE]
- 2. dislocation (need fault tests)
+ 2. dislocation [DONE]
2-D
1. axial/shear (DirichletBC)
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-06 20:35:55 UTC (rev 14901)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-06 20:37:51 UTC (rev 14902)
@@ -93,6 +93,7 @@
assert(0 != upDir);
assert(0 != normalDir);
assert(0 != _quadrature);
+ assert(0 != _normalizer);
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
@@ -123,12 +124,16 @@
topology::Field<topology::SubMesh>& slip = _fields->get("slip");
slip.newSection(vertices, cs->spaceDim());
slip.allocate();
+ slip.vectorFieldType(topology::FieldBase::VECTOR);
+ slip.scale(_normalizer->lengthScale());
// Allocate cumulative slip field
_fields->add("cumulative slip", "cumulative_slip");
topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
cumSlip.newSection(slip);
cumSlip.allocate();
+ cumSlip.vectorFieldType(topology::FieldBase::VECTOR);
+ cumSlip.scale(_normalizer->lengthScale());
// Setup pseudo-stiffness of cohesive cells to improve conditioning
// of Jacobian matrix
@@ -149,6 +154,9 @@
// Create empty tractions field for change in fault tractions.
_fields->add("tractions", "tractions_change");
+ topology::Field<topology::SubMesh>& tractions = _fields->get("tractions");
+ tractions.vectorFieldType(topology::FieldBase::VECTOR);
+ tractions.scale(_normalizer->pressureScale());
} // initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/pylith/tests/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/tests/Fault.py 2009-05-06 20:35:55 UTC (rev 14901)
+++ short/3D/PyLith/trunk/pylith/tests/Fault.py 2009-05-06 20:37:51 UTC (rev 14902)
@@ -16,7 +16,7 @@
import numpy
-def check_info(testcase, filename, mesh, fieldNames):
+def check_vertex_fields(testcase, filename, mesh, fieldNames):
"""
Check properties.
"""
@@ -36,6 +36,51 @@
tolerance = 1.0e-5
for name in fieldNames:
+ valuesE = testcase.calcFaultField(name, data['vertices'])
+ values = data['vertex_fields'][name]
+
+ (nverticesE, dimE) = valuesE.shape
+ if 1 == dimE:
+ values = values.reshape( (nvertices, dimE) )
+ (nvertices, dim) = values.shape
+
+ testcase.assertEqual(nverticesE, nvertices)
+ testcase.assertEqual(dimE, dim)
+
+ for i in xrange(dim):
+ ratio = numpy.abs(1.0 - values[:,i]/valuesE[:,i])
+ diff = numpy.abs(values[:,i] - valuesE[:,i])
+ mask = valuesE[:,i] != 0.0
+ okay = mask*(ratio < tolerance) + ~mask*(diff < tolerance)
+ if numpy.sum(okay) != nvertices:
+ print "Error in component %d of field '%s'." % (i, name)
+ print "Expected values:",valuesE
+ print "Output values:",values
+ testcase.assertEqual(numpy.sum(okay), nvertices)
+
+ return
+
+
+def check_data(testcase, filename, mesh, fieldNames):
+ """
+ Check properties.
+ """
+ data = testcase.reader.read(filename)
+
+ # Check cells
+ (ncells, ncorners) = data['cells'].shape
+ testcase.assertEqual(mesh['ncells'], ncells)
+ testcase.assertEqual(mesh['ncorners'], ncorners)
+
+ # Check vertices
+ (nvertices, spaceDim) = data['vertices'].shape
+ testcase.assertEqual(mesh['nvertices'], nvertices)
+ testcase.assertEqual(mesh['spaceDim'], spaceDim)
+
+ # Check fault information
+ tolerance = 1.0e-5
+
+ for name in fieldNames:
valuesE = testcase.calcFaultInfo(name, data['vertices'])
values = data['vertex_fields'][name]
Modified: short/3D/PyLith/trunk/pylith/tests/StateVariables.py
===================================================================
--- short/3D/PyLith/trunk/pylith/tests/StateVariables.py 2009-05-06 20:35:55 UTC (rev 14901)
+++ short/3D/PyLith/trunk/pylith/tests/StateVariables.py 2009-05-06 20:37:51 UTC (rev 14902)
@@ -39,9 +39,13 @@
valuesE = testcase.calcStateVar(name, data['vertices'], data['cells'])
values = data['cell_fields'][name]
- (ncellsE, dim) = valuesE.shape
- values = values.reshape( (ncells, dim) )
+ (ncellsE, dimE) = valuesE.shape
+ if 1 == dimE:
+ values = values.reshape( (ncells, dimE) )
+ (ncells, dim) = values.shape
+
testcase.assertEqual(ncellsE, ncells)
+ testcase.assertEqual(dimE, dim)
for i in xrange(dim):
ratio = numpy.abs(1.0 - values[:,i]/valuesE[:,i])
Modified: short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py 2009-05-06 20:35:55 UTC (rev 14901)
+++ short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py 2009-05-06 20:37:51 UTC (rev 14902)
@@ -18,8 +18,8 @@
from TestLine2 import TestLine2
from pylith.utils.VTKDataReader import has_vtk
from pylith.utils.VTKDataReader import VTKDataReader
+from pylith.tests.Fault import check_vertex_fields
-
# Local version of PyLithApp
from pylith.apps.PyLithApp import PyLithApp
class DislocationApp(PyLithApp):
@@ -75,13 +75,26 @@
return
filename = "%s-fault_info.vtk" % self.outputRoot
- from pylith.tests.Fault import check_info
fields = ["normal_dir", "final_slip", "slip_time"]
- check_info(self, filename, self.faultMesh, fields)
+ check_vertex_fields(self, filename, self.faultMesh, fields)
return
+ def test_fault_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_t0000000.vtk" % self.outputRoot
+ fields = ["cumulative_slip", "traction_change"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
def calcDisplacements(self, vertices):
"""
Calculate displacement field given coordinates of vertices.
@@ -123,7 +136,7 @@
return stateVar
- def calcFaultInfo(self, name, vertices):
+ def calcFaultField(self, name, vertices):
"""
Calculate fault info.
"""
@@ -132,10 +145,16 @@
finalSlip = -0.5
slipTime = 0.0
+ slip = 1.0
+ exx = -0.025
+ lp2m = self.density*self.vp**2
+ traction = -exx*lp2m
+
nvertices = self.faultMesh['nvertices']
if name == "normal_dir":
- field = normalDir*numpy.ones( (nvertices, 1), dtype=numpy.float64)
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir
elif name == "final_slip":
field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
@@ -144,8 +163,16 @@
elif name == "slip_time":
field = slipTime*numpy.ones( (nvertices, 1), dtype=numpy.float64)
+ elif name == "cumulative_slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = finalSlip
+
+ elif name == "traction_change":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = traction
+
else:
- raise ValueError("Unknown fault info field '%s'." % name)
+ raise ValueError("Unknown fault field '%s'." % name)
return field
More information about the CIG-COMMITS
mailing list