commit 638f2fff672dc9e39577379ecada3612975a0dc1
parent 53552af739157da351f96a3a8eb1b4d2688cc90e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 12 Jan 2021 11:41:25 +0100
Begin the implementation of the "build octree" process
Diffstat:
7 files changed, 654 insertions(+), 20 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -26,9 +26,10 @@ option(NO_TEST "Do not build tests" OFF)
find_package(AtrTP 0.0 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
find_package(RCMake 0.4 REQUIRED)
-find_package(RSys 0.10 REQUIRED)
+find_package(RSys 0.12 REQUIRED)
find_package(StarTetraHedra 0.0 REQUIRED)
find_package(StarUVM 0.0 REQUIRED)
+find_package(StarVX 0.2 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
@@ -50,10 +51,12 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(ATRSTM_FILES_SRC
atrstm.c
+ atrstm_build_octrees.c
atrstm_log.c
atrstm_partition.c
atrstm_setup_uvm.c)
set(ATRSTM_FILES_INC
+ atrstm_build_octrees.h
atrstm_c.h
atrstm_log.h
atrstm_partition.h)
@@ -69,7 +72,7 @@ rcmake_prepend_path(ATRSTM_FILES_INC_API ${ATRSTM_SOURCE_DIR})
rcmake_prepend_path(ATRSTM_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(atrstm SHARED ${ATRSTM_FILES_SRC} ${ATRSTM_FILES_INC} ${ATRSTM_FILES_INC_API})
-target_link_libraries(atrstm AtrTP RSys StarTetraHedra StarUVM)
+target_link_libraries(atrstm AtrTP RSys StarTetraHedra StarUVM StarVX)
if(CMAKE_COMPILER_IS_GNUCC)
set_target_properties(atrstm PROPERTIES LINK_FLAGS "${OpenMP_C_FLAGS}")
diff --git a/src/atrstm.c b/src/atrstm.c
@@ -21,6 +21,7 @@
#include <rsys/mem_allocator.h>
#include <star/suvm.h>
+#include <star/svx.h>
#include <omp.h>
@@ -53,6 +54,7 @@ release_atrstm(ref_T* ref)
atrstm = CONTAINER_OF(ref, struct atrstm, ref);
if(atrstm->atrtp) ATRTP(ref_put(atrstm->atrtp));
if(atrstm->suvm) SUVM(device_ref_put(atrstm->suvm));
+ if(atrstm->svx) SVX(device_ref_put(atrstm->svx));
if(atrstm->logger == &atrstm->logger__) logger_release(&atrstm->logger__);
str_release(&atrstm->name);
MEM_RM(atrstm->allocator, atrstm);
@@ -138,6 +140,10 @@ atrstm_create
res = atrtp_load(atrstm->atrtp, args->atrtp_filename);
if(res != RES_OK) goto error;
+ /* Create the Star-VoXel device */
+ res = svx_device_create
+ (atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->svx);
+
exit:
if(out_atrstm) *out_atrstm = atrstm;
return res;
diff --git a/src/atrstm_build_octrees.c b/src/atrstm_build_octrees.c
@@ -0,0 +1,528 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "atrstm_c.h"
+#include "atrstm_build_octrees.h"
+#include "atrstm_log.h"
+#include "atrstm_svx.h"
+#include "atrstm_partition.h"
+
+#include <star/suvm.h>
+#include <star/svx.h>
+
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_size_t.h>
+#include <rsys/morton.h>
+
+#include <omp.h>
+
+struct build_octree_context {
+ struct atrstm* atrstm;
+ struct pool* pools; /* Per thread pool of voxel partitions */
+
+ /* Optical thickness threshold criteria for the merge process */
+ double tau_threshold;
+};
+#define BUILD_OCTREE_CONTEXT_NULL__ { NULL, NULL, 0 }
+static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
+ BUILD_OCTREE_CONTEXT_NULL__;
+
+/* Define the darray of list of primitive ids */
+#define DARRAY_NAME prims_list
+#define DARRAY_DATA struct darray_size_t
+#define DARRAY_FUNCTOR_INIT darray_size_t_init
+#define DARRAY_FUNCTOR_RELEASE darray_size_t_release
+#define DARRAY_FUNCTOR_COPY darray_size_t_copy
+#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_size_t_copy_and_release
+#include <rsys/dynamic_array.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+check_build_octrees_args(const struct build_octrees_args* args)
+{
+ return args
+ && args->grid_max_definition[0]
+ && args->grid_max_definition[1]
+ && args->grid_max_definition[2]
+ && args->optical_thickness_threshold
+ && (!args->cache.stream || args->cache.name);
+}
+
+static INLINE void
+register_primitive
+ (const struct suvm_primitive* prim,
+ const double low[3],
+ const double upp[3],
+ void* context)
+{
+ struct darray_size_t* prims = context;
+ ASSERT(prim && low && upp && context);
+ ASSERT(low[0] < upp[0]);
+ ASSERT(low[1] < upp[1]);
+ ASSERT(low[2] < upp[2]);
+ (void)low, (void)upp;
+ CHK(darray_size_t_push_back(prims, &prim->iprim) == RES_OK);
+}
+
+static res_T
+voxelize_partition
+ (struct atrstm* atrstm,
+ const struct darray_size_t* prims, /* Primitives intersecting the partition */
+ const double part_low[3], /* Lower bound of the partition */
+ const double part_upp[3], /* Upper bound of the partition */
+ const double vxsz[3], /* Size of a partition voxel */
+ struct part* part)
+{
+ size_t iprim;
+ res_T res = RES_OK;
+ ASSERT(atrstm && prims && part_low && part_upp && vxsz && part);
+ ASSERT(part_low[0] < part_upp[0]);
+ ASSERT(part_low[1] < part_upp[1]);
+ ASSERT(part_low[2] < part_upp[2]);
+ ASSERT(vxsz[0] > 0 && vxsz[1] > 0 && vxsz[2] > 0);
+
+ part_clear_voxels(part);
+
+ FOR_EACH(iprim, 0, darray_size_t_size_get(prims)) {
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ double poly_low[3];
+ double poly_upp[3];
+ uint32_t ivxl_low[3];
+ uint32_t ivxl_upp[3];
+ uint32_t ivxl[3];
+
+ /* Retrieve the primitive data and setup its polyhedron */
+ SUVM(volume_get_primitive(atrstm->volume, iprim, &prim));
+ SUVM(primitive_setup_polyhedron(&prim, &poly));
+
+ /* Clamp the poly AABB to the partition boundaries */
+ poly_low[0] = MMAX(poly.lower[0], part_low[0]);
+ poly_low[1] = MMAX(poly.lower[1], part_low[1]);
+ poly_low[2] = MMAX(poly.lower[2], part_low[2]);
+ poly_upp[0] = MMIN(poly.upper[0], part_upp[0]);
+ poly_upp[1] = MMIN(poly.upper[1], part_upp[1]);
+ poly_upp[2] = MMIN(poly.upper[2], part_upp[2]);
+
+ /* Transform the polyhedron AABB in partition voxel space */
+ ivxl_low[0] = (uint32_t)((poly_low[0] - part_low[0]) / vxsz[0]);
+ ivxl_low[1] = (uint32_t)((poly_low[1] - part_low[1]) / vxsz[1]);
+ ivxl_low[2] = (uint32_t)((poly_low[2] - part_low[2]) / vxsz[2]);
+ ivxl_upp[0] = (uint32_t)ceil((poly_upp[0] - part_low[0]) / vxsz[0]);
+ ivxl_upp[1] = (uint32_t)ceil((poly_upp[1] - part_low[1]) / vxsz[1]);
+ ivxl_upp[2] = (uint32_t)ceil((poly_upp[2] - part_low[2]) / vxsz[2]);
+ ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
+ ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
+ ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
+
+ /* Iterate over the partition voxels intersected by the primitive AABB */
+ FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) {
+ float vxl_low[3]; /* Voxel lower bound */
+ float vxl_upp[3]; /* Voxel upper bound */
+ uint64_t mcode[3]; /* Morton code cache */
+
+ vxl_low[2] = (float)((double)ivxl[2] * vxsz[2] + part_low[2]);
+ vxl_upp[2] = vxl_low[2] + (float)vxsz[2];
+ mcode[2] = morton3D_encode_u21(ivxl[2]) << 0;
+
+ FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) {
+ vxl_low[1] = (float)((double)ivxl[1] * vxsz[1] + part_low[1]);
+ vxl_upp[1] = vxl_low[1] + (float)vxsz[1];
+ mcode[1] = (morton3D_encode_u21(ivxl[1]) << 1) | mcode[2];
+
+ FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) {
+ enum suvm_intersection_type intersect;
+ vxl_low[0] = (float)((double)ivxl[0] * vxsz[0] + part_low[0]);
+ vxl_upp[0] = vxl_low[0] + (float)vxsz[0];
+
+ /* NOTE mcode[0] store the morton index of the voxel */
+ mcode[0] = (morton3D_encode_u21(ivxl[0]) << 2) | mcode[1];
+
+ intersect = suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
+ if(intersect != SUVM_INTERSECT_NONE) {
+ float* vox = part_get_voxel(part, mcode[2]);
+ /* TODO Register the optical properties of the primitive against
+ * the current voxel. */
+ (void)vox;
+ }
+ }
+ }
+ }
+ }
+
+ return res;
+}
+
+static res_T
+voxelize_volumetric_mesh
+ (struct atrstm* atrstm,
+ const struct build_octrees_args* args,
+ struct pool* pools)
+{
+ struct darray_prims_list prims_list; /* Per thread list of primitives */
+ const size_t part_def = PARTITION_DEFINITION;
+ size_t nparts[3]; /* #partitions along the 3 axis */
+ size_t nparts_adjusted;
+ double vol_low[3];
+ double vol_upp[3];
+ double vxsz[3];
+ uint64_t ipart;
+ ATOMIC res = RES_OK;
+ ASSERT(atrstm && args && pools);
+
+ darray_prims_list_init(atrstm->allocator, &prims_list);
+
+ /* Allocate the per thread primitive lists */
+ res = darray_prims_list_resize(&prims_list, atrstm->nthreads);
+ if(res != RES_OK) goto error;
+
+ /* Fetch the volume AABB and compute the size of one voxel */
+ SUVM(volume_get_aabb(atrstm->volume, vol_low, vol_upp));
+ vxsz[0] = (vol_upp[0] - vol_low[0]) / (double)args->grid_max_definition[0];
+ vxsz[1] = (vol_upp[1] - vol_low[1]) / (double)args->grid_max_definition[1];
+ vxsz[2] = (vol_upp[2] - vol_low[2]) / (double)args->grid_max_definition[2];
+
+ /* Compute the number of partitions */
+ nparts[0] = (args->grid_max_definition[0] + (part_def-1)/*ceil*/) / part_def;
+ nparts[1] = (args->grid_max_definition[1] + (part_def-1)/*ceil*/) / part_def;
+ nparts[2] = (args->grid_max_definition[2] + (part_def-1)/*ceil*/) / part_def;
+
+ /* Compute the overall #partitions allowing Morton indexing */
+ nparts_adjusted = MMAX(nparts[0], MMAX(nparts[1], nparts[2]));
+ nparts_adjusted = round_up_pow2(nparts_adjusted);
+ nparts_adjusted = nparts_adjusted*nparts_adjusted*nparts_adjusted;
+
+ /* Iterate over the partition according to their Morton order and voxelize
+ * the primitives that intersect it */
+ omp_set_num_threads((int)atrstm->nthreads);
+ #pragma omp parallel for schedule(static, 1/*chunk size*/)
+ for(ipart = 0; ipart < nparts_adjusted; ++ipart) {
+ struct darray_size_t* prims = NULL;
+ struct pool* pool = NULL;
+ double part_low[3];
+ double part_upp[3];
+ uint32_t part_ids[3];
+ const int ithread = omp_get_thread_num();
+ struct part* part = NULL;
+ res_T res_local = RES_OK;
+
+ /* Handle error */
+ if(ATOMIC_GET(&res) != RES_OK) continue;
+
+ /* Fetch the per thread list of primitives and partition pool */
+ prims = darray_prims_list_data_get(&prims_list)+ithread;
+ darray_size_t_clear(prims);
+ pool = pools + ithread;
+
+ morton_xyz_decode_u21(ipart, part_ids);
+
+ /* Check that the partition is not out of bound due to Morton indexing */
+ if(part_ids[0] >= nparts[0]
+ || part_ids[1] >= nparts[1]
+ || part_ids[2] >= nparts[2])
+ continue;
+
+ /* Compute the partition AABB */
+ part_low[0] = (double)(part_ids[0] * part_def) * vxsz[0] + vol_low[0];
+ part_low[1] = (double)(part_ids[1] * part_def) * vxsz[1] + vol_low[1];
+ part_low[2] = (double)(part_ids[2] * part_def) * vxsz[2] + vol_low[2];
+ part_upp[0] = part_low[0] + (double)part_def * vxsz[0];
+ part_upp[1] = part_low[1] + (double)part_def * vxsz[1];
+ part_upp[2] = part_low[2] + (double)part_def * vxsz[2];
+
+ /* Find the list of primitives that intersects the partition */
+ res_local = suvm_volume_intersect_aabb(atrstm->volume, part_low, part_upp,
+ register_primitive, prims);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ part = pool_new_partition(pool, ipart);
+ if(!part) { /* An error occurs */
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ res = voxelize_partition(atrstm, prims, part_low, part_upp, vxsz, part);
+ if(res != RES_OK) {
+ pool_free_partition(pool, part);
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ pool_commit_partition(pool, part);
+ }
+
+exit:
+ darray_prims_list_release(&prims_list);
+ return (res_T)res;
+error:
+ goto exit;
+}
+
+static void
+voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
+{
+ const struct build_octree_context* ctx = context;
+ struct pool* pool = NULL;
+ struct part* part = NULL;
+ float* vxl = NULL;
+ uint64_t ivxl, ipart, ipool;
+ ASSERT(xyz && dst && ctx);
+ (void)xyz;
+
+ /* Retrieve the partition morton id and the per partition voxel morton id
+ * from the overall voxel morton code. As voxels, partitions are sorted in
+ * morton order and thus the partition morton ID is encoded in the MSB of the
+ * morton code while the voxels morton ID is stored in its LSB. */
+ ipart = (mcode >> (LOG2_PARTITION_DEFINITION*3));
+ ivxl = (mcode & BIT_U64(LOG2_PARTITION_DEFINITION*3));
+
+ /* Compute the pool index containing the partition. Partitions are
+ * alternatively stored into the per thread pool. Consequentlu the i^th
+ * partition is stored in the (i % #thread)^th pool. */
+ ipool = ipart / ctx->atrstm->nthreads;
+
+ /* Fetch the pool storing the partition */
+ pool = ctx->pools + ipool;
+
+ /* Fetch the partition storing the voxel */
+ part = pool_fetch_partition(pool, ipart);
+
+ if(!part) { /* An error occurs */
+ memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
+ } else {
+ vxl = part_get_voxel(part, ivxl); /* Fetch the voxel data */
+
+ /* Copy voxel data */
+ memcpy(dst, vxl, NFLOATS_PER_VOXEL * sizeof(float));
+
+ /* Free the partition if the fetch voxel was its last one */
+ if(ivxl == PARTITION_NVOXELS - 1) {
+ pool_free_partition(pool, part);
+ }
+ }
+}
+
+static void
+voxels_merge_component
+ (void* dst,
+ const enum atrstm_component cpnt,
+ const void* voxels[],
+ const size_t nvoxs)
+{
+ float ka[2] = {FLT_MAX,-FLT_MAX};
+ float ks[2] = {FLT_MAX,-FLT_MAX};
+ float kext[2] = {FLT_MAX,-FLT_MAX};
+ size_t ivox;
+ ASSERT(dst && voxs && nvoxs);
+
+ FOR_EACH(ivox, 0, nvoxs) {
+ const float* vox = voxels[ivox];
+ ka[0] = MMIN(ka[0], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN));
+ ka[1] = MMAX(ka[1], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX));
+ ks[0] = MMIN(ks[0], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN));
+ ks[1] = MMAX(ks[1], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX));
+ kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
+ kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
+ }
+ vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, ka[0]);
+ vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, ka[1]);
+ vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, ks[0]);
+ vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, ks[1]);
+ vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, kext[0]);
+ vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, kext[1]);
+}
+
+static void
+voxels_merge
+ (void* dst,
+ const void* voxels[],
+ const size_t nvoxs,
+ void* ctx)
+{
+ (void)ctx;
+ voxels_merge_component(dst, ATRSTM_CPNT_GAS, voxels, nvoxs);
+ voxels_merge_component(dst, ATRSTM_CPNT_SOOT, voxels, nvoxs);
+}
+
+static INLINE void
+voxels_component_compute_kext_range
+ (const enum atrstm_component cpnt,
+ const struct svx_voxel voxels[],
+ const size_t nvoxs,
+ float kext[2])
+{
+ size_t ivox;
+ ASSERT(voxels && nvoxs && kext);
+
+ kext[0] = FLT_MAX;
+ kext[1] =-FLT_MAX;
+
+ /* Compute the Kext range of the submitted voxels */
+ FOR_EACH(ivox, 0, nvoxs) {
+ const float* vox = voxels[ivox].data;
+ kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
+ kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
+ }
+}
+
+static int
+voxels_challenge_merge
+ (const struct svx_voxel voxels[],
+ const size_t nvoxs,
+ void* context)
+{
+ const struct build_octree_context* ctx = context;
+ double lower[3] = { DBL_MAX, DBL_MAX, DBL_MAX};
+ double upper[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX};
+ double sz_max;
+ size_t ivox;
+ float kext_gas[2] = {0, 0};
+ float kext_soot[2] = {0, 0};
+ int merge_gas = 0;
+ int merge_soot = 0;
+ ASSERT(voxels && nvoxs && context);
+
+ /* Compute the AABB of the group of voxels and its maximum size */
+ FOR_EACH(ivox, 0, nvoxs) {
+ lower[0] = MMIN(voxels[ivox].lower[0], lower[0]);
+ lower[1] = MMIN(voxels[ivox].lower[1], lower[1]);
+ lower[2] = MMIN(voxels[ivox].lower[2], lower[2]);
+ upper[0] = MMAX(voxels[ivox].upper[0], upper[0]);
+ upper[1] = MMAX(voxels[ivox].upper[1], upper[1]);
+ upper[2] = MMAX(voxels[ivox].upper[2], upper[2]);
+ }
+ sz_max = upper[0] - lower[0];
+ sz_max = MMAX(upper[1] - lower[1], sz_max);
+ sz_max = MMAX(upper[2] - lower[2], sz_max);
+
+ /* Compute the Kext range for each component */
+ voxels_component_compute_kext_range(ATRSTM_CPNT_GAS, voxels, nvoxs, kext_gas);
+ voxels_component_compute_kext_range(ATRSTM_CPNT_SOOT, voxels, nvoxs, kext_soot);
+ ASSERT(kext_gas[1] >= kext_gas[0]);
+ ASSERT(kext_soot[1] >= kext_soot[0]);
+
+ /* Check if the voxels are candidates to the merge process by evaluating its
+ * optical thickness regarding the maximum size of their AABB and a user
+ * defined threshold */
+ merge_gas = (kext_gas[1] - kext_gas[0])*sz_max <= ctx->tau_threshold;
+ merge_soot = (kext_soot[1] - kext_soot[0])*sz_max <= ctx->tau_threshold;
+ return merge_gas && merge_soot;
+}
+
+static res_T
+build_octree
+ (struct atrstm* atrstm,
+ const struct build_octrees_args* args,
+ struct pool* pools)
+{
+ struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL;
+ struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
+ double lower[3];
+ double upper[3];
+ size_t definition[3];
+ res_T res = RES_OK;
+ ASSERT(atrstm && args && pools);
+
+ /* Setup the build octree context */
+ ctx.atrstm = atrstm;
+ ctx.pools = pools;
+ ctx.tau_threshold = args->optical_thickness_threshold;
+
+ /* Setup the voxel descriptor */
+ vox_desc.get = voxel_get;
+ vox_desc.merge = voxels_merge;
+ vox_desc.challenge_merge = voxels_challenge_merge;
+ vox_desc.context = &ctx;
+ vox_desc.size = NFLOATS_PER_VOXEL * sizeof(float);
+
+ SUVM(volume_get_aabb(atrstm->volume, lower, upper));
+
+ /* Create the octree */
+ definition[0] = (size_t)args->grid_max_definition[0];
+ definition[1] = (size_t)args->grid_max_definition[1];
+ definition[2] = (size_t)args->grid_max_definition[2];
+ res = svx_octree_create
+ (atrstm->svx, lower, upper, definition, &vox_desc, &atrstm->octree);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+build_octrees(struct atrstm* atrstm, const struct build_octrees_args* args)
+{
+ struct pool* pools = NULL;
+ size_t i;
+ ATOMIC res = RES_OK;
+ ASSERT(atrstm && check_build_octrees_args(args));
+
+ /* Allocate the per thread pool of partition */
+ pools = MEM_CALLOC(atrstm->allocator, atrstm->nthreads, sizeof(*pools));
+ if(!pools) { res = RES_MEM_ERR; goto error; }
+ FOR_EACH(i, 0, atrstm->nthreads) {
+ res = pool_init(atrstm->allocator, pools+i);
+ if(res != RES_OK) goto error;
+ }
+
+ omp_set_nested(1); /* Enable nested threads for voxelize_volumetric_mesh */
+ #pragma omp parallel sections num_threads(2)
+ {
+ #pragma omp section
+ {
+ res_T res_local = voxelize_volumetric_mesh(atrstm, args, pools);
+ if(res_local != RES_OK) {
+ size_t ipool;
+ log_err(atrstm, "Error voxelizing the volumetric mesh -- %s\n",
+ res_to_cstr(res_local));
+ FOR_EACH(ipool, 0, atrstm->nthreads) {
+ pool_invalidate(pools+ipool);
+ }
+ ATOMIC_SET(&res, res_local);
+ }
+ }
+
+ #pragma omp section
+ {
+ res_T res_local = build_octree(atrstm, args, pools);
+ if(res_local != RES_OK) {
+ size_t ipool;
+ log_err(atrstm, "Error buildin the octree -- %s\n",
+ res_to_cstr(res_local));
+ FOR_EACH(ipool, 0, atrstm->nthreads) {
+ pool_invalidate(pools+ipool);
+ }
+ ATOMIC_SET(&res, res_local);
+ }
+ }
+ }
+
+exit:
+ if(pools) {
+ FOR_EACH(i, 0, atrstm->nthreads) pool_release(pools+i);
+ MEM_RM(atrstm->allocator, pools);
+ }
+ return (res_T)res;
+error:
+ goto exit;
+}
diff --git a/src/atrstm_build_octrees.h b/src/atrstm_build_octrees.h
@@ -0,0 +1,51 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef ATRSTM_BUILD_OCTREES_H
+#define ATRSTM_BUILD_OCTREES_H
+
+#include <rsys/rsys.h>
+
+/* Forward declaration of external data type */
+struct sth;
+
+struct cache {
+ FILE* stream; /* Stream where data are stored */
+ const char* name; /* Name of the stream where data are stored */
+ int force_update; /* Force overwrite of the cached data */
+};
+#define CACHE_NULL__ {NULL, NULL, 0}
+static const struct cache CACHE_NULL = CACHE_NULL__;
+
+struct build_octrees_args {
+ unsigned grid_max_definition[3]; /* Max definition of the Octrees */
+ double optical_thickness_threshold; /* Threshold used to build the octrees */
+ struct cache cache; /* Cache of the acceleration data structures */
+};
+#define BUILD_OCTREES_ARGS_NULL__ {{0,0,0}, 0, CACHE_NULL__}
+static const struct build_octrees_args BUILD_OCTREES_ARGS_NULL =
+ BUILD_OCTREES_ARGS_NULL__;
+
+extern LOCAL_SYM res_T
+build_octrees
+ (struct atrstm* atstmm,
+ const struct build_octrees_args* args);
+
+extern LOCAL_SYM void
+octrees_clean
+ (struct atrstm* atrstm);
+
+#endif /* ATRSTM_BUILD_OCTREES_H */
+
diff --git a/src/atrstm_c.h b/src/atrstm_c.h
@@ -32,6 +32,9 @@ struct atrstm {
struct suvm_device* suvm;
struct suvm_volume* volume;
+ struct svx_device* svx;
+ struct svx_tree* octree;
+
unsigned nthreads; /* #nthreads */
struct str name; /* Name of the gas mixture */
diff --git a/src/atrstm_partition.c b/src/atrstm_partition.c
@@ -19,8 +19,6 @@
#include <rsys/mem_allocator.h>
#include <rsys/mutex.h>
-#include <string.h>
-
/******************************************************************************
* Local functions
******************************************************************************/
@@ -92,7 +90,7 @@ pool_release(struct pool* pool)
}
struct part*
-pool_new_partition(struct pool* pool)
+pool_new_partition(struct pool* pool, const size_t ipart)
{
struct list_node* node = NULL;
struct part* part = NULL;
@@ -102,11 +100,17 @@ pool_new_partition(struct pool* pool)
if(is_list_empty(&pool->parts_free)) {
cond_wait(pool->cond_new, pool->mutex);
}
- node = list_head(&pool->parts_free);
- list_del(node);
- mutex_unlock(pool->mutex);
+ if(pool->error) {
+ part = NULL;
+ mutex_unlock(pool->mutex);
+ } else {
+ node = list_head(&pool->parts_free);
+ list_del(node);
+ mutex_unlock(pool->mutex);
- part = CONTAINER_OF(node, struct part, node);
+ part = CONTAINER_OF(node, struct part, node);
+ part->id = ipart;
+ }
return part;
}
@@ -147,7 +151,7 @@ pool_fetch_partition(struct pool* pool, size_t part_id)
ASSERT(pool);
mutex_lock(pool->mutex);
- for(;;) {
+ while(!pool->error) {
/* Search for a partition that match the submitted id */
LIST_FOR_EACH(node, &pool->parts_full) {
struct part* part = CONTAINER_OF(node, struct part, node);
@@ -173,3 +177,15 @@ pool_fetch_partition(struct pool* pool, size_t part_id)
return found_part;
}
+void
+pool_invalidate(struct pool* pool)
+{
+ ASSERT(pool);
+ mutex_lock(pool->mutex);
+ pool->error = 1;
+ mutex_unlock(pool->mutex);
+
+ cond_broadcast(pool->cond_new);
+ cond_broadcast(pool->cond_fetch);
+}
+
diff --git a/src/atrstm_partition.h b/src/atrstm_partition.h
@@ -18,9 +18,15 @@
#include "atrstm_svx.h"
#include <rsys/list.h>
+#include <string.h>
/* Definition of a partition along the 4 axis */
-#define PARTITION_DEFINITION 32
+#define LOG2_PARTITION_DEFINITION 5
+#define PARTITION_DEFINITION BIT(LOG2_PARTITION_DEFINITION) /*32*/
+#define PARTITION_NVOXELS \
+ ( PARTITION_DEFINITION \
+ * PARTITION_DEFINITION \
+ * PARTITION_DEFINITION)
/* Max #partitions stored into the partition pool */
#define POOL_MAX_NPARTITIONS 16
@@ -32,11 +38,7 @@ struct mutex;
struct part {
/* Morton ordered list of voxels */
- float voxels
- [ PARTITION_DEFINITION
- * PARTITION_DEFINITION
- * PARTITION_DEFINITION
- * NFLOATS_PER_VOXEL];
+ float voxels[PARTITION_NVOXELS * NFLOATS_PER_VOXEL];
struct list_node node;
size_t id; /* Unique identifier of the partition */
};
@@ -45,7 +47,7 @@ struct pool {
/* Linked list of pre-allocated partitions */
struct list_node parts_free;
- /* List of available partition sorted in ascending order wrt partition id */
+ /* List of available partition sorted in ascending order wrt partition id */
struct list_node parts_full;
struct mutex* mutex;
@@ -53,8 +55,27 @@ struct pool {
struct cond* cond_fetch;
struct mem_allocator* allocator;
+
+ ATOMIC error;
};
+static FINLINE void
+part_clear_voxels(struct part* part)
+{
+ ASSERT(part);
+ memset(part, 0, sizeof(part->voxels));
+}
+
+static FINLINE float*
+part_get_voxel
+ (struct part* part,
+ const size_t ivoxel) /* Morton code of the voxel */
+{
+ ASSERT(part);
+ ASSERT(ivoxel < PARTITION_NVOXELS);
+ return part->voxels + ivoxel * NFLOATS_PER_VOXEL;
+}
+
extern LOCAL_SYM res_T
pool_init
(struct mem_allocator* allocator, /* May be NULL <=> default allocator */
@@ -64,10 +85,12 @@ extern LOCAL_SYM void
pool_release
(struct pool* pool);
-/* Return a free partition. Wait until a free partition is available */
+/* Return a free partition. Wait until a free partition is available. Return
+ * NULL on error. */
extern LOCAL_SYM struct part*
pool_new_partition
- (struct pool* pool);
+ (struct pool* pool,
+ const size_t ipart); /* Identifier of the partition */
/* Register the partition as a free partition against the pool */
extern LOCAL_SYM void
@@ -82,10 +105,14 @@ pool_commit_partition
struct part* part);
/* Return the partition whose id is `part_id`. Wait until the expected
- * partition is available */
+ * partition is available. Return NULL on error. */
extern LOCAL_SYM struct part*
pool_fetch_partition
(struct pool* pool,
const size_t part_id);
+extern LOCAL_SYM void
+pool_invalidate
+ (struct pool* pool);
+
#endif /* ATRSTM_PARTITION_H */
