ShapeFunctions

Create Lagrange finite element shape functions.

• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 

==================== SEG2 ====================
N1 = lambda r : -(r - 1)/2
N2 = lambda r : (r + 1)/2

dN1 = [lambda r : -1/2]
dN2 = [lambda r : 1/2]

ddN1 = [lambda r : 0]
ddN2 = [lambda r : 0]

dddN1 = [lambda r : 0]
dddN2 = [lambda r : 0]

ddddN1 = [lambda r : 0]
ddddN2 = [lambda r : 0]



============= SEG2 EulerBernoulli =============
N1 = lambda r : (r - 1)**2*(r + 2)/4
N2 = lambda r : (r - 1)**2*(r + 1)/8
N3 = lambda r : -(r - 2)*(r + 1)**2/4
N4 = lambda r : (r - 1)*(r + 1)**2/8

dN1 = [lambda r : (r - 1)**2/4 + (r + 2)*(2*r - 2)/4]
dN2 = [lambda r : (r - 1)**2/8 + (r + 1)*(2*r - 2)/8]
dN3 = [lambda r : -(r - 2)*(2*r + 2)/4 - (r + 1)**2/4]
dN4 = [lambda r : (r - 1)*(2*r + 2)/8 + (r + 1)**2/8]

ddN1 = [lambda r : 3*r/2]
ddN2 = [lambda r : 3*r/4 - 1/4]
ddN3 = [lambda r : -3*r/2]
ddN4 = [lambda r : 3*r/4 + 1/4]

dddN1 = [lambda r : 3/2]
dddN2 = [lambda r : 3/4]
dddN3 = [lambda r : -3/2]
dddN4 = [lambda r : 3/4]

ddddN1 = [lambda r : 0]
ddddN2 = [lambda r : 0]
ddddN3 = [lambda r : 0]
ddddN4 = [lambda r : 0]



==================== SEG3 ====================
N1 = lambda r : r*(r - 1)/2
N2 = lambda r : r*(r + 1)/2
N3 = lambda r : -(r - 1)*(r + 1)

dN1 = [lambda r : r - 1/2]
dN2 = [lambda r : r + 1/2]
dN3 = [lambda r : -2*r]

ddN1 = [lambda r : 1]
ddN2 = [lambda r : 1]
ddN3 = [lambda r : -2]

dddN1 = [lambda r : 0]
dddN2 = [lambda r : 0]
dddN3 = [lambda r : 0]

ddddN1 = [lambda r : 0]
ddddN2 = [lambda r : 0]
ddddN3 = [lambda r : 0]



============= SEG3 EulerBernoulli =============
N1 = lambda r : r**2*(r - 1)**2*(3*r + 4)/4
N2 = lambda r : r**2*(r - 1)**2*(r + 1)/8
N3 = lambda r : -r**2*(r + 1)**2*(3*r - 4)/4
N4 = lambda r : r**2*(r - 1)*(r + 1)**2/8
N5 = lambda r : (r - 1)**2*(r + 1)**2
N6 = lambda r : r*(r - 1)**2*(r + 1)**2/2

dN1 = [lambda r : 3*r**2*(r - 1)**2/4 + r**2*(2*r - 2)*(3*r + 4)/4 + r*(r - 1)**2*(3*r + 4)/2]
dN2 = [lambda r : r**2*(r - 1)**2/8 + r**2*(r + 1)*(2*r - 2)/8 + r*(r - 1)**2*(r + 1)/4]
dN3 = [lambda r : -3*r**2*(r + 1)**2/4 - r**2*(2*r + 2)*(3*r - 4)/4 - r*(r + 1)**2*(3*r - 4)/2]
dN4 = [lambda r : r**2*(r - 1)*(2*r + 2)/8 + r**2*(r + 1)**2/8 + r*(r - 1)*(r + 1)**2/4]
dN5 = [lambda r : (r - 1)**2*(2*r + 2) + (r + 1)**2*(2*r - 2)]
dN6 = [lambda r : r*(r - 1)**2*(2*r + 2)/2 + r*(r + 1)**2*(2*r - 2)/2 + (r - 1)**2*(r + 1)**2/2]

ddN1 = [lambda r : 3*r**2*(2*r - 2)/2 + r**2*(3*r + 4)/2 + 3*r*(r - 1)**2 + r*(2*r - 2)*(3*r + 4) + (r - 1)**2*(3*r + 4)/2]
ddN2 = [lambda r : r**2*(r + 1)/4 + r**2*(2*r - 2)/4 + r*(r - 1)**2/2 + r*(r + 1)*(2*r - 2)/2 + (r - 1)**2*(r + 1)/4]
ddN3 = [lambda r : -3*r**2*(2*r + 2)/2 - r**2*(3*r - 4)/2 - 3*r*(r + 1)**2 - r*(2*r + 2)*(3*r - 4) - (r + 1)**2*(3*r - 4)/2]
ddN4 = [lambda r : r**2*(r - 1)/4 + r**2*(2*r + 2)/4 + r*(r - 1)*(2*r + 2)/2 + r*(r + 1)**2/2 + (r - 1)*(r + 1)**2/4]
ddN5 = [lambda r : 2*(r - 1)**2 + 2*(r + 1)**2 + 2*(2*r - 2)*(2*r + 2)]
ddN6 = [lambda r : r*(r - 1)**2 + r*(r + 1)**2 + r*(2*r - 2)*(2*r + 2) + (r - 1)**2*(2*r + 2) + (r + 1)**2*(2*r - 2)]

dddN1 = [lambda r : 9*r**2/2 + 9*r*(2*r - 2) + 3*r*(3*r + 4) + 9*(r - 1)**2/2 + 3*(2*r - 2)*(3*r + 4)/2]
dddN2 = [lambda r : 3*r**2/4 + 3*r*(r + 1)/2 + 3*r*(2*r - 2)/2 + 3*(r - 1)**2/4 + 3*(r + 1)*(2*r - 2)/4]
dddN3 = [lambda r : -9*r**2/2 - 9*r*(2*r + 2) - 3*r*(3*r - 4) - 9*(r + 1)**2/2 - 3*(2*r + 2)*(3*r - 4)/2]
dddN4 = [lambda r : 3*r**2/4 + 3*r*(r - 1)/2 + 3*r*(2*r + 2)/2 + 3*(r - 1)*(2*r + 2)/4 + 3*(r + 1)**2/4]
dddN5 = [lambda r : 24*r]
dddN6 = [lambda r : 3*r*(2*r - 2) + 3*r*(2*r + 2) + 3*(r - 1)**2 + 3*(r + 1)**2 + 3*(2*r - 2)*(2*r + 2)]

ddddN1 = [lambda r : 90*r - 12]
ddddN2 = [lambda r : 15*r - 3]
ddddN3 = [lambda r : -90*r - 12]
ddddN4 = [lambda r : 15*r + 3]
ddddN5 = [lambda r : 24]
ddddN6 = [lambda r : 60*r]



==================== SEG4 ====================
N1 = lambda r : -9*r**3/16 + 9*r**2/16 + r/16 - 1/16
N2 = lambda r : 9*r**3/16 + 9*r**2/16 - r/16 - 1/16
N3 = lambda r : 27*r**3/16 - 9*r**2/16 - 27*r/16 + 9/16
N4 = lambda r : -27*r**3/16 - 9*r**2/16 + 27*r/16 + 9/16

dN1 = [lambda r : -27*r**2/16 + 9*r/8 + 1/16]
dN2 = [lambda r : 27*r**2/16 + 9*r/8 - 1/16]
dN3 = [lambda r : 81*r**2/16 - 9*r/8 - 27/16]
dN4 = [lambda r : -81*r**2/16 - 9*r/8 + 27/16]

ddN1 = [lambda r : 9/8 - 27*r/8]
ddN2 = [lambda r : 27*r/8 + 9/8]
ddN3 = [lambda r : 81*r/8 - 9/8]
ddN4 = [lambda r : -81*r/8 - 9/8]

dddN1 = [lambda r : -27/8]
dddN2 = [lambda r : 27/8]
dddN3 = [lambda r : 81/8]
dddN4 = [lambda r : -81/8]

ddddN1 = [lambda r : 0]
ddddN2 = [lambda r : 0]
ddddN3 = [lambda r : 0]
ddddN4 = [lambda r : 0]



============= SEG4 EulerBernoulli =============
N1 = lambda r : 891*r**7/512 - 729*r**6/512 - 1413*r**5/512 + 1215*r**4/512 + 137207031249999*r**3/250000000000000 - 237304687500001*r**2/500000000000000 - 292968749999991*r/10000000000000000 + 253906250000009/10000000000000000
N2 = lambda r : 81*r**7/512 - 81*r**6/512 - 99*r**5/512 + 99*r**4/512 + 92773437499999*r**3/2500000000000000 - 185546875000001*r**2/5000000000000000 - r/512 + 1/512
N3 = lambda r : -891*r**7/512 - 729*r**6/512 + 275976562499999*r**5/100000000000000 + 1215*r**4/512 - 109765624999999*r**3/200000000000000 - 237304687499999*r**2/500000000000000 + 146484374999991*r/5000000000000000 + 63476562499999/2500000000000000
N4 = lambda r : 81*r**7/512 + 81*r**6/512 - 99*r**5/512 - 99*r**4/512 + 92773437499999*r**3/2500000000000000 + 371093749999999*r**2/10000000000000000 - r/512 - 1/512
N5 = lambda r : 2187*r**7/512 + 729*r**6/512 - 5589*r**5/512 - 1215*r**4/512 + 450878906249999*r**3/50000000000000 + 94921874999999*r**2/200000000000000 - 1215*r/512 + 474609375000003/1000000000000000
N6 = lambda r : 729*r**7/512 - 243*r**6/512 - 1539*r**5/512 + 513*r**4/512 + 891*r**3/512 - 290039062500001*r**2/500000000000000 - 158203124999999*r/1000000000000000 + 527343750000009/10000000000000000
N7 = lambda r : -2187*r**7/512 + 729*r**6/512 + 5589*r**5/512 - 1215*r**4/512 - 450878906249999*r**3/50000000000000 + 474609375000007*r**2/1000000000000000 + 1215*r/512 + 118652343749999/250000000000000
N8 = lambda r : 729*r**7/512 + 243*r**6/512 - 300585937499999*r**5/100000000000000 - 513*r**4/512 + 174023437499999*r**3/100000000000000 + 290039062499999*r**2/500000000000000 - 79101562499999*r/500000000000000 - 527343749999993/10000000000000000

dN1 = [lambda r : 6237*r**6/512 - 2187*r**5/256 - 7065*r**4/512 + 1215*r**3/128 + 411621093749997*r**2/250000000000000 - 237304687500001*r/250000000000000 - 292968749999991/10000000000000000]
dN2 = [lambda r : 567*r**6/512 - 243*r**5/256 - 495*r**4/512 + 99*r**3/128 + 278320312499997*r**2/2500000000000000 - 185546875000001*r/2500000000000000 - 1/512]
dN3 = [lambda r : -6237*r**6/512 - 2187*r**5/256 + 275976562499999*r**4/20000000000000 + 1215*r**3/128 - 329296874999997*r**2/200000000000000 - 237304687499999*r/250000000000000 + 146484374999991/5000000000000000]
dN4 = [lambda r : 567*r**6/512 + 243*r**5/256 - 495*r**4/512 - 99*r**3/128 + 278320312499997*r**2/2500000000000000 + 371093749999999*r/5000000000000000 - 1/512]
dN5 = [lambda r : 15309*r**6/512 + 2187*r**5/256 - 27945*r**4/512 - 1215*r**3/128 + 1352636718749997*r**2/50000000000000 + 94921874999999*r/100000000000000 - 1215/512]
dN6 = [lambda r : 5103*r**6/512 - 729*r**5/256 - 7695*r**4/512 + 513*r**3/128 + 2673*r**2/512 - 290039062500001*r/250000000000000 - 158203124999999/1000000000000000]
dN7 = [lambda r : -15309*r**6/512 + 2187*r**5/256 + 27945*r**4/512 - 1215*r**3/128 - 1352636718749997*r**2/50000000000000 + 474609375000007*r/500000000000000 + 1215/512]
dN8 = [lambda r : 5103*r**6/512 + 729*r**5/256 - 300585937499999*r**4/20000000000000 - 513*r**3/128 + 522070312499997*r**2/100000000000000 + 290039062499999*r/250000000000000 - 79101562499999/500000000000000]

ddN1 = [lambda r : 18711*r**5/256 - 10935*r**4/256 - 7065*r**3/128 + 3645*r**2/128 + 411621093749997*r/125000000000000 - 237304687500001/250000000000000]
ddN2 = [lambda r : 1701*r**5/256 - 1215*r**4/256 - 495*r**3/128 + 297*r**2/128 + 278320312499997*r/1250000000000000 - 185546875000001/2500000000000000]
ddN3 = [lambda r : -18711*r**5/256 - 10935*r**4/256 + 275976562499999*r**3/5000000000000 + 3645*r**2/128 - 329296874999997*r/100000000000000 - 237304687499999/250000000000000]
ddN4 = [lambda r : 1701*r**5/256 + 1215*r**4/256 - 495*r**3/128 - 297*r**2/128 + 278320312499997*r/1250000000000000 + 371093749999999/5000000000000000]
ddN5 = [lambda r : 45927*r**5/256 + 10935*r**4/256 - 27945*r**3/128 - 3645*r**2/128 + 1352636718749997*r/25000000000000 + 94921874999999/100000000000000]
ddN6 = [lambda r : 15309*r**5/256 - 3645*r**4/256 - 7695*r**3/128 + 1539*r**2/128 + 2673*r/256 - 290039062500001/250000000000000]
ddN7 = [lambda r : -45927*r**5/256 + 10935*r**4/256 + 27945*r**3/128 - 3645*r**2/128 - 1352636718749997*r/25000000000000 + 474609375000007/500000000000000]
ddN8 = [lambda r : 15309*r**5/256 + 3645*r**4/256 - 300585937499999*r**3/5000000000000 - 1539*r**2/128 + 522070312499997*r/50000000000000 + 290039062499999/250000000000000]

dddN1 = [lambda r : 93555*r**4/256 - 10935*r**3/64 - 21195*r**2/128 + 3645*r/64 + 411621093749997/125000000000000]
dddN2 = [lambda r : 8505*r**4/256 - 1215*r**3/64 - 1485*r**2/128 + 297*r/64 + 278320312499997/1250000000000000]
dddN3 = [lambda r : -93555*r**4/256 - 10935*r**3/64 + 827929687499997*r**2/5000000000000 + 3645*r/64 - 329296874999997/100000000000000]
dddN4 = [lambda r : 8505*r**4/256 + 1215*r**3/64 - 1485*r**2/128 - 297*r/64 + 278320312499997/1250000000000000]
dddN5 = [lambda r : 229635*r**4/256 + 10935*r**3/64 - 83835*r**2/128 - 3645*r/64 + 1352636718749997/25000000000000]
dddN6 = [lambda r : 76545*r**4/256 - 3645*r**3/64 - 23085*r**2/128 + 1539*r/64 + 2673/256]
dddN7 = [lambda r : -229635*r**4/256 + 10935*r**3/64 + 83835*r**2/128 - 3645*r/64 - 1352636718749997/25000000000000]
dddN8 = [lambda r : 76545*r**4/256 + 3645*r**3/64 - 901757812499997*r**2/5000000000000 - 1539*r/64 + 522070312499997/50000000000000]

ddddN1 = [lambda r : 93555*r**3/64 - 32805*r**2/64 - 21195*r/64 + 3645/64]
ddddN2 = [lambda r : 8505*r**3/64 - 3645*r**2/64 - 1485*r/64 + 297/64]
ddddN3 = [lambda r : -93555*r**3/64 - 32805*r**2/64 + 827929687499997*r/2500000000000 + 3645/64]
ddddN4 = [lambda r : 8505*r**3/64 + 3645*r**2/64 - 1485*r/64 - 297/64]
ddddN5 = [lambda r : 229635*r**3/64 + 32805*r**2/64 - 83835*r/64 - 3645/64]
ddddN6 = [lambda r : 76545*r**3/64 - 10935*r**2/64 - 23085*r/64 + 1539/64]
ddddN7 = [lambda r : -229635*r**3/64 + 32805*r**2/64 + 83835*r/64 - 3645/64]
ddddN8 = [lambda r : 76545*r**3/64 + 10935*r**2/64 - 901757812499997*r/2500000000000 - 1539/64]



==================== SEG5 ====================
N1 = lambda r : 2*r**4/3 - 2*r**3/3 - r**2/6 + r/6
N2 = lambda r : 2*r**4/3 + 2*r**3/3 - r**2/6 - r/6
N3 = lambda r : -8*r**4/3 + 4*r**3/3 + 8*r**2/3 - 4*r/3
N4 = lambda r : 4*r**4 - 5*r**2 + 1
N5 = lambda r : -8*r**4/3 - 4*r**3/3 + 8*r**2/3 + 4*r/3

dN1 = [lambda r : 8*r**3/3 - 2*r**2 - r/3 + 1/6]
dN2 = [lambda r : 8*r**3/3 + 2*r**2 - r/3 - 1/6]
dN3 = [lambda r : -32*r**3/3 + 4*r**2 + 16*r/3 - 4/3]
dN4 = [lambda r : 16*r**3 - 10*r]
dN5 = [lambda r : -32*r**3/3 - 4*r**2 + 16*r/3 + 4/3]

