[cig-commits] r6879 - in short/3D/PyLith/trunk: libsrc/feassemble
unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Mon May 14 21:44:16 PDT 2007
Author: brad
Date: 2007-05-14 21:44:14 -0700 (Mon, 14 May 2007)
New Revision: 6879
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXY.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXYZ.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearYZ.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData3DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py
Log:
Fixed bugs related to having tranpose of Jacobian instead of Jacobian. All C++ quadrature unit tests now pass.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -117,68 +117,68 @@
_jacobianDet[iQuadPt] = det;
// Compute inverse of Jacobian at quadrature point
- const double d10 =
+ const double d01 =
_jacobian[i00]*_jacobian[i11] -
- _jacobian[i01]*_jacobian[i10];
- const double d21 =
+ _jacobian[i10]*_jacobian[i01];
+ const double d12 =
_jacobian[i10]*_jacobian[i21] -
- _jacobian[i11]*_jacobian[i20];
- const double d20 =
+ _jacobian[i20]*_jacobian[i11];
+ const double d02 =
_jacobian[i00]*_jacobian[i21] -
- _jacobian[i01]*_jacobian[i20];
- if (fabs(d20) > _minJacobian) {
- // Jinv00 = 1/d10 * J11
- // Jinv10 = 1/d10 * -J10
- // Jinv02 = 1/d10 * -J01
- // Jinv11 = 1/d10 * J00
- _jacobianInv[iJ+0] = _jacobian[i11] / d10; // Jinv00
- _jacobianInv[iJ+1] = -_jacobian[i10] / d10; // Jinv10
- _jacobianInv[iJ+2] = -_jacobian[i01] / d10; // Jinv01
- _jacobianInv[iJ+3] = _jacobian[i00] / d10; // Jinv11
- if (fabs(d21) > _minJacobian) {
- // Jinv02 = 1/d21 -J11
- // Jinv12 = 1/d21 J10
- _jacobianInv[iJ+4] = -_jacobian[i11] / d21; // Jinv02
- _jacobianInv[iJ+5] = _jacobian[i10] / d21; // Jinv12
+ _jacobian[i20]*_jacobian[i01];
+ if (fabs(d01) > _minJacobian) {
+ // Jinv00 = 1/d01 * J11
+ // Jinv01 = 1/d01 * -J01
+ // Jinv10 = 1/d01 * -J10
+ // Jinv11 = 1/d01 * J00
+ _jacobianInv[iJ+0] = _jacobian[i11] / d01; // Jinv00
+ _jacobianInv[iJ+1] = -_jacobian[i01] / d01; // Jinv01
+ _jacobianInv[iJ+3] = -_jacobian[i10] / d01; // Jinv10
+ _jacobianInv[iJ+4] = _jacobian[i00] / d01; // Jinv11
+ if (fabs(d12) > _minJacobian) {
+ // Jinv02 = 1/d12 -J11
+ // Jinv12 = 1/d12 J10
+ _jacobianInv[iJ+2] = -_jacobian[i11] / d12; // Jinv02
+ _jacobianInv[iJ+5] = _jacobian[i10] / d12; // Jinv12
- } else if (fabs(d20) > _minJacobian) {
- // Jinv02 = 1/d20 -J01
- // Jinv12 = 1/d20 J00
- _jacobianInv[iJ+4] = -_jacobian[i01] / d20; // Jinv02
- _jacobianInv[iJ+5] = _jacobian[i00] / d20; // Jinv12
+ } else if (fabs(d02) > _minJacobian) {
+ // Jinv02 = 1/d02 -J01
+ // Jinv12 = 1/d02 J00
+ _jacobianInv[iJ+2] = -_jacobian[i01] / d02; // Jinv02
+ _jacobianInv[iJ+5] = _jacobian[i00] / d02; // Jinv12
} else {
- _jacobianInv[iJ+4] = 0.0; // Jinv02
+ _jacobianInv[iJ+2] = 0.0; // Jinv02
_jacobianInv[iJ+5] = 0.0; // Jinv12
} // if/else
- } else if (fabs(d20) > _minJacobian) {
- // Jinv00 = 1/d20 * J21
- // Jinv01 = 1/d20 * -J20
- // Jinv20 = 1/d20 * -J01
- // Jinv21 = 1/d20 * J00
- _jacobianInv[iJ+0] = _jacobian[i21] / d20; // Jinv00
- _jacobianInv[iJ+1] = -_jacobian[i20] / d20; // Jinv10
- _jacobianInv[iJ+4] = -_jacobian[i01] / d20; // Jinv02
- _jacobianInv[iJ+5] = _jacobian[i00] / d20; // Jinv12
- if (fabs(d21) > _minJacobian) {
- // Jinv10 = 1/d21 J21
- // Jinv11 = 1/d21 -J20
- _jacobianInv[iJ+2] = -_jacobian[i21] / d21; // Jinv01
- _jacobianInv[iJ+3] = _jacobian[i20] / d21; // Jinv11
+ } else if (fabs(d02) > _minJacobian) {
+ // Jinv00 = 1/d02 * J21
+ // Jinv02 = 1/d02 * -J01
+ // Jinv10 = 1/d02 * -J20
+ // Jinv12 = 1/d02 * J00
+ _jacobianInv[iJ+0] = _jacobian[i21] / d02; // Jinv00
+ _jacobianInv[iJ+2] = -_jacobian[i01] / d02; // Jinv02
+ _jacobianInv[iJ+3] = -_jacobian[i20] / d02; // Jinv10
+ _jacobianInv[iJ+5] = _jacobian[i00] / d02; // Jinv12
+ if (fabs(d12) > _minJacobian) {
+ // Jinv01 = 1/d12 J21
+ // Jinv11 = 1/d12 -J20
+ _jacobianInv[iJ+1] = -_jacobian[i21] / d12; // Jinv01
+ _jacobianInv[iJ+4] = _jacobian[i20] / d12; // Jinv11
} else {
- _jacobianInv[iJ+2] = 0.0; // Jinv01
- _jacobianInv[iJ+3] = 0.0; // Jinv11
+ _jacobianInv[iJ+1] = 0.0; // Jinv01
+ _jacobianInv[iJ+4] = 0.0; // Jinv11
} // if/else
- } else if (fabs(d21) > _minJacobian) {
+ } else if (fabs(d12) > _minJacobian) {
_jacobianInv[iJ+0] = 0.0; // Jinv00
- _jacobianInv[iJ+1] = 0.0; // Jinv10
- // Jinv01 = 1/d21 J21
- // Jinv11 = 1/d21 -J20
- // Jinv02 = 1/d21 -J11
- // Jinv12 = 1/d21 J10
- _jacobianInv[iJ+2] = _jacobian[i21] / d21; // Jinv01
- _jacobianInv[iJ+3] = -_jacobian[i20] / d21; // Jin11
- _jacobianInv[iJ+4] = -_jacobian[i11] / d21; // Jinv02
- _jacobianInv[iJ+5] = _jacobian[i10] / d21; // Jinv12
+ _jacobianInv[iJ+3] = 0.0; // Jinv10
+ // Jinv01 = 1/d12 J21
+ // Jinv02 = 1/d12 -J11
+ // Jinv11 = 1/d12 -J20
+ // Jinv12 = 1/d12 J10
+ _jacobianInv[iJ+1] = _jacobian[i21] / d12; // Jinv01
+ _jacobianInv[iJ+2] = -_jacobian[i11] / d12; // Jinv02
+ _jacobianInv[iJ+4] = -_jacobian[i20] / d12; // Jinv11
+ _jacobianInv[iJ+5] = _jacobian[i10] / d12; // Jinv12
} else
throw std::runtime_error("Could not invert Jacobian.");
} // for
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXY.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXY.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXY.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -71,7 +71,7 @@
const double pylith::feassemble::QuadratureData2Din3DLinearXY::_jacobianInv[] = {
1.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- 1.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00,
};
pylith::feassemble::QuadratureData2Din3DLinearXY::QuadratureData2Din3DLinearXY(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXYZ.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXYZ.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXYZ.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -70,8 +70,8 @@
};
const double pylith::feassemble::QuadratureData2Din3DLinearXYZ::_jacobianInv[] = {
- 2.90909091e-01, -1.81818182e-01, 1.09090909e-01,
- 1.81818182e-01, -4.32432432e-01, 2.70270270e-01,
+ 2.90909091e-01, 1.09090909e-01, -4.32432432e-01,
+ -1.81818182e-01, 1.81818182e-01, 2.70270270e-01,
};
pylith::feassemble::QuadratureData2Din3DLinearXYZ::QuadratureData2Din3DLinearXYZ(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -70,8 +70,8 @@
};
const double pylith::feassemble::QuadratureData2Din3DLinearXZ::_jacobianInv[] = {
- -1.00000000e+00, 0.00000000e+00, 0.00000000e+00,
- 0.00000000e+00, -0.00000000e+00, 1.00000000e+00,
+ -1.00000000e+00, 0.00000000e+00, -0.00000000e+00,
+ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
};
pylith::feassemble::QuadratureData2Din3DLinearXZ::QuadratureData2Din3DLinearXZ(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearYZ.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearYZ.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearYZ.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -70,8 +70,8 @@
};
const double pylith::feassemble::QuadratureData2Din3DLinearYZ::_jacobianInv[] = {
+ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
- 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
};
pylith::feassemble::QuadratureData2Din3DLinearYZ::QuadratureData2Din3DLinearYZ(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DQuadratic.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DQuadratic.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -106,12 +106,12 @@
};
const double pylith::feassemble::QuadratureData2Din3DQuadratic::_jacobianInv[] = {
- 4.66101695e-02, -4.53389831e-01, 3.00847458e-01,
- -1.99152542e-01, -4.45344130e-02, 4.33198381e-01,
- 7.20338983e-02, -4.27966102e-01, 3.26271186e-01,
- -1.73728814e-01, -6.71936759e-02, 3.99209486e-01,
- 4.91071429e-02, -4.50892857e-01, 3.16964286e-01,
- -1.83035714e-01, -4.68085106e-02, 4.29787234e-01,
+ 4.66101695e-02, 3.00847458e-01, -4.45344130e-02,
+ -4.53389831e-01, -1.99152542e-01, 4.33198381e-01,
+ 7.20338983e-02, 3.26271186e-01, -6.71936759e-02,
+ -4.27966102e-01, -1.73728814e-01, 3.99209486e-01,
+ 4.91071429e-02, 3.16964286e-01, -4.68085106e-02,
+ -4.50892857e-01, -1.83035714e-01, 4.29787234e-01,
};
pylith::feassemble::QuadratureData2Din3DQuadratic::QuadratureData2Din3DQuadratic(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData3DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData3DQuadratic.cc 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData3DQuadratic.cc 2007-05-15 04:44:14 UTC (rev 6879)
@@ -130,7 +130,7 @@
6.33333333e-01, 1.26666667e+00, 2.96666667e+00,
2.36666667e+00, 1.40000000e+00, 5.66666667e-01,
-3.00000000e-01, 3.03333333e+00, 2.66666667e+00,
- 6.33333333e-01, 1.00000000e-00, 2.70000000e+00,
+ 6.33333333e-01, 1.00000000e+00, 2.70000000e+00,
2.63333333e+00, 1.66666667e+00, 8.33333333e-01,
-3.33333333e-02, 3.30000000e+00, 2.66666667e+00,
3.66666667e-01, 7.33333333e-01, 2.96666667e+00,
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py 2007-05-15 00:51:39 UTC (rev 6878)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py 2007-05-15 04:44:14 UTC (rev 6879)
@@ -30,7 +30,7 @@
quadrature.spaceDim, quadrature.cellDim),
dtype=numpy.float64)
jacobianInv = numpy.zeros( (quadrature.numQuadPts,
- quadrature.spaceDim, quadrature.cellDim),
+ quadrature.cellDim, quadrature.spaceDim),
dtype=numpy.float64)
jacobianDet = numpy.zeros( (quadrature.numQuadPts,), dtype=numpy.float64)
@@ -50,52 +50,46 @@
jacobianDet[iQuad] = det
if 1 == quadrature.cellDim:
- jacobianInv[iQuad] = 1.0 / j
+ jacobianInv[iQuad] = 1.0 / j.transpose()
elif 2 == quadrature.cellDim:
minJacobian = 1.0e-06
- jj10 = j[[0,1],:]
- jj21 = j[[1,2],:]
- jj20 = j[[0,2],:]
- det10 = numpy.linalg.det(jj10)
- det21 = numpy.linalg.det(jj21)
- det20 = numpy.linalg.det(jj20)
- if abs(det10) > minJacobian:
- ij10 = numpy.linalg.inv(jj10)
- if abs(det21) > minJacobian:
- ij21 = numpy.linalg.inv(jj21)
- jacobianInv[iQuad] = numpy.array([ [ij10[0,0], ij10[1,0]],
- [ij10[0,1], ij10[1,1]],
- [ij21[0,1], ij21[1,1]] ],
+ jj01 = j[[0,1],:]
+ jj12 = j[[1,2],:]
+ jj02 = j[[0,2],:]
+ det01 = numpy.linalg.det(jj01)
+ det12 = numpy.linalg.det(jj12)
+ det02 = numpy.linalg.det(jj02)
+ if abs(det01) > minJacobian:
+ ij01 = numpy.linalg.inv(jj01)
+ if abs(det12) > minJacobian:
+ ij12 = numpy.linalg.inv(jj12)
+ jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1], ij12[0,1]],
+ [ij01[1,0], ij01[1,1], ij12[1,1]] ],
dtype=numpy.float64)
- elif abs(det20) > minJacobian:
- ij20 = numpy.linalg.inv(jj20)
- jacobianInv[iQuad] = numpy.array([ [ij10[0,0], ij10[1,0]],
- [ij10[0,1], ij10[1,1]],
- [ij20[0,1], ij20[1,1]] ],
+ elif abs(det02) > minJacobian:
+ ij02 = numpy.linalg.inv(jj02)
+ jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,0], ij02[0,1]],
+ [ij01[1,0], ij11[1,1], ij02[1,1]] ],
dtype=numpy.float64)
else:
- jacobianInv[iQuad] = numpy.array([ [ij10[0,0], ij10[0,1]],
- [ij10[1,0], ij10[1,1]],
- [ 0.0, 0.0] ],
+ jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1], 0.0],
+ [ij01[1,0], ij01[1,1], 0.0] ],
dtype=numpy.float64)
- elif abs(det20) > minJacobian:
- ij20 = numpy.linalg.inv(jj20)
- if abs(det21) > minJacobian:
- ij21 = numpy.linalg.inv(jj21)
- jacobianInv[iQuad] = numpy.array([ [ij20[0,0], ij20[1,0]],
- [ij21[0,0], ij21[1,0]],
- [ij20[0,1], ij20[1,1]] ],
+ elif abs(det02) > minJacobian:
+ ij02 = numpy.linalg.inv(jj02)
+ if abs(det12) > minJacobian:
+ ij12 = numpy.linalg.inv(jj12)
+ jacobianInv[iQuad] = numpy.array([ [ij02[0,0], ij12[0,0], ij02[0,1]],
+ [ij02[1,0], ij12[1,0], ij02[1,1]] ],
dtype=numpy.float64)
else:
- jacobianInv[iQuad] = numpy.array([ [ij20[0,0], ij20[1,0]],
- [ 0.0, 0.0],
- [ij20[0,1], ij20[1,1]] ],
+ jacobianInv[iQuad] = numpy.array([ [ij02[0,0], 0.0, ij02[0,1]],
+ [ij02[1,0], 0.0, ij02[1,1]] ],
dtype=numpy.float64)
- elif abs(det21) > minJacobian:
- ij21 = numpy.linalg.inv(jj21)
- jacobianInv[iQuad] = numpy.array([ [ 0.0, 0.0],
- [ij21[0,0], ij21[1,0]],
- [ij21[0,1], ij21[1,1]] ],
+ elif abs(det12) > minJacobian:
+ ij12 = numpy.linalg.inv(jj12)
+ jacobianInv[iQuad] = numpy.array([ [0.0, ij12[0,0], ij12[0,1]],
+ [0.0, ij12[1,0], ij12[1,1]] ],
dtype=numpy.float64)
else:
raise ValueError("Could not find inverse of Jacobian.")
More information about the cig-commits
mailing list