In [ ]:
from sympy import sin, cos, tan, exp, sqrt, init_printing, symbols, fcode, Eq
from sympy.utilities.codegen import codegen
from sympy.matrices import *
from __future__ import print_function
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
init_printing()
In [ ]:
# Define symbols
x1,y1,z1 = symbols('x1 y1 z1', real=True)
x2,y2,z2 = symbols('x2 y2 z2', real=True)
x3,y3,z3 = symbols('x3 y3 z3', real=True)
x4,y4,z4 = symbols('x4 y4 z4', real=True)
In [ ]:
# B element (vector)
r1 = Matrix([x1,y1,z1])
r2 = Matrix([x2,y2,z2])
r12 = r2 - r1
e12 = r12/r12.norm()
s1 = -e12
s2 = e12
Bstre = Matrix([s1,s2])
In [ ]:
# B derivative (matrix)
derBstre = Matrix([Bstre.diff(x1),Bstre.diff(y1),Bstre.diff(z1),
Bstre.diff(x2),Bstre.diff(y2),Bstre.diff(z2)]).reshape(6,6)
In [ ]:
# Generate subroutine and print
Bder = MatrixSymbol('Bder',6,6)
subroutine = codegen(('derBstre', Eq(Bder,derBstre)), 'f95',None,'internal')
print(subroutine[0][1])
#f = open(subroutine[0][0],'w')
#f.write(subroutine[0][1])
#f.close()
In [ ]:
# B element (vector)
r1 = Matrix([x1,y1,z1])
r2 = Matrix([x2,y2,z2])
r3 = Matrix([x3,y3,z3])
r31 = r1 - r3
r32 = r2 - r3
e31 = r31/r31.norm()
e32 = r32/r32.norm()
sinphi = e31.cross(e32).norm()
cosphi = e31.dot(e32)
r31 = r31.norm()
r32 = r32.norm()
s1 = (cosphi*e31-e32)/r31/sinphi
s2 = (cosphi*e32-e31)/r32/sinphi
s3 = ((r31-r32*cosphi)*e31+(r32-r31*cosphi)*e32)/r31/r32/sinphi
Bbend = Matrix([s1,s2,s3])
In [ ]:
# B derivative (matrix)
derBbend = Matrix([Bbend.diff(x1),Bbend.diff(y1),Bbend.diff(z1),
Bbend.diff(x2),Bbend.diff(y2),Bbend.diff(z2),
Bbend.diff(x3),Bbend.diff(y3),Bbend.diff(z3)]).reshape(9,9)
In [ ]:
# Generate subroutine and print
Bder = MatrixSymbol('Bder',9,9)
subroutine = codegen(('derBbend', Eq(Bder,derBbend)), 'f95',None,'internal')
f = open(subroutine[0][0],'w')
f.write(subroutine[0][1])
f.close()
In [279]:
# B element (vector)
r1 = Matrix([x1,y1,z1])
r2 = Matrix([x2,y2,z2])
r3 = Matrix([x3,y3,z3])
r4 = Matrix([x4,y4,z4])
r12 = r2 - r1
r21 = r1 - r2
r23 = r3 - r2
r32 = r2 - r3
r43 = r3 - r4
r34 = r4 - r3
e12 = r12/r12.norm()
e21 = r21/r21.norm()
e23 = r23/r23.norm()
e32 = r32/r32.norm()
e43 = r43/r43.norm()
e34 = r34/r34.norm()
sinphi2 = e21.cross(e23).norm()
sinphi3 = e32.cross(e34).norm()
cosphi2 = e21.dot(e23)
cosphi3 = e32.dot(e34)
r12 = r12.norm()
r23 = r23.norm()
r43 = r43.norm()
#
s1 = -e12.cross(e23)/r12/sinphi2**2
s2 = (r23-r12*cosphi2)*e12.cross(e23)/r23/r12/sinphi2**2 + cosphi3*e43.cross(e32)/r23/sinphi3**2
s3 = (r23-r43*cosphi3)*e43.cross(e32)/r23/r43/sinphi3**2 + cosphi2*e12.cross(e23)/r23/sinphi2**2
s4 = -e43.cross(e32)/r43/sinphi3**2
Bdihe = Matrix([s1,s2,s3,s4])
In [280]:
# B derivative (matrix)
derBdihe = Matrix([Bdihe.diff(x1),Bdihe.diff(y1),Bdihe.diff(z1),
Bdihe.diff(x2),Bdihe.diff(y2),Bdihe.diff(z2),
Bdihe.diff(x3),Bdihe.diff(y3),Bdihe.diff(z3),
Bdihe.diff(x4),Bdihe.diff(y4),Bdihe.diff(z4)]).reshape(12,12)
In [281]:
# Generate subroutine and print
Bder = MatrixSymbol('Bder',12,12)
subroutine = codegen(('derBdihe', Eq(Bder,derBdihe)), 'f95',None,'internal')
f = open(subroutine[0][0],'w')
f.write(subroutine[0][1])
f.close()
In [ ]: