Fixed bug in exchange class; problem was in the argument order in the test

Author: davidcortesortuno

fidimag/common/c_clib.pyx

@@ -238,7 +238,7 @@ def init_vector_func_fast(m0, mesh, double[:] field, norm=False, *args):
 #
 # - Methods that are redefined in inherited classes are only specified in the
 # base class
-# - These can be declared in a .pxd file: 
+# - These can be declared in a .pxd file:
 
 cdef extern from "c_energy.h":
 
@@ -250,27 +250,29 @@ cdef extern from "c_energy.h":
         void compute_field(double t)
         double compute_energy()
         void setup(int nx, int ny, int nz, double dx, double dy, double dz,
-                   double unit_length, double *spin,
-                   double *Ms, double *Ms_inv)
-
-        bool set_up
-        int nx, ny, nz, n
-        double dx, dy, dz
-        double unit_length
-        double *spin
-        double *Ms
-        double *Ms_inv
-        double *field
-        double *energy
-        double *coordinates
-        int *ngbs
+                   double unit_length, double *spin, double *Ms, double *Ms_inv,
+                   double *coordinates, int *ngbs,
+                   double *energy, double *field
+                   )
+
+        # bool set_up
+        # int nx, ny, nz, n
+        # double dx, dy, dz
+        # double unit_length
+        # double *spin
+        # double *Ms
+        # double *Ms_inv
+        # double *field
+        # double *energy
+        # double *coordinates
+        # int *ngbs
 
     cdef cppclass ExchangeEnergy(Energy):
         ExchangeEnergy() except +
 
         void init(double *A)
 
-        double *A
+        # double *A
 
 # cdef extern from "c_energy.cpp":
 #     pass
@@ -285,17 +287,17 @@ cdef extern from "c_energy.h":
 # __cinit__ and __dealloc__ methods which are guaranteed to be called exactly
 # once upon creation and deletion of the Python instance.
 
-cdef class PyEnergy:
-    cdef Energy *thisptr
-    # Try cinit:
-    def __cinit__(self):
-        # Should we allocate?
-        # self.thisptr = new ExchangeEnergy()
-        print("In Python A")
-
-# We need to inherit from a cdef-ined base class to make the inheritance to
-# work. This way the constructors are called properly
-cdef class PyExchangeEnergy(PyEnergy):
+#cdef class PyEnergy:
+#    cdef Energy *thisptr
+#    # Try cinit:
+#    def __cinit__(self):
+#        # Should we allocate?
+#        # self.thisptr = new Energy()
+#        print("In Python A")
+
+# No need to replicate C++ class structure; we make the ExchangeEnergy a
+# base class in Python (which we know inherits from Energy in C++)
+cdef class PyExchangeEnergy(object):
     cdef ExchangeEnergy *derivedptr
     # Try cinit:
     def __cinit__(self, double [:] A):
@@ -304,9 +306,8 @@ cdef class PyExchangeEnergy(PyEnergy):
         self.derivedptr = new ExchangeEnergy()
         self.derivedptr.init(&A[0])
 
-        if self.thisptr:
-            print("in B: deallocating old A")
-            del self.thisptr
+    def __dealloc__(self):
+        del self.derivedptr
 
     # DEBUG: check contents of the A array
     # def printA(self):
@@ -315,21 +316,20 @@ cdef class PyExchangeEnergy(PyEnergy):
     #         lst.append(self.derivedptr.A[i])
     #     print(lst)
 
-    # We could use another constructor if we use this method:
-    # def __cinit__(self):
-    #     if type(self) is PyExchangeEnergy:
-    #         self.thisptr = new ExchangeEnergy()
-    # def __dealloc__(self):
-    #     if type(self) is PyExchangeEnergy:
-    #         del self.thisptr
-
     # Necessary?:
     def compute_field(self, t):
         self.derivedptr.compute_field(t)
+
     def compute_energy(self, time):
         return self.derivedptr.compute_energy()
+
     def setup(self, nx, ny, nz, dx, dy, dz, unit_length,
-              double [:] spin, double [:] Ms, double [:] Ms_inv):
-        return self.derivedptr.setup(nx, ny, nz, dx, dy, dz, unit_length,
-                                     &spin[0], &Ms[0], &Ms_inv[0])
+              double [:] spin, double [:] Ms, double [:] Ms_inv,
+              double [:, :] coordinates, int [:, :] neighbours,
+              double [:] energy, double [:] field):
 
