Commit b29ef89c863fcdf7b689b2b0438a365b43268065 - libctl

Import Debian changes 4.4.0-1 libctl (4.4.0-1) experimental; urgency=medium * New upstream release Thorsten Alteholz 4 years ago

7 changed file(s) with 162 addition(s) and 55 deletion(s). Raw diff Collapse all Expand all

-0

NEWS.md less more

0	0	# Libctl Release Notes
	1
	2	## libctl 4.4.0
	3
	4	11/12/19
	5
	6	* `geom_object_volume` function to get the volume of a 3d object
	7	(accelerates `box_overlap_with_object` for objects completely within a box) (#45).
	8
	9	* Bugfix to geometry tree search for empty dimensions.
1	10
2	11	## libctl 4.3.0
3	12

-2

configure.ac less more

0	0	# Process this file with autoconf to produce a configure script.
1		AC_INIT(libctl, 4.3.0, stevenj@alum.mit.edu)
	1	AC_INIT(libctl, 4.4.0, stevenj@alum.mit.edu)
2	2	AC_CONFIG_SRCDIR([src/ctl.c])
3	3	AC_CONFIG_HEADERS([config.h src/ctl.h])
4	4	AC_CONFIG_MACRO_DIR([m4])

7	7	# Shared-library version number; indicates api compatibility, and is
8	8	# not the same as the "public" version number. (Don't worry about this
9	9	# except for public releases.)
10		SHARED_VERSION_INFO="8:0:1" # CURRENT:REVISION:AGE
	10	SHARED_VERSION_INFO="9:0:2" # CURRENT:REVISION:AGE
11	11
12	12	AM_INIT_AUTOMAKE([foreign])
13	13	AC_SUBST(SHARED_VERSION_INFO)

-0

debian/changelog less more

	0	libctl (4.4.0-1) experimental; urgency=medium
	1
	2	* New upstream release
	3
	4	-- Thorsten Alteholz <debian@alteholz.de> Wed, 20 Nov 2019 19:25:46 +0100
	5
0	6	libctl (4.3.0-2) experimental; urgency=medium
1	7
2	8	* deactivate tests for i386 and powerpc for now

-0

debian/libctl7.symbols less more

209	209	geom_fix_objects@Base 4.1.0
210	210	geom_get_bounding_box@Base 4.1.0
211	211	geom_initialize@Base 4.1.0
	212	geom_object_volume@Base 4.4.0
212	213	geom_tree_search@Base 4.1.0
213	214	geom_tree_search_next@Base 4.1.0
214	215	geometric_object_copy@Base 4.1.0

-0

utils/ctlgeom.h less more

67	67	extern void geom_fix_lattice0(LATTICE *L);
68	68	extern void geom_cartesian_lattice(void);
69	69	extern void geom_cartesian_lattice0(LATTICE *L);
	70	extern double geom_object_volume(GEOMETRIC_OBJECT o);
70	71	extern boolean point_in_objectp(vector3 p, GEOMETRIC_OBJECT o);
71	72	extern boolean point_in_periodic_objectp(vector3 p, GEOMETRIC_OBJECT o);
72	73	extern boolean point_in_fixed_objectp(vector3 p, GEOMETRIC_OBJECT o);

+57

-6

utils/geom.c less more

573	573	for (kkk = -1; kkk <= 1; ++kkk) { \
574	574	shiftby.z = kkk * geometry_lattice.size.z; \
575	575	body; \
	576	if (geometry_lattice.size.z == 0) break; \
576	577	} \
	578	if (geometry_lattice.size.y == 0) break; \
577	579	} \
	580	if (geometry_lattice.size.x == 0) break; \
578	581	} \
579	582	break; \
580	583	} \

