Changeset 28136

Timestamp:

03/08/24 06:15:26 (13 months ago)

Author:

Mathieu Morlighem

Message:

CHG: allow to swap edges

File:

• issm/trunk-jpl/src/m/mesh/squaremesh.m (modified) (5 diffs)

issm/trunk-jpl/src/m/mesh/squaremesh.m

@@ -1 @@
-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 +11 @@
 
 %process options
-if nargin == 5,
+if nargin == 5
         computeconnectivity = 1;
+        orientation = 1;
+elseif nargin==6
+        orientation = 1;
 end
 
@@ -19 +22 @@
 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 +58 @@
 %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 +89 @@
 
 %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);