rnatm

commit 2847f14eaa3231788b44b3609ed77f3a624ba414
parent 1c670d11cbca39f4546a55fa3a0fa4eee655a01e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  8 Nov 2022 11:12:26 +0100

Fix rnatm_cell_get_radcoef precision inaccuracy

When a cell has the same radiative coefficient for all its vertices,
interpolation introduced an imprecision that made the result not
strictly equal to the radiative coefficient shared by the vertices. As a
result, this result may not be included in the k_min, k_max of the
corresponding node of the acceleration structure. In this commit, we
avoid interpolating the radiative coefficients of the vertices of the
cell if they are equal.

Diffstat:
M
src/rnatm.h
 | 
7
++++---
M
src/rnatm_properties.c
 | 
29
++++++++++++++++++-----------

2 files changed, 22 insertions(+), 14 deletions(-)
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -143,8 +143,9 @@ struct rnatm_cell_get_radcoef_args {
 
   enum rnatm_radcoef radcoef;
 
-  /* For debug: check that the retrieved radcoef is in [k_min, k_max] */
-  double k_min;
+  /* For debug: check that the retrieved radcoef is < k_max. Do not provide
+   * k_min because the radiative coefficient of a component may be less than
+   * the radiative coefficient of the mixture */
   double k_max;
 };
 #define RNATM_CELL_GET_RADCOEF_ARGS_NULL__ {                                   \
@@ -152,7 +153,7 @@ struct rnatm_cell_get_radcoef_args {
   {0, 0, 0, 0}, /* Barycentric coordinates */                                  \
   0, 0, /* Spectral data (band, quadrature pointh) */                          \
   RNATM_RADCOEFS_COUNT__,                                                      \
-  -DBL_MAX, DBL_MAX /* For debug: Radcoef range */                             \
+  DBL_MAX /* For debug: Radcoef range */                                       \
 }
 static const struct rnatm_cell_get_radcoef_args
 RNATM_CELL_GET_RADCOEF_ARGS_NULL = RNATM_CELL_GET_RADCOEF_ARGS_NULL__;
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -151,7 +151,7 @@ check_rnatm_cell_get_radcoef_args
     return RES_BAD_ARG;
 
   /* Invalid K range */
-  if(args->k_min > args->k_max)
+  if(args->k_max < 0)
     return RES_BAD_ARG;
 
   return RES_OK;
@@ -1111,7 +1111,7 @@ rnatm_cell_get_radcoef
    const struct rnatm_cell_get_radcoef_args* args,
    double* out_k)
 {
-  float per_vertex_k[4];
+  float vtx_k[4];
   double k = NaN;
   res_T res = RES_OK;
 
@@ -1121,15 +1121,22 @@ rnatm_cell_get_radcoef
 
   /* Get the radiative coefficients of the tetrahedron */
   tetra_get_radcoef(atm, &args->cell.prim, args->cell.component, args->iband,
-    args->iquad, args->radcoef, per_vertex_k);
-
-  /* Calculate the radiative coefficient by linearly interpolating the
-   * coefficients defined by tetrahedron vertex */
-  k = per_vertex_k[0] * args->barycentric_coords[0]
-    + per_vertex_k[1] * args->barycentric_coords[1]
-    + per_vertex_k[2] * args->barycentric_coords[2]
-    + per_vertex_k[3] * args->barycentric_coords[3];
-  ASSERT(args->k_min <= k && k <= args->k_max);
+    args->iquad, args->radcoef, vtx_k);
+
+  if(vtx_k[0] == vtx_k[1]
+  && vtx_k[0] == vtx_k[2]
+  && vtx_k[0] == vtx_k[3]) {
+    /* Avoid precision issue when iterpolating the same value */
+    k = vtx_k[0];
+  } else {
+    /* Calculate the radiative coefficient by linearly interpolating the
+     * coefficients defined by tetrahedron vertex */
+    k = vtx_k[0] * args->barycentric_coords[0]
+      + vtx_k[1] * args->barycentric_coords[1]
+      + vtx_k[2] * args->barycentric_coords[2]
+      + vtx_k[3] * args->barycentric_coords[3];
+  }
+  ASSERT(k <= args->k_max);
 
 exit:
   *out_k = k;