From 3048dfb1776b76e384f08952a1b79644c335f8cc Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 2 Oct 2019 18:13:04 +0800 Subject: [PATCH 1/5] geom_object_volume function --- utils/ctlgeom.h | 1 + utils/geom.c | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+) 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..0296332 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 (4/3 * 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 * 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? */ + } +} + /**************************************************************************/ /**************************************************************************/ From 5848216a4a047390f5e76a2376d48627c8d1c42c Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 22 Oct 2019 21:12:23 -0400 Subject: [PATCH 2/5] more tests --- utils/geom.c | 5 +++ utils/geomtst.c | 104 ++++++++++++++++++++++++++++-------------------- 2 files changed, 66 insertions(+), 43 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 0296332..07e3d6e 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -1394,6 +1394,11 @@ number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o, unsigned i; geom_get_bounding_box(o, &bb); + if (!is_ellipsoid && + 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); /* 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..05e855b 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,11 @@ 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); - +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 +198,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); } else printf("Got %dd overlap %g vs. %g with tol = %e\n", 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 += (1 + 1e-10) * (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) @@ -249,15 +264,18 @@ int main(void) geom_initialize(); 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); } + printf("**** whole box overlap: ****\n"); + for (i = 0; i < ntest; ++i) + test_volume(tol); return 0; } From e6f8f4add4dab760386572c091f405170ade0fa4 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 29 Oct 2019 21:24:23 -0400 Subject: [PATCH 3/5] some bug fixes --- utils/geom.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 07e3d6e..fd8851c 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -1025,7 +1025,7 @@ double geom_object_volume(GEOMETRIC_OBJECT o) switch (o.which_subclass) { case GEOM SPHERE: { number radius = o.subclass.sphere_data->radius; - return (4/3 * K_PI) * radius*radius*radius; + return (1.333333333333333333 * K_PI) * radius*radius*radius; } case GEOM CYLINDER: { number height = o.subclass.cylinder_data->height; @@ -1039,7 +1039,7 @@ double geom_object_volume(GEOMETRIC_OBJECT o) } case GEOM BLOCK: { vector3 size = o.subclass.block_data->size; - double vol = size.x * size.y * size.z * matrix3x3_determinant(geometry_lattice.basis) / matrix3x3_determinant(o.subclass.block_data->projection_matrix); + 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: { @@ -1398,7 +1398,7 @@ number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o, 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); /* o is completely contained within b */ + return geom_object_volume(o) / V0; /* 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) From 47ff10a7c137292f596526beca1e35185f871bf8 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 29 Oct 2019 21:33:43 -0400 Subject: [PATCH 4/5] more fixes --- utils/geom.c | 2 +- utils/geomtst.c | 35 ++++++++++++++++++++++++++++------- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index fd8851c..8f59691 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -1398,7 +1398,7 @@ number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o, 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; /* o is completely contained within b */ + 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 05e855b..b2644f2 100644 --- a/utils/geomtst.c +++ b/utils/geomtst.c @@ -178,6 +178,27 @@ geometric_object random_object_and_lattice(void) return o; } +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)) { @@ -198,11 +219,11 @@ void check_overlap(double tol, double olap0, double olap, int dim, geometric_obj b.low.x, b.low.y, b.low.z, b.high.x, b.high.y, b.high.z); display_geometric_object_info(2, o); - 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, @@ -244,7 +265,7 @@ static void test_volume(double tol) geom_get_bounding_box(o, &b); olap1 = box_overlap_with_object(b, o, tol/100, 10000/tol); - b.low.x += (1 + 1e-10) * (b.high.x - b.low.x); /* b no longer contains o */ + 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); @@ -263,6 +284,9 @@ 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, @@ -273,9 +297,6 @@ int main(void) ellipsoid_overlap_with_object, simple_ellip_overlap); } - printf("**** whole box overlap: ****\n"); - for (i = 0; i < ntest; ++i) - test_volume(tol); return 0; } From 4a616ae6974f607fbe36d4f5637d9fc1bdae2a41 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 5 Nov 2019 21:07:06 -0500 Subject: [PATCH 5/5] don't use geom_object_volume in 1d and 2d cases --- utils/geom.c | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/geom.c b/utils/geom.c index 8f59691..9df47a2 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -1395,6 +1395,7 @@ number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o, 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)