diff --git a/apps/kida/kida.c b/apps/kida/kida.c
index 4b9d2486dfe34da0f4eca0ea1dbe2586451b8798..2ac43f1b4749f5c24fd1025cc5b9256393e73a94 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 d9d17445d1e9bc9ac966403f91be3c673f803001..558159a7c59c5d79bce29664b0e1f7857356ece6 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 fb9b732e8b2a6d8f06df59cae0ba8b052299093c..b644fc1f9f8773af6d6459f6e6ba0d4a1a9411d1 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 494b7c56364ace04454f57fc2dc7a6dfcee7e087..dded18761ccf0e118e842e0963e88f856ab4e9a2 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 09762b4df7dc4fc89ec87c9ccc4ea1699e7a0742..81d9805eae3987a9c97d27e4d4f19bc72b33678b 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 20d8830fc1c2e13d48abce20233d7eb077889489..44067505766feee6dd270a116a4acdd5b135f618 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 19ec1972a6b4b0504a1e0c1dee5f97044c83a20d..f3416aae4bbbc8eafe3ba314adb2bfdc433ca4f1 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