add wrappers for c arrays; star wrapping the tGrid struct

pyloc3d.py

@@ -3,9 +3,91 @@
 #  import os
 #  os.environ['LD_LIBRARY_PATH'] = "/ssd/du23lag/python/ahloc3d"
 
+
+from cffi import FFI
+import numpy as np
 from _ahloc_3d import ffi, lib
 from ctypes import create_string_buffer
 
+
+ffi = FFI()
+def ctype(ptr):
+    return ffi.getctype(ffi.typeof(ptr).item)
+
+def wrap_in_numpy_array(ptr, nn):
+    # taken from https://gist.github.com/yig/77667e676163bbfc6c44af02657618a6
+    T = ctype(ptr)
+    assert(T == 'double')
+    return np.frombuffer(ffi.buffer(ptr,nn*ffi.sizeof(T)), np.dtype('f8'))
+
+class Grid:
+    def __init__(self, g):
+        assert(ctype(g) == "tGrid")
+        self.ng = g.ng # grid number
+        self.nn = g.nn # nu*nv*nw
+        self.nu = g.nu
+        self.nv = g.nv
+        self.nw = g.nw
+        # local coordinates
+        self.umax = g.umax
+        self.umin = g.umin
+        self.vmax = g.vmax
+        self.vmin = g.vmin
+        self.wmax = g.wmax
+        self.wmin = g.wmin
+        # 3d cartesian coordinates
+        self.xcart = wrap_in_numpy_array(g.xcart, g.nn)
+        self.ycart = wrap_in_numpy_array(g.ycart, g.nn)
+        self.zcart = wrap_in_numpy_array(g.zcart, g.nn)
+        # 3d local coordinates
+        self.ucoord = wrap_in_numpy_array(g.ucoord, g.nn)
+        self.vcoord = wrap_in_numpy_array(g.vcoord, g.nn)
+        self.wcoord = wrap_in_numpy_array(g.wcoord, g.nn)
+        # 1d local coordinates
+        self.ucoord_1d = wrap_in_numpy_array(g.ucoord_1d, g.nu)
+        self.vcoord_1d = wrap_in_numpy_array(g.vcoord_1d, g.nv)
+        self.wcoord_1d = wrap_in_numpy_array(g.wcoord_1d, g.nw)
+        # 1d spectral weights
+        self.uweigh_1d = wrap_in_numpy_array(g.uweigh_1d, g.nu)
+        self.vweigh_1d = wrap_in_numpy_array(g.vweigh_1d, g.nv)
+        self.wweigh_1d = wrap_in_numpy_array(g.wweigh_1d, g.nw)
+        # Jacobian matrix
+        self.J_du_dx = wrap_in_numpy_array(g.J_du_dx, g.nn)
+        self.J_du_dy = wrap_in_numpy_array(g.J_du_dy, g.nn)
+        self.J_du_dz = wrap_in_numpy_array(g.J_du_dz, g.nn)
+        self.J_dv_dx = wrap_in_numpy_array(g.J_dv_dx, g.nn)
+        self.J_dv_dy = wrap_in_numpy_array(g.J_dv_dy, g.nn)
+        self.J_dv_dz = wrap_in_numpy_array(g.J_dv_dz, g.nn)
+        self.J_dw_dx = wrap_in_numpy_array(g.J_dw_dx, g.nn)
+        self.J_dw_dy = wrap_in_numpy_array(g.J_dw_dy, g.nn)
+        self.J_dw_dz = wrap_in_numpy_array(g.J_dw_dz, g.nn)
+        # derivative matrices
+        self.DM_u = wrap_in_numpy_array(g.DM_u, g.nu**2)
+        self.DM_v = wrap_in_numpy_array(g.DM_v, g.nv**2)
+        self.DM_w = wrap_in_numpy_array(g.DM_w, g.nw**2)
+
+        # TODO Reshape arrays
+        # TODO Check if dimensions are correct
+
+        # CC=0, CS=1, SS=2, SC=3, CU=4, DF=5
+        # grid type
+        self.type = g.type
+
+        # XP=0, YP=1, ZP=2, XM=3, YM=4, ZM=5, CENTER=6
+        # grid orientation
+        self.orientation = g.orientation
+
+        # fields
+        all_fields_cart = wrap_in_numpy_array(g.fields_cart, g.nn*lib.n_fields_cart)
+        self.fields_cart = np.split(all_fields_cart, lib.n_fields_cart)
+
+        # FULL3D=0, AXISYM=1
+        ndim = 3 if lib.sym_mode == 0 else 2
+        # derivatives of fields
+        all_derivatives = wrap_in_numpy_array(g.derivatives, g.nn*lib.n_derivatives*ndim)
+        self.derivatives = np.split(all_fields_cart, lib.n_derivatives*ndim)
+
+
 parfile = "./test.par"
 cstr_parfile = parfile.encode('utf8')
 lib.setup_output(cstr_parfile)
@@ -16,6 +98,12 @@ lib.stdout_verbose = 0
 #  print("stdout_verbose = ", lib.stdout_verbose)
 n_work = lib.n_work
 work = lib.work
+grids = []
 for i in range(n_work):
-    lib.load_data_and_derivatives(work[i])
+    w = work[i]
+    lib.load_data_and_derivatives(w)
+    ngrids = w.ngrids
+    for n in range(ngrids):
+        grids.append(Grid(w.grids[n]))
 
+print("wrapped number of grids = ", len(grids))