[cig-commits] r14945 - in short/3D/PyLith/trunk: libsrc/faults pylith/problems unittests/libtests/faults unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Sat May 9 13:49:01 PDT 2009
Author: brad
Date: 2009-05-09 13:49:00 -0700 (Sat, 09 May 2009)
New Revision: 14945
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
Log:
Fixed bug in not updating fault geometry when integrating residual (confinged to SWIG branch?). Fixed bugs associated with switching explicit time stepping to incremental form.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-09 19:29:28 UTC (rev 14944)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-09 20:49:00 UTC (rev 14945)
@@ -191,16 +191,18 @@
const int numQuadPts = _quadrature->numQuadPts();
const double_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
- const double_array& basis = _quadrature->basis();
- const double_array& jacobianDet = _quadrature->jacobianDet();
// Allocate vectors for cell values
double_array orientationCell(numConstraintVert*orientationSize);
double_array stiffnessCell(numConstraintVert);
double_array solutionCell(numCorners*spaceDim);
double_array residualCell(numCorners*spaceDim);
+
+ // Tributary area for the current for each vertex.
+ double_array areaVertexCell(numConstraintVert);
+
+ // Total fault area associated with each vertex (assembled over all cells).
double_array areaCell(numConstraintVert);
- double_array areaAssembledCell(numConstraintVert);
// Get cohesive cells
const ALE::Obj<SieveMesh>& sieveMesh = residual.mesh().sieveMesh();
@@ -237,8 +239,7 @@
_fields->get("area").section();
assert(!areaSection.isNull());
topology::Mesh::RestrictVisitor areaVisitor(*areaSection,
- areaAssembledCell.size(),
- &areaAssembledCell[0]);
+ areaCell.size(), &areaCell[0]);
topology::Field<topology::Mesh>& solution = fields->solution();
const ALE::Obj<RealSection>& solutionSection = solution.section();
@@ -255,15 +256,20 @@
c_iter != cellsCohesiveEnd;
++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
- areaCell = 0.0;
+ _quadrature->retrieveGeometry(c_fault);
+ areaVertexCell = 0.0;
residualCell = 0.0;
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
// Compute contributory area for cell (to weight contributions)
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double dArea = wt*basis[iQuad*numBasis+iBasis];
- areaCell[iBasis] += dArea;
+ areaVertexCell[iBasis] += dArea;
} // for
} // for
@@ -291,9 +297,9 @@
const int indexK = iConstraint + 2*numConstraintVert;
const double pseudoStiffness = stiffnessCell[iConstraint];
- assert(areaAssembledCell[iConstraint] > 0);
+ assert(areaCell[iConstraint] > 0);
const double wt = pseudoStiffness *
- areaCell[iConstraint] / areaAssembledCell[iConstraint];
+ areaVertexCell[iConstraint] / areaCell[iConstraint];
// Get orientation at constraint vertex
const double* orientationVertex =
@@ -1052,6 +1058,7 @@
// Allocate area field.
_fields->add("area", "area");
+
topology::Field<topology::SubMesh>& area = _fields->get("area");
area.newSection(topology::FieldBase::VERTICES_FIELD, 1);
area.allocate();
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-09 19:29:28 UTC (rev 14944)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-09 20:49:00 UTC (rev 14945)
@@ -118,9 +118,9 @@
logEvent = "%sprestep" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- dispTpdt = self.fields.get("disp(t+dt)")
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
for constraint in self.constraints:
- constraint.setFieldIncr(t, t+dt, dispTpdt)
+ constraint.setFieldIncr(t, t+dt, dispIncr)
needNewJacobian = False
for integrator in self.integratorsMesh + self.integratorsSubMesh:
@@ -145,7 +145,7 @@
self._info.log("Solving equations.")
residual = self.fields.get("residual")
- dispIncr = self.fields.get("dispIncr(t->t+dt")
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
self.solver.solve(dispIncr, self.jacobian, residual)
self._logger.eventEnd(logEvent)
@@ -160,11 +160,11 @@
self._logger.eventBegin(logEvent)
dispIncr = self.fields.get("dispIncr(t->t+dt)")
- disp = self.fields.get("disp(t)")
+ dispT = self.fields.get("disp(t)")
dispTmdt = self.fields.get("disp(t-dt)")
dispTmdt.copy(dispT)
- disp += dispIncr
+ dispT += dispIncr
dispIncr.zero()
Formulation.poststep(self, t, dt)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-05-09 19:29:28 UTC (rev 14944)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-05-09 20:49:00 UTC (rev 14945)
@@ -283,7 +283,7 @@
for (int i=0; i < fiberDimE; ++i) {
const int index = iVertex*spaceDim+i;
const double valE = valsE[index] * pseudoStiffness;
- if (valE > tolerance)
+ if (fabs(valE) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
@@ -301,7 +301,7 @@
fault.useSolnIncr(true);
fault.integrateResidual(residual, t, &fields);
- //residual.view("RESIDUAL"); // DEBUGGING
+ residual.view("RESIDUAL"); // DEBUGGING
// Check values
const double* valsE = _data->valsResidualIncr;
@@ -322,7 +322,7 @@
for (int i=0; i < fiberDimE; ++i) {
const int index = iVertex*spaceDim+i;
const double valE = valsE[index] * pseudoStiffness;
- if (valE > tolerance)
+ if (fabs(valE) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
@@ -478,7 +478,7 @@
for (int i=0; i < fiberDimE; ++i) {
const int index = iVertex*spaceDim+i;
const double valE = valsE[index] * pseudoStiffness;
- if (valE > tolerance)
+ if (fabs(valE) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
@@ -517,7 +517,7 @@
for (int i=0; i < fiberDimE; ++i) {
const int index = iVertex*spaceDim+i;
const double valE = valsE[index] * pseudoStiffness;
- if (valE > tolerance)
+ if (fabs(valE) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2009-05-09 19:29:28 UTC (rev 14944)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2009-05-09 20:49:00 UTC (rev 14945)
@@ -131,7 +131,7 @@
integrator.material(&material);
CPPUNIT_ASSERT_EQUAL(false, integrator._useSolnIncr);
try {
- integrator.useSolnIncr(true);
+ integrator.useSolnIncr(false);
// Should have thrown exception, so don't make it here.
CPPUNIT_ASSERT(false);
More information about the CIG-COMMITS
mailing list