[cig-commits] r16188 - in short/3D/PyLith/trunk: libsrc/faults unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Wed Jan 27 16:14:40 PST 2010
Author: brad
Date: 2010-01-27 16:14:37 -0800 (Wed, 27 Jan 2010)
New Revision: 16188
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py
short/3D/PyLith/trunk/unittests/libtests/faults/data/slipth.py
Log:
More work on unit tests for adjustSolnLumped().
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-01-28 00:14:37 UTC (rev 16188)
@@ -894,20 +894,15 @@
switch(spaceDim)
{ // switch
case 1 : {
- const double d00 = 1.0/Ai[0] + 1.0/Aj[0];
- const double dinv00 = 1.0 / d00;
+ const double Spp = 1.0/Ai[0] + 1.0/Aj[0];
+ const double Sinvpp = 1.0 / Spp;
// Aru = A_i^{-1} r_i - A_j^{-1} r_j + u_i - u_j
const double Aru = ri[0]/Ai[0] - rj[0]/Aj[0] + ui[0] - uj[0];
// dl_k = D^{-1}( C_{ki} Aru - d_k)
const double Aruslip = Aru - slipVertex[0];
- const double dlp = dinv00*Aruslip;
- std::cout << "Aru: " << Aru
- << ", Aruslip: " << Aruslip
- << ", dinv00: " << dinv00
- << ", dlp: " << dlp
- << std::endl;
+ const double dlp = Sinvpp * Aruslip;
// Update displacements at node I
solutionCell[indexI*spaceDim+0] = wt * -1.0/Ai[0] * dlp;
@@ -925,8 +920,6 @@
const double Cpy = orientationVertex[1];
const double Cqx = orientationVertex[2];
const double Cqy = orientationVertex[3];
- const double Crx = orientationVertex[4];
- const double Cry = orientationVertex[5];
const double Spp =
Cpx*Cpx*(1.0/Ai[0] + 1.0/Aj[0]) +
@@ -934,6 +927,7 @@
const double Spq =
Cpx*Cqx*(1.0/Ai[0] + 1.0/Aj[0]) +
Cpy*Cqy*(1.0/Ai[1] + 1.0/Aj[1]);
+ assert(fabs(Spq) < 1.0e-06);
const double Sqq =
Cqx*Cqx*(1.0/Ai[0] + 1.0/Aj[0]) +
Cqy*Cqy*(1.0/Ai[1] + 1.0/Aj[1]);
@@ -954,6 +948,22 @@
const double dlp = Sinvpp*Arupslip + Sinvpq*Aruqslip;
const double dlq = Sinvpq*Arupslip + Sinvqq*Aruqslip;
+ std::cout << "slip: " << slipVertex[0]
+ << ", ri: " << ri[0]
+ << ", rj: " << rj[0]
+ << ", ui: " << ui[0]
+ << ", uj: " << uj[0]
+ << ", Arup: " << Arup
+ << ", Aruq: " << Aruq
+ << ", Arupslip: " << Arupslip
+ << ", Aruqslip: " << Aruqslip
+ << ", Sinvpp: " << Sinvpp
+ << ", Sinvqq: " << Sinvqq
+ << ", dlp: " << dlp
+ << ", dlq: " << dlq
+ << ", wt: " << wt
+ << std::endl;
+
// Update displacements at node I
solutionCell[indexI*spaceDim+0] =
wt * -1.0/Ai[0] * (Cpx*dlp + Cqx*dlq);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2010-01-28 00:14:37 UTC (rev 16188)
@@ -676,6 +676,15 @@
} // for
} // setup disp
+ // compute residual so that slip and residual are setup
+ const double t = 2.134;
+ const double dt = 0.01;
+ fault.timeStep(dt);
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
+ fault.integrateResidual(residual, t, &fields);
+ residual.complete();
+ fault.integrateResidualAssembled(residual, t, &fields);
+
{ // setup disp increment
const ALE::Obj<RealSection>& dispIncrSection = fields.get("dispIncr(t->t+dt)").section();
CPPUNIT_ASSERT(!dispIncrSection.isNull());
@@ -705,15 +714,6 @@
} // setup disp
jacobian.complete();
- // compute residual so that slip and residual are setup
- const double t = 2.134;
- const double dt = 0.01;
- fault.timeStep(dt);
- topology::Field<topology::Mesh>& residual = fields.get("residual");
- fault.integrateResidual(residual, t, &fields);
- residual.complete();
- fault.integrateResidualAssembled(residual, t, &fields);
-
const topology::Field<topology::Mesh>& solution = fields.get("dispIncr(t->t+dt)");
#if 1 // DEBUGGING
solution.view("SOLUTION BEFORE ADJUSTMENT");
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh 2010-01-28 00:14:37 UTC (rev 16188)
@@ -42,6 +42,7 @@
CPPUNIT_TEST( testIntegrateJacobian );
CPPUNIT_TEST( testIntegrateJacobianAssembled );
CPPUNIT_TEST( testIntegrateJacobianAssembledLumped );
+ CPPUNIT_TEST( testAdjustSolnLumped );
CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testCalcTractionsChange );
CPPUNIT_TEST( testSplitField );
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2010-01-28 00:14:37 UTC (rev 16188)
@@ -135,10 +135,10 @@
const double pylith::faults::CohesiveKinDataLine2::_fieldIncrAdjusted[] = {
1.1,
- 1.2, // 3
+ -1.11580527225, // 3
1.3,
- 1.4, // 5
- 1.5, // 6
+ 3.52282149957, // 5
+ 6.59477159896, // 6
};
pylith::faults::CohesiveKinDataLine2::CohesiveKinDataLine2(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2010-01-28 00:14:37 UTC (rev 16188)
@@ -106,8 +106,34 @@
8.8, 9.8, // 9
};
+const double pylith::faults::CohesiveKinDataTri3::_fieldIncr[] = {
+ 3.1, 4.1,
+ 3.2, 4.2, // 3
+ 3.3, 4.3, // 4
+ 3.4, 4.4,
+ 3.5, 4.5, // 6
+ 3.7, 4.7, // 7
+ 3.6, 4.6, // 8
+ 3.8, 4.8, // 9
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_jacobianLumped[] = {
+ 1.1, 1.1,
+ 1.2, 1.2, // 3
+ 1.3, 1.3, // 4
+ 1.4, 1.4,
+ 1.5, 1.5, // 6
+ 1.7, 1.7, // 7
+ 1.6, 1.6, // 8
+ 1.8, 1.8, // 9
+};
+
const int pylith::faults::CohesiveKinDataTri3::_numConstraintVert = 2;
+const int pylith::faults::CohesiveKinDataTri3::_constraintVertices[] = {
+ 8, 9
+};
+
const double pylith::faults::CohesiveKinDataTri3::_orientation[] = {
0.0, -1.0, -1.0, 0.0,
0.0, -1.0, -1.0, 0.0
@@ -118,10 +144,6 @@
1.0,
};
-const int pylith::faults::CohesiveKinDataTri3::_constraintVertices[] = {
- 8, 9
-};
-
const double pylith::faults::CohesiveKinDataTri3::_residual[] = {
0.0, 0.0,
-9.6, -8.6, // 3
@@ -277,6 +299,17 @@
0.0, 0.0,
};
+const double pylith::faults::CohesiveKinDataTri3::_fieldIncrAdjusted[] = {
+ 3.1, 4.1,
+ 3.2, 4.2, // 3
+ 3.3, 4.3, // 4
+ 3.4, 4.4,
+ 3.5, 4.5, // 6
+ 3.7, 4.7, // 7
+ 3.6, 4.6, // 8
+ 3.8, 4.8, // 9
+};
+
pylith::faults::CohesiveKinDataTri3::CohesiveKinDataTri3(void)
{ // constructor
meshFilename = const_cast<char*>(_meshFilename);
@@ -296,12 +329,15 @@
riseTimeFilename = const_cast<char*>(_riseTimeFilename);
matPropsFilename = const_cast<char*>(_matPropsFilename);
fieldT = const_cast<double*>(_fieldT);
+ fieldIncr = const_cast<double*>(_fieldIncr);
+ jacobianLumped = const_cast<double*>(_jacobianLumped);
orientation = const_cast<double*>(_orientation);
area = const_cast<double*>(_area);
- constraintVertices = const_cast<int*>(_constraintVertices);
+ residualIncr = const_cast<double*>(_residualIncr);
residual = const_cast<double*>(_residual);
- residualIncr = const_cast<double*>(_residualIncr);
jacobian = const_cast<double*>(_jacobian);
+ fieldIncrAdjusted = const_cast<double*>(_fieldIncrAdjusted);
+ constraintVertices = const_cast<int*>(_constraintVertices);
numConstraintVert = _numConstraintVert;
} // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh 2010-01-28 00:14:37 UTC (rev 16188)
@@ -57,14 +57,17 @@
static const char* _matPropsFilename; ///< Name of db for bulk mat properties.
//@}
- static const double _fieldT[]; ///< Solution field at time t.
+ static const double _fieldT[]; ///< Field over domain at time t.
+ static const double _fieldIncr[]; ///< Solution increment field over domain at time t.
+ static const double _jacobianLumped[]; ///< Lumped Jacobian.
static const double _orientation[]; ///< Expected values for fault orientation.
static const double _area[]; ///< Expected values for fault area.
- static const int _constraintVertices[]; ///< Expected points for constraint vertices
static const double _residual[]; ///< Expected values from residual calculation.
- static const double _residualIncr[]; ///< Expected values from residual calculation using solution increment.
+ static const double _residualIncr[]; ///< Expected values from residual calculation with solution increment.
static const double _jacobian[]; ///< Expected values from Jacobian calculation.
+ static const double _fieldIncrAdjusted[]; ///< Expected values for colution increment after adjustment.
+ static const int _constraintVertices[]; ///< Expected points for constraint vertices
static const int _numConstraintVert; ///< Number of constraint vertices
};
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py 2010-01-28 00:14:37 UTC (rev 16188)
@@ -1,15 +1,76 @@
-# Line2
-C = 1.0
-dk = 1.89546413727; dl = 1.5
-Ai = 2.2; ri = 7.5; ui = 7.2; dui = 1.2
-Aj = 2.4; rj = -7.5; uj = 7.4; duj = 1.4
+test = "tri3"
-Si = (Ai * Aj) / (Ai + Aj)
-ddl = Si * (C*(ri/Ai - rj/Aj + ui - uj) - dk)
+# ----------------------------------------------------------------------
+if test == "line2":
+ C = 1.0
+ dk = 1.89546413727; l = 1.5
+ Ai = 2.2; ri = 7.5; ui = 7.2; dui = 1.2
+ Aj = 2.4; rj = -7.5; uj = 7.4; duj = 1.4
-ddui = -C / Ai * ddl
-dduj = +C / Aj * ddl
+ Si = (Ai * Aj) / (Ai + Aj)
+ Aru = ri/Ai - rj/Aj + ui - uj
+ Aruslip = C*Aru - dk
+ dlp = Si * Aruslip
-print ddl
-print ddui
-print dui+ddui,duj+dduj,dl+ddl
+ ddui = -C / Ai * dlp
+ dduj = +C / Aj * dlp
+
+ #print "Aru",Aru
+ #print "Aruslip",Aruslip
+ #print "Si",Si
+ #print "dlip",dlp
+
+ print dui+ddui,duj+dduj,l+dlp
+
+# ----------------------------------------------------------------------
+elif test == "tri3":
+
+ # Lagrange vertex 8, vertex i: 3, vertex j: 6
+ Cpx = 0.0
+ Cpy = -1.0
+ Cqx = -1.0
+ Cqy = 0.0
+
+ dkp = 1.89546413727; lp = 1.5
+ dkq = 0.08241148423; lq = 1.5
+
+ # vertex i
+ Aix = 1.2; Aiy = 1.2
+ rix = -9.6; riy = -8.6
+ uix = 8.2; uiy = 9.2
+ duix = 3.2; duiy = 4.2
+
+ # vertex j
+ Ajx = 1.5; Ajy = 1.5
+ rjx = +9.6; rjy = +8.6
+ ujx = 8.5; ujy = 9.5
+ dujx = 3.5; dujy = 4.5
+
+ Sppi = Aix*Ajx / (Aix + Ajx)
+ Sqqi = Aix*Ajx / (Aix + Ajx)
+
+ Arux = rix / Aix - rjx / Ajx + uix - ujx
+ Aruy = riy / Aiy - rjy / Ajy + uiy - ujy
+
+ Arup = Cpx*Arux + Cpy*Aruy
+ Aruq = Cqx*Arux + Cqy*Aruy
+ Arupslip = Arup - dkp
+ Aruqslip = Aruq - dkq
+
+ dlp = Sppi * Arupslip
+ dlq = Sqqi * Aruqslip
+
+ dduix = -1.0/Aix * (Cpx*dlp + Cqx*dlq)
+ dduiy = -1.0/Aiy * (Cpy*dlp + Cqy*dlq)
+
+ ddujx = +1.0/Ajx * (Cpx*dlp + Cqx*dlq)
+ ddujy = +1.0/Ajy * (Cpy*dlp + Cqy*dlq)
+
+ #print "Aru",Aru
+ #print "Aruslip",Aruslip
+ #print "Si",Si
+ #print "dlip",dlp
+
+ print duix+dduix,duiy+dduiy
+ print dujx+ddujx,dujy+ddujy
+ print lp+dlp,lq+dlq
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/slipth.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/slipth.py 2010-01-27 20:35:42 UTC (rev 16187)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/slipth.py 2010-01-28 00:14:37 UTC (rev 16188)
@@ -1,6 +1,6 @@
from math import *
-finalSlip = (2.3, 0.0, 0.0)
+finalSlip = (2.3, 0.1, 0.0)
riseTime = 1.4
slipMag = (finalSlip[0]**2+finalSlip[1]**2+finalSlip[2]**2)**0.5
More information about the CIG-COMMITS
mailing list