From bcee34e596f593df0d69569f65719dc91ef97a73 Mon Sep 17 00:00:00 2001
From: Jan Mandel <jan.mandel@gmail.com>
Date: Sat, 12 Dec 2020 15:27:55 -0700
Subject: [PATCH] progress by 5% in sparse_assemblt

---
 femwind/femwind_test.m    | 87 +++++++++++++++++++++++++----------------------
 femwind/sparse_assembly.m | 15 +++++---
 2 files changed, 58 insertions(+), 44 deletions(-)

diff --git a/femwind/femwind_test.m b/femwind/femwind_test.m
index f051d2b..0fbb3a0 100644
--- a/femwind/femwind_test.m
+++ b/femwind/femwind_test.m
@@ -1,9 +1,11 @@
 disp('femwind_test')
 
+do_plot=0
+
 % dimensions in elements
 sc=1; % mesh scale
 n = sc*[20,10,5];
-sc2=2;
+sc2=10;
 n(1:2)=n(1:2)*sc2
 h = [1,1,1]/sc;
 fprintf('linear array of %ix%ix%i cells\n',n(1),n(2),n(3))
@@ -17,30 +19,33 @@ X = regular_mesh(n,h,1.5^(1/sc));
 X = add_terrain_to_mesh(X,'hill','squash',0.3);
 CX = center_mesh(X);
 
-% show mesh
-hold off
-figure
-plot_mesh_3d(X), hold on, 
-title('The wind mesh, wind vector in centers, lambda in corners')
-
 % initial wind at the centers of the elements
 U0={ones(n),zeros(n),zeros(n)};
 
-% show initial wind
-figure 
-plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
-plot_wind_3d(CX,U0)
-hold off
-axis equal
-title('Initial wind')
+if do_plot
+
+    % show mesh
+    hold off
+    figure
+    plot_mesh_3d(X), hold on, 
+    title('The wind mesh, wind vector in centers, lambda in corners')
+
+    % show initial wind
+    figure 
+    plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
+    plot_wind_3d(CX,U0)
+    hold off
+    axis equal
+    title('Initial wind')
 
-% show initial wind
-figure 
-plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
-plot_wind_3d(CX,U0,1)
-hold off
-axis equal
-title('Initial wind lowest layer')
+    % show initial wind
+    figure 
+    plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
+    plot_wind_3d(CX,U0,1)
+    hold off
+    axis equal
+    title('Initial wind lowest layer')
+end
 
 % assemble sparse system matrix
 [K,F,~] = sparse_assembly(A,X,U0,lambda);
@@ -54,26 +59,28 @@ lambda = sparse_solve(K,F,X,'s');
 % assemble final wind
 [~,~,W] = sparse_assembly(A,X,U0,lambda);
 
-% plot resulting wind
-figure
-plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
-plot_wind_3d(CX,W)
-hold off
-axis equal
-title(['Final wind a=',string_diag_A])
+if do_plot
+    % plot resulting wind
+    figure
+    plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
+    plot_wind_3d(CX,W)
+    hold off
+    axis equal
+    title(['Final wind a=',string_diag_A])
 
-figure
-plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
-plot_wind_3d(CX,W,1)
-hold off
-axis equal
-title(['Final wind lowest layer a=',string_diag_A])
+    figure
+    plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
+    plot_wind_3d(CX,W,1)
+    hold off
+    axis equal
+    title(['Final wind lowest layer a=',string_diag_A])
 
-figure
-plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
-plot_wind_3d(CX,W,1:2)
-hold off
-axis equal
-title(['Final wind lowest layers a=',string_diag_A])
+    figure
+    plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on, 
+    plot_wind_3d(CX,W,1:2)
+    hold off
+    axis equal
+    title(['Final wind lowest layers a=',string_diag_A])
+end
 % condition_number=scond(K)
 n
diff --git a/femwind/sparse_assembly.m b/femwind/sparse_assembly.m
index dbbad02..97e7c3b 100644
--- a/femwind/sparse_assembly.m
+++ b/femwind/sparse_assembly.m
@@ -28,10 +28,11 @@ W = {zeros(n-1),zeros(n-1),zeros(n-1)};
 tstart=tic;
 Xloc = zeros(3,8);
 kglo=zeros(1,8);
-for i3=1:n(3)-1
-    fprintf(' %g %g %g',100*i3/(n(3)-1))
-    for i2=1:n(2)-1
-        for i1=1:n(1)-1  % loop over elements
+m = n-1;
+done_last=0;
+for i3=1:m(3)
+    for i2=1:m(2)
+        for i1=1:m(1)  % loop over elements
             % now in element (i1,i2,i3)
             for j3=0:1 % loop over corners of the element
                 for j2=0:1
@@ -55,6 +56,12 @@ for i3=1:n(3)-1
                 W{i}(i1,i2,i3)=grad(i);
             end                         
         end
+        done = 100*((i3-1)*m(2)+i2)/(m(3)*m(2));
+        done = round(done);
+        if done>done_last+5
+            fprintf(' %g%% ',done)
+            done_last=done;
+        end
     end
 end
 for i=1:3
-- 
2.11.4.GIT