finish Poisson3RT0

data/mixedPossiondata.m

@@ -0,0 +1,32 @@
+function pde = mixedPossiondata
+%% MIXBCDATA3 mix boundary condition data for Poisson equation in 3-D
+%
+%     u = sin(pi*x)^2*sin(pi*y)^2*sin(pi*z)^2;
+%     f = - div Du
+%
+% Created by Yongke Wu.
+
+pde.Du = @Du;
+pde.f = @f;
+pde.u = @u;
+
+    function s = u(p)
+        x = p(:,1); y = p(:,2); z = p(:,3);
+        s = (sin(pi.*x)).^2.*(sin(pi.*y)).^2.*(sin(pi.*z)).^2;
+    end
+
+    function s = Du(p)
+        x = p(:,1); y = p(:,2); z = p(:,3);
+        s = 2*pi.*[sin(pi.*x).*cos(pi.*x).*(sin(pi.*y)).^2.*(sin(pi.*z)).^2,...
+                   (sin(pi.*x)).^2.*sin(pi.*y).*cos(pi.*y).*(sin(pi.*z)).^2,...
+                   (sin(pi.*x)).^2.*(sin(pi.*y)).^2.*sin(pi.*z).*cos(pi.*z)];
+    end
+
+    function s = f(p)
+        x = p(:,1); y = p(:,2); z = p(:,3);
+        s = - 2*pi^2*((cos(pi*x)).^2.*(sin(pi*y)).^2.*(sin(pi*z)).^2 ...
+                  + (cos(pi*y)).^2.*(sin(pi*x)).^2.*(sin(pi*z)).^2 ...
+                  + (cos(pi*z)).^2.*(sin(pi*x)).^2.*(sin(pi*y)).^2 ...
+                  - 3*(sin(pi*x)).^2.*(sin(pi*y)).^2.*(sin(pi.*z)).^2);
+    end
+end

dof/dof3face.m

@@ -22,4 +22,4 @@
                    elem(:,[1 2 4]); elem(:,[1 3 2])]); % induced ordering
 [face, ~, j] = myunique(sort(totalFace,2));               
 elem2face = uint32(reshape(j,NT,4));
-elem2faceSign = int8(reshape(sum(sign(diff(totalFace(:,[1:3,1]),1,2)),2),NT,4));      
+elem2faceSign = reshape(sum(sign(diff(totalFace(:,[1:3,1]),1,2)),2),NT,4);      

dof/sortelem.m

@@ -3,11 +3,18 @@
 %
 % [elem,bdFlag] = SORTELEM(elem,bdFlag) sorts the elem such that
 % elem(t,1)< elem(t,2)< elem(t,3). A simple sort(elem,2) cannot
-% sort bdFlag.
+% switch bdFlag.
+%
+% See also  sortelem3
 %
 % Copyright (C) Long Chen. See COPYRIGHT.txt for details.
 
-%% Step 1: make elem(:,3) to be the biggest one
+%% Dimension check
+if size(elem,2) >= 4 % 3D case
+    [elem,bdFlag] = sortelem3(elem,bdFlag);
+end
+
+%% Step 1: make elem(:,3) to the largest one
 [tempvar,idx] = max(elem,[],2); %#ok<ASGLU>
 elem(idx==1,1:3) = elem(idx==1,[2 3 1]);
 elem(idx==2,1:3) = elem(idx==2,[3 1 2]);
@@ -16,7 +23,7 @@
     bdFlag(idx==2,1:3) = bdFlag(idx==2,[3 1 2]);
 end
 
-%% Step 2: sort first two vertices such that elem(:,1)<elem(:2)
+%% Step 2: swtich the first two vertices such that elem(:,1)<elem(:2)
 idx = (elem(:,2) < elem(:,1));
 elem(idx,[1 2]) = elem(idx,[2 1]);
 if exist('bdFlag','var') && ~isempty(bdFlag)

dof/sortelem3.m