ddN1 = [lambda r : 8*r**2 - 4*r - 1/3]
ddN2 = [lambda r : 8*r**2 + 4*r - 1/3]
ddN3 = [lambda r : -32*r**2 + 8*r + 16/3]
ddN4 = [lambda r : 48*r**2 - 10]
ddN5 = [lambda r : -32*r**2 - 8*r + 16/3]

dddN1 = [lambda r : 16*r - 4]
dddN2 = [lambda r : 16*r + 4]
dddN3 = [lambda r : 8 - 64*r]
dddN4 = [lambda r : 96*r]
dddN5 = [lambda r : -64*r - 8]

ddddN1 = [lambda r : 16]
ddddN2 = [lambda r : 16]
ddddN3 = [lambda r : -64]
ddddN4 = [lambda r : 96]
ddddN5 = [lambda r : -64]



============= SEG5 EulerBernoulli =============
N1 = lambda r : 100*r**9/27 - 13037037037037*r**8/4000000000000 - 58*r**7/9 + 23111111111111*r**6/4000000000000 + 91*r**5/36 - 9111111111111*r**4/4000000000000 - 31*r**3/108 + 7*r**2/27
N2 = lambda r : 2*r**9/9 - 2*r**8/9 - r**7/3 + r**6/3 + r**5/8 - r**4/8 - r**3/72 + r**2/72
N3 = lambda r : -100*r**9/27 - 325925925925927*r**8/100000000000000 + 58*r**7/9 + 577777777777781*r**6/100000000000000 - 91*r**5/36 - 11388888888889*r**4/5000000000000 + 31*r**3/108 + 7*r**2/27
N4 = lambda r : 2*r**9/9 + 2*r**8/9 - r**7/3 - r**6/3 + r**5/8 + r**4/8 - r**3/72 - r**2/72
N5 = lambda r : 640*r**9/27 - 474074074074067*r**8/100000000000000 - 544*r**7/9 + 7111111111111*r**6/500000000000 + 448*r**5/9 - 142222222222221*r**4/10000000000000 - 352*r**3/27 + 474074074074069*r**2/100000000000000
N6 = lambda r : 32*r**9/9 - 16*r**8/9 - 8*r**7 + 4*r**6 + 16*r**5/3 - 8*r**4/3 - 8*r**3/9 + 4*r**2/9
N7 = lambda r : 16*r**8 - 40*r**6 + 33*r**4 - 10*r**2 + 1
N8 = lambda r : 8*r**9 - 20*r**7 + 33*r**5/2 - 5*r**3 + r/2
N9 = lambda r : -640*r**9/27 - 474074074074081*r**8/100000000000000 + 544*r**7/9 + 8888888888889*r**6/625000000000 - 448*r**5/9 - 8888888888889*r**4/625000000000 + 352*r**3/27 + 237037037037039*r**2/50000000000000
N10 = lambda r : 32*r**9/9 + 16*r**8/9 - 8*r**7 - 4*r**6 + 16*r**5/3 + 8*r**4/3 - 8*r**3/9 - 4*r**2/9

dN1 = [lambda r : 100*r**8/3 - 13037037037037*r**7/500000000000 - 406*r**6/9 + 69333333333333*r**5/2000000000000 + 455*r**4/36 - 9111111111111*r**3/1000000000000 - 31*r**2/36 + 14*r/27]
dN2 = [lambda r : 2*r**8 - 16*r**7/9 - 7*r**6/3 + 2*r**5 + 5*r**4/8 - r**3/2 - r**2/24 + r/36]
dN3 = [lambda r : -100*r**8/3 - 325925925925927*r**7/12500000000000 + 406*r**6/9 + 1733333333333343*r**5/50000000000000 - 455*r**4/36 - 11388888888889*r**3/1250000000000 + 31*r**2/36 + 14*r/27]
dN4 = [lambda r : 2*r**8 + 16*r**7/9 - 7*r**6/3 - 2*r**5 + 5*r**4/8 + r**3/2 - r**2/24 - r/36]
dN5 = [lambda r : 640*r**8/3 - 474074074074067*r**7/12500000000000 - 3808*r**6/9 + 21333333333333*r**5/250000000000 + 2240*r**4/9 - 142222222222221*r**3/2500000000000 - 352*r**2/9 + 474074074074069*r/50000000000000]
dN6 = [lambda r : 32*r**8 - 128*r**7/9 - 56*r**6 + 24*r**5 + 80*r**4/3 - 32*r**3/3 - 8*r**2/3 + 8*r/9]
dN7 = [lambda r : 128*r**7 - 240*r**5 + 132*r**3 - 20*r]
dN8 = [lambda r : 72*r**8 - 140*r**6 + 165*r**4/2 - 15*r**2 + 1/2]
dN9 = [lambda r : -640*r**8/3 - 474074074074081*r**7/12500000000000 + 3808*r**6/9 + 26666666666667*r**5/312500000000 - 2240*r**4/9 - 8888888888889*r**3/156250000000 + 352*r**2/9 + 237037037037039*r/25000000000000]
dN10 = [lambda r : 32*r**8 + 128*r**7/9 - 56*r**6 - 24*r**5 + 80*r**4/3 + 32*r**3/3 - 8*r**2/3 - 8*r/9]

ddN1 = [lambda r : 800*r**7/3 - 91259259259259*r**6/500000000000 - 812*r**5/3 + 69333333333333*r**4/400000000000 + 455*r**3/9 - 27333333333333*r**2/1000000000000 - 31*r/18 + 14/27]
ddN2 = [lambda r : 16*r**7 - 112*r**6/9 - 14*r**5 + 10*r**4 + 5*r**3/2 - 3*r**2/2 - r/12 + 1/36]
ddN3 = [lambda r : -800*r**7/3 - 2281481481481489*r**6/12500000000000 + 812*r**5/3 + 1733333333333343*r**4/10000000000000 - 455*r**3/9 - 34166666666667*r**2/1250000000000 + 31*r/18 + 14/27]
ddN4 = [lambda r : 16*r**7 + 112*r**6/9 - 14*r**5 - 10*r**4 + 5*r**3/2 + 3*r**2/2 - r/12 - 1/36]
ddN5 = [lambda r : 5120*r**7/3 - 3318518518518469*r**6/12500000000000 - 7616*r**5/3 + 21333333333333*r**4/50000000000 + 8960*r**3/9 - 426666666666663*r**2/2500000000000 - 704*r/9 + 474074074074069/50000000000000]
ddN6 = [lambda r : 256*r**7 - 896*r**6/9 - 336*r**5 + 120*r**4 + 320*r**3/3 - 32*r**2 - 16*r/3 + 8/9]
ddN7 = [lambda r : 896*r**6 - 1200*r**4 + 396*r**2 - 20]
ddN8 = [lambda r : 576*r**7 - 840*r**5 + 330*r**3 - 30*r]
ddN9 = [lambda r : -5120*r**7/3 - 3318518518518567*r**6/12500000000000 + 7616*r**5/3 + 26666666666667*r**4/62500000000 - 8960*r**3/9 - 26666666666667*r**2/156250000000 + 704*r/9 + 237037037037039/25000000000000]
ddN10 = [lambda r : 256*r**7 + 896*r**6/9 - 336*r**5 - 120*r**4 + 320*r**3/3 + 32*r**2 - 16*r/3 - 8/9]

dddN1 = [lambda r : 5600*r**6/3 - 273777777777777*r**5/250000000000 - 4060*r**4/3 + 69333333333333*r**3/100000000000 + 455*r**2/3 - 27333333333333*r/500000000000 - 31/18]
dddN2 = [lambda r : 112*r**6 - 224*r**5/3 - 70*r**4 + 40*r**3 + 15*r**2/2 - 3*r - 1/12]
dddN3 = [lambda r : -5600*r**6/3 - 6844444444444467*r**5/6250000000000 + 4060*r**4/3 + 1733333333333343*r**3/2500000000000 - 455*r**2/3 - 34166666666667*r/625000000000 + 31/18]
dddN4 = [lambda r : 112*r**6 + 224*r**5/3 - 70*r**4 - 40*r**3 + 15*r**2/2 + 3*r - 1/12]
dddN5 = [lambda r : 35840*r**6/3 - 9955555555555407*r**5/6250000000000 - 38080*r**4/3 + 21333333333333*r**3/12500000000 + 8960*r**2/3 - 426666666666663*r/1250000000000 - 704/9]
dddN6 = [lambda r : 1792*r**6 - 1792*r**5/3 - 1680*r**4 + 480*r**3 + 320*r**2 - 64*r - 16/3]
dddN7 = [lambda r : 5376*r**5 - 4800*r**3 + 792*r]
dddN8 = [lambda r : 4032*r**6 - 4200*r**4 + 990*r**2 - 30]
dddN9 = [lambda r : -35840*r**6/3 - 9955555555555701*r**5/6250000000000 + 38080*r**4/3 + 26666666666667*r**3/15625000000 - 8960*r**2/3 - 26666666666667*r/78125000000 + 704/9]
dddN10 = [lambda r : 1792*r**6 + 1792*r**5/3 - 1680*r**4 - 480*r**3 + 320*r**2 + 64*r - 16/3]

ddddN1 = [lambda r : 11200*r**5 - 273777777777777*r**4/50000000000 - 16240*r**3/3 + 207999999999999*r**2/100000000000 + 910*r/3 - 27333333333333/500000000000]
ddddN2 = [lambda r : 672*r**5 - 1120*r**4/3 - 280*r**3 + 120*r**2 + 15*r - 3]
ddddN3 = [lambda r : -11200*r**5 - 6844444444444467*r**4/1250000000000 + 16240*r**3/3 + 5200000000000029*r**2/2500000000000 - 910*r/3 - 34166666666667/625000000000]
ddddN4 = [lambda r : 672*r**5 + 1120*r**4/3 - 280*r**3 - 120*r**2 + 15*r + 3]
ddddN5 = [lambda r : 71680*r**5 - 9955555555555407*r**4/1250000000000 - 152320*r**3/3 + 63999999999999*r**2/12500000000 + 17920*r/3 - 426666666666663/1250000000000]
ddddN6 = [lambda r : 10752*r**5 - 8960*r**4/3 - 6720*r**3 + 1440*r**2 + 640*r - 64]
ddddN7 = [lambda r : 26880*r**4 - 14400*r**2 + 792]
ddddN8 = [lambda r : 24192*r**5 - 16800*r**3 + 1980*r]
ddddN9 = [lambda r : -71680*r**5 - 9955555555555701*r**4/1250000000000 + 152320*r**3/3 + 80000000000001*r**2/15625000000 - 17920*r/3 - 26666666666667/78125000000]
ddddN10 = [lambda r : 10752*r**5 + 8960*r**4/3 - 6720*r**3 - 1440*r**2 + 640*r + 64]



==================== TRI3 ====================
N1 = lambda r, s : -r - s + 1
N2 = lambda r, s : r
N3 = lambda r, s : s

dN1 = [lambda r, s : -1, lambda r, s : -1]
dN2 = [lambda r, s : 1, lambda r, s : 0]
dN3 = [lambda r, s : 0, lambda r, s : 1]

ddN1 = [lambda r, s : 0, lambda r, s : 0]
ddN2 = [lambda r, s : 0, lambda r, s : 0]
ddN3 = [lambda r, s : 0, lambda r, s : 0]

dddN1 = [lambda r, s : 0, lambda r, s : 0]
dddN2 = [lambda r, s : 0, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 0]

ddddN1 = [lambda r, s : 0, lambda r, s : 0]
ddddN2 = [lambda r, s : 0, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 0]



==================== TRI6 ====================
N1 = lambda r, s : (r + s - 1)*(2*r + 2*s - 1)
N2 = lambda r, s : r*(2*r - 1)
N3 = lambda r, s : s*(2*s - 1)
N4 = lambda r, s : -4*r*(r + s - 1)
N5 = lambda r, s : 4*r*s
N6 = lambda r, s : -4*s*(r + s - 1)

dN1 = [lambda r, s : 4*r + 4*s - 3, lambda r, s : 4*r + 4*s - 3]
dN2 = [lambda r, s : 4*r - 1, lambda r, s : 0]
dN3 = [lambda r, s : 0, lambda r, s : 4*s - 1]
dN4 = [lambda r, s : -8*r - 4*s + 4, lambda r, s : -4*r]
dN5 = [lambda r, s : 4*s, lambda r, s : 4*r]
dN6 = [lambda r, s : -4*s, lambda r, s : -4*r - 8*s + 4]

ddN1 = [lambda r, s : 4, lambda r, s : 4]
ddN2 = [lambda r, s : 4, lambda r, s : 0]
ddN3 = [lambda r, s : 0, lambda r, s : 4]
ddN4 = [lambda r, s : -8, lambda r, s : 0]
ddN5 = [lambda r, s : 0, lambda r, s : 0]
ddN6 = [lambda r, s : 0, lambda r, s : -8]

dddN1 = [lambda r, s : 0, lambda r, s : 0]
dddN2 = [lambda r, s : 0, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 0]
dddN4 = [lambda r, s : 0, lambda r, s : 0]
dddN5 = [lambda r, s : 0, lambda r, s : 0]
dddN6 = [lambda r, s : 0, lambda r, s : 0]

ddddN1 = [lambda r, s : 0, lambda r, s : 0]
ddddN2 = [lambda r, s : 0, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 0]
ddddN4 = [lambda r, s : 0, lambda r, s : 0]
ddddN5 = [lambda r, s : 0, lambda r, s : 0]
ddddN6 = [lambda r, s : 0, lambda r, s : 0]



==================== TRI10 ====================
N1 = lambda r, s : -9*r**3/2 - 27*r**2*s/2 + 9*r**2 - 27*r*s**2/2 + 18*r*s - 11*r/2 - 9*s**3/2 + 9*s**2 - 11*s/2 + 1
N2 = lambda r, s : 9*r**3/2 - 9*r**2/2 + r
N3 = lambda r, s : 9*s**3/2 - 9*s**2/2 + s
N4 = lambda r, s : 27*r**3/2 + 27*r**2*s - 45*r**2/2 + 27*r*s**2/2 - 45*r*s/2 + 9*r
N5 = lambda r, s : -27*r**3/2 - 27*r**2*s/2 + 18*r**2 + 9*r*s/2 - 9*r/2
N6 = lambda r, s : 27*r**2*s/2 - 9*r*s/2
N7 = lambda r, s : 27*r*s**2/2 - 9*r*s/2
N8 = lambda r, s : -27*r*s**2/2 + 9*r*s/2 - 27*s**3/2 + 18*s**2 - 9*s/2
N9 = lambda r, s : 27*r**2*s/2 + 27*r*s**2 - 45*r*s/2 + 27*s**3/2 - 45*s**2/2 + 9*s
N10 = lambda r, s : -27*r**2*s - 27*r*s**2 + 27*r*s

dN1 = [lambda r, s : -27*r**2/2 - 27*r*s + 18*r - 27*s**2/2 + 18*s - 11/2, lambda r, s : -27*r**2/2 - 27*r*s + 18*r - 27*s**2/2 + 18*s - 11/2]
dN2 = [lambda r, s : 27*r**2/2 - 9*r + 1, lambda r, s : 0]
dN3 = [lambda r, s : 0, lambda r, s : 27*s**2/2 - 9*s + 1]
dN4 = [lambda r, s : 81*r**2/2 + 54*r*s - 45*r + 27*s**2/2 - 45*s/2 + 9, lambda r, s : 27*r**2 + 27*r*s - 45*r/2]
dN5 = [lambda r, s : -81*r**2/2 - 27*r*s + 36*r + 9*s/2 - 9/2, lambda r, s : -27*r**2/2 + 9*r/2]
dN6 = [lambda r, s : 27*r*s - 9*s/2, lambda r, s : 27*r**2/2 - 9*r/2]
dN7 = [lambda r, s : 27*s**2/2 - 9*s/2, lambda r, s : 27*r*s - 9*r/2]
dN8 = [lambda r, s : -27*s**2/2 + 9*s/2, lambda r, s : -27*r*s + 9*r/2 - 81*s**2/2 + 36*s - 9/2]
dN9 = [lambda r, s : 27*r*s + 27*s**2 - 45*s/2, lambda r, s : 27*r**2/2 + 54*r*s - 45*r/2 + 81*s**2/2 - 45*s + 9]
dN10 = [lambda r, s : -54*r*s - 27*s**2 + 27*s, lambda r, s : -27*r**2 - 54*r*s + 27*r]

ddN1 = [lambda r, s : -27*r - 27*s + 18, lambda r, s : -27*r - 27*s + 18]
ddN2 = [lambda r, s : 27*r - 9, lambda r, s : 0]
ddN3 = [lambda r, s : 0, lambda r, s : 27*s - 9]
ddN4 = [lambda r, s : 81*r + 54*s - 45, lambda r, s : 27*r]
ddN5 = [lambda r, s : -81*r - 27*s + 36, lambda r, s : 0]
ddN6 = [lambda r, s : 27*s, lambda r, s : 0]
ddN7 = [lambda r, s : 0, lambda r, s : 27*r]
ddN8 = [lambda r, s : 0, lambda r, s : -27*r - 81*s + 36]
ddN9 = [lambda r, s : 27*s, lambda r, s : 54*r + 81*s - 45]
ddN10 = [lambda r, s : -54*s, lambda r, s : -54*r]

