## Generating finite element sampling functions with sympy 

This notebook describes how to use `sympy` to generate symbolic expressions describing sampling of finite elements of different order. It starts with some simple 1D examples and then moves to 2D and 3D and ends with a sampling method for the `VTK_LAGRANGE_HEXAHEDRON` VTK cell type. 


```python
import sympy as sy
import numpy as np 

%matplotlib notebook
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d 
```

## 1D elements

Let's start with the 1D lagrange polynomial basis functions! 

https://www.longqi.cf/python/2014/03/24/implement-of-lagrange-polynomial-in-sympy/ has a nice implementation already: 


```python
def LagrangPoly(x,order,i,xi=None):
    if xi==None:
        xi=sy.symbols('x:%d'%(order+1))
    index = range(order+1)    
    return sy.prod([(x-xi[j])/(xi[i]-xi[j]) for j in index if j != i])
```

where 𝑥 is a sympy symbol, order is the polynomial order (1 for linear, 2 for quadratic, 3 for cubic, etc.), 𝑖 the node index and 𝑥𝑖 are the node locations. If we run with order=2 and i=0 without node locations, we'll get the polynomial expression for 𝐿0: 


```python
x=sy.symbols('x')
LagrangPoly(x,2,0)
```




$\displaystyle \frac{\left(x - x_{1}\right) \left(x - x_{2}\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)}$



For the rest of the notebook, we're going to use [-1,0,1] for the node locations in each dimension, so let's supply those to get a simpler expression: 


```python
LP = LagrangPoly(x,2,0,[-1,0,1])
LP
```




$\displaystyle - x \left(\frac{1}{2} - \frac{x}{2}\right)$




```python
sy.simplify(LP)
```




$\displaystyle \frac{x \left(x - 1\right)}{2}$



The 2 basis functions for a 1D linear element are 


```python
sy.pprint(
    (sy.simplify(LagrangPoly(x,1,0,[-1,1])),
    sy.simplify(LagrangPoly(x,1,1,[-1,1])))
)

```

    ⎛1   x  x   1⎞
    ⎜─ - ─, ─ + ─⎟
    ⎝2   2  2   2⎠


and then the 3 basis functions for a 1D quadratic element are 


```python
sy.pprint(
    (sy.simplify(LagrangPoly(x,2,0,[-1,0,1])),
    sy.simplify(LagrangPoly(x,2,1,[-1,0,1])),
    sy.simplify(LagrangPoly(x,2,2,[-1,0,1]))
    )
)

```

    ⎛x⋅(x - 1)       2  x⋅(x + 1)⎞
    ⎜─────────, 1 - x , ─────────⎟
    ⎝    2                  2    ⎠


To get the functional form for evaluating at a point off a node, we also need the known values of the element at the nodes:


```python
vals = [sy.symbols('f'+str(i)) for i in range(3)]
vals
```




    [f0, f1, f2]



So our in-cell interpolation function for the linear element is 


```python
shape_funcs = [
    sy.simplify(LagrangPoly(x,1,0,[-1,1])),
    sy.simplify(LagrangPoly(x,1,1,[-1,1]))
]
sample_expression = sum([vals[i]*shape_funcs[i] for i in range(2)])
sample_expression
```




$\displaystyle f_{0} \left(\frac{1}{2} - \frac{x}{2}\right) + f_{1} \left(\frac{x}{2} + \frac{1}{2}\right)$



And for the quadratic: 


```python
shape_funcs = [
    sy.simplify(LagrangPoly(x,2,0,[-1,0,1])),
    sy.simplify(LagrangPoly(x,2,1,[-1,0,1])),
    sy.simplify(LagrangPoly(x,2,2,[-1,0,1]))
]
sample_expression = sum([vals[i]*shape_funcs[i] for i in range(3)])
sample_expression
```




$\displaystyle \frac{f_{0} x \left(x - 1\right)}{2} + f_{1} \left(1 - x^{2}\right) + \frac{f_{2} x \left(x + 1\right)}{2}$



