extract line

src/ref_iso.c

@@ -876,6 +876,86 @@ REF_FCN REF_STATUS ref_iso_boom_zone(FILE *file, REF_GRID ref_grid,
   return REF_SUCCESS;
 }
 
+REF_FCN REF_STATUS ref_iso_boomray(const char *filename, REF_GRID ref_grid,
+                                   REF_DBL *field, REF_INT ldim,
+                                   const char **scalar_names, REF_DBL *segment0,
+                                   REF_DBL *segment1) {
+  FILE *file;
+  REF_DBL ds[3], dt[3];
+  REF_GRID ray_grid;
+  REF_DBL *ray_field;
+  REF_DBL *local_xyzf, *xyzf, *t;
+  REF_NODE ref_node;
+  REF_MPI ref_mpi = ref_grid_mpi(ref_grid);
+  REF_INT node, local_n, n, *source, i, *order;
+  file = fopen(filename, "w");
+  if (NULL == file) printf("unable to open %s\n", filename);
+  RNS(file, "unable to open file");
+  fprintf(file, "title=\"tecplot refine gather\"\n");
+  fprintf(file, "variables = \"x\" \"y\" \"z\"");
+  if (NULL != scalar_names) {
+    for (i = 0; i < ldim; i++) fprintf(file, " \"%s\"", scalar_names[i]);
+  } else {
+    for (i = 0; i < ldim; i++) fprintf(file, " \"V%d\"", i + 1);
+  }
+  fprintf(file, "\n");
+  for (i = 0; i < 3; i++) ds[i] = segment1[i] - segment0[i];
+  RSS(ref_math_normalize(ds), "segment unit vector");
+
+  RSS(ref_iso_cast(&ray_grid, &ray_field, ref_grid, field, ldim, segment0,
+                   segment1),
+      "cast");
+  ref_node = ref_grid_node(ray_grid);
+  local_n = 0;
+  each_ref_node_valid_node(ref_node, node) {
+    if (ref_node_owned(ref_node, node)) {
+      local_n++;
+    }
+  }
+  ref_malloc(local_xyzf, local_n * (3 + ldim), REF_DBL);
+  local_n = 0;
+  each_ref_node_valid_node(ref_node, node) {
+    if (ref_node_owned(ref_node, node)) {
+      for (i = 0; i < 3; i++)
+        local_xyzf[i + (3 + ldim) * local_n] = ref_node_xyz(ref_node, i, node);
+      for (i = 0; i < ldim; i++)
+        local_xyzf[3 + i + (3 + ldim) * local_n] = ray_field[i + node * ldim];
+      local_n++;
+    }
+  }
+  ref_free(ray_field);
+  ref_grid_free(ray_grid);
+
+  RSS(ref_mpi_allconcat(ref_mpi, 3 + ldim, local_n, local_xyzf, &n, &source,
+                        (void **)(&xyzf), REF_DBL_TYPE),
+      "concat");
+  ref_free(local_xyzf);
+  if (ref_mpi_once(ref_mpi)) {
+    ref_free(source);
+    ref_malloc(t, n, REF_DBL);
+
+    for (node = 0; node < n; node++) {
+      for (i = 0; i < 3; i++) dt[i] = xyzf[i + (3 + ldim) * node] - segment0[i];
+      t[node] = ref_math_dot(dt, ds);
+    }
+    ref_malloc(order, n, REF_INT);
+    RSS(ref_sort_heap_dbl(n, t, order), "sort t");
+    fprintf(file, " zone t=\"((%.2f,%.2f,%.2f),(%.2f,%.2f,%.2f))\"\n",
+            segment0[0], segment0[1], segment0[2], segment1[0], segment1[1],
+            segment1[2]);
+    for (node = 0; node < n; node++) {
+      for (i = 0; i < 3 + ldim; i++) {
+        fprintf(file, " %.15e", xyzf[i + (3 + ldim) * order[node]]);
+      }
+      fprintf(file, "\n");
+    }
+    ref_free(order);
+    ref_free(t);
+    ref_free(xyzf);
+  }
+  return REF_SUCCESS;
+}
+
 REF_FCN REF_STATUS ref_iso_slice(REF_GRID *iso_grid, REF_GRID ref_grid,
                                  REF_DBL *normal, REF_DBL offset, REF_INT ldim,
                                  REF_DBL *in, REF_DBL **out) {

src/ref_iso.h

@@ -62,6 +62,10 @@ REF_FCN REF_STATUS ref_iso_boom_zone(FILE *file, REF_GRID ref_grid,
                                      REF_DBL *field, REF_INT ldim,
                                      REF_DBL *center, REF_DBL aoa, REF_DBL phi,
                                      REF_DBL h);
+REF_FCN REF_STATUS ref_iso_boomray(const char *filename, REF_GRID ref_grid,
+                                   REF_DBL *field, REF_INT ldim,
+                                   const char **scalar_names, REF_DBL *xyz0,
+                                   REF_DBL *xyz1);
 
 REF_FCN REF_STATUS ref_iso_slice(REF_GRID *iso_grid, REF_GRID ref_grid,
                                  REF_DBL *normal, REF_DBL offset, REF_INT ldim,

src/ref_subcommand.c

@@ -320,6 +320,10 @@ static void visualize_help(const char *name) {
   printf(
       "   --slice <nx> <ny> <nz> <offset> <slice.extension> "
       "extracts a slice.\n");
+  printf(
+      "   --boomray <x0> <y0> <z0> <x1> <y1> <z1> <ray.tec> "
+      "extracts a ray of dp/pinf defined by two points.\n");
+
   printf("\n");
 }
 
@@ -4560,6 +4564,50 @@ static REF_STATUS visualize(REF_MPI ref_mpi, int argc, char *argv[]) {
     return REF_SUCCESS;
   }
 
+  {
+    REF_BOOL boom_ray = REF_FALSE;
+    REF_DBL *dp_pinf = NULL;
+    for (pos = 0; pos < argc - 1; pos++) {
+      if (strcmp(argv[pos], "--boomray") == 0) {
+        REF_DBL xyz0[3], xyz1[3];
+
+        char *boomray_filename;
+        const char *vars[] = {"dp/pinf"};
+        RAS(pos < argc - 7,
+            "not enough arguments for --boomray <x0> <y0> <z0> <x1> <y1> <z1> "
+            "<ray.tec>");
+        boom_ray = REF_TRUE;
+        if (NULL == dp_pinf) {
+          REF_INT node;
+          ref_malloc(dp_pinf, ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
+          each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
+            REF_INT pressure_index = 4;
+            REF_DBL gamma = 1.4;
+            dp_pinf[node] =
+                (field[pressure_index + ldim * node] - 1.0 / gamma) * gamma;
+          }
+        }
+        xyz0[0] = atof(argv[pos + 1]);
+        xyz0[1] = atof(argv[pos + 2]);
+        xyz0[2] = atof(argv[pos + 3]);
+        xyz1[0] = atof(argv[pos + 4]);
+        xyz1[1] = atof(argv[pos + 5]);
+        xyz1[2] = atof(argv[pos + 6]);
+        boomray_filename = argv[pos + 7];
+        RSS(ref_iso_boomray(boomray_filename, ref_grid, dp_pinf, 1, vars, xyz0,
+                            xyz1),
+            "boomray");
+        ref_free(dp_pinf);
+      }
+    }
+    if (boom_ray) {
+      ref_free(dp_pinf);
+      ref_free(field);
+      RSS(ref_grid_free(ref_grid), "free grid");
+      return REF_SUCCESS;
+    }
+  }
+
   RXS(ref_args_find(argc, argv, "--subtract", &pos), REF_NOT_FOUND,
       "arg search");
   if (REF_EMPTY != pos && pos < argc - 1) {