From 4c1dd61aada288cdccc9632ca83dcd7f79ed850a Mon Sep 17 00:00:00 2001
From: Orestis <orestis.malaspinas@pm.me>
Date: Sun, 21 Feb 2021 20:36:18 +0100
Subject: [PATCH] added mtetras for isosurface stuff

---
 apps/kida/kida.c                 | 14 ++++++++++++--
 src/c_interface/triangle_field.c | 17 ++++++++++++-----
 src/c_interface/triangle_field.h |  4 +++-
 src/futhark.pkg                  |  2 +-
 src/liblbm_f32.fut               |  4 +++-
 src/liblbm_f64.fut               |  4 +++-
 src/polygonise.fut               | 24 ++++++++++++++++++++++--
 7 files changed, 56 insertions(+), 13 deletions(-)

diff --git a/apps/kida/kida.c b/apps/kida/kida.c
index 4b9d248..2ac43f1 100644
--- a/apps/kida/kida.c
+++ b/apps/kida/kida.c
@@ -109,14 +109,24 @@ int main(int argc, char **argv) {
             vtk_file_close(&vtk);
             free(fname);
 
-            triangle_field tri = triangle_field_polygonise(density, ctx, 1.0);
+            triangle_field tri = triangle_field_polygonise(density, ctx, 1.0, mtetras); // mtetras or mcubes
 
-            char *stl_fname = stl_create_fname(prefix, "iso_surface", it);
+            char *stl_fname = stl_create_fname(prefix, "iso_surface_tetras", it);
             stl_file stl = stl_file_create(stl_fname);
             stl_file_write_triangle_field(&stl, tri);
             stl_file_close(&stl);
             free(stl_fname);
             triangle_field_deallocate(tri);
+
+            tri = triangle_field_polygonise(density, ctx, 1.0, mcubes); // mtetras or mcubes
+
+            stl_fname = stl_create_fname(prefix, "iso_surface_cubes", it);
+            stl = stl_file_create(stl_fname);
+            stl_file_write_triangle_field(&stl, tri);
+            stl_file_close(&stl);
+            free(stl_fname);
+            triangle_field_deallocate(tri);
+
         }
 
         f = lbm_collide_and_stream(ctx, &f, vel_ini, d, units.omega, delta_iter);
diff --git a/src/c_interface/triangle_field.c b/src/c_interface/triangle_field.c
index d9d1744..558159a 100644
--- a/src/c_interface/triangle_field.c
+++ b/src/c_interface/triangle_field.c
@@ -56,16 +56,23 @@ triangle_field triangle_field_from_futhark_f32_2d(futhark_2d *tri, struct futhar
     return field;
 }
 
