added the ability to calculate higher order mimetic divergence

src/matlab/alt_pascal.m

@@ -0,0 +1,13 @@
+function [ A ] = alt_pascal( n )
+%   Calculates the alternating pascal triangle up to n rows.
+%   This is the change of basis matrix from the basis {i in {0,1,...,n-1}:(e^(h*D)-1)^i} to 
+%   the standard basis. {i in {0,1,...,n-1}:(e^(i*h*D))}
+    A=eye(n);
+    A(:,1)=2*rem(1:n,2) - 1; %Alternating vector [1,-1,1,-1,...]
+    for i=2:n
+        A(i,2:end)=A(i-1,1:end-1)-A(i-1,2:end);
+    end
+
+
+
+end

src/matlab/basisMatrixForGradient.m

@@ -0,0 +1,43 @@
+function [A]=basisMatrixForGradient(k)
+% Returns a k+1 by k+1 triangular matrix.
+% This is the change of basis matrix from the basis {1}U{f_0(x),f_1(x),...,f_k(x)} 
+% to {1,x,x^3,...,x^(2k+1)}.Where f_i(x)=∫_1^x (t^2-1)^idt. 
+% This is the corresponding change of basis matrix for
+% extensionOperatorGradient like altPascal was for extensionOperatorDivergence.
+
+    A=zeros(k+1);
+    A(1,1)=1;
+    currentRow=zeros([1,k+1]);
+    %This is f_i(x), 
+    %Initial row set to f_0(x) which is just x-1
+    currentRow(1)=-1;
+    currentRow(2)=1;
+    alternatingBinom=zeros([1,k+1]); %this list represents x*(x^2-1)^i, which if expanded out
+    % has coordinate vector [0,(-1)^(i),(-1)^(i-1)binom(i,1),...,1]
+    alternatingBinom(2)=-1;
+    alternatingBinom(3)=1;
+    alternatingBinomCopy=zeros([1,k+1]);
+    i=1;
+    leftBoundaryBinom=1; %(-1)^(i), the first term of our alternatingBinom
+    while 1
+        A(i+1,:)=currentRow;
+        if (i==k)
+            break
+        end
+        currentRow=(1/(2*i+1)*(alternatingBinom+currentRow)-currentRow);
+        %f_i(x)=(x(x^2-1)^i)/(2i+1)-(2i)/(2i+1)f_{i-1}(x)
+        alternatingBinomCopy(2)=leftBoundaryBinom;
+        for j=3:(i+2)
+            alternatingBinomCopy(j)=alternatingBinom(j-1)-alternatingBinom(j);
+        end
+        alternatingBinomCopy(i+3)=1;
+        for j=2:(i+3)
+            alternatingBinom(j)=alternatingBinomCopy(j);
+        end
+        leftBoundaryBinom=-leftBoundaryBinom;
+        i=i+1;
+    end
+
+
+
+end

src/matlab/calculateInteriorRow.m

