commit 0daf7eac33f62de458131fe9c47320bfdb4ad1e0
parent e131f6e67e45fee616d525b19a7b3d66d2a896a5
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 7 Aug 2018 14:13:27 +0200
Add and test the htcp_desc_get_voxel_aabb function
Diffstat:
4 files changed, 116 insertions(+), 7 deletions(-)
diff --git a/src/htcp.c b/src/htcp.c
@@ -37,6 +37,7 @@ struct htcp {
double lower[3];
double vxsz[2]; /* Size of the voxels in X and Y */
struct darray_double vxsz_z; /* Size of the voxels along the Z dimension */
+ struct darray_double coord_z; /* Lower coordinate of the voxel along Z */
double* RCT; /* Mapped memory */
double* RVT; /* Mapped memory */
@@ -96,6 +97,7 @@ reset_htcp(struct htcp* htcp)
htcp->vxsz[0] = -1;
htcp->vxsz[1] = -1;
darray_double_clear(&htcp->vxsz_z);
+ darray_double_clear(&htcp->coord_z);
if(htcp->RCT) munmap(htcp->RCT, htcp->RCT_length);
if(htcp->RVT) munmap(htcp->RVT, htcp->RVT_length);
if(htcp->PABST) munmap(htcp->PABST, htcp->PABST_length);
@@ -165,12 +167,26 @@ load_stream(struct htcp* htcp, FILE* stream, const char* stream_name)
log_err(htcp,
"%s: could not allocate memory to store the size of the voxels in Z.\n",
stream_name);
- res = RES_MEM_ERR;
goto error;
}
READ(darray_double_data_get(&htcp->vxsz_z), nz, "Z voxel size(s)");
#undef READ
+ if(htcp->irregular_z) {
+ const double* size = NULL;
+ double* coord = NULL;
+ size_t i;
+ res = darray_double_resize(&htcp->coord_z, nz);
+ if(res != RES_OK) {
+ log_err(htcp, "%s: could not allocate memory to store the lower "
+ "coordinate of the voxels in Z.\n", FUNC_NAME);
+ goto error;
+ }
+ size = darray_double_cdata_get(&htcp->vxsz_z);
+ coord = darray_double_data_get(&htcp->coord_z);
+ FOR_EACH(i, 0, nz) coord[i] = i ? coord[i-1] + size[i-1] : htcp->lower[2];
+ }
+
map_len =
(size_t)htcp->definition[X]
* (size_t)htcp->definition[Y]
@@ -224,6 +240,7 @@ release_htcp(ref_T* ref)
htcp = CONTAINER_OF(ref, struct htcp, ref);
reset_htcp(htcp);
darray_double_release(&htcp->vxsz_z);
+ darray_double_release(&htcp->coord_z);
MEM_RM(htcp->allocator, htcp);
}
@@ -265,6 +282,7 @@ htcp_create
htcp->verbose = verbose;
htcp->pagesize_os = (size_t)sysconf(_SC_PAGESIZE);
darray_double_init(htcp->allocator, &htcp->vxsz_z);
+ darray_double_init(htcp->allocator, &htcp->coord_z);
exit:
if(out_htcp) *out_htcp = htcp;
@@ -345,6 +363,7 @@ htcp_get_desc(const struct htcp* htcp, struct htcp_desc* desc)
desc->vxsz_x = htcp->vxsz[0];
desc->vxsz_y = htcp->vxsz[1];
desc->vxsz_z = darray_double_cdata_get(&htcp->vxsz_z);
+ desc->coord_z = htcp->irregular_z ? darray_double_cdata_get(&htcp->coord_z) : NULL;
desc->RCT = htcp->RCT;
desc->RVT = htcp->RVT;
desc->PABST = htcp->PABST;
diff --git a/src/htcp.h b/src/htcp.h
@@ -56,13 +56,15 @@ struct htcp_desc {
double vxsz_x; /* Voxel size in X */
double vxsz_y; /* Voxel size in Y */
const double* vxsz_z; /* Voxel size along Z */
+ const double* coord_z; /* Voxel lower bound along Z. NULL if !irregular_z */
const double* RCT;
const double* RVT;
const double* PABST; /* Pressure */
const double* T; /* Temperature */
};
-#define HTCP_DESC_NULL__ {0,-1,{0,0,0},0,{-1,-1,-1},-1,-1,NULL,NULL,NULL,NULL,NULL}
+#define HTCP_DESC_NULL__ \
+ {0,-1,{0,0,0},0,{-1,-1,-1},-1,-1,NULL,NULL,NULL,NULL,NULL,NULL}
static const struct htcp_desc HTCP_DESC_NULL = HTCP_DESC_NULL__;
BEGIN_DECLS
@@ -105,7 +107,7 @@ static FINLINE double
htcp_dblgrid4D_at__
(const double* grid,
const struct htcp_desc* desc,
- size_t x, size_t y, size_t z, size_t t)
+ const size_t x, const size_t y, const size_t z, const size_t t)
{
size_t row, slice, array;
ASSERT(desc && grid);
@@ -122,7 +124,7 @@ htcp_dblgrid4D_at__
static FINLINE double
htcp_desc_RCT_at
(const struct htcp_desc* desc,
- size_t x, size_t y, size_t z, size_t t)
+ const size_t x, const size_t y, const size_t z, const size_t t)
{
return htcp_dblgrid4D_at__(desc->RCT, desc, x, y, z, t);
}
@@ -131,7 +133,7 @@ htcp_desc_RCT_at
static FINLINE double
htcp_desc_RVT_at
(const struct htcp_desc* desc,
- size_t x, size_t y, size_t z, size_t t)
+ const size_t x, const size_t y, const size_t z, const size_t t)
{
return htcp_dblgrid4D_at__(desc->RVT, desc, x, y, z, t);
}
@@ -139,7 +141,7 @@ htcp_desc_RVT_at
static FINLINE double
htcp_desc_PABST_at
(const struct htcp_desc* desc,
- size_t x, size_t y, size_t z, size_t t)
+ const size_t x, const size_t y, const size_t z, const size_t t)
{
return htcp_dblgrid4D_at__(desc->PABST, desc, x, y, z, t);
}
@@ -147,11 +149,35 @@ htcp_desc_PABST_at
static FINLINE double
htcp_desc_T_at
(const struct htcp_desc* desc,
- size_t x, size_t y, size_t z, size_t t)
+ const size_t x, const size_t y, const size_t z, const size_t t)
{
return htcp_dblgrid4D_at__(desc->T, desc, x, y, z, t);
}
+static FINLINE void
+htcp_desc_get_voxel_aabb
+ (const struct htcp_desc* desc,
+ const size_t x, const size_t y, const size_t z,
+ double lower[3],
+ double upper[3])
+{
+ ASSERT(desc && lower && upper);
+ ASSERT(x < desc->spatial_definition[0]);
+ ASSERT(y < desc->spatial_definition[1]);
+ ASSERT(z < desc->spatial_definition[2]);
+ lower[0] = (double)x * desc->vxsz_x;
+ lower[1] = (double)y * desc->vxsz_y;
+ upper[0] = lower[0] + desc->vxsz_x;
+ upper[1] = lower[1] + desc->vxsz_y;
+ if(!desc->irregular_z) {
+ lower[2] = (double)z * desc->vxsz_z[0];
+ upper[2] = lower[2] + desc->vxsz_z[0];
+ } else {
+ lower[2] = desc->coord_z[z];
+ upper[2] = lower[2] + desc->vxsz_z[z];
+ }
+}
+
END_DECLS
#endif /* HTCP_H */
diff --git a/src/test_htcp_load.c b/src/test_htcp_load.c
@@ -16,6 +16,8 @@
#include "htcp.h"
#include "test_htcp_utils.h"
+#include <rsys/double3.h>
+
#include <unistd.h>
int
@@ -113,6 +115,26 @@ main(int argc, char** argv)
CHK(desc.vxsz_x == 1);
CHK(desc.vxsz_y == 2);
CHK(desc.vxsz_z[0] == 3);
+ CHK(desc.coord_z == NULL);
+
+ FOR_EACH(x, 0, desc.spatial_definition[0]) {
+ double low[3];
+ low[0] = (double)x * desc.vxsz_x;
+ FOR_EACH(y, 0, desc.spatial_definition[1]) {
+ low[1] = (double)y * desc.vxsz_y;
+ FOR_EACH(z, 0, desc.spatial_definition[2]) {
+ double vox_low[3], vox_upp[3];
+ double upp[3];
+ low[2] = (double)z * desc.vxsz_z[0];
+ upp[0] = low[0] + desc.vxsz_x;
+ upp[1] = low[1] + desc.vxsz_y;
+ upp[2] = low[2] + desc.vxsz_z[0];
+ htcp_desc_get_voxel_aabb(&desc, x, y, z, vox_low, vox_upp);
+ CHK(d3_eq_eps(vox_low, low, 1.e-6));
+ CHK(d3_eq_eps(vox_upp, upp, 1.e-6));
+ }
+ }
+ }
n = desc.spatial_definition[0]
* desc.spatial_definition[1]
diff --git a/src/test_htcp_load_from_file.c b/src/test_htcp_load_from_file.c
@@ -18,6 +18,7 @@
#include "htcp.h"
#include "test_htcp_utils.h"
+#include <rsys/double3.h>
#include <rsys/math.h>
#include <string.h>
@@ -96,6 +97,45 @@ check_variable
CHK(fclose(fp) == 0);
}
+static void
+check_misc(const struct htcp* htcp)
+{
+ struct htcp_desc desc;
+ double low[3], upp[3];
+ size_t x, y, z;
+
+ printf("Check miscellaneous\n");
+
+ CHK(htcp_get_desc(htcp, &desc) == RES_OK);
+ CHK(desc.irregular_z || desc.coord_z == NULL);
+
+ FOR_EACH(x, 0, desc.spatial_definition[0]) {
+ low[0] = (double)x * desc.vxsz_x;
+ upp[0] = low[0] + desc.vxsz_x;
+ FOR_EACH(y, 0, desc.spatial_definition[1]) {
+ low[1] = (double)y * desc.vxsz_y;
+ upp[1] = low[1] + desc.vxsz_y;
+ FOR_EACH(z, 0, desc.spatial_definition[2]) {
+ double vox_low[3];
+ double vox_upp[3];
+ if(!desc.irregular_z) {
+ low[2] = (double)z * desc.vxsz_z[0];
+ upp[2] = low[2] + desc.vxsz_z[0];
+ } else {
+ size_t i;
+ low[2] = desc.lower[2];
+ FOR_EACH(i, 0, z) low[2] += desc.vxsz_z[i];
+ CHK(eq_eps(low[2], desc.coord_z[z], 1.e-6));
+ upp[2] = low[2] + desc.vxsz_z[z];
+ }
+ htcp_desc_get_voxel_aabb(&desc, x, y, z, vox_low, vox_upp);
+ CHK(d3_eq_eps(vox_low, low, 1.e-6));
+ CHK(d3_eq_eps(vox_upp, upp, 1.e-6));
+ }
+ }
+ }
+}
+
int
main(int argc, char** argv)
{
@@ -132,7 +172,9 @@ main(int argc, char** argv)
* desc.spatial_definition[2]
* desc.time_definition;
+ printf("Irregular Z: %i\n", desc.irregular_z);
check_descriptor(&desc, path, base);
+ check_misc(htcp);
check_variable(desc.RVT, n, "RVT", path, base, NULL, NULL);
check_variable(desc.RCT, n, "RCT", path, base, NULL, NULL);
check_variable(desc.PABST, n, "PABST", path, base, NULL, NULL);