+        return self.derivedptr.setup(nx, ny, nz, dx, dy, dz, unit_length,
+                                     &spin[0], &Ms[0], &Ms_inv[0],
+                                     &coordinates[0, 0], &neighbours[0, 0],
+                                     &energy[0], &field[0]
+                                     )

native/include/c_energy.h

@@ -4,7 +4,7 @@
 class Energy {
 public:
     Energy() {std::cout << "In A; at " << this << "\n";};
-    ~Energy() {std::cout << "Killing A\n";};
+    virtual ~Energy() {std::cout << "Killing A\n";};
     bool set_up;
     int nx, ny, nz, n;
     double dx, dy, dz;
@@ -19,7 +19,7 @@ class Energy {
     double compute_energy();
     void setup(int nx, int ny, int nz, double dx, double dy, double dz,
                double unit_length, double *spin, double *Ms, double *Ms_inv,
-               double *coordinates, double *ngbs, 
+               double *coordinates, int *ngbs, 
                double *energy, double *field
                );
     virtual void compute_field(double t) {};
@@ -28,6 +28,7 @@ class Energy {
 class ExchangeEnergy : public Energy {
 public:
     ExchangeEnergy() {std::cout << "In B; at " << this << "\n";};
+    ~ExchangeEnergy() {std::cout << "Killing B\n";};
     void init(double *A) {
         this->set_up = false;
         this->A = A;

native/src/c_energy.cpp

@@ -12,12 +12,14 @@ double Energy::compute_energy() {
 
 void Energy::setup(int nx, int ny, int nz, double dx, double dy, double dz,
                    double unit_length, double *spin, double *Ms, double *Ms_inv,
-                   double *coordinates, double *ngbs, 
+                   double *coordinates, int *ngbs, 
                    double *energy, double *field
-                   ); {
+                   ) {
+
     this->nx = nx;
     this->ny = ny;
     this->nz = nz;
+    this->n = nx * ny * nz;
     this->dx = dx;
     this->dy = dy;
     this->dz = dz;
@@ -29,7 +31,7 @@ void Energy::setup(int nx, int ny, int nz, double dx, double dy, double dz,
     this->Ms_inv = Ms_inv;
     this->coordinates = coordinates;
     this->ngbs = ngbs;
-    this ->energy->energy;
+    this->energy = energy;
     this->field = field;
 
     set_up = true;
@@ -190,5 +192,6 @@ void ExchangeEnergy::compute_field(double t) {
         field[3 * i]     = fx * Ms_inv[i] * MU_0_INV;
         field[3 * i + 1] = fy * Ms_inv[i] * MU_0_INV;
         field[3 * i + 2] = fz * Ms_inv[i] * MU_0_INV;
+
     }
 }

tests/test_C_classes.py

@@ -1,6 +1,41 @@
 from fidimag.extensions.c_clib import PyExchangeEnergy
+from fidimag.common import CuboidMesh
+
+import time
 import numpy as np
 
 A = np.ones(9, dtype=np.double)
 print(A)
 Exch = PyExchangeEnergy(A)
+# A[2] = 4
+# Exch.printA()
+
+nx, ny, nz = 3, 3, 1
+n = nx * ny * nz
+dx, dy, dz = 1, 1, 1
+unit_length = 1.
+
+mesh = CuboidMesh(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)
+
+spin = np.ones(3 * n, dtype=np.double)
+spin[::3] = 0
+spin[1::3] = 0
+spin[3 * 4:3 * 4 + 3] = [0, 1, 0]
+Ms = np.ones(n)
+Ms_inv = 1 / Ms
+
+field = np.zeros(3 * n, dtype=np.double)
+energy = np.zeros(n, dtype=np.double)
+
+Exch.setup(mesh.nx, mesh.ny, mesh.nz, mesh.dx, mesh.dy, mesh.dz,
+           mesh.unit_length,
+           spin, Ms, Ms_inv,
+           mesh.coordinates, mesh.neighbours,
+           energy, field
+           )
+
+Exch.compute_field(0)
+print(field)
+print(energy)
+
+# time.sleep(5)