discretized_laplacian_taylor_expansion.py

import sympy 

h = sympy.Symbol("h")
f = sympy.Symbol("f")
e = sympy.Symbol("e")
c = sympy.Symbol("c")

def taylor_expansion(a, b, c):
    n = 4
    derivatives = [0] * n
    steps = [0] * n

    result = sympy.IndexedBase("p")

    for order in range(1, n+1):
        for i in range(0, 3**order):
            for j in range(0, order):
                derivatives[j] = int((i / (3**(order-1-j))) % 3 + 1)
                step = [a, b, c][derivatives[j] - 1]
                steps[j] = step

            total_step = 1
            symbol_name = ""
            for j in range(0, order):
                total_step *= (steps[j] * h)
                component = derivatives[j] - 1
                symbol_name += ["x", "y", "z"][component]
            symbol_name = ''.join(sorted(symbol_name))

            result += (1 / sympy.factorial(order)) * total_step * sympy.Symbol(symbol_name)

    return result;

def discretized_laplacian(face, edge, corner):
    # Note: any linear combination of the face, edge and corner coefficients will
    # produce a valid discretized Laplacian, except for when normalization == 0
    center = face * 6 + edge * 12 + corner * 8
    normalization = face + 4*edge + 4*corner

    expr = (
        - center * taylor_expansion(0, 0, 0)

        + face * taylor_expansion( 0, -1,  0)
        + face * taylor_expansion( 0,  0,  1)
        + face * taylor_expansion( 0,  0, -1)
        + face * taylor_expansion( 1,  0,  0)
        + face * taylor_expansion(-1,  0,  0)
        + face * taylor_expansion( 0,  1,  0)

        +  edge * taylor_expansion( 1,  1,  0)
        +  edge * taylor_expansion( 1, -1,  0)
        +  edge * taylor_expansion( 1,  0,  1)
        +  edge * taylor_expansion( 1,  0, -1)
        +  edge * taylor_expansion(-1,  1,  0)
        +  edge * taylor_expansion(-1, -1,  0)
        +  edge * taylor_expansion(-1,  0,  1)
        +  edge * taylor_expansion(-1,  0, -1)
        +  edge * taylor_expansion( 0,  1,  1)
        +  edge * taylor_expansion( 0,  1, -1)
        +  edge * taylor_expansion( 0, -1,  1)
        +  edge * taylor_expansion( 0, -1, -1)

        +  corner * taylor_expansion( 1,  1,  1)
        +  corner * taylor_expansion( 1,  1, -1)
        +  corner * taylor_expansion( 1, -1,  1)
        +  corner * taylor_expansion( 1, -1, -1)
        +  corner * taylor_expansion(-1,  1,  1)
        +  corner * taylor_expansion(-1,  1, -1)
        +  corner * taylor_expansion(-1, -1,  1)
        +  corner * taylor_expansion(-1, -1, -1)
    ) / (normalization*h*h)
    print("27-point {} {} {} {} {}:     ".format(center, face, edge, corner, normalization), expr.simplify())
    return expr

L7  = discretized_laplacian(1, 0, 0)
L13 = discretized_laplacian(0, 1, 0)
L9  = discretized_laplacian(0, 0, 1)

L = f*L7 + e*L13 + c*L9

print(L7)
print(L13)
print(L9)
print(L)