Index: /issm/trunk-jpl/src/m/mesh/squaremesh.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/squaremesh.m	(revision 28135)
+++ /issm/trunk-jpl/src/m/mesh/squaremesh.m	(revision 28136)
@@ -1,3 +1,3 @@
-function md=squaremesh(md,Lx,Ly,nx,ny,computeconnectivity)
+function md=squaremesh(md,Lx,Ly,nx,ny,computeconnectivity,orientation)
 %SQUAREMESH - create a structured square mesh 
 %
@@ -11,6 +11,9 @@
 
 %process options
-if nargin == 5,
+if nargin == 5
 	computeconnectivity = 1;
+	orientation = 1;
+elseif nargin==6
+	orientation = 1;
 end
 
@@ -19,62 +22,29 @@
 nods=nx*ny;
 
-%Old method
-if 0,
-	%initialization
-	index=zeros(nel,3);
-	x=zeros(nx*ny,1);
-	y=zeros(nx*ny,1);
+%prepare coordinates of vertices
+x = repmat(linspace(0,Lx,nx),[ny 1]);
+x = reshape(x,[nx*ny 1]);
+y = repmat(linspace(0,Ly,ny)',[1 nx]);
+y = reshape(y,[nx*ny 1]);
 
-	%create coordinates
-	for n=1:nx,
-		for m=1:ny,
-			x((n-1)*ny+m)=(n-1.);
-			y((n-1)*ny+m)=(m-1.);
-		end
-	end
+%do first column of elements first
+nels1 = 2*(ny-1);
+index = ones(nels1,3);
+%First column
+index(2:2:nels1,1) = 2:ny;
+index(3:2:nels1,1) = 2:ny-1;
+%2d column
+index(1:2:nels1,2) = ny+1:2*ny-1;
+index(2:2:nels1,2) = ny+1:2*ny-1;
+%3rd column
+index(1:2:nels1,3) = 2:ny;
+index(2:2:nels1,3) = ny+2:2*ny;
 
-	%create index
-	for n=1:(nx-1)
-		for m=1:(ny-1),
-			A=(n-1)*ny+m;
-			B=A+1;
-			C=n*ny+m;
-			D=C+1;
-			index((n-1)*(ny-1)*2+2*(m-1)+1,:)=[A C B];
-			index((n-1)*(ny-1)*2+2*m,:)=[B C D];
-		end
-	end
-
-	%Scale  x and y
-	x=x/max(x)*Lx;
-	y=y/max(y)*Ly;
-else
-	%New method (faster!)
-	x = repmat(linspace(0,Lx,nx),[ny 1]);
-	x = reshape(x,[nx*ny 1]);
-	y = repmat(linspace(0,Ly,ny)',[1 nx]);
-	y = reshape(y,[nx*ny 1]);
-
-	%do first column of elements first
-	nels1 = 2*(ny-1);
-	index = ones(nels1,3);
-	%First column
-	index(2:2:nels1,1) = 2:ny;
-	index(3:2:nels1,1) = 2:ny-1;
-	%2d column
-	index(1:2:nels1,2) = ny+1:2*ny-1;
-	index(2:2:nels1,2) = ny+1:2*ny-1;
-	%3rd column
-	index(1:2:nels1,3) = 2:ny;
-	index(2:2:nels1,3) = ny+2:2*ny;
-
-	%Now copy column and offset with ny, nx times
-	index = repmat(index,[nx-1 1]);
-	offset = repmat([0:ny:(nx-2)*ny],[nels1 1]);
-	offset = reshape(offset,[(nx-1)*nels1,1]);
-	offset = repmat(offset,[1,3]);
-	index = index + offset;
-end
-
+%Now copy column and offset with ny, nx times
+index = repmat(index,[nx-1 1]);
+offset = repmat([0:ny:(nx-2)*ny],[nels1 1]);
+offset = reshape(offset,[(nx-1)*nels1,1]);
+offset = repmat(offset,[1,3]);
+index = index + offset;
 
 %create segments
@@ -88,4 +58,21 @@
 %back edge
 segments(2*(ny-1)+(nx-1)+1:2*(nx-1)+2*(ny-1),:)=[[1:ny:(nx-2)*ny+1]' [ny+1:ny:ny*(nx-1)+1]' [1:2*(ny-1):2*(nx-2)*(ny-1)+1]'];
+
+%Do we need to change the orientation of the diagonals?
+if orientation==1
+	%nothing to do
+elseif orientation==2
+	%switch diagonals
+	indexold = index;
+	index(1:2:end,3) = indexold(2:2:end,3);
+	index(2:2:end,2) = indexold(1:2:end,1);
+elseif orientation==3
+	%alternate diagonals
+	indexold = index;
+	index(1:4:end,3) = indexold(2:4:end,3);
+	index(2:4:end,2) = indexold(1:4:end,1);
+else
+	error('not supported');
+end
 
 %plug coordinates and nodes
@@ -102,5 +89,5 @@
 
 %Now, build the connectivity tables for this mesh.
-if computeconnectivity,
+if computeconnectivity
 	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);