[MOAB-dev] r2838 - in MOAB/trunk: . tools/iMesh tools/mbcoupler

kraftche at cae.wisc.edu kraftche at cae.wisc.edu
Fri Apr 17 14:59:41 CDT 2009


Author: kraftche
Date: 2009-04-17 14:59:41 -0500 (Fri, 17 Apr 2009)
New Revision: 2838

Added:
   MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
   MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
Removed:
   MOAB/trunk/MBElemUtil.cpp
   MOAB/trunk/MBElemUtil.hpp
Modified:
   MOAB/trunk/MBGeomUtil.cpp
   MOAB/trunk/MBGeomUtil.hpp
   MOAB/trunk/Makefile.am
   MOAB/trunk/tools/iMesh/testc_cbind.c
   MOAB/trunk/tools/mbcoupler/Makefile.am
Log:
o Revert previouis change because MBElemUtils contains a second implementation
  of point_in_trilinear_hex that depends on a bunch of C code in mbcoupler

o Copy my original point_in_trilinear_hex implementation back into MBGeomUtil



Deleted: MOAB/trunk/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/MBElemUtil.cpp	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/MBElemUtil.cpp	2009-04-17 19:59:41 UTC (rev 2838)
@@ -1,297 +0,0 @@
-#include <iostream>
-#include <limits>
-#include <assert.h>
-
-#include "MBElemUtil.hpp"
-#include "MBMatrix3.hpp"
-#include "types.h"
-
-namespace MBElemUtil {
-
-/**\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 ;
-};
-
-bool VolMap::solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const
-{
-  const double error_tol_sqr = tol*tol;
-  double det;
-  xi = center_xi();
-  MBCartVect delta = evaluate(xi) - x;
-  MBMatrix3 J;
-  while (delta % delta > error_tol_sqr) {
-    J = jacobian(xi);
-    det = J.determinant();
-    if (det < std::numeric_limits<double>::epsilon())
-      return false;
-    xi -= J.inverse(1.0/det) * delta;
-    delta = evaluate( xi ) - x;
-  }
-  return true;
-}
-
-/**\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];
-};
-
-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); }
-
-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;
-}
-
-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;
-}
-
-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 nat_coords_trilinear_hex2(const MBCartVect hex[8], 
-                               const MBCartVect& xyz,
-                               MBCartVect &ncoords,
-                               double etol)       
-
-{
-  const int ndim = 3;
-  const int nverts = 8;
-  const int vertMap[nverts] = {0,1,3,2, 4,5,7,6}; //Map from nat to lex ordering
-
-  const int n = 2; //linear
-  real coords[ndim*nverts]; //buffer
-
-  real *xm[ndim];
-  for(int i=0; i<ndim; i++)
-    xm[i] = coords + i*nverts;
-    
-  //stuff hex into coords
-  for(int i=0; i<nverts; i++){
-    real vcoord[ndim];
-    hex[i].get(vcoord);
-   
-    for(int d=0; d<ndim; d++)
-      coords[d*nverts + vertMap[i]] = vcoord[d];
-    
-  }
-
-  double dist = 0.0;
-  MBElemUtil::hex_findpt(xm, n, xyz, ncoords, dist);
-  if (3*EPS < dist) {
-      // outside element, set extremal values to something outside range
-    for (int j = 0; j < 3; j++) {
-      if (ncoords[j] < (-1.0-etol) || ncoords[j] > (1.0+etol))
-        ncoords[j] *= 10;
-    }
-  }
-  
-}
-
-bool point_in_trilinear_hex(const MBCartVect *hex, 
-                            const MBCartVect& xyz,
-                            double etol) 
-{
-  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;
-}
-
-
-bool point_in_trilinear_hex(const MBCartVect *hex, 
-                            const MBCartVect& xyz, 
-                            const MBCartVect& box_min, 
-                            const MBCartVect& box_max,
-                            double etol) 
-{
-    // 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 );
-}
-
-
-
-// Wrapper to James Lottes' findpt routines
-// hex_findpt
-// Find the parametric coordinates of an xyz point inside
-// a 3d hex spectral element with n nodes per dimension
-// xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
-// n: number of nodes per dimension -- n=2 for a linear element
-// xyz: input, point to find
-// rst: output: parametric coords of xyz inside the element. If xyz is outside the element, rst will be the coords of the closest point
-// dist: output: distance between xyz and the point with parametric coords rst
-extern "C"{
-#include "types.h"
-#include "poly.h"
-#include "tensor.h"
-#include "findpt.h"
-#include "extrafindpt.h"
-#include "errmem.h"
-}
-
-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
-  real *z[3];
-  lagrange_data ld[3];
-  opt_data_3 data;
-
-  //triplicates
-  for(int d=0; d<3; d++){
-    z[d] = tmalloc(real, n);
-    lobatto_nodes(z[d], n); 
-    lagrange_setup(&ld[d], z[d], n);
-  }
-
-  opt_alloc_3(&data, ld);
-
-  //find nearest point
-  real x_star[3];
-  xyz.get(x_star);
-
-  real r[3] = {0, 0, 0 }; // initial guess for parametric coords
-  unsigned c = opt_no_constraints_3;
-  dist = opt_findpt_3(&data, (const real **)xm, x_star, r, &c);
-  //c tells us if we landed inside the element or exactly on a face, edge, or node
-
-  //copy parametric coords back
-  rst = r;
-
-  //Clean-up (move to destructor if we decide to cache)
-  opt_free_3(&data);  
-  for(int d=0; d<3; ++d) 
-    lagrange_free(&ld[d]);
-  for(int d=0; d<3; ++d) 
-    free(z[d]);
-}
-
-
-
-
-// hex_eval
-// Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given parametric coordinates
-// field: field values for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
-// n: number of nodes per dimension -- n=2 for a linear element
-// rst: input: parametric coords of the point where we want to evaluate the field
-// value: output: value of field at rst
-
-void hex_eval(real *field,
-	      int n,
-	      MBCartVect rstCartVec,
-	      double &value)       
-{
-  int d;
-  real rst[3];
-  rstCartVec.get(rst);
-
-  //can cache stuff below
-  lagrange_data ld[3]; 
-  real *z[3];
-  for(d=0;d<3;++d){
-    z[d] = tmalloc(real, n);
-    lobatto_nodes(z[d], n);
-    lagrange_setup(&ld[d], z[d], n);
-  } 
-
-  //cut and paste -- see findpt.c
-  const unsigned 
-    nf = n*n,
-    ne = n,
-    nw = 2*n*n + 3*n;
-  real *od_work = tmalloc(real, 6*nf + 9*ne + nw);
-
-  //piece that we shouldn't want to cache
-  for(d=0; d<3; d++){
-    lagrange_0(&ld[d], rst[d]);
-  }
-  
-  value = tensor_i3(ld[0].J,ld[0].n,
-		    ld[1].J,ld[1].n,
-		    ld[2].J,ld[2].n,
-		    field,
-		    od_work);
-
-  //all this could be cached
-  for(d=0; d<3; d++){
-    free(z[d]);
-    lagrange_free(&ld[d]); 
-  }
-  free(od_work);
-}
-
-} // namespace MBElemUtil
-

Deleted: MOAB/trunk/MBElemUtil.hpp
===================================================================
--- MOAB/trunk/MBElemUtil.hpp	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/MBElemUtil.hpp	2009-04-17 19:59:41 UTC (rev 2838)
@@ -1,38 +0,0 @@
-#include "MBCore.hpp"
-#include "MBCartVect.hpp"
-
-namespace MBElemUtil {
-
-  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);
-  
-  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
-  void nat_coords_trilinear_hex2(const MBCartVect* hex_corners, 
-                                 const MBCartVect& x, 
-                                 MBCartVect& xi,
-                                 double til);
-
-
-
-  void hex_findpt(double *xm[3],
-                  int n,
-                  MBCartVect xyz, 
-                  MBCartVect& rst,
-                  double& dist);
-
-  void hex_eval(double *field,
-		int n,
-		MBCartVect rst,
-		double &value);
-} // namespace MBElemUtil

Modified: MOAB/trunk/MBGeomUtil.cpp
===================================================================
--- MOAB/trunk/MBGeomUtil.cpp	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/MBGeomUtil.cpp	2009-04-17 19:59:41 UTC (rev 2838)
@@ -21,6 +21,7 @@
 #include "MBCartVect.hpp"
 #include "MBCN.hpp"
 #include "MBGeomUtil.hpp"
+#include "MBMatrix3.hpp"
 #include <cmath>
 #include <algorithm>
 #include <assert.h>
@@ -1129,4 +1130,117 @@
   return closest % closest < tolerance * tolerance;
 }
 
+
+/**\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 ;
+};
+
+bool VolMap::solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const
+{
+  const double error_tol_sqr = tol*tol;
+  double det;
+  xi = center_xi();
+  MBCartVect delta = evaluate(xi) - x;
+  MBMatrix3 J;
+  while (delta % delta > error_tol_sqr) {
+    J = jacobian(xi);
+    det = J.determinant();
+    if (det < std::numeric_limits<double>::epsilon())
+      return false;
+    xi -= J.inverse(1.0/det) * delta;
+    delta = evaluate( xi ) - x;
+  }
+  return true;
+}
+
+/**\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];
+};
+
+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); }
+
+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;
+}
+
+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;
+}
+
+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 );
+}
+
+
+bool point_in_trilinear_hex(const MBCartVect *hex, 
+                            const MBCartVect& xyz,
+                            double etol) 
+{
+  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;
+}
+
+
+
 } // namespace MBGeoemtry

Modified: MOAB/trunk/MBGeomUtil.hpp
===================================================================
--- MOAB/trunk/MBGeomUtil.hpp	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/MBGeomUtil.hpp	2009-04-17 19:59:41 UTC (rev 2838)
@@ -254,6 +254,18 @@
                       const MBCartVect& box_center,
                       const MBCartVect& box_dims);
 
+//
+// point_in_trilinear_hex
+// Tests if a point in xyz space is within a hex element defined with
+// its eight vertex points forming a trilinear basis function.  Computes
+// the natural coordinates with respect to the hex of the xyz point 
+// and checks if each are between +/-1.  If anyone is outside the range
+// the function returns false, otherwise it returns true.
+//
+bool point_in_trilinear_hex(const MBCartVect *hex, 
+                            const MBCartVect& xyz,
+                            double etol);
+
 } // namespace MBGeoemtry
 
 #endif

Modified: MOAB/trunk/Makefile.am
===================================================================
--- MOAB/trunk/Makefile.am	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/Makefile.am	2009-04-17 19:59:41 UTC (rev 2838)
@@ -139,7 +139,6 @@
   MBCN.cpp \
   MBCNArrays.hpp \
   MBCartVect.cpp \
-  MBElemUtil.cpp \
   MBMatrix3.cpp \
   MBMatrix3.hpp \
   MBCore.cpp \
@@ -246,7 +245,6 @@
   MBCN_protos.h \
   MBCartVect.hpp \
   MBCore.hpp \
-  MBElemUtil.hpp \
   MBEntityType.h \
   MBEntityHandle.h \
   MBError.hpp \

Modified: MOAB/trunk/tools/iMesh/testc_cbind.c
===================================================================
--- MOAB/trunk/tools/iMesh/testc_cbind.c	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/tools/iMesh/testc_cbind.c	2009-04-17 19:59:41 UTC (rev 2838)
@@ -404,6 +404,65 @@
   return TRUE;
 }
 
+int adjacent_indices_test(iMesh_Instance mesh)
+{
+  int i, count[5], idxcount[5], err, idxsum = 0;
+  iBase_EntityHandle *ent_arr, *adj_arr;
+  int *off_arr, *idx_arr;
+  int ent_arr_size, adj_arr_size, off_arr_size, idx_arr_size;
+  int junk1, junk2, junk3, junk4;
+  
+  count[4] = 0;
+  for (i = 0; i < 4; ++i) {
+    iMesh_getNumOfType( mesh, 0, i, count+i, &err );
+    count[4] += count[i];
+  }
+  
+  
+  for (i = 0; i < 5; ++i) {
+    ent_arr = adj_arr = NULL;
+    off_arr = idx_arr = NULL;
+    junk1 = junk2 = junk3 = junk4 = 0;
+    iMesh_getAdjEntIndices( mesh, 0, i, iMesh_ALL_TOPOLOGIES, iBase_VERTEX,
+                            &ent_arr, &junk1, &ent_arr_size, 
+                            &adj_arr, &junk2, &adj_arr_size,
+                            &idx_arr, &junk3, &idx_arr_size,
+                            &off_arr, &junk4, &off_arr_size,
+                            &err );
+    free( ent_arr );
+    free( adj_arr );
+    free( idx_arr );
+    free( off_arr );
+    
+    if (err) {
+      printf( "%s:%d : iMesh_getAdjEntIndices failed with dimension = %d.  Error code = %d\n",
+              __FILE__, __LINE__, i, err );
+      return FALSE;
+    }
+                            
+    if (count[i] != ent_arr_size) {
+      printf( "%s:%d : Invalid number of source entities returned from iMesh_getAdjEntIndices with dimension = %d.  Expected %d, got %d\n",
+              __FILE__, __LINE__, i, count[i], ent_arr_size );
+      return FALSE;
+    }
+    
+    idxcount[i] = idx_arr_size;
+  }
+  
+  idxsum = 0;
+  for (i =0; i < 4; ++i)
+    idxsum += idxcount[i];
+  if (idxsum != idxcount[4]) {
+    printf( "%s:%d : Invalid adjacent entity count returned from iMesh_getAdjEntIndices for iBase_ALL_TYPES.  Expected %d (%d + %d + %d + %d), got %d\n",
+            __FILE__, __LINE__, idxsum, idxcount[0], idxcount[1], idxcount[2], 
+                                        idxcount[3], idxcount[4] );
+    return FALSE;
+  }
+  
+  return TRUE;
+}
+
+
 int qsort_comp_handles( const void* h1, const void* h2 )
 {
   return *(char**)h1 - *(char**)h2;
@@ -2121,6 +2180,15 @@
   number_tests++;
   printf("\n");
 
+    /* adjacent_indices_test */
+  printf("   adjacent_indices_test: ");
+  result = adjacent_indices_test(mesh);
+  handle_error_code(result, &number_tests_failed,
+                    &number_tests_not_implemented,
+                    &number_tests_successful);
+  number_tests++;
+  printf("\n");
+
     /* entity connectivity test */
   printf("   entity_connectivity_test: ");
   result = entity_connectivity_test(mesh);

Copied: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp (from rev 2835, MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp)
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp	                        (rev 0)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp	2009-04-17 19:59:41 UTC (rev 2838)
@@ -0,0 +1,297 @@
+#include <iostream>
+#include <limits>
+#include <assert.h>
+
+#include "MBElemUtil.hpp"
+#include "MBMatrix3.hpp"
+#include "types.h"
+
+namespace MBElemUtil {
+
+/**\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 ;
+};
+
+bool VolMap::solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const
+{
+  const double error_tol_sqr = tol*tol;
+  double det;
+  xi = center_xi();
+  MBCartVect delta = evaluate(xi) - x;
+  MBMatrix3 J;
+  while (delta % delta > error_tol_sqr) {
+    J = jacobian(xi);
+    det = J.determinant();
+    if (det < std::numeric_limits<double>::epsilon())
+      return false;
+    xi -= J.inverse(1.0/det) * delta;
+    delta = evaluate( xi ) - x;
+  }
+  return true;
+}
+
+/**\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];
+};
+
+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); }
+
+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;
+}
+
+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;
+}
+
+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 nat_coords_trilinear_hex2(const MBCartVect hex[8], 
+                               const MBCartVect& xyz,
+                               MBCartVect &ncoords,
+                               double etol)       
+
+{
+  const int ndim = 3;
+  const int nverts = 8;
+  const int vertMap[nverts] = {0,1,3,2, 4,5,7,6}; //Map from nat to lex ordering
+
+  const int n = 2; //linear
+  real coords[ndim*nverts]; //buffer
+
+  real *xm[ndim];
+  for(int i=0; i<ndim; i++)
+    xm[i] = coords + i*nverts;
+    
+  //stuff hex into coords
+  for(int i=0; i<nverts; i++){
+    real vcoord[ndim];
+    hex[i].get(vcoord);
+   
+    for(int d=0; d<ndim; d++)
+      coords[d*nverts + vertMap[i]] = vcoord[d];
+    
+  }
+
+  double dist = 0.0;
+  MBElemUtil::hex_findpt(xm, n, xyz, ncoords, dist);
+  if (3*EPS < dist) {
+      // outside element, set extremal values to something outside range
+    for (int j = 0; j < 3; j++) {
+      if (ncoords[j] < (-1.0-etol) || ncoords[j] > (1.0+etol))
+        ncoords[j] *= 10;
+    }
+  }
+  
+}
+
+bool point_in_trilinear_hex(const MBCartVect *hex, 
+                            const MBCartVect& xyz,
+                            double etol) 
+{
+  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;
+}
+
+
+bool point_in_trilinear_hex(const MBCartVect *hex, 
+                            const MBCartVect& xyz, 
+                            const MBCartVect& box_min, 
+                            const MBCartVect& box_max,
+                            double etol) 
+{
+    // 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 );
+}
+
+
+
+// Wrapper to James Lottes' findpt routines
+// hex_findpt
+// Find the parametric coordinates of an xyz point inside
+// a 3d hex spectral element with n nodes per dimension
+// xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
+// n: number of nodes per dimension -- n=2 for a linear element
+// xyz: input, point to find
+// rst: output: parametric coords of xyz inside the element. If xyz is outside the element, rst will be the coords of the closest point
+// dist: output: distance between xyz and the point with parametric coords rst
+extern "C"{
+#include "types.h"
+#include "poly.h"
+#include "tensor.h"
+#include "findpt.h"
+#include "extrafindpt.h"
+#include "errmem.h"
+}
+
+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
+  real *z[3];
+  lagrange_data ld[3];
+  opt_data_3 data;
+
+  //triplicates
+  for(int d=0; d<3; d++){
+    z[d] = tmalloc(real, n);
+    lobatto_nodes(z[d], n); 
+    lagrange_setup(&ld[d], z[d], n);
+  }
+
+  opt_alloc_3(&data, ld);
+
+  //find nearest point
+  real x_star[3];
+  xyz.get(x_star);
+
+  real r[3] = {0, 0, 0 }; // initial guess for parametric coords
+  unsigned c = opt_no_constraints_3;
+  dist = opt_findpt_3(&data, (const real **)xm, x_star, r, &c);
+  //c tells us if we landed inside the element or exactly on a face, edge, or node
+
+  //copy parametric coords back
+  rst = r;
+
+  //Clean-up (move to destructor if we decide to cache)
+  opt_free_3(&data);  
+  for(int d=0; d<3; ++d) 
+    lagrange_free(&ld[d]);
+  for(int d=0; d<3; ++d) 
+    free(z[d]);
+}
+
+
+
+
+// hex_eval
+// Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given parametric coordinates
+// field: field values for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
+// n: number of nodes per dimension -- n=2 for a linear element
+// rst: input: parametric coords of the point where we want to evaluate the field
+// value: output: value of field at rst
+
+void hex_eval(real *field,
+	      int n,
+	      MBCartVect rstCartVec,
+	      double &value)       
+{
+  int d;
+  real rst[3];
+  rstCartVec.get(rst);
+
+  //can cache stuff below
+  lagrange_data ld[3]; 
+  real *z[3];
+  for(d=0;d<3;++d){
+    z[d] = tmalloc(real, n);
+    lobatto_nodes(z[d], n);
+    lagrange_setup(&ld[d], z[d], n);
+  } 
+
+  //cut and paste -- see findpt.c
+  const unsigned 
+    nf = n*n,
+    ne = n,
+    nw = 2*n*n + 3*n;
+  real *od_work = tmalloc(real, 6*nf + 9*ne + nw);
+
+  //piece that we shouldn't want to cache
+  for(d=0; d<3; d++){
+    lagrange_0(&ld[d], rst[d]);
+  }
+  
+  value = tensor_i3(ld[0].J,ld[0].n,
+		    ld[1].J,ld[1].n,
+		    ld[2].J,ld[2].n,
+		    field,
+		    od_work);
+
+  //all this could be cached
+  for(d=0; d<3; d++){
+    free(z[d]);
+    lagrange_free(&ld[d]); 
+  }
+  free(od_work);
+}
+
+} // namespace MBElemUtil
+


Property changes on: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Copied: MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp (from rev 2835, MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp)
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp	                        (rev 0)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp	2009-04-17 19:59:41 UTC (rev 2838)
@@ -0,0 +1,38 @@
+#include "MBCore.hpp"
+#include "MBCartVect.hpp"
+
+namespace MBElemUtil {
+
+  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);
+  
+  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
+  void nat_coords_trilinear_hex2(const MBCartVect* hex_corners, 
+                                 const MBCartVect& x, 
+                                 MBCartVect& xi,
+                                 double til);
+
+
+
+  void hex_findpt(double *xm[3],
+                  int n,
+                  MBCartVect xyz, 
+                  MBCartVect& rst,
+                  double& dist);
+
+  void hex_eval(double *field,
+		int n,
+		MBCartVect rst,
+		double &value);
+} // namespace MBElemUtil


Property changes on: MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Modified: MOAB/trunk/tools/mbcoupler/Makefile.am
===================================================================
--- MOAB/trunk/tools/mbcoupler/Makefile.am	2009-04-17 19:48:52 UTC (rev 2837)
+++ MOAB/trunk/tools/mbcoupler/Makefile.am	2009-04-17 19:59:41 UTC (rev 2838)
@@ -14,6 +14,7 @@
 endif
 
 libmbcoupler_la_SOURCES = \
+   MBElemUtil.cpp \
    poly.c   \
    poly.h   \
    tensor.c \
@@ -28,6 +29,7 @@
    MBCoupler.cpp
 
 libmbcoupler_la_include_HEADERS = \
+   MBElemUtil.hpp \
    MBCoupler.hpp
 
 libmbcoupler_la_includedir = $(includedir)



More information about the moab-dev mailing list