# Code generation for barycentric coordinates

This notebook consists of three main parts, all used in `ocellaris/utils/taylor_basis.py` for Lagrange to Taylor and Taylor to Lagrange transformations for discontinuous Galerkin shape functions on simplices.

**1: Code generation for quadratic barycentric shape functions**

Given a triangle/tetrahedron with barycentric coordinates, find the linear and quadratic shape functions. We here work only in barycentric coordinates. The linear basis functions should be the coordinates themselves and the interesting part is the quadratic shape functions

**2: Code generation for the cartesian to barycentric transformation**

Given cartesian coordinates of all cell vertices, compute the barycentric coordinate transformation for an arbitrary point inside the cell. This is done for 2D (fast) and 3D (dog slow)

**3: Code generation for the transformation matrices**

Produce the actual transformation matrices for the Lagrange to Taylor transform and some required parts of the matrices for the Taylor to Lagrange transform


```python
import time
from IPython.display import display, HTML
import numpy
from matplotlib import pyplot
import sympy
from sympy import *
sympy.init_printing(use_latex=True)
%matplotlib inline
```

## Barycentric basis functions in 2D


```python
barycentric_basis = {}

for P in (1, 2):
    display(HTML('<b>Barycentric basis in 2D of order %d</b>' % P))
    D = 2
    N = D + 1
    M = N*P

    # Define the barycentric coordinates L
    Ls = list(symbols('L_(1:%d)' % N))
    Ls.append(1 - sum(Ls))

    # Define the coeffcients of the basis functions, C
    Cs = symbols('C_(1:%d)' % (1 + M*M))
    
    # DOF locations
    w = 1/S(2)
    locs = [(1, 0, 0), (0, 1, 0), (0, 0, 1), (0, w, w), (w, 0, w), (w, w, 0)][:M]

    # Create the linear basis
    Bs = []
    for i in range(M):
        cs = Cs[i*M:(i + 1)*M]
        if P == 1:
            Bs.append(cs[0] + cs[1]*Ls[0] + cs[2]*Ls[1])
        else:
            Bs.append(cs[0] + cs[1]*Ls[0] + cs[2]*Ls[1] + cs[3]*Ls[0]**2 + cs[4]*Ls[1]**2 + cs[5]*Ls[0]*Ls[1])

    # Conditions on the basis, Bi(xj) = delta_j(xj)
    eqs = []
    for iloc, loc in enumerate(locs):
        subs = {Ls[i]: loc[i] for i in range(N)}
        for jb in range(M):
            rhs = 1 if iloc == jb else 0
            eqs.append((Bs[jb] - rhs).subs(subs))

    # Find the coefficients of the basis
    sol = solve(eqs, Cs)
    Ln = Symbol('L_%d' % N)
    Lnsubs = {Ls[-1]: Ln}
    Bsol = []
    for i in range(M):
        Bi = Bs[i]
        Bi = factor(Bi).subs(Lnsubs).simplify().subs(Lnsubs)
        Bi = Bi.subs(sol)
        Bi = factor(Bi).subs(Lnsubs).simplify().subs(Lnsubs)
        Bi = Bi.subs(2*Ls[-1] - 1, 2*Ln -1)
        Bsol.append(Bi.simplify())
        display(Eq(Symbol('B_%d' % i), Bsol[-1]))
    
    barycentric_basis[(D, P)] = Bsol, Ls[:-1] + [Ln]

```

## Barycentric basis functions in 3D


```python
for P in (1, 2):
    display(HTML('<b>Barycentric basis in 3D of order %d</b>' % P))
    D = 3
    N = D + 1
    M = 4 if P == 1 else 10

    # Define the barycentric coordinates L
    Ls = list(symbols('L_(1:%d)' % N))
    Ls.append(1 - sum(Ls))

    # Define the coeffcients of the basis functions, C
    Cs = symbols('C_(1:%d)' % (1 + M*M))
    
    # DOF locations
    w = 1/S(2)
    locs = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1), (0, 0, w, w),
            (0, w, 0, w), (0, w, w, 0), (w, 0, 0, w), (w, 0, w, 0), (w, w, 0, 0)][:M]

    # Define the basis in terms of the unknown coefficients, C_i
    Bs = []
    for i in range(M):
        cs = Cs[i*M:(i + 1)*M]
        if P == 1:
            Bs.append(cs[0] + cs[1]*Ls[0] + cs[2]*Ls[1] + cs[3]*Ls[2])
        else:
            Bs.append(cs[0] + cs[1]*Ls[0] + cs[2]*Ls[1] + cs[3]*Ls[2] +
                      cs[4]*Ls[0]**2 + cs[5]*Ls[1]**2 + cs[6]*Ls[2]**2 +
                      cs[7]*Ls[0]*Ls[1] + cs[8]*Ls[0]*Ls[2] + cs[9]*Ls[1]*Ls[2])

    # Conditions on the basis, Bi(xj) = delta_j(xj)
    eqs = []
    for iloc, loc in enumerate(locs):
        subs = {Ls[i]: loc[i] for i in range(N)}
        for jb in range(M):
            rhs = 1 if iloc == jb else 0
            eqs.append((Bs[jb] - rhs).subs(subs))

    # Find the coefficients of the basis
    sol = solve(eqs, Cs)
    Ln = Symbol('L_%d' % N)
    Lnsubs = {Ls[-1]: Ln}
    Bsol = []
    for i in range(M):
        Bi = Bs[i]
        Bi = factor(Bi).subs(Lnsubs).simplify().subs(Lnsubs)
        Bi = Bi.subs(sol)
        Bi = factor(Bi).subs(Lnsubs).simplify().subs(Lnsubs)
        Bi = Bi.subs(2*Ls[-1] - 1, 2*Ln -1)
        Bsol.append(Bi.simplify())
        display(Eq(Symbol('B_%d' % i), Bsol[-1]))
    
    barycentric_basis[(D, P)] = Bsol, Ls[:-1] + [Ln]

```

## Code generation


```python
def tocode(Bi):
    s = str(Bi)
    s = s.replace('L_1', 'L1')
    s = s.replace('L_2', 'L2')
    s = s.replace('L_3', 'L3')
    s = s.replace('L_4', 'L4')
    return s

for k, (Bs, Ls) in barycentric_basis.items():
    display(HTML('<b>Code for %dD, order %d' % k))
    s = '\n'.join('phi[%d] = %s' % (i, tocode(Bi)) for i, Bi in enumerate(Bs))
    display(HTML('<pre>' + s + '</pre>'))
```


<b>Code for 2D, order 1



<pre>phi[0] = L1
phi[1] = L2
phi[2] = L3</pre>



<b>Code for 2D, order 2



<pre>phi[0] = L1*(2*L1 - 1)
phi[1] = L2*(2*L2 - 1)
phi[2] = L3*(2*L3 - 1)
phi[3] = 4*L2*L3
phi[4] = 4*L1*L3
phi[5] = 4*L1*L2</pre>



<b>Code for 3D, order 1



<pre>phi[0] = L1
phi[1] = L2
phi[2] = L3
phi[3] = L4</pre>



<b>Code for 3D, order 2



<pre>phi[0] = L1*(2*L1 - 1)
phi[1] = L2*(2*L2 - 1)
phi[2] = L3*(2*L3 - 1)
phi[3] = L4*(2*L4 - 1)
phi[4] = 4*L3*L4
phi[5] = 4*L2*L4
phi[6] = 4*L2*L3
phi[7] = 4*L1*L4
phi[8] = 4*L1*L3
phi[9] = 4*L1*L2</pre>


## Compute barycentric coordinates from cartesian vertex coordinates


```python
barycentric_transform = {}

xs = x1, x2, x3, x4 = symbols('x(1:5)', real=True)
ys = y1, y2, y3, y4 = symbols('y(1:5)', real=True)
zs = z1, z2, z3, z4 = symbols('z(1:5)', real=True)
x, y, z = symbols('x y z', real=True)

for D in (2, 3):
    display(HTML('<h3>Cartesian to barycentric %dD</h3>' % D))
    t1 = time.time()
    N = D + 1
    
    # Coefficients
    Cs = symbols('C_(:%d)' % (N*D))

    # Barycentric coordinates
    if D == 2:
        Ls = [Cs[0] + Cs[1]*x + Cs[2]*y,
              Cs[3] + Cs[4]*x + Cs[5]*y]
        Ls.append(1 - Ls[0] - Ls[1])
    else:
        Ls = [Cs[0] + Cs[1]*x + Cs[2]*y + Cs[3]*z,
              Cs[4] + Cs[5]*x + Cs[6]*y + Cs[7]*z,
              Cs[8] + Cs[9]*x + Cs[10]*y + Cs[11]*z,]
        Ls.append(1 - Ls[0] - Ls[1] - Ls[2])

    # Build equation system
    eqs = []
    # Loop over corners, i
    for i in range(N):
        vertex_pos = {x: xs[i], y: ys[i], z: zs[i]}
        # Loop over barycentric coordinates, j
        for j in range(N-1):
            rhs = 1 if i == j else 0
            eqs.append(Ls[j].subs(vertex_pos) - rhs)

    display(Eq(Symbol('N_\\text{eq}'), len(eqs)))
    display(Eq(Symbol('N_\\text{Cs}'), len(Cs)))

    # Solve the equation system
    display(HTML('<b style="color:#999900">This is going to take some time, especially in 3D ...</b>'))
    sol = solve(eqs, Cs)
    
    F = Symbol('F')
    if D == 2:
        fac = x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2
        #diff = fac - (y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3)
        #display(expand(diff).simplify())
    else:
        fac = x1*y2*z3 - x1*y2*z4 - x1*y3*z2 + x1*y3*z4 + x1*y4*z2 - x1*y4*z3 - x2*y1*z3 + x2*y1*z4 + x2*y3*z1 - x2*y3*z4 - x2*y4*z1 + x2*y4*z3 + x3*y1*z2 - x3*y1*z4 - x3*y2*z1 + x3*y2*z4 + x3*y4*z1 - x3*y4*z2 - x4*y1*z2 + x4*y1*z3 + x4*y2*z1 - x4*y2*z3 - x4*y3*z1 + x4*y3*z2
    display(Eq(F, factor(fac)))
    
    # Show the solution
    Lsol = []
    for i in range(N):
        Li = (Ls[i] * fac).subs(sol).simplify()
        Lsol.append(Li/F)
        display(Eq(Symbol('L_%d' % i), factor(Lsol[-1])))
    
    barycentric_transform[D] = Lsol
    display(HTML('<b style="color:#009900">Done in %.2f seconds</b>' % (time.time() - t1)))

```

## Code generation


```python
def tocode(Li):
    return str(Li)

for D, Ls in barycentric_transform.items():
    display(HTML('<b>Code for %dD' % D))
    s = '\n'.join('L%df = %s' % (i + 1, tocode(Li)) for i, Li in enumerate(Ls))
    display(HTML('<pre>' + s + '</pre>'))

```


<b>Code for 2D



<pre>L1f = (x*(y2 - y3) + x2*y3 - x3*y2 - y*(x2 - x3))/F
L2f = (-x*(y1 - y3) - x1*y3 + x3*y1 + y*(x1 - x3))/F
L3f = (x*y1 - x*y2 - x1*y + x1*y2 + x2*y - x2*y1)/F</pre>



<b>Code for 3D



<pre>L1f = (x*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) - x2*y3*z4 + x2*y4*z3 + x3*y2*z4 - x3*y4*z2 - x4*y2*z3 + x4*y3*z2 - y*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3) + z*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3))/F
L2f = (-x*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + x1*y3*z4 - x1*y4*z3 - x3*y1*z4 + x3*y4*z1 + x4*y1*z3 - x4*y3*z1 + y*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3) - z*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3))/F
L3f = (x*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2) - x1*y2*z4 + x1*y4*z2 + x2*y1*z4 - x2*y4*z1 - x4*y1*z2 + x4*y2*z1 - y*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2) + z*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2))/F
L4f = (-x*y1*z2 + x*y1*z3 + x*y2*z1 - x*y2*z3 - x*y3*z1 + x*y3*z2 + x1*y*z2 - x1*y*z3 - x1*y2*z + x1*y2*z3 + x1*y3*z - x1*y3*z2 - x2*y*z1 + x2*y*z3 + x2*y1*z - x2*y1*z3 - x2*y3*z + x2*y3*z1 + x3*y*z1 - x3*y*z2 - x3*y1*z + x3*y1*z2 + x3*y2*z - x3*y2*z1)/F</pre>


### Test implementation against known answers


```python
# First we check for linearity
cell = {x1: 1, y1: 1, z1: 1,
        x2: 2, y2: 1, z2: 2,
        x3: 4, y3: 2, z3: 3,
        x4: 7, y4: 2, z4: 4}
Fs = {2: x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2,
      3: x1*y2*z3 - x1*y2*z4 - x1*y3*z2 + x1*y3*z4 + x1*y4*z2 - x1*y4*z3 - x2*y1*z3 + x2*y1*z4 + x2*y3*z1 - x2*y3*z4 - x2*y4*z1 + x2*y4*z3 + x3*y1*z2 - x3*y1*z4 - x3*y2*z1 + x3*y2*z4 + x3*y4*z1 - x3*y4*z2 - x4*y1*z2 + x4*y1*z3 + x4*y2*z1 - x4*y2*z3 - x4*y3*z1 + x4*y3*z2}
for D, Ls in barycentric_transform.items():
    display(HTML('<b>Code for %dD' % D))
    display(Eq(Symbol('F'), Fs[D].subs(cell)))
    for i, L in enumerate(Ls):
        LF = (L*F).subs(cell).subs(F, Fs[D].subs(cell)).simplify()
        display(Eq(Symbol('L_{%d}\\,F' % i), LF))
```