dddN1 = [lambda r, s : -27, lambda r, s : -27]
dddN2 = [lambda r, s : 27, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 27]
dddN4 = [lambda r, s : 81, lambda r, s : 0]
dddN5 = [lambda r, s : -81, lambda r, s : 0]
dddN6 = [lambda r, s : 0, lambda r, s : 0]
dddN7 = [lambda r, s : 0, lambda r, s : 0]
dddN8 = [lambda r, s : 0, lambda r, s : -81]
dddN9 = [lambda r, s : 0, lambda r, s : 81]
dddN10 = [lambda r, s : 0, lambda r, s : 0]

ddddN1 = [lambda r, s : 0, lambda r, s : 0]
ddddN2 = [lambda r, s : 0, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 0]
ddddN4 = [lambda r, s : 0, lambda r, s : 0]
ddddN5 = [lambda r, s : 0, lambda r, s : 0]
ddddN6 = [lambda r, s : 0, lambda r, s : 0]
ddddN7 = [lambda r, s : 0, lambda r, s : 0]
ddddN8 = [lambda r, s : 0, lambda r, s : 0]
ddddN9 = [lambda r, s : 0, lambda r, s : 0]
ddddN10 = [lambda r, s : 0, lambda r, s : 0]



==================== TRI15 ====================
N1 = lambda r, s : 32*r**4/3 + 128*r**3*s/3 - 80*r**3/3 + 64*r**2*s**2 - 80*r**2*s + 70*r**2/3 + 128*r*s**3/3 - 80*r*s**2 + 140*r*s/3 - 25*r/3 + 32*s**4/3 - 80*s**3/3 + 70*s**2/3 - 25*s/3 + 1
N2 = lambda r, s : 32*r**4/3 - 16*r**3 + 22*r**2/3 - r
N3 = lambda r, s : 32*s**4/3 - 16*s**3 + 22*s**2/3 - s
N4 = lambda r, s : -128*r**4/3 - 128*r**3*s + 96*r**3 - 128*r**2*s**2 + 192*r**2*s - 208*r**2/3 - 128*r*s**3/3 + 96*r*s**2 - 208*r*s/3 + 16*r
N5 = lambda r, s : 64*r**4 + 128*r**3*s - 128*r**3 + 64*r**2*s**2 - 144*r**2*s + 76*r**2 - 16*r*s**2 + 28*r*s - 12*r
N6 = lambda r, s : -128*r**4/3 - 128*r**3*s/3 + 224*r**3/3 + 32*r**2*s - 112*r**2/3 - 16*r*s/3 + 16*r/3
N7 = lambda r, s : 128*r**3*s/3 - 32*r**2*s + 16*r*s/3
N8 = lambda r, s : 64*r**2*s**2 - 16*r**2*s - 16*r*s**2 + 4*r*s
N9 = lambda r, s : 128*r*s**3/3 - 32*r*s**2 + 16*r*s/3
N10 = lambda r, s : -128*r*s**3/3 + 32*r*s**2 - 16*r*s/3 - 128*s**4/3 + 224*s**3/3 - 112*s**2/3 + 16*s/3
N11 = lambda r, s : 64*r**2*s**2 - 16*r**2*s + 128*r*s**3 - 144*r*s**2 + 28*r*s + 64*s**4 - 128*s**3 + 76*s**2 - 12*s
N12 = lambda r, s : -128*r**3*s/3 - 128*r**2*s**2 + 96*r**2*s - 128*r*s**3 + 192*r*s**2 - 208*r*s/3 - 128*s**4/3 + 96*s**3 - 208*s**2/3 + 16*s
N13 = lambda r, s : 128*r**3*s + 256*r**2*s**2 - 224*r**2*s + 128*r*s**3 - 224*r*s**2 + 96*r*s
N14 = lambda r, s : -128*r**3*s - 128*r**2*s**2 + 160*r**2*s + 32*r*s**2 - 32*r*s
N15 = lambda r, s : -128*r**2*s**2 + 32*r**2*s - 128*r*s**3 + 160*r*s**2 - 32*r*s

dN1 = [lambda r, s : 128*r**3/3 + 128*r**2*s - 80*r**2 + 128*r*s**2 - 160*r*s + 140*r/3 + 128*s**3/3 - 80*s**2 + 140*s/3 - 25/3, lambda r, s : 128*r**3/3 + 128*r**2*s - 80*r**2 + 128*r*s**2 - 160*r*s + 140*r/3 + 128*s**3/3 - 80*s**2 + 140*s/3 - 25/3]
dN2 = [lambda r, s : 128*r**3/3 - 48*r**2 + 44*r/3 - 1, lambda r, s : 0]
dN3 = [lambda r, s : 0, lambda r, s : 128*s**3/3 - 48*s**2 + 44*s/3 - 1]
dN4 = [lambda r, s : -512*r**3/3 - 384*r**2*s + 288*r**2 - 256*r*s**2 + 384*r*s - 416*r/3 - 128*s**3/3 + 96*s**2 - 208*s/3 + 16, lambda r, s : -128*r**3 - 256*r**2*s + 192*r**2 - 128*r*s**2 + 192*r*s - 208*r/3]
dN5 = [lambda r, s : 256*r**3 + 384*r**2*s - 384*r**2 + 128*r*s**2 - 288*r*s + 152*r - 16*s**2 + 28*s - 12, lambda r, s : 128*r**3 + 128*r**2*s - 144*r**2 - 32*r*s + 28*r]
dN6 = [lambda r, s : -512*r**3/3 - 128*r**2*s + 224*r**2 + 64*r*s - 224*r/3 - 16*s/3 + 16/3, lambda r, s : -128*r**3/3 + 32*r**2 - 16*r/3]
dN7 = [lambda r, s : 128*r**2*s - 64*r*s + 16*s/3, lambda r, s : 128*r**3/3 - 32*r**2 + 16*r/3]
dN8 = [lambda r, s : 128*r*s**2 - 32*r*s - 16*s**2 + 4*s, lambda r, s : 128*r**2*s - 16*r**2 - 32*r*s + 4*r]
dN9 = [lambda r, s : 128*s**3/3 - 32*s**2 + 16*s/3, lambda r, s : 128*r*s**2 - 64*r*s + 16*r/3]
dN10 = [lambda r, s : -128*s**3/3 + 32*s**2 - 16*s/3, lambda r, s : -128*r*s**2 + 64*r*s - 16*r/3 - 512*s**3/3 + 224*s**2 - 224*s/3 + 16/3]
dN11 = [lambda r, s : 128*r*s**2 - 32*r*s + 128*s**3 - 144*s**2 + 28*s, lambda r, s : 128*r**2*s - 16*r**2 + 384*r*s**2 - 288*r*s + 28*r + 256*s**3 - 384*s**2 + 152*s - 12]
dN12 = [lambda r, s : -128*r**2*s - 256*r*s**2 + 192*r*s - 128*s**3 + 192*s**2 - 208*s/3, lambda r, s : -128*r**3/3 - 256*r**2*s + 96*r**2 - 384*r*s**2 + 384*r*s - 208*r/3 - 512*s**3/3 + 288*s**2 - 416*s/3 + 16]
dN13 = [lambda r, s : 384*r**2*s + 512*r*s**2 - 448*r*s + 128*s**3 - 224*s**2 + 96*s, lambda r, s : 128*r**3 + 512*r**2*s - 224*r**2 + 384*r*s**2 - 448*r*s + 96*r]
dN14 = [lambda r, s : -384*r**2*s - 256*r*s**2 + 320*r*s + 32*s**2 - 32*s, lambda r, s : -128*r**3 - 256*r**2*s + 160*r**2 + 64*r*s - 32*r]
dN15 = [lambda r, s : -256*r*s**2 + 64*r*s - 128*s**3 + 160*s**2 - 32*s, lambda r, s : -256*r**2*s + 32*r**2 - 384*r*s**2 + 320*r*s - 32*r]

ddN1 = [lambda r, s : 128*r**2 + 256*r*s - 160*r + 128*s**2 - 160*s + 140/3, lambda r, s : 128*r**2 + 256*r*s - 160*r + 128*s**2 - 160*s + 140/3]
ddN2 = [lambda r, s : 128*r**2 - 96*r + 44/3, lambda r, s : 0]
ddN3 = [lambda r, s : 0, lambda r, s : 128*s**2 - 96*s + 44/3]
ddN4 = [lambda r, s : -512*r**2 - 768*r*s + 576*r - 256*s**2 + 384*s - 416/3, lambda r, s : -256*r**2 - 256*r*s + 192*r]
ddN5 = [lambda r, s : 768*r**2 + 768*r*s - 768*r + 128*s**2 - 288*s + 152, lambda r, s : 128*r**2 - 32*r]
ddN6 = [lambda r, s : -512*r**2 - 256*r*s + 448*r + 64*s - 224/3, lambda r, s : 0]
ddN7 = [lambda r, s : 256*r*s - 64*s, lambda r, s : 0]
ddN8 = [lambda r, s : 128*s**2 - 32*s, lambda r, s : 128*r**2 - 32*r]
ddN9 = [lambda r, s : 0, lambda r, s : 256*r*s - 64*r]
ddN10 = [lambda r, s : 0, lambda r, s : -256*r*s + 64*r - 512*s**2 + 448*s - 224/3]
ddN11 = [lambda r, s : 128*s**2 - 32*s, lambda r, s : 128*r**2 + 768*r*s - 288*r + 768*s**2 - 768*s + 152]
ddN12 = [lambda r, s : -256*r*s - 256*s**2 + 192*s, lambda r, s : -256*r**2 - 768*r*s + 384*r - 512*s**2 + 576*s - 416/3]
ddN13 = [lambda r, s : 768*r*s + 512*s**2 - 448*s, lambda r, s : 512*r**2 + 768*r*s - 448*r]
ddN14 = [lambda r, s : -768*r*s - 256*s**2 + 320*s, lambda r, s : -256*r**2 + 64*r]
ddN15 = [lambda r, s : -256*s**2 + 64*s, lambda r, s : -256*r**2 - 768*r*s + 320*r]

dddN1 = [lambda r, s : 256*r + 256*s - 160, lambda r, s : 256*r + 256*s - 160]
dddN2 = [lambda r, s : 256*r - 96, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 256*s - 96]
dddN4 = [lambda r, s : -1024*r - 768*s + 576, lambda r, s : -256*r]
dddN5 = [lambda r, s : 1536*r + 768*s - 768, lambda r, s : 0]
dddN6 = [lambda r, s : -1024*r - 256*s + 448, lambda r, s : 0]
dddN7 = [lambda r, s : 256*s, lambda r, s : 0]
dddN8 = [lambda r, s : 0, lambda r, s : 0]
dddN9 = [lambda r, s : 0, lambda r, s : 256*r]
dddN10 = [lambda r, s : 0, lambda r, s : -256*r - 1024*s + 448]
dddN11 = [lambda r, s : 0, lambda r, s : 768*r + 1536*s - 768]
dddN12 = [lambda r, s : -256*s, lambda r, s : -768*r - 1024*s + 576]
dddN13 = [lambda r, s : 768*s, lambda r, s : 768*r]
dddN14 = [lambda r, s : -768*s, lambda r, s : 0]
dddN15 = [lambda r, s : 0, lambda r, s : -768*r]

ddddN1 = [lambda r, s : 256, lambda r, s : 256]
ddddN2 = [lambda r, s : 256, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 256]
ddddN4 = [lambda r, s : -1024, lambda r, s : 0]
ddddN5 = [lambda r, s : 1536, lambda r, s : 0]
ddddN6 = [lambda r, s : -1024, lambda r, s : 0]
ddddN7 = [lambda r, s : 0, lambda r, s : 0]
ddddN8 = [lambda r, s : 0, lambda r, s : 0]
ddddN9 = [lambda r, s : 0, lambda r, s : 0]
ddddN10 = [lambda r, s : 0, lambda r, s : -1024]
ddddN11 = [lambda r, s : 0, lambda r, s : 1536]
ddddN12 = [lambda r, s : 0, lambda r, s : -1024]
ddddN13 = [lambda r, s : 0, lambda r, s : 0]
ddddN14 = [lambda r, s : 0, lambda r, s : 0]
ddddN15 = [lambda r, s : 0, lambda r, s : 0]



==================== QUAD4 ====================
N1 = lambda r, s : (r - 1)*(s - 1)/4
N2 = lambda r, s : -(r + 1)*(s - 1)/4
N3 = lambda r, s : (r + 1)*(s + 1)/4
N4 = lambda r, s : -(r - 1)*(s + 1)/4

dN1 = [lambda r, s : s/4 - 1/4, lambda r, s : r/4 - 1/4]
dN2 = [lambda r, s : 1/4 - s/4, lambda r, s : -r/4 - 1/4]
dN3 = [lambda r, s : s/4 + 1/4, lambda r, s : r/4 + 1/4]
dN4 = [lambda r, s : -s/4 - 1/4, lambda r, s : 1/4 - r/4]

ddN1 = [lambda r, s : 0, lambda r, s : 0]
ddN2 = [lambda r, s : 0, lambda r, s : 0]
ddN3 = [lambda r, s : 0, lambda r, s : 0]
ddN4 = [lambda r, s : 0, lambda r, s : 0]

dddN1 = [lambda r, s : 0, lambda r, s : 0]
dddN2 = [lambda r, s : 0, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 0]
dddN4 = [lambda r, s : 0, lambda r, s : 0]

ddddN1 = [lambda r, s : 0, lambda r, s : 0]
ddddN2 = [lambda r, s : 0, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 0]
ddddN4 = [lambda r, s : 0, lambda r, s : 0]



==================== QUAD8 ====================
N1 = lambda r, s : -(r - 1)*(s - 1)*(r + s + 1)/4
N2 = lambda r, s : -(r + 1)*(s - 1)*(r - s - 1)/4
N3 = lambda r, s : (r + 1)*(s + 1)*(r + s - 1)/4
N4 = lambda r, s : (r - 1)*(s + 1)*(r - s + 1)/4
N5 = lambda r, s : (r - 1)*(r + 1)*(s - 1)/2
N6 = lambda r, s : -(r + 1)*(s - 1)*(s + 1)/2
N7 = lambda r, s : -(r - 1)*(r + 1)*(s + 1)/2
N8 = lambda r, s : (r - 1)*(s - 1)*(s + 1)/2

dN1 = [lambda r, s : -(r - 1)*(s - 1)/4 - (s - 1)*(r + s + 1)/4, lambda r, s : -(r - 1)*(s - 1)/4 - (r - 1)*(r + s + 1)/4]
dN2 = [lambda r, s : -(r + 1)*(s - 1)/4 - (s - 1)*(r - s - 1)/4, lambda r, s : (r + 1)*(s - 1)/4 - (r + 1)*(r - s - 1)/4]
dN3 = [lambda r, s : (r + 1)*(s + 1)/4 + (s + 1)*(r + s - 1)/4, lambda r, s : (r + 1)*(s + 1)/4 + (r + 1)*(r + s - 1)/4]
dN4 = [lambda r, s : (r - 1)*(s + 1)/4 + (s + 1)*(r - s + 1)/4, lambda r, s : -(r - 1)*(s + 1)/4 + (r - 1)*(r - s + 1)/4]
dN5 = [lambda r, s : (r - 1)*(s - 1)/2 + (r + 1)*(s - 1)/2, lambda r, s : (r - 1)*(r + 1)/2]
dN6 = [lambda r, s : -(s - 1)*(s + 1)/2, lambda r, s : -(r + 1)*(s - 1)/2 - (r + 1)*(s + 1)/2]
dN7 = [lambda r, s : -(r - 1)*(s + 1)/2 - (r + 1)*(s + 1)/2, lambda r, s : -(r - 1)*(r + 1)/2]
dN8 = [lambda r, s : (s - 1)*(s + 1)/2, lambda r, s : (r - 1)*(s - 1)/2 + (r - 1)*(s + 1)/2]

ddN1 = [lambda r, s : 1/2 - s/2, lambda r, s : 1/2 - r/2]
ddN2 = [lambda r, s : 1/2 - s/2, lambda r, s : r/2 + 1/2]
ddN3 = [lambda r, s : s/2 + 1/2, lambda r, s : r/2 + 1/2]
ddN4 = [lambda r, s : s/2 + 1/2, lambda r, s : 1/2 - r/2]
ddN5 = [lambda r, s : s - 1, lambda r, s : 0]
ddN6 = [lambda r, s : 0, lambda r, s : -r - 1]
ddN7 = [lambda r, s : -s - 1, lambda r, s : 0]
ddN8 = [lambda r, s : 0, lambda r, s : r - 1]

dddN1 = [lambda r, s : 0, lambda r, s : 0]
dddN2 = [lambda r, s : 0, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 0]
dddN4 = [lambda r, s : 0, lambda r, s : 0]
dddN5 = [lambda r, s : 0, lambda r, s : 0]
dddN6 = [lambda r, s : 0, lambda r, s : 0]
dddN7 = [lambda r, s : 0, lambda r, s : 0]
dddN8 = [lambda r, s : 0, lambda r, s : 0]

