Source updates to make python-toast compile in Windows

arch.h

@@ -5,7 +5,7 @@
 #ifndef __ARCH_H
 #define __ARCH_H
 
-#if defined(WIN32)||defined(WIN64)
+#if defined(_WIN32)
 
 #include <process.h>
 #include <direct.h>
@@ -18,7 +18,7 @@
 #define isnan(arg) _isnan((arg))
 
 // avoid annoying warnings
-#define hypot(x,y) _hypot((x),(y))
+#define hypot _hypot
 #define isnan(x) _isnan((x))
 #define getpid() _getpid()
 #define getcwd(buffer,maxlen) _getcwd(buffer,maxlen)
@@ -27,6 +27,6 @@
 inline double drand48(void) { return (double)rand()/(double)RAND_MAX; }
 // quick fix of missing function. this should be improved upon!
 
-#endif
+#endif //_WIN32
 
 #endif // !__ARCH_H

script/python/tutorials/fwd2.py

@@ -1,78 +1,78 @@
-# This pytoast example solves the forward problem
-# for a homogeneous 2D circular problem
-
-# Import various modules
-import os
-import numpy as np
-from numpy import matrix
-from scipy import sparse
-from scipy.sparse import linalg
-from numpy.random import rand
-import matplotlib.pyplot as plt
-import matplotlib.animation as animation
-from matplotlib import ion
-
-ion()
-
-# PyToast environment
-execfile(os.getenv("TOASTDIR") + "/ptoast_install.py")
-
-# Import the toast modules
-from toast import mesh
-from toast import raster
-
-# Set the file paths
-meshdir = os.path.expandvars("$TOASTDIR/test/2D/meshes/")
-meshfile = meshdir + "circle25_32.msh"
-qmfile = meshdir + "circle25_32x32.qm"
-
-# Load the mesh and source/detector specs
-hmesh = mesh.Read(meshfile)
-mesh.ReadQM(hmesh,qmfile)
-
-# Extract mesh geometry
-nlist,elist,eltp = mesh.Data (hmesh)
-nlen = nlist.shape[0]
-
-grd = np.array([64,64])
-hraster = raster.Make (hmesh,grd);
-
-# Homogeneous parameter distributions
-mua = np.ones ((1,nlen)) * 0.025
-mus = np.ones ((1,nlen)) * 2.0
-ref = np.ones ((1,nlen)) * 1.4
-freq = 100
-
-# Set up the linear system
-qvec = mesh.Qvec (hmesh)
-mvec = mesh.Mvec (hmesh)
-nq = qvec.shape[0]
-
-phi = mesh.Fields(hmesh,hraster,qvec,mvec,mua,mus,ref,freq)
-
-fig1 = plt.figure()
-ims = []
-for q in range(nq):
-    phi = np.log (phi[q,:])
-    bphi = raster.MapBasis (hraster, 'S->B', phi)
-    bphi = np.reshape (bphi, grd)
-    ims.append ((plt.imshow(bphi.real),))
-
-im_ani = animation.ArtistAnimation(fig1, ims, interval=50, repeat_delay=3000, blit=True)
-
-plt.show()
-
-
-# Display fields
-plt.subplot(1,2,1)
-im = plt.imshow(bphi.real)
-plt.title('log amplitude')
-plt.colorbar()
-#plt.draw()
-#plt.show()
-plt.subplot(1,2,2)
-im = plt.imshow(bphi.imag)
-plt.title('phase')
-plt.colorbar()
-plt.show()
-
+# This pytoast example solves the forward problem
+# for a homogeneous 2D circular problem
+
+# Import various modules
+import os
+import numpy as np
+from numpy import matrix
+from scipy import sparse
+from scipy.sparse import linalg
+from numpy.random import rand
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+#from matplotlib import ion
+
+plt.ion()
+
+# PyToast environment
+execfile(os.getenv("TOASTDIR") + "/ptoast_install.py")
+
+# Import the toast modules
+from toast import mesh
+from toast import raster
+
+# Set the file paths
+meshdir = os.path.expandvars("$TOASTDIR/test/2D/meshes/")
+meshfile = meshdir + "circle25_32.msh"
+qmfile = meshdir + "circle25_32x32.qm"
+
+# Load the mesh and source/detector specs
+hmesh = mesh.Read(meshfile)
+mesh.ReadQM(hmesh,qmfile)
+
+# Extract mesh geometry
+nlist,elist,eltp = mesh.Data (hmesh)
+nlen = nlist.shape[0]
+
+grd = np.array([64,64])
+hraster = raster.Make (hmesh,grd);
+
+# Homogeneous parameter distributions
+mua = np.ones ((1,nlen)) * 0.025
+mus = np.ones ((1,nlen)) * 2.0
+ref = np.ones ((1,nlen)) * 1.4
+freq = 100
+
+# Set up the linear system
+qvec = mesh.Qvec (hmesh)
+mvec = mesh.Mvec (hmesh)
+nq = qvec.shape[0]
+
+phi = mesh.Fields(hmesh,hraster,qvec,mvec,mua,mus,ref,freq)
+
+fig1 = plt.figure()
+ims = []
+for q in range(nq):
+    lphi = np.log (phi[q,:])
+    bphi = raster.MapBasis (hraster, 'S->B', lphi)
+    bphi = np.reshape (bphi, grd)
+    ims.append ((plt.imshow(bphi.real),))
+
+im_ani = animation.ArtistAnimation(fig1, ims, interval=50, repeat_delay=3000, blit=True)
+
+plt.show()
+
+
+# Display fields
+plt.subplot(1,2,1)
+im = plt.imshow(bphi.real)
+plt.title('log amplitude')
+plt.colorbar()
+#plt.draw()
+#plt.show()
+plt.subplot(1,2,2)
+im = plt.imshow(bphi.imag)
+plt.title('phase')
+plt.colorbar()
+plt.show()
+