@@ -0,0 +1,45 @@
+function output=calculateInteriorRow(k)
+% Returns the interior row for the kth order mimetic divergence or gradient.
+%
+% Parameters:
+%                k : Order of accuracy
+    output=zeros([1,k]); 
+    % Represents ∑_{i=0}^{k/2-1}
+    % (∏_{j=1}^{i}(0.125/j-0.25))/(2i+1))(e^(h/2*d/dx)-e^(-h/2*d/dx))^(2i+1).
+    % Our answer. 
+    r=idivide(k,int32(2));
+    output(r)=-1.0;
+    output(r+1)=1.0; %our starting value for this series is (e^(h/2*d/dx)-e^(-h/2*d/dx))
+    if (k>2)
+        nextTerm=zeros([1,k]);% NextTerm in the series.
+        %Initially set to (e^(h/2*d/dx)-e^(-h/2*d/dx))^(3)
+        nextTermCopy=zeros([1,k]); 
+        nextTerm(r-1)=-1.0;
+        nextTerm(r)=3.0;
+        nextTerm(r+1)=-3.0;
+        nextTerm(r+2)=1.0;
+        scalefactor=1; %this number is ∏_{j=1}^{i}(0.125/j-0.25).
+        i=1;
+        while 1
+            scalefactor=scalefactor*(0.125/i-0.25);
+            output=output+scalefactor/(2.0*i+1.0)*nextTerm;
+            if (i==r-1)
+                break
+            end
+            nextTermCopy(r-i-1)=-1.0; %Updates next term from (e^(h/2*d/dx)-e^(-h/2*d/dx))^(2i+1)
+             %to (e^(h/2*d/dx)-e^(-h/2*d/dx))^(2i+3)
+            nextTermCopy(r-i)=1.0*(2.0*i+3.0);
+            nextTermCopy(r+i+1)=-1.0*(2.0*i+3.0);
+            nextTermCopy(r+i+2)=1.0;
+            for j=(r-i+1):(r+i)
+                nextTermCopy(j)=nextTerm(j+1)-2*nextTerm(j)+nextTerm(j-1);
+            end
+            for j=(r-i-1):(r+i+2)
+                nextTerm(j)=nextTermCopy(j);
+            end
+            i=i+1;
+
+
+        end
+    end
+end

src/matlab/div.m

@@ -26,20 +26,22 @@
     assert(k >= 2, 'k >= 2');
     assert(mod(k, 2) == 0, 'k % 2 = 0');
     assert(m >= 2*k+1, ['m >= ' num2str(2*k+1) ' for k = ' num2str(k)]);
+    n_rows=m+2;
+    n_cols=m+1;
 
-    D = sparse(m+2, m+1);
+    D = sparse(n_rows,n_cols);
 
     switch k
         case 2
-            for i = 2:m+1
+            for i = 2:n_cols
                D(i, i-1:i) = [-1 1];
             end
 
         case 4
             A = [-11/12 17/24 3/8 -5/24 1/24];
             D(2, 1:5) = A;
             D(m+1, m-3:end) = -fliplr(A);
-            for i = 3:m
+            for i = 3:n_cols-1
                D(i, i-2:i+1) = [1/24 -9/8 9/8 -1/24];
             end
 
@@ -48,7 +50,7 @@
                     31/960  -687/640 129/128   19/192 -3/32    21/640  -3/640];
             D(2:3, 1:7) = A;
             D(m:m+1, m-5:end) = -rot90(A,2);
-            for i = 4:m-1
+            for i = 4:n_cols-2
                 D(i, i-3:i+2) = [-3/640 25/384 -75/64 75/64 -25/384 3/640];
             end
 
@@ -58,9 +60,60 @@
                    -59/17920    1175/21504 -1165/1024   1135/1024    25/3072  -251/5120   25/1024   -45/7168     5/7168];
             D(2:4, 1:9) = A;
             D(m-1:m+1, m-7:end) = -rot90(A,2);
-            for i = 5:m-2
+            for i = 5:n_cols-3
                 D(i, i-4:i+3) = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168];
             end
+        otherwise 
+            neighbors = zeros(1, k); % Bandwidth = k
+            neighbors(1) = 1/2 - k/2;
+            for i = 2 : k
+                neighbors(i) = neighbors(i-1)+1;
+            end
+
+
+            A = invvander(neighbors);
+
+
+            b = zeros(k, 1);
+            b(2) = 1;
+
+
+            coeffs = A*b;
+
+            j = 1;
+            for i = k/2+1 : n_rows - k/2
+                D(i, j:j+k-1) = coeffs;
+                j = j + 1;
+            end
+
+            p = k/2-1;
+            q = k+1;
+            A = sparse(p, q);
+            for i = 1 : p % For each row of A
+                neighbors = zeros(1, q); % k+1 points are used for the boundaries
+                neighbors(1) = 1/2 - i; % Shifting the stencil to the right
+                for j = 2 : q
+                    neighbors(j) = neighbors(j-1)+1;
+                end
+                V = invvander(neighbors);
+                b = zeros(q, 1);
+                b(2) = 1;
+                coeffs = V*b;
+                A(i, 1:q) = coeffs;
+            end
+
+            D(2:p+1, 1:q) = A;
+
+
+
+            Pp = fliplr(speye(p));
+            Pq = fliplr(speye(q));
+
+            A = -Pp*A*Pq;
+
+
+           D(n_rows-p:n_rows-1, n_cols-q+1:n_cols) = A;
     end
     D = (1/dx).*D;
