Update generate, coordinates storing.

Author: rpep

fidimag/atomistic/demag.py

@@ -6,7 +6,7 @@
 from fidimag.atomistic.energy import Energy
 import fidimag.extensions.fmm as fmm
 import time
-
+import sys
 
 
 class Demag(Energy):
@@ -47,7 +47,7 @@ def __init__(self, calc_every=0, name='Demag'):
 
     def setup(self, mesh, spin, mu_s, mu_s_inv):
         super(Demag, self).setup(mesh, spin, mu_s, mu_s_inv)
-        
+
         # Ryan Pepper 04/04/2019
         # We *do not* need to scale by mesh.unit_length**3 here!
         # This is because in the base energy class, dx, dy and dz
@@ -102,7 +102,7 @@ def compute_energy(self):
         return energy / self.scale
 
 
-class DemagFMM(Energy): 
+class DemagFMM(Energy):
     def __init__(self, order, ncrit, theta, name="DemagFMM"):
         self.name = name
         assert order > 0, "Order must be 1 or higher"
@@ -116,14 +116,16 @@ def __init__(self, order, ncrit, theta, name="DemagFMM"):
     def setup(self, mesh, spin, mu_s, mu_s_inv):
         super(DemagFMM, self).setup(mesh, spin, mu_s, mu_s_inv)
         self.n = mesh.n
-        print(mesh.coordinates)
         self.m_temp = spin.copy()
         self.m_temp[0::3] *= self.mu_s
         self.m_temp[1::3] *= self.mu_s
         self.m_temp[2::3] *= self.mu_s
+        self.coords = mesh.coordinates * mesh.unit_length
+        print(np.min(self.coords[:, 0]))
+        print(np.max(self.coords[:, 0]))
         self.fmm = fmm.FMM(self.n, self.ncrit, self.theta,
                            self.order,
-                           mesh.coordinates * mesh.unit_length,
+                           self.coords,
                            self.m_temp)
 
     def compute_field(self, t=0, spin=None):
@@ -134,6 +136,6 @@ def compute_field(self, t=0, spin=None):
 
         self.field[:] = 0.0
         #self.fmm.set(self.m_temp)
-        self.fmm.compute_field(self.theta, self.field)
+        self.fmm.compute_field(self.field)
         self.field *= 1e-7
         return self.field

fidimag/atomistic/fmmlib/calculate.cpp

@@ -144,6 +144,7 @@ void interact_dehnen(size_t A, size_t B, std::vector<Cell> &cells, std::vector<P
   }
 }
 
+
 void interact_dehnen_lazy(const size_t A, const size_t B,
                           const std::vector<Cell> &cells,
                           const std::vector<Particle> &particles,
@@ -155,6 +156,7 @@ void interact_dehnen_lazy(const size_t A, const size_t B,
   const double dy = cells[A].y - cells[B].y;
   const double dz = cells[A].z - cells[B].z;
   const double R = sqrt(dx*dx + dy*dy + dz*dz);
+  // std::cout << "cells["<<A<<"].rmax = " << cells[A].rmax << "cells["<<B<<"].rmax = " << cells[B].rmax << std::endl;
 
   if (R*theta > (cells[A].rmax + cells[B].rmax)) {
     //if (cells[A].nleaf < ncrit && cells[B].nleaf < ncrit) {
@@ -305,7 +307,6 @@ void evaluate_approx_lazy(std::vector<Particle> &particles, std::vector<Cell> &c
       #pragma omp barrier
       evaluate_L2P(particles, cells, F, ncrit, order);
     }
-
 }
 
 
@@ -315,15 +316,17 @@ void build_interaction_lists(std::vector<std::pair<size_t, size_t>> &M2L_list,
                              std::vector<Particle> &particles,
                              double theta,
                              size_t order,
-                             size_t ncrit
-    ) {
+                             size_t ncrit) {
+
+
+    interact_dehnen_lazy(0, 0, cells, particles, theta, order, ncrit, M2L_list, P2P_list);
+
+    std::cout << "M2L_list.size() = " << M2L_list.size() << std::endl;
+    std::cout << "P2P_list.size() = " << P2P_list.size() << std::endl;
 
     std::sort(M2L_list.begin(), M2L_list.end(),
            [](std::pair<size_t, size_t> &left, std::pair<size_t, size_t> &right) {
                 return left.first < right.first;
                }
            );
-
-    std::cout << "M2L_list size = " << M2L_list.size() << std::endl;
-    std::cout << "P2P_list size = " << P2P_list.size() << std::endl;
-}
+}