src/libfe/ellist.cc

@@ -104,7 +104,7 @@ Element *ElementList::SideNeighbour (int el, int side)
     return pel->sdnbhr[side];
 }
 
-int ElementList::EdgeAdjacentElement (int el, int side)
+int ElementList::EdgeAdjacentElement (int el, int side) const
 {
     const int max_sidenode = 10;
     int i, i1, i2, nn, nd[max_sidenode];

src/libfe/ellist.h

@@ -69,7 +69,7 @@ class FELIB ElementList {
     // Mesh::PopulateNeighbourLists
     // Note: This replaces 'EdgeAdjacentElement'
 
-    int EdgeAdjacentElement (int el, int side);
+    int EdgeAdjacentElement (int el, int side) const;
     // returns the number of the element adjacent to 'el' at side 'side'
     // returns -1 if no element is found, i.e. if 'side' is a boundary
 

src/libfe/mesh.cc

@@ -1038,7 +1038,7 @@ void Mesh::SetupElementMatrices (void)
 */
 // ***********************************************
 
-void Mesh::SetupNeighbourList ()
+void Mesh::SetupNeighbourList () const
 {
     int el, side;
     if (nbhrs) return;   // list exists already

src/libfe/mesh.h

@@ -80,7 +80,7 @@ class FELIB Mesh {
     ElementList elist;   // list of mesh elements
     ParameterList plist; // list of mesh parameters (per node)
 
-    int **nbhrs;
+    mutable int **nbhrs;
     // list of edge-adjacent neighbours for each element
     // this is only defined after a call to SetupNeighbourList
 
@@ -229,7 +229,7 @@ class FELIB Mesh {
 	int &side, double sub = 0) const;
     // This contains the old algorithm which is used whenever the above fails
 
-    void SetupNeighbourList ();
+    void SetupNeighbourList () const;
     // sets up a list of edge-adjacent elements for each element in the mesh
     // and stores it in nbhrs
 

src/libfe/pix4.cc

@@ -579,6 +579,48 @@ double Pixel4::BndIntPFF (int i, int j, const RVector &P) const
     return res;
 }
 
+double Pixel4::Size() const
+{
+	return size;
+}
+
+double Pixel4::IntF (int i) const
+{
+	return 0.25*dx*dy;
+}
+
+double Pixel4::IntFF (int i, int j) const
+{
+	return intff(i,j);
+}
+
+RSymMatrix Pixel4::IntFF() const
+{
+	return intff;
+}
+
+RSymMatrix Pixel4::IntDD () const
+{
+	return intdd;
+}
+
+double Pixel4::IntDD (int i, int j) const
+{
+	return intdd(i,j);
+}
+
+double Pixel4::IntFDD (int i, int j, int k) const
+{
+	return intfdd[i](j,k);
+}
+
+double Pixel4::IntPDD (int i, int j, const RVector &P) const
+{
+	double res = 0.0;
+    for (int k = 0; k < 4; k++) res += P[Node[k]] * intfdd[k](i,j);
+    return res;
+}
+
 double Pixel4::dx = 0.0;
 double Pixel4::dy = 0.0;
 double Pixel4::size = 0.0;

src/libfe/pix4.h

@@ -50,7 +50,7 @@ class FELIB Pixel4: public Element_Structured_2D {
     inline int nSideNode (int /*side*/) const { return 2; }
     int SideNode (int side, int node) const;
 
-    inline double Size() const { return size; }
+    double Size() const;
 
     Point Local (const NodeList &nlist, const Point &glob) const
     { return Local (glob); }
@@ -73,23 +73,24 @@ class FELIB Pixel4: public Element_Structured_2D {
     RDenseMatrix GlobalShapeD (const NodeList &nlist, const Point &glob) const
     { return GlobalShapeD (glob); }
 
-    inline double IntF (int i) const
-    { return 0.25*dx*dy; }
-    inline double IntFF (int i, int j) const
-    { return intff(i,j); }
-    inline RSymMatrix IntFF() const
-    { return intff; }
+	double IntF (int i) const;
+
+	double IntFF (int i, int j) const;
+
+    RSymMatrix IntFF() const;
+
     double IntFFF (int i, int j, int k) const;
     double IntPFF (int i, int j, const RVector &P) const;
     RSymMatrix IntPFF (const RVector& P) const;
-    inline RSymMatrix IntDD () const { return intdd; }
-    inline double IntDD (int i, int j) const { return intdd(i,j); }
-    inline double IntFDD (int i, int j, int k) const { return intfdd[i](j,k); }
-    inline double IntPDD (int i, int j, const RVector &P) const
-    { double res = 0.0;
-      for (int k = 0; k < 4; k++) res += P[Node[k]] * intfdd[k](i,j);
-      return res;
-    }
+
+    RSymMatrix IntDD () const;
+
+    double IntDD (int i, int j) const;
+
+    double IntFDD (int i, int j, int k) const;
+
+    double IntPDD (int i, int j, const RVector &P) const;
+
     inline RSymMatrix IntPDD (const RVector& P) const
     { RSymMatrix pdd(4);
       for (int i = 0; i < 4; i++)

src/libfe/vox8.cc

@@ -175,6 +175,46 @@ RDenseMatrix Voxel8::LocalShapeD (const Point &loc) const
     return der;
 }
 
+double Voxel8::Size() const
+{
+	return size;
+}
+
+double Voxel8::IntF (int i) const
+{
+	return intf;
+}
+
+double Voxel8::IntFF (int i, int j) const
+{
+	return intff(i,j);
+}
+
+RSymMatrix Voxel8::IntFF() const
+{
+	return intff;
+}
+
+double Voxel8::IntFFF (int i, int j, int k) const
+{
+	return intfff[i](j,k);
+}
+
+RSymMatrix Voxel8::IntDD () const
+{
+	return intdd;
+}
+
+double Voxel8::IntDD (int i, int j) const
+{
+	return intdd(i,j);
+}
+
+double Voxel8::IntFDD (int i, int j, int k) const
+{
+	return intfdd[i](j,k);
+}
+
 double Voxel8::IntPFF (int i, int j, const RVector &P) const
 {
     return (intfff[0](i,j) * P[Node[0]] +