Global to local indexing

DD/SchurNKP/demoSchurNKP.m

@@ -1,4 +1,34 @@
-function [lam] = demoSchurNKP(shape, ref, m, k)
+function [lam] = demoSchurNKP(m, shape, ref)
+% m = Number of gridpoint in 1D per quadrilateral element
+% shape = 1, Equilateral triangle
+% shape = 2, Right triangle
+% shape = 3, Isoceles triangle
+% shape = 4, Scalene triangle
+% shape = 5, L shaped membrane
+% shape = 6, Unit disc
+% shape = 7, Newtonian binary stars
+% ref = Level of refinement (0 = no refinement)
+
+function F=myrhs(x,y)
+    F=1+0*x+0*y;
+end
+
+R1=1;
+R2=1;
+L=6;
+function F=binrhs(z,rho)
+    z1= L;
+    z2=-L;
+    r1=hypot(rho,z-z1);
+    r2=hypot(rho,z-z2);
+    F1=(r1<=R1).*sinc(pi*r1/R1);
+    F2=(r2<=R2).*sinc(pi*r2/R2);
+    F=F1+F2;
+end
+
+RHS=@myrhs;
+metric = @(x,y) 1+0*x+0*y;
+
 % Triangle Vertices
 tri=zeros(3,4);
 tri(:,1)=[1i;exp(1i*pi*7/6);exp(-1i*pi*1/6)]; % Equilateral
@@ -22,21 +52,19 @@
         case 2
             [z0,quads,curv]=coregrid(1);
         case 3
-            [z0,quads,curv]=bingrid(1,1,3);
+            [z0,quads,curv]=bingrid(R1,R2,L);
+            RHS=@binrhs;
+            metric = @(z,rho) abs(rho)+0*z;
     end
 end
 
 % Mesh refinement
 for j=1:ref
-    [~, adj, ~, ~, bnd] = meshtopo(quads);
+    [~, adj, bnd] = meshtopo(quads);
     [z0,quads,curv]=quadmeshrefine(z0,quads,curv,adj,bnd);
 end
 
-if nargin>3
-    lam = meshSchurNKP(z0, quads, curv, m, k);
-else
-    lam = meshSchurNKP(z0, quads, curv, m);
-end
+lam = meshSchurNKP(m, z0, quads, curv, 0, RHS, metric);
 
 % delete(findall(gcf,'Type','light'));
 % shading faceted;

DD/SchurNKP/diffgeom.m

@@ -18,5 +18,5 @@
 G11=reshape(gg(1,1,:), m,n);
 G12=reshape(gg(1,2,:), m,n);
 G22=reshape(gg(2,2,:), m,n);
-jac=reshape(sqrt(max(detG{0},0)), m,n);
+jac=reshape(sqrt(detG{0}), m,n);
 end

DD/SchurNKP/feedSchurNKP.m

@@ -1,6 +1,6 @@
-function [S,nkp,gf,kd,gb] = feedSchurNKP(S,net,A1,B1,A2,B2,C1,C2)
+function [ischur,jschur,eschur,nkp,gf]=feedSchurNKP(GID,A1,B1,A2,B2,C1,C2)
 % Updates SchurNKP S given domain topology net (EWNS) and NKP (A1,B1,A2,B2)
-% Returns preconditioned Green's function gf and constrained basis gb = kd'
+% Returns preconditioned Green's function gf
 m=size(A1,1);
 n=size(B2,1);
 kd1=2:m-1; rd1=[1,m];
@@ -15,25 +15,19 @@
 
 % Green's function
 LL=1./((L1.*D1)*D2'+D1*(L2.*D2)');
