[MOAB-dev] r1936 - in MOAB/trunk: . tools/mbcoupler
kraftche at mcs.anl.gov
kraftche at mcs.anl.gov
Thu Jun 26 19:09:29 CDT 2008
Author: kraftche
Date: 2008-06-26 19:09:29 -0500 (Thu, 26 Jun 2008)
New Revision: 1936
Added:
MOAB/trunk/tools/mbcoupler/ElemUtilTest.cpp
Modified:
MOAB/trunk/MBCartVect.hpp
MOAB/trunk/MBMatrix3.hpp
MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
MOAB/trunk/tools/mbcoupler/Makefile.am
Log:
o Change MBElemUtil from a class full of static functions (the old way)
to a namespace (the new way.)
o Add unit tests for nat_coords_trilinear_hex
o Re-implement nat_coords_trilinear_hex correctly. New implementation
uses simple Newton's method solver with a generic mapping function
class. Extensible to any volume element type by providing a mapping
function class.
Modified: MOAB/trunk/MBCartVect.hpp
===================================================================
--- MOAB/trunk/MBCartVect.hpp 2008-06-26 20:09:48 UTC (rev 1935)
+++ MOAB/trunk/MBCartVect.hpp 2008-06-27 00:09:29 UTC (rev 1936)
@@ -69,6 +69,8 @@
inline double length() const; //!< vector length
+ inline double length_squared() const;
+
inline void normalize(); //!< make unit length
inline void flip(); //!< flip direction
@@ -112,6 +114,9 @@
inline double MBCartVect::length() const
{ return std::sqrt( *this % *this ); }
+inline double MBCartVect::length_squared() const
+ { return d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; }
+
inline void MBCartVect::normalize()
{ *this /= length(); }
Modified: MOAB/trunk/MBMatrix3.hpp
===================================================================
--- MOAB/trunk/MBMatrix3.hpp 2008-06-26 20:09:48 UTC (rev 1935)
+++ MOAB/trunk/MBMatrix3.hpp 2008-06-27 00:09:29 UTC (rev 1936)
@@ -147,6 +147,7 @@
inline MBMatrix3 transpose() const;
+ inline bool invert();
};
inline MBMatrix3 operator+( const MBMatrix3& a, const MBMatrix3& b )
@@ -217,6 +218,23 @@
i * (d[0] * d[4] - d[3] * d[1]) );
}
+inline bool MBMatrix3::invert()
+{
+ 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]) );
+ return true;
+}
+
inline MBMatrix3 MBMatrix3::transpose() const
{
return MBMatrix3( d[0], d[3], d[6],
Added: MOAB/trunk/tools/mbcoupler/ElemUtilTest.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/ElemUtilTest.cpp (rev 0)
+++ MOAB/trunk/tools/mbcoupler/ElemUtilTest.cpp 2008-06-27 00:09:29 UTC (rev 1936)
@@ -0,0 +1,83 @@
+#include "TestUtil.hpp"
+#include "MBElemUtil.hpp"
+#include <iostream>
+
+void test_hex_nat_coords();
+
+int main()
+{
+ int rval = 0;
+ rval += RUN_TEST(test_hex_nat_coords);
+ return rval;
+}
+
+const MBCartVect cube_corners[8] = { MBCartVect( 0, 0, 0 ),
+ MBCartVect( 1, 0, 0 ),
+ MBCartVect( 1, 1, 0 ),
+ MBCartVect( 0, 1, 0 ),
+ MBCartVect( 0, 0, 1 ),
+ MBCartVect( 1, 0, 1 ),
+ MBCartVect( 1, 1, 1 ),
+ MBCartVect( 0, 1, 1 ) };
+
+
+const MBCartVect hex_corners[8] = { MBCartVect( 1.0, 0.0, 0.0 ),
+ MBCartVect( 1.0, 1.0, 0.3 ),
+ MBCartVect( 0.0, 2.0, 0.6 ),
+ MBCartVect( 0.2, 1.1, 0.4 ),
+ MBCartVect( 1.5, 0.3, 1.0 ),
+ MBCartVect( 1.5, 1.3, 1.0 ),
+ MBCartVect( 0.5, 2.3, 1.0 ),
+ MBCartVect( 0.7, 1.4, 1.0 ) };
+
+/** shape function for trilinear hex */
+MBCartVect hex_map( const MBCartVect& xi, const MBCartVect* corners )
+{
+ MBCartVect x(0.0);
+ x += (1 - xi[0]) * (1 - xi[1]) * (1 - xi[2]) * corners[0];
+ x += (1 + xi[0]) * (1 - xi[1]) * (1 - xi[2]) * corners[1];
+ x += (1 + xi[0]) * (1 + xi[1]) * (1 - xi[2]) * corners[2];
+ x += (1 - xi[0]) * (1 + xi[1]) * (1 - xi[2]) * corners[3];
+ x += (1 - xi[0]) * (1 - xi[1]) * (1 + xi[2]) * corners[4];
+ x += (1 + xi[0]) * (1 - xi[1]) * (1 + xi[2]) * corners[5];
+ x += (1 + xi[0]) * (1 + xi[1]) * (1 + xi[2]) * corners[6];
+ x += (1 - xi[0]) * (1 + xi[1]) * (1 + xi[2]) * corners[7];
+ return x *= 0.125;
+}
+
+
+void test_hex_nat_coords()
+{
+ MBCartVect xi, result_xi;
+ bool valid;
+ const double EPS = 1e-6;
+
+ // first test with cube because it's easier to debug failures
+ for (xi[0] = -1; xi[0] <= 1; xi[0] += 0.2) {
+ for (xi[1] = -1; xi[1] <= 1; xi[1] += 0.2) {
+ for (xi[2] = -1; xi[2] <= 1; xi[2] += 0.2) {
+ const MBCartVect pt = hex_map(xi, cube_corners);
+ valid = MBElemUtil::nat_coords_trilinear_hex( cube_corners, pt, result_xi, EPS/10 );
+ CHECK(valid);
+ CHECK_REAL_EQUAL( xi[0], result_xi[0], EPS );
+ CHECK_REAL_EQUAL( xi[1], result_xi[1], EPS );
+ CHECK_REAL_EQUAL( xi[2], result_xi[2], EPS );
+ }
+ }
+ }
+
+ // now test with distorted hex
+ for (xi[0] = -1; xi[0] <= 1; xi[0] += 0.2) {
+ for (xi[1] = -1; xi[1] <= 1; xi[1] += 0.2) {
+ for (xi[2] = -1; xi[2] <= 1; xi[2] += 0.2) {
+ const MBCartVect pt = hex_map(xi, hex_corners);
+ valid = MBElemUtil::nat_coords_trilinear_hex( hex_corners, pt, result_xi, EPS/10 );
+ CHECK(valid);
+ CHECK_REAL_EQUAL( xi[0], result_xi[0], EPS );
+ CHECK_REAL_EQUAL( xi[1], result_xi[1], EPS );
+ CHECK_REAL_EQUAL( xi[2], result_xi[2], EPS );
+ }
+ }
+ }
+}
+
Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2008-06-26 20:09:48 UTC (rev 1935)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2008-06-27 00:09:29 UTC (rev 1936)
@@ -2,153 +2,118 @@
#include <assert.h>
#include "MBElemUtil.hpp"
+#include "MBMatrix3.hpp"
#include "types.h"
-//
-// nat_coords_trilinear_hex
-// Given an MBEntityHandle defining a MBHEX element defined by a
-// trilinear basis and a set of xyz points in space, it finds the
-// cooresponding natural coordinates and returns them to the pointer
-// nat. The method uses an iterative technique, guessing initial
-// coordinates {0, 0, 0} and calculating new coordinates based on the
-// ones just calculated.
-//
-void MBElemUtil::nat_coords_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- MBCartVect &ncoords,
- double etol)
-{
+namespace MBElemUtil {
- // short-circuit for now...
- ncoords[0] = ncoords[1] = ncoords[2] = 0.5;
- return;
-
- MBCartVect nat(0.);
- double A[8]; double B[8]; double C[8]; double D[8];
- double Ax, By, Cz;
- double err = 1.e37;
+/**\brief Class representing a 3-D mapping function (e.g. shape function for volume element) */
+class VolMap {
+ public:
+ /**\brief Return $\vec \xi$ corresponding to logical center of element */
+ virtual MBCartVect center_xi() const = 0;
+ /**\brief Evaluate mapping function (calculate $\vec x = F($\vec \xi)$ )*/
+ virtual MBCartVect evaluate( const MBCartVect& xi ) const = 0;
+ /**\brief Evaluate Jacobian of mapping function */
+ virtual MBMatrix3 jacobian( const MBCartVect& xi ) const = 0;
+ /**\brief Evaluate inverse of mapping function (calculate $\vec \xi = F^-1($\vec x)$ )*/
+ bool solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const ;
+};
- double xt, yt, zt;
- double nxi, neta, nmu;
- double pxi, peta, pmu;
- double tmp;
-
- // Iterative estimate of natural coordinates
- while ( err > etol ) {
+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;
+ xi = center_xi();
+ MBCartVect delta = evaluate(xi) - x;
+ MBMatrix3 J;
+ while (delta.length_squared() > error_tol_sqr) {
+ J = jacobian(xi);
+ if (!J.invert() || ++iterations > max_iterations)
+ return false;
+ xi -= J * delta;
+ delta = evaluate( xi ) - x;
+ }
+ return true;
+}
- // Estimate the xi-coordinate
- A[0] = (1 - nat[1]) * (1 - nat[2]);
- A[1] = A[0];
- A[2] = (1 + nat[1]) * (1 - nat[2]);
- A[3] = A[2];
- A[4] = (1 - nat[1]) * (1 + nat[2]);
- A[5] = A[4];
- A[6] = (1 + nat[1]) * (1 + nat[2]);
- A[7] = A[6];
+/**\brief Shape function for trilinear hexahedron */
+class LinearHexMap : public VolMap {
+ public:
+ LinearHexMap( const MBCartVect* corner_coords ) : corners(corner_coords) {}
+ virtual MBCartVect center_xi() const;
+ virtual MBCartVect evaluate( const MBCartVect& xi ) const;
+ virtual MBMatrix3 jacobian( const MBCartVect& xi ) const;
+ private:
+ const MBCartVect* corners;
+ static const double corner_xi[8][3];
+};
- Ax = 0;
- for (unsigned int i = 0; i < 8; i++) {
- Ax = Ax + A[i]*hex[i][0];
- }
- double tmp_d = -A[0]*hex[0][0] + A[1]*hex[1][0] + A[2]*hex[2][0] - A[3]*hex[3][0]
- -A[4]*hex[4][0] + A[5]*hex[5][0] + A[6]*hex[6][0] - A[7]*hex[7][0];
- if (0.0 == tmp_d) nat[0] = 2.0;
- else nat[0] = (8*xyz[0] - Ax ) / tmp_d;
+const double LinearHexMap::corner_xi[8][3] = { { -1, -1, -1 },
+ { 1, -1, -1 },
+ { 1, 1, -1 },
+ { -1, 1, -1 },
+ { -1, -1, 1 },
+ { 1, -1, 1 },
+ { 1, 1, 1 },
+ { -1, 1, 1 } };
+MBCartVect LinearHexMap::center_xi() const
+ { return MBCartVect(0.0); }
- // Estimate the eta-coordinate
- B[0] = (1 - nat[0]) * (1 - nat[2]);
- B[1] = (1 + nat[0]) * (1 - nat[2]);
- B[2] = B[1];
- B[3] = B[0];
- B[4] = (1 - nat[0]) * (1 + nat[2]);
- B[5] = (1 + nat[0]) * (1 + nat[2]);
- B[6] = B[5];
- B[7] = B[4];
+MBCartVect LinearHexMap::evaluate( const MBCartVect& xi ) const
+{
+ MBCartVect x(0.0);
+ for (unsigned i = 0; i < 8; ++i) {
+ const double N_i = (1 + xi[0]*corner_xi[i][0])
+ * (1 + xi[1]*corner_xi[i][1])
+ * (1 + xi[2]*corner_xi[i][2]);
+ x += N_i * corners[i];
+ }
+ x *= 0.125;
+ return x;
+}
- By = 0;
- for (unsigned int i = 0; i < 8; i++) {
- By = By + B[i]*hex[i][1];
- }
- tmp_d = -B[0]*hex[0][1] - B[1]*hex[1][1] + B[2]*hex[2][1] + B[3]*hex[3][1]
- -B[4]*hex[4][1] - B[5]*hex[5][1] + B[6]*hex[6][1] + B[7]*hex[7][1];
-
- if (0.0 == tmp_d) nat[1] = 2.0;
- else nat[1] = (8*xyz[1] - By ) / tmp_d;
+MBMatrix3 LinearHexMap::jacobian( const MBCartVect& xi ) const
+{
+ MBMatrix3 J(0.0);
+ for (unsigned i = 0; i < 8; ++i) {
+ const double xi_p = 1 + xi[0]*corner_xi[i][0];
+ const double eta_p = 1 + xi[1]*corner_xi[i][1];
+ const double zeta_p = 1 + xi[2]*corner_xi[i][2];
+ const double dNi_dxi = corner_xi[i][0] * eta_p * zeta_p;
+ const double dNi_deta = corner_xi[i][1] * xi_p * zeta_p;
+ const double dNi_dzeta = corner_xi[i][2] * xi_p * eta_p;
+ J(0,0) += dNi_dxi * corners[i][0];
+ J(1,0) += dNi_dxi * corners[i][1];
+ J(2,0) += dNi_dxi * corners[i][2];
+ J(0,1) += dNi_deta * corners[i][0];
+ J(1,1) += dNi_deta * corners[i][1];
+ J(2,1) += dNi_deta * corners[i][2];
+ J(0,2) += dNi_dzeta * corners[i][0];
+ J(1,2) += dNi_dzeta * corners[i][1];
+ J(2,2) += dNi_dzeta * corners[i][2];
+ }
+ return J *= 0.125;
+}
- // Estimate the mu-coordinate
- C[0] = (1 - nat[0]) * (1 - nat[1]);
- C[1] = (1 + nat[0]) * (1 - nat[1]);
- C[2] = (1 + nat[0]) * (1 + nat[1]);
- C[3] = (1 - nat[0]) * (1 + nat[1]);
- C[4] = C[0];
- C[5] = C[1];
- C[6] = C[2];
- C[7] = C[3];
-
- Cz = 0;
- for (unsigned int i = 0; i < 8; i++) {
- Cz = Cz + C[i]*hex[i][2];
- }
- tmp_d = -C[0]*hex[0][2] - C[1]*hex[1][2] - C[2]*hex[2][2] - C[3]*hex[3][2]
- +C[4]*hex[4][2] + C[5]*hex[5][2] + C[6]*hex[6][2] + C[7]*hex[7][2];
-
- if (0.0 == tmp_d) nat[2] = 2.0;
- else nat[2] = (8*xyz[2] - Cz ) / tmp_d;
-
- // Shortcut variables...
- nxi = 1 - nat[0];
- neta = 1 - nat[1];
- nmu = 1 - nat[2];
- pxi = 1 + nat[0];
- peta = 1 + nat[1];
- pmu = 1 + nat[2];
- D[0] = nxi * neta * nmu;
- D[1] = pxi * neta * nmu;
- D[2] = pxi * peta * nmu;
- D[3] = nxi * peta * nmu;
- D[4] = nxi * neta * pmu;
- D[5] = pxi * neta * pmu;
- D[6] = pxi * peta * pmu;
- D[7] = nxi * peta * pmu;
-
- // Compute corresponding estimates for x, y, and z to check
- xt = 0.125 * ( D[0] * hex[0][0] + D[1] * hex[1][0] +
- D[2] * hex[2][0] + D[3] * hex[3][0] +
- D[4] * hex[4][0] + D[5] * hex[5][0] +
- D[6] * hex[6][0] + D[7] * hex[7][0] );
-
- yt = 0.125 * ( D[0] * hex[0][1] + D[1] * hex[1][1] +
- D[2] * hex[2][1] + D[3] * hex[3][1] +
- D[4] * hex[4][1] + D[5] * hex[5][1] +
- D[6] * hex[6][1] + D[7] * hex[7][1] );
-
- zt = 0.125 * ( D[0] * hex[0][2] + D[1] * hex[1][2] +
- D[2] * hex[2][2] + D[3] * hex[3][2] +
- D[4] * hex[4][2] + D[5] * hex[5][2] +
- D[6] * hex[6][2] + D[7] * hex[7][2] );
-
- // Compute error
- err = fabs(xt - xyz[0]);
- tmp = fabs(yt - xyz[1]);
- if (tmp > err) err = tmp;
- tmp = fabs(zt - xyz[2]);
- if (tmp > err) err = tmp;
-
- }
-
- ncoords = nat;
- return;
-
+bool nat_coords_trilinear_hex( const MBCartVect* corner_coords,
+ const MBCartVect& x,
+ MBCartVect& xi,
+ double tol )
+{
+ return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
}
+
//
// nat_coords_trilinear_hex2
// Duplicate functionality of nat_coords_trilinear_hex using hex_findpt
//
-void MBElemUtil::nat_coords_trilinear_hex2(MBCartVect hex[8],
- MBCartVect xyz,
- MBCartVect &ncoords,
- double etol)
+void nat_coords_trilinear_hex2(const MBCartVect hex[8],
+ const MBCartVect& xyz,
+ MBCartVect &ncoords,
+ double etol)
{
const int ndim = 3;
@@ -184,9 +149,9 @@
}
-bool MBElemUtil::point_in_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- double etol)
+bool point_in_trilinear_hex(const MBCartVect *hex,
+ const MBCartVect& xyz,
+ double etol)
{
const double one = 1.000001;
@@ -203,10 +168,11 @@
}
-bool MBElemUtil::point_in_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- MBCartVect box_min, MBCartVect box_max,
- double etol)
+bool point_in_trilinear_hex(const MBCartVect *hex,
+ const MBCartVect& xyz,
+ const MBCartVect& box_min,
+ const MBCartVect& box_max,
+ double etol)
{
const double one = 1.000001;
@@ -247,11 +213,11 @@
#include "errmem.h"
}
-void MBElemUtil::hex_findpt(real *xm[3],
- int n,
- MBCartVect xyz,
- MBCartVect &rst,
- double &dist)
+void hex_findpt(real *xm[3],
+ int n,
+ MBCartVect xyz,
+ MBCartVect &rst,
+ double &dist)
{
//compute stuff that only depends on the order -- could be cached
@@ -287,3 +253,5 @@
for(int d=0; d<3; ++d)
free(z[d]);
}
+
+} // namespace MBElemUtil
Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp 2008-06-26 20:09:48 UTC (rev 1935)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp 2008-06-27 00:09:29 UTC (rev 1936)
@@ -1,35 +1,33 @@
#include "MBCore.hpp"
#include "MBCartVect.hpp"
-class MBElemUtil {
+namespace MBElemUtil {
-public:
-
- static void nat_coords_trilinear_hex(MBCartVect*,
- MBCartVect,
- MBCartVect&,
- double);
- static bool point_in_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- double etol);
+ bool nat_coords_trilinear_hex(const MBCartVect* hex_corners,
+ const MBCartVect& x,
+ MBCartVect& xi,
+ double tol);
+ bool point_in_trilinear_hex(const MBCartVect *hex_corners,
+ const MBCartVect& xyz,
+ double etol);
- static bool point_in_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- MBCartVect box_min,
- MBCartVect box_max,
- double etol);
+ bool point_in_trilinear_hex(const MBCartVect *hex_corners,
+ const MBCartVect& xyz,
+ const MBCartVect& box_min,
+ const MBCartVect& box_max,
+ double etol);
//wrapper to hex_findpt
- static void nat_coords_trilinear_hex2(MBCartVect*,
- MBCartVect,
- MBCartVect&,
- double);
+ void nat_coords_trilinear_hex2(const MBCartVect* hex_corners,
+ const MBCartVect& x,
+ MBCartVect& xi,
+ double til);
- static void hex_findpt(double *xm[3],
- int n,
- MBCartVect xyz,
- MBCartVect& rst,
- double& dist);
-};
+ void hex_findpt(double *xm[3],
+ int n,
+ MBCartVect xyz,
+ MBCartVect& rst,
+ double& dist);
+} // namespace MBElemUtil
Modified: MOAB/trunk/tools/mbcoupler/Makefile.am
===================================================================
--- MOAB/trunk/tools/mbcoupler/Makefile.am 2008-06-26 20:09:48 UTC (rev 1935)
+++ MOAB/trunk/tools/mbcoupler/Makefile.am 2008-06-27 00:09:29 UTC (rev 1936)
@@ -34,10 +34,12 @@
# check that only libraries are going in $(libdir)
cfgdir = $(libdir)
-TESTS = findpt_test
+TESTS = findpt_test elem_util_test
check_PROGRAMS = $(TESTS)
findpt_test_SOURCES = findpt_test.cpp
findpt_test_LDADD = libmbcoupler.la $(top_builddir)/libMOAB.la
+elem_util_test_SOURCES = ElemUtilTest.cpp
+elem_util_test_LDADD = libmbcoupler.la $(top_builddir)/libMOAB.la
if PARALLEL
TESTS += mbcoupler_test
More information about the moab-dev
mailing list