```python
def test_cell_2D(vertices):
    print('Checking cell', vertices, '...', end=' ')
    x1, y1 = vertices[0]
    x2, y2 = vertices[1]
    x3, y3 = vertices[2]
    
    def our(x, y):
        "Our solution"
        F = x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2
        L1f = (x*(y2 - y3) + x2*y3 - x3*y2 - y*(x2 - x3))/F
        L2f = (-x*(y1 - y3) - x1*y3 + x3*y1 + y*(x1 - x3))/F
        L3f = (x*y1 - x*y2 - x1*y + x1*y2 + x2*y - x2*y1)/F
        return L1f, L2f, L3f
    
    def ref(x, y):
        D = (y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3)
        L1f = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3))/D
        L2f = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3))/D
        L3f = 1 - L1f - L2f
        return L1f, L2f, L3f
    
    # Generate some positions to check
    other_positions = []
    for A, B in [(1/3, 1/3), (1/2, 1/4), (1/4, 1/2),
                 (3/5, 1/5), (1/5, 3/5), (1/5, 1/5)]:
        C = 1 - A - B
        other_positions.append((A*x1 + B*x2 + C*x3,
                                A*y1 + B*y2 + C*y3))
    
    OK = True
    for v in vertices + other_positions:
        #print('Checking', v)
        o1, o2, o3 = our(v[0], v[1])
        r1, r2, r3 = ref(v[0], v[1])
        for o, r in zip((o1, o2, o3), (r1, r2, r3)):
            if not abs(o - r) < 1e-15:
                print(o, '-', r, '=', o-r, '   <-- ERROR!!!')
                OK = False
    
    print('OK' if OK else 'ERROR')

print('Testing 2D cells')
test_cell_2D([(0, 0), (1, 0), (0, 1)])
test_cell_2D([(1, 1), (1, 0), (0, 1)])
test_cell_2D([(cell[x1], cell[y1]),
              (cell[x2], cell[y2]),
              (cell[x3], cell[y3])])
test_cell_2D([(-3, 4), (8, -1), (-15, -1)])

def test_cell_3D(vertices):
    print('Checking cell', vertices, '...', end=' ')
    x1, y1, z1 = vertices[0]
    x2, y2, z2 = vertices[1]
    x3, y3, z3 = vertices[2]
    x4, y4, z4 = vertices[3]
    
    def our(x, y, z):
        "Our solution"
        F = x1*y2*z3 - x1*y2*z4 - x1*y3*z2 + x1*y3*z4 + x1*y4*z2 - x1*y4*z3 - x2*y1*z3 + x2*y1*z4 + x2*y3*z1 - x2*y3*z4 - x2*y4*z1 + x2*y4*z3 + x3*y1*z2 - x3*y1*z4 - x3*y2*z1 + x3*y2*z4 + x3*y4*z1 - x3*y4*z2 - x4*y1*z2 + x4*y1*z3 + x4*y2*z1 - x4*y2*z3 - x4*y3*z1 + x4*y3*z2
        L1f = (x*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) - x2*y3*z4 + x2*y4*z3 + x3*y2*z4 - x3*y4*z2 - x4*y2*z3 + x4*y3*z2 - y*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3) + z*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3))/F
        L2f = (-x*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + x1*y3*z4 - x1*y4*z3 - x3*y1*z4 + x3*y4*z1 + x4*y1*z3 - x4*y3*z1 + y*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3) - z*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3))/F
        L3f = (x*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2) - x1*y2*z4 + x1*y4*z2 + x2*y1*z4 - x2*y4*z1 - x4*y1*z2 + x4*y2*z1 - y*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2) + z*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2))/F
        L4f = (-x*y1*z2 + x*y1*z3 + x*y2*z1 - x*y2*z3 - x*y3*z1 + x*y3*z2 + x1*y*z2 - x1*y*z3 - x1*y2*z + x1*y2*z3 + x1*y3*z - x1*y3*z2 - x2*y*z1 + x2*y*z3 + x2*y1*z - x2*y1*z3 - x2*y3*z + x2*y3*z1 + x3*y*z1 - x3*y*z2 - x3*y1*z + x3*y1*z2 + x3*y2*z - x3*y2*z1)/F
        return L1f, L2f, L3f, L4f
    
    def ref(x, y, z):
        if (x, y, z) == vertices[0]:
            return 1, 0, 0, 0
        elif (x, y, z) == vertices[1]:
            return 0, 1, 0, 0
        elif (x, y, z) == vertices[2]:
            return 0, 0, 1, 0
        elif (x, y, z) == vertices[3]:
            return 0, 0, 0, 1
        elif (x, y, z) == other_positions[0]: # centroid
            return 1/4, 1/4, 1/4, 1/4
    
    ## Generate some positions to check
    other_positions = []
    for A, B, C in [(1/4, 1/4, 1/4)]:
        D = 1 - A - B - C
        other_positions.append((A*x1 + B*x2 + C*x3 + D*x4,
                                A*y1 + B*y2 + C*y3 + D*y4,
                                A*z1 + B*z2 + C*z3 + D*z4))
    
    OK = True
    for v in vertices + other_positions:
        #print('Checking', v)
        os = our(v[0], v[1], v[2])
        rs = ref(v[0], v[1], v[2])
        for i in range(4):
            o, r = os[i], rs[i]
            #print('    %d: % .2f - % .2f = % .2f' % (i, o, r, o-r),
            #      '   <-- ERROR!!!' if not abs(o - r) < 1e-15 else '')
            if not abs(o - r) < 1e-15:
                OK = False
    
    print('OK' if OK else 'ERROR')

print('Testing 3D cells')
test_cell_3D([(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)])
test_cell_3D([(0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 1)])
test_cell_3D([(cell[x1], cell[y1], cell[z1]),
              (cell[x2], cell[y2], cell[z2]),
              (cell[x3], cell[y3], cell[z3]),
              (cell[x4], cell[y4], cell[z4])])
test_cell_3D([(-4, 2, 33), (-2, -1, -50), (-10, 1, -30), (55, 22, -30)])
```

    Testing 2D cells
    Checking cell [(0, 0), (1, 0), (0, 1)] ... OK
    Checking cell [(1, 1), (1, 0), (0, 1)] ... OK
    Checking cell [(1, 1), (2, 1), (4, 2)] ... OK
    Checking cell [(-3, 4), (8, -1), (-15, -1)] ... OK
    Testing 3D cells
    Checking cell [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)] ... OK
    Checking cell [(0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 1)] ... OK
    Checking cell [(1, 1, 1), (2, 1, 2), (4, 2, 3), (7, 2, 4)] ... OK
    Checking cell [(-4, 2, 33), (-2, -1, -50), (-10, 1, -30), (55, 22, -30)] ... OK


