[MOAB-dev] r1938 - in MOAB/trunk: . tools/mbcoupler
kraftche at mcs.anl.gov
kraftche at mcs.anl.gov
Fri Jun 27 10:05:52 CDT 2008
Author: kraftche
Date: 2008-06-27 10:05:52 -0500 (Fri, 27 Jun 2008)
New Revision: 1938
Modified:
MOAB/trunk/MBMatrix3.hpp
MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
Log:
o Terminate calculation of natural coordinates as soon as the Jacobian
is negative (much faster that waiting for iteration limit when point
in an invalid region outside of the element.)
o Check return value of nat_coords_trilinear_hex in point_in_trilinear_hex.
NOTE: point_in_trilinear_hex will return false for inverted elements
for points "inside" the element.
o Some minor cleanups to the two point_in_trilinear_hex functions.
Modified: MOAB/trunk/MBMatrix3.hpp
===================================================================
--- MOAB/trunk/MBMatrix3.hpp 2008-06-27 01:14:12 UTC (rev 1937)
+++ MOAB/trunk/MBMatrix3.hpp 2008-06-27 15:05:52 UTC (rev 1938)
@@ -146,9 +146,20 @@
inline MBMatrix3 inverse() const;
+ // invert matrix without re-calculating the
+ // reciprocal of the determinant.
+ inline MBMatrix3 inverse( double inverse_det ) const;
+
+ inline bool positive_definite() const;
+ inline bool positive_definite( double& determinant_out ) const;
+
inline MBMatrix3 transpose() const;
inline bool invert();
+
+ // Calculate determinant of 2x2 submatrix composed of the
+ // elements not in the passed row or column.
+ inline double subdet( int r, int c ) const;
};
inline MBMatrix3 operator+( const MBMatrix3& a, const MBMatrix3& b )
@@ -205,9 +216,23 @@
- d[2] * d[4] * d[6];
}
-inline MBMatrix3 MBMatrix3::inverse() const
+inline bool MBMatrix3::positive_definite( double& det ) const
{
- double i = 1.0 / determinant();
+ double subdet6 = d[1]*d[5]-d[2]*d[4];
+ double subdet7 = d[2]*d[3]-d[0]*d[5];
+ double subdet8 = d[0]*d[4]-d[1]*d[3];
+ det = d[6]*subdet6 + d[7]*subdet7 + d[8]*subdet8;
+ return d[0] > 0 && subdet8 > 0 && det > 0;
+}
+
+inline bool MBMatrix3::positive_definite() const
+{
+ double d;
+ return positive_definite(d);
+}
+
+inline MBMatrix3 MBMatrix3::inverse( double i ) const
+{
return MBMatrix3( i * (d[4] * d[8] - d[5] * d[7]),
i * (d[2] * d[7] - d[8] * d[1]),
i * (d[1] * d[5] - d[4] * d[2]),
@@ -217,6 +242,11 @@
i * (d[3] * d[7] - d[6] * d[4]),
i * (d[1] * d[6] - d[7] * d[0]),
i * (d[0] * d[4] - d[3] * d[1]) );
+}
+
+inline MBMatrix3 MBMatrix3::inverse() const
+{
+ return inverse( 1.0 / determinant() );
}
inline bool MBMatrix3::invert()
@@ -224,15 +254,7 @@
double i = 1.0 / determinant();
if (!finite(i) || fabs(i) < std::numeric_limits<double>::epsilon())
return false;
- *this= MBMatrix3( i * (d[4] * d[8] - d[5] * d[7]),
- i * (d[2] * d[7] - d[8] * d[1]),
- i * (d[1] * d[5] - d[4] * d[2]),
- i * (d[5] * d[6] - d[8] * d[3]),
- i * (d[0] * d[8] - d[6] * d[2]),
- i * (d[2] * d[3] - d[5] * d[0]),
- i * (d[3] * d[7] - d[6] * d[4]),
- i * (d[1] * d[6] - d[7] * d[0]),
- i * (d[0] * d[4] - d[3] * d[1]) );
+ *this = inverse( i );
return true;
}
@@ -243,6 +265,12 @@
d[2], d[5], d[8] );
}
+inline double MBMatrix3::subdet( int r, int c ) const
+{
+ const int r1 = (r+1)%3, r2 = (r+2)%3;
+ const int c1 = (c+1)%3, c2 = (c+2)%3;
+ return d[3*r1+c1]*d[3*r2+c2] - d[3*r1+c2]*d[3*r2+c1];
+}
MBErrorCode EigenDecomp( const MBMatrix3& a,
double Eigenvalues[3],
Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2008-06-27 01:14:12 UTC (rev 1937)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2008-06-27 15:05:52 UTC (rev 1938)
@@ -1,4 +1,5 @@
#include <iostream>
+#include <limits>
#include <assert.h>
#include "MBElemUtil.hpp"
@@ -23,16 +24,16 @@
bool VolMap::solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const
{
const double error_tol_sqr = tol*tol;
- const int max_iterations = 10000;
- int iterations = 0;
+ double det;
xi = center_xi();
MBCartVect delta = evaluate(xi) - x;
MBMatrix3 J;
- while (delta.length_squared() > error_tol_sqr) {
+ while (delta % delta > error_tol_sqr) {
J = jacobian(xi);
- if (!J.invert() || ++iterations > max_iterations)
+ det = J.determinant();
+ if (det < std::numeric_limits<double>::epsilon())
return false;
- xi -= J * delta;
+ xi -= J.inverse(1.0/det) * delta;
delta = evaluate( xi ) - x;
}
return true;
@@ -153,18 +154,11 @@
const MBCartVect& xyz,
double etol)
{
-
- const double one = 1.000001;
-
- MBCartVect nat(0.);
- nat_coords_trilinear_hex(hex, xyz, nat, etol);
-
- for (unsigned int i = 0; i < 3; i++) {
- if ((nat[i] > one) || (nat[i] < -one)) return false;
- }
-
- return true;
-
+ MBCartVect xi;
+ return nat_coords_trilinear_hex( hex, xyz, xi, etol )
+ && fabs(xi[0])-1 < etol
+ && fabs(xi[1])-1 < etol
+ && fabs(xi[2])-1 < etol;
}
@@ -174,22 +168,14 @@
const MBCartVect& box_max,
double etol)
{
-
- const double one = 1.000001;
-
- if ((xyz[0] < box_min[0]) || (xyz[0] > box_max[0])) return false;
- if ((xyz[1] < box_min[1]) || (xyz[1] > box_max[1])) return false;
- if ((xyz[2] < box_min[2]) || (xyz[2] > box_max[2])) return false;
-
- MBCartVect nat(0.);
- nat_coords_trilinear_hex(hex, xyz, nat, etol);
-
- for (unsigned int i = 0; i < 3; i++) {
- if ((nat[i] > one) || (nat[i] < -one)) return false;
- }
-
- return true;
-
+ // all values scaled by 2 (eliminates 3 flops)
+ const MBCartVect mid = box_max + box_min;
+ const MBCartVect dim = box_max - box_min;
+ const MBCartVect pt = 2*xyz - mid;
+ return fabs(pt[0]) - dim[0] < etol &&
+ fabs(pt[1]) - dim[1] < etol &&
+ fabs(pt[2]) - dim[2] < etol &&
+ point_in_trilinear_hex( hex, xyz, etol );
}
More information about the moab-dev
mailing list