Merge branch 'develop' of github.com:ifilot/pytessel

.gitignore

@@ -149,3 +149,4 @@ test.py
 nose_debug
 test.ply
 .vscode
+*.ply

example/build_metaballs_icosahedron_rectangular_grid.py

@@ -0,0 +1,60 @@
+from pytessel import PyTessel
+import numpy as np
+
+def main():
+    """
+    Build 6 .ply files of increasing quality
+    """
+    pytessel = PyTessel()
+
+    sz = 3
+
+    x = np.linspace(-sz, sz, 30)
+    y = np.linspace(-sz, sz, 40)
+    z = np.linspace(-sz, sz, 50)
+
+    xx, yy, zz, field = icosahedron_field(x,y,z)
+    print(field.shape)
+
+    unitcell = np.diag(np.ones(3) * sz * 2)
+    pytessel = PyTessel()
+    isovalue = 3.75
+    vertices, normals, indices = pytessel.marching_cubes(field.flatten(), field.shape, unitcell.flatten(), isovalue)
+
+    pytessel.write_ply('icosahedron_metaballs.ply', vertices, normals, indices)
+
+def icosahedron_field(x,y,z):
+    """
+    Produce a scalar field for the icosahedral metaballs
+    """
+    phi = (1 + np.sqrt(5)) / 2
+    vertices = [
+        [0,1,phi],
+        [0,-1,-phi],
+        [0,1,-phi],
+        [0,-1,phi],
+        [1,phi,0],
+        [-1,-phi,0],
+        [1,-phi,0],
+        [-1,phi,0],
+        [phi,0,1],
+        [-phi,0,-1],
+        [phi,0,-1],
+        [-phi,0,1]
+    ]
+
+    zz,yy,xx = np.meshgrid(z,y,x, indexing='ij')
+    field = np.zeros_like(xx)
+    for v in vertices:
+        field += f(xx,yy,zz,v[0], v[1],v[2])
+
+    return xx, yy, zz, field
+
+def f(x,y,z,X0,Y0,Z0):
+    """
+    Single metaball function
+    """
+    return 1 / ((x - X0)**2 + (y - Y0)**2 + (z - Z0)**2)
+
+if __name__ == '__main__':
+    main()

pytessel/scalar_field.cpp

@@ -60,6 +60,11 @@ float ScalarField::get_value_interp(float x, float y, float z) const {
     // cast the input to the nearest grid point
     Vec3 r = this->realspace_to_grid(x,y,z);
 
+    // recast
+    if(r.x < 0.0) r.x += (float)this->grid_dimensions[0];
+    if(r.y < 0.0) r.y += (float)this->grid_dimensions[1];
+    if(r.z < 0.0) r.z += (float)this->grid_dimensions[2];
+
     // calculate value using trilinear interpolation
     float xd = remainderf(r.x, 1.0);
     float yd = remainderf(r.y, 1.0);
@@ -69,12 +74,12 @@ float ScalarField::get_value_interp(float x, float y, float z) const {
     if(yd < 0) yd += 1.0;
     if(zd < 0) zd += 1.0;
 
-    float x0 = floor(r.x);
-    float x1 = ceil(r.x);
-    float y0 = floor(r.y);
-    float y1 = ceil(r.y);
-    float z0 = floor(r.z);
-    float z1 = ceil(r.z);
+    float x0 = fmod(floor(r.x), this->grid_dimensions[0]);
+    float x1 = fmod(ceil(r.x), this->grid_dimensions[0]);
+    float y0 = fmod(floor(r.y), this->grid_dimensions[1]);
+    float y1 = fmod(ceil(r.y), this->grid_dimensions[1]);
+    float z0 = fmod(floor(r.z), this->grid_dimensions[2]);
+    float z1 = fmod(ceil(r.z), this->grid_dimensions[2]);
 
     return
     this->get_value(x0, y0, z0) * (1.0 - xd) * (1.0 - yd) * (1.0 - zd) +
@@ -169,9 +174,9 @@ Vec3 ScalarField::realspace_to_direct(float x, float y, float z) const {
 Vec3 ScalarField::realspace_to_grid(float i, float j, float k) const {
     Vec3 g = this->realspace_to_direct(i, j, k);
 
-    g.x *= float(this->grid_dimensions[0]-1);
-    g.y *= float(this->grid_dimensions[1]-1);
-    g.z *= float(this->grid_dimensions[2]-1);
+    g.x *= float(this->grid_dimensions[0]);
+    g.y *= float(this->grid_dimensions[1]);
+    g.z *= float(this->grid_dimensions[2]);
 
     return g;
 }
@@ -208,3 +213,27 @@ void ScalarField::inverse(const mat33& mat, mat33* invmat) {
     (*invmat)[2][1] = (mat[2][0] * mat[0][1] - mat[0][0] * mat[2][1]) * invdet;
     (*invmat)[2][2] = (mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]) * invdet;
 }
+
+std::vector<float> ScalarField::get_unitcell_vf() const {
+    std::vector<float> ans(9,0.0);
+
+    for(unsigned int i=0; i<3; i++) {
+        for(unsigned int j=0; j<3; j++) {
+            ans[i*3 + j] = this->unitcell[i][j];
+        }
+    }
+
+    return ans;
+}
+
+std::vector<float> ScalarField::get_unitcell_inverse() const {
+    std::vector<float> ans(9,0.0);
+
+    for(unsigned int i=0; i<3; i++) {
+        for(unsigned int j=0; j<3; j++) {
+            ans[i*3 + j] = this->unitcell_inverse[i][j];
+        }
+    }
+
+    return ans;
+}

pytessel/scalar_field.h

@@ -92,5 +92,13 @@ class ScalarField{
         return this->unitcell;
     }
 
+    inline const mat33& get_mat_inverse() const {
+        return this->unitcell_inverse;
+    }
+
+    std::vector<float> get_unitcell_vf() const;
+
+    std::vector<float> get_unitcell_inverse() const;
+
 private:
 };

setup.py

@@ -59,14 +59,12 @@ def find_windows_versions():
         os.environ['PATH'] = ";".join(newpaths)
 
 if os.name == 'posix' and sys.platform != 'darwin':
-    os.environ['CFLAGS'] = '-I/usr/include/glm'
     extra_compile_args = ["-Wno-date-time", "-fopenmp", "-fPIC"]
     extra_link_args = ["-fopenmp"]
 elif os.name == 'nt':
     extra_compile_args = ["/wd4244"]
     extra_link_args = []
 elif sys.platform == 'darwin':
-    os.environ['CFLAGS'] = '-I/usr/local/Cellar/boost/1.81.0_1/include -I /usr/local/Cellar/glm/0.9.9.8/include'
     extra_compile_args = ["-Wno-date-time", "-fPIC", "-std=c++14"]
     extra_link_args = []