@@ -3,11 +3,13 @@
 %
 % [elem,bdFlag] = sortelem3(elem,bdFlag) sorts the elem such that
 % elem(t,1)< elem(t,2)< elem(t,3)<elem(t,4). A simple sort(elem,2) cannot
-% sort bdFlag.
+% switch bdFlag.
+% 
+% See also  sortelem
 %
 % Copyright (C) Long Chen. See COPYRIGHT.txt for details.
 
-%% Step 1: elem(:,4) is the biggest one
+%% Step 1: elem(:,4) is the largest one
 [tempvar,idx] = max(elem,[],2);  %#ok<*ASGLU>
 elem(idx==1,1:4) = elem(idx==1,[2 4 3 1]);
 elem(idx==2,1:4) = elem(idx==2,[3 4 1 2]);
@@ -17,6 +19,7 @@
     bdFlag(idx==2,1:4) = bdFlag(idx==2,[3 4 1 2]);
     bdFlag(idx==3,1:4) = bdFlag(idx==3,[4 2 1 3]);
 end
+
 %% Step 2: elem(:,1) is the smallest one
 [tempvar,idx] = min(elem(:,1:3),[],2);
 % elem(idx==1,1:3) = elem(idx==1,[1 2 3]);

equation/Poisson3RT0.m

@@ -25,7 +25,7 @@
 
 %% Data structure
 elemold = elem;
-[elem,bdFlag] = sortelem(elem,bdFlag);  % ascend ordering
+[elem,bdFlag] = sortelem3(elem,bdFlag);  % ascend ordering
 [elem2face,face] = dof3face(elem);
 [Dlambda,volume,elemSign] = gradbasis3(node,elem);
 
@@ -46,7 +46,7 @@
 % C. zero matrix.
 C = sparse(Nu,Nu);
 
