[MOAB-dev] r5969 - in MOAB/trunk: test/parallel tools/mbcslam

[iulian at mcs.anl.gov]

Author: iulian
Date: 2013-02-01 16:58:16 -0600 (Fri, 01 Feb 2013)
New Revision: 5969

Added:
   MOAB/trunk/tools/mbcslam/spherical_area_test.cpp
Modified:
   MOAB/trunk/test/parallel/par_intx_sph.cpp
   MOAB/trunk/tools/mbcslam/CslamUtils.cpp
   MOAB/trunk/tools/mbcslam/CslamUtils.hpp
   MOAB/trunk/tools/mbcslam/Makefile.am
Log:
add some tests for spherical area estimations
so far, the relative errors are about 1.e-14;
still need to test how accurate is Girard's theorem


Modified: MOAB/trunk/test/parallel/par_intx_sph.cpp
===================================================================
--- MOAB/trunk/test/parallel/par_intx_sph.cpp	2013-02-01 20:14:52 UTC (rev 5968)
+++ MOAB/trunk/test/parallel/par_intx_sph.cpp	2013-02-01 22:58:16 UTC (rev 5969)
@@ -264,5 +264,9 @@
   std::stringstream outf;
   outf<<"intersect" << rank<<".h5m";
   rval = mb.write_file(outf.str().c_str(), 0, 0, &outputSet, 1);
+  double intx_area = area_on_sphere(&mb, outputSet, radius);
+  double arrival_area = area_on_sphere(&mb, euler_set, radius) ;
+  std::cout<< "On rank : " << rank << " arrival area: " << arrival_area<<
+      "  intersection area:" << intx_area << " rel error: " << fabs((intx_area-arrival_area)/arrival_area) << "\n";
   CHECK_ERR(rval);
 }

Modified: MOAB/trunk/tools/mbcslam/CslamUtils.cpp
===================================================================
--- MOAB/trunk/tools/mbcslam/CslamUtils.cpp	2013-02-01 20:14:52 UTC (rev 5968)
+++ MOAB/trunk/tools/mbcslam/CslamUtils.cpp	2013-02-01 22:58:16 UTC (rev 5969)
@@ -665,4 +665,34 @@
   return Radius*Radius *correction;
 
 }
+
+double area_on_sphere(Interface * mb, EntityHandle set, double R)
+{
+  // get all entities of dimension 2
+  // then get the connectivity, etc
+  Range inputRange;
+  ErrorCode rval = mb->get_entities_by_dimension(set, 2, inputRange);
+  if (MB_SUCCESS != rval)
+    return -1;
+
+  // compare total area with 4*M_PI * R^2
+  double total_area = 0.;
+  for (Range::iterator eit=inputRange.begin(); eit!= inputRange.end(); eit++ )
+  {
+    EntityHandle eh=*eit;
+    // get the nodes, then the coordinates
+    const EntityHandle * verts;
+    int num_nodes;
+    rval = mb->get_connectivity(eh, verts, num_nodes);
+    if (MB_SUCCESS != rval)
+      return -1;
+    std::vector<double> coords(3*num_nodes);
+    // get coordinates
+    rval = mb->get_coords(verts, num_nodes, &coords[0]);
+    if (MB_SUCCESS != rval)
+      return -1;
+    total_area+=area_spherical_polygon (&coords[0], num_nodes, R);
+  }