From 2365c97c9968a56bce324fa36855344f4c208599 Mon Sep 17 00:00:00 2001
From: Shuhao Cao <scao.math@gmail.com>
Date: Tue, 4 Aug 2020 00:40:34 -0500
Subject: [PATCH] Update mg solver for linear Nedelec elements.

---
 equation/Maxwell1saddle.m                |  4 +---
 example/fem/Maxwell/cubeMaxwellSaddle.m  |  4 ++--
 example/fem/Maxwell/cubeMaxwellSaddle1.m |  5 ++---
 fem/getmassmat3.m                        |  4 ++--
 solver/mgMaxwellsaddle.m                 | 14 ++++++++------
 5 files changed, 15 insertions(+), 16 deletions(-)

diff --git a/equation/Maxwell1saddle.m b/equation/Maxwell1saddle.m
index be9864d..eb3717d 100644
--- a/equation/Maxwell1saddle.m
+++ b/equation/Maxwell1saddle.m
@@ -337,9 +337,7 @@
 end
 
 % multigrid options 
-% option.mg.isFreeEdge = isFreeEdge; % needed in mg solver
-option.mg.isFreeEdge = isFreeEdge;
-% option.mg.isFreeNode = isFreeNode; 
+option.mg.isFreeEdge = isFreeEdge; % needed in mg solver
 option.mg.isFreeNode = isFreeNode; 
 switch method
     case 'DIRECT'
diff --git a/example/fem/Maxwell/cubeMaxwellSaddle.m b/example/fem/Maxwell/cubeMaxwellSaddle.m
index a016583..ae7e50f 100644
--- a/example/fem/Maxwell/cubeMaxwellSaddle.m
+++ b/example/fem/Maxwell/cubeMaxwellSaddle.m
@@ -22,7 +22,7 @@
 %%
 bdFlag = setboundary3(node,elem,'Dirichlet');
 option.printlevel = 0;
-option.solver = 'direct';
+option.solver = 'mg';
 
 %% Parameters
 maxIt = 3; 
@@ -37,7 +37,7 @@
     [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);    
     [soln,eqn,info] = Maxwellsaddle(node,elem,bdFlag,pde,option);    
     u = soln.u;
-    fprintf('\n # of DoFs = %d \n',length(u));
+    fprintf('\n\n # of DoFs = %d \n',length(u));
     % compute error
     uI = edgeinterpolate(pde.exactu,node,eqn.edge);
     energyErr(k) = getHcurlerror3ND(node,elem,pde.curlu,u);
diff --git a/example/fem/Maxwell/cubeMaxwellSaddle1.m b/example/fem/Maxwell/cubeMaxwellSaddle1.m
index 5419d11..140ed50 100644
--- a/example/fem/Maxwell/cubeMaxwellSaddle1.m
+++ b/example/fem/Maxwell/cubeMaxwellSaddle1.m
@@ -7,7 +7,6 @@
 
 %% Defacult setting
 [node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
-[node,elem] = uniformrefine3(node,elem);
 %%
 pde.J = @(p) [sin(p(:,1)).*cos(p(:,2)).*sin(p(:,3)), ...
               cos(p(:,1)).*sin(p(:,2)).*sin(p(:,3)), ...
@@ -39,14 +38,14 @@
     [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);     
     [soln,eqn,info] = Maxwell1saddle(node,elem,bdFlag,pde,option);  
     u = soln.u;
-    fprintf('\n # of DoFs = %d \n',length(u));
+    fprintf('\n\n # of DoFs = %d \n',length(u));
     % compute error
     uI = edgeinterpolate1(pde.exactu,node,eqn.edge);
     energyErr(k) = getHcurlerror3ND1(node,elem,pde.curlu,u);
     L2Err(k) = getL2error3ND1(node,elem,pde.exactu,u);
     uIuhErr(k) = sqrt((u-uI)'*(eqn.A)*(u-uI));        
 %     L2Err(k) = sqrt((u-uI)'*(eqn.M)*(u-uI));        
-    fprintf('||curl(u-u_h)|| is %g \n',energyErr(k))
+    fprintf('\n ||curl(u-u_h)|| is %g \n',energyErr(k))
     N(k) = length(u);
     h(k) = 1./(size(node,1)^(1/3)-1);   
 end
diff --git a/fem/getmassmat3.m b/fem/getmassmat3.m
index cd956da..4253edb 100644
--- a/fem/getmassmat3.m
+++ b/fem/getmassmat3.m
@@ -5,8 +5,8 @@
 % space specified by elemType.  
 %
 % The type can be: 
-% - P1: full mass matrix for P1 element
-% - lump: lumped mass matrix for P1 element
+% - 'P1': full mass matrix for P1 element
+% - 'lump': lumped mass matrix for P1 element
 
 N = size(node,1);
 NT = length(volume);
diff --git a/solver/mgMaxwellsaddle.m b/solver/mgMaxwellsaddle.m
index e395c2e..ecc81d6 100644
--- a/solver/mgMaxwellsaddle.m
+++ b/solver/mgMaxwellsaddle.m
@@ -61,14 +61,15 @@
     else % linear ND (quadratic Lagrange for p)
         elem2dof = dof3P2(elem);
         M = getmassmat3(node,elem2dof,volume);
-        Mvlump = sum(M,2);
+%         Mvlump = sum(M,2);
+        Mvlump = diag(M);
     end
 end
 if N>= Ng
     DMinv = spdiags(1./Mvlump(option.isFreeNode),0,Ng,Ng);
 else
-    isFreeDofu = [option.isFreeNode; option.isFreeEdge];
-    DMinv = spdiags(1./Mvlump(isFreeDofu),0,Ng,Ng);
+    isFreeDofp = [option.isFreeNode; option.isFreeEdge];
+    DMinv = spdiags(1./Mvlump(isFreeDofp),0,Ng,Ng);
 end
 f = f + G*(DMinv*g);  % add second equation to the first one
 Abar = A + G*DMinv*G'; % Hodge Laplacian
@@ -85,6 +86,7 @@
 else
     HB = [];
 end
+
 if N>= Ng
     [u,info] = mgHodgeLapE(Abar,f,node,elem,bdFlag,option); % lowest order
 else
@@ -93,14 +95,14 @@
 
 Apoption.x0 = p0;
 Apoption.printlevel = 1;
-Apoption.freeDof = option.isFreeNode;
 rg = g-G'*u;
 
 if N>=Ng
+    Apoption.freeDof = option.isFreeNode;
     [p,info2] = mg(Ap,rg,elem,Apoption,HB);
 else
-    [~,edge] = dof3edge(elem);
-    [p,info2] = mg(Ap,rg,elem,Apoption,HB,edge);
+    Apoption.freeDof = [option.isFreeNode; option.isFreeEdge];
+    [p,info2] = mg(Ap,rg,elem,Apoption,HB);
 end
 % [p,info2] = amg(Ap,rg,Apmgoption);