geometric_object_volume (#45)

* geom_object_volume function

* more tests

* some bug fixes

* more fixes

* don't use geom_object_volume in 1d and 2d cases

Author: stevengj

utils/ctlgeom.h

@@ -68,6 +68,7 @@ extern void geom_fix_lattice(void);
 extern void geom_fix_lattice0(LATTICE *L);
 extern void geom_cartesian_lattice(void);
 extern void geom_cartesian_lattice0(LATTICE *L);
+extern double geom_object_volume(GEOMETRIC_OBJECT o);
 extern boolean point_in_objectp(vector3 p, GEOMETRIC_OBJECT o);
 extern boolean point_in_periodic_objectp(vector3 p, GEOMETRIC_OBJECT o);
 extern boolean point_in_fixed_objectp(vector3 p, GEOMETRIC_OBJECT o);

utils/geom.c

@@ -1016,6 +1016,47 @@ matrix3x3 CTLIO square_basis(matrix3x3 basis, vector3 size)
   return matrix3x3_mult(matrix3x3_inverse(basis), square);
 }
 
+/**************************************************************************/
+
+/* compute the 3d volume enclosed by a geometric object o. */
+
+double geom_object_volume(GEOMETRIC_OBJECT o)
+{
+     switch (o.which_subclass) {
+	 case GEOM SPHERE: {
+	     number radius = o.subclass.sphere_data->radius;
+          return (1.333333333333333333 * K_PI) * radius*radius*radius;
+	 }
+	 case GEOM CYLINDER: {
+	     number height = o.subclass.cylinder_data->height;
+	     number radius = o.subclass.cylinder_data->radius;
+	     number radius2 = o.subclass.cylinder_data->which_subclass == CYL CONE ? o.subclass.cylinder_data->subclass.cone_data->radius2 : radius;
+          double vol = height * (K_PI/3) * (radius*radius + radius*radius2 + radius2*radius2);
+          if (o.subclass.cylinder_data->which_subclass == CYL WEDGE)
+               return vol * fabs(o.subclass.cylinder_data->subclass.wedge_data->wedge_angle) / (2*K_PI);
+          else
+               return vol;
+	 }
+	 case GEOM BLOCK: {
+          vector3 size = o.subclass.block_data->size;
+          double vol = size.x * size.y * size.z * fabs(matrix3x3_determinant(geometry_lattice.basis) / matrix3x3_determinant(o.subclass.block_data->projection_matrix));
+          return o.subclass.block_data->which_subclass == BLK BLOCK_SELF ? vol : vol * (K_PI/6);
+	 }
+	 case GEOM PRISM: {
+           vector3_list vertices = o.subclass.prism_data->vertices_p;
+           double area = 0;
+           int i;
+           for (i = 0; i < vertices.num_items; ++i) {
+                int i1 = (i + 1) % vertices.num_items;
+                area += 0.5 * (vertices.items[i1].x - vertices.items[i].x) * (vertices.items[i1].y + vertices.items[i].y);
+           }
+           return fabs(area) * o.subclass.prism_data->height;
+	 }
+	 default:
+	      return 0; /* unsupported object types? */
+     }
+}
+
 /**************************************************************************/
 /**************************************************************************/
 
@@ -1353,6 +1394,12 @@ number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o,
      unsigned i;
 
      geom_get_bounding_box(o, &bb);
+     if (!is_ellipsoid &&
+         !empty_x && !empty_y && !empty_z && /* todo: optimize 1d and 2d cases */
+         bb.low.x >= b.low.x && bb.high.x <= b.high.x &&
+         bb.low.y >= b.low.y && bb.high.y <= b.high.y &&
+         bb.low.z >= b.low.z && bb.high.z <= b.high.z)
+       return geom_object_volume(o) / (V0 * fabs(matrix3x3_determinant(geometry_lattice.basis))); /* o is completely contained within b */
      geom_box_intersection(&bb, &b, &bb);
      if (bb.low.x > bb.high.x || bb.low.y > bb.high.y || bb.low.z > bb.high.z
 	 || (!empty_x && bb.low.x == bb.high.x)

utils/geomtst.c

@@ -165,43 +165,42 @@ static double simple_ellip_overlap(geom_box b, geometric_object o, double tol)
      return olap0;
 }
 
-static void test_overlap(double tol,
-			 number (*box_overlap_with_object)
-			 (geom_box b, geometric_object o,
-			  number tol, integer maxeval),
-			 double (*simple_overlap)
-			 (geom_box b, geometric_object o, double tol))
+geometric_object random_object_and_lattice(void)
 {
      geometric_object o = random_object();
-     vector3 dir = random_unit_vector3();
-     geom_box b;
-     double d, olap, olap0;
-     int dim;
-
 #if 1
      geometry_lattice.basis1 = random_unit_vector3();
      geometry_lattice.basis2 = random_unit_vector3();
      geometry_lattice.basis3 = random_unit_vector3();
      geom_fix_lattice();
      geom_fix_object_ptr(&o);
 #endif
+     return o;
+}
 
