sphere shell experiments

optimized_stokes_laplacian.m

@@ -38,7 +38,6 @@
     % get function handle from filename
     Aeq_fun = str2func(filename);
   end
-  Aeq_fun
   [Cb,Q,I,J,JE,JF,JT] = partial_dual_cell(V,T,'voronoi');
   nt = size(T,1);
   nf = max(JF(:));
@@ -99,11 +98,15 @@
 
   Aeq_all = [Aeq*A3;A4];
   b = [zeros(size(Aeq,1),1);ones(nf,1)];
-  %params = optimset('Display','off');
-  B = quadprog(H,f,[],[],Aeq_all,b,zeros(3*nf,1),ones(3*nf,1));
-  if isempty(B)
-    keyboard
+  params = optimset('Display','off');
+
+  if any(B(:)<-1e-12)
+    B = quadprog(H,f,[],[],Aeq_all,b,zeros(3*nf,1),ones(3*nf,1),[],params);
+  else
+    B = B(:);
+    fprintf('hooray! yur mesh is noice, skipping solve.\n');
   end
+
   Cq_f = reshape(A3*B,[],3);
   C = Cb;
   C(nv+ne+(1:nf),:) = Cq_f;

stokes_laplacian.m

@@ -21,9 +21,9 @@
   %tsurf(Q,C,'Tets',0,'CData',T(sub2ind(size(T),I,J)));
   KK = T(sub2ind(size(T),II,JJ));
 
-  [~,~,U] = unique(sort(FF,2),'rows');
-  counts = accumarray(U,1);
-  once = counts(U) == 1;
+  %[~,~,U] = unique(sort(FF,2),'rows');
+  %counts = accumarray(U,1);
+  %once = counts(U) == 1;
 
   nt = size(T,1);
   pattern = repmat([false(3*nt,1);true(3*nt,1)],4,1);
@@ -34,6 +34,7 @@
 
   % Area-weighted normals
   N = normals(C,FF);
+
   D = [];
   for d = 1:size(N,2)
     Nd = N(:,d);

test_sphere_shell.m

@@ -0,0 +1,90 @@
+% target edge length
+meths = {'dual','barycentric','alexa','optimized'};
+hs = [0.12 0.1 0.08 0.06 0.04 0.02];
+D = zeros(numel(hs),3,numel(meths));
+clf;
+hold on;
+for hi = 1:numel(hs)
+  h = hs(hi);
+
+  %%N = [100 1000 10000];
+  %%H = [];
+  %%for n = N
+  %%  [V,F] = lloyd_sphere(n);
+  %%  H = [H avgedge(V,F)];
+  %%end
+  %%%N == round(pi/2*H.^-2);
+  %
+  %VV = [];
+  %FF = [];
+  %RR = [];
+  %for R = [1.0 0.75 0.5]
+  %  [V,F] = sphere_with_target_edge_length(R,h);
+  %  FF = [FF;F+size(VV,1)];
+  %  VV = [VV;V];
+  %  RR = [RR;repmat(R,size(V,1),1)];
+  %end
+  %[TV,TT] = tetrahedralize(VV,FF,'Flags',sprintf('-Yq1a%0.17f',h^3));
+  %BC = barycenter(TV,TT);
+  %W = winding_number(V,F,BC);
+  %TT = TT(abs(W)<0.5,:);
+  %[TV,~,~,TT] = remove_unreferenced(TV,TT);
+  %I050 = find(RR==0.50);
+  %I075 = find(RR==0.75);
+  %I100 = find(RR==1.00);
+  %
+  %%writeMESH(sprintf('sphere_%0.2f.mesh',h),TV,TT,boundary_faces(TT));
+  %save(sprintf('sphere_%0.2f.mat',h),'TV','TT','I050','I075','I100','h');
+  load(sprintf('sphere_%0.2f.mat',h));
+
+  for mi = 1:numel(meths)
+    method = meths{mi};
+    switch method
+    case 'dual'
+      [L,M] = dual_laplacian(TV,TT);
+    case 'barycentric'
+      L = cotmatrix(TV,TT);
+      M = massmatrix(TV,TT,'barycentric');
+    case {'alexa','optimized'}
+      [L, ~, ~, C, Q, I, J] = stokes_laplacian(TV, TT, method);
+      M = stokes_massmatrix(TV, TT, C, Q, I, J);
+    end
+
+    Z = min_quad_with_fixed(-L,[],[I050;I100],[ones(numel(I050),1);zeros(numel(I100),1)]);
+    D(hi,:,mi) = [avgedge(TV,TT) var(Z(I075)) size(TV,1)];
+    D(hi,2,mi)
+
+    if mi == numel(meths)
+      clf;
+      hold on;
+      for ki = 1:numel(meths)
+        loglog(D(1:hi,1,ki),D(1:hi,2,ki),'LineWidth',2);
+      end
+      hold off;
+      set(gca,'YScale','log','XScale','log');
+      set(gca,'XDir','reverse');
+      legend(meths,'FontSize',20);
+      drawnow;
+    end
+  end
+
+end
+
+
+
+function [V,F] = sphere_with_target_edge_length(R,h)
+  ns = (pi*(h/R).^-2).*[3 6];
+  % binary search on avgedge length 
+  for iter = 1:7
+    n = round(mean(ns));
+    [V,F] = lloyd_sphere(n);
+    V = V*R;
+    h_iter = avgedge(V,F);
+    if h_iter > h
+      ns(1) = n;
+    else
+      ns(2) = n;
+    end
+  end
+
+end