diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h index b5e9bf4..10bf357 100644 --- a/utils/ctlgeom.h +++ b/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); diff --git a/utils/geom.c b/utils/geom.c index 2b510ad..9df47a2 100644 --- a/utils/geom.c +++ b/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) diff --git a/utils/geomtst.c b/utils/geomtst.c index 7303116..b2644f2 100644 --- a/utils/geomtst.c +++ b/utils/geomtst.c @@ -165,19 +165,9 @@ 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(); @@ -185,23 +175,32 @@ static void test_overlap(double tol, 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;