-     b.low = make_vector3(myurand(-1,0), myurand(-1,0), myurand(-1,0));
-     b.high = make_vector3(myurand(0,1), myurand(0,1), myurand(0,1));
-
-     d = find_edge(o, dir, 10, tol);
-     b.low = vector3_plus(b.low, vector3_scale(d, dir));
-     b.high = vector3_plus(b.high, vector3_scale(d, dir));
-
-     dim = rand() % 3 + 1;
-     if (dim < 3)
-	  b.low.z = b.high.z = 0;
-     if (dim < 2)
-	  b.low.y = b.high.y = 0;
-
-     olap = box_overlap_with_object(b, o, tol/100, 10000/tol);
-     olap0 = simple_overlap(b, o, tol/2);
+static const char *object_name(geometric_object o)
+{
+     switch (o.which_subclass) {
+	 case CYLINDER:
+	      switch (o.subclass.cylinder_data->which_subclass) {
+		  case WEDGE: return "wedge";
+		  case CONE: return "cone";
+		  case CYLINDER_SELF: return "cylinder";
+	      }
+	 case SPHERE: return "sphere";
+	 case BLOCK:
+	      switch (o.subclass.block_data->which_subclass) {
+		  case ELLIPSOID: return "ellipsoid";
+		  case BLOCK_SELF: return "block";
+	      }
+	 case PRISM: return "prism";
+	 case COMPOUND_GEOMETRIC_OBJECT: return "compound object";
+	 default: return "geometric object";
+     }
+}
 
+void check_overlap(double tol, double olap0, double olap, int dim, geometric_object o, geom_box b)
+{
      if (fabs(olap0 - olap) > 2 * tol * fabs(olap)) {
 	  fprintf(stderr, "Large error %e in overlap (%g vs. %g) for:\n"
 		  "  lattice = (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)\n"
@@ -220,22 +219,59 @@ static void test_overlap(double tol,
 		  b.low.x, b.low.y, b.low.z,
 		  b.high.x, b.high.y, b.high.z);
 	  display_geometric_object_info(2, o);
-#if 1
-	  while (1) {
-	       tol /= sqrt(2.0);
-	       fprintf(stderr, "olap = %g, olap0 = %g (with tol = %e)\n",
-		       box_overlap_with_object(b, o, tol/100, 10000/tol),
-		       simple_overlap(b, o, tol/2), tol);
-	  }
-#endif
-	  exit(1);
+	  /* exit(1); */
      }
      else
-	  printf("Got %dd overlap %g vs. %g with tol = %e\n",
-		 dim,olap,olap0,tol);
+	  printf("Got %s %dd overlap %g vs. %g with tol = %e\n",
+		 object_name(o), dim,olap,olap0,tol);
+}
+
+static void test_overlap(double tol,
+			 number (*box_overlap_with_object)
+			 (geom_box b, geometric_object o,
+			  number tol, integer maxeval),
+			 double (*simple_overlap)
+			 (geom_box b, geometric_object o, double tol))
+{
+     geometric_object o = random_object_and_lattice();
+     vector3 dir = random_unit_vector3();
+     geom_box b;
+     double d, olap, olap0;
+     int dim;
+
+     b.low = make_vector3(myurand(-1,0), myurand(-1,0), myurand(-1,0));
+     b.high = make_vector3(myurand(0,1), myurand(0,1), myurand(0,1));
+     d = find_edge(o, dir, 10, tol);
+     b.low = vector3_plus(b.low, vector3_scale(d, dir));
+     b.high = vector3_plus(b.high, vector3_scale(d, dir));
+
+     dim = rand() % 3 + 1;
+     if (dim < 3)
+	  b.low.z = b.high.z = 0;
+     if (dim < 2)
+	  b.low.y = b.high.y = 0;
+
+     olap = box_overlap_with_object(b, o, tol/100, 10000/tol);
+     olap0 = simple_overlap(b, o, tol/2);
+     check_overlap(tol, olap0, olap, dim, o, b);
+     geometric_object_destroy(o);
+}
+
+static void test_volume(double tol)
+{
+     geometric_object o = random_object_and_lattice();
+     geom_box b;
+     double olap1, olap2;
+
+     geom_get_bounding_box(o, &b);
+     olap1 = box_overlap_with_object(b, o, tol/100, 10000/tol);
+     b.low.x += 1e-7 * (b.high.x - b.low.x); /* b no longer contains o */
+     olap2 = box_overlap_with_object(b, o, tol/100, 10000/tol);
+     check_overlap(tol, olap1, olap2, 3, o, b);
      geometric_object_destroy(o);
 }
 
+
 /************************************************************************/
 
 int main(void)
@@ -248,15 +284,18 @@ int main(void)
 
      geom_initialize();
 
+     printf("**** whole box overlap: ****\n");
+     for (i = 0; i < ntest; ++i)
+          test_volume(tol);
      for (i = 0; i < ntest; ++i) {
-	  printf("**** box overlap: ****\n");
-	  test_overlap(tol,
-		       box_overlap_with_object,
-		       simple_overlap);
-	  printf("**** ellipsoid overlap: ****\n");
-	  test_overlap(tol,
-		       ellipsoid_overlap_with_object,
-		       simple_ellip_overlap);
+          printf("**** box overlap: ****\n");
+          test_overlap(tol,
+                    box_overlap_with_object,
+                    simple_overlap);
+          printf("**** ellipsoid overlap: ****\n");
+          test_overlap(tol,
+                    ellipsoid_overlap_with_object,
+                    simple_ellip_overlap);
      }
 
      return 0;