1013	1016	}
1014	1017
1015	1018	/**************************************************************************/
	1019
	1020	/* compute the 3d volume enclosed by a geometric object o. */
	1021
	1022	double geom_object_volume(GEOMETRIC_OBJECT o)
	1023	{
	1024	switch (o.which_subclass) {
	1025	case GEOM SPHERE: {
	1026	number radius = o.subclass.sphere_data->radius;
	1027	return (1.333333333333333333 * K_PI) * radiusradiusradius;
	1028	}
	1029	case GEOM CYLINDER: {
	1030	number height = o.subclass.cylinder_data->height;
	1031	number radius = o.subclass.cylinder_data->radius;
	1032	number radius2 = o.subclass.cylinder_data->which_subclass == CYL CONE ? o.subclass.cylinder_data->subclass.cone_data->radius2 : radius;
	1033	double vol = height * (K_PI/3) * (radiusradius + radiusradius2 + radius2*radius2);
	1034	if (o.subclass.cylinder_data->which_subclass == CYL WEDGE)
	1035	return vol * fabs(o.subclass.cylinder_data->subclass.wedge_data->wedge_angle) / (2*K_PI);
	1036	else
	1037	return vol;
	1038	}
	1039	case GEOM BLOCK: {
	1040	vector3 size = o.subclass.block_data->size;
	1041	double vol = size.x * size.y * size.z * fabs(matrix3x3_determinant(geometry_lattice.basis) / matrix3x3_determinant(o.subclass.block_data->projection_matrix));
	1042	return o.subclass.block_data->which_subclass == BLK BLOCK_SELF ? vol : vol * (K_PI/6);
	1043	}
	1044	case GEOM PRISM: {
	1045	vector3_list vertices = o.subclass.prism_data->vertices_p;
	1046	double area = 0;
	1047	int i;
	1048	for (i = 0; i < vertices.num_items; ++i) {
	1049	int i1 = (i + 1) % vertices.num_items;
	1050	area += 0.5 * (vertices.items[i1].x - vertices.items[i].x) * (vertices.items[i1].y + vertices.items[i].y);
	1051	}
	1052	return fabs(area) * o.subclass.prism_data->height;
	1053	}
	1054	default:
	1055	return 0; /* unsupported object types? */
	1056	}
	1057	}
	1058
	1059	/**************************************************************************/
1016	1060	/**************************************************************************/
1017	1061
1018	1062	/* Fast geometry routines */

1349	1393	unsigned i;
1350	1394
1351	1395	geom_get_bounding_box(o, &bb);
	1396	if (!is_ellipsoid &&
	1397	!empty_x && !empty_y && !empty_z && /* todo: optimize 1d and 2d cases */
	1398	bb.low.x >= b.low.x && bb.high.x <= b.high.x &&
	1399	bb.low.y >= b.low.y && bb.high.y <= b.high.y &&
	1400	bb.low.z >= b.low.z && bb.high.z <= b.high.z)
	1401	return geom_object_volume(o) / (V0 * fabs(matrix3x3_determinant(geometry_lattice.basis))); /* o is completely contained within b */
