Changeset 22502

Timestamp:

03/06/18 14:46:56 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: sped up squaremesh (needed for BedMachine Ant.)

File:

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

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

@@ -1 @@
-function md=squaremesh(md,Lx,Ly,nx,ny)
+function md=squaremesh(md,Lx,Ly,nx,ny,computeconnectivity)
 %SQUAREMESH - create a structured square mesh 
 %
@@ -10 +10 @@
 %      [md]=squaremesh(md,Lx,Ly,nx,ny)
 
+%process options
+if nargin == 5,
+        computeconnectivity = 1;
+end
+
 %get number of elements and number of nodes
 nel=(nx-1)*(ny-1)*2;
 nods=nx*ny;
 
-%initialization
-index=zeros(nel,3);
-x=zeros(nx*ny,1);
-y=zeros(nx*ny,1);
+%Old method
+if 0,
+        %initialization
+        index=zeros(nel,3);
+        x=zeros(nx*ny,1);
+        y=zeros(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.);
+        %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
+
+        %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
 
-%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;
 
 %create segments
@@ -67 +102 @@
 
 %Now, build the connectivity tables for this mesh.
-md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
-md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+if computeconnectivity,
+        md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+        md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+end