From 849a91473815c141080f3572dfb7d1457cea5e2a Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@visiblegeology.com>
Date: Sun, 4 Aug 2013 07:54:33 -0700
Subject: [PATCH] LOM grid ding and visualization code.

---
 SimPEG/LomView.py     | 61 +++++++++++++++++++++++++++++++++++++++++++
 SimPEG/exampleGrid.py | 25 ++++++++++++++++++
 2 files changed, 86 insertions(+)
 create mode 100644 SimPEG/LomView.py
 create mode 100644 SimPEG/exampleGrid.py

diff --git a/SimPEG/LomView.py b/SimPEG/LomView.py
new file mode 100644
index 00000000..1596dfd9
--- /dev/null
+++ b/SimPEG/LomView.py
@@ -0,0 +1,61 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+from mpl_toolkits.mplot3d import Axes3D
+from utils import mkvc
+
+
+class LomView(object):
+    """
+    Provides viewing functions for TensorMesh
+
+    This class is inherited by TensorMesh
+    """
+    def __init__(self):
+        pass
+
+    def plotGrid(self):
+        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions."""
+        NN = self.r(self.gridN, 'N', 'N', 'M')
+        if self.dim == 2:
+            fig = plt.figure(2)
+            fig.clf()
+            ax = plt.subplot(111)
+            X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
+            Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
+
+            X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
+            Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
+
+            X = np.r_[X1, X2]
+            Y = np.r_[Y1, Y2]
+
+            plt.plot(X, Y)
+        elif self.dim == 3:
+            fig = plt.figure(3)
+            fig.clf()
+            ax = fig.add_subplot(111, projection='3d')
+            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
+            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
+            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
+
+            X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
+            Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
+            Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
+
+            X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
+            Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
+            Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
+
+            X = np.r_[X1, X2, X3]
+            Y = np.r_[Y1, Y2, Y3]
+            Z = np.r_[Z1, Z2, Z3]
+
+            plt.plot(X, Y, 'b', zs=Z)
+            ax.set_zlabel('x3')
+
+        ax.grid(True)
+        ax.hold(False)
+        ax.set_xlabel('x1')
+        ax.set_ylabel('x2')
+        fig.show()
diff --git a/SimPEG/exampleGrid.py b/SimPEG/exampleGrid.py
new file mode 100644
index 00000000..a49326df
--- /dev/null
+++ b/SimPEG/exampleGrid.py
@@ -0,0 +1,25 @@
+import numpy as np
+from utils import mkvc, ndgrid
+
+
+def exampleLomGird(nC, exType):
+    assert type(nC) == list, "nC must be a list containing the number of nodes"
+    assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions"
+    exType = exType.lower()
+
+    possibleTypes = ['rect', 'rotate']
+    assert exType in possibleTypes, "Not a possible example type."
+
+    if exType == 'rect':
+        return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
+    elif exType == 'rotate':
+        if len(nC) == 2:
+            X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
+            amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2)
+            amt[amt < 0] = 0
+            return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt
+        elif len(nC) == 3:
+            X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
+            amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2)
+            amt[amt < 0] = 0
+            return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt