solely irregular, winding, or non-constant coefficient case handling

Author: ajozefiak

src/derivative_operators/derivative_operator_functions.jl

@@ -594,35 +594,26 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqComposi
 
         # Handle first case non-additively
         N = diff_axis(A.ops[1])
+
         if N == 1
-            if operating_dims[2] == 1
-                mul!(x_temp,A.ops[1],view(M,1:x_temp_1+2,1:x_temp_2))
-            else
-                mul!(x_temp,A.ops[1],M)
-            end
+            mul!(x_temp, A.ops[1], view(M,1:x_temp_1+2,1:x_temp_2,1:x_temp_3))
+        elseif N == 2
+             mul!(x_temp, A.ops[1], view(M,1:x_temp_1,1:x_temp_2+2,1:x_temp_3))
         else
-            if operating_dims[1] == 1
-                mul!(x_temp,A.ops[1],view(M,1:x_temp_1,1:x_temp_2+2))
-            else
-                mul!(x_temp,A.ops[1],M)
-            end
+            mul!(x_temp, A.ops[1], view(M,1:x_temp_1,1:x_temp_2,1:x_temp_3+2))
         end
 
         for L in A.ops[2:end]
             N = diff_axis(L)
+
             if N == 1
-                if operating_dims[2] == 1
-                    mul!(x_temp,L,view(M,1:x_temp_1+2,1:x_temp_2), overwrite = false)
-                else
-                    mul!(x_temp,L,M, overwrite = false)
-                end
+                mul!(x_temp, L, view(M,1:x_temp_1+2,1:x_temp_2,1:x_temp_3), overwrite = false)
+            elseif N == 2
+                 mul!(x_temp, L, view(M,1:x_temp_1,1:x_temp_2+2,1:x_temp_3), overwrite = false)
             else
-                if operating_dims[1] == 1
-                    mul!(x_temp,L,view(M,1:x_temp_1,1:x_temp_2+2), overwrite = false)
-                else
-                    mul!(x_temp,L,M, overwrite = false)
-                end
+                mul!(x_temp, L, view(M,1:x_temp_1,1:x_temp_2,1:x_temp_3+2), overwrite = false)
             end
+            
         end
     end
 end

test/2D_3D_fast_multiplication.jl

@@ -941,3 +941,67 @@ end
      + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]
 
 end
+
+@testset "x, y, and z are both irregular grids" begin
+
+    N = 100
+    dx = cumsum(rand(N+2))
+    dy = cumsum(rand(N+2))
+    dz = cumsum(rand(N+2))
+    M = zeros(N+2,N+2,N+2)
+
+    for i in 1:N+2
+        for j in 1:N+2
+            for k in 1:N+2
+                M[i,j,k] = cos(dx[i])+sin(dy[j]) + exp(dz[k])
+            end
+        end
+    end
+
+    # Lx2 has 0 boundary points
+    Lx2 = CenteredDifference{1}(2,2,dx,N)
+    # Lx3 has 1 boundary point
+    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
+    # Lx4 has 2 boundary points
+    Lx4 = CenteredDifference{1}(4,4,dx,N)
+
+    # Ly2 has 0 boundary points
+    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
+    # Ly3 has 1 boundary point
+    Ly3 = CenteredDifference{2}(3,3,dy,N)
+    # Ly4 has 2 boundary points
+    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)
+
+    # Lz2 has 0 boundary points
+    Lz2 = CenteredDifference{3}(2,2,dz,N)
+    # Lz3 has 1 boundary point
+    Lz3 = CenteredDifference{3}(3,3,dz,N)
+    # Lz4 has 2 boundary points
+    Lz4 = CenteredDifference{3}(4,4,dz,N)
+
+    # Test composition of all first-dimension operators
+    A = Lx2+Lx3+Lx4
+    M_temp = zeros(N,N+2,N+2)
+    mul!(M_temp, A, M)
+    @test M_temp ≈ (Lx2*M + Lx3*M + Lx4*M)
+
+    # Test composition of all second-dimension operators
+    A = Ly2+Ly3+Ly4
+    M_temp = zeros(N+2,N,N+2)
+    mul!(M_temp, A, M)
+    @test M_temp ≈ (Ly2*M + Ly3*M + Ly4*M)
+
+    # Test composition of all second-dimension operators
+    A = Lz2+Lz3+Lz4
+    M_temp = zeros(N+2,N+2,N)
+    mul!(M_temp, A, M)
+    @test M_temp ≈ (Lz2*M + Lz3*M + Lz4*M)
+
+    # Test composition of all operators
+    A = Lx2+Lx3+Lx4+Ly2+Ly3+Ly4+Lz2+Lz3+Lz4
+    M_temp = zeros(N,N,N)
+    mul!(M_temp, A, M)
+    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
+     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]
+
+end