ddddN1 = [lambda r, s : 0, lambda r, s : 0]
ddddN2 = [lambda r, s : 0, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 0]
ddddN4 = [lambda r, s : 0, lambda r, s : 0]
ddddN5 = [lambda r, s : 0, lambda r, s : 0]
ddddN6 = [lambda r, s : 0, lambda r, s : 0]
ddddN7 = [lambda r, s : 0, lambda r, s : 0]
ddddN8 = [lambda r, s : 0, lambda r, s : 0]



==================== QUAD9 ====================
N1 = lambda r, s : r*s*(r - 1)*(s - 1)/4
N2 = lambda r, s : r*s*(r + 1)*(s - 1)/4
N3 = lambda r, s : r*s*(r + 1)*(s + 1)/4
N4 = lambda r, s : r*s*(r - 1)*(s + 1)/4
N5 = lambda r, s : -s*(r - 1)*(r + 1)*(s - 1)/2
N6 = lambda r, s : -r*(r + 1)*(s - 1)*(s + 1)/2
N7 = lambda r, s : -s*(r - 1)*(r + 1)*(s + 1)/2
N8 = lambda r, s : -r*(r - 1)*(s - 1)*(s + 1)/2
N9 = lambda r, s : (r - 1)*(r + 1)*(s - 1)*(s + 1)

dN1 = [lambda r, s : r*s*(s - 1)/4 + s*(r - 1)*(s - 1)/4, lambda r, s : r*s*(r - 1)/4 + r*(r - 1)*(s - 1)/4]
dN2 = [lambda r, s : r*s*(s - 1)/4 + s*(r + 1)*(s - 1)/4, lambda r, s : r*s*(r + 1)/4 + r*(r + 1)*(s - 1)/4]
dN3 = [lambda r, s : r*s*(s + 1)/4 + s*(r + 1)*(s + 1)/4, lambda r, s : r*s*(r + 1)/4 + r*(r + 1)*(s + 1)/4]
dN4 = [lambda r, s : r*s*(s + 1)/4 + s*(r - 1)*(s + 1)/4, lambda r, s : r*s*(r - 1)/4 + r*(r - 1)*(s + 1)/4]
dN5 = [lambda r, s : -s*(r - 1)*(s - 1)/2 - s*(r + 1)*(s - 1)/2, lambda r, s : -s*(r - 1)*(r + 1)/2 - (r - 1)*(r + 1)*(s - 1)/2]
dN6 = [lambda r, s : -r*(s - 1)*(s + 1)/2 - (r + 1)*(s - 1)*(s + 1)/2, lambda r, s : -r*(r + 1)*(s - 1)/2 - r*(r + 1)*(s + 1)/2]
dN7 = [lambda r, s : -s*(r - 1)*(s + 1)/2 - s*(r + 1)*(s + 1)/2, lambda r, s : -s*(r - 1)*(r + 1)/2 - (r - 1)*(r + 1)*(s + 1)/2]
dN8 = [lambda r, s : -r*(s - 1)*(s + 1)/2 - (r - 1)*(s - 1)*(s + 1)/2, lambda r, s : -r*(r - 1)*(s - 1)/2 - r*(r - 1)*(s + 1)/2]
dN9 = [lambda r, s : (r - 1)*(s - 1)*(s + 1) + (r + 1)*(s - 1)*(s + 1), lambda r, s : (r - 1)*(r + 1)*(s - 1) + (r - 1)*(r + 1)*(s + 1)]

ddN1 = [lambda r, s : s*(s - 1)/2, lambda r, s : r*(r - 1)/2]
ddN2 = [lambda r, s : s*(s - 1)/2, lambda r, s : r*(r + 1)/2]
ddN3 = [lambda r, s : s*(s + 1)/2, lambda r, s : r*(r + 1)/2]
ddN4 = [lambda r, s : s*(s + 1)/2, lambda r, s : r*(r - 1)/2]
ddN5 = [lambda r, s : -s*(s - 1), lambda r, s : -(r - 1)*(r + 1)]
ddN6 = [lambda r, s : -(s - 1)*(s + 1), lambda r, s : -r*(r + 1)]
ddN7 = [lambda r, s : -s*(s + 1), lambda r, s : -(r - 1)*(r + 1)]
ddN8 = [lambda r, s : -(s - 1)*(s + 1), lambda r, s : -r*(r - 1)]
ddN9 = [lambda r, s : 2*(s - 1)*(s + 1), lambda r, s : 2*(r - 1)*(r + 1)]

dddN1 = [lambda r, s : 0, lambda r, s : 0]
dddN2 = [lambda r, s : 0, lambda r, s : 0]
dddN3 = [lambda r, s : 0, lambda r, s : 0]
dddN4 = [lambda r, s : 0, lambda r, s : 0]
dddN5 = [lambda r, s : 0, lambda r, s : 0]
dddN6 = [lambda r, s : 0, lambda r, s : 0]
dddN7 = [lambda r, s : 0, lambda r, s : 0]
dddN8 = [lambda r, s : 0, lambda r, s : 0]
dddN9 = [lambda r, s : 0, lambda r, s : 0]

ddddN1 = [lambda r, s : 0, lambda r, s : 0]
ddddN2 = [lambda r, s : 0, lambda r, s : 0]
ddddN3 = [lambda r, s : 0, lambda r, s : 0]
ddddN4 = [lambda r, s : 0, lambda r, s : 0]
ddddN5 = [lambda r, s : 0, lambda r, s : 0]
ddddN6 = [lambda r, s : 0, lambda r, s : 0]
ddddN7 = [lambda r, s : 0, lambda r, s : 0]
ddddN8 = [lambda r, s : 0, lambda r, s : 0]
ddddN9 = [lambda r, s : 0, lambda r, s : 0]



=================== TETRA4 ===================
N1 = lambda r, s, t : -r - s - t + 1
N2 = lambda r, s, t : r
N3 = lambda r, s, t : s
N4 = lambda r, s, t : t

dN1 = [lambda r, s, t : -1, lambda r, s, t : -1, lambda r, s, t : -1]
dN2 = [lambda r, s, t : 1, lambda r, s, t : 0, lambda r, s, t : 0]
dN3 = [lambda r, s, t : 0, lambda r, s, t : 1, lambda r, s, t : 0]
dN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 1]

ddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



=================== TETRA10 ===================
N1 = lambda r, s, t : (r + s + t - 1)*(2*r + 2*s + 2*t - 1)
N2 = lambda r, s, t : r*(2*r - 1)
N3 = lambda r, s, t : s*(2*s - 1)
N4 = lambda r, s, t : t*(2*t - 1)
N5 = lambda r, s, t : -4*r*(r + s + t - 1)
N6 = lambda r, s, t : 4*r*s
N7 = lambda r, s, t : -4*s*(r + s + t - 1)
N8 = lambda r, s, t : -4*t*(r + s + t - 1)
N9 = lambda r, s, t : 4*s*t
N10 = lambda r, s, t : 4*r*t

dN1 = [lambda r, s, t : 4*r + 4*s + 4*t - 3, lambda r, s, t : 4*r + 4*s + 4*t - 3, lambda r, s, t : 4*r + 4*s + 4*t - 3]
dN2 = [lambda r, s, t : 4*r - 1, lambda r, s, t : 0, lambda r, s, t : 0]
dN3 = [lambda r, s, t : 0, lambda r, s, t : 4*s - 1, lambda r, s, t : 0]
dN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 4*t - 1]
dN5 = [lambda r, s, t : -8*r - 4*s - 4*t + 4, lambda r, s, t : -4*r, lambda r, s, t : -4*r]
dN6 = [lambda r, s, t : 4*s, lambda r, s, t : 4*r, lambda r, s, t : 0]
dN7 = [lambda r, s, t : -4*s, lambda r, s, t : -4*r - 8*s - 4*t + 4, lambda r, s, t : -4*s]
dN8 = [lambda r, s, t : -4*t, lambda r, s, t : -4*t, lambda r, s, t : -4*r - 4*s - 8*t + 4]
dN9 = [lambda r, s, t : 0, lambda r, s, t : 4*t, lambda r, s, t : 4*s]
dN10 = [lambda r, s, t : 4*t, lambda r, s, t : 0, lambda r, s, t : 4*r]

ddN1 = [lambda r, s, t : 4, lambda r, s, t : 4, lambda r, s, t : 4]
ddN2 = [lambda r, s, t : 4, lambda r, s, t : 0, lambda r, s, t : 0]
ddN3 = [lambda r, s, t : 0, lambda r, s, t : 4, lambda r, s, t : 0]
ddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 4]
ddN5 = [lambda r, s, t : -8, lambda r, s, t : 0, lambda r, s, t : 0]
ddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN7 = [lambda r, s, t : 0, lambda r, s, t : -8, lambda r, s, t : 0]
ddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : -8]
ddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



==================== HEXA8 ====================
N1 = lambda r, s, t : -(r - 1)*(s - 1)*(t - 1)/8
N2 = lambda r, s, t : (r + 1)*(s - 1)*(t - 1)/8
N3 = lambda r, s, t : -(r + 1)*(s + 1)*(t - 1)/8
N4 = lambda r, s, t : (r - 1)*(s + 1)*(t - 1)/8
N5 = lambda r, s, t : (r - 1)*(s - 1)*(t + 1)/8
N6 = lambda r, s, t : -(r + 1)*(s - 1)*(t + 1)/8
N7 = lambda r, s, t : (r + 1)*(s + 1)*(t + 1)/8
N8 = lambda r, s, t : -(r - 1)*(s + 1)*(t + 1)/8

dN1 = [lambda r, s, t : -(s - 1)*(t - 1)/8, lambda r, s, t : -(r - 1)*(t - 1)/8, lambda r, s, t : -(r - 1)*(s - 1)/8]
dN2 = [lambda r, s, t : (s - 1)*(t - 1)/8, lambda r, s, t : (r + 1)*(t - 1)/8, lambda r, s, t : (r + 1)*(s - 1)/8]
dN3 = [lambda r, s, t : -(s + 1)*(t - 1)/8, lambda r, s, t : -(r + 1)*(t - 1)/8, lambda r, s, t : -(r + 1)*(s + 1)/8]
dN4 = [lambda r, s, t : (s + 1)*(t - 1)/8, lambda r, s, t : (r - 1)*(t - 1)/8, lambda r, s, t : (r - 1)*(s + 1)/8]
dN5 = [lambda r, s, t : (s - 1)*(t + 1)/8, lambda r, s, t : (r - 1)*(t + 1)/8, lambda r, s, t : (r - 1)*(s - 1)/8]
dN6 = [lambda r, s, t : -(s - 1)*(t + 1)/8, lambda r, s, t : -(r + 1)*(t + 1)/8, lambda r, s, t : -(r + 1)*(s - 1)/8]
dN7 = [lambda r, s, t : (s + 1)*(t + 1)/8, lambda r, s, t : (r + 1)*(t + 1)/8, lambda r, s, t : (r + 1)*(s + 1)/8]
dN8 = [lambda r, s, t : -(s + 1)*(t + 1)/8, lambda r, s, t : -(r - 1)*(t + 1)/8, lambda r, s, t : -(r - 1)*(s + 1)/8]

ddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



=================== HEXA20 ===================
N1 = lambda r, s, t : (r - 1)*(s - 1)*(t - 1)*(r + s + t + 2)/8
N2 = lambda r, s, t : (r + 1)*(s - 1)*(t - 1)*(r - s - t - 2)/8
N3 = lambda r, s, t : -(r + 1)*(s + 1)*(t - 1)*(r + s - t - 2)/8
N4 = lambda r, s, t : -(r - 1)*(s + 1)*(t - 1)*(r - s + t + 2)/8
N5 = lambda r, s, t : -(r - 1)*(s - 1)*(t + 1)*(r + s - t + 2)/8
N6 = lambda r, s, t : -(r + 1)*(s - 1)*(t + 1)*(r - s + t - 2)/8
N7 = lambda r, s, t : (r + 1)*(s + 1)*(t + 1)*(r + s + t - 2)/8
N8 = lambda r, s, t : (r - 1)*(s + 1)*(t + 1)*(r - s - t + 2)/8
N9 = lambda r, s, t : -(r - 1)*(r + 1)*(s - 1)*(t - 1)/4
N10 = lambda r, s, t : -(r - 1)*(s - 1)*(s + 1)*(t - 1)/4
N11 = lambda r, s, t : -(r - 1)*(s - 1)*(t - 1)*(t + 1)/4
N12 = lambda r, s, t : (r + 1)*(s - 1)*(s + 1)*(t - 1)/4
N13 = lambda r, s, t : (r + 1)*(s - 1)*(t - 1)*(t + 1)/4
N14 = lambda r, s, t : (r - 1)*(r + 1)*(s + 1)*(t - 1)/4
N15 = lambda r, s, t : -(r + 1)*(s + 1)*(t - 1)*(t + 1)/4
N16 = lambda r, s, t : (r - 1)*(s + 1)*(t - 1)*(t + 1)/4
N17 = lambda r, s, t : (r - 1)*(r + 1)*(s - 1)*(t + 1)/4
N18 = lambda r, s, t : (r - 1)*(s - 1)*(s + 1)*(t + 1)/4
N19 = lambda r, s, t : -(r + 1)*(s - 1)*(s + 1)*(t + 1)/4
N20 = lambda r, s, t : -(r - 1)*(r + 1)*(s + 1)*(t + 1)/4