## Code generation for transformation matrices, 2D

Use sympy to derive the coefficients needed to convert from DG polynomials
on a triangle to a Taylor basis by expressing the derivatives at the centre
point of the triangle


```python
%reset -f

from sympy import Symbol, symbols, S
x1, x2, x3, y1, y2, y3 = symbols('x1 x2 x3 y1 y2 y3')
x, y, D = symbols('x y D')

# Barycentric coordinates from cartesian corner coordinates
L1f = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3))/D
L2f = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3))/D
L3f = 1 - L1f - L2f
L1, L2, L3 = Symbol('L1')(x, y), Symbol('L2')(x, y), Symbol('L3')(x, y)

# Construct replacements
Ls = [L1, L2, L3]
Lfs = [L1f, L2f, L3f]
replacements = []
for LX, LXf in zip(Ls, Lfs):
    replacements.append((LX.diff(x, x), LXf.diff(x, x)))
    replacements.append((LX.diff(y, y), LXf.diff(y, y)))
    replacements.append((LX.diff(x, y), LXf.diff(x, y)))
    replacements.append((LX.diff(x), LXf.diff(x)))
    replacements.append((LX.diff(y), LXf.diff(y)))
    replacements.append((LX, 1/S(3)))
replacements.append((x, (x1 + x2 + x3)/3))
replacements.append((y, (y1 + y2 + y3)/3))

def print_code(name, func, index):
    code = ['# %s' % name]
    for i, phii in enumerate(phi):
        v = func(phii)
        v = v.subs(replacements)
        v = v.simplify()
        code.append('A[icell, %d, %d] = %s' % (index, i, v))
    code.append('')
    print('\n'.join(code))

print('########################################')
print('Degree 1 in 2D:\n')

# Barycentric linear shape functions on a triangle
phi = [L1, L2, L3]

print_code('Value at xc (also the cell average value)', lambda q: q, 0)
print_code('d/dx', lambda q: q.diff(x), 1)
print_code('d/dy', lambda q: q.diff(y), 2)

print('########################################')
print('Degree 2 in 2D:\n')

# Barycentric quadratic shape functions on a triangle
# see documentation/notebooks/barycentric.ipynb
phi = [None]*6
phi[0] = L1*(2*L1 - 1)
phi[1] = L2*(2*L2 - 1)
phi[2] = L3*(2*L3 - 1)
phi[3] = 4*L2*L3
phi[4] = 4*L1*L3
phi[5] = 4*L1*L2

#print_code('Value at xc (MODIFY ME TO GET THE CELL AVERAGE)', lambda q: q, 0)
print('# Cell average value')
print('A[icell, 0, 3] = 1/3')
print('A[icell, 0, 4] = 1/3')
print('A[icell, 0, 5] = 1/3\n')
print_code('d/dx', lambda q: q.diff(x), 1)
print_code('d/dy', lambda q: q.diff(y), 2)
print_code('d/dx^2', lambda q: q.diff(x, x), 3)
print_code('d/dy^2', lambda q: q.diff(y, y), 4)
print_code('d/dx*dy', lambda q: q.diff(x, y), 5)

# Cell average terms
BT_QUADRATURE_TABLE = {}
BT_QUADRATURE_TABLE[1] = ([[1/S(3), 1/S(3)]],
                          [1.0])
BT_QUADRATURE_TABLE[2] = ([[2/S(3), 1/S(6)], [1/S(6), 2/S(3)], [1/S(6), 1/S(6)]],
                          [1/S(3), 1/S(3), 1/S(3)])
BT_QUADRATURE_TABLE[3] = ([[1/S(3), 1/S(3)], [S(3)/5, S(1)/5], [S(1)/5, S(3)/5], [S(1)/5, S(1)/5]],
                          [-S(9)/16, S(25)/48, S(25)/48, S(25)/48])

xv = x1*L1 + x2*L2 + x3*L3
yv = y1*L1 + y2*L2 + y3*L3
xc = xv.subs({L1: 1/S(3), L2: 1/S(3), L3: 1/S(3)})
yc = yv.subs({L1: 1/S(3), L2: 1/S(3), L3: 1/S(3)})

exprs = [('bar_xx', (xv-xc)**2/2),
         ('bar_yy', (yv-yc)**2/2),
         ('bar_xy', (xv-xc)*(yv-yc))]
for i in range(6):
    exprs.append(('A[icell, 0, %d]' % i, phi[i]))

order = 3
points, weights = BT_QUADRATURE_TABLE[order]
for name, expr in exprs:
    v = 0
    for point, weight in zip(points, weights):
        loc = {L1: point[0],
               L2: point[1],
               L3: 1 - point[0] - point[1]}
        v += weight*expr.subs(loc)
    v = v.simplify().simplify()
    v = v.subs({x1: Symbol('x1'), x2: Symbol('x2'), x3: Symbol('x3'),
                y1: Symbol('y1'), y2: Symbol('y2'), y3: Symbol('y3')})
    print('%s = %s' % (name, v))

```

    ########################################
    Degree 1 in 2D:
    
    # Value at xc (also the cell average value)
    A[icell, 0, 0] = 1/3
    A[icell, 0, 1] = 1/3
    A[icell, 0, 2] = 1/3
    
    # d/dx
    A[icell, 1, 0] = (y2 - y3)/D
    A[icell, 1, 1] = (-y1 + y3)/D
    A[icell, 1, 2] = (y1 - y2)/D
    
    # d/dy
    A[icell, 2, 0] = (-x2 + x3)/D
    A[icell, 2, 1] = (x1 - x3)/D
    A[icell, 2, 2] = (-x1 + x2)/D
    
    ########################################
    Degree 2 in 2D:
    
    # Cell average value
    A[icell, 0, 3] = 1/3
    A[icell, 0, 4] = 1/3
    A[icell, 0, 5] = 1/3
    
    # d/dx
    A[icell, 1, 0] = (y2 - y3)/(3*D)
    A[icell, 1, 1] = (-y1 + y3)/(3*D)
    A[icell, 1, 2] = (y1 - y2)/(3*D)
    A[icell, 1, 3] = 4*(-y2 + y3)/(3*D)
    A[icell, 1, 4] = 4*(y1 - y3)/(3*D)
    A[icell, 1, 5] = 4*(-y1 + y2)/(3*D)
    
    # d/dy
    A[icell, 2, 0] = (-x2 + x3)/(3*D)
    A[icell, 2, 1] = (x1 - x3)/(3*D)
    A[icell, 2, 2] = (-x1 + x2)/(3*D)
    A[icell, 2, 3] = 4*(x2 - x3)/(3*D)
    A[icell, 2, 4] = 4*(-x1 + x3)/(3*D)
    A[icell, 2, 5] = 4*(x1 - x2)/(3*D)
    
    # d/dx^2
    A[icell, 3, 0] = 4*(y2 - y3)**2/D**2
    A[icell, 3, 1] = 4*(y1 - y3)**2/D**2
    A[icell, 3, 2] = 4*(-y1 + y2)**2/D**2
    A[icell, 3, 3] = 8*(-y1 + y2)*(y1 - y3)/D**2
    A[icell, 3, 4] = 8*(y1 - y2)*(y2 - y3)/D**2
    A[icell, 3, 5] = -8*(y1 - y3)*(y2 - y3)/D**2
    
    # d/dy^2
    A[icell, 4, 0] = 4*(x2 - x3)**2/D**2
    A[icell, 4, 1] = 4*(x1 - x3)**2/D**2
    A[icell, 4, 2] = 4*(x1 - x2)**2/D**2
    A[icell, 4, 3] = -8*(x1 - x2)*(x1 - x3)/D**2
    A[icell, 4, 4] = 8*(x1 - x2)*(x2 - x3)/D**2
    A[icell, 4, 5] = -8*(x1 - x3)*(x2 - x3)/D**2
    
    # d/dx*dy
    A[icell, 5, 0] = -4*(x2 - x3)*(y2 - y3)/D**2
    A[icell, 5, 1] = -4*(x1 - x3)*(y1 - y3)/D**2
    A[icell, 5, 2] = 4*(x1 - x2)*(-y1 + y2)/D**2
    A[icell, 5, 3] = 4*((x1 - x2)*(y1 - y3) + (x1 - x3)*(y1 - y2))/D**2
    A[icell, 5, 4] = -(4*(x1 - x2)*(y2 - y3) + 4*(x2 - x3)*(y1 - y2))/D**2
    A[icell, 5, 5] = 4*((x1 - x3)*(y2 - y3) + (x2 - x3)*(y1 - y3))/D**2
    
    bar_xx = x1**2/36 - x1*x2/36 - x1*x3/36 + x2**2/36 - x2*x3/36 + x3**2/36
    bar_yy = y1**2/36 - y1*y2/36 - y1*y3/36 + y2**2/36 - y2*y3/36 + y3**2/36
    bar_xy = x1*y1/18 - x1*y2/36 - x1*y3/36 - x2*y1/36 + x2*y2/18 - x2*y3/36 - x3*y1/36 - x3*y2/36 + x3*y3/18
    A[icell, 0, 0] = 0
    A[icell, 0, 1] = 0
    A[icell, 0, 2] = 0
    A[icell, 0, 3] = 1/3
    A[icell, 0, 4] = 1/3
    A[icell, 0, 5] = 1/3


