Fix init_prism centroid; make test-prism's parallel loops thread-safe

utils/geom.c

@@ -247,7 +247,13 @@ void geom_initialize(void) {
    initialized.
 
    point_in_fixed_objectp additionally requires that geom_fix_object
-   has been called on o (if the lattice basis is non-orthogonal).  */
+   has been called on o (if the lattice basis is non-orthogonal).
+
+   NOT THREAD-SAFE: calls geom_fix_object_ptr, which mutates per-object
+   cached state (and for prisms, frees+reallocates shared arrays).
+   Callers that query the same object from multiple threads must
+   geom_fix_object_ptr once up-front and then call
+   point_in_fixed_objectp / point_in_fixed_pobjectp directly. */
 
 boolean CTLIO point_in_objectp(vector3 p, geometric_object o) {
   geom_fix_object_ptr(&o);
@@ -395,7 +401,13 @@ vector3 from_geom_object_coords(vector3 p, geometric_object o) {
    in lattice coordinates, using the surface of the object that the
    point is "closest" to for some definition of "closest" that is
    reasonable (at least for points near to the object). The length and
-   sign of the normal vector are arbitrary. */
+   sign of the normal vector are arbitrary.
+
+   NOT THREAD-SAFE: calls geom_fix_object_ptr, which mutates per-object
+   cached state (and for prisms, frees+reallocates shared arrays).
+   Callers that query the same object from multiple threads must
+   geom_fix_object_ptr once up-front and then call
+   normal_to_fixed_object directly. */
 
 vector3 CTLIO normal_to_object(vector3 p, geometric_object o) {
   geom_fix_object_ptr(&o);
@@ -523,7 +535,13 @@ vector3 normal_to_fixed_object(vector3 p, geometric_object o) {
 /**************************************************************************/
 
 /* Like point_in_objectp, but also checks the object shifted
-   by the lattice vectors: */
+   by the lattice vectors.
+
+   NOT THREAD-SAFE: calls geom_fix_object_ptr, which mutates per-object
+   cached state (and for prisms, frees+reallocates shared arrays).
+   Callers that query the same object from multiple threads must
+   geom_fix_object_ptr once up-front and then call
+   point_in_periodic_fixed_objectp directly. */
 
 boolean CTLIO point_in_periodic_objectp(vector3 p, geometric_object o) {
   geom_fix_object_ptr(&o);
@@ -2600,12 +2618,14 @@ void init_prism(geometric_object *o) {
     vertices = (vector3 *)realloc(vertices, num_vertices * sizeof(vector3));
   }
 
-  // compute centroid of vertices
+  // compute centroid of vertices; prsm->centroid is assigned below, AFTER
+  // the optional shift to o->center is applied (otherwise prsm->centroid
+  // would be stale relative to the shifted vertices).
   vector3 centroid = {0.0, 0.0, 0.0};
   int nv;
   for (nv = 0; nv < num_vertices; nv++)
     centroid = vector3_plus(centroid, vertices[nv]);
-  prsm->centroid = centroid = vector3_scale(1.0 / ((double)num_vertices), centroid);
+  centroid = vector3_scale(1.0 / ((double)num_vertices), centroid);
 
   // make sure all vertices lie in a plane, i.e. that the normal
   // vectors to all triangles (v_n, v_{n+1}, centroid) agree.
@@ -2656,6 +2676,7 @@ void init_prism(geometric_object *o) {
       vertices[nv] = vector3_plus(vertices[nv], shift);
     centroid = vector3_plus(centroid, shift);
   }
+  prsm->centroid = centroid;
 
   // compute rotation matrix that operates on a vector of cartesian coordinates
   // to yield the coordinates of the same point in the prism coordinate system.

utils/test-prism.c

@@ -280,11 +280,15 @@ int test_point_inclusion(geometric_object the_block, geometric_object the_prism,
   int num_failed = 0, num_adjusted = 0, n;
   char *exclude_env = getenv("LIBCTL_EXCLUDE_BOUNDARIES");
   boolean libctl_include_boundaries = !(exclude_env && exclude_env[0] == '1');
+  // Use point_in_fixed_objectp instead of point_in_objectp because the
+  // latter calls geom_fix_object_ptr internally and is not thread-safe.
+  // the_block and the_prism were fixed up in run_unit_tests above
+  // before this parallel loop runs.
 #pragma omp parallel for schedule(dynamic) reduction(+ : num_failed, num_adjusted)
   for (n = 0; n < num_tests; n++) {
     vector3 p = random_point_in_box(min_corner, max_corner);
-    boolean in_block = point_in_objectp(p, the_block);
-    boolean in_prism = point_in_objectp(p, the_prism);
+    boolean in_block = point_in_fixed_objectp(p, the_block);
+    boolean in_prism = point_in_fixed_objectp(p, the_prism);
 
     if (in_block != in_prism) {
       // retry with boundary exclusion/inclusion reversed
@@ -330,8 +334,10 @@ int test_normal_to_object(geometric_object the_block, geometric_object the_prism
     vector3 p = (urand(0.0, 1.0) < PFACE) ? random_point_on_prism(the_prism)
                                           : random_point_in_box(min_corner, max_corner);
 
-    vector3 nhat_block = standardize(normal_to_object(p, the_block));
-    vector3 nhat_prism = standardize(normal_to_object(p, the_prism));
+    // normal_to_fixed_object instead of normal_to_object: the latter is
+    // not thread-safe (calls geom_fix_object_ptr internally).
+    vector3 nhat_block = standardize(normal_to_fixed_object(p, the_block));
+    vector3 nhat_prism = standardize(normal_to_fixed_object(p, the_prism));
     if (!vector3_nearly_equal(nhat_block, nhat_prism, tolerance)) num_failed++;
 
     if (f) {
@@ -1150,6 +1156,11 @@ int run_unit_tests() {
     vector3 shift = vector3_scale(urand(0.0, 1.0), random_unit_vector3());
     the_block.center = vector3_plus(the_block.center, shift);
     the_prism.center = vector3_plus(the_prism.center, shift);
+    // Sync each object's internal cache (block projection_matrix, prism
+    // vertices_p / centroid) to the new center, once, single-threaded,
+    // before the parallel test loops below run.
+    geom_fix_object_ptr(&the_block);
+    geom_fix_object_ptr(&the_prism);
   }
 
   char *s = getenv("LIBCTL_TEST_PRISM_LOG");