dN1 = [lambda r, s, t : (r - 1)*(s - 1)*(t - 1)/8 + (s - 1)*(t - 1)*(r + s + t + 2)/8, lambda r, s, t : (r - 1)*(s - 1)*(t - 1)/8 + (r - 1)*(t - 1)*(r + s + t + 2)/8, lambda r, s, t : (r - 1)*(s - 1)*(t - 1)/8 + (r - 1)*(s - 1)*(r + s + t + 2)/8]
dN2 = [lambda r, s, t : (r + 1)*(s - 1)*(t - 1)/8 + (s - 1)*(t - 1)*(r - s - t - 2)/8, lambda r, s, t : -(r + 1)*(s - 1)*(t - 1)/8 + (r + 1)*(t - 1)*(r - s - t - 2)/8, lambda r, s, t : -(r + 1)*(s - 1)*(t - 1)/8 + (r + 1)*(s - 1)*(r - s - t - 2)/8]
dN3 = [lambda r, s, t : -(r + 1)*(s + 1)*(t - 1)/8 - (s + 1)*(t - 1)*(r + s - t - 2)/8, lambda r, s, t : -(r + 1)*(s + 1)*(t - 1)/8 - (r + 1)*(t - 1)*(r + s - t - 2)/8, lambda r, s, t : (r + 1)*(s + 1)*(t - 1)/8 - (r + 1)*(s + 1)*(r + s - t - 2)/8]
dN4 = [lambda r, s, t : -(r - 1)*(s + 1)*(t - 1)/8 - (s + 1)*(t - 1)*(r - s + t + 2)/8, lambda r, s, t : (r - 1)*(s + 1)*(t - 1)/8 - (r - 1)*(t - 1)*(r - s + t + 2)/8, lambda r, s, t : -(r - 1)*(s + 1)*(t - 1)/8 - (r - 1)*(s + 1)*(r - s + t + 2)/8]
dN5 = [lambda r, s, t : -(r - 1)*(s - 1)*(t + 1)/8 - (s - 1)*(t + 1)*(r + s - t + 2)/8, lambda r, s, t : -(r - 1)*(s - 1)*(t + 1)/8 - (r - 1)*(t + 1)*(r + s - t + 2)/8, lambda r, s, t : (r - 1)*(s - 1)*(t + 1)/8 - (r - 1)*(s - 1)*(r + s - t + 2)/8]
dN6 = [lambda r, s, t : -(r + 1)*(s - 1)*(t + 1)/8 - (s - 1)*(t + 1)*(r - s + t - 2)/8, lambda r, s, t : (r + 1)*(s - 1)*(t + 1)/8 - (r + 1)*(t + 1)*(r - s + t - 2)/8, lambda r, s, t : -(r + 1)*(s - 1)*(t + 1)/8 - (r + 1)*(s - 1)*(r - s + t - 2)/8]
dN7 = [lambda r, s, t : (r + 1)*(s + 1)*(t + 1)/8 + (s + 1)*(t + 1)*(r + s + t - 2)/8, lambda r, s, t : (r + 1)*(s + 1)*(t + 1)/8 + (r + 1)*(t + 1)*(r + s + t - 2)/8, lambda r, s, t : (r + 1)*(s + 1)*(t + 1)/8 + (r + 1)*(s + 1)*(r + s + t - 2)/8]
dN8 = [lambda r, s, t : (r - 1)*(s + 1)*(t + 1)/8 + (s + 1)*(t + 1)*(r - s - t + 2)/8, lambda r, s, t : -(r - 1)*(s + 1)*(t + 1)/8 + (r - 1)*(t + 1)*(r - s - t + 2)/8, lambda r, s, t : -(r - 1)*(s + 1)*(t + 1)/8 + (r - 1)*(s + 1)*(r - s - t + 2)/8]
dN9 = [lambda r, s, t : -(r - 1)*(s - 1)*(t - 1)/4 - (r + 1)*(s - 1)*(t - 1)/4, lambda r, s, t : -(r - 1)*(r + 1)*(t - 1)/4, lambda r, s, t : -(r - 1)*(r + 1)*(s - 1)/4]
dN10 = [lambda r, s, t : -(s - 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -(r - 1)*(s - 1)*(t - 1)/4 - (r - 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -(r - 1)*(s - 1)*(s + 1)/4]
dN11 = [lambda r, s, t : -(s - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(s - 1)*(t - 1)/4 - (r - 1)*(s - 1)*(t + 1)/4]
dN12 = [lambda r, s, t : (s - 1)*(s + 1)*(t - 1)/4, lambda r, s, t : (r + 1)*(s - 1)*(t - 1)/4 + (r + 1)*(s + 1)*(t - 1)/4, lambda r, s, t : (r + 1)*(s - 1)*(s + 1)/4]
dN13 = [lambda r, s, t : (s - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : (r + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : (r + 1)*(s - 1)*(t - 1)/4 + (r + 1)*(s - 1)*(t + 1)/4]
dN14 = [lambda r, s, t : (r - 1)*(s + 1)*(t - 1)/4 + (r + 1)*(s + 1)*(t - 1)/4, lambda r, s, t : (r - 1)*(r + 1)*(t - 1)/4, lambda r, s, t : (r - 1)*(r + 1)*(s + 1)/4]
dN15 = [lambda r, s, t : -(s + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -(r + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -(r + 1)*(s + 1)*(t - 1)/4 - (r + 1)*(s + 1)*(t + 1)/4]
dN16 = [lambda r, s, t : (s + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(s + 1)*(t - 1)/4 + (r - 1)*(s + 1)*(t + 1)/4]
dN17 = [lambda r, s, t : (r - 1)*(s - 1)*(t + 1)/4 + (r + 1)*(s - 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(r + 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(r + 1)*(s - 1)/4]
dN18 = [lambda r, s, t : (s - 1)*(s + 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(s - 1)*(t + 1)/4 + (r - 1)*(s + 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(s - 1)*(s + 1)/4]
dN19 = [lambda r, s, t : -(s - 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -(r + 1)*(s - 1)*(t + 1)/4 - (r + 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -(r + 1)*(s - 1)*(s + 1)/4]
dN20 = [lambda r, s, t : -(r - 1)*(s + 1)*(t + 1)/4 - (r + 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(r + 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(r + 1)*(s + 1)/4]

ddN1 = [lambda r, s, t : (s - 1)*(t - 1)/4, lambda r, s, t : (r - 1)*(t - 1)/4, lambda r, s, t : (r - 1)*(s - 1)/4]
ddN2 = [lambda r, s, t : (s - 1)*(t - 1)/4, lambda r, s, t : -(r + 1)*(t - 1)/4, lambda r, s, t : -(r + 1)*(s - 1)/4]
ddN3 = [lambda r, s, t : -(s + 1)*(t - 1)/4, lambda r, s, t : -(r + 1)*(t - 1)/4, lambda r, s, t : (r + 1)*(s + 1)/4]
ddN4 = [lambda r, s, t : -(s + 1)*(t - 1)/4, lambda r, s, t : (r - 1)*(t - 1)/4, lambda r, s, t : -(r - 1)*(s + 1)/4]
ddN5 = [lambda r, s, t : -(s - 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(t + 1)/4, lambda r, s, t : (r - 1)*(s - 1)/4]
ddN6 = [lambda r, s, t : -(s - 1)*(t + 1)/4, lambda r, s, t : (r + 1)*(t + 1)/4, lambda r, s, t : -(r + 1)*(s - 1)/4]
ddN7 = [lambda r, s, t : (s + 1)*(t + 1)/4, lambda r, s, t : (r + 1)*(t + 1)/4, lambda r, s, t : (r + 1)*(s + 1)/4]
ddN8 = [lambda r, s, t : (s + 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(t + 1)/4, lambda r, s, t : -(r - 1)*(s + 1)/4]
ddN9 = [lambda r, s, t : -(s - 1)*(t - 1)/2, lambda r, s, t : 0, lambda r, s, t : 0]
ddN10 = [lambda r, s, t : 0, lambda r, s, t : -(r - 1)*(t - 1)/2, lambda r, s, t : 0]
ddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : -(r - 1)*(s - 1)/2]
ddN12 = [lambda r, s, t : 0, lambda r, s, t : (r + 1)*(t - 1)/2, lambda r, s, t : 0]
ddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : (r + 1)*(s - 1)/2]
ddN14 = [lambda r, s, t : (s + 1)*(t - 1)/2, lambda r, s, t : 0, lambda r, s, t : 0]
ddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : -(r + 1)*(s + 1)/2]
ddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : (r - 1)*(s + 1)/2]
ddN17 = [lambda r, s, t : (s - 1)*(t + 1)/2, lambda r, s, t : 0, lambda r, s, t : 0]
ddN18 = [lambda r, s, t : 0, lambda r, s, t : (r - 1)*(t + 1)/2, lambda r, s, t : 0]
ddN19 = [lambda r, s, t : 0, lambda r, s, t : -(r + 1)*(t + 1)/2, lambda r, s, t : 0]
ddN20 = [lambda r, s, t : -(s + 1)*(t + 1)/2, lambda r, s, t : 0, lambda r, s, t : 0]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN17 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN19 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN20 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN17 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN19 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN20 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



=================== HEXA27 ===================
N1 = lambda r, s, t : r*s*t*(r - 1)*(s - 1)*(t - 1)/8
N2 = lambda r, s, t : r*s*t*(r + 1)*(s - 1)*(t - 1)/8
N3 = lambda r, s, t : r*s*t*(r + 1)*(s + 1)*(t - 1)/8
N4 = lambda r, s, t : r*s*t*(r - 1)*(s + 1)*(t - 1)/8
N5 = lambda r, s, t : r*s*t*(r - 1)*(s - 1)*(t + 1)/8
N6 = lambda r, s, t : r*s*t*(r + 1)*(s - 1)*(t + 1)/8
N7 = lambda r, s, t : r*s*t*(r + 1)*(s + 1)*(t + 1)/8
N8 = lambda r, s, t : r*s*t*(r - 1)*(s + 1)*(t + 1)/8
N9 = lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s - 1)*(t - 1)/4
N10 = lambda r, s, t : -r*t*(r - 1)*(s - 1)*(s + 1)*(t - 1)/4
N11 = lambda r, s, t : -r*s*(r - 1)*(s - 1)*(t - 1)*(t + 1)/4
N12 = lambda r, s, t : -r*t*(r + 1)*(s - 1)*(s + 1)*(t - 1)/4
N13 = lambda r, s, t : -r*s*(r + 1)*(s - 1)*(t - 1)*(t + 1)/4
N14 = lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s + 1)*(t - 1)/4
N15 = lambda r, s, t : -r*s*(r + 1)*(s + 1)*(t - 1)*(t + 1)/4
N16 = lambda r, s, t : -r*s*(r - 1)*(s + 1)*(t - 1)*(t + 1)/4
N17 = lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s - 1)*(t + 1)/4
N18 = lambda r, s, t : -r*t*(r - 1)*(s - 1)*(s + 1)*(t + 1)/4
N19 = lambda r, s, t : -r*t*(r + 1)*(s - 1)*(s + 1)*(t + 1)/4
N20 = lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s + 1)*(t + 1)/4
N21 = lambda r, s, t : t*(r - 1)*(r + 1)*(s - 1)*(s + 1)*(t - 1)/2
N22 = lambda r, s, t : s*(r - 1)*(r + 1)*(s - 1)*(t - 1)*(t + 1)/2
N23 = lambda r, s, t : r*(r - 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1)/2
N24 = lambda r, s, t : r*(r + 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1)/2
N25 = lambda r, s, t : s*(r - 1)*(r + 1)*(s + 1)*(t - 1)*(t + 1)/2
N26 = lambda r, s, t : t*(r - 1)*(r + 1)*(s - 1)*(s + 1)*(t + 1)/2
N27 = lambda r, s, t : -(r - 1)*(r + 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1)

dN1 = [lambda r, s, t : r*s*t*(s - 1)*(t - 1)/8 + s*t*(r - 1)*(s - 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r - 1)*(t - 1)/8 + r*t*(r - 1)*(s - 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r - 1)*(s - 1)/8 + r*s*(r - 1)*(s - 1)*(t - 1)/8]
dN2 = [lambda r, s, t : r*s*t*(s - 1)*(t - 1)/8 + s*t*(r + 1)*(s - 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r + 1)*(t - 1)/8 + r*t*(r + 1)*(s - 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r + 1)*(s - 1)/8 + r*s*(r + 1)*(s - 1)*(t - 1)/8]
dN3 = [lambda r, s, t : r*s*t*(s + 1)*(t - 1)/8 + s*t*(r + 1)*(s + 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r + 1)*(t - 1)/8 + r*t*(r + 1)*(s + 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r + 1)*(s + 1)/8 + r*s*(r + 1)*(s + 1)*(t - 1)/8]
dN4 = [lambda r, s, t : r*s*t*(s + 1)*(t - 1)/8 + s*t*(r - 1)*(s + 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r - 1)*(t - 1)/8 + r*t*(r - 1)*(s + 1)*(t - 1)/8, lambda r, s, t : r*s*t*(r - 1)*(s + 1)/8 + r*s*(r - 1)*(s + 1)*(t - 1)/8]
dN5 = [lambda r, s, t : r*s*t*(s - 1)*(t + 1)/8 + s*t*(r - 1)*(s - 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r - 1)*(t + 1)/8 + r*t*(r - 1)*(s - 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r - 1)*(s - 1)/8 + r*s*(r - 1)*(s - 1)*(t + 1)/8]
dN6 = [lambda r, s, t : r*s*t*(s - 1)*(t + 1)/8 + s*t*(r + 1)*(s - 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r + 1)*(t + 1)/8 + r*t*(r + 1)*(s - 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r + 1)*(s - 1)/8 + r*s*(r + 1)*(s - 1)*(t + 1)/8]
dN7 = [lambda r, s, t : r*s*t*(s + 1)*(t + 1)/8 + s*t*(r + 1)*(s + 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r + 1)*(t + 1)/8 + r*t*(r + 1)*(s + 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r + 1)*(s + 1)/8 + r*s*(r + 1)*(s + 1)*(t + 1)/8]
dN8 = [lambda r, s, t : r*s*t*(s + 1)*(t + 1)/8 + s*t*(r - 1)*(s + 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r - 1)*(t + 1)/8 + r*t*(r - 1)*(s + 1)*(t + 1)/8, lambda r, s, t : r*s*t*(r - 1)*(s + 1)/8 + r*s*(r - 1)*(s + 1)*(t + 1)/8]
dN9 = [lambda r, s, t : -s*t*(r - 1)*(s - 1)*(t - 1)/4 - s*t*(r + 1)*(s - 1)*(t - 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(t - 1)/4 - t*(r - 1)*(r + 1)*(s - 1)*(t - 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s - 1)/4 - s*(r - 1)*(r + 1)*(s - 1)*(t - 1)/4]
dN10 = [lambda r, s, t : -r*t*(s - 1)*(s + 1)*(t - 1)/4 - t*(r - 1)*(s - 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -r*t*(r - 1)*(s - 1)*(t - 1)/4 - r*t*(r - 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -r*t*(r - 1)*(s - 1)*(s + 1)/4 - r*(r - 1)*(s - 1)*(s + 1)*(t - 1)/4]
dN11 = [lambda r, s, t : -r*s*(s - 1)*(t - 1)*(t + 1)/4 - s*(r - 1)*(s - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r - 1)*(t - 1)*(t + 1)/4 - r*(r - 1)*(s - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r - 1)*(s - 1)*(t - 1)/4 - r*s*(r - 1)*(s - 1)*(t + 1)/4]
dN12 = [lambda r, s, t : -r*t*(s - 1)*(s + 1)*(t - 1)/4 - t*(r + 1)*(s - 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -r*t*(r + 1)*(s - 1)*(t - 1)/4 - r*t*(r + 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -r*t*(r + 1)*(s - 1)*(s + 1)/4 - r*(r + 1)*(s - 1)*(s + 1)*(t - 1)/4]
dN13 = [lambda r, s, t : -r*s*(s - 1)*(t - 1)*(t + 1)/4 - s*(r + 1)*(s - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r + 1)*(t - 1)*(t + 1)/4 - r*(r + 1)*(s - 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r + 1)*(s - 1)*(t - 1)/4 - r*s*(r + 1)*(s - 1)*(t + 1)/4]
dN14 = [lambda r, s, t : -s*t*(r - 1)*(s + 1)*(t - 1)/4 - s*t*(r + 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(t - 1)/4 - t*(r - 1)*(r + 1)*(s + 1)*(t - 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s + 1)/4 - s*(r - 1)*(r + 1)*(s + 1)*(t - 1)/4]
dN15 = [lambda r, s, t : -r*s*(s + 1)*(t - 1)*(t + 1)/4 - s*(r + 1)*(s + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r + 1)*(t - 1)*(t + 1)/4 - r*(r + 1)*(s + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r + 1)*(s + 1)*(t - 1)/4 - r*s*(r + 1)*(s + 1)*(t + 1)/4]
dN16 = [lambda r, s, t : -r*s*(s + 1)*(t - 1)*(t + 1)/4 - s*(r - 1)*(s + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r - 1)*(t - 1)*(t + 1)/4 - r*(r - 1)*(s + 1)*(t - 1)*(t + 1)/4, lambda r, s, t : -r*s*(r - 1)*(s + 1)*(t - 1)/4 - r*s*(r - 1)*(s + 1)*(t + 1)/4]
dN17 = [lambda r, s, t : -s*t*(r - 1)*(s - 1)*(t + 1)/4 - s*t*(r + 1)*(s - 1)*(t + 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(t + 1)/4 - t*(r - 1)*(r + 1)*(s - 1)*(t + 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s - 1)/4 - s*(r - 1)*(r + 1)*(s - 1)*(t + 1)/4]
dN18 = [lambda r, s, t : -r*t*(s - 1)*(s + 1)*(t + 1)/4 - t*(r - 1)*(s - 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -r*t*(r - 1)*(s - 1)*(t + 1)/4 - r*t*(r - 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -r*t*(r - 1)*(s - 1)*(s + 1)/4 - r*(r - 1)*(s - 1)*(s + 1)*(t + 1)/4]
dN19 = [lambda r, s, t : -r*t*(s - 1)*(s + 1)*(t + 1)/4 - t*(r + 1)*(s - 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -r*t*(r + 1)*(s - 1)*(t + 1)/4 - r*t*(r + 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -r*t*(r + 1)*(s - 1)*(s + 1)/4 - r*(r + 1)*(s - 1)*(s + 1)*(t + 1)/4]
dN20 = [lambda r, s, t : -s*t*(r - 1)*(s + 1)*(t + 1)/4 - s*t*(r + 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(t + 1)/4 - t*(r - 1)*(r + 1)*(s + 1)*(t + 1)/4, lambda r, s, t : -s*t*(r - 1)*(r + 1)*(s + 1)/4 - s*(r - 1)*(r + 1)*(s + 1)*(t + 1)/4]
dN21 = [lambda r, s, t : t*(r - 1)*(s - 1)*(s + 1)*(t - 1)/2 + t*(r + 1)*(s - 1)*(s + 1)*(t - 1)/2, lambda r, s, t : t*(r - 1)*(r + 1)*(s - 1)*(t - 1)/2 + t*(r - 1)*(r + 1)*(s + 1)*(t - 1)/2, lambda r, s, t : t*(r - 1)*(r + 1)*(s - 1)*(s + 1)/2 + (r - 1)*(r + 1)*(s - 1)*(s + 1)*(t - 1)/2]
dN22 = [lambda r, s, t : s*(r - 1)*(s - 1)*(t - 1)*(t + 1)/2 + s*(r + 1)*(s - 1)*(t - 1)*(t + 1)/2, lambda r, s, t : s*(r - 1)*(r + 1)*(t - 1)*(t + 1)/2 + (r - 1)*(r + 1)*(s - 1)*(t - 1)*(t + 1)/2, lambda r, s, t : s*(r - 1)*(r + 1)*(s - 1)*(t - 1)/2 + s*(r - 1)*(r + 1)*(s - 1)*(t + 1)/2]
dN23 = [lambda r, s, t : r*(s - 1)*(s + 1)*(t - 1)*(t + 1)/2 + (r - 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : r*(r - 1)*(s - 1)*(t - 1)*(t + 1)/2 + r*(r - 1)*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : r*(r - 1)*(s - 1)*(s + 1)*(t - 1)/2 + r*(r - 1)*(s - 1)*(s + 1)*(t + 1)/2]
dN24 = [lambda r, s, t : r*(s - 1)*(s + 1)*(t - 1)*(t + 1)/2 + (r + 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : r*(r + 1)*(s - 1)*(t - 1)*(t + 1)/2 + r*(r + 1)*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : r*(r + 1)*(s - 1)*(s + 1)*(t - 1)/2 + r*(r + 1)*(s - 1)*(s + 1)*(t + 1)/2]
dN25 = [lambda r, s, t : s*(r - 1)*(s + 1)*(t - 1)*(t + 1)/2 + s*(r + 1)*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : s*(r - 1)*(r + 1)*(t - 1)*(t + 1)/2 + (r - 1)*(r + 1)*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : s*(r - 1)*(r + 1)*(s + 1)*(t - 1)/2 + s*(r - 1)*(r + 1)*(s + 1)*(t + 1)/2]
dN26 = [lambda r, s, t : t*(r - 1)*(s - 1)*(s + 1)*(t + 1)/2 + t*(r + 1)*(s - 1)*(s + 1)*(t + 1)/2, lambda r, s, t : t*(r - 1)*(r + 1)*(s - 1)*(t + 1)/2 + t*(r - 1)*(r + 1)*(s + 1)*(t + 1)/2, lambda r, s, t : t*(r - 1)*(r + 1)*(s - 1)*(s + 1)/2 + (r - 1)*(r + 1)*(s - 1)*(s + 1)*(t + 1)/2]
dN27 = [lambda r, s, t : -(r - 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1) - (r + 1)*(s - 1)*(s + 1)*(t - 1)*(t + 1), lambda r, s, t : -(r - 1)*(r + 1)*(s - 1)*(t - 1)*(t + 1) - (r - 1)*(r + 1)*(s + 1)*(t - 1)*(t + 1), lambda r, s, t : -(r - 1)*(r + 1)*(s - 1)*(s + 1)*(t - 1) - (r - 1)*(r + 1)*(s - 1)*(s + 1)*(t + 1)]

ddN1 = [lambda r, s, t : s*t*(s - 1)*(t - 1)/4, lambda r, s, t : r*t*(r - 1)*(t - 1)/4, lambda r, s, t : r*s*(r - 1)*(s - 1)/4]
ddN2 = [lambda r, s, t : s*t*(s - 1)*(t - 1)/4, lambda r, s, t : r*t*(r + 1)*(t - 1)/4, lambda r, s, t : r*s*(r + 1)*(s - 1)/4]
ddN3 = [lambda r, s, t : s*t*(s + 1)*(t - 1)/4, lambda r, s, t : r*t*(r + 1)*(t - 1)/4, lambda r, s, t : r*s*(r + 1)*(s + 1)/4]
ddN4 = [lambda r, s, t : s*t*(s + 1)*(t - 1)/4, lambda r, s, t : r*t*(r - 1)*(t - 1)/4, lambda r, s, t : r*s*(r - 1)*(s + 1)/4]
ddN5 = [lambda r, s, t : s*t*(s - 1)*(t + 1)/4, lambda r, s, t : r*t*(r - 1)*(t + 1)/4, lambda r, s, t : r*s*(r - 1)*(s - 1)/4]
ddN6 = [lambda r, s, t : s*t*(s - 1)*(t + 1)/4, lambda r, s, t : r*t*(r + 1)*(t + 1)/4, lambda r, s, t : r*s*(r + 1)*(s - 1)/4]
ddN7 = [lambda r, s, t : s*t*(s + 1)*(t + 1)/4, lambda r, s, t : r*t*(r + 1)*(t + 1)/4, lambda r, s, t : r*s*(r + 1)*(s + 1)/4]
ddN8 = [lambda r, s, t : s*t*(s + 1)*(t + 1)/4, lambda r, s, t : r*t*(r - 1)*(t + 1)/4, lambda r, s, t : r*s*(r - 1)*(s + 1)/4]
ddN9 = [lambda r, s, t : -s*t*(s - 1)*(t - 1)/2, lambda r, s, t : -t*(r - 1)*(r + 1)*(t - 1)/2, lambda r, s, t : -s*(r - 1)*(r + 1)*(s - 1)/2]
ddN10 = [lambda r, s, t : -t*(s - 1)*(s + 1)*(t - 1)/2, lambda r, s, t : -r*t*(r - 1)*(t - 1)/2, lambda r, s, t : -r*(r - 1)*(s - 1)*(s + 1)/2]
ddN11 = [lambda r, s, t : -s*(s - 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*(r - 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*s*(r - 1)*(s - 1)/2]
ddN12 = [lambda r, s, t : -t*(s - 1)*(s + 1)*(t - 1)/2, lambda r, s, t : -r*t*(r + 1)*(t - 1)/2, lambda r, s, t : -r*(r + 1)*(s - 1)*(s + 1)/2]
ddN13 = [lambda r, s, t : -s*(s - 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*(r + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*s*(r + 1)*(s - 1)/2]
ddN14 = [lambda r, s, t : -s*t*(s + 1)*(t - 1)/2, lambda r, s, t : -t*(r - 1)*(r + 1)*(t - 1)/2, lambda r, s, t : -s*(r - 1)*(r + 1)*(s + 1)/2]
ddN15 = [lambda r, s, t : -s*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*(r + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*s*(r + 1)*(s + 1)/2]
ddN16 = [lambda r, s, t : -s*(s + 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*(r - 1)*(t - 1)*(t + 1)/2, lambda r, s, t : -r*s*(r - 1)*(s + 1)/2]
ddN17 = [lambda r, s, t : -s*t*(s - 1)*(t + 1)/2, lambda r, s, t : -t*(r - 1)*(r + 1)*(t + 1)/2, lambda r, s, t : -s*(r - 1)*(r + 1)*(s - 1)/2]
ddN18 = [lambda r, s, t : -t*(s - 1)*(s + 1)*(t + 1)/2, lambda r, s, t : -r*t*(r - 1)*(t + 1)/2, lambda r, s, t : -r*(r - 1)*(s - 1)*(s + 1)/2]
ddN19 = [lambda r, s, t : -t*(s - 1)*(s + 1)*(t + 1)/2, lambda r, s, t : -r*t*(r + 1)*(t + 1)/2, lambda r, s, t : -r*(r + 1)*(s - 1)*(s + 1)/2]
ddN20 = [lambda r, s, t : -s*t*(s + 1)*(t + 1)/2, lambda r, s, t : -t*(r - 1)*(r + 1)*(t + 1)/2, lambda r, s, t : -s*(r - 1)*(r + 1)*(s + 1)/2]
ddN21 = [lambda r, s, t : t*(s - 1)*(s + 1)*(t - 1), lambda r, s, t : t*(r - 1)*(r + 1)*(t - 1), lambda r, s, t : (r - 1)*(r + 1)*(s - 1)*(s + 1)]
ddN22 = [lambda r, s, t : s*(s - 1)*(t - 1)*(t + 1), lambda r, s, t : (r - 1)*(r + 1)*(t - 1)*(t + 1), lambda r, s, t : s*(r - 1)*(r + 1)*(s - 1)]
ddN23 = [lambda r, s, t : (s - 1)*(s + 1)*(t - 1)*(t + 1), lambda r, s, t : r*(r - 1)*(t - 1)*(t + 1), lambda r, s, t : r*(r - 1)*(s - 1)*(s + 1)]
ddN24 = [lambda r, s, t : (s - 1)*(s + 1)*(t - 1)*(t + 1), lambda r, s, t : r*(r + 1)*(t - 1)*(t + 1), lambda r, s, t : r*(r + 1)*(s - 1)*(s + 1)]
ddN25 = [lambda r, s, t : s*(s + 1)*(t - 1)*(t + 1), lambda r, s, t : (r - 1)*(r + 1)*(t - 1)*(t + 1), lambda r, s, t : s*(r - 1)*(r + 1)*(s + 1)]
ddN26 = [lambda r, s, t : t*(s - 1)*(s + 1)*(t + 1), lambda r, s, t : t*(r - 1)*(r + 1)*(t + 1), lambda r, s, t : (r - 1)*(r + 1)*(s - 1)*(s + 1)]
ddN27 = [lambda r, s, t : -2*(s - 1)*(s + 1)*(t - 1)*(t + 1), lambda r, s, t : -2*(r - 1)*(r + 1)*(t - 1)*(t + 1), lambda r, s, t : -2*(r - 1)*(r + 1)*(s - 1)*(s + 1)]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN17 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN19 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN20 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN21 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN22 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN23 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN24 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN25 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN26 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN27 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN17 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN19 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN20 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN21 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN22 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN23 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN24 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN25 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN26 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN27 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



=================== PRISM6 ===================
N1 = lambda r, s, t : (t - 1)*(r + s - 1)/2
N2 = lambda r, s, t : -r*(t - 1)/2
N3 = lambda r, s, t : -s*(t - 1)/2
N4 = lambda r, s, t : -(t + 1)*(r + s - 1)/2
N5 = lambda r, s, t : r*(t + 1)/2
N6 = lambda r, s, t : s*(t + 1)/2

dN1 = [lambda r, s, t : t/2 - 1/2, lambda r, s, t : t/2 - 1/2, lambda r, s, t : r/2 + s/2 - 1/2]
dN2 = [lambda r, s, t : 1/2 - t/2, lambda r, s, t : 0, lambda r, s, t : -r/2]
dN3 = [lambda r, s, t : 0, lambda r, s, t : 1/2 - t/2, lambda r, s, t : -s/2]
dN4 = [lambda r, s, t : -t/2 - 1/2, lambda r, s, t : -t/2 - 1/2, lambda r, s, t : -r/2 - s/2 + 1/2]
dN5 = [lambda r, s, t : t/2 + 1/2, lambda r, s, t : 0, lambda r, s, t : r/2]
dN6 = [lambda r, s, t : 0, lambda r, s, t : t/2 + 1/2, lambda r, s, t : s/2]

ddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



=================== PRISM15 ===================
N1 = lambda r, s, t : -(t - 1)*(r + s - 1)*(2*r + 2*s + t)/2
N2 = lambda r, s, t : -r*(t - 1)*(2*r - t - 2)/2
N3 = lambda r, s, t : -s*(t - 1)*(2*s - t - 2)/2
N4 = lambda r, s, t : (t + 1)*(r + s - 1)*(2*r + 2*s - t)/2
N5 = lambda r, s, t : r*(t + 1)*(2*r + t - 2)/2
N6 = lambda r, s, t : s*(t + 1)*(2*s + t - 2)/2
N7 = lambda r, s, t : 2*r*(t - 1)*(r + s - 1)
N8 = lambda r, s, t : 2*s*(t - 1)*(r + s - 1)
N9 = lambda r, s, t : (t - 1)*(t + 1)*(r + s - 1)
N10 = lambda r, s, t : -2*r*s*(t - 1)
N11 = lambda r, s, t : -r*(t - 1)*(t + 1)
N12 = lambda r, s, t : -s*(t - 1)*(t + 1)
N13 = lambda r, s, t : -2*r*(t + 1)*(r + s - 1)
N14 = lambda r, s, t : -2*s*(t + 1)*(r + s - 1)
N15 = lambda r, s, t : 2*r*s*(t + 1)

dN1 = [lambda r, s, t : -(t - 1)*(r + s - 1) - (t - 1)*(2*r + 2*s + t)/2, lambda r, s, t : -(t - 1)*(r + s - 1) - (t - 1)*(2*r + 2*s + t)/2, lambda r, s, t : -(t - 1)*(r + s - 1)/2 - (r + s - 1)*(2*r + 2*s + t)/2]
dN2 = [lambda r, s, t : -r*(t - 1) - (t - 1)*(2*r - t - 2)/2, lambda r, s, t : 0, lambda r, s, t : r*(t - 1)/2 - r*(2*r - t - 2)/2]
dN3 = [lambda r, s, t : 0, lambda r, s, t : -s*(t - 1) - (t - 1)*(2*s - t - 2)/2, lambda r, s, t : s*(t - 1)/2 - s*(2*s - t - 2)/2]
dN4 = [lambda r, s, t : (t + 1)*(r + s - 1) + (t + 1)*(2*r + 2*s - t)/2, lambda r, s, t : (t + 1)*(r + s - 1) + (t + 1)*(2*r + 2*s - t)/2, lambda r, s, t : -(t + 1)*(r + s - 1)/2 + (r + s - 1)*(2*r + 2*s - t)/2]
dN5 = [lambda r, s, t : r*(t + 1) + (t + 1)*(2*r + t - 2)/2, lambda r, s, t : 0, lambda r, s, t : r*(t + 1)/2 + r*(2*r + t - 2)/2]
dN6 = [lambda r, s, t : 0, lambda r, s, t : s*(t + 1) + (t + 1)*(2*s + t - 2)/2, lambda r, s, t : s*(t + 1)/2 + s*(2*s + t - 2)/2]
dN7 = [lambda r, s, t : 2*r*(t - 1) + 2*(t - 1)*(r + s - 1), lambda r, s, t : 2*r*(t - 1), lambda r, s, t : 2*r*(r + s - 1)]
dN8 = [lambda r, s, t : 2*s*(t - 1), lambda r, s, t : 2*s*(t - 1) + 2*(t - 1)*(r + s - 1), lambda r, s, t : 2*s*(r + s - 1)]
dN9 = [lambda r, s, t : (t - 1)*(t + 1), lambda r, s, t : (t - 1)*(t + 1), lambda r, s, t : (t - 1)*(r + s - 1) + (t + 1)*(r + s - 1)]
dN10 = [lambda r, s, t : -2*s*(t - 1), lambda r, s, t : -2*r*(t - 1), lambda r, s, t : -2*r*s]
dN11 = [lambda r, s, t : -(t - 1)*(t + 1), lambda r, s, t : 0, lambda r, s, t : -r*(t - 1) - r*(t + 1)]
dN12 = [lambda r, s, t : 0, lambda r, s, t : -(t - 1)*(t + 1), lambda r, s, t : -s*(t - 1) - s*(t + 1)]
dN13 = [lambda r, s, t : -2*r*(t + 1) - 2*(t + 1)*(r + s - 1), lambda r, s, t : -2*r*(t + 1), lambda r, s, t : -2*r*(r + s - 1)]
dN14 = [lambda r, s, t : -2*s*(t + 1), lambda r, s, t : -2*s*(t + 1) - 2*(t + 1)*(r + s - 1), lambda r, s, t : -2*s*(r + s - 1)]
dN15 = [lambda r, s, t : 2*s*(t + 1), lambda r, s, t : 2*r*(t + 1), lambda r, s, t : 2*r*s]

ddN1 = [lambda r, s, t : 2 - 2*t, lambda r, s, t : 2 - 2*t, lambda r, s, t : -r - s + 1]
ddN2 = [lambda r, s, t : 2 - 2*t, lambda r, s, t : 0, lambda r, s, t : r]
ddN3 = [lambda r, s, t : 0, lambda r, s, t : 2 - 2*t, lambda r, s, t : s]
ddN4 = [lambda r, s, t : 2*t + 2, lambda r, s, t : 2*t + 2, lambda r, s, t : -r - s + 1]
ddN5 = [lambda r, s, t : 2*t + 2, lambda r, s, t : 0, lambda r, s, t : r]
ddN6 = [lambda r, s, t : 0, lambda r, s, t : 2*t + 2, lambda r, s, t : s]
ddN7 = [lambda r, s, t : 4*t - 4, lambda r, s, t : 0, lambda r, s, t : 0]
ddN8 = [lambda r, s, t : 0, lambda r, s, t : 4*t - 4, lambda r, s, t : 0]
ddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 2*r + 2*s - 2]
ddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : -2*r]
ddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : -2*s]
ddN13 = [lambda r, s, t : -4*t - 4, lambda r, s, t : 0, lambda r, s, t : 0]
ddN14 = [lambda r, s, t : 0, lambda r, s, t : -4*t - 4, lambda r, s, t : 0]
ddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]



=================== PRISM18 ===================
N1 = lambda r, s, t : t*(t - 1)*(r + s - 1)*(2*r + 2*s - 1)/2
N2 = lambda r, s, t : r*t*(2*r - 1)*(t - 1)/2
N3 = lambda r, s, t : s*t*(2*s - 1)*(t - 1)/2
N4 = lambda r, s, t : t*(t + 1)*(r + s - 1)*(2*r + 2*s - 1)/2
N5 = lambda r, s, t : r*t*(2*r - 1)*(t + 1)/2
N6 = lambda r, s, t : s*t*(2*s - 1)*(t + 1)/2
N7 = lambda r, s, t : -2*r*t*(t - 1)*(r + s - 1)
N8 = lambda r, s, t : -2*s*t*(t - 1)*(r + s - 1)
N9 = lambda r, s, t : -(t - 1)*(t + 1)*(r + s - 1)*(2*r + 2*s - 1)
N10 = lambda r, s, t : 2*r*s*t*(t - 1)
N11 = lambda r, s, t : -r*(2*r - 1)*(t - 1)*(t + 1)
N12 = lambda r, s, t : -s*(2*s - 1)*(t - 1)*(t + 1)
N13 = lambda r, s, t : -2*r*t*(t + 1)*(r + s - 1)
N14 = lambda r, s, t : -2*s*t*(t + 1)*(r + s - 1)
N15 = lambda r, s, t : 2*r*s*t*(t + 1)
N16 = lambda r, s, t : 4*r*(t - 1)*(t + 1)*(r + s - 1)
N17 = lambda r, s, t : 4*s*(t - 1)*(t + 1)*(r + s - 1)
N18 = lambda r, s, t : -4*r*s*(t - 1)*(t + 1)

dN1 = [lambda r, s, t : t*(t - 1)*(r + s - 1) + t*(t - 1)*(2*r + 2*s - 1)/2, lambda r, s, t : t*(t - 1)*(r + s - 1) + t*(t - 1)*(2*r + 2*s - 1)/2, lambda r, s, t : t*(r + s - 1)*(2*r + 2*s - 1)/2 + (t - 1)*(r + s - 1)*(2*r + 2*s - 1)/2]
dN2 = [lambda r, s, t : r*t*(t - 1) + t*(2*r - 1)*(t - 1)/2, lambda r, s, t : 0, lambda r, s, t : r*t*(2*r - 1)/2 + r*(2*r - 1)*(t - 1)/2]
dN3 = [lambda r, s, t : 0, lambda r, s, t : s*t*(t - 1) + t*(2*s - 1)*(t - 1)/2, lambda r, s, t : s*t*(2*s - 1)/2 + s*(2*s - 1)*(t - 1)/2]
dN4 = [lambda r, s, t : t*(t + 1)*(r + s - 1) + t*(t + 1)*(2*r + 2*s - 1)/2, lambda r, s, t : t*(t + 1)*(r + s - 1) + t*(t + 1)*(2*r + 2*s - 1)/2, lambda r, s, t : t*(r + s - 1)*(2*r + 2*s - 1)/2 + (t + 1)*(r + s - 1)*(2*r + 2*s - 1)/2]
dN5 = [lambda r, s, t : r*t*(t + 1) + t*(2*r - 1)*(t + 1)/2, lambda r, s, t : 0, lambda r, s, t : r*t*(2*r - 1)/2 + r*(2*r - 1)*(t + 1)/2]
dN6 = [lambda r, s, t : 0, lambda r, s, t : s*t*(t + 1) + t*(2*s - 1)*(t + 1)/2, lambda r, s, t : s*t*(2*s - 1)/2 + s*(2*s - 1)*(t + 1)/2]
dN7 = [lambda r, s, t : -2*r*t*(t - 1) - 2*t*(t - 1)*(r + s - 1), lambda r, s, t : -2*r*t*(t - 1), lambda r, s, t : -2*r*t*(r + s - 1) - 2*r*(t - 1)*(r + s - 1)]
dN8 = [lambda r, s, t : -2*s*t*(t - 1), lambda r, s, t : -2*s*t*(t - 1) - 2*t*(t - 1)*(r + s - 1), lambda r, s, t : -2*s*t*(r + s - 1) - 2*s*(t - 1)*(r + s - 1)]
dN9 = [lambda r, s, t : -2*(t - 1)*(t + 1)*(r + s - 1) - (t - 1)*(t + 1)*(2*r + 2*s - 1), lambda r, s, t : -2*(t - 1)*(t + 1)*(r + s - 1) - (t - 1)*(t + 1)*(2*r + 2*s - 1), lambda r, s, t : -(t - 1)*(r + s - 1)*(2*r + 2*s - 1) - (t + 1)*(r + s - 1)*(2*r + 2*s - 1)]
dN10 = [lambda r, s, t : 2*s*t*(t - 1), lambda r, s, t : 2*r*t*(t - 1), lambda r, s, t : 2*r*s*t + 2*r*s*(t - 1)]
dN11 = [lambda r, s, t : -2*r*(t - 1)*(t + 1) - (2*r - 1)*(t - 1)*(t + 1), lambda r, s, t : 0, lambda r, s, t : -r*(2*r - 1)*(t - 1) - r*(2*r - 1)*(t + 1)]
dN12 = [lambda r, s, t : 0, lambda r, s, t : -2*s*(t - 1)*(t + 1) - (2*s - 1)*(t - 1)*(t + 1), lambda r, s, t : -s*(2*s - 1)*(t - 1) - s*(2*s - 1)*(t + 1)]
dN13 = [lambda r, s, t : -2*r*t*(t + 1) - 2*t*(t + 1)*(r + s - 1), lambda r, s, t : -2*r*t*(t + 1), lambda r, s, t : -2*r*t*(r + s - 1) - 2*r*(t + 1)*(r + s - 1)]
dN14 = [lambda r, s, t : -2*s*t*(t + 1), lambda r, s, t : -2*s*t*(t + 1) - 2*t*(t + 1)*(r + s - 1), lambda r, s, t : -2*s*t*(r + s - 1) - 2*s*(t + 1)*(r + s - 1)]
dN15 = [lambda r, s, t : 2*s*t*(t + 1), lambda r, s, t : 2*r*t*(t + 1), lambda r, s, t : 2*r*s*t + 2*r*s*(t + 1)]
dN16 = [lambda r, s, t : 4*r*(t - 1)*(t + 1) + 4*(t - 1)*(t + 1)*(r + s - 1), lambda r, s, t : 4*r*(t - 1)*(t + 1), lambda r, s, t : 4*r*(t - 1)*(r + s - 1) + 4*r*(t + 1)*(r + s - 1)]
dN17 = [lambda r, s, t : 4*s*(t - 1)*(t + 1), lambda r, s, t : 4*s*(t - 1)*(t + 1) + 4*(t - 1)*(t + 1)*(r + s - 1), lambda r, s, t : 4*s*(t - 1)*(r + s - 1) + 4*s*(t + 1)*(r + s - 1)]
dN18 = [lambda r, s, t : -4*s*(t - 1)*(t + 1), lambda r, s, t : -4*r*(t - 1)*(t + 1), lambda r, s, t : -4*r*s*(t - 1) - 4*r*s*(t + 1)]

ddN1 = [lambda r, s, t : 2*t*(t - 1), lambda r, s, t : 2*t*(t - 1), lambda r, s, t : (r + s - 1)*(2*r + 2*s - 1)]
ddN2 = [lambda r, s, t : 2*t*(t - 1), lambda r, s, t : 0, lambda r, s, t : r*(2*r - 1)]
ddN3 = [lambda r, s, t : 0, lambda r, s, t : 2*t*(t - 1), lambda r, s, t : s*(2*s - 1)]
ddN4 = [lambda r, s, t : 2*t*(t + 1), lambda r, s, t : 2*t*(t + 1), lambda r, s, t : (r + s - 1)*(2*r + 2*s - 1)]
ddN5 = [lambda r, s, t : 2*t*(t + 1), lambda r, s, t : 0, lambda r, s, t : r*(2*r - 1)]
ddN6 = [lambda r, s, t : 0, lambda r, s, t : 2*t*(t + 1), lambda r, s, t : s*(2*s - 1)]
ddN7 = [lambda r, s, t : -4*t*(t - 1), lambda r, s, t : 0, lambda r, s, t : -4*r*(r + s - 1)]
ddN8 = [lambda r, s, t : 0, lambda r, s, t : -4*t*(t - 1), lambda r, s, t : -4*s*(r + s - 1)]
ddN9 = [lambda r, s, t : -4*(t - 1)*(t + 1), lambda r, s, t : -4*(t - 1)*(t + 1), lambda r, s, t : -2*(r + s - 1)*(2*r + 2*s - 1)]
ddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 4*r*s]
ddN11 = [lambda r, s, t : -4*(t - 1)*(t + 1), lambda r, s, t : 0, lambda r, s, t : -2*r*(2*r - 1)]
ddN12 = [lambda r, s, t : 0, lambda r, s, t : -4*(t - 1)*(t + 1), lambda r, s, t : -2*s*(2*s - 1)]
ddN13 = [lambda r, s, t : -4*t*(t + 1), lambda r, s, t : 0, lambda r, s, t : -4*r*(r + s - 1)]
ddN14 = [lambda r, s, t : 0, lambda r, s, t : -4*t*(t + 1), lambda r, s, t : -4*s*(r + s - 1)]
ddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 4*r*s]
ddN16 = [lambda r, s, t : 8*(t - 1)*(t + 1), lambda r, s, t : 0, lambda r, s, t : 8*r*(r + s - 1)]
ddN17 = [lambda r, s, t : 0, lambda r, s, t : 8*(t - 1)*(t + 1), lambda r, s, t : 8*s*(r + s - 1)]
ddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : -8*r*s]

dddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN17 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
dddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

ddddN1 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN2 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN3 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN4 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN5 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN6 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN7 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN8 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN9 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN10 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN11 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN12 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN13 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN14 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN15 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN16 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN17 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]
ddddN18 = [lambda r, s, t : 0, lambda r, s, t : 0, lambda r, s, t : 0]

import numpy as np
import matplotlib.pyplot as plt
from EasyFEA import Display

try:
    import sympy
except ModuleNotFoundError:
    raise Exception("sympy must be installed!")

Display.Clear()

# SEGMENTS
# TRIANGLES
# QUADRANGLES
# TETRAHEDRON
# HEXAHEDRON
# PRISM

# ----------------------------------------------
# Options
# ----------------------------------------------

plot_dN = True
plot_ddN = True
plot_dddN = True
plot_ddddN = True

# coords = sympy.symbols("x, y, z")
coords = sympy.symbols("r, s, t")

# ----------------------------------------------
# Public functions
# ----------------------------------------------


def Compute_and_Print(
    polynom, *args, useSimplify=True, useFactor=True, useEulerBernoulli=False
):
    """Compute and print shape functions and their derivatives for a given polynom.

    Parameters
    ----------
    polynom : function
        Polynomial basis function taking 1, 2, or 3 arguments.
    args : tuple
        Coordinates (1, 2, or 3 lists of x, y, z).
    useSimplify : bool, optional
        Simplify the shape functions. Default is True.
    useFactor : bool, optional
        Factor the shape functions. Default is True.
    useEulerBernoulli : bool, optional
        useEulerBernoulli shape functions. Default is False.
    """

    local_coords, dim = __Get_local_coords_and_dim(*args)

    shape_functions = __Get_shape_functions(
        polynom, local_coords, dim, useSimplify, useFactor, useEulerBernoulli
    )

    __Print_functions(shape_functions, dim, "N")

    # derivative_shape_functions
    if plot_dN:
        dN_functions = __Get_derivative_functions(shape_functions, dim, 1)
        __Print_functions(dN_functions, dim, "dN")

    if plot_ddN:
        ddN_functions = __Get_derivative_functions(shape_functions, dim, 2)
        __Print_functions(ddN_functions, dim, "ddN")

    if plot_dddN:
        dddN_functions = __Get_derivative_functions(shape_functions, dim, 3)
        __Print_functions(dddN_functions, dim, "dddN")

    if plot_ddddN:
        ddddN_functions = __Get_derivative_functions(shape_functions, dim, 4)
        __Print_functions(ddddN_functions, dim, "ddddN")


def Plot_Nodes(title: str, *args):
    local_coords, dim = __Get_local_coords_and_dim(*args)
    nPe = local_coords.shape[0]

    list_x, list_y, list_z = local_coords.T

    if dim == 3:
        ax = Display.Init_Axes(3, elev=16, azim=37)
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_zlabel("z")
        ax.set_title(title)
        ax.axis("equal")

        ax.scatter(list_x, list_y, list_z)
        [ax.text(list_x[i], list_y[i], list_z[i], i) for i in range(nPe)]
    else:
        ax = Display.Init_Axes(2)
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_title(title)
        ax.axis("equal")
        ax.grid(zorder=-10)

        ax.scatter(list_x, list_y)
        [ax.text(list_x[i], list_y[i], i + 1) for i in range(nPe)]


# ----------------------------------------------
# Private functions do not touch
# ----------------------------------------------


def __Get_local_coords_and_dim(*args):
    # Get dim
    dim = len(args)
    assert dim in [1, 2, 3], "The number of lists in args must be 1, 2, or 3."

    # Get coordinates
    list_x = args[0]
    nPe = len(list_x)

    list_y = args[1] if dim > 1 else [0] * nPe
    assert (
        len(list_y) == nPe
    ), "The length of list_y must be equal to the length of list_x."

    list_z = args[2] if dim > 2 else [0] * nPe
    assert (
        len(list_z) == nPe
    ), "The length of list_z must be equal to the length of list_x."

    local_coords = np.array([list_x, list_y, list_z]).T

    return local_coords, dim


def __Get_shape_functions(
    polynom,
    local_coords: np.ndarray,
    dim: int,
    useSimplify=True,
    useFactor=True,
    useEulerBernoulli=False,
) -> list:
    nPe = local_coords.shape[0]

    nF = __Get_functions_per_node(polynom, local_coords, dim)

    if useEulerBernoulli:
        assert nF == 2, "euler bernoulli shape functions use 2 functions per node!"

    # construct matrix a
    matrix_A = __Get_matrix_A(polynom, local_coords, dim)

    # Get symbols and coords symbols
    symbols = sympy.symbols(f"x0:{nPe * nF}")

    functions = []

    for i in range(nPe * nF):
        # construct vector b
        vector_b = np.zeros(nPe * nF)
        if useEulerBernoulli:
            vector_b[i] = 1 if i % 2 == 0 else 1 / 2
            coefs = polynom(*coords[:dim])[0]
        else:
            vector_b[i] = 1
            coefs = polynom(*coords[:dim])

        # solve x from A x = b
        vector_x = np.linalg.solve(matrix_A, vector_b)
        # check that A x = b
        assert np.linalg.norm(matrix_A @ vector_x - vector_b) <= 1e-12

        # construct shape function
        function = sum(coef * term for coef, term in zip(symbols, coefs))
        # apply values
        function = function.subs({key: value for key, value in zip(symbols, vector_x)})

        # apply function display properties
        function = __chop(function)
        if useSimplify:
            function = function.nsimplify()
        if useFactor:
            function = function.factor()

        functions.append(function)

    return functions


def __Get_matrix_A(polynom, local_coords: np.ndarray, dim: int):
    nPe = local_coords.shape[0]

    nF = __Get_functions_per_node(polynom, local_coords, dim)
    indexes = np.arange(nPe * nF)
    if nF > 1:
        indexes = indexes.reshape(-1, 2)

    matrix_A = np.zeros((nPe * nF, nPe * nF))

    for n in range(nPe):
        matrix_A[indexes[n], :] = polynom(*local_coords[n, :dim])

    return matrix_A


def __Get_functions_per_node(polynom, local_coords: np.ndarray, dim: int) -> int:
    eval = np.array(polynom(*local_coords[0, :dim]))
    if len(eval.shape) == 2:
        return eval.shape[0]
    elif len(eval.shape) == 1:
        return 1
    else:
        raise Exception("polynom must be a function or a list of functions.")


def __Get_derivative_functions(functions, dim, order):
    assert isinstance(functions, list), "functions must be a list"
    assert dim in [1, 2, 3]
    assert order >= 1 and isinstance(order, int), "order must be >= 1"

    derivative_functions = []

    for function in functions:
        functions_per_dim = []

        # loop on dimensions
        for coord in coords[:dim]:
            func = function  # copy of `function`
            # loop to derive the func function
            for _ in range(order):
                func = func.diff(coord)
            functions_per_dim.append(func)

        derivative_functions.append(functions_per_dim)

    return derivative_functions


def __Print_functions(functions: list, dim: int, name="", printArray=False):
    # lamba string (e.g. lamda r, s)
    lambda_str = f"lambda {', '.join(str(coord) for coord in coords[:dim])}"

    # print each functions
    for i, function in enumerate(functions):
        if isinstance(function, list):
            print(
                f"{name}{i + 1} = [{', '.join(f'{lambda_str} : {func}' for func in function)}]"
            )
        else:
            print(f"{name}{i + 1} = {lambda_str} : {function}")

    print()
    if printArray:
        end = ".reshape(-1, 1)" if len(np.shape(functions)) == 1 else ""
        nF = len(functions)
        print(
            f"{name} = np.array([{', '.join(f'{name}{i + 1}' for i in range(nF))}]){end}\n"
        )


def __chop(expr):
    """Recursive function that replaces small values in the sympy expression with zero."""
    if isinstance(expr, sympy.Float) and abs(expr) < 1e-12:
        return sympy.S.Zero
    elif isinstance(expr, sympy.Add):
        return sympy.Add(*[__chop(arg) for arg in expr.args])
    elif isinstance(expr, sympy.Mul):
        return sympy.Mul(*[__chop(arg) for arg in expr.args])
    elif isinstance(expr, sympy.Pow):
        return sympy.Pow(__chop(expr.base), expr.exp)
    else:
        return expr


# ----------------------------------------------
# SEGMENTS
# ----------------------------------------------


def Do_Segments():
    # ----------------------------------------------
    # SEG 2
    # ----------------------------------------------
    #       v
    #       ^
    #       |
    #       |
    #  0----+----1 --> u
    # ----------------------------------------------

    name = "SEG2"
    Display.Section(name)

    list_x = [-1, 1]
    Plot_Nodes(name, list_x)

    polynom = lambda x: [x, 1]
    Compute_and_Print(polynom, list_x)

    Display.Section(name + " EulerBernoulli")
    polynom = lambda x: [[1, x, x**2, x**3], [0, 1, 2 * x, 3 * x**2]]
    Compute_and_Print(polynom, list_x, useEulerBernoulli=True)

    # ----------------------------------------------
    # SEG 3
    # ----------------------------------------------
    #       v
    #       ^
    #       |
    #       |
    #  0----2----1 --> u
    # ----------------------------------------------

    name = "SEG3"
    Display.Section(name)

    list_x = [-1, 1, 0]
    Plot_Nodes(name, list_x)

    polynom = lambda x: [x**2, x, 1]
    Compute_and_Print(polynom, list_x)

    Display.Section(name + " EulerBernoulli")
    polynom = lambda x: [
        [1, x, x**2, x**3, x**4, x**5],
        [0, 1, 2 * x, 3 * x**2, 4 * x**3, 5 * x**4],
    ]
    Compute_and_Print(polynom, list_x, useEulerBernoulli=True)

    # ----------------------------------------------
    # SEG 4
    # ----------------------------------------------
    #        v
    #        ^
    #        |
    #        |
    #  0---2-+-3---1 --> u
    # ----------------------------------------------

    name = "SEG4"
    Display.Section(name)

    list_x = [-1, 1, -1 / 3, 1 / 3]
    Plot_Nodes(name, list_x)

    polynom = lambda x: [x**3, x**2, x, 1]
    Compute_and_Print(polynom, list_x, useFactor=False)

    Display.Section(name + " EulerBernoulli")
    polynom = lambda x: [
        [1, x, x**2, x**3, x**4, x**5, x**6, x**7],
        [0, 1, 2 * x, 3 * x**2, 4 * x**3, 5 * x**4, 6 * x**5, 7 * x**6],
    ]
    Compute_and_Print(polynom, list_x, useEulerBernoulli=True, useFactor=False)

    # ----------------------------------------------
    # SEG 5
    # ----------------------------------------------
    #        v
    #        ^
    #        |
    #        |
    #  0--2--3--4--1 --> u
    # ----------------------------------------------

    name = "SEG5"
    Display.Section(name)

    list_x = [-1, 1, -1 / 2, 0, 1 / 2]
    Plot_Nodes(name, list_x)

    polynom = lambda x: [x**4, x**3, x**2, x, 1]

    Compute_and_Print(polynom, list_x, useFactor=False)

    Display.Section(name + " EulerBernoulli")
    polynom = lambda x: [
        [1, x, x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9],
        [
            0,
            1,
            2 * x,
            3 * x**2,
            4 * x**3,
            5 * x**4,
            6 * x**5,
            7 * x**6,
            8 * x**7,
            9 * x**8,
        ],
    ]
    Compute_and_Print(polynom, list_x, useEulerBernoulli=True, useFactor=False)


# ----------------------------------------------
# TRIANGLES
# ----------------------------------------------


def Do_Triangles():
    # ----------------------------------------------
    # TRI3
    # ----------------------------------------------
    # v
    # ^
    # |
    # 2
    # |`\
    # |  `\
    # |    `\
    # |      `\
    # |        `\
    # 0----------1 --> u
    # ----------------------------------------------

    name = "TRI3"
    Display.Section(name)

    list_x = [0, 1, 0]
    list_y = [0, 0, 1]
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [x, y, 1]

    Compute_and_Print(polynom, list_x, list_y)

    # ----------------------------------------------
    # TRI6
    # ----------------------------------------------
    # v
    # ^
    # |
    # 2
    # |`\
    # |  `\
    # 5    `4
    # |      `\
    # |        `\
    # 0----3-----1 --> u
    # ----------------------------------------------

    name = "TRI6"
    Display.Section(name)

    list_x = [0, 1, 0, 0.5, 0.5, 0]
    list_y = [0, 0, 1, 0, 0.5, 0.5]
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [x**2, y**2, x * y, x, y, 1]

    Compute_and_Print(polynom, list_x, list_y)

    # ----------------------------------------------
    # TRI10
    # ----------------------------------------------
    # v
    # ^
    # |
    # 2
    # | \
    # 7   6
    # |     \
    # 8  (9)  5
    # |         \
    # 0---3---4---1 --> u
    # ----------------------------------------------

    name = "TRI10"
    Display.Section(name)

    # fmt: off
    list_x = [
        0,1,0,
        1/3,2/3,
        2/3, 1/3,
        0,0,
        1/3
    ]
    list_y = [
        0,0,1,
        0,0,
        1/3, 2/3,
        2/3,1/3,
        1/3
    ]
    # fmt: off
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [x**3, y**3, x**2 * y, x * y**2, x**2, y**2, x * y, x, y, 1]

    Compute_and_Print(polynom, list_x, list_y, useFactor=False)

    # ----------------------------------------------
    # TRI15
    # ----------------------------------------------
    #
    # 2
    # | \
    # 9   8
    # |     \
    # 10 (14)  7
    # |         \
    # 11 (12) (13) 6
    # |             \
    # 0---3---4---5---1
    # ----------------------------------------------

    name = "TRI15"
    Display.Section(name)

    # fmt: off
    list_x = [
        0,1,0,
        1/4,1/2,3/4,
        3/4,1/2,1/4,
        0,0,0,
        1/4,1/2,1/4
    ]
    list_y = [
        0,0,1,
        0,0,0,
        1/4,1/2,3/4,
        3/4,1/2,1/4,
        1/4,1/4,1/2
    ]
    # fmt: on
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [
        x**4,
        x**3 * y,
        x**2 * y**2,
        x * y**3,
        y**4,
        x**3,
        x**2 * y,
        x * y**2,
        y**3,
        x**2,
        x * y,
        y**2,
        x,
        y,
        1,
    ]

    Compute_and_Print(polynom, list_x, list_y, useFactor=False)


# ----------------------------------------------
# QUADRANGLES
# ----------------------------------------------


def Do_Quadrangles():
    # ----------------------------------------------
    # QUAD4
    # ----------------------------------------------
    #       v
    #       ^
    #       |
    # 3-----------2
    # |     |     |
    # |     |     |
    # |     +---- | --> u
    # |           |
    # |           |
    # 0-----------1
    # ----------------------------------------------

    name = "QUAD4"
    Display.Section(name)

    list_x = [-1, 1, 1, -1]
    list_y = [-1, -1, 1, 1]
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [x * y, x, y, 1]

    Compute_and_Print(polynom, list_x, list_y)

    # ----------------------------------------------
    # QUAD8
    # ----------------------------------------------
    #       v
    #       ^
    #       |
    # 3-----6-----2
    # |     |     |
    # |     |     |
    # 7     +---- 5 --> u
    # |           |
    # |           |
    # 0-----4-----1
    # ----------------------------------------------

    name = "QUAD8"
    Display.Section(name)

    list_x = [-1, 1, 1, -1, 0, 1, 0, -1]
    list_y = [-1, -1, 1, 1, -1, 0, 1, 0]
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [x**2 * y, y**2 * x, x**2, y**2, x * y, x, y, 1]

    Compute_and_Print(polynom, list_x, list_y, useFactor=True)

    # ----------------------------------------------
    # QUAD9
    # ----------------------------------------------
    #       v
    #       ^
    #       |
    # 3-----6-----2
    # |     |     |
    # |     |     |
    # 7     8---- 5 --> u
    # |           |
    # |           |
    # 0-----4-----1
    # ----------------------------------------------

    name = "QUAD9"
    Display.Section(name)

    list_x = [-1, 1, 1, -1, 0, 1, 0, -1, 0]
    list_y = [-1, -1, 1, 1, -1, 0, 1, 0, 0]
    Plot_Nodes(name, list_x, list_y)

    polynom = lambda x, y: [x**2 * y**2, x**2 * y, y**2 * x, x**2, y**2, x * y, x, y, 1]

    Compute_and_Print(polynom, list_x, list_y, useFactor=True)


# ----------------------------------------------
# TETRAHEDRON
# ----------------------------------------------


def Do_Tetrahedron():
    # ----------------------------------------------
    # TETRA4
    # ----------------------------------------------
    #                    v
    #                  .
    #                ,/
    #               /
    #            2
    #          ,/|`\
    #        ,/  |  `\
    #      ,/    '.   `\
    #    ,/       |     `\
    #  ,/         |       `\
    # 0-----------'.--------1 --> u
    #  `\.         |      ,/
    #     `\.      |    ,/
    #        `\.   '. ,/
    #           `\. |/
    #              `3
    #                 `\.
    #                    ` w
    # ----------------------------------------------

    name = "TETRA4"
    Display.Section(name)

    list_x = [0, 1, 0, 0]
    list_y = [0, 0, 1, 0]
    list_z = [0, 0, 0, 1]
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [x, y, z, 1]

    Compute_and_Print(polynom, list_x, list_y, list_z)

    # ----------------------------------------------
    # TETRA10
    # ----------------------------------------------
    #                    v
    #                  .
    #                ,/
    #               /
    #            2
    #          ,/|`\
    #        ,/  |  `\
    #      ,6    '.   `5
    #    ,/       8     `\
    #  ,/         |       `\
    # 0--------4--'.--------1 --> u
    #  `\.         |      ,/
    #     `\.      |    ,9
    #        `7.   '. ,/
    #           `\. |/
    #              `3
    #                 `\.
    #                    ` w
    # ----------------------------------------------

    name = "TETRA10"
    Display.Section(name)

    list_x = [0, 1, 0, 0, 0.5, 0.5, 0, 0, 0, 0.5]
    list_y = [0, 0, 1, 0, 0, 0.5, 0.5, 0, 0.5, 0]
    list_z = [0, 0, 0, 1, 0, 0, 0, 0.5, 0.5, 0.5]
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [x**2, y**2, z**2, x * y, x * z, y * z, x, y, z, 1]

    Compute_and_Print(polynom, list_x, list_y, list_z)


# ----------------------------------------------
# HEXAHEDRON
# ----------------------------------------------


def Do_Hexahedron():
    # ----------------------------------------------
    # HEXA8
    # ----------------------------------------------
    #        v
    # 3----------2
    # |\     ^   |\
    # | \    |   | \
    # |  \   |   |  \
    # |   7------+---6
    # |   |  +-- |-- | -> u
    # 0---+---\--1   |
    #  \  |    \  \  |
    #   \ |     \  \ |
    #    \|      w  \|
    #     4----------5
    # ----------------------------------------------

    name = "HEXA8"
    Display.Section(name)

    list_x = [-1, 1, 1, -1, -1, 1, 1, -1]
    list_y = [-1, -1, 1, 1, -1, -1, 1, 1]
    list_z = [-1, -1, -1, -1, 1, 1, 1, 1]
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [x * y * z, x * y, x * z, y * z, x, y, z, 1]

    Compute_and_Print(polynom, list_x, list_y, list_z)

    # ----------------------------------------------
    # HEXA20
    # ----------------------------------------------
    #        v
    # 3----13----2
    # |\     ^   |\
    # | 15   |   | 14
    # 9  \   |   11 \
    # |   7----19+---6
    # |   |  +-- |-- | -> u
    # 0---+-8-\--1   |
    #  \  17   \  \  18
    #  10 |     \  12|
    #    \|      w  \|
    #     4----16----5
    # ----------------------------------------------

    name = "HEXA20"
    Display.Section(name)

    # fmt: off
    list_x = [-1,1,1,-1,
            -1,1,1,-1,
            0,-1,-1,1,
            1,0,1,-1,
            0,-1,1,0]
    list_y = [-1,-1,1,1,
            -1,-1,1,1,
            -1,0,-1,0,
            -1,1,1,1,
            -1,0,0,1]
    list_z = [-1,-1,-1,-1,
            1,1,1,1,
            -1,-1,0,-1,
            0,-1,0,0,
            1,1,1,1]
    # fmt: on
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [
        x**2 * y * z,
        y**2 * x * z,
        z**2 * x * y,
        x**2 * y,
        y**2 * x,
        z**2 * x,
        x**2 * z,
        y**2 * z,
        z**2 * y,
        x**2,
        y**2,
        z**2,
        x * y,
        x * z,
        y * z,
        x,
        y,
        z,
        x * y * z,
        1,
    ]

    Compute_and_Print(polynom, list_x, list_y, list_z, useSimplify=True, useFactor=True)

    # ----------------------------------------------
    # HEXA27
    # ----------------------------------------------
    #
    # 3----13----2
    # |\         |\
    # |15    24  | 14
    # 9  \ 20    11 \
    # |   7----19+---6
    # |22 |  26  | 23|
    # 0---+-8----1   |
    #  \ 17    25 \  18
    #  10 |  21    12|
    #    \|         \|
    #     4----16----5
    # ----------------------------------------------

    name = "HEXA27"
    Display.Section(name)

    # fmt: off
    list_x = [
        -1,1,1,-1,
        -1,1,1,-1,
        0,-1,-1,1,
        1,0,1,-1,
        0,-1,1,0,
        0,0,-1,1,
        0,0,0
    ]
    list_y = [
        -1,-1,1,1,
        -1,-1,1,1,
        -1,0,-1,0,
        -1,1,1,1,
        -1,0,0,1,
        0,-1,0,0,
        1,0,0
    ]
    list_z = [
        -1,-1,-1,-1,
        1,1,1,1,
        -1,-1,0,-1,
        0,-1,0,0,
        1,1,1,1,
        -1,0,0,0,
        0,1,0
    ]
    # fmt: on
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [
        x**2 * z**2 * y,
        x**2 * y**2 * z,
        y**2 * z**2 * x,
        x**2 * z**2 * y**2,
        x**2 * y**2,
        x**2 * z**2,
        y**2 * z**2,
        x**2 * y * z,
        y**2 * x * z,
        z**2 * x * y,
        x**2 * y,
        y**2 * x,
        z**2 * x,
        x**2 * z,
        y**2 * z,
        z**2 * y,
        x**2,
        y**2,
        z**2,
        x * y,
        x * z,
        y * z,
        x,
        y,
        z,
        x * y * z,
        1,
    ]

    Compute_and_Print(polynom, list_x, list_y, list_z, useSimplify=True, useFactor=True)


# ----------------------------------------------
# PRISM
# ----------------------------------------------


def Do_Prism():
    # ----------------------------------------------
    # PRISM6
    # ----------------------------------------------
    #            w
    #            ^
    #            |
    #            3
    #          ,/|`\
    #        ,/  |  `\
    #      ,/    |    `\
    #     4------+------5
    #     |      |      |
    #     |    ,/|`\    |
    #     |  ,/  |  `\  |
    #     |,/    |    `\|
    #    ,|      |      |\
    #  ,/ |      0      | `\
    # u   |    ,/ `\    |    v
    #     |  ,/     `\  |
    #     |,/         `\|
    #     1-------------2
    # ----------------------------------------------

    name = "PRISM6"
    Display.Section(name)

    list_x = [0, 1, 0, 0, 1, 0]
    list_y = [0, 0, 1, 0, 0, 1]
    list_z = [-1, -1, -1, 1, 1, 1]
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [x * z, y * z, x, y, z, 1]

    Compute_and_Print(polynom, list_x, list_y, list_z)

    # ----------------------------------------------
    # PRISM15
    # ----------------------------------------------
    #            w
    #            ^
    #            |
    #            3
    #          ,/|`\
    #        12  |  13
    #      ,/    |    `\
    #     4------14-----5
    #     |      8      |
    #     |    ,/|`\    |
    #     |  ,/  |  `\  |
    #     |,/    |    `\|
    #    ,10     |      11
    #  ,/ |      0      | \
    # u   |    ,/ `\    |   v
    #     |  ,6     `7  |
    #     |,/         `\|
    #     1------9------2
    # ----------------------------------------------

    name = "PRISM15"
    Display.Section(name)

    list_x = [0, 1, 0, 0, 1, 0, 0.5, 0, 0, 0.5, 1, 0, 0.5, 0, 0.5]
    list_y = [0, 0, 1, 0, 0, 1, 0, 0.5, 0, 0.5, 0, 1, 0, 0.5, 0.5]
    list_z = [-1, -1, -1, 1, 1, 1, -1, -1, 0, -1, 0, 0, 1, 1, 1]
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [
        x**2 * z,
        y**2 * z,
        z**2 * x,
        z**2 * y,
        x * y * z,
        x**2,
        y**2,
        z**2,
        x * y,
        x * z,
        y * z,
        x,
        y,
        z,
        1,
    ]

    Compute_and_Print(polynom, list_x, list_y, list_z)

    # ----------------------------------------------
    # PRISM18
    # ----------------------------------------------
    #            w
    #            ^
    #            |
    #            3
    #          ,/|`\
    #        12  |  13
    #      ,/    |    `\
    #     4------14-----5
    #     |      8      |
    #     |    ,/|`\    |
    #     |  15  |  16  |
    #     |,/    |    `\|
    #    ,10-----17-----11
    #  ,/ |      0      | \
    # u   |    ,/ `\    |   v
    #     |  ,6     `7  |
    #     |,/         `\|
    #     1------9------2
    # ----------------------------------------------

    name = "PRISM18"
    Display.Section(name)

    list_x = [0, 1, 0, 0, 1, 0, 0.5, 0, 0, 0.5, 1, 0, 0.5, 0, 0.5, 0.5, 0, 0.5]
    list_y = [0, 0, 1, 0, 0, 1, 0, 0.5, 0, 0.5, 0, 1, 0, 0.5, 0.5, 0, 0.5, 0.5]
    list_z = [-1, -1, -1, 1, 1, 1, -1, -1, 0, -1, 0, 0, 1, 1, 1, 0, 0, 0]
    Plot_Nodes(name, list_x, list_y, list_z)

    polynom = lambda x, y, z: [
        x**2,
        y**2,
        x * y,
        x,
        y,
        1,
        x**2 * z,
        y**2 * z,
        x * y * z,
        x * z,
        y * z,
        z,
        x**2 * z**2,
        y**2 * z**2,
        x * y * z**2,
        x * z**2,
        y * z**2,
        z**2,
    ]

    Compute_and_Print(polynom, list_x, list_y, list_z)


# ----------------------------------------------
# MAIN
# ----------------------------------------------

if __name__ == "__main__":
    Do_Segments()

    Do_Triangles()

    Do_Quadrangles()

    Do_Tetrahedron()

    Do_Hexahedron()

    Do_Prism()

    plt.show()