Merged in vijaysm/add-eigen-support (pull request #248)

Add support for Eigen library for Dense matrix/vector algebra in MOAB

src/GeomUtil.cpp

@@ -1449,7 +1449,7 @@ bool VolMap::solve_inverse( const CartVect& x, CartVect& xi, double tol ) const
     det = J.determinant();
     if (det < std::numeric_limits<double>::epsilon())
       return false;
-    xi -= J.inverse(1.0/det) * delta;
+    xi -= J.inverse() * delta;
     delta = evaluate( xi ) - x;
   }
   return true;

src/LocalDiscretization/ElemEvaluator.cpp

@@ -61,7 +61,7 @@ namespace moab {
         }
 
           // new params tries to eliminate residual
-        *cvparams -= J.inverse(1.0/det) * res;
+        *cvparams -= J.inverse() * res;
 
           // get the new forward-evaluated position, and its difference from the target pt
         rval = (*eval)(params, verts, ndim, 

src/LocalDiscretization/LinearTet.cpp

@@ -35,7 +35,7 @@ namespace moab
     }
 
     ErrorCode LinearTet::evalFcn(const double *params, const double *field, const int /*ndim*/, const int num_tuples, 
-                                 double */*work*/, double *result) {
+                                 double* /*work*/, double *result) {
       assert(params && field && num_tuples > 0);
       std::vector<double> f0(num_tuples);
       std::copy(field, field+num_tuples, f0.begin());
@@ -50,7 +50,7 @@ namespace moab
       return MB_SUCCESS;
     }
 
-    ErrorCode LinearTet::integrateFcn(const double *field, const double */*verts*/, const int nverts, const int /*ndim*/, const int num_tuples,
+    ErrorCode LinearTet::integrateFcn(const double *field, const double* /*verts*/, const int nverts, const int /*ndim*/, const int num_tuples,
                                       double *work, double *result) 
     {
       assert(field && num_tuples > 0);
@@ -122,9 +122,11 @@ namespace moab
       CartVect res = new_pos - *cvposn;
       Matrix3 J;
       rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
+#ifndef NDEBUG
       double det = J.determinant();
       assert(det > std::numeric_limits<double>::epsilon());
-      Matrix3 Ji = J.inverse(1.0/det);
+#endif
+      Matrix3 Ji = J.inverse();
 
       int iters=0;
         // while |res| larger than tol

src/LocalDiscretization/LinearTri.cpp

@@ -119,9 +119,11 @@ namespace moab
       CartVect res = new_pos - *cvposn;
       Matrix3 J;
       rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
+#ifndef NDEBUG
       double det = J.determinant();
       assert(det > std::numeric_limits<double>::epsilon());
-      Matrix3 Ji = J.inverse(1.0/det);
+#endif
+      Matrix3 Ji = J.inverse();
 
       int iters=0;
         // while |res| larger than tol

src/Makefile.am

@@ -14,6 +14,8 @@ if PARALLEL
   AM_CPPFLAGS += -I$(srcdir)/parallel
 endif
 
+include moab/EigenHeaders
+
 SUBDIRS += io LocalDiscretization verdict RefineMesh .
 libMOAB_la_LIBADD += io/libmoabio.la LocalDiscretization/libLocalDiscretization.la verdict/libmoabverdict.la RefineMesh/libRefineMesh.la
 
@@ -219,6 +221,7 @@ nobase_libMOAB_la_include_HEADERS = \
   moab/Util.hpp \
   moab/WriteUtilIface.hpp \
   moab/WriterIface.hpp \
+  $(EIGEN_INST_HDRS) \
   MBEntityType.h \
   MBCN.h \
   MBCN_protos.h \

src/OrientedBox.cpp

@@ -35,8 +35,6 @@
 #include "moab/CN.hpp"
 #include "moab/OrientedBox.hpp"
 #include "moab/Range.hpp"
-#include "moab/Matrix3.hpp"
-#include "moab/Util.hpp"
 #include <ostream>
 #include <assert.h>
 #include <limits>
@@ -47,17 +45,17 @@ std::ostream& operator<<( std::ostream& s, const OrientedBox& b )
 {
   return s << b.center 
            << " + " 
-           << b.axis[0] 
+           << b.axes.col(0) 
 #if MB_ORIENTED_BOX_UNIT_VECTORS
            << ":" << b.length[0] 
 #endif
            << " x " 
-           << b.axis[1] 
+           << b.axes.col(1) 
 #if MB_ORIENTED_BOX_UNIT_VECTORS
            << ":" << b.length[1] 
 #endif
            << " x " 
-           << b.axis[2]
+           << b.axes.col(2)
 #if MB_ORIENTED_BOX_UNIT_VECTORS
            << ":" << b.length[2] 
 #endif
@@ -86,47 +84,59 @@ static double point_perp( const CartVect& p,   // closest to this point
 #endif
   return Util::is_finite(t) ? t : 0.0;
 }
-
-OrientedBox::OrientedBox( const CartVect axes[3], const CartVect& mid )
-  : center(mid)
-{
-    // re-order axes by length
-  CartVect len( axes[0].length(), axes[1].length(), axes[2].length() );
-  axis[0] = axes[0];
-  axis[1] = axes[1];
-  axis[2] = axes[2];
   
+void OrientedBox::order_axes_by_length(double ax1_len, double ax2_len, double ax3_len)
+{
+
+  CartVect len( ax1_len, ax2_len, ax3_len );
+ 
   if (len[2] < len[1])
-  {
-    if (len[2] < len[0]) {
-      std::swap( len[0], len[2] );
-      std::swap( axis[0], axis[2] );
+    {
+      if (len[2] < len[0]) {
+	std::swap( len[0], len[2] );
+	axes.swapcol( 0, 2 );
+      }
     }
-  }
   else if (len[1] < len[0]) {
     std::swap( len[0], len[1] );
-    std::swap( axis[0], axis[1] );
+    axes.swapcol( 0, 1 );
   }
   if (len[1] > len[2]) {
     std::swap( len[1], len[2] );
-    std::swap( axis[1], axis[2] );
+    axes.swapcol( 1, 2 );
   }
   
 #if MB_ORIENTED_BOX_UNIT_VECTORS
-  this->length = len;
-  if (len[0] > 0.0)
-    axis[0] /= len[0];
-  if (len[1] > 0.0)
-    axis[1] /= len[1];
-  if (len[2] > 0.0)
-    axis[2] /= len[2];
+    length = len;
+    if (len[0] > 0.0)
+      axes.colscale(0, 1.0/len[0]);
+    if (len[1] > 0.0)
+      axes.colscale(1, 1.0/len[1]);
+    if (len[2] > 0.0)
+      axes.colscale(2 ,1.0/len[2]);
 #endif
 
 #if MB_ORIENTED_BOX_OUTER_RADIUS
-  radius = len.length();
+    radius = len.length();
 #endif
 }
 
+OrientedBox::OrientedBox( const CartVect axes_in[3], const CartVect& mid )
+  : center(mid)
+{
+
+  axes = Matrix3( axes_in[0], axes_in[1], axes_in[2], false );
+
+  order_axes_by_length( axes_in[0].length(),axes_in[1].length(),axes_in[2].length() );
+  
+}
+
+OrientedBox::OrientedBox( const Matrix3& axes_mat, const CartVect& mid )
+  : center(mid), axes(axes_mat)
+{
+  order_axes_by_length( axes.col(0).length(), axes.col(1).length(), axes.col(2).length() );
+}
+
 ErrorCode OrientedBox::tag_handle( Tag& handle_out,
                                        Interface* instance,
                                        const char* name)
@@ -174,13 +184,11 @@ static ErrorCode box_from_axes( OrientedBox& result,
   for (Range::iterator i = points.begin(); i != points.end(); ++i)
   {
     CartVect coords;
-    rval = instance->get_coords( &*i, 1, coords.array() );
-    if (MB_SUCCESS != rval)
-      return rval;
-    
+    rval = instance->get_coords( &*i, 1, coords.array() );MB_CHK_ERR(rval);
+
     for (int d = 0; d < 3; ++d)
     {
-      double t = point_perp( coords, result.center, result.axis[d] );
+      const double t = point_perp( coords, result.center, result.axes.col(d) );
       if (t < min[d])
         min[d] = t;
       if (t > max[d])
@@ -190,14 +198,14 @@ static ErrorCode box_from_axes( OrientedBox& result,
   
     // We now have a box defined by three orthogonal line segments
     // that intersect at the center of the box.  Each line segment
-    // is defined as result.center + t * result.axis[i], where the
+    // is defined as result.center + t * result.axes[i], where the
     // range of t is [min[i], max[i]].
   
     // Calculate new center
-  CartVect mid = 0.5 * (min + max);
-  result.center += mid[0] * result.axis[0] +
-                   mid[1] * result.axis[1] +
-                   mid[2] * result.axis[2];
+  const CartVect mid = 0.5 * (min + max);
+  result.center += mid[0] * result.axes.col(0) +
+                   mid[1] * result.axes.col(1) +
+                   mid[2] * result.axes.col(2);
   
     // reorder axes by length
   CartVect range = 0.5 * (max - min);
@@ -205,25 +213,25 @@ static ErrorCode box_from_axes( OrientedBox& result,
   {
     if (range[2] < range[0]) {
       std::swap( range[0], range[2] );
-      std::swap( result.axis[0], result.axis[2] );
+      result.axes.swapcol( 0, 2 );
     }
   }
   else if (range[1] < range[0]) {
     std::swap( range[0], range[1] );
-    std::swap( result.axis[0], result.axis[1] );
+    result.axes.swapcol( 0, 1 );
   }
   if (range[1] > range[2]) {
     std::swap( range[1], range[2] );
-    std::swap( result.axis[1], result.axis[2] );
+    result.axes.swapcol( 1, 2 );
   }
 
     // scale axis to encompass all points, divide in half
 #if MB_ORIENTED_BOX_UNIT_VECTORS
   result.length = range;
 #else
-  result.axis[0] *= range[0];
-  result.axis[1] *= range[1];
-  result.axis[2] *= range[2];
+  result.axes.colscale(0, range[0]);
+  result.axes.colscale(1, range[1]);
+  result.axes.colscale(2, range[2]);
 #endif
 
 #if MB_ORIENTED_BOX_OUTER_RADIUS
@@ -269,8 +277,10 @@ ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result,
   a /= count;
 
     // Get axes (Eigenvectors) from covariance matrix
-  double lambda[3];
-  moab::Matrix::EigenDecomp( a, lambda, result.axis );
+  CartVect lambda;
+  a.eigen_decomposition(lambda, result.axes);
+
+  // moab::Matrix::EigenDecomp( a, lambda, result.axes );
   
     // Calculate center and extents of box given orientation defined by axes
   return box_from_axes( result, instance, vertices );
@@ -351,8 +361,8 @@ ErrorCode OrientedBox::compute_from_covariance_data(
                                                 const Range& vertices )
 {
   if (data.area <= 0.0) {
-    CartVect axis[3] = { CartVect(0.), CartVect(0.), CartVect(0.) };
-    result = OrientedBox( axis, CartVect(0.) );
+    Matrix3 empty_axes (0.0);
+    result = OrientedBox( empty_axes, CartVect(0.) );
     return MB_SUCCESS;
   }
 
@@ -364,8 +374,8 @@ ErrorCode OrientedBox::compute_from_covariance_data(
   data.matrix -= outer_product( result.center, result.center );
 
     // get axes (Eigenvectors) from covariance matrix
-  double lamda[3];
-  moab::Matrix::EigenDecomp( data.matrix, lamda, result.axis );
+  CartVect lambda;
+  data.matrix.eigen_decomposition(lambda, result.axes);
 
     // We now have only the axes.  Calculate proper center
     // and extents for enclosed points.
@@ -376,13 +386,13 @@ bool OrientedBox::contained( const CartVect& point, double tol ) const
 {
   CartVect from_center = point - center;
 #if MB_ORIENTED_BOX_UNIT_VECTORS
-  return fabs(from_center % axis[0]) - length[0] <= tol &&
-         fabs(from_center % axis[1]) - length[1] <= tol &&
-         fabs(from_center % axis[2]) - length[2] <= tol ;
+  return fabs(from_center % axes.col(0)) - length[0] <= tol &&
+         fabs(from_center % axes.col(1)) - length[1] <= tol &&
+         fabs(from_center % axes.col(2)) - length[2] <= tol ;
 #else
   for (int i = 0; i < 3; ++i) {
-    double length = axis[i].length();
-    if (fabs(from_center % axis[i]) - length*length > length*tol)
+    double length = axes.col(i).length();
+    if (fabs(from_center % axes.col(i)) - length*length > length*tol)
       return false;
   }
   return true;
@@ -418,13 +428,13 @@ ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
 //      {
 //        CartVect corner( center );
 //#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
-//        corner += i * box.length[0] * box.axis[0];
-//        corner += j * box.length[1] * box.axis[1];
-//        corner += k * box.length[2] * box.axis[2];
+//        corner += i * box.length[0] * box.axis.col(0);
+//        corner += j * box.length[1] * box.axis.col(1);
+//        corner += k * box.length[2] * box.axis.col(2);
 //#else
-//        corner += i * box.axis[0];
-//        corner += j * box.axis[1];
-//        corner += k * box.axis[2];
+//        corner += i * box.axis.col(0);
+//        corner += j * box.axis.col(1);
+//        corner += k * box.axis.col(2);
 //#endif
 //        if (!contained( corner, tol ))
 //          return false;
@@ -552,26 +562,17 @@ bool OrientedBox::intersect_ray( const CartVect& ray_origin,
     }
   }
   
-    // get transpose of axes
-    // Note: if axes were stored as a matrix, could skip
-    // transpose and just switch order of operands in
-    // matrix-vector multiplies below. - J.K.
-  //Matrix3 B( axis[0][0], axis[1][0], axis[2][0],
-  //             axis[0][1], axis[1][1], axis[2][1],
-  //             axis[0][2], axis[1][2], axis[2][2] );
-  Matrix3 B( axis[0][0], axis[0][1], axis[0][2],
-               axis[1][0], axis[1][1], axis[1][2],
-               axis[2][0], axis[2][1], axis[2][2] );
-  //CartVect T = B * -center;
+  // get transpose of axes
+  Matrix3 B = Matrix::transpose(axes);
+
   
-    // transform ray to box coordintae system
-  //CartVect par_pos = T + B * b;
+  // transform ray to box coordintae system
   CartVect par_pos = B * (ray_origin - center);
   CartVect par_dir = B * ray_direction;
 
   // Fast Rejection Test: Ray will not intersect if it is going away from the box.
   // This will not work for rays with neg_ray_len. length[0] is half of box width 
-  // along axis[0].
+  // along axes.col(0).
   const double half_x = length[0] + reps;
   const double half_y = length[1] + reps;
   const double half_z = length[2] + reps;
@@ -658,9 +659,9 @@ ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
     CartVect coords(center);
     for (int j = 0; j < 3; ++j){
 #if MB_ORIENTED_BOX_UNIT_VECTORS
-      coords += signs[i][j] * (axis[j]*length[j]);
+      coords += signs[i][j] * (axes.col(j)*length[j]);
 #else
-      coords += signs[i][j] * axis[j];
+      coords += signs[i][j] * axes.col(j);
 #endif
     }
     EntityHandle handle;
@@ -689,9 +690,9 @@ void OrientedBox::closest_location_in_box(
   const CartVect from_center = input_position - center;
 
 #if MB_ORIENTED_BOX_UNIT_VECTORS
-  CartVect local( from_center % axis[0],
-                    from_center % axis[1],
-                    from_center % axis[2] );
+  CartVect local( from_center % axes.col(0),
+                    from_center % axes.col(1),
+                    from_center % axes.col(2) );
 
   for (int i = 0; i < 3; ++i) {
     if (local[i] < -length[i])
@@ -700,9 +701,9 @@ void OrientedBox::closest_location_in_box(
       local[i] =  length[i];
   }
 #else
-  CartVect local( (from_center % axis[0]) / (axis[0] % axis[0]),
-                    (from_center % axis[1]) / (axis[1] % axis[1]),
-                    (from_center % axis[2]) / (axis[2] % axis[2]) );
+  CartVect local( (from_center % axes.col(0)) / (axes.col(0) % axes.col(0)),
+                    (from_center % axes.col(1)) / (axes.col(1) % axes.col(1)),
+                    (from_center % axes.col(2)) / (axes.col(2) % axes.col(2)) );
 
   for (int i = 0; i < 3; ++i) {
     if (local[i] < -1.0)
@@ -713,9 +714,9 @@ void OrientedBox::closest_location_in_box(
 #endif
 
   output_position = center
-                  + local[0] * axis[0] 
-                  + local[1] * axis[1]
-                  + local[2] * axis[2];
+                  + local[0] * axes.col(0) 
+                  + local[1] * axes.col(1)
+                  + local[2] * axes.col(2);
 }
   
 } // namespace moab

src/OrientedBoxTreeTool.cpp

@@ -361,7 +361,7 @@ static ErrorCode split_box( Interface* instance,
       centroid += coords[j];
     centroid /= conn_len;
     
-    if ((box.axis[axis] % (centroid - box.center)) < 0.0)
+    if ((box.axis(axis) % (centroid - box.center)) < 0.0)
       left_list.insert( *i );
     else
       right_list.insert( *i );
@@ -380,7 +380,7 @@ ErrorCode OrientedBoxTreeTool::build_tree( const Range& entities,
   ErrorCode rval;
   
   if (entities.empty()) {
-    CartVect axis[3] = { CartVect(0.), CartVect(0.), CartVect(0.) };
+    Matrix3 axis;
     tmp_box = OrientedBox( axis, CartVect(0.) );
   }
   else {
@@ -472,7 +472,7 @@ static ErrorCode split_sets( Interface* ,
   std::list<OrientedBoxTreeTool::SetData>::const_iterator i;
   for (i = sets.begin(); i != sets.end(); ++i) {
     CartVect centroid(i->box_data.center / i->box_data.area);
-    if ((box.axis[axis] % (centroid - box.center)) < 0.0)
+    if ((box.axis(axis) % (centroid - box.center)) < 0.0)
       left.push_back( *i );
     else
       right.push_back( *i );
@@ -1879,9 +1879,9 @@ ErrorCode TreeNodePrinter::print_geometry( EntityHandle node )
   
   outputStream << box.center << "  Radius: " 
                << box.inner_radius() << " - " << box.outer_radius() << std::endl
-               << '+' << box.axis[0] << " : " << length[0] << std::endl
-               << 'x' << box.axis[1] << " : " << length[1] << std::endl
-               << 'x' << box.axis[2] << " : " << length[2] << std::endl;
+               << '+' << box.axis(0) << " : " << length[0] << std::endl
+               << 'x' << box.axis(1) << " : " << length[1] << std::endl
+               << 'x' << box.axis(2) << " : " << length[2] << std::endl;
   return MB_SUCCESS;
 }
 

src/SmoothFace.cpp

 #include <iostream>
 #include <fstream>
-#include "SmoothFace.hpp"
-
 #include <algorithm>
 #include <iomanip>
+#include <cassert>
+#include <limits>
+#include "SmoothFace.hpp"
 
-#include "assert.h"
-// included in the header now
-// #include "Range.hpp"
-// #include "CartVect.hpp"
-
-// some defines from CUBIT
 #define GEOMETRY_RESABS 1.e-6
-#define CUBIT_DBL_MAX 1.e+30
-//#define DBL_EPSILON  1.e-8
-#include <float.h>
+#define mbsqr(a) ((a)*(a))
+#define mbcube(a) (mbsqr(a) * (a))
+#define mbquart(a) (mbsqr(a) * mbsqr(a))
 
 namespace moab {
 
@@ -845,19 +840,19 @@ ErrorCode SmoothFace::eval_bezier_patch(EntityHandle tri, CartVect &areacoord,
   ctrl_pts[4] = vN[1]; //
 
   //i=4; j=0; k=0;
-  double B = quart(areacoord[0]);
+  double B = mbquart(areacoord[0]);
   pt += B * ctrl_pts[0];
 
   //i=3; j=1; k=0;
-  B = 4.0 * cube(areacoord[0]) * areacoord[1];
+  B = 4.0 * mbcube(areacoord[0]) * areacoord[1];
   pt += B * ctrl_pts[1];
 
   //i=2; j=2; k=0;
-  B = 6.0 * sqr(areacoord[0]) * sqr(areacoord[1]);
+  B = 6.0 * mbsqr(areacoord[0]) * mbsqr(areacoord[1]);
   pt += B * ctrl_pts[2];
 
   //i=1; j=3; k=0;
-  B = 4.0 * areacoord[0] * cube(areacoord[1]);
+  B = 4.0 * areacoord[0] * mbcube(areacoord[1]);
   pt += B * ctrl_pts[3];
 
   //edge = facet->edge(0);
@@ -869,19 +864,19 @@ ErrorCode SmoothFace::eval_bezier_patch(EntityHandle tri, CartVect &areacoord,
   ctrl_pts[4] = vN[2]; //
 
   //i=0; j=4; k=0;
-  B = quart(areacoord[1]);
+  B = mbquart(areacoord[1]);
   pt += B * ctrl_pts[0];
 
   //i=0; j=3; k=1;
-  B = 4.0 * cube(areacoord[1]) * areacoord[2];
+  B = 4.0 * mbcube(areacoord[1]) * areacoord[2];
   pt += B * ctrl_pts[1];
 
   //i=0; j=2; k=2;
-  B = 6.0 * sqr(areacoord[1]) * sqr(areacoord[2]);
+  B = 6.0 * mbsqr(areacoord[1]) * mbsqr(areacoord[2]);
   pt += B * ctrl_pts[2];
 
   //i=0; j=1; k=3;
-  B = 4.0 * areacoord[1] * cube(areacoord[2]);
+  B = 4.0 * areacoord[1] * mbcube(areacoord[2]);
   pt += B * ctrl_pts[3];
 
   //edge = facet->edge(1);
@@ -893,31 +888,31 @@ ErrorCode SmoothFace::eval_bezier_patch(EntityHandle tri, CartVect &areacoord,
   ctrl_pts[4] = vN[0]; //
 
   //i=0; j=0; k=4;
-  B = quart(areacoord[2]);
+  B = mbquart(areacoord[2]);
   pt += B * ctrl_pts[0];
 
   //i=1; j=0; k=3;
-  B = 4.0 * areacoord[0] * cube(areacoord[2]);
+  B = 4.0 * areacoord[0] * mbcube(areacoord[2]);
   pt += B * ctrl_pts[1];
 
   //i=2; j=0; k=2;
-  B = 6.0 * sqr(areacoord[0]) * sqr(areacoord[2]);
+  B = 6.0 * mbsqr(areacoord[0]) * mbsqr(areacoord[2]);
   pt += B * ctrl_pts[2];
 
   //i=3; j=0; k=1;
-  B = 4.0 * cube(areacoord[0]) * areacoord[2];
+  B = 4.0 * mbcube(areacoord[0]) * areacoord[2];
   pt += B * ctrl_pts[3];
 
   //i=2; j=1; k=1;
-  B = 12.0 * sqr(areacoord[0]) * areacoord[1] * areacoord[2];
+  B = 12.0 * mbsqr(areacoord[0]) * areacoord[1] * areacoord[2];
   pt += B * P_facet[0];
 
   //i=1; j=2; k=1;
-  B = 12.0 * areacoord[0] * sqr(areacoord[1]) * areacoord[2];
+  B = 12.0 * areacoord[0] * mbsqr(areacoord[1]) * areacoord[2];
   pt += B * P_facet[1];
 
   //i=1; j=1; k=2;
-  B = 12.0 * areacoord[0] * areacoord[1] * sqr(areacoord[2]);
+  B = 12.0 * areacoord[0] * areacoord[1] * mbsqr(areacoord[2]);
   pt += B * P_facet[2];
 
   return MB_SUCCESS;
@@ -1001,7 +996,7 @@ void SmoothFace::facet_area_coordinate(EntityHandle facet,
         p[2][1], p[2][2]);
     if (fabs(area2) < tol)