## higher dimensions 

To construct the basis functions for higher dimensions, we simply multiply on another lagrange basis function for the new dimension. 

### linear rectangular element 
$2^2$ total nodes (4 corners): 


```python

y=sy.symbols('y')
shape_funcs = []

for x_i in range(2):
    for y_i in range(2):
        LP1 = LagrangPoly(x,1,x_i,[-1,1])
        LP2 = LagrangPoly(y,1,y_i,[-1,1])
        shape_funcs.append(sy.simplify(LP2*LP1))
        
sy.pprint(shape_funcs)
```

    ⎡(x - 1)⋅(y - 1)  -(x - 1)⋅(y + 1)   -(x + 1)⋅(y - 1)   (x + 1)⋅(y + 1)⎤
    ⎢───────────────, ─────────────────, ─────────────────, ───────────────⎥
    ⎣       4                 4                  4                 4       ⎦


### quadratic rectangular element 

$3^3$ total nodes (4 corners, 4 edge-centers and 1 area-center):


```python
shape_funcs = []

for x_i in range(3):
    for y_i in range(3):
        LP1 = LagrangPoly(x,2,x_i,[-1,0,1])
        LP2 = LagrangPoly(y,2,y_i,[-1,0,1])
        shape_funcs.append(sy.simplify(LP2*LP1))

for sf in shape_funcs:
    sy.pprint(sf)
```

    x⋅y⋅(x - 1)⋅(y - 1)
    ───────────────────
             4         
    -x⋅(x - 1)⋅(y - 1)⋅(y + 1) 
    ───────────────────────────
                 2             
    x⋅y⋅(x - 1)⋅(y + 1)
    ───────────────────
             4         
    -y⋅(x - 1)⋅(x + 1)⋅(y - 1) 
    ───────────────────────────
                 2             
    (x - 1)⋅(x + 1)⋅(y - 1)⋅(y + 1)
    -y⋅(x - 1)⋅(x + 1)⋅(y + 1) 
    ───────────────────────────
                 2             
    x⋅y⋅(x + 1)⋅(y - 1)
    ───────────────────
             4         
    -x⋅(x + 1)⋅(y - 1)⋅(y + 1) 
    ───────────────────────────
                 2             
    x⋅y⋅(x + 1)⋅(y + 1)
    ───────────────────
             4         


### 3D linear hexahedral element 

$2^3=8$ total nodes (8 corner vertices): 


```python
z=sy.symbols('z')
shape_funcs = []

for z_i in range(2):
    for y_i in range(2):
        for x_i in range(2):        
            LP1 = LagrangPoly(x,1,x_i,[-1,1])
            LP2 = LagrangPoly(y,1,y_i,[-1,1])
            LP3 = LagrangPoly(z,1,z_i,[-1,1])
            shape_funcs.append(sy.simplify(LP1 * LP2 * LP3))
        

for sf in shape_funcs:
    sy.pprint(sf) 
```

    -(x - 1)⋅(y - 1)⋅(z - 1) 
    ─────────────────────────
                8            
    (x + 1)⋅(y - 1)⋅(z - 1)
    ───────────────────────
               8           
    (x - 1)⋅(y + 1)⋅(z - 1)
    ───────────────────────
               8           
    -(x + 1)⋅(y + 1)⋅(z - 1) 
    ─────────────────────────
                8            
    (x - 1)⋅(y - 1)⋅(z + 1)
    ───────────────────────
               8           
    -(x + 1)⋅(y - 1)⋅(z + 1) 
    ─────────────────────────
                8            
    -(x - 1)⋅(y + 1)⋅(z + 1) 
    ─────────────────────────
                8            
    (x + 1)⋅(y + 1)⋅(z + 1)
    ───────────────────────
               8           


### 3D quadratic hexahedral element 

$3^3=27$ total nodes (8 corner vertices, 12 edge-centers, 6 face-centers, 1 volume-center): 


