From 1b063325038b906b872c91b85d72b11e1636e918 Mon Sep 17 00:00:00 2001
From: Ryan Pepper <ryan.pepper@soton.ac.uk>
Date: Wed, 15 Jan 2020 21:40:57 +0000
Subject: [PATCH] Add atomistic demag to master

---
 doc/ipynb/steepest_descent_atomistic.ipynb |  2 +-
 fidimag/atomistic/demag.py                 | 35 ++++++++++++++++++++++
 2 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/doc/ipynb/steepest_descent_atomistic.ipynb b/doc/ipynb/steepest_descent_atomistic.ipynb
index 800a051a..190f6b82 100644
--- a/doc/ipynb/steepest_descent_atomistic.ipynb
+++ b/doc/ipynb/steepest_descent_atomistic.ipynb
@@ -1014,7 +1014,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.5"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,
diff --git a/fidimag/atomistic/demag.py b/fidimag/atomistic/demag.py
index a5928220..8f35d0cf 100644
--- a/fidimag/atomistic/demag.py
+++ b/fidimag/atomistic/demag.py
@@ -94,3 +94,38 @@ def compute_energy(self):
         self.energy /= self.scale
 
         return energy / self.scale
+
+
+class DemagFMM(Energy):
+    def __init__(self, order, ncrit, theta, name="DemagFMM", type='fmm'):
+        self.type = type
+        if self.type == 'fmm':
+            self._type = 0
+        elif self.type == 'bh':
+            self._type = 1
+
+        self.name = name
+        assert order > 0, "Order must be 1 or higher"
+        assert order < 11, "Order bust be < 11"
+        self.order = order
+        assert ncrit >= 2, "ncrit must be greater than 1."
+        self.ncrit = ncrit
+        assert theta >= 0.0, "theta must be >= 0.0"
+        self.theta = theta
+
+    def setup(self, mesh, spin, mu_s, mu_s_inv):
+        super(DemagFMM, self).setup(mesh, spin, mu_s, mu_s_inv)
+        self.coords = mesh.coordinates * mesh.unit_length
+        self.m_temp = np.zeros_like(spin)
+        self.fmm = fmm.FMM(self.n, self.ncrit, self.theta,
+                           self.order,
+                           self.coords,
+                           self.m_temp,
+                           self.mu_s, self._type
+                           )
+
+    def compute_field(self, t=0, spin=None):
+        self.m_temp[:] = spin if spin is not None else self.spin
+        self.fmm.compute_field(self.field)
+        self.field *= -1e-7
+        return self.field