+
 end

src/matlab/div1dfft.m

@@ -0,0 +1,11 @@
+function [du]=div1dfft(u,k,dx)
+    m=size(u)
+    m=m(1)
+    A=extensionOperatorDivergence(k,m-1);
+    deriv=calculateInteriorRow(k)'; %size is k-1
+
+    uExt=A*u; %size is m+k-2
+    du=-ifft(fft(deriv,m+2*k-4) .* fft(uExt,m+2*k-4))/dx;
+    du=du(k:(k+m-2))
+    du=[0;du;0]
+end

src/matlab/divergencehelper.m

@@ -0,0 +1,25 @@
+function D=divergencehelper(k,m)
+%Calculates the mimetic divergence operator of order k,
+%   barring dividing by dx. 
+%   k :order of accuracy.
+%   m :number of cells. 
+    interior=calculateInteriorRow(k); %The interior row of our matrix 
+    numNonZeros=m*k; % The number of nonZeroElements in the interior stencil matrix
+    rowList=zeros([1,numNonZeros]);
+    columnList=zeros([1,numNonZeros]);
+    valueList=zeros([1,numNonZeros]);
+    elementsInserted=1;
+    for i=2:(m+1)
+        for j=1:(k)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=i+j-2;
+            valueList(elementsInserted)=interior(j);
+            elementsInserted=elementsInserted+1;
+        end
+    end
+    interiorScheme=sparse(rowList,columnList,valueList,m+2,m+k-1);
+    % This matrix is just using the interior scheme.
+    % Since we've extended our vector field beyond the boundaries it'd normally have, we don't need to use
+    % seperate power series expansions for the stencils near the boundary. 
+    D=interiorScheme*extensionOperatorDivergence(k,m);
+end

src/matlab/extensionOperatorDivergence.m