-function uu=greenF(F,b1,b2,b0)
-    uu=zeros(m,n);
-    uu(rd1,rd2)=b0;
-    uu(rd1,kd2)=b1;
-    uu(kd1,rd2)=b2;
+function uu=greenF(F,uu)
     uu=E1*uu*E2';
     rhs=F-E1(:,kd1)'*(A1*uu*B1'+A2*uu*B2')*E2(:,kd2);
     uu=uu+V1(:,kd1)*(LL.*(V1(kd1,kd1)'*rhs*V2(kd2,kd2)))*V2(:,kd2)';
 end
-gf=@(F,b1,b2,b0) greenF(F,b1,b2,b0);
-kd=@(uu) reshape(E1(:,kd1)'*uu*E2(:,kd2),[],1);
-gb=@(uu) E1(:,kd1)*reshape(uu,length(kd1),length(kd2))*E2(:,kd2)';
+gf=@(F,uu) greenF(F,uu);
 
 % Schur complement
-blocks=size(S,1)/m;
-mask=@(x,y) sparse(x,y,1,blocks,blocks);
-
-T1=net(1:2);
-T2=net(3:4);
+ischur=zeros(m*m,16);
+jschur=zeros(m*m,16);
+eschur=zeros(m*m,16);
+T1=rd1(any(GID(rd1,:),2));
+T2=rd2(any(GID(:,rd2),1));
 
 % Building blocks
 MV1=E1'*A2*V1(:,kd1);
@@ -42,8 +36,7 @@
 KV2=E2'*B2*V2(:,kd2);
 
 % [E W] x [E W]
-[x,y]=meshgrid(T1(T1>0), T1(T1>0) );
-[I,J]=meshgrid(rd1(T1>0),rd1(T1>0));
+[I,J]=meshgrid(T1,T1);
 X11=(KV1(I(:),:).*KV1(J(:),:))*LL;
 X12=(KV1(I(:),:).*MV1(J(:),:))*LL;
 X21=(MV1(I(:),:).*KV1(J(:),:))*LL;
@@ -54,12 +47,15 @@
     AXX=AXX-MV2*diag(X12(i,:))*KV2';
     AXX=AXX-KV2*diag(X21(i,:))*MV2';
     AXX=AXX-KV2*diag(X22(i,:))*KV2';
-    S=S+kron(mask(x(i),y(i)), (AXX+AXX')/2);
+    [x,y]=ndgrid(GID(I(i),:), GID(J(i),:));
+
+    eschur(:,i)=AXX(:);
+    ischur(:,i)=x(:);
+    jschur(:,i)=y(:);
 end
 
 % [N S] x [N S]
-[x,y]=meshgrid(T2(T2>0) ,T2(T2>0) );
-[I,J]=meshgrid(rd2(T2>0),rd2(T2>0));
+[I,J]=meshgrid(T2,T2);
 Y11=(KV2(I(:),:).*KV2(J(:),:))*LL';
 Y12=(KV2(I(:),:).*MV2(J(:),:))*LL';
 Y21=(MV2(I(:),:).*KV2(J(:),:))*LL';
@@ -70,19 +66,32 @@
     AYY=AYY-MV1*diag(Y12(i,:))*KV1';
     AYY=AYY-KV1*diag(Y21(i,:))*MV1';
     AYY=AYY-KV1*diag(Y22(i,:))*KV1';
-    S=S+kron(mask(x(i),y(i)), (AYY+AYY')/2);
+    [x,y]=ndgrid(GID(:,I(i)), GID(:,J(i)));
+
+    eschur(:,i+4)=AYY(:);
+    ischur(:,i+4)=x(:);
+    jschur(:,i+4)=y(:);
 end
 
 % [E W] x [N S]
-[x,y]=meshgrid(T1(T1>0), T2(T2>0) );
-[I,J]=meshgrid(rd1(T1>0),rd2(T2>0));
+[I,J]=meshgrid(T1,T2);
 for i=1:numel(I)
     AXY=B1(:,J(i))*A1(I(i),:)+B2(:,J(i))*A2(I(i),:);
     AXY=AXY-MV2*((MV2(J(i),:)'*KV1(I(i),:)).*LL')*KV1';
     AXY=AXY-KV2*((MV2(J(i),:)'*MV1(I(i),:)).*LL')*KV1';
     AXY=AXY-MV2*((KV2(J(i),:)'*KV1(I(i),:)).*LL')*MV1';
     AXY=AXY-KV2*((KV2(J(i),:)'*MV1(I(i),:)).*LL')*MV1';
-    S=S+kron(mask(x(i),y(i)), AXY)+kron(mask(y(i),x(i)), AXY');
+
+    [x,y]=ndgrid(GID(I(i),:), GID(:,J(i)));
+    eschur(:,i+8)=AXY(:);
+    ischur(:,i+8)=x(:);
+    jschur(:,i+8)=y(:);
+
+    AXY=AXY';
+    [x,y]=ndgrid(GID(:,J(i)), GID(I(i),:));
+    eschur(:,i+12)=AXY(:);
+    ischur(:,i+12)=x(:);
+    jschur(:,i+12)=y(:);
 end
 end
 

DD/SchurNKP/lapGalerkin.m

@@ -15,12 +15,9 @@
 
 function vv=stiffFun(uu)
     uu=reshape(uu,[m,m]);
-    vv=zeros(size(uu));
-    vv=vv+(ED'*(C11.*(ED*uu*E'))*E);
-    vv=vv+(ED'*(C12.*(E*uu*ED'))*E);
-    vv=vv+(E'*(C12.*(ED*uu*E'))*ED);
-    vv=vv+(E'*(C22.*(E*uu*ED'))*ED);
-    vv=reshape(vv, size(uu));
+    ux=ED*uu*E';
+    uy=E*uu*ED';
+    vv=ED'*(C11.*ux+C12.*uy)*E+E'*(C12.*ux+C22.*uy)*ED;
 end
 stiff=@stiffFun;
 
@@ -33,14 +30,12 @@
 
 % This is so crucial, and I don't know why
 P11=lowrank(C11,1);
-P12=lowrank(C12,1);
+P12=C12;
 P22=lowrank(C22,1);
 
 % Block-to-row permuted operator (Van Loan, 1993)
 % Here we use some nifty algebra to optimize the computation
 % v11=(ED1'*diag(P11*diag(E2*X*E2'))*ED1);
-
-
 function b=stiff_hat(x, tflag)
     X=reshape(x,floor(sqrt(numel(x))),[]).';
     b=zeros(size(X));

DD/SchurNKP/loc2glob.m

@@ -0,0 +1,53 @@
+function [GID,nid] = loc2glob(m, corners, adj)
+% Local to global indexing
+GID=zeros(m,m,size(corners,1));
+rd=[1,m];
+se=m-2;
+
+% Corners
+b=(corners==0);
+corners=corners+se*size(adj,1);
+nid=max(corners(:));
+corners(b)=0;
+
+for c=1:size(corners,1)
+    GID([1,end],[1,end],c)=reshape(corners(c,:),[2,2]);
+end
+
+% Edges
+for e=1:size(adj,1)
+    eid=se*(e-1)+1:se*e;
+
+    % Quad 1
+    s=adj(e,1);
+    q=adj(e,2);
+    if(s>2)
+        % Vertical edge
+        ns=rd(s-2);
+        dr=1-2*(GID(1,ns,q)>GID(m,ns,q));
+        GID(2:end-1,ns,q)=eid(abs(min(dr,dr*se)):dr:abs(max(dr,dr*se)));
+    else
+        % Horizontal edge
+        ew=rd(s);
+        dr=1-2*(GID(ew,1,q)>GID(ew,m,q));
+        GID(ew,2:end-1,q)=eid(abs(min(dr,dr*se)):dr:abs(max(dr,dr*se)));
+    end
+
+    % Quad 2
+    s=adj(e,3);
+    q=adj(e,4);
+    if(s>2)
+        % Vertical edge
+        ns=rd(s-2);
+        dr=1-2*(GID(1,ns,q)>GID(m,ns,q));
+        GID(2:end-1,ns,q)=eid(abs(min(dr,dr*se)):dr:abs(max(dr,dr*se)));
+    else
+        % Horizontal edge
+        ew=rd(s);
+        dr=1-2*(GID(ew,1,q)>GID(ew,m,q));
+        GID(ew,2:end-1,q)=eid(abs(min(dr,dr*se)):dr:abs(max(dr,dr*se)));
+    end
+
+end
+
+end