-triangle_field triangle_field_polygonise(scalar_field rho, struct futhark_context *ctx, T isovalue) {
-    scalar_field_write_file(rho, "./tmp/", "rho", 0);
-    printf("sizeof %ld %d %d %d\n", sizeof(*rho.data), rho.nx, rho.ny, rho.nz);
-
+triangle_field triangle_field_polygonise(scalar_field rho, struct futhark_context *ctx, T isovalue, poly_method pm ) {
     futhark_3d *field = scalar_field_to_futhark_f32_3d(rho, ctx);
     assert(NULL != field && "Problem with rho to field.");
 
     futhark_2d *triangles = NULL;
     int64_t out0 = -1;
-    futhark_entry_main_polygonise(ctx, &out0, &triangles, field, isovalue);
+    switch (pm) {
+        case mcubes:
+            futhark_entry_main_mcubes_polygonise(ctx, &out0, &triangles, field, isovalue);        
+            break;
+        case mtetras:
+            futhark_entry_main_mtetras_polygonise(ctx, &out0, &triangles, field, isovalue);        
+            break;
+        default:
+            assert(false && "Should not have a default here.");
+    }
+
     char *str = futhark_context_get_error(ctx);
     if (str != NULL) {
         printf("%s\n", str);
diff --git a/src/c_interface/triangle_field.h b/src/c_interface/triangle_field.h
index fb9b732..b644fc1 100644
--- a/src/c_interface/triangle_field.h
+++ b/src/c_interface/triangle_field.h
@@ -4,6 +4,8 @@
 #include "scalar_field.h"
 #include "triangle.h"
 
+typedef enum { mcubes, mtetras } poly_method;
+
 typedef struct {
     fint num_triangles;
     triangle *t;
@@ -19,6 +21,6 @@ void triangle_field_from_futhark_f32_2d_inplace(futhark_2d *tri, struct futhark_
 
 triangle_field triangle_field_from_futhark_f32_2d(futhark_2d *tri, struct futhark_context *ctx);
 
-triangle_field triangle_field_polygonise(scalar_field rho, struct futhark_context *ctx, T isovalue);
+triangle_field triangle_field_polygonise(scalar_field rho, struct futhark_context *ctx, T isovalue, poly_method pm);
 
 #endif
\ No newline at end of file
diff --git a/src/futhark.pkg b/src/futhark.pkg
index 494b7c5..dded187 100644
--- a/src/futhark.pkg
+++ b/src/futhark.pkg
@@ -1,4 +1,4 @@
 require {
   github.com/athas/vector 0.7.0 #f5aaa219a02949d1eb2978f81419cf5614abfa3d
-  github.com/omalaspinas/futcubes 0.0.0-20210220002337+2055ab4b878406e92ef6eeea32b210ab429b0050 #2055ab4b878406e92ef6eeea32b210ab429b0050
+  github.com/omalaspinas/futcubes 0.0.0-20210221184312+bf89c12957be3be7d5f6c5a4efe75bef260654a5 #bf89c12957be3be7d5f6c5a4efe75bef260654a5
 }
diff --git a/src/liblbm_f32.fut b/src/liblbm_f32.fut
index 09762b4..81d9805 100644
--- a/src/liblbm_f32.fut
+++ b/src/liblbm_f32.fut
@@ -11,7 +11,9 @@ module sim_d3q19_reg_bfl_f32 = mk_sim d3q19_reg_bfl_f32
 
 entry main_d3q19_reg = sim_d3q19_reg_bfl_f32.run
 
-entry main_polygonise = polygonise_f32
+entry main_mcubes_polygonise = mcubes_polygonise_f32
+
+entry main_mtetras_polygonise = mtetras_polygonise_f32
 
 -- Test the sub1 function.
 -- ==
diff --git a/src/liblbm_f64.fut b/src/liblbm_f64.fut
index 20d8830..4406750 100644
--- a/src/liblbm_f64.fut
+++ b/src/liblbm_f64.fut
@@ -11,7 +11,9 @@ module sim_d3q19_reg_bfl_f64 = mk_sim d3q19_reg_bfl_f64
 
 entry main_d3q19_reg_f64 = sim_d3q19_reg_bfl_f64.run
 
-entry main_polygonise = polygonise_f64
+entry main_mcubes_polygonise = mcubes_polygonise_f64
+
+entry main_mtetras_polygonise = mtetras_polygonise_f64
 
 
 -- -- executes Lattice-Boltzmann it times and
diff --git a/src/polygonise.fut b/src/polygonise.fut
index 19ec197..f3416aa 100644
--- a/src/polygonise.fut
+++ b/src/polygonise.fut
@@ -1,9 +1,10 @@
 import "lib/github.com/omalaspinas/futcubes/mcubes"
+import "lib/github.com/omalaspinas/futcubes/mtetras"
 
-let polygonise_f64 [nx][ny][nz] (rho: [nx][ny][nz]f64) (isolevel: f64) =
+let mcubes_polygonise_f64 [nx][ny][nz] (rho: [nx][ny][nz]f64) (isolevel: f64) =
     map (mcubes.triangle_to_array) (mcubes.polygonise_field rho isolevel)
 
-let polygonise_f32 [nx][ny][nz] (rho: [nx][ny][nz]f32) (isolevel: f32) =
+let mcubes_polygonise_f32 [nx][ny][nz] (rho: [nx][ny][nz]f32) (isolevel: f32) =
     let rho' = 
         map (\rx -> 
             map (\rxy -> 
@@ -14,6 +15,25 @@ let polygonise_f32 [nx][ny][nz] (rho: [nx][ny][nz]f32) (isolevel: f32) =
         ) rho
     let tri = map (mcubes.triangle_to_array) (mcubes.polygonise_field rho' (f64.f32 isolevel))
     let si = length tri
+    in 
+        (si , map (\tx -> 
+            map(\txi -> f32.f64 txi) tx
+        ) tri ) 
+
+let mtetras_polygonise_f64 [nx][ny][nz] (rho: [nx][ny][nz]f64) (isolevel: f64) =
+    map (mtetras.triangle_to_array) (mtetras.polygonise_field rho isolevel)
+
+let mtetras_polygonise_f32 [nx][ny][nz] (rho: [nx][ny][nz]f32) (isolevel: f32) =
+    let rho' = 
+        map (\rx -> 
+            map (\rxy -> 
+                map (\rxyz -> 
+                    f64.f32 rxyz
+                ) rxy
+            ) rx 
+        ) rho
+    let tri = map (mtetras.triangle_to_array) (mtetras.polygonise_field rho' (f64.f32 isolevel))
+    let si = length tri
     in 
         (si , map (\tx -> 
             map(\txi -> f32.f64 txi) tx
-- 
GitLab