[cig-commits] r20627 - in short/3D/PyLith/branches/v1.7-trunk: libsrc/pylith/feassemble unittests/libtests/feassemble unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Fri Aug 24 16:05:59 PDT 2012
Author: brad
Date: 2012-08-24 16:05:58 -0700 (Fri, 24 Aug 2012)
New Revision: 20627
Modified:
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DLinear.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DQuadratic.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DLinear.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DQuadratic.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DLinear.py
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DQuadratic.py
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DLinear.py
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DQuadratic.py
Log:
Finished work on C++ unit tests for calculating stable explicit time step.
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -1020,8 +1020,7 @@
// ----------------------------------------------------------------------
// Compute volume of tetrahedral cell.
PylithScalar
-pylith::feassemble::ElasticityExplicitTet4::_volume(
- const scalar_array& coordinatesCell) const
+pylith::feassemble::ElasticityExplicitTet4::_volume(const scalar_array& coordinatesCell) const
{ // __volume
assert(12 == coordinatesCell.size());
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -277,11 +277,133 @@
PetscLogFlops(numEdges*9);
- // TODO: Add radius of inscribed sphere
- // r = V / (3*A)
+ // Radius of inscribed sphere
+ const PylithScalar v = volume(coordinatesCell);
+ const PylithScalar a = faceArea(coordinatesCell, 0) +
+ faceArea(coordinatesCell, 1) +
+ faceArea(coordinatesCell, 2) +
+ faceArea(coordinatesCell, 3);
+
+ const PylithScalar r = v / (3.0*a);
+ if (r < minWidth) {
+ minWidth = r;
+ } // if
+ PetscLogFlops(5);
+
return minWidth;
} // minCellWidth
+// ----------------------------------------------------------------------
+// Compute cell volume.
+PylithScalar
+pylith::feassemble::GeometryTet3D::volume(const scalar_array& coordinatesCell) const
+{ // volume
+ assert(12 == coordinatesCell.size() ||
+ 30 == coordinatesCell.size()); // :KLUDGE: allow quadratic
+
+ const PylithScalar x0 = coordinatesCell[0];
+ const PylithScalar y0 = coordinatesCell[1];
+ const PylithScalar z0 = coordinatesCell[2];
+
+ const PylithScalar x1 = coordinatesCell[3];
+ const PylithScalar y1 = coordinatesCell[4];
+ const PylithScalar z1 = coordinatesCell[5];
+
+ const PylithScalar x2 = coordinatesCell[6];
+ const PylithScalar y2 = coordinatesCell[7];
+ const PylithScalar z2 = coordinatesCell[8];
+
+ const PylithScalar x3 = coordinatesCell[9];
+ const PylithScalar y3 = coordinatesCell[10];
+ const PylithScalar z3 = coordinatesCell[11];
+ const PylithScalar det =
+ x1*(y2*z3-y3*z2)-y1*(x2*z3-x3*z2)+(x2*y3-x3*y2)*z1 -
+ x0*((y2*z3-y3*z2)-y1*(z3-z2)+(y3-y2)*z1) +
+ y0*((x2*z3-x3*z2)-x1*(z3-z2)+(x3-x2)*z1) -
+ z0*((x2*y3-x3*y2)-x1*(y3-y2)+(x3-x2)*y1);
+
+ const PylithScalar v = det / 6.0;
+ PetscLogFlops(48);
+
+ return v;
+} // volume
+
+// ----------------------------------------------------------------------
+// Compute area of face.
+PylithScalar
+pylith::feassemble::GeometryTet3D::faceArea(const scalar_array& coordinatesCell,
+ const int face) const
+{ // faceArea
+ assert(12 == coordinatesCell.size() ||
+ 30 == coordinatesCell.size()); // :KLUDGE: allow quadratic
+
+ const PylithScalar x0 = coordinatesCell[0];
+ const PylithScalar y0 = coordinatesCell[1];
+ const PylithScalar z0 = coordinatesCell[2];
+
+ const PylithScalar x1 = coordinatesCell[3];
+ const PylithScalar y1 = coordinatesCell[4];
+ const PylithScalar z1 = coordinatesCell[5];
+
+ const PylithScalar x2 = coordinatesCell[6];
+ const PylithScalar y2 = coordinatesCell[7];
+ const PylithScalar z2 = coordinatesCell[8];
+
+ const PylithScalar x3 = coordinatesCell[9];
+ const PylithScalar y3 = coordinatesCell[10];
+ const PylithScalar z3 = coordinatesCell[11];
+
+ PylithScalar a[3];
+ PylithScalar b[3];
+ switch (face) {
+ case 0:
+ a[0] = x3-x1;
+ a[1] = y3-y1;
+ a[2] = z3-z1;
+ b[0] = x2-x1;
+ b[1] = y2-y1;
+ b[2] = z2-z1;
+ break;
+ case 1:
+ a[0] = x3-x2;
+ a[1] = y3-y2;
+ a[2] = z3-z2;
+ b[0] = x0-x2;
+ b[1] = y0-y2;
+ b[2] = z0-z2;
+ break;
+ case 2:
+ a[0] = x3-x0;
+ a[1] = y3-y0;
+ a[2] = z3-z0;
+ b[0] = x1-x0;
+ b[1] = y1-y0;
+ b[2] = z1-z0;
+ break;
+ case 3:
+ a[0] = x1-x0;
+ a[1] = y1-y0;
+ a[2] = z1-z0;
+ b[0] = x2-x0;
+ b[1] = y2-y0;
+ b[2] = z2-z0;
+ break;
+ default:
+ assert(0);
+ throw std::logic_error("Unknown face.");
+ } // switch
+
+ const PylithScalar areaX = a[1]*b[2] - a[2]*b[1];
+ const PylithScalar areaY = a[2]*b[0] - a[0]*b[2];
+ const PylithScalar areaZ = a[0]*b[1] - a[1]*b[0];
+
+ const PylithScalar area = 0.5*sqrt(areaX*areaX + areaY*areaY + areaZ*areaZ);
+ PetscLogFlops(22);
+
+ return area;
+} // faceArea
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh 2012-08-24 23:05:58 UTC (rev 20627)
@@ -116,6 +116,22 @@
*/
PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+ /** Compute cell volume.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Volume of cell.
+ */
+ PylithScalar volume(const scalar_array& coordinatesCell) const;
+
+ /** Compute area of face.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @param face Index of vertex across from face.
+ * @returns Area of cell face.
+ */
+ PylithScalar faceArea(const scalar_array& coordinatesCell,
+ const int face) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -215,7 +215,6 @@
PetscLogFlops(numEdges*6);
-#if 1
// Ad-hoc to account for distorted cells.
// Radius of inscribed circle.
const PylithScalar xA = coordinatesCell[0];
@@ -234,7 +233,6 @@
minWidth = r;
PetscLogFlops(3*6 + 3 + 8);
-#endif
return minWidth;
} // minCellWidth
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -180,7 +180,7 @@
const int size = residualSection->sizeWithBC();
CPPUNIT_ASSERT_EQUAL(sizeE, size);
-#if 0
+#if 0 // DEBUGGING
residual.view("RESIDUAL");
std::cout << "EXPECTED RESIDUAL" << std::endl;
for (int i=0; i < size; ++i)
@@ -220,7 +220,7 @@
const int size = residualSection->sizeWithBC();
CPPUNIT_ASSERT_EQUAL(sizeE, size);
-#if 0
+#if 0 // DEBUGGING
residual.view("RESIDUAL");
std::cout << "EXPECTED RESIDUAL" << std::endl;
for (int i=0; i < size; ++i)
@@ -372,7 +372,8 @@
_initialize(&mesh, &integrator, &fields);
const PylithScalar dtStable = integrator.stableTimeStep(mesh);
- CPPUNIT_ASSERT_EQUAL(1.0, dtStable/_data->dtStableExplicit);
+ const PylithScalar tolerance = 1.0e-6;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dtStable/_data->dtStableExplicit, tolerance);
} // testStableTimeStep
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitData2DLinear::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitData2DLinear::_dtStableExplicit = 4.45596817e-05;
+const PylithScalar pylith::feassemble::ElasticityExplicitData2DLinear::_dtStableExplicit = 5.09461158e-05;
const PylithScalar pylith::feassemble::ElasticityExplicitData2DLinear::_gravityVec[] = {
0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitData2DQuadratic::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitData2DQuadratic::_dtStableExplicit = 7.34858411e-05;
+const PylithScalar pylith::feassemble::ElasticityExplicitData2DQuadratic::_dtStableExplicit = 9.32914833e-05;
const PylithScalar pylith::feassemble::ElasticityExplicitData2DQuadratic::_gravityVec[] = {
0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitData3DLinear::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitData3DLinear::_dtStableExplicit = 1.11111111e-04;
+const PylithScalar pylith::feassemble::ElasticityExplicitData3DLinear::_dtStableExplicit = 3.78899352e-06;
const PylithScalar pylith::feassemble::ElasticityExplicitData3DLinear::_gravityVec[] = {
0.00000000e+00, 0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitData3DQuadratic::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitData3DQuadratic::_dtStableExplicit = 1.11111111e-04;
+const PylithScalar pylith::feassemble::ElasticityExplicitData3DQuadratic::_dtStableExplicit = 9.81816592e-06;
const PylithScalar pylith::feassemble::ElasticityExplicitData3DQuadratic::_gravityVec[] = {
0.00000000e+00, 0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DLinear.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DLinear.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DLinear.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DLinear::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DLinear::_dtStableExplicit = 4.45596817e-05;
+const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DLinear::_dtStableExplicit = 5.09461158e-05;
const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DLinear::_gravityVec[] = {
0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DQuadratic.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DQuadratic.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData2DQuadratic.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DQuadratic::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DQuadratic::_dtStableExplicit = 7.34858411e-05;
+const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DQuadratic::_dtStableExplicit = 9.32914833e-05;
const PylithScalar pylith::feassemble::ElasticityExplicitGravData2DQuadratic::_gravityVec[] = {
0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DLinear.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DLinear.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DLinear.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DLinear::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DLinear::_dtStableExplicit = 1.11111111e-04;
+const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DLinear::_dtStableExplicit = 3.78899352e-06;
const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DLinear::_gravityVec[] = {
0.00000000e+00, 0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DQuadratic.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DQuadratic.cc 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/ElasticityExplicitGravData3DQuadratic.cc 2012-08-24 23:05:58 UTC (rev 20627)
@@ -43,7 +43,7 @@
const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DQuadratic::_dt = 1.00000000e-02;
-const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DQuadratic::_dtStableExplicit = 1.11111111e-04;
+const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DQuadratic::_dtStableExplicit = 9.81816592e-06;
const PylithScalar pylith::feassemble::ElasticityExplicitGravData3DQuadratic::_gravityVec[] = {
0.00000000e+00, 0.00000000e+00, -1.00000000e+08,
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DLinear.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DLinear.py 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DLinear.py 2012-08-24 23:05:58 UTC (rev 20627)
@@ -60,7 +60,7 @@
a = (0.1**2 + 0.9**2)**0.5
b = (1.3**2 + 0.7**2)**0.5
- c = (1.0**2 + 0.2**2)**0.5
+ c = (1.2**2 + 0.2**2)**0.5
k = 0.5 * (a + b + c)
r = (k*(k-a)*(k-b)*(k-c))**0.5 / k
self.minCellWidth = r
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DQuadratic.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DQuadratic.py 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh2DQuadratic.py 2012-08-24 23:05:58 UTC (rev 20627)
@@ -66,7 +66,7 @@
dtype=numpy.float64)
a = (2.0**2 + 1.2**2)**0.5
- b = (1.5**2 + 0.3**2)**0.5
+ b = (2.5**2 + 0.3**2)**0.5
c = (0.5**2 + 1.5**2)**0.5
k = 0.5 * (a + b + c)
r = (k*(k-a)*(k-b)*(k-c))**0.5 / k
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DLinear.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DLinear.py 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DLinear.py 2012-08-24 23:05:58 UTC (rev 20627)
@@ -64,9 +64,25 @@
v1 = self.vertices[1,:]
v2 = self.vertices[2,:]
v3 = self.vertices[3,:]
- vol = 1.0
- area = 0.5
- r = vol / (3*area)
+ vol = 1.0/6.0*numpy.linalg.det(numpy.array([[1.0, v0[0], v0[1], v0[2]],
+ [1.0, v1[0], v1[1], v1[2]],
+ [1.0, v2[0], v2[1], v2[2]],
+ [1.0, v3[0], v3[1], v3[2]]],
+ dtype=numpy.float64))
+ cross012 = numpy.cross(v1-v0, v2-v0)
+ area012 = 0.5*(numpy.dot(cross012, cross012))**0.5
+
+ cross013 = numpy.cross(v1-v0, v3-v0)
+ area013 = 0.5*(numpy.dot(cross013, cross013))**0.5
+
+ cross123 = numpy.cross(v2-v1, v3-v1)
+ area123 = 0.5*(numpy.dot(cross123, cross123))**0.5
+
+ cross203 = numpy.cross(v0-v2, v3-v2)
+ area203 = 0.5*(numpy.dot(cross203, cross203))**0.5
+
+ area = area012 + area013 + area123 + area203;
+ r = vol / (3.0*area)
self.minCellWidth = r
return
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DQuadratic.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DQuadratic.py 2012-08-24 19:48:13 UTC (rev 20626)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/data/Mesh3DQuadratic.py 2012-08-24 23:05:58 UTC (rev 20627)
@@ -78,9 +78,26 @@
v1 = self.vertices[1,:]
v2 = self.vertices[2,:]
v3 = self.vertices[3,:]
- vol = 1.0
- area = 0.5
- r = vol / (3*area)
+
+ vol = 1.0/6.0*numpy.linalg.det(numpy.array([[1.0, v0[0], v0[1], v0[2]],
+ [1.0, v1[0], v1[1], v1[2]],
+ [1.0, v2[0], v2[1], v2[2]],
+ [1.0, v3[0], v3[1], v3[2]]],
+ dtype=numpy.float64))
+ cross012 = numpy.cross(v1-v0, v2-v0)
+ area012 = 0.5*(numpy.dot(cross012, cross012))**0.5
+
+ cross013 = numpy.cross(v1-v0, v3-v0)
+ area013 = 0.5*(numpy.dot(cross013, cross013))**0.5
+
+ cross123 = numpy.cross(v2-v1, v3-v1)
+ area123 = 0.5*(numpy.dot(cross123, cross123))**0.5
+
+ cross203 = numpy.cross(v0-v2, v3-v2)
+ area203 = 0.5*(numpy.dot(cross203, cross203))**0.5
+
+ area = area012 + area013 + area123 + area203;
+ r = vol / (3.0*area)
self.minCellWidth = r
return
More information about the CIG-COMMITS
mailing list