[cig-commits] r8287 - in short/3D/PyLith/trunk: libsrc/bc pylith/bc
unittests/libtests/bc/data
brad at geodynamics.org
brad at geodynamics.org
Tue Nov 13 10:55:28 PST 2007
Author: brad
Date: 2007-11-13 10:55:28 -0800 (Tue, 13 Nov 2007)
New Revision: 8287
Modified:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
Log:
Fixed integration of residual term for absorbing dampers (does not depend on disp at t+dt). Updated unit tests accordingly.
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-11-13 18:55:28 UTC (rev 8287)
@@ -288,8 +288,7 @@
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] +=
dampingConstsCell[iQuad*spaceDim+iDim] *
- valIJ * (-dispTpdtCell[jBasis*spaceDim+iDim]
- + dispTmdtCell[jBasis*spaceDim+iDim]);
+ valIJ * (dispTmdtCell[jBasis*spaceDim+iDim]);
} // for
} // for
Modified: short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py 2007-11-13 18:55:28 UTC (rev 8287)
@@ -86,12 +86,13 @@
Verify compatibility of configuration.
"""
BoundaryCondition.verifyConfiguration(self)
- if self.quadrature.spaceDim != self.mesh.dimension-1:
+ if self.quadrature.cellDim != self.mesh.dimension()-1:
raise ValueError, \
"Quadrature scheme and mesh are incompatible.\n" \
"Dimension for quadrature: %d\n" \
"Dimension of mesh boundary '%s': %d" % \
- (self.quadrature.spaceDim, self.label, self.mesh.dimension-1)
+ (self.quadrature.cellDim,
+ self.label, self.mesh.dimension()-1)
return
@@ -102,6 +103,7 @@
Setup members using inventory.
"""
BoundaryCondition._configure(self)
+ self.quadrature = self.inventory.quadrature
return
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc 2007-11-13 18:55:28 UTC (rev 8287)
@@ -122,18 +122,18 @@
7.5e+06, 1.25e+07, 7.5e+06,
};
const double pylith::bc::AbsorbingDampersDataHex8::_valsResidual[] = {
- -5.50000000e+06, 1.29166667e+07, -1.10000000e+07,
- 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- -1.50000000e+07, 2.25000000e+07, -3.00000000e+07,
- 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- -9.50000000e+06, 9.58333333e+06, -1.90000000e+07,
- 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- -8.50000000e+06, 1.04166667e+07, -1.70000000e+07,
- 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- -2.10000000e+07, 1.75000000e+07, -4.20000000e+07,
- 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- -1.25000000e+07, 7.08333333e+06, -2.50000000e+07,
- 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 9.50000021e+06, 2.33333337e+07, 2.65000006e+07,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 2.10000004e+07, 4.00000007e+07, 5.70000012e+07,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 1.15000002e+07, 1.66666670e+07, 3.05000006e+07,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 1.10000002e+07, 1.83333337e+07, 2.95000006e+07,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 2.40000005e+07, 3.00000007e+07, 6.30000012e+07,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 1.30000002e+07, 1.16666670e+07, 3.35000006e+07,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
};
const double pylith::bc::AbsorbingDampersDataHex8::_valsJacobian[] = {
3.33333333e+06, 0.0, 0.0, // 0x
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc 2007-11-13 18:55:28 UTC (rev 8287)
@@ -67,9 +67,9 @@
17.5e+6,
};
const double pylith::bc::AbsorbingDampersDataLine2::_valsResidual[] = {
- -12.5e+6*0.2/0.5,
+ 12.5e+6*1.0/0.5,
0.0,
- -17.5e+6*0.6/0.5,
+ 17.5e+6*1.2/0.5,
};
const double pylith::bc::AbsorbingDampersDataLine2::_valsJacobian[] = {
12.5e+6/0.5, 0.0, 0.0,
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc 2007-11-13 18:55:28 UTC (rev 8287)
@@ -78,12 +78,12 @@
1.25e+07, 7.5e+06,
};
const double pylith::bc::AbsorbingDampersDataQuad4::_valsResidual[] = {
- -7.5e+06, 1.5e+06,
- -7.5e+06, 1.5e+06,
+ 2.625e+07, 3.15e+07,
+ 2.625e+07, 3.15e+07,
0.0, 0.0,
0.0, 0.0,
- -2.75e+07, -1.35e+07,
- -2.75e+07, -1.35e+07,
+ 3.625e+07, 3.0e+07,
+ 3.625e+07, 3.0e+07,
};
const double pylith::bc::AbsorbingDampersDataQuad4::_valsJacobian[] = {
1.25e+07, 0.0, // 0x
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc 2007-11-13 18:55:28 UTC (rev 8287)
@@ -73,10 +73,10 @@
1.25e+07, 7.5e+06, 7.5e+06
};
const double pylith::bc::AbsorbingDampersDataTet4::_valsResidual[] = {
- -2.22222222e+06, 8.33333333e+05, -2.66666667e+06,
+ 4.861111111e+06, 5.83333333e+06, 8.33333333e+06,
0.0, 0.0, 0.0,
- -2.22222222e+06, 8.33333333e+05, -2.66666667e+06,
- -2.22222222e+06, 8.33333333e+05, -2.66666667e+06,
+ 4.861111111e+06, 5.83333333e+06, 8.33333333e+06,
+ 4.861111111e+06, 5.83333333e+06, 8.33333333e+06,
};
const double pylith::bc::AbsorbingDampersDataTet4::_valsJacobian[] = {
1.38888889e+06, 0.0, 0.0, // 0x
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc 2007-11-13 18:55:28 UTC (rev 8287)
@@ -72,9 +72,9 @@
};
const double pylith::bc::AbsorbingDampersDataTri3::_valsResidual[] = {
0.0, 0.0,
- -1.2e+07, -2.0e+06,
+ 2.4e+07, 1.0e+07,
0.0, 0.0,
- -1.2e+07, -2.0e+06,
+ 2.4e+07, 1.0e+07,
};
const double pylith::bc::AbsorbingDampersDataTri3::_valsJacobian[] = {
0.0, 0.0, // 0x
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py 2007-11-13 18:32:00 UTC (rev 8286)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py 2007-11-13 18:55:28 UTC (rev 8287)
@@ -31,16 +31,16 @@
edgeLen = 2.0**0.5
N0 = 0.5
N1 = 0.5
- velocityX = 0.5*(N0+N1)*(N0*(1.5-1.1)+N1*(2.1-1.3)) / (2.0*dt)
- velocityY = 0.5*(N0+N1)*(N0*(2.4-1.8)+N1*(2.4-2.2)) / (2.0*dt)
+ dispX = 0.5*(N0+N1)*(N0*(1.1)+N1*(1.3)) / (2.0*dt)
+ dispY = 0.5*(N0+N1)*(N0*(1.8)+N1*(2.2)) / (2.0*dt)
normal = [0.5**0.5, -0.5**0.5]
constNormal = density*vp
constTangential = density*vs
dampingConsts = [abs(constNormal*normal[0] - constTangential*normal[1]),
abs(constNormal*normal[1] + constTangential*normal[0])]
- residualX = -dampingConsts[0]*velocityX*edgeLen
- residualY = -dampingConsts[1]*velocityY*edgeLen
+ residualX = dampingConsts[0]*dispX*edgeLen
+ residualY = dampingConsts[1]*dispY*edgeLen
residual = [residualX, residualY,
residualX, residualY]
@@ -78,12 +78,12 @@
dt = 0.25
N0 = 0.5
N1 = 0.5
- velocity1X = 0.5*(N0+N1)*(N0*(1.5-1.1) + N1*(1.2-1.0)) / (2.0*dt)
- velocity1Y = 0.5*(N0+N1)*(N0*(2.4-1.8) + N1*(1.6-2.4)) / (2.0*dt)
+ disp1X = 0.5*(N0+N1)*(N0*(1.1) + N1*(1.0)) / (2.0*dt)
+ disp1Y = 0.5*(N0+N1)*(N0*(1.8) + N1*(2.4)) / (2.0*dt)
normal1 = [-1.0, 0.0]
- velocity2X = (N0*(2.4-1.4) + N1*(2.7-1.5)) / (2.0*dt)
- velocity2Y = (N0*(2.4-2.4) + N1*(3.4-1.6)) / (2.0*dt)
+ disp2X = 0.5*(N0+N1)*(N0*(1.4) + N1*(1.5)) / (2.0*dt)
+ disp2Y = 0.5*(N0+N1)*(N0*(2.4) + N1*(1.6)) / (2.0*dt)
normal2 = [1.0, 0.0]
edgeLen = 2.0
@@ -94,10 +94,10 @@
abs(constNormal*normal1[1] + constTangential*normal1[0]),
abs(constNormal*normal2[0] - constTangential*normal2[1]),
abs(constNormal*normal2[1] + constTangential*normal2[0])]
- residual1X = -dampingConsts[0]*velocity1X*edgeLen
- residual1Y = -dampingConsts[1]*velocity1Y*edgeLen
- residual2X = -dampingConsts[2]*velocity2X*edgeLen
- residual2Y = -dampingConsts[3]*velocity2Y*edgeLen
+ residual1X = dampingConsts[0]*disp1X*edgeLen
+ residual1Y = dampingConsts[1]*disp1Y*edgeLen
+ residual2X = dampingConsts[2]*disp2X*edgeLen
+ residual2Y = dampingConsts[3]*disp2Y*edgeLen
residual = [residual1X, residual1Y,
residual1X, residual1Y,
residual2X, residual2Y,
@@ -119,9 +119,9 @@
print "damping constants:"
for v in dampingConsts:
print " %16.8e" % v
- print "velocity:"
- print " velocity1: ",velocity1X," ",velocity1Y
- print " velocity2: ",velocity2X," ",velocity2Y
+ print "disp:"
+ print " disp1: ",disp1X," ",disp1Y
+ print " disp2: ",disp2X," ",disp2Y
print "values for residual:"
for v in residual:
print " %16.8e" % v
@@ -145,12 +145,9 @@
N0 = 1.0/3.0
N1 = 1.0/3.0
N2 = 1.0/3.0
- velocityX = (N0+N1+N2)/3.0 * \
- (N0*(2.1-1.3)+N1*(1.8-1.2)+N2*(1.2-1.0)) / (2.0*dt)
- velocityY = (N0+N1+N2)/3.0 * \
- (N0*(2.4-2.2)+N1*(2.0-2.4)+N2*(1.6-2.4)) / (2.0*dt)
- velocityZ = (N0+N1+N2)/3.0 * \
- (N0*(5.2-3.6)+N1*(4.6-3.4)+N2*(3.4-3.0)) / (2.0*dt)
+ dispX = (N0+N1+N2)/3.0 * (N0*(1.3)+N1*(1.2)+N2*(1.0)) / (2.0*dt)
+ dispY = (N0+N1+N2)/3.0 * (N0*(2.2)+N1*(2.4)+N2*(2.4)) / (2.0*dt)
+ dispZ = (N0+N1+N2)/3.0 * (N0*(3.6)+N1*(3.4)+N2*(3.0)) / (2.0*dt)
normal = [-1.0, 0.0, 0.0]
tangent1 = [0.0, -1.0, 0.0]
tangent2 = [0.0, 0.0, 1.0]
@@ -166,9 +163,9 @@
abs(constNormal*normal[2] +
constTangential*tangent1[2] +
constTangential*tangent2[2])]
- residualX = -dampingConsts[0]*velocityX*area
- residualY = -dampingConsts[1]*velocityY*area
- residualZ = -dampingConsts[2]*velocityZ*area
+ residualX = dampingConsts[0]*dispX*area
+ residualY = dampingConsts[1]*dispY*area
+ residualZ = dampingConsts[2]*dispZ*area
residual = [residualX, residualY, residualZ,
residualX, residualY, residualZ,
residualX, residualY, residualZ]
@@ -224,18 +221,6 @@
[0.16666667, 0.0446582, 0.16666667, 0.62200847]]
cells = [[2, 8, 6, 0],
[4, 10, 8, 2]]
- dispTpdt = [[1.2, 1.1, 3.4],
- [1.5, 1.0, 4.0],
- [1.8, 0.9, 4.6],
- [2.1, 0.8, 5.2],
- [2.4, 0.7, 5.8],
- [2.7, 0.6, 6.4],
- [3.0, 0.5, 7.0],
- [3.3, 0.4, 7.6],
- [3.6, 0.3, 8.2],
- [3.9, 0.2, 8.8],
- [4.2, 0.1, 9.4],
- [4.5, 0.0, 10.0]]
dispTmdt = [[1.0, 2.4, 3.0],
[1.1, 2.2, 3.2],
[1.2, 2.0, 3.4],
@@ -272,25 +257,25 @@
N1 = b[1]
N2 = b[2]
N3 = b[3]
- v0x = (dispTpdt[cell[0]][0]-dispTmdt[cell[0]][0])/(2.0*dt)
- v0y = (dispTpdt[cell[0]][1]-dispTmdt[cell[0]][1])/(2.0*dt)
- v0z = (dispTpdt[cell[0]][2]-dispTmdt[cell[0]][2])/(2.0*dt)
- v1x = (dispTpdt[cell[1]][0]-dispTmdt[cell[1]][0])/(2.0*dt)
- v1y = (dispTpdt[cell[1]][1]-dispTmdt[cell[1]][1])/(2.0*dt)
- v1z = (dispTpdt[cell[1]][2]-dispTmdt[cell[1]][2])/(2.0*dt)
- v2x = (dispTpdt[cell[2]][0]-dispTmdt[cell[2]][0])/(2.0*dt)
- v2y = (dispTpdt[cell[2]][1]-dispTmdt[cell[2]][1])/(2.0*dt)
- v2z = (dispTpdt[cell[2]][2]-dispTmdt[cell[2]][2])/(2.0*dt)
- v3x = (dispTpdt[cell[3]][0]-dispTmdt[cell[3]][0])/(2.0*dt)
- v3y = (dispTpdt[cell[3]][1]-dispTmdt[cell[3]][1])/(2.0*dt)
- v3z = (dispTpdt[cell[3]][2]-dispTmdt[cell[3]][2])/(2.0*dt)
- velocityX = N0*v0x + N1*v1x + N2*v2x + N3*v3x
- velocityY = N0*v0y + N1*v1y + N2*v2y + N3*v3y
- velocityZ = N0*v0z + N1*v1z + N2*v2z + N3*v3z
+ d0x = dispTmdt[cell[0]][0]/(2.0*dt)
+ d0y = dispTmdt[cell[0]][1]/(2.0*dt)
+ d0z = dispTmdt[cell[0]][2]/(2.0*dt)
+ d1x = dispTmdt[cell[1]][0]/(2.0*dt)
+ d1y = dispTmdt[cell[1]][1]/(2.0*dt)
+ d1z = dispTmdt[cell[1]][2]/(2.0*dt)
+ d2x = dispTmdt[cell[2]][0]/(2.0*dt)
+ d2y = dispTmdt[cell[2]][1]/(2.0*dt)
+ d2z = dispTmdt[cell[2]][2]/(2.0*dt)
+ d3x = dispTmdt[cell[3]][0]/(2.0*dt)
+ d3y = dispTmdt[cell[3]][1]/(2.0*dt)
+ d3z = dispTmdt[cell[3]][2]/(2.0*dt)
+ dispX = N0*d0x + N1*d1x + N2*d2x + N3*d3x
+ dispY = N0*d0y + N1*d1y + N2*d2y + N3*d3y
+ dispZ = N0*d0z + N1*d1z + N2*d2z + N3*d3z
- residualX = -dampingConsts[0] * velocityX * area * jacobianDet
- residualY = -dampingConsts[1] * velocityY * area * jacobianDet
- residualZ = -dampingConsts[2] * velocityZ * area * jacobianDet
+ residualX = dampingConsts[0] * dispX * area * jacobianDet
+ residualY = dampingConsts[1] * dispY * area * jacobianDet
+ residualZ = dampingConsts[2] * dispZ * area * jacobianDet
residual[cell[0],:] += N0*numpy.array([residualX,residualY,residualZ])
residual[cell[1],:] += N1*numpy.array([residualX,residualY,residualZ])
residual[cell[2],:] += N2*numpy.array([residualX,residualY,residualZ])
@@ -359,7 +344,7 @@
#calcTri3()
#calcQuad4()
#calcTet4()
-#calcHex8()
+calcHex8()
# End of file
More information about the cig-commits
mailing list