star-uvm

Spatial structuring of unstructured volumetric meshes
git clone git://git.meso-star.fr/star-uvm.git
Log | Files | Refs | README | LICENSE

test_suvm_volume.c (31173B)


      1 /* Copyright (C) 2020-2023 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "suvm.h"
     17 #include "test_suvm_utils.h"
     18 #include "test_suvm_ball.h"
     19 #include "test_suvm_box.h"
     20 
     21 #include <rsys/dynamic_array.h>
     22 #include <rsys/double3.h>
     23 #include <rsys/float3.h>
     24 
     25 #define DARRAY_NAME prim
     26 #define DARRAY_DATA struct suvm_primitive
     27 #include <rsys/dynamic_array.h>
     28 
     29 struct mesh {
     30   const double* vertices; /* List of double[3] */
     31   size_t nvertices;
     32   const size_t* tetrahedra; /* List of size_t[4] */
     33   size_t ntetrahedra;
     34   size_t tetrahedron_data_alignment;
     35   size_t vertex_data_alignment;
     36 };
     37 static const struct mesh MESH_NULL;
     38 
     39 /*******************************************************************************
     40  * Geometry functions
     41  ******************************************************************************/
     42 static void
     43 get_indices(const size_t itetra, size_t ids[4], void* ctx)
     44 {
     45   struct mesh* msh = ctx;
     46   CHK(itetra < msh->ntetrahedra);
     47   CHK(ids != NULL);
     48   ids[0] = msh->tetrahedra[itetra*4+0];
     49   ids[1] = msh->tetrahedra[itetra*4+1];
     50   ids[2] = msh->tetrahedra[itetra*4+2];
     51   ids[3] = msh->tetrahedra[itetra*4+3];
     52 }
     53 
     54 static void
     55 get_position(const size_t ivert, double pos[3], void* ctx)
     56 {
     57   struct mesh* msh = ctx;
     58   CHK(ctx != NULL);
     59   CHK(pos != NULL);
     60   CHK(ivert < msh->nvertices);
     61   pos[0] = msh->vertices[ivert*3+0];
     62   pos[1] = msh->vertices[ivert*3+1];
     63   pos[2] = msh->vertices[ivert*3+2];
     64 }
     65 
     66 /* Use tetra indices as tetra data */
     67 static void
     68 get_tetra_data(const size_t itetra, void* data, void* ctx)
     69 {
     70   struct mesh* msh = ctx;
     71   size_t* ids = data;
     72   CHK(ctx != NULL);
     73   CHK(data != NULL);
     74   CHK(itetra < msh->ntetrahedra);
     75   CHK(IS_ALIGNED(data, msh->tetrahedron_data_alignment));
     76   ids[0] = msh->tetrahedra[itetra*4+0];
     77   ids[1] = msh->tetrahedra[itetra*4+1];
     78   ids[2] = msh->tetrahedra[itetra*4+2];
     79   ids[3] = msh->tetrahedra[itetra*4+3];
     80 }
     81 
     82 /* Use vertex position as vertex data */
     83 static void
     84 get_vert_data(const size_t ivert, void* data, void* ctx)
     85 {
     86   struct mesh* msh = ctx;
     87   double* pos = data;
     88   CHK(ctx != NULL);
     89   CHK(data != NULL);
     90   CHK(ivert < msh->nvertices);
     91   CHK(IS_ALIGNED(data, msh->vertex_data_alignment));
     92   pos[0] = msh->vertices[ivert*3+0];
     93   pos[1] = msh->vertices[ivert*3+1];
     94   pos[2] = msh->vertices[ivert*3+2];
     95 }
     96 
     97 /*******************************************************************************
     98  * Helper functions
     99  ******************************************************************************/
    100 static INLINE int
    101 cmp_prim(const void* a, const void* b)
    102 {
    103   const struct suvm_primitive* prim0 = a;
    104   const struct suvm_primitive* prim1 = b;
    105   CHK(prim0 && prim1);
    106   return prim0->iprim < prim1->iprim ? -1
    107        : prim0->iprim > prim1->iprim ? +1 : 0;
    108 }
    109 
    110 static void
    111 check_volume_polyhedra(const struct suvm_volume* vol, struct mesh* msh)
    112 {
    113   size_t iprim;
    114   size_t nprims;
    115 
    116   CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK);
    117   CHK(nprims == msh->ntetrahedra);
    118 
    119   FOR_EACH(iprim, 0, nprims) {
    120     struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
    121     struct suvm_polyhedron poly;
    122     double v[4][3];
    123     double E0[3];
    124     double E1[3];
    125     double N[3];
    126     size_t ids[4];
    127     float f[3];
    128     CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
    129     CHK(!SUVM_PRIMITIVE_NONE(&prim));
    130     CHK(prim.iprim == iprim);
    131     CHK(prim.nvertices == 4);
    132     CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
    133 
    134     get_indices(iprim, ids, msh);
    135     get_position(ids[0], v[0], msh);
    136     get_position(ids[1], v[1], msh);
    137     get_position(ids[2], v[2], msh);
    138     get_position(ids[3], v[3], msh);
    139 
    140     CHK(f3_eq(poly.v[0], f3_set_d3(f, v[0])));
    141     CHK(f3_eq(poly.v[1], f3_set_d3(f, v[1])));
    142     CHK(f3_eq(poly.v[2], f3_set_d3(f, v[2])));
    143     CHK(f3_eq(poly.v[3], f3_set_d3(f, v[3])));
    144 
    145     d3_sub(E0, v[1], v[0]);
    146     d3_sub(E1, v[2], v[0]);
    147     d3_cross(N, E0, E1);
    148     f3_normalize(f, f3_set_d3(f, N));
    149     CHK(f3_eq_eps(poly.N[0], f, 1.e-4f));
    150 
    151     d3_sub(E0, v[3], v[0]);
    152     d3_sub(E1, v[1], v[0]);
    153     d3_cross(N, E0, E1);
    154     f3_normalize(f, f3_set_d3(f, N));
    155     CHK(f3_eq_eps(poly.N[1], f, 1.e-4f));
    156 
    157     d3_sub(E0, v[3], v[1]);
    158     d3_sub(E1, v[2], v[1]);
    159     d3_cross(N, E0, E1);
    160     f3_normalize(f, f3_set_d3(f, N));
    161     CHK(f3_eq_eps(poly.N[2], f, 1.e-4f));
    162 
    163     d3_sub(E0, v[3], v[2]);
    164     d3_sub(E1, v[0], v[2]);
    165     d3_cross(N, E0, E1);
    166     f3_normalize(f, f3_set_d3(f, N));
    167     CHK(f3_eq_eps(poly.N[3], f, 1.e-4f));
    168   }
    169 }
    170 
    171 static void
    172 check_volume_mesh(const struct suvm_volume* vol, struct mesh* msh)
    173 {
    174   struct suvm_mesh_desc desc = SUVM_MESH_DESC_NULL;
    175   size_t i;
    176 
    177   CHK(suvm_volume_get_mesh_desc(vol, &desc) == RES_OK);
    178   CHK(desc.nvertices == msh->nvertices);
    179   CHK(desc.nprimitives == msh->ntetrahedra);
    180   CHK(desc.dvertex == 3);
    181   CHK(desc.dprimitive == 4);
    182 
    183   FOR_EACH(i, 0, desc.nvertices) {
    184     CHK(desc.positions[i*3+0] == (float)msh->vertices[i*3+0]);
    185     CHK(desc.positions[i*3+1] == (float)msh->vertices[i*3+1]);
    186     CHK(desc.positions[i*3+2] == (float)msh->vertices[i*3+2]);
    187   }
    188   FOR_EACH(i, 0, desc.nprimitives) {
    189     CHK(desc.indices[i*4+0] == (uint32_t)msh->tetrahedra[i*4+0]);
    190     CHK(desc.indices[i*4+1] == (uint32_t)msh->tetrahedra[i*4+1]);
    191     CHK(desc.indices[i*4+2] == (uint32_t)msh->tetrahedra[i*4+2]);
    192     CHK(desc.indices[i*4+3] == (uint32_t)msh->tetrahedra[i*4+3]);
    193   }
    194 }
    195 
    196 static void
    197 test_tetrahedral_mesh_creation(struct suvm_device* dev)
    198 {
    199   struct mesh msh = MESH_NULL;
    200   struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
    201   struct suvm_volume* vol = NULL;
    202   struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
    203   struct suvm_mesh_desc desc = SUVM_MESH_DESC_NULL;
    204   struct suvm_polyhedron poly;
    205   size_t nprims;
    206 
    207   args.ntetrahedra = box_ntetras;
    208   args.nvertices = box_nverts;
    209   args.get_indices = get_indices;
    210   args.get_position = get_position;
    211   args.tetrahedron_data = SUVM_DATA_NULL;
    212   args.vertex_data = SUVM_DATA_NULL;
    213   args.context = &msh;
    214 
    215   msh.vertices = box_vertices;
    216   msh.nvertices = box_nverts;
    217   msh.tetrahedra = box_indices;
    218   msh.ntetrahedra = box_ntetras;
    219 
    220   CHK(suvm_tetrahedral_mesh_create(NULL, &args, &vol) == RES_BAD_ARG);
    221   CHK(suvm_tetrahedral_mesh_create(dev, NULL, &vol) == RES_BAD_ARG);
    222   CHK(suvm_tetrahedral_mesh_create(dev, &args, NULL) == RES_BAD_ARG);
    223   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    224 
    225   CHK(suvm_volume_ref_get(NULL) == RES_BAD_ARG);
    226   CHK(suvm_volume_ref_get(vol) == RES_OK);
    227   CHK(suvm_volume_ref_put(NULL) == RES_BAD_ARG);
    228   CHK(suvm_volume_ref_put(vol) == RES_OK);
    229   CHK(suvm_volume_ref_put(vol) == RES_OK);
    230 
    231   args.ntetrahedra = 0;
    232   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    233   args.ntetrahedra = box_ntetras;
    234   args.nvertices = 0;
    235   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    236   args.nvertices = box_nverts;
    237   args.get_indices = NULL;
    238   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    239   args.get_indices = get_indices;
    240   args.get_position = NULL;
    241   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    242   args.get_position = get_position;
    243   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    244   CHK(suvm_volume_ref_put(vol) == RES_OK);
    245 
    246   args.precompute_normals = 1;
    247   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    248   CHK(suvm_volume_ref_put(vol) == RES_OK);
    249 
    250   args.tetrahedron_data.get = get_tetra_data;
    251   args.tetrahedron_data.size = sizeof(size_t[4]);
    252   args.tetrahedron_data.alignment = 64;
    253   args.vertex_data.get = get_vert_data;
    254   args.vertex_data.size = sizeof(double[3]);
    255   args.vertex_data.alignment = 32;
    256   msh.tetrahedron_data_alignment = args.tetrahedron_data.alignment;
    257   msh.vertex_data_alignment = args.vertex_data.alignment;
    258   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    259   CHK(suvm_volume_ref_put(vol) == RES_OK);
    260 
    261   args.tetrahedron_data.size = 0;
    262   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    263   args.tetrahedron_data.size = sizeof(size_t[4]);
    264   args.tetrahedron_data.alignment = 0;
    265   msh.tetrahedron_data_alignment = 0;
    266   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    267   args.tetrahedron_data.alignment = 3;
    268   msh.tetrahedron_data_alignment = 3;
    269   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    270   args.tetrahedron_data.alignment = 64;
    271   msh.tetrahedron_data_alignment = 64;
    272   args.vertex_data.size = 0;
    273   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    274   args.vertex_data.size = sizeof(double[3]);
    275   args.vertex_data.alignment = 5;
    276   msh.vertex_data_alignment = 5;
    277   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    278   args.vertex_data.alignment = 0;
    279   msh.vertex_data_alignment = 0;
    280   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_BAD_ARG);
    281   args.vertex_data.alignment = 32;
    282   msh.vertex_data_alignment = 32;
    283   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    284 
    285   CHK(suvm_volume_get_primitives_count(NULL, &nprims) == RES_BAD_ARG);
    286   CHK(suvm_volume_get_primitives_count(vol, NULL) == RES_BAD_ARG);
    287   CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK);
    288   CHK(nprims == msh.ntetrahedra);
    289 
    290   CHK(suvm_volume_get_primitive(NULL, 0, &prim) == RES_BAD_ARG);
    291   CHK(suvm_volume_get_primitive(vol, nprims, &prim) == RES_BAD_ARG);
    292   CHK(suvm_volume_get_primitive(vol, 0, NULL) == RES_BAD_ARG);
    293   CHK(suvm_volume_get_primitive(vol, 0, &prim) == RES_OK);
    294 
    295   CHK(suvm_primitive_setup_polyhedron(NULL, &poly) == RES_BAD_ARG);
    296   CHK(suvm_primitive_setup_polyhedron(&prim, NULL) == RES_BAD_ARG);
    297   CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
    298 
    299   check_volume_polyhedra(vol, &msh);
    300 
    301   CHK(suvm_volume_get_mesh_desc(NULL, &desc) == RES_BAD_ARG);
    302   CHK(suvm_volume_get_mesh_desc(vol, NULL) == RES_BAD_ARG);
    303   CHK(suvm_volume_get_mesh_desc(vol, &desc) == RES_OK);
    304 
    305   check_volume_mesh(vol, &msh);
    306 
    307   CHK(suvm_volume_ref_put(vol) == RES_OK);
    308 }
    309 
    310 static void
    311 check_prim
    312   (const struct suvm_primitive* prim,
    313    const struct suvm_volume* vol,
    314    struct mesh* msh)
    315 {
    316   size_t ids[4];
    317   double verts[4][3];
    318   const size_t* data;
    319   const double* vert_data[4];
    320 
    321   /* Check primitive data */
    322   CHK(prim != NULL);
    323   CHK(vol != NULL);
    324   CHK(prim->volume__ == vol);
    325   CHK(!SUVM_PRIMITIVE_NONE(prim));
    326   CHK(prim->data != NULL);
    327   CHK(prim->data != NULL);
    328   CHK(prim->vertex_data[0] != NULL);
    329   CHK(prim->vertex_data[1] != NULL);
    330   CHK(prim->vertex_data[2] != NULL);
    331   CHK(prim->vertex_data[3] != NULL);
    332   CHK(prim->nvertices == 4);
    333 
    334   /* Fetch tetrahedron vertices */
    335   get_indices(prim->iprim, ids, msh);
    336   CHK(prim->indices[0] == ids[0]);
    337   CHK(prim->indices[1] == ids[1]);
    338   CHK(prim->indices[2] == ids[2]);
    339   CHK(prim->indices[3] == ids[3]);
    340   get_position(ids[0], verts[0], msh);
    341   get_position(ids[1], verts[1], msh);
    342   get_position(ids[2], verts[2], msh);
    343   get_position(ids[3], verts[3], msh);
    344 
    345   /* Check per primitive data */
    346   data = prim->data;
    347   CHK(data[0] == ids[0]);
    348   CHK(data[1] == ids[1]);
    349   CHK(data[2] == ids[2]);
    350   CHK(data[3] == ids[3]);
    351 
    352   /* Check per vertex data */
    353   vert_data[0] =  prim->vertex_data[0];
    354   vert_data[1] =  prim->vertex_data[1];
    355   vert_data[2] =  prim->vertex_data[2];
    356   vert_data[3] =  prim->vertex_data[3];
    357   CHK(d3_eq(vert_data[0], verts[0]));
    358   CHK(d3_eq(vert_data[1], verts[1]));
    359   CHK(d3_eq(vert_data[2], verts[2]));
    360   CHK(d3_eq(vert_data[3], verts[3]));
    361 }
    362 
    363 static void
    364 check_prim_at
    365   (const struct suvm_primitive* prim,
    366    const struct suvm_volume* vol,
    367    const double pos[3],
    368    const double bcoords[4],
    369    struct mesh* msh)
    370 {
    371   size_t ids[4];
    372   double verts[4][3];
    373   double E0[3], E1[3];
    374   double N[4][3];
    375   double dst[4];
    376   double p[3];
    377   const double* vert_data[4];
    378 
    379   check_prim(prim, vol, msh);
    380 
    381   /* Fetch tetrahedron vertices */
    382   get_indices(prim->iprim, ids, msh);
    383   get_position(ids[0], verts[0], msh);
    384   get_position(ids[1], verts[1], msh);
    385   get_position(ids[2], verts[2], msh);
    386   get_position(ids[3], verts[3], msh);
    387 
    388   /* Compute tetrahdron normals */
    389   d3_sub(E0, verts[1], verts[0]);
    390   d3_sub(E1, verts[2], verts[0]);
    391   d3_cross(N[0], E0, E1);
    392   d3_sub(E0, verts[1], verts[3]);
    393   d3_sub(E1, verts[0], verts[3]);
    394   d3_cross(N[1], E0, E1);
    395   d3_sub(E0, verts[2], verts[3]);
    396   d3_sub(E1, verts[1], verts[3]);
    397   d3_cross(N[2], E0, E1);
    398   d3_sub(E0, verts[0], verts[3]);
    399   d3_sub(E1, verts[2], verts[3]);
    400   d3_cross(N[3], E0, E1);
    401 
    402   /* Evaluate the distance from pos to the tetrahedron planes */
    403   dst[0] = d3_dot(N[0], d3_sub(p, pos, verts[0]));
    404   dst[1] = d3_dot(N[1], d3_sub(p, pos, verts[1]));
    405   dst[2] = d3_dot(N[2], d3_sub(p, pos, verts[2]));
    406   dst[3] = d3_dot(N[3], d3_sub(p, pos, verts[3]));
    407 
    408   /* Check that pos lies into the tetrahedron */
    409   CHK(dst[0] >= 0);
    410   CHK(dst[1] >= 0);
    411   CHK(dst[2] >= 0);
    412   CHK(dst[3] >= 0);
    413 
    414   vert_data[0] =  prim->vertex_data[0];
    415   vert_data[1] =  prim->vertex_data[1];
    416   vert_data[2] =  prim->vertex_data[2];
    417   vert_data[3] =  prim->vertex_data[3];
    418 
    419   /* Check interpolation parameter */
    420   p[0] =
    421     vert_data[0][0] * bcoords[0]
    422   + vert_data[1][0] * bcoords[1]
    423   + vert_data[2][0] * bcoords[2]
    424   + vert_data[3][0] * bcoords[3];
    425   p[1] =
    426     vert_data[0][1] * bcoords[0]
    427   + vert_data[1][1] * bcoords[1]
    428   + vert_data[2][1] * bcoords[2]
    429   + vert_data[3][1] * bcoords[3];
    430   p[2] =
    431     vert_data[0][2] * bcoords[0]
    432   + vert_data[1][2] * bcoords[1]
    433   + vert_data[2][2] * bcoords[2]
    434   + vert_data[3][2] * bcoords[3];
    435   CHK(d3_eq_eps(p, pos, 1.e-6));
    436 }
    437 
    438 static void
    439 check_prims_intersect_aabb
    440   (struct darray_prim* primitives,
    441    const struct suvm_volume* vol,
    442    const double low[3],
    443    const double upp[3],
    444    struct mesh* msh)
    445 {
    446   struct suvm_primitive* prims = NULL;
    447   size_t nprims;
    448   size_t iprim;
    449   size_t iprim_challenge;
    450   size_t ids[4];
    451   double vtx[4][3];
    452   double prim_low[3];
    453   double prim_upp[3];
    454   float lowf[3];
    455   float uppf[3];
    456 
    457   CHK(primitives);
    458   prims = darray_prim_data_get(primitives);
    459   nprims = darray_prim_size_get(primitives);
    460 
    461   FOR_EACH(iprim, 0, nprims) {
    462     check_prim(prims + iprim, vol, msh);
    463 
    464     /* Fetch tetrahedron vertices */
    465     get_indices(prims[iprim].iprim, ids, msh);
    466     get_position(ids[0], vtx[0], msh);
    467     get_position(ids[1], vtx[1], msh);
    468     get_position(ids[2], vtx[2], msh);
    469     get_position(ids[3], vtx[3], msh);
    470 
    471     /* Compute the primitive AABB */
    472     prim_low[0] = MMIN(MMIN(vtx[0][0], vtx[1][0]), MMIN(vtx[2][0], vtx[3][0]));
    473     prim_low[1] = MMIN(MMIN(vtx[0][1], vtx[1][1]), MMIN(vtx[2][1], vtx[3][1]));
    474     prim_low[2] = MMIN(MMIN(vtx[0][2], vtx[1][2]), MMIN(vtx[2][2], vtx[3][2]));
    475     prim_upp[0] = MMAX(MMAX(vtx[0][0], vtx[1][0]), MMAX(vtx[2][0], vtx[3][0]));
    476     prim_upp[1] = MMAX(MMAX(vtx[0][1], vtx[1][1]), MMAX(vtx[2][1], vtx[3][1]));
    477     prim_upp[2] = MMAX(MMAX(vtx[0][2], vtx[1][2]), MMAX(vtx[2][2], vtx[3][2]));
    478 
    479     /* Check AABB intersections */
    480     CHK(prim_low[0] < upp[0] && prim_upp[0] > low[0]);
    481     CHK(prim_low[1] < upp[1] && prim_upp[1] > low[1]);
    482     CHK(prim_low[2] < upp[2] && prim_upp[2] > low[2]);
    483   }
    484 
    485   /* Sort the primitives in ascending order wrt its identifier */
    486   qsort(prims, nprims, sizeof(*prims), cmp_prim);
    487 
    488   /* Exhaustively check all the mesh primitives against the AABB and ensure
    489    * that the input primitives are effectively the same of the ones detected by
    490    * this brute force procedure */
    491   iprim_challenge = 0;
    492   f3_set_d3(lowf, low);
    493   f3_set_d3(uppf, upp);
    494   FOR_EACH(iprim, 0, msh->ntetrahedra) {
    495     struct suvm_primitive prim;
    496     struct suvm_polyhedron tetra;
    497     enum suvm_intersection_type intersect;
    498 
    499     CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
    500     CHK(suvm_primitive_setup_polyhedron(&prim, &tetra) == RES_OK);
    501 
    502     /* Check that the primitive intersects the AABB */
    503     intersect = suvm_polyhedron_intersect_aabb(&tetra, lowf, uppf);
    504     if(intersect != SUVM_INTERSECT_NONE) {
    505       CHK(iprim == prims[iprim_challenge].iprim);
    506       ++iprim_challenge;
    507     }
    508   }
    509   /* All submitted primitives were exhaustively challenged */
    510   CHK(iprim_challenge == nprims);
    511 }
    512 
    513 static void
    514 check_hash
    515   (const struct suvm_volume* vol,
    516    const struct mesh* msh,
    517    const int has_prim_data,
    518    const int has_vert_data)
    519 {
    520   struct suvm_mesh_desc msh_desc;
    521   void* mem = NULL;
    522   float* pos = NULL;
    523   uint32_t* ids = NULL;
    524   void* prims = NULL;
    525   void* verts = NULL;
    526   size_t sz_pos = 0;
    527   size_t sz_ids = 0;
    528   size_t sz_prims = 0;
    529   size_t sz_verts = 0;
    530   size_t i;
    531   hash256_T hash0, hash1;
    532 
    533   CHK(suvm_volume_compute_hash(NULL, 0, hash0) == RES_BAD_ARG);
    534   CHK(suvm_volume_compute_hash(vol, 0, NULL) == RES_BAD_ARG);
    535   CHK(suvm_volume_compute_hash(vol, 0, hash0) == RES_OK);
    536 
    537   hash_sha256(NULL, 0, hash1);
    538   CHK(hash256_eq(hash0, hash1));
    539 
    540   msh_desc = SUVM_MESH_DESC_NULL;
    541   CHK(suvm_mesh_desc_compute_hash(NULL, hash0) == RES_BAD_ARG);
    542   CHK(suvm_mesh_desc_compute_hash(&msh_desc, NULL) == RES_BAD_ARG);
    543   CHK(suvm_mesh_desc_compute_hash(&msh_desc, hash0) == RES_OK);
    544 
    545   CHK(suvm_volume_get_mesh_desc(vol, &msh_desc) == RES_OK);
    546   CHK(suvm_mesh_desc_compute_hash(&msh_desc, hash1) == RES_OK);
    547   CHK(!hash256_eq(hash0, hash1));
    548 
    549   /* Compute data size to hash */
    550   sz_pos = msh->nvertices*sizeof(float[3]);
    551   sz_ids = msh->ntetrahedra*sizeof(uint32_t[4]);
    552   if(has_prim_data) {
    553      sz_prims = msh->ntetrahedra*sizeof(size_t[4]);
    554   }
    555   if(has_vert_data) {
    556     sz_verts = msh->nvertices*sizeof(double[3]);
    557   }
    558   mem = mem_calloc(1, sz_pos + sz_ids + sz_prims + sz_verts);
    559   CHK(mem != NULL);
    560 
    561   /* Copy data to hash into the allocated memory block. Be carefull the memory
    562    * layout used by SUVM to hash the volume. First the vertices, then the
    563    * indices followed by the per prim data and finally the per vertex data */
    564 
    565   /* SUVM stores the vertices in single precision. Convert vertex coordinates
    566    * in float to ensure bit correspondance */
    567   pos = mem;
    568   FOR_EACH(i, 0, msh->nvertices) {
    569     pos[i*3 + 0] = (float)msh->vertices[i*3+0];
    570     pos[i*3 + 1] = (float)msh->vertices[i*3+1];
    571     pos[i*3 + 2] = (float)msh->vertices[i*3+2];
    572   }
    573 
    574   /* SUVM stores the tetrahedron indices as 32 bits unsigned integers. Convert
    575    * indices to ensure bit correspondance */
    576   ids = (uint32_t*)((char*)mem + sz_pos);
    577   FOR_EACH(i, 0, msh->ntetrahedra) {
    578     ids[i*4 + 0] = (uint32_t)msh->tetrahedra[i*4+0];
    579     ids[i*4 + 1] = (uint32_t)msh->tetrahedra[i*4+1];
    580     ids[i*4 + 2] = (uint32_t)msh->tetrahedra[i*4+2];
    581     ids[i*4 + 3] = (uint32_t)msh->tetrahedra[i*4+3];
    582   }
    583 
    584   /* Copy per primitive data */
    585   if(has_prim_data) {
    586     prims = ((char*)mem + sz_pos + sz_ids);
    587     memcpy(prims, msh->tetrahedra, msh->ntetrahedra*sizeof(size_t[4]));
    588   }
    589 
    590   /* Copy per vertex data */
    591   if(has_vert_data) {
    592     verts = ((char*)mem + sz_pos + sz_ids + sz_prims);
    593     memcpy(verts, msh->vertices, msh->nvertices*sizeof(double[3]));
    594   }
    595 
    596   CHK(suvm_volume_compute_hash(vol, SUVM_POSITIONS, hash0) == RES_OK);
    597   hash_sha256(pos, sz_pos, hash1);
    598   CHK(hash256_eq(hash0, hash1));
    599 
    600   CHK(suvm_volume_compute_hash(vol, SUVM_INDICES, hash0) == RES_OK);
    601   hash_sha256(ids, sz_ids, hash1);
    602   CHK(hash256_eq(hash0, hash1));
    603 
    604   CHK(suvm_volume_compute_hash(vol, SUVM_PRIMITIVE_DATA, hash0) == RES_OK);
    605   hash_sha256(prims, sz_prims, hash1);
    606   CHK(hash256_eq(hash0, hash1));
    607 
    608   CHK(suvm_volume_compute_hash(vol, SUVM_VERTEX_DATA, hash0) == RES_OK);
    609   hash_sha256(verts, sz_verts, hash1);
    610   CHK(hash256_eq(hash0, hash1));
    611 
    612   CHK(suvm_volume_compute_hash(vol, SUVM_POSITIONS|SUVM_INDICES, hash0) == RES_OK);
    613   hash_sha256(mem, sz_pos + sz_ids, hash1);
    614   CHK(hash256_eq(hash0, hash1));
    615 
    616   CHK(suvm_volume_compute_hash
    617     (vol, SUVM_POSITIONS|SUVM_INDICES|SUVM_PRIMITIVE_DATA, hash0) == RES_OK);
    618   hash_sha256(mem, sz_pos + sz_ids + sz_prims, hash1);
    619   CHK(hash256_eq(hash0, hash1));
    620 
    621   CHK(suvm_volume_compute_hash(vol,
    622      SUVM_POSITIONS|SUVM_INDICES|SUVM_PRIMITIVE_DATA|SUVM_VERTEX_DATA,
    623      hash0) == RES_OK);
    624   hash_sha256(mem, sz_pos + sz_ids + sz_prims + sz_verts, hash1);
    625   CHK(hash256_eq(hash0, hash1));
    626 
    627   mem_rm(mem);
    628 }
    629 
    630 static void
    631 test_volume_hash(struct suvm_device* dev, struct mesh* msh)
    632 {
    633   struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
    634   struct suvm_volume* vol = NULL;
    635 
    636   args.ntetrahedra = msh->ntetrahedra;
    637   args.nvertices = msh->nvertices;
    638   args.get_indices = get_indices;
    639   args.get_position = get_position;
    640   args.context = msh;
    641   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    642 
    643   check_hash(vol, msh, 0, 0);
    644   CHK(suvm_volume_ref_put(vol) == RES_OK);
    645 
    646   args.vertex_data.get = get_vert_data;
    647   args.vertex_data.size = sizeof(double[3]);
    648   args.vertex_data.alignment = ALIGNOF(double[3]);
    649   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    650 
    651   check_hash(vol, msh, 0, 1);
    652   CHK(suvm_volume_ref_put(vol) == RES_OK);
    653 
    654   args.tetrahedron_data.get = get_tetra_data;
    655   args.tetrahedron_data.size = sizeof(size_t[4]);
    656   args.tetrahedron_data.alignment = ALIGNOF(size_t[4]);
    657   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    658   check_hash(vol, msh, 1, 1);
    659 
    660   CHK(suvm_volume_ref_put(vol) == RES_OK);
    661 }
    662 
    663 static void
    664 prim_intersect_aabb
    665   (const struct suvm_primitive* prim,
    666    const double low[3],
    667    const double upp[3],
    668    void* context)
    669 {
    670   struct darray_prim* prims = context;
    671   CHK(prim && low && upp && context);
    672   CHK(low[0] < upp[0]);
    673   CHK(low[1] < upp[1]);
    674   CHK(low[2] < upp[2]);
    675   CHK(darray_prim_push_back(prims, prim) == RES_OK);
    676 }
    677 
    678 static void
    679 test_volume_at_cube(struct suvm_device* dev)
    680 {
    681   struct darray_prim prims;
    682   struct mesh msh = MESH_NULL;
    683   struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
    684   struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
    685   struct suvm_volume* vol = NULL;
    686   double bcoords[4];
    687   double pos[3];
    688   double tmp[3];
    689   double low[3], upp[3];
    690   const size_t nsamps = 100;
    691   size_t i;
    692 
    693   darray_prim_init(NULL, &prims);
    694 
    695   args.ntetrahedra = box_ntetras;
    696   args.nvertices = box_nverts;
    697   args.get_indices = get_indices;
    698   args.get_position = get_position;
    699   args.tetrahedron_data.get = get_tetra_data;
    700   args.tetrahedron_data.size = sizeof(size_t[4]);
    701   args.tetrahedron_data.alignment = ALIGNOF(size_t[4]);
    702   args.vertex_data.get = get_vert_data;
    703   args.vertex_data.size = sizeof(double[3]);
    704   args.vertex_data.alignment = ALIGNOF(double[3]);
    705   args.context = &msh;
    706 
    707   msh.vertices = box_vertices;
    708   msh.nvertices = box_nverts;
    709   msh.tetrahedra = box_indices;
    710   msh.ntetrahedra = box_ntetras;
    711   msh.tetrahedron_data_alignment = args.tetrahedron_data.alignment;
    712   msh.vertex_data_alignment = args.vertex_data.alignment;
    713 
    714   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    715   check_volume_mesh(vol, &msh);
    716 
    717   pos[0] = 0.25;
    718   pos[1] = 0.4;
    719   pos[2] = 0.4;
    720 
    721   CHK(suvm_volume_at(NULL, pos, &prim, bcoords) == RES_BAD_ARG);
    722   CHK(suvm_volume_at(vol, NULL, &prim, bcoords) == RES_BAD_ARG);
    723   CHK(suvm_volume_at(vol, pos, NULL, bcoords) == RES_BAD_ARG);
    724   CHK(suvm_volume_at(vol, pos, &prim, NULL) == RES_BAD_ARG);
    725   CHK(suvm_volume_at(vol, pos, &prim, bcoords) == RES_OK);
    726   CHK(!SUVM_PRIMITIVE_NONE(&prim));
    727   CHK(prim.iprim == 2);
    728   check_prim_at(&prim, vol, pos, bcoords, &msh);
    729 
    730   CHK(suvm_volume_get_aabb(NULL, low, upp) == RES_BAD_ARG);
    731   CHK(suvm_volume_get_aabb(vol, NULL, upp) == RES_BAD_ARG);
    732   CHK(suvm_volume_get_aabb(vol, low, NULL) == RES_BAD_ARG);
    733   CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK);
    734   CHK(d3_eq(low, d3_splat(tmp, 0)));
    735   CHK(d3_eq(upp, d3_splat(tmp, 1)));
    736   FOR_EACH(i, 0, nsamps) {
    737     double samp[3];
    738     samp[0] = low[0] + rand_canonic()*(upp[0] - low[0]);
    739     samp[1] = low[1] + rand_canonic()*(upp[1] - low[1]);
    740     samp[2] = low[2] + rand_canonic()*(upp[2] - low[2]);
    741     CHK(suvm_volume_at(vol, samp, &prim, bcoords) == RES_OK);
    742     check_prim_at(&prim, vol, samp, bcoords, &msh);
    743   }
    744 
    745   pos[0] = upp[0] + 0.1;
    746   pos[1] = upp[1] + 0.1;
    747   pos[2] = upp[2] + 0.1;
    748   CHK(suvm_volume_at(vol, pos, &prim, bcoords) == RES_OK);
    749   CHK(SUVM_PRIMITIVE_NONE(&prim));
    750 
    751   d3_splat(low, 0.25);
    752   d3_splat(upp, 0.75);
    753 
    754   CHK(suvm_volume_intersect_aabb
    755     (NULL, low, upp, prim_intersect_aabb, &prims) == RES_BAD_ARG);
    756   CHK(suvm_volume_intersect_aabb
    757     (vol, NULL, upp, prim_intersect_aabb, &prims) == RES_BAD_ARG);
    758   CHK(suvm_volume_intersect_aabb
    759     (vol, low, NULL, prim_intersect_aabb, &prims) == RES_BAD_ARG);
    760   CHK(suvm_volume_intersect_aabb
    761     (vol, low, upp, NULL, &prims) == RES_BAD_ARG);
    762   CHK(suvm_volume_intersect_aabb
    763     (vol, low, upp, prim_intersect_aabb, &prims) == RES_OK);
    764 
    765   CHK(darray_prim_size_get(&prims) == box_ntetras);
    766   check_prims_intersect_aabb(&prims, vol, low, upp, &msh);
    767 
    768   CHK(suvm_volume_intersect_aabb
    769     (vol, upp, low, prim_intersect_aabb, &prims) == RES_BAD_ARG);
    770   CHK(suvm_volume_intersect_aabb
    771     (vol, low, low, prim_intersect_aabb, &prims) == RES_BAD_ARG);
    772 
    773   darray_prim_clear(&prims);
    774   d3(low,-1, 0, 0);
    775   d3(upp, 0, 1, 1);
    776   CHK(suvm_volume_intersect_aabb
    777     (vol, low, upp, prim_intersect_aabb, &prims) == RES_OK);
    778   CHK(darray_prim_size_get(&prims) == 0);
    779 
    780   d3(low, 1, 0, 0);
    781   d3(upp, 2, 1, 1);
    782   CHK(suvm_volume_intersect_aabb
    783     (vol, low, upp, prim_intersect_aabb, &prims) == RES_OK);
    784   CHK(darray_prim_size_get(&prims) == 0);
    785 
    786   d3(low, 0.9, 0, 0);
    787   d3(upp, 2, 1, 1);
    788   CHK(suvm_volume_intersect_aabb
    789     (vol, low, upp, prim_intersect_aabb, &prims) == RES_OK);
    790   CHK(darray_prim_size_get(&prims) == 10);
    791   check_prims_intersect_aabb(&prims, vol, low, upp, &msh);
    792 
    793   CHK(suvm_volume_ref_put(vol) == RES_OK);
    794   darray_prim_release(&prims);
    795 }
    796 
    797 static void
    798 test_volume_hash_cube(struct suvm_device* dev)
    799 {
    800   struct mesh msh = MESH_NULL;
    801 
    802   msh.vertices = box_vertices;
    803   msh.nvertices = box_nverts;
    804   msh.tetrahedra = box_indices;
    805   msh.ntetrahedra = box_ntetras;
    806   msh.tetrahedron_data_alignment = ALIGNOF(size_t[4]);
    807   msh.vertex_data_alignment = ALIGNOF(double[3]);
    808   test_volume_hash(dev, &msh);
    809 }
    810 
    811 static void
    812 test_volume_at_ball(struct suvm_device* dev)
    813 {
    814   struct darray_prim prims;
    815   struct mesh msh = MESH_NULL;
    816   struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
    817   struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
    818   struct suvm_volume* vol = NULL;
    819   double bcoords[4];
    820   double aabb_low[3], aabb_upp[3], aabb_sz;
    821   double low[3], upp[3];
    822   double vec[3];
    823   size_t ivert;
    824   const size_t nsamps = 10000;
    825   const size_t naabbs = 100;
    826   size_t i;
    827   int check_at = 0;
    828   int check_intersect_aabb = 0;
    829 
    830   darray_prim_init(NULL, &prims);
    831 
    832   args.precompute_normals = 1;
    833   args.ntetrahedra = ball_ntetrahedra;
    834   args.nvertices = ball_nvertices;
    835   args.get_indices = get_indices;
    836   args.get_position = get_position;
    837   args.tetrahedron_data.get = get_tetra_data;
    838   args.tetrahedron_data.size = sizeof(size_t[4]);
    839   args.tetrahedron_data.alignment = ALIGNOF(size_t[4]);
    840   args.vertex_data.get = get_vert_data;
    841   args.vertex_data.size = sizeof(double[3]);
    842   args.vertex_data.alignment = ALIGNOF(double[3]);
    843   args.context = &msh;
    844 
    845   msh.vertices = ball_vertices;
    846   msh.nvertices = ball_nvertices;
    847   msh.tetrahedra = ball_tetrahedra;
    848   msh.ntetrahedra = ball_ntetrahedra;
    849   msh.tetrahedron_data_alignment = args.tetrahedron_data.alignment;
    850   msh.vertex_data_alignment = args.vertex_data.alignment;
    851 
    852   /* Compute the ball aabb */
    853   d3_splat(aabb_low, DBL_MAX);
    854   d3_splat(aabb_upp,-DBL_MAX);
    855   FOR_EACH(ivert, 0, ball_nvertices) {
    856     aabb_low[0] = MMIN(aabb_low[0], ball_vertices[ivert*3+0]);
    857     aabb_low[1] = MMIN(aabb_low[1], ball_vertices[ivert*3+1]);
    858     aabb_low[2] = MMIN(aabb_low[2], ball_vertices[ivert*3+2]);
    859     aabb_upp[0] = MMAX(aabb_upp[0], ball_vertices[ivert*3+0]);
    860     aabb_upp[1] = MMAX(aabb_upp[1], ball_vertices[ivert*3+1]);
    861     aabb_upp[2] = MMAX(aabb_upp[2], ball_vertices[ivert*3+2]);
    862   }
    863   aabb_sz = d3_len(d3_sub(vec, aabb_upp, aabb_low));
    864 
    865   CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
    866   check_volume_mesh(vol, &msh);
    867 
    868   check_volume_polyhedra(vol, &msh);
    869 
    870   /* Check volume AABB */
    871   CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK);
    872   CHK(d3_eq_eps(aabb_low, low, aabb_sz*1.e-6));
    873   CHK(d3_eq_eps(aabb_upp, upp, aabb_sz*1.e-6));
    874 
    875   /* Check at operator */
    876   check_at = 0;
    877   FOR_EACH(i, 0, nsamps) {
    878     double samp[3];
    879     samp[0] = low[0] + rand_canonic()*(upp[0] - low[0]);
    880     samp[1] = low[1] + rand_canonic()*(upp[1] - low[1]);
    881     samp[2] = low[2] + rand_canonic()*(upp[2] - low[2]);
    882     CHK(suvm_volume_at(vol, samp, &prim, bcoords) == RES_OK);
    883 
    884     if(!SUVM_PRIMITIVE_NONE(&prim)) {
    885       check_at += 1;
    886       check_prim_at(&prim, vol, samp, bcoords, &msh);
    887     }
    888   }
    889   CHK(check_at != 0);
    890 
    891   /* Check intersect AABB */
    892   check_intersect_aabb = 0;
    893   FOR_EACH(i, 0, naabbs) {
    894     double pt0[3];
    895     double pt1[3];
    896 
    897     /* Sample 2 random points into the mesh AABB */
    898     pt0[0] = low[0] + rand_canonic()*(upp[0] - low[0]);
    899     pt0[1] = low[1] + rand_canonic()*(upp[1] - low[1]);
    900     pt0[2] = low[2] + rand_canonic()*(upp[2] - low[2]);
    901     pt1[0] = low[0] + rand_canonic()*(upp[0] - low[0]);
    902     pt1[1] = low[1] + rand_canonic()*(upp[1] - low[1]);
    903     pt1[2] = low[2] + rand_canonic()*(upp[2] - low[2]);
    904 
    905     /* Build the bounding box of the sampled point */
    906     aabb_low[0] = MMIN(pt0[0], pt1[0]);
    907     aabb_low[1] = MMIN(pt0[1], pt1[1]);
    908     aabb_low[2] = MMIN(pt0[2], pt1[2]);
    909     aabb_upp[0] = MMAX(pt0[0], pt1[0]);
    910     aabb_upp[1] = MMAX(pt0[1], pt1[1]);
    911     aabb_upp[2] = MMAX(pt0[2], pt1[2]);
    912 
    913     CHK(suvm_volume_intersect_aabb
    914       (vol, aabb_low, aabb_upp, prim_intersect_aabb, &prims) == RES_OK);
    915 
    916     check_intersect_aabb = darray_prim_size_get(&prims) != 0;
    917     check_prims_intersect_aabb(&prims, vol, aabb_low, aabb_upp, &msh);
    918     darray_prim_clear(&prims);
    919   }
    920   CHK(check_intersect_aabb != 0);
    921 
    922   CHK(suvm_volume_ref_put(vol) == RES_OK);
    923   darray_prim_release(&prims);
    924 }
    925 
    926 static void
    927 test_volume_hash_ball(struct suvm_device* dev)
    928 {
    929   struct mesh msh = MESH_NULL;
    930 
    931   msh.vertices = ball_vertices;
    932   msh.nvertices = ball_nvertices;
    933   msh.tetrahedra = ball_tetrahedra;
    934   msh.ntetrahedra = ball_ntetrahedra;
    935   msh.tetrahedron_data_alignment = ALIGNOF(size_t[4]);
    936   msh.vertex_data_alignment = ALIGNOF(double[3]);
    937   test_volume_hash(dev, &msh);
    938 }
    939 
    940 /*******************************************************************************
    941  * Main function
    942  ******************************************************************************/
    943 int
    944 main(int argc, char** argv)
    945 {
    946   struct suvm_device* dev = NULL;
    947   (void)argc, (void)argv;
    948 
    949   CHK(suvm_device_create(NULL, &mem_default_allocator, 1, &dev) == RES_OK);
    950 
    951   test_tetrahedral_mesh_creation(dev);
    952   test_volume_at_cube(dev);
    953   test_volume_at_ball(dev);
    954   test_volume_hash_cube(dev);
    955   test_volume_hash_ball(dev);
    956 
    957   CHK(suvm_device_ref_put(dev) == RES_OK);
    958 
    959   check_memory_allocator(&mem_default_allocator);
    960   CHK(mem_allocated_size() == 0);
    961   return 0;
    962 }
    963