@@ -0,0 +1,53 @@
+function E=extensionOperatorDivergence(k,m)
+% Returns a m+k-1 by m+1 one matrix. Which pads an m+1` point vector field
+% With k/2-1 additional points on each side. These points are
+% approximations of
+% f(-(k/2-1)h),f(-(k/2-2)h),....,f(-h),f(1+h),f(1+2h),...,f(1+(k/2-1)h),
+% with O(h^(k+1)) error on each of the approximations. 
+% 
+% 
+% Parameters:
+%                k : Order of accuracy
+%                m : Number of cells
+    numOfRows=m+k-1;
+    numOfColumns=m+1;
+    numNonZeros=(k-2)*(k+1)+m+1;%The number of nonzero elements in E.
+    rowList=zeros([1,numNonZeros]);%Row indices for non zero elements
+    columnList=zeros([1,numNonZeros]); %Column indices for non zero elements
+    valueList=zeros([1,numNonZeros]); %Values of non zero elements
+    elementsInserted=1; %Number of elements inserted into rowList,ColumnList
+    %and ValueList
+    r=idivide(k,int32(2));
+    R=zeros(r-1,k+1);% Calculates the matrix R_{i,j}=binom(i-k/2,j-1) 
+    R(:,1)=ones([r-1,1]);
+    currentTerm=(-1);
+    for i=2:(k+1)
+        R(r-1,i)=currentTerm;
+        currentTerm=currentTerm*(-1);
+    end
+    for i=(r-2):(-1):1
+        for j=2:(k+1)
+            R(i,j)=R(i+1,j)-R(i,j-1);
+        end
+    end
+    R=R*alt_pascal(k+1); %Converts R into the standard basis 
+    for i=1:(r-1) %
+        for j=1:(k+1)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=j;
+            valueList(elementsInserted)=R(i,j);
+            rowList(elementsInserted+1)=m+k-i;
+            columnList(elementsInserted+1)=m+2-j;
+            valueList(elementsInserted+1)=R(i,j);
+            elementsInserted=elementsInserted+2;
+        end
+    end
+    for i=r:m+r
+        rowList(elementsInserted)=i;
+        columnList(elementsInserted)=i+1-r;
+        valueList(elementsInserted)=1;
+        elementsInserted=elementsInserted+1;
+    end
+    E=sparse(rowList,columnList,valueList,numOfRows,numOfColumns);
+
+end

src/matlab/extensionOperatorGradient.m

@@ -0,0 +1,48 @@
+function E=extensionOperatorGradient(k,m)
+    % Returns a m+k by m+2 matrix, 
+    % Which turns the list [f(0),f(h/2),f(3h/2),...,f(1-h/2),f(1)]
+    % Into
+    % [f(-(k/2-1/2)h),f(-(k/2-3/2)h),...,f(-h/2),
+    % f(h/2),...f(1-h/2),f(1+h/2),...f(1+(k/2-1/2)*h)]
+    numOfRows=m+k;
+    numOfColumns=m+2;
+    r=idivide(k,int32(2));
+    numNonZeros=m+k*(k+1); 
+    % The number of nonZeroElements in the interior stencil matrix
+    rowList=zeros([1,numNonZeros]);
+    columnList=zeros([1,numNonZeros]);
+    valueList=zeros([1,numNonZeros]);
+    elementsInserted=1;
+    R=zeros([r,k+1]);
+    R(:,2)=ones([1,r]);
+    R(r,2:end)=2*rem(1:k,2) - 1;
+    for i=(r-1):(-1):1
+        for j=3:(k+1)
+            R(i,j)=R(i+1,j)-R(i,j-1);
+        end
+    end
+
+    Scalefactor=diag(-(k-1):2:-1);
+    R=Scalefactor*R;
+    R(:,1)=ones([1,r]);
+    R=R*basisMatrixForGradient(k);
+    for i=1:r
+        for j=1:(k+1)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=j;
+            valueList(elementsInserted)=R(i,j);
+            rowList(elementsInserted+1)=m+k+1-i;
+            columnList(elementsInserted+1)=m+3-j;
+            valueList(elementsInserted+1)=R(i,j);
+            elementsInserted=elementsInserted+2;
+        end
+    end
+    for j=2:(m+1)
+        rowList(elementsInserted)=r+j-1;
+        columnList(elementsInserted)=j;
+        valueList(elementsInserted)=1;
+        elementsInserted=elementsInserted+1;
+    end
+    E=sparse(rowList,columnList,valueList,numOfRows,numOfColumns);
+
+end

src/matlab/grad.m

@@ -67,6 +67,8 @@
             for i = 5:m-3
                 G(i, i-3:i+4) = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168];
             end
+        otherwise 
+            G=gradienthelper(k,m);
     end
     G = (1/dx).*G;
 end

src/matlab/gradienthelper.m

@@ -0,0 +1,19 @@
+function G=gradienthelper(k,m)
+    interior=calculateInteriorRow(k);
+    numNonZeros=(m+1)*k;
+    rowList=zeros([1,numNonZeros]);
+    columnList=zeros([1,numNonZeros]);
+    valueList=zeros([1,numNonZeros]);
+    elementsInserted=1;
+    for i=1:(m+1)
+        for j=1:(k)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=i+j-1;
+            valueList(elementsInserted)=interior(j);
+            elementsInserted=elementsInserted+1;
+        end
+    end
+    interiorScheme=sparse(rowList,columnList,valueList);
+    G=interiorScheme*extensionOperatorGradient(k,m);
+
+end