```python
shape_funcs = []

for z_i in range(3):
    for y_i in range(3):
        for x_i in range(3):        
            LP1 = LagrangPoly(x,2,x_i,[-1,0,1])
            LP2 = LagrangPoly(y,2,y_i,[-1,0,1])
            LP3 = LagrangPoly(z,2,z_i,[-1,0,1])
            shape_funcs.append(sy.simplify(LP1 * LP2 * LP3))

for sf in shape_funcs:
    sy.pprint(sf)

```

    x⋅y⋅z⋅(x - 1)⋅(y - 1)⋅(z - 1)
    ─────────────────────────────
                  8              
    -y⋅z⋅(x - 1)⋅(x + 1)⋅(y - 1)⋅(z - 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x + 1)⋅(y - 1)⋅(z - 1)
    ─────────────────────────────
                  8              
    -x⋅z⋅(x - 1)⋅(y - 1)⋅(y + 1)⋅(z - 1) 
    ─────────────────────────────────────
                      4                  
    z⋅(x - 1)⋅(x + 1)⋅(y - 1)⋅(y + 1)⋅(z - 1)
    ─────────────────────────────────────────
                        2                    
    -x⋅z⋅(x + 1)⋅(y - 1)⋅(y + 1)⋅(z - 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x - 1)⋅(y + 1)⋅(z - 1)
    ─────────────────────────────
                  8              
    -y⋅z⋅(x - 1)⋅(x + 1)⋅(y + 1)⋅(z - 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x + 1)⋅(y + 1)⋅(z - 1)
    ─────────────────────────────
                  8              
    -x⋅y⋅(x - 1)⋅(y - 1)⋅(z - 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    y⋅(x - 1)⋅(x + 1)⋅(y - 1)⋅(z - 1)⋅(z + 1)
    ─────────────────────────────────────────
                        2                    
    -x⋅y⋅(x + 1)⋅(y - 1)⋅(z - 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    x⋅(x - 1)⋅(y - 1)⋅(y + 1)⋅(z - 1)⋅(z + 1)
    ─────────────────────────────────────────
                        2                    
    -(x - 1)⋅(x + 1)⋅(y - 1)⋅(y + 1)⋅(z - 1)⋅(z + 1)
    x⋅(x + 1)⋅(y - 1)⋅(y + 1)⋅(z - 1)⋅(z + 1)
    ─────────────────────────────────────────
                        2                    
    -x⋅y⋅(x - 1)⋅(y + 1)⋅(z - 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    y⋅(x - 1)⋅(x + 1)⋅(y + 1)⋅(z - 1)⋅(z + 1)
    ─────────────────────────────────────────
                        2                    
    -x⋅y⋅(x + 1)⋅(y + 1)⋅(z - 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x - 1)⋅(y - 1)⋅(z + 1)
    ─────────────────────────────
                  8              
    -y⋅z⋅(x - 1)⋅(x + 1)⋅(y - 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x + 1)⋅(y - 1)⋅(z + 1)
    ─────────────────────────────
                  8              
    -x⋅z⋅(x - 1)⋅(y - 1)⋅(y + 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    z⋅(x - 1)⋅(x + 1)⋅(y - 1)⋅(y + 1)⋅(z + 1)
    ─────────────────────────────────────────
                        2                    
    -x⋅z⋅(x + 1)⋅(y - 1)⋅(y + 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x - 1)⋅(y + 1)⋅(z + 1)
    ─────────────────────────────
                  8              
    -y⋅z⋅(x - 1)⋅(x + 1)⋅(y + 1)⋅(z + 1) 
    ─────────────────────────────────────
                      4                  
    x⋅y⋅z⋅(x + 1)⋅(y + 1)⋅(z + 1)
    ─────────────────────────────
                  8              


## node numbering

different FEM packages can use slightly different node number conventions. For example, the 3D quadratic hexahedral interpolation function would be given by 


```python
vals = [sy.symbols('f'+str(i)) for i in range(27)]
sample_expression = sum([vals[i]*shape_funcs[i] for i in range(27)])
sample_expression
```




$\displaystyle \frac{f_{0} x y z \left(x - 1\right) \left(y - 1\right) \left(z - 1\right)}{8} - \frac{f_{1} y z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(z - 1\right)}{4} + \frac{f_{10} y \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(z - 1\right) \left(z + 1\right)}{2} - \frac{f_{11} x y \left(x + 1\right) \left(y - 1\right) \left(z - 1\right) \left(z + 1\right)}{4} + \frac{f_{12} x \left(x - 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{2} - f_{13} \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right) + \frac{f_{14} x \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{2} - \frac{f_{15} x y \left(x - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{4} + \frac{f_{16} y \left(x - 1\right) \left(x + 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{2} - \frac{f_{17} x y \left(x + 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{4} + \frac{f_{18} x y z \left(x - 1\right) \left(y - 1\right) \left(z + 1\right)}{8} - \frac{f_{19} y z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(z + 1\right)}{4} + \frac{f_{2} x y z \left(x + 1\right) \left(y - 1\right) \left(z - 1\right)}{8} + \frac{f_{20} x y z \left(x + 1\right) \left(y - 1\right) \left(z + 1\right)}{8} - \frac{f_{21} x z \left(x - 1\right) \left(y - 1\right) \left(y + 1\right) \left(z + 1\right)}{4} + \frac{f_{22} z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z + 1\right)}{2} - \frac{f_{23} x z \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z + 1\right)}{4} + \frac{f_{24} x y z \left(x - 1\right) \left(y + 1\right) \left(z + 1\right)}{8} - \frac{f_{25} y z \left(x - 1\right) \left(x + 1\right) \left(y + 1\right) \left(z + 1\right)}{4} + \frac{f_{26} x y z \left(x + 1\right) \left(y + 1\right) \left(z + 1\right)}{8} - \frac{f_{3} x z \left(x - 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right)}{4} + \frac{f_{4} z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right)}{2} - \frac{f_{5} x z \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right)}{4} + \frac{f_{6} x y z \left(x - 1\right) \left(y + 1\right) \left(z - 1\right)}{8} - \frac{f_{7} y z \left(x - 1\right) \left(x + 1\right) \left(y + 1\right) \left(z - 1\right)}{4} + \frac{f_{8} x y z \left(x + 1\right) \left(y + 1\right) \left(z - 1\right)}{8} - \frac{f_{9} x y \left(x - 1\right) \left(y - 1\right) \left(z - 1\right) \left(z + 1\right)}{4}$



and while the first 8 nodes ($f_0$ to $f_7$) generally correspond to the 8 corner vertices, there is no set convention across FEM implementations. 

page 68 of https://fenicsproject.org/pub/documents/ufc/ufc-user-manual/ufc-user-manual.pdf describes the convention for UFC (used by fenics). 

## VTK_LAGRANGE_HEXAHEDRON node numbering convention 

The node numbering for the VTK_LAGRANGE_HEXAHEDRON (VTK type 72) is buried in the VTK source code: 

https://gitlab.kitware.com/vtk/vtk/-/blob/7a0b92864c96680b1f42ee84920df556fc6ebaa3/Common/DataModel/vtkHigherOrderInterpolation.cxx


Some other useful links:
* https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
* https://github.com/Kitware/VTK/blob/0ce0d74e67927fd964a27c045d68e2f32b5f65f7/Common/DataModel/vtkCellType.h#L112
* https://github.com/ju-kreber/paraview-scripts
* https://discourse.paraview.org/t/about-high-order-non-traditional-lagrange-finite-element/1577/4

here's the VTK ordering: 

```
corners:                   edges:               

       z
7----------6            .----14----.      
|\     ^   |\           |\         |\     
| \    |   | \          | 15       | 13   
|  \   |   |  \         19 \       18 \   
|   4------+---5        |   .----12+---.  
|   |  +-- |-- | -> x   |   |      |   |  
3---+---\--2   |        .---+-10---.   |   
 \  |    \  \  |         \  16      \  17  
  \ |     \  \ |         11 |        9 |   
   \|      y  \|           \|         \|   
    0----------1            .----8-----.   
    
   
center-face node numbers 

y-z plane at x = -1 :   20 
y-z plane at x = +1 :   21
x-z plane at y = -1 :   22
x-z plane at y = +1 :   24
x-y plane at z = -1 :   23
x-y plane at z = +1 :   25

volume-center point node number: 26 
```

Note that edge numbers 18 and 19 were switched by this VTK PR https://gitlab.kitware.com/vtk/vtk/-/commit/7a0b92864c96680b1f42ee84920df556fc6ebaa3  The above numbering is for after VTK 9.0, before VTK 9.0 would require a check and switch. 
 


```python
corner_coords = [
    [-1,-1,-1],
    [ 1,-1,-1],    
    [ 1, 1,-1],
    [-1, 1,-1],
    [-1,-1, 1],
    [ 1,-1, 1],
    [ 1, 1, 1],
    [-1, 1, 1],
]

# the corner nodes defining the edges
edge_nodes = [
    [0, 1], 
    [1, 2],
    [2, 3],
    [0, 3],
    [4, 5],
    [5, 6],
    [6, 7],
    [4, 7], 
    [0, 4],
    [1, 5],
    [2, 6],
    [3, 7],
]

# the corner nodes defning the faces
face_nodes = [ 
    [0,3,4,7],
    [1,2,5,6],
    [2,3,6,7],
    [0,1,2,3],
    [0,1,4,5],
    [4,5,6,7],
]


```


```python
edge_coords = []
for edge in edge_nodes:
    edge_center = (np.array(corner_coords[edge[0]]) + np.array(corner_coords[edge[1]]))/2
    edge_coords.append(edge_center.tolist())
    
face_coords = []    
for face in face_nodes:
    coord = np.array([0,0,0])    
    for i in range(0,4):
        coord += np.array(corner_coords[face[i]]) 
    face_coords.append(coord/4)    

vol_center_coords=[np.array(corner_coords).sum(axis=0)/8]
```


```python
corner_coords=np.array(corner_coords)
edge_coords=np.array(edge_coords)
face_coords=np.array(face_coords)
vol_center_coords=np.array(vol_center_coords)
```


```python
fig=plt.figure()
ax = fig.add_subplot(111, projection='3d')

node_num=0
clrs=['r','g','b','k']
type_num=0
all_coords = []
for coords in [corner_coords, edge_coords, face_coords, vol_center_coords]:   
    
    for node in coords:
        ax.plot(node[0],node[1],node[2],marker='.',color=clrs[type_num])
        ax.text(node[0],node[1],node[2], str(node_num),fontsize=12,color=clrs[type_num])
        all_coords.append(node)
        node_num+=1 
    type_num+=1

all_coords = np.array(all_coords)
lncol=[0,0,0,.4]
for xyz in [1,-1]:
    ax.plot([-1,1],[xyz,xyz],[xyz,xyz],color=lncol)
    ax.plot([xyz,xyz],[-1,1],[xyz,xyz],color=lncol)
    ax.plot([xyz,xyz],[xyz,xyz],[-1,1],color=lncol)    
    ax.plot([-1,1],[-xyz,-xyz],[xyz,xyz],color=lncol)
    ax.plot([-xyz,-xyz],[-1,1],[xyz,xyz],color=lncol)
    ax.plot([-xyz,-xyz],[xyz,xyz],[-1,1],color=lncol)    

ax.plot([-1,1],[0.,0.],[0.,0.],'--',color=lncol)    
ax.plot([0.,0.],[-1,1],[0.,0.],'--',color=lncol)    
ax.plot([0.,0.],[0.,0.],[-1,1],'--',color=lncol)
```


    <IPython.core.display.Javascript object>









    [<mpl_toolkits.mplot3d.art3d.Line3D at 0x7fe0cb2e7b80>]



**So to build the proper interpolation function, we need to account for this node numbering:**


```python
# corresponding quadratic lagrange poly 
shape_funcs = []
vtk_node_num = []
crd=[-1,0,1]
for z_i in range(3):
    for y_i in range(3):
        for x_i in range(3):        
            LP1 = LagrangPoly(x,2,x_i,crd)
            LP2 = LagrangPoly(y,2,y_i,crd)
            LP3 = LagrangPoly(z,2,z_i,crd)
            
            # find the VTK node number             
            indx = np.where((all_coords[:,0]==crd[x_i]) & (all_coords[:,1]==crd[y_i]) & (all_coords[:,2]==crd[z_i]))[0][0]                                   
            vtk_node_num.append(indx)
            shape_funcs.append(sy.simplify(LP1 * LP2 * LP3))
            
```


```python
vals = [sy.symbols('f'+str(i)) for i in vtk_node_num]
sample_expression = sum([vals[i]*shape_funcs[i] for i in range(27)])
sample_expression
```




$\displaystyle \frac{f_{0} x y z \left(x - 1\right) \left(y - 1\right) \left(z - 1\right)}{8} + \frac{f_{1} x y z \left(x + 1\right) \left(y - 1\right) \left(z - 1\right)}{8} - \frac{f_{10} y z \left(x - 1\right) \left(x + 1\right) \left(y + 1\right) \left(z - 1\right)}{4} - \frac{f_{11} x z \left(x - 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right)}{4} - \frac{f_{12} y z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(z + 1\right)}{4} - \frac{f_{13} x z \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z + 1\right)}{4} - \frac{f_{14} y z \left(x - 1\right) \left(x + 1\right) \left(y + 1\right) \left(z + 1\right)}{4} - \frac{f_{15} x z \left(x - 1\right) \left(y - 1\right) \left(y + 1\right) \left(z + 1\right)}{4} - \frac{f_{16} x y \left(x - 1\right) \left(y - 1\right) \left(z - 1\right) \left(z + 1\right)}{4} - \frac{f_{17} x y \left(x + 1\right) \left(y - 1\right) \left(z - 1\right) \left(z + 1\right)}{4} - \frac{f_{18} x y \left(x + 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{4} - \frac{f_{19} x y \left(x - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{4} + \frac{f_{2} x y z \left(x + 1\right) \left(y + 1\right) \left(z - 1\right)}{8} + \frac{f_{20} x \left(x - 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{2} + \frac{f_{21} x \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{2} + \frac{f_{22} y \left(x - 1\right) \left(x + 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right)}{2} + \frac{f_{23} z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right)}{2} + \frac{f_{24} y \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(z - 1\right) \left(z + 1\right)}{2} + \frac{f_{25} z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z + 1\right)}{2} - f_{26} \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right) \left(z + 1\right) + \frac{f_{3} x y z \left(x - 1\right) \left(y + 1\right) \left(z - 1\right)}{8} + \frac{f_{4} x y z \left(x - 1\right) \left(y - 1\right) \left(z + 1\right)}{8} + \frac{f_{5} x y z \left(x + 1\right) \left(y - 1\right) \left(z + 1\right)}{8} + \frac{f_{6} x y z \left(x + 1\right) \left(y + 1\right) \left(z + 1\right)}{8} + \frac{f_{7} x y z \left(x - 1\right) \left(y + 1\right) \left(z + 1\right)}{8} - \frac{f_{8} y z \left(x - 1\right) \left(x + 1\right) \left(y - 1\right) \left(z - 1\right)}{4} - \frac{f_{9} x z \left(x + 1\right) \left(y - 1\right) \left(y + 1\right) \left(z - 1\right)}{4}$



As an extra aside, if we don't want to assume the node positions: 


```python
# corresponding quadratic lagrange poly 
shape_funcs = []
vtk_node_num = []
crd=[-1,0,1]
for z_i in range(3):
    for y_i in range(3):
        for x_i in range(3):        
            LP1 = LagrangPoly(x,2,x_i)
            LP2 = LagrangPoly(y,2,y_i)
            LP3 = LagrangPoly(z,2,z_i)
            
            # find the VTK node number             
            indx = np.where((all_coords[:,0]==crd[x_i]) & (all_coords[:,1]==crd[y_i]) & (all_coords[:,2]==crd[z_i]))[0][0]                                   
            vtk_node_num.append(indx)
            shape_funcs.append(sy.simplify(LP1 * LP2 * LP3))
            
```


```python
vals = [sy.symbols('f'+str(i)) for i in vtk_node_num]
sample_expression = sum([vals[i]*shape_funcs[i] for i in range(27)])
sample_expression
```




$\displaystyle \frac{f_{0} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right)^{3}} + \frac{f_{1} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)} - \frac{f_{10} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{2}} - \frac{f_{11} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)} - \frac{f_{12} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{2}} - \frac{f_{13} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - z\right) \left(x_{2} - y\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{3}} - \frac{f_{14} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{1} - z\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{3}} - \frac{f_{15} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - z\right) \left(x_{2} - y\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{2}} - \frac{f_{16} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)} - \frac{f_{17} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{2}} - \frac{f_{18} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{3}} - \frac{f_{19} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{2}} + \frac{f_{2} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - y\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)^{2}} + \frac{f_{20} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right) \left(x_{1} - x_{2}\right)^{2}} + \frac{f_{21} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right) \left(x_{1} - x_{2}\right)^{3}} + \frac{f_{22} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right) \left(x_{1} - x_{2}\right)^{3}} + \frac{f_{23} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right) \left(x_{1} - x_{2}\right)^{2}} + \frac{f_{24} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right) \left(x_{1} - x_{2}\right)^{2}} + \frac{f_{25} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - z\right) \left(x_{2} - y\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right) \left(x_{1} - x_{2}\right)^{3}} - \frac{f_{26} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{1} - x_{2}\right)^{3}} + \frac{f_{3} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)} + \frac{f_{4} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)} + \frac{f_{5} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)^{2}} + \frac{f_{6} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{1} - z\right)}{\left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)^{3}} + \frac{f_{7} \left(x - x_{1}\right) \left(x - x_{2}\right) \left(x_{0} - y\right) \left(x_{0} - z\right) \left(x_{1} - y\right) \left(x_{1} - z\right)}{\left(x_{0} - x_{1}\right) \left(x_{0} - x_{2}\right)^{3} \left(x_{1} - x_{2}\right)^{2}} - \frac{f_{8} \left(x - x_{0}\right) \left(x - x_{2}\right) \left(x_{1} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{3} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)} - \frac{f_{9} \left(x - x_{0}\right) \left(x - x_{1}\right) \left(x_{0} - y\right) \left(x_{1} - z\right) \left(x_{2} - y\right) \left(x_{2} - z\right)}{\left(x_{0} - x_{1}\right)^{2} \left(x_{0} - x_{2}\right)^{2} \left(x_{1} - x_{2}\right)^{2}}$




```python

```

cf601dbcc5b26ca79fe9af922e71a1fc34950126
ipynb
Jupyter Notebook
ASPECT_VTK_quad_hex_mapping.ipynb
chrishavlin/unstructured_vtu
b2b2babee2537cf0b633b2f95167f5ae8e78c4ea
BSD-3-Clause
2020-08-05T19:37:01.000Z
2020-08-05T19:37:01.000Z
ASPECT_VTK_quad_hex_mapping.ipynb
chrishavlin/unstructured_vtu
b2b2babee2537cf0b633b2f95167f5ae8e78c4ea
BSD-3-Clause
2020-08-15T15:37:59.000Z
2020-08-15T15:37:59.000Z
ASPECT_VTK_quad_hex_mapping.ipynb
chrishavlin/unstructured_vtu
b2b2babee2537cf0b633b2f95167f5ae8e78c4ea
BSD-3-Clause
2020-06-23T14:35:12.000Z
2020-06-23T14:35:12.000Z
Qwen/Qwen-72B
1. YES
2. YES