-A = [M B';B C];
+A = [M B';B 1e-10*C];
 
 %% Assemble right hand side.% if ~isempty(bdEdge)
 fu = zeros(Nu,1);

example/fem/mixedPoisson/cubePoisson3RT0.m

@@ -1,25 +1,31 @@
-[node,elem] = cubemesh([-1,1,-1,1,-1,1],1); 
+
+[node,elem,~] = cubemesh([0,1,0,1,0,1],1);
 bdFlag = setboundary(node,elem,'Dirichlet'); 
-pde = mixBCdata3;
-maxIt = 3;
+% pde = mixBCdata3;
+pde = mixedPossiondata;
+maxIt = 4;
 option.solver = 'uzawapcg';
+option.fquadorder = 5;
 
 %% Solve and plot
-err = zeros(maxIt,3); N = zeros(maxIt,1);
+err = zeros(maxIt,2); 
+N = zeros(maxIt,1);
 for i =1:maxIt
-%     temperr = mixedPoisson(node,elem,pde,bdFlag,option);
+%     err(i,1:2) = mixedPoisson(node,elem,pde,bdFlag,option);
     [u,sigma,eqn] = Poisson3RT0(node,elem,bdFlag,pde,option);
 %     err(i,1) = temperr(1);
 %     err(i,2) = temperr(2);
     err(i,1) = getL2error3RT0(node,elem,pde.Du,sigma);
     sigmaI = faceinterpolate3(pde.Du,node,elem);
     err(i,2) = getL2error3RT0(node,elem,pde.Du,sigmaI);
 % %     err(i,3) = getL2error3RT0(node,elem,pde.Du,sigma,[]);
-    err(i,3) = sqrt((sigma-sigmaI)'*eqn.M*(sigma-sigmaI));
+%     err(i,3) = sqrt((sigma-sigmaI)'*eqn.M*(sigma-sigmaI));
 %     err(i,4) = getHdiverror3RT0(node,elem,pde.f,-sigma,[]);
     N(i) = size(elem,1);
     [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
 end
 figure(1); hold on; clf;
-r3 = showrate3(N,err(:,1),1,'r-','L2sigma',N,err(:,2),1,'k-','L2u',...
-               N,err(:,3),1,'b-','\sigma_I - \sigma_h');
+showrate2(N,err(:,1),1,'r-','L2sigma',N,err(:,2),1,'k-','L2u');
+% 
+% r3 = showrate3(N,err(:,1),1,'r-','L2sigma',N,err(:,2),1,'k-','L2u',...
+%                N,err(:,3),1,'b-','\sigma_I - \sigma_h');

fem/faceinterpolate3.m

@@ -2,11 +2,12 @@
 %% FACEINTERPOLATE3 interpolate to face elements RT0.
 %
 % uI = faceinterpolate3(u,node,face) interpolates a given function u
-% into the lowesr order RT0. The coefficient is given by the face integral 
+% into the lowest order RT0. The coefficient is given by the face integral 
 % int_f u*n ds. 
 %
 % uI = faceinterpolate3(u,node,elem) when |face| is not given, the function
-% will generate the face by |[elem2face,face] = dof3face(elem)|.
+% will generate the face by |[elem2face,face] = dof3face(elem)| and the
+% ascend order is used |elem = sortelem3(elem)|.
 %
 % uI = faceinterpolate3(u,node,face,6) the last input is the quadrature
 % order; see quadpts. 
@@ -39,7 +40,7 @@
     face = elem;
 else
     elem = sortelem3(elem); 
-    [elem2face,face] = dof3face(elem); 
+    [~,face] = dof3face(elem); 
 end
 
 %% normal vector of each face

fem/getL2error3RT0.m

@@ -1,10 +1,14 @@
 function err = getL2error3RT0(node,elem,sigma,sigmah,markedElem)
 %% GETL2ERROR3RT0  L2 norm of RT0 element in 3D.
 % 
-%  The input sigma can now just be function boundle. sigmah is the 
-%  flux through faces. 
 %  err = getL2error3RT0(node,elem,sigma,sigmah,markedElem)
-% 
+%
+%  The input sigma can just be a function boundle. sigmah is the flux
+%  through faces. markedElem is used to compute the error in certain region
+%  only.
+%  
+%  Note that the ascend ordering of elem is used. 
+%   
 % Example
 %   
 %   [node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
@@ -38,15 +42,17 @@
 for p = 1:nQuad
     % quadrature points in the x-y-z coordinate
     pxyz = lambda(p,1)*node(elem(:,1),:) ...
-        + lambda(p,2)*node(elem(:,2),:) ... 
-        + lambda(p,3)*node(elem(:,3),:) ...
-        + lambda(p,4)*node(elem(:,4),:);
+         + lambda(p,2)*node(elem(:,2),:) ... 
+         + lambda(p,3)*node(elem(:,3),:) ...
+         + lambda(p,4)*node(elem(:,4),:);
     sigmap = sigma(pxyz);
     sigmahp = zeros(NT,3);
     for l = 1:4 % for each basis
         i = locFace(l,1); j = locFace(l,2); k = locFace(l,3);
-        % phi_k = lambda_iRot_j - lambda_jRot_i;
-        sigmahp = sigmahp + repmat(sigmah(elem2face(:,l)),1,3)*2.*...
+        % phi_l = 2(lambda_i Dlambda_j x Dlambda_k + 
+        %           lambda_j Dlambda_k x Dlambda_i + 
+        %           lambda_k Dlambda_i x Dlambda_j)
+        sigmahp = sigmahp + repmat(2*sigmah(elem2face(:,l)),1,3).*...
                    (lambda(p,i)*mycross(Dlambda(:,:,j),Dlambda(:,:,k)) + ...
                     lambda(p,j)*mycross(Dlambda(:,:,k),Dlambda(:,:,i)) + ...    
                     lambda(p,k)*mycross(Dlambda(:,:,i),Dlambda(:,:,j)));

solver/mg.m

@@ -387,7 +387,8 @@
    end
 end
 if condest(Ai{1}) > 1e16 % Ai{1} is singular
-    coarsegridsolver = 'pcg';
+    Ai{1} = Ai{1} + 1e-12*speye(size(Ai{1}));
+%     coarsegridsolver = 'pcg';
 end
 
 %% No coarsening or coarsened nodes is small