Added a class for nonspinning puncture data

AHFpy/Puncture.py

+import os.path
+import numpy as np
+import h5py
+
+from AHFpy.data import ab_indexes, ij_indexes
+from AHFpy.data import write_h5_dset
+
+import matplotlib.pyplot as plt
+
+class Puncture():
+    """
+    Puncture data in Cartesian coordinates
+    """
+    
+    def __init__(self,
+                 M = 1., # Mass parameter
+                 x = np.arange(-5+0.125/2,5,0.25/2), # better staggered!
+                 y = np.arange(-5+0.125/2,5,0.25/2),
+                 z = np.arange(-5+0.125/2,5,0.25/2)):
+
+        self.M = M
+
+        self.x = x
+        self.y = y
+        self.z = z
+        
+        self.grid = None        
+        self.grid = self.CartGrid()
+        
+    def CartGrid(self):
+        """
+        Cartesian mesh grid
+        """
+        #if not self.grid == None:
+        #    return self.grid
+        return np.array(np.meshgrid(self.x,self.y,self.z,
+                           indexing="ij"))
+        
+    def Psi4(self):
+        """
+        Conformal factor
+        """
+        x,y,z = self.CartGrid()
+        r = np.sqrt(x**2 + y**2 + z**2)
+        return (1+0.5*self.M/r)**4
+            
+    def FlatMetric(self):
+        """
+        Minkowski metric (-,+,+,+)
+        """
+        g = {}
+        for (a,b) in ab_indexes:
+            g[(a,b)] = +1 if a == b else 0.
+        g[(0,0)] = -1
+        return g
+    
+    def Metric(self):
+        """
+        Metric $g_{ab}$
+        """
+        eta = self.FlatMetric()
+        Psi4 = self.Psi4()
+        g = {}
+        for (a,b) in ab_indexes:
+            g[(a,b)] = Psi4 * eta[(a,b)]
+        return g
+
+    def ExtrinsicCurvature(self):
+        """
+        Extrinsic cuvrvature $K_{ij}$
+        """
+        K = {}
+        for (i,j) in ij_indexes:
+            K[(i,j)] = 0.
+        return K
+        
+    def Output(self,fname):
+        """
+        Write HDF5 data
+        """
+        g = self.Metric()
+        K = self.ExtrinsicCurvature()
+        points = self.x, self.y, self.z
+        return write_h5_dset(fname,
+                             time = 0, points = points,
+                             g = g, K = K)
+    
+    
+if __name__ == "__main__":
+
+
+    dfile = "puncture.h5"
+    p = Puncture()
+    p.Output(dfile)
+