star-line

commit fe3d4473bae17697db5c6d87aa84b21dc5c2d8ee
parent c0b032e352a1d7b736f0daef3230e6c9318ae4e5
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  6 Mar 2026 15:29:35 +0100

Update how the upper mesh of a line is calculated

Until now, the line was meshed by evaluating the value of the line at
each selected wave number. Then, only if the mesh was to be an upper
bound of its associated line, the value of the mesh vertex was changed
to the value of the previous vertex (only the upper half of the line,
strictly decreasing, being meshed).

While this worked reasonably well with a fine mesh, it is no longer the
case with a coarse mesh, as has been the situation since validation
a8b999b. Hence this validation, which ensures an upper bound mesh for
the line by calculating its value directly from a wave number slightly
ahead of the actual wave number. The offset to be applied to the
original wave number is dictated by the parameter on the number of
vertices to be used to mesh each line, so that the user has some control
over this numerical parameter.

Diffstat:
M
src/sln_line.c
 | 
107
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------

1 file changed, 81 insertions(+), 26 deletions(-)
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -246,7 +246,7 @@ error:
 
 /* Calculate line values for a set of wave numbers */
 static res_T
-eval_mesh
+eval_mesh_fit
   (const struct sln_tree* tree,
    const struct line* line,
    const struct darray_double* wavenumbers,
@@ -254,7 +254,8 @@ eval_mesh
 {
   const double* nu = NULL;
   double* ha = NULL;
-  size_t ivertex, nvertices;
+  size_t ivertex = 0;
+  size_t nvertices = 0;
   res_T res = RES_OK;
   ASSERT(tree && line && wavenumbers && values);
 
@@ -277,34 +278,88 @@ error:
   goto exit;
 }
 
-static void
-snap_mesh_to_upper_bound
-  (const struct darray_double* wavenumbers,
+/* Calculate an upper bound for the line values for a set of wave numbers */
+static res_T
+eval_mesh_upper
+  (const struct sln_tree* tree,
+   const struct line* line,
+   const struct darray_double* wavenumbers,
+  /* #vertices used to mesh the line on its first interval, from the center of
+   * the line to the width of the line at mid-height. This parameter is used to
+   * calculate an upper bound mesh for the line in accordance with its width.
+   * This user-defined discretization parameter controls the desired mesh
+   * approximation. Thus, the mesh will be upper bounding “but not too much”
+   * with respect to this numerical parameter */
+   const size_t nvertices_around_line_center,
    struct darray_double* values)
 {
+  const double* nu = NULL;
   double* ha = NULL;
-  size_t ivertex, nvertices;
-  ASSERT(wavenumbers && values);
-  ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values));
-  (void)wavenumbers;
+  double nu_offset = 0;
+  size_t ivertex = 0;
+  size_t nvertices = 0;
+  res_T res = RES_OK;
+  ASSERT(tree && line && wavenumbers && values && nvertices_around_line_center);
+
+  /* The line is meshed on its upper half, which is a strictly decreasing
+   * function. To ensure that the mesh is an upper bound of the line, it is
+   * therefore sufficient to evaluate its value at a wave number slightly lower
+   * than the current wave number.
+   *
+   * This refers to the behavior of the following nu_offset. It corresponds to
+   * an offset to be applied to the current nu so as to evaluate the line
+   * slightly before the actual nu, and thus have a line value greater than its
+   * actual value. */
+  nu_offset = line->gamma_l / (double)nvertices_around_line_center;
 
-  ha = darray_double_data_get(values);
   nvertices = darray_double_size_get(wavenumbers);
+  ASSERT(nvertices);
 
-  /* Ensure that the stored vertex value is an exclusive upper bound of the
-   * original value. We do this by storing a value in single precision that is
-   * strictly greater than its encoding in double precision */
-  if(ha[0] != (float)ha[0]) {
-    ha[0] = nextafterf((float)ha[0], FLT_MAX);
-  }
+  res = darray_double_resize(values, nvertices);
+  if(res != RES_OK) goto error;
 
-  /* We have meshed the upper half of the line which is a strictly decreasing
-   * function. To ensure that the mesh is an upper limit of this function,
-   * simply align the value of each vertex with the value of the preceding
-   * vertex */
-  FOR_EACH_REVERSE(ivertex, nvertices-1, 0) {
-    ha[ivertex] = ha[ivertex-1];
+  nu = darray_double_cdata_get(wavenumbers);
+  ha = darray_double_data_get(values);
+
+  FOR_EACH(ivertex, 0, nvertices) {
+    double nu_adjusted = nu[ivertex];
+
+   /* The first node is actually the center of the line, so it is the absolute
+    * upper limit of the line. No offset should therefore be applied to its wave
+    * number to ensure that the resulting mesh is an upper bound of the line */
+    if(ivertex == 0) {
+      ASSERT(line->wavenumber == nu_adjusted);
+    } else {
+      nu_adjusted -= nu_offset;
+
+      if(nu_adjusted == nu[ivertex]) {
+        /* Offset was too small regarding the numeric precision */
+        nu_adjusted = nextafter(-DBL_MAX, nu_adjusted);
+        nu_adjusted = MMIN(nu_adjusted, line->wavenumber);
+      }
+    }
+
+    ha[ivertex] = line_value(tree, line, nu_adjusted);
+
+#ifndef NDEBUG
+    if(ivertex != 0) {
+      double ha_ref = line_value(tree, line, nu[ivertex]);
+      ASSERT(ha_ref < ha[ivertex]);
+    }
+#endif
+
+    /* Ensure that the value of the vertex, stored in single precision, is
+     * indeed an exclusive upper bound of the original value. */
+    if(ha[ivertex] != (float)ha[ivertex]) {
+      ha[ivertex] = nextafterf((float)ha[ivertex], FLT_MAX);
+    }
   }
+
+exit:
+  return res;
+error:
+  ERROR(tree->sln, "Error evaluating the line mesh -- %s.\n", res_to_cstr(res));
+  goto exit;
 }
 
 static INLINE int
@@ -534,13 +589,13 @@ line_mesh
 
   /* Evaluate the mesh vertices, i.e. define the line value for the list of
    * wavenumbers */
-  eval_mesh(tree, &line, &wavenumbers, &values);
-
   switch(tree->args.mesh_type) {
     case SLN_MESH_UPPER:
-      snap_mesh_to_upper_bound(&wavenumbers, &values);
+      eval_mesh_upper(tree, &line, &wavenumbers, nvertices_adjusted, &values);
+      break;
+    case SLN_MESH_FIT:
+      eval_mesh_fit(tree, &line, &wavenumbers, &values);
       break;
-    case SLN_MESH_FIT: /* Do nothing */ break;
     default: FATAL("Unreachable code.\n"); break;
   }