1352	1402	geom_box_intersection(&bb, &b, &bb);
1353	1403	if (bb.low.x > bb.high.x \|\| bb.low.y > bb.high.y \|\| bb.low.z > bb.high.z
1354	1404	\|\| (!empty_x && bb.low.x == bb.high.x)

1606	1656	{
1607	1657	int division_nobjects[3][2] = {{0,0},{0,0},{0,0}};
1608	1658	number division_point[3];
1609		int best = 0;
	1659	int best = -1;
1610	1660	int i, j, n1, n2;
1611	1661
1612	1662	if (!t)

1624	1674	number of objects in the partitioned boxes and finding
1625	1675	the best partition. */
1626	1676	for (i = 0; i < dimensions; ++i) {
	1677	if (VEC_I(t->b.high, i) == VEC_I(t->b.low, i)) continue; /* skip empty dimensions */
1627	1678	find_best_partition(t->nobjects, t->objects, i, &division_point[i],
1628	1679	&division_nobjects[i][0],
1629	1680	&division_nobjects[i][1]);
1630		if (MAX(division_nobjects[i][0], division_nobjects[i][1]) <
	1681	if (best < 0 \|\|
	1682	MAX(division_nobjects[i][0], division_nobjects[i][1]) <
1631	1683	MAX(division_nobjects[best][0], division_nobjects[best][1]))
1632	1684	best = i;
1633	1685	}
1634	1686
1635	1687	/* don't do anything if division makes the worst case worse or if
1636	1688	it fails to improve the best case: */
1637		if (MAX(division_nobjects[best][0], division_nobjects[best][1]) + 1 >
1638		t->nobjects \|\|
1639		MIN(division_nobjects[best][0], division_nobjects[best][1]) + 1 >=
1640		t->nobjects)
	1689	if (best < 0 \|\|
	1690	MAX(division_nobjects[best][0], division_nobjects[best][1]) + 1 > t->nobjects \|\|
	1691	MIN(division_nobjects[best][0], division_nobjects[best][1]) + 1 >= t->nobjects)
1641	1692	return; /* division didn't help us */
1642	1693
1643	1694	divide_geom_box(&t->b, best, division_point[best], &t->b1, &t->b2);

+86

-47

utils/geomtst.c less more

164	164	return olap0;
165	165	}
166	166
167		static void test_overlap(double tol,
168		number (*box_overlap_with_object)
169		(geom_box b, geometric_object o,
170		number tol, integer maxeval),
171		double (*simple_overlap)
172		(geom_box b, geometric_object o, double tol))
	167	geometric_object random_object_and_lattice(void)
173	168	{
174	169	geometric_object o = random_object();
175		vector3 dir = random_unit_vector3();
176		geom_box b;
177		double d, olap, olap0;
178		int dim;
179
180	170	#if 1
181	171	geometry_lattice.basis1 = random_unit_vector3();
182	172	geometry_lattice.basis2 = random_unit_vector3();

184	174	geom_fix_lattice();
185	175	geom_fix_object_ptr(&o);
186	176	#endif
187
188		b.low = make_vector3(myurand(-1,0), myurand(-1,0), myurand(-1,0));
189		b.high = make_vector3(myurand(0,1), myurand(0,1), myurand(0,1));
190
191		d = find_edge(o, dir, 10, tol);
192		b.low = vector3_plus(b.low, vector3_scale(d, dir));
193		b.high = vector3_plus(b.high, vector3_scale(d, dir));
194
195		dim = rand() % 3 + 1;
196		if (dim < 3)
197		b.low.z = b.high.z = 0;
198		if (dim < 2)
199		b.low.y = b.high.y = 0;
200
201		olap = box_overlap_with_object(b, o, tol/100, 10000/tol);
202		olap0 = simple_overlap(b, o, tol/2);
203
	177	return o;
	178	}
	179
	180	static const char *object_name(geometric_object o)
	181	{
	182	switch (o.which_subclass) {
	183	case CYLINDER:
	184	switch (o.subclass.cylinder_data->which_subclass) {
	185	case WEDGE: return "wedge";
	186	case CONE: return "cone";
	187	case CYLINDER_SELF: return "cylinder";
	188	}
	189	case SPHERE: return "sphere";
	190	case BLOCK:
	191	switch (o.subclass.block_data->which_subclass) {
	192	case ELLIPSOID: return "ellipsoid";
	193	case BLOCK_SELF: return "block";
	194	}
	195	case PRISM: return "prism";
	196	case COMPOUND_GEOMETRIC_OBJECT: return "compound object";
	197	default: return "geometric object";
	198	}
	199	}
	200
	201	void check_overlap(double tol, double olap0, double olap, int dim, geometric_object o, geom_box b)
	202	{
204	203	if (fabs(olap0 - olap) > 2 * tol * fabs(olap)) {
205	204	fprintf(stderr, "Large error %e in overlap (%g vs. %g) for:\n"
206	205	" lattice = (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)\n"

219	218	b.low.x, b.low.y, b.low.z,
220	219	b.high.x, b.high.y, b.high.z);
221	220	display_geometric_object_info(2, o);
222		#if 1
223		while (1) {
224		tol /= sqrt(2.0);
225		fprintf(stderr, "olap = %g, olap0 = %g (with tol = %e)\n",
226		box_overlap_with_object(b, o, tol/100, 10000/tol),
227		simple_overlap(b, o, tol/2), tol);
228		}
229		#endif
230		exit(1);
	221	/* exit(1); */
231	222	}
232	223	else
233		printf("Got %dd overlap %g vs. %g with tol = %e\n",
234		dim,olap,olap0,tol);
	224	printf("Got %s %dd overlap %g vs. %g with tol = %e\n",
	225	object_name(o), dim,olap,olap0,tol);
	226	}
	227
	228	static void test_overlap(double tol,
	229	number (*box_overlap_with_object)
	230	(geom_box b, geometric_object o,
	231	number tol, integer maxeval),
	232	double (*simple_overlap)
	233	(geom_box b, geometric_object o, double tol))
	234	{
	235	geometric_object o = random_object_and_lattice();
	236	vector3 dir = random_unit_vector3();
	237	geom_box b;
	238	double d, olap, olap0;
	239	int dim;
	240
	241	b.low = make_vector3(myurand(-1,0), myurand(-1,0), myurand(-1,0));
	242	b.high = make_vector3(myurand(0,1), myurand(0,1), myurand(0,1));
	243	d = find_edge(o, dir, 10, tol);
	244	b.low = vector3_plus(b.low, vector3_scale(d, dir));
	245	b.high = vector3_plus(b.high, vector3_scale(d, dir));
	246
	247	dim = rand() % 3 + 1;
	248	if (dim < 3)
	249	b.low.z = b.high.z = 0;
	250	if (dim < 2)
	251	b.low.y = b.high.y = 0;
	252
	253	olap = box_overlap_with_object(b, o, tol/100, 10000/tol);
	254	olap0 = simple_overlap(b, o, tol/2);
	255	check_overlap(tol, olap0, olap, dim, o, b);
235	256	geometric_object_destroy(o);
236	257	}
	258
	259	static void test_volume(double tol)
	260	{
	261	geometric_object o = random_object_and_lattice();
	262	geom_box b;
	263	double olap1, olap2;
	264
	265	geom_get_bounding_box(o, &b);
	266	olap1 = box_overlap_with_object(b, o, tol/100, 10000/tol);
	267	b.low.x += 1e-7 * (b.high.x - b.low.x); /* b no longer contains o */
	268	olap2 = box_overlap_with_object(b, o, tol/100, 10000/tol);
	269	check_overlap(tol, olap1, olap2, 3, o, b);
	270	geometric_object_destroy(o);
	271	}
	272
237	273
238	274	/************************************************************************/
239	275

247	283
248	284	geom_initialize();
249	285
	286	printf("** whole box overlap: **\n");
	287	for (i = 0; i < ntest; ++i)
	288	test_volume(tol);
250	289	for (i = 0; i < ntest; ++i) {
251		printf("** box overlap: **\n");
252		test_overlap(tol,
253		box_overlap_with_object,
254		simple_overlap);
255		printf("** ellipsoid overlap: **\n");
256		test_overlap(tol,
257		ellipsoid_overlap_with_object,
258		simple_ellip_overlap);
	290	printf("** box overlap: **\n");
	291	test_overlap(tol,
	292	box_overlap_with_object,
	293	simple_overlap);
	294	printf("** ellipsoid overlap: **\n");
	295	test_overlap(tol,
	296	ellipsoid_overlap_with_object,
	297	simple_ellip_overlap);
259	298	}
260	299
261	300	return 0;