## Code generation for transformation matrices, 3D

Use sympy to derive the coefficients needed to convert from DG polynomials
on a tetrahedron to a Taylor basis by expressing the derivatives at the centre
point of the triangle


```python
%reset -f

from sympy import Symbol, symbols, S, sqrt
x1, x2, x3, x4 = symbols('x1 x2 x3 x4', real=True)
y1, y2, y3, y4 = symbols('y1 y2 y3 y4', real=True)
z1, z2, z3, z4 = symbols('z1 z2 z3 z4', real=True)
x, y, z, F = symbols('x y z F', real=True)

#F = x1*y2*z3 - x1*y2*z4 - x1*y3*z2 + x1*y3*z4 + x1*y4*z2 - x1*y4*z3 - x2*y1*z3 + x2*y1*z4 + x2*y3*z1 - x2*y3*z4 - x2*y4*z1 + x2*y4*z3 + x3*y1*z2 - x3*y1*z4 - x3*y2*z1 + x3*y2*z4 + x3*y4*z1 - x3*y4*z2 - x4*y1*z2 + x4*y1*z3 + x4*y2*z1 - x4*y2*z3 - x4*y3*z1 + x4*y3*z2
L1f = (x*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) - x2*y3*z4 + x2*y4*z3 + x3*y2*z4 - x3*y4*z2 - x4*y2*z3 + x4*y3*z2 - y*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3) + z*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3))/F
L2f = (-x*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + x1*y3*z4 - x1*y4*z3 - x3*y1*z4 + x3*y4*z1 + x4*y1*z3 - x4*y3*z1 + y*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3) - z*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3))/F
L3f = (x*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2) - x1*y2*z4 + x1*y4*z2 + x2*y1*z4 - x2*y4*z1 - x4*y1*z2 + x4*y2*z1 - y*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2) + z*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2))/F
L4f = (-x*y1*z2 + x*y1*z3 + x*y2*z1 - x*y2*z3 - x*y3*z1 + x*y3*z2 + x1*y*z2 - x1*y*z3 - x1*y2*z + x1*y2*z3 + x1*y3*z - x1*y3*z2 - x2*y*z1 + x2*y*z3 + x2*y1*z - x2*y1*z3 - x2*y3*z + x2*y3*z1 + x3*y*z1 - x3*y*z2 - x3*y1*z + x3*y1*z2 + x3*y2*z - x3*y2*z1)/F 
L1, L2, L3, L4 = Symbol('L1')(x, y, z), Symbol('L2')(x, y, z), Symbol('L3')(x, y, z), Symbol('L4')(x, y, z)

# Construct replacements
Ls = [L1, L2, L3, L4]
Lfs = [L1f, L2f, L3f, L4f]
replacements = []
for LX, LXf in zip(Ls, Lfs):
    replacements.append((LX.diff(x, x), LXf.diff(x, x)))
    replacements.append((LX.diff(y, y), LXf.diff(y, y)))
    replacements.append((LX.diff(z, z), LXf.diff(z, z)))
    replacements.append((LX.diff(x, y), LXf.diff(x, y)))
    replacements.append((LX.diff(x, z), LXf.diff(x, z)))
    replacements.append((LX.diff(y, z), LXf.diff(y, z)))
    replacements.append((LX.diff(x), LXf.diff(x)))
    replacements.append((LX.diff(y), LXf.diff(y)))
    replacements.append((LX.diff(z), LXf.diff(z)))
    replacements.append((LX, 1/S(4)))
replacements.append((x, (x1 + x2 + x3 + x4)/4))
replacements.append((y, (y1 + y2 + y3 + y4)/4))
replacements.append((z, (z1 + z2 + z3 + z4)/4))

def print_code(name, func, index):
    code = ['# %s' % name]
    for i, phii in enumerate(phi):
        v = func(phii)
        v = v.subs(replacements)
        v = v.simplify()
        vstr = str(v).replace('Abs', 'abs')
        code.append('A[icell, %d, %d] = %s' % (index, i, vstr))
    code.append('')
    print('\n'.join(code))

print('########################################')
print('Degree 1 in 3D:\n')

# Barycentric linear shape functions on a tetrahedron
phi = [L1, L2, L3, L4]

print_code('Value at xc (also the cell average value)', lambda q: q, 0)
print_code('d/dx', lambda q: q.diff(x), 1)
print_code('d/dy', lambda q: q.diff(y), 2)
print_code('d/dz', lambda q: q.diff(z), 3)

print('########################################')
print('Degree 2 in 3D:\n')

# Barycentric quadratic shape functions on a tetrahedron
# see documentation/notebooks/barycentric.ipynb
phi = [None]*10
phi[0] = L1*(2*L1 - 1)
phi[1] = L2*(2*L2 - 1)
phi[2] = L3*(2*L3 - 1)
phi[3] = L4*(2*L4 - 1)
phi[4] = 4*L3*L4
phi[5] = 4*L2*L4
phi[6] = 4*L2*L3
phi[7] = 4*L1*L4
phi[8] = 4*L1*L3
phi[9] = 4*L1*L2

#print_code('Value at xc', lambda q: q, 0)
print('# Cell average value')
print('A[icell, 0, 0] = -1/20')
print('A[icell, 0, 1] = -1/20')
print('A[icell, 0, 2] = -1/20')
print('A[icell, 0, 3] = -1/20')
print('A[icell, 0, 4] = 1/5')
print('A[icell, 0, 5] = 1/5')
print('A[icell, 0, 6] = 1/5')
print('A[icell, 0, 7] = 1/5')
print('A[icell, 0, 8] = 1/5')
print('A[icell, 0, 9] = 1/5\n')

print_code('d/dx', lambda q: q.diff(x), 1)
print_code('d/dy', lambda q: q.diff(y), 2)
print_code('d/dz', lambda q: q.diff(z), 3)
print_code('d/dx^2', lambda q: q.diff(x, x), 4)
print_code('d/dy^2', lambda q: q.diff(y, y), 5)
print_code('d/dz^2', lambda q: q.diff(z, z), 6)
print_code('d/dx*dy', lambda q: q.diff(x, y), 7)
print_code('d/dx*dz', lambda q: q.diff(x, z), 8)
print_code('d/dy*dz', lambda q: q.diff(y, z), 9)

xv = x1*L1 + x2*L2 + x3*L3 + x4*L4
yv = y1*L1 + y2*L2 + y3*L3 + y4*L4
zv = z1*L1 + z2*L2 + z3*L3 + z4*L4
xc = xv.subs({L1: 1/S(4), L2: 1/S(4), L3: 1/S(4), L4: 1/S(4)})
yc = yv.subs({L1: 1/S(4), L2: 1/S(4), L3: 1/S(4), L4: 1/S(4)})
zc = zv.subs({L1: 1/S(4), L2: 1/S(4), L3: 1/S(4), L4: 1/S(4)})

exprs = [('bar_xx', (xv-xc)**2/2),
         ('bar_yy', (yv-yc)**2/2),
         ('bar_zz', (zv-zc)**2/2),
         ('bar_xy', (xv-xc)*(yv-yc)),
         ('bar_xz', (xv-xc)*(zv-zc)),
         ('bar_yz', (yv-yc)*(zv-zc))]
for i in range(10):
    exprs.append(('A[icell, 0, %d]' % i, phi[i]))

# see, e.g., https://github.com/libMesh/libmesh/blob/master/src/quadrature/quadrature_gauss_3D.C
BT_QUADRATURE_TABLE = {}
I = S(1)
a = 1/(3*sqrt(5) - 5)
b = 1/(1*sqrt(5) + 5)
BT_QUADRATURE_TABLE[1] = ([[I/4, I/4, I/4]], [I/6])
BT_QUADRATURE_TABLE[2] = ([[a, b, b], [b, a, b], [b, b, a], [b, b, b]], [I/24, I/24, I/24, I/24])
BT_QUADRATURE_TABLE[3] = ([[I/4, I/4, I/4],
                           [I/2, I/6, I/6],
                           [I/6, I/2, I/6],
                           [I/6, I/6, I/2],
                           [I/6, I/6, I/6]],
                          [S(-2)/15, S(3)/40, S(3)/40, S(3)/40, S(3)/40])

order = 3
points, weights = BT_QUADRATURE_TABLE[order]
for name, expr in exprs:
    v = 0
    for point, weight in zip(points, weights):
        loc = {L1: point[0],
               L2: point[1],
               L3: point[2],
               L4: 1 - point[0] - point[1] - point[2]}
        v += 6*weight*expr.subs(loc).simplify()
    v = v.simplify().simplify()
    v = v.subs({x1: Symbol('x1'), x2: Symbol('x2'), x3: Symbol('x3'), x4: Symbol('x4'),
                y1: Symbol('y1'), y2: Symbol('y2'), y3: Symbol('y3'), y4: Symbol('y4'),
                z1: Symbol('z1'), z2: Symbol('z2'), z3: Symbol('z3'), z4: Symbol('z4')})
    print('%s = %s' % (name, v))
```

    ########################################
    Degree 1 in 3D:
    
    # Value at xc (also the cell average value)
    A[icell, 0, 0] = 1/4
    A[icell, 0, 1] = 1/4
    A[icell, 0, 2] = 1/4
    A[icell, 0, 3] = 1/4
    
    # d/dx
    A[icell, 1, 0] = (y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)/F
    A[icell, 1, 1] = (-y1*z3 + y1*z4 + y3*z1 - y3*z4 - y4*z1 + y4*z3)/F
    A[icell, 1, 2] = (y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)/F
    A[icell, 1, 3] = (-y1*z2 + y1*z3 + y2*z1 - y2*z3 - y3*z1 + y3*z2)/F
    
    # d/dy
    A[icell, 2, 0] = (-x2*z3 + x2*z4 + x3*z2 - x3*z4 - x4*z2 + x4*z3)/F
    A[icell, 2, 1] = (x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)/F
    A[icell, 2, 2] = (-x1*z2 + x1*z4 + x2*z1 - x2*z4 - x4*z1 + x4*z2)/F
    A[icell, 2, 3] = (x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)/F
    
    # d/dz
    A[icell, 3, 0] = (x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)/F
    A[icell, 3, 1] = (-x1*y3 + x1*y4 + x3*y1 - x3*y4 - x4*y1 + x4*y3)/F
    A[icell, 3, 2] = (x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)/F
    A[icell, 3, 3] = (-x1*y2 + x1*y3 + x2*y1 - x2*y3 - x3*y1 + x3*y2)/F
    
    ########################################
    Degree 2 in 3D:
    
    # Cell average value
    A[icell, 0, 0] = -1/20
    A[icell, 0, 1] = -1/20
    A[icell, 0, 2] = -1/20
    A[icell, 0, 3] = -1/20
    A[icell, 0, 4] = 1/5
    A[icell, 0, 5] = 1/5
    A[icell, 0, 6] = 1/5
    A[icell, 0, 7] = 1/5
    A[icell, 0, 8] = 1/5
    A[icell, 0, 9] = 1/5
    
    # d/dx
    A[icell, 1, 0] = 0
    A[icell, 1, 1] = 0
    A[icell, 1, 2] = 0
    A[icell, 1, 3] = 0
    A[icell, 1, 4] = (y1*z3 - y1*z4 - y2*z3 + y2*z4 - y3*z1 + y3*z2 + y4*z1 - y4*z2)/F
    A[icell, 1, 5] = (-y1*z2 + y1*z4 + y2*z1 - y2*z3 + y3*z2 - y3*z4 - y4*z1 + y4*z3)/F
    A[icell, 1, 6] = (y1*z2 - y1*z3 - y2*z1 + y2*z4 + y3*z1 - y3*z4 - y4*z2 + y4*z3)/F
    A[icell, 1, 7] = (-y1*z2 + y1*z3 + y2*z1 - y2*z4 - y3*z1 + y3*z4 + y4*z2 - y4*z3)/F
    A[icell, 1, 8] = (y1*z2 - y1*z4 - y2*z1 + y2*z3 - y3*z2 + y3*z4 + y4*z1 - y4*z3)/F
    A[icell, 1, 9] = (-y1*z3 + y1*z4 + y2*z3 - y2*z4 + y3*z1 - y3*z2 - y4*z1 + y4*z2)/F
    
    # d/dy
    A[icell, 2, 0] = 0
    A[icell, 2, 1] = 0
    A[icell, 2, 2] = 0
    A[icell, 2, 3] = 0
    A[icell, 2, 4] = (-x1*z3 + x1*z4 + x2*z3 - x2*z4 + x3*z1 - x3*z2 - x4*z1 + x4*z2)/F
    A[icell, 2, 5] = (x1*z2 - x1*z4 - x2*z1 + x2*z3 - x3*z2 + x3*z4 + x4*z1 - x4*z3)/F
    A[icell, 2, 6] = (-x1*z2 + x1*z3 + x2*z1 - x2*z4 - x3*z1 + x3*z4 + x4*z2 - x4*z3)/F
    A[icell, 2, 7] = (x1*z2 - x1*z3 - x2*z1 + x2*z4 + x3*z1 - x3*z4 - x4*z2 + x4*z3)/F
    A[icell, 2, 8] = (-x1*z2 + x1*z4 + x2*z1 - x2*z3 + x3*z2 - x3*z4 - x4*z1 + x4*z3)/F
    A[icell, 2, 9] = (x1*z3 - x1*z4 - x2*z3 + x2*z4 - x3*z1 + x3*z2 + x4*z1 - x4*z2)/F
    
    # d/dz
    A[icell, 3, 0] = 0
    A[icell, 3, 1] = 0
    A[icell, 3, 2] = 0
    A[icell, 3, 3] = 0
    A[icell, 3, 4] = (x1*y3 - x1*y4 - x2*y3 + x2*y4 - x3*y1 + x3*y2 + x4*y1 - x4*y2)/F
    A[icell, 3, 5] = (-x1*y2 + x1*y4 + x2*y1 - x2*y3 + x3*y2 - x3*y4 - x4*y1 + x4*y3)/F
    A[icell, 3, 6] = (x1*y2 - x1*y3 - x2*y1 + x2*y4 + x3*y1 - x3*y4 - x4*y2 + x4*y3)/F
    A[icell, 3, 7] = (-x1*y2 + x1*y3 + x2*y1 - x2*y4 - x3*y1 + x3*y4 + x4*y2 - x4*y3)/F
    A[icell, 3, 8] = (x1*y2 - x1*y4 - x2*y1 + x2*y3 - x3*y2 + x3*y4 + x4*y1 - x4*y3)/F
    A[icell, 3, 9] = (-x1*y3 + x1*y4 + x2*y3 - x2*y4 + x3*y1 - x3*y2 - x4*y1 + x4*y2)/F
    
    # d/dx^2
    A[icell, 4, 0] = 4*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)**2/F**2
    A[icell, 4, 1] = 4*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3)**2/F**2
    A[icell, 4, 2] = 4*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)**2/F**2
    A[icell, 4, 3] = 4*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2)**2/F**2
    A[icell, 4, 4] = -8*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)/F**2
    A[icell, 4, 5] = 8*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3)/F**2
    A[icell, 4, 6] = -8*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3)/F**2
    A[icell, 4, 7] = -8*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)/F**2
    A[icell, 4, 8] = 8*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)/F**2
    A[icell, 4, 9] = -8*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)/F**2
    
    # d/dy^2
    A[icell, 5, 0] = 4*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)**2/F**2
    A[icell, 5, 1] = 4*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)**2/F**2
    A[icell, 5, 2] = 4*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)**2/F**2
    A[icell, 5, 3] = 4*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)**2/F**2
    A[icell, 5, 4] = -8*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)/F**2
    A[icell, 5, 5] = 8*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)/F**2
    A[icell, 5, 6] = -8*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)/F**2
    A[icell, 5, 7] = -8*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)/F**2
    A[icell, 5, 8] = 8*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)/F**2
    A[icell, 5, 9] = -8*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)/F**2
    
    # d/dz^2
    A[icell, 6, 0] = 4*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)**2/F**2
    A[icell, 6, 1] = 4*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)**2/F**2
    A[icell, 6, 2] = 4*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)**2/F**2
    A[icell, 6, 3] = 4*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)**2/F**2
    A[icell, 6, 4] = -8*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)/F**2
    A[icell, 6, 5] = 8*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)/F**2
    A[icell, 6, 6] = -8*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)/F**2
    A[icell, 6, 7] = -8*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)/F**2
    A[icell, 6, 8] = 8*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)/F**2
    A[icell, 6, 9] = -8*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)/F**2
    
    # d/dx*dy
    A[icell, 7, 0] = -4*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)/F**2
    A[icell, 7, 1] = -4*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3)/F**2
    A[icell, 7, 2] = -4*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)/F**2
    A[icell, 7, 3] = -4*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2)/F**2
    A[icell, 7, 4] = 4*((x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2) + (x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2))/F**2
    A[icell, 7, 5] = -(4*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + 4*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2))/F**2
    A[icell, 7, 6] = 4*((x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + (x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2))/F**2
    A[icell, 7, 7] = 4*((x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) + (x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2))/F**2
    A[icell, 7, 8] = -(4*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) + 4*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2))/F**2
    A[icell, 7, 9] = 4*((x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) + (x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3))/F**2
    
    # d/dx*dz
    A[icell, 8, 0] = 4*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3)/F**2
    A[icell, 8, 1] = 4*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3)/F**2
    A[icell, 8, 2] = 4*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2)/F**2
    A[icell, 8, 3] = 4*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2)/F**2
    A[icell, 8, 4] = -(4*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2) + 4*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2))/F**2
    A[icell, 8, 5] = 4*((x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + (x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2))/F**2
    A[icell, 8, 6] = -(4*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3) + 4*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2))/F**2
    A[icell, 8, 7] = -(4*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) + 4*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)*(y1*z2 - y1*z3 - y2*z1 + y2*z3 + y3*z1 - y3*z2))/F**2
    A[icell, 8, 8] = 4*((x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) + (x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)*(y1*z2 - y1*z4 - y2*z1 + y2*z4 + y4*z1 - y4*z2))/F**2
    A[icell, 8, 9] = -(4*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(y2*z3 - y2*z4 - y3*z2 + y3*z4 + y4*z2 - y4*z3) + 4*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)*(y1*z3 - y1*z4 - y3*z1 + y3*z4 + y4*z1 - y4*z3))/F**2
    
    # d/dy*dz
    A[icell, 9, 0] = -4*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3)/F**2
    A[icell, 9, 1] = -4*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)/F**2
    A[icell, 9, 2] = -4*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)/F**2
    A[icell, 9, 3] = -4*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)/F**2
    A[icell, 9, 4] = 4*((x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2) + (x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2))/F**2
    A[icell, 9, 5] = -(4*(x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3) + 4*(x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2))/F**2
    A[icell, 9, 6] = 4*((x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3) + (x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2))/F**2
    A[icell, 9, 7] = 4*((x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3) + (x1*z2 - x1*z3 - x2*z1 + x2*z3 + x3*z1 - x3*z2)*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3))/F**2
    A[icell, 9, 8] = -(4*(x1*y2 - x1*y4 - x2*y1 + x2*y4 + x4*y1 - x4*y2)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3) + 4*(x1*z2 - x1*z4 - x2*z1 + x2*z4 + x4*z1 - x4*z2)*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3))/F**2
    A[icell, 9, 9] = 4*((x1*y3 - x1*y4 - x3*y1 + x3*y4 + x4*y1 - x4*y3)*(x2*z3 - x2*z4 - x3*z2 + x3*z4 + x4*z2 - x4*z3) + (x1*z3 - x1*z4 - x3*z1 + x3*z4 + x4*z1 - x4*z3)*(x2*y3 - x2*y4 - x3*y2 + x3*y4 + x4*y2 - x4*y3))/F**2
    
    bar_xx = (-3*x1 + x2 + x3 + x4)**2/640 + (x1 - 3*x2 + x3 + x4)**2/640 + (x1 + x2 - 3*x3 + x4)**2/640 + (x1 + x2 + x3 - 3*x4)**2/640
    bar_yy = (-3*y1 + y2 + y3 + y4)**2/640 + (y1 - 3*y2 + y3 + y4)**2/640 + (y1 + y2 - 3*y3 + y4)**2/640 + (y1 + y2 + y3 - 3*y4)**2/640
    bar_zz = (-3*z1 + z2 + z3 + z4)**2/640 + (z1 - 3*z2 + z3 + z4)**2/640 + (z1 + z2 - 3*z3 + z4)**2/640 + (z1 + z2 + z3 - 3*z4)**2/640
    bar_xy = 3*x1*y1/80 - x1*y2/80 - x1*y3/80 - x1*y4/80 - x2*y1/80 + 3*x2*y2/80 - x2*y3/80 - x2*y4/80 - x3*y1/80 - x3*y2/80 + 3*x3*y3/80 - x3*y4/80 - x4*y1/80 - x4*y2/80 - x4*y3/80 + 3*x4*y4/80
    bar_xz = 3*x1*z1/80 - x1*z2/80 - x1*z3/80 - x1*z4/80 - x2*z1/80 + 3*x2*z2/80 - x2*z3/80 - x2*z4/80 - x3*z1/80 - x3*z2/80 + 3*x3*z3/80 - x3*z4/80 - x4*z1/80 - x4*z2/80 - x4*z3/80 + 3*x4*z4/80
    bar_yz = 3*y1*z1/80 - y1*z2/80 - y1*z3/80 - y1*z4/80 - y2*z1/80 + 3*y2*z2/80 - y2*z3/80 - y2*z4/80 - y3*z1/80 - y3*z2/80 + 3*y3*z3/80 - y3*z4/80 - y4*z1/80 - y4*z2/80 - y4*z3/80 + 3*y4*z4/80
    A[icell, 0, 0] = -1/20
    A[icell, 0, 1] = -1/20
    A[icell, 0, 2] = -1/20
    A[icell, 0, 3] = -1/20
    A[icell, 0, 4] = 1/5
    A[icell, 0, 5] = 1/5
    A[icell, 0, 6] = 1/5
    A[icell, 0, 7] = 1/5
    A[icell, 0, 8] = 1/5
    A[icell, 0, 9] = 1/5



```python

```

649be102524024e7ebbbe9c99b28114e56d827f8
ipynb
Jupyter Notebook
documentation/notebooks/barycentric.ipynb
TormodLandet/Ocellaris
6b4b2515fb881b1ed8d8fd8d8c23a8e1990ada58
Apache-2.0
2017-11-07T12:19:44.000Z
2017-11-07T12:19:44.000Z
documentation/notebooks/barycentric.ipynb
TormodLandet/Ocellaris
6b4b2515fb881b1ed8d8fd8d8c23a8e1990ada58
Apache-2.0
documentation/notebooks/barycentric.ipynb
TormodLandet/Ocellaris
6b4b2515fb881b1ed8d8fd8d8c23a8e1990ada58
Apache-2.0
2018-05-02T17:17:01.000Z
2019-03-11T13:09:40.000Z
Qwen/Qwen-72B
1. YES
2. YES