00001
00002 from swiginac import *
00003 from SyFi import *
00004 from sfc.integral_tools import geometry_mapping
00005 from sfc.symbolic_utils import symbol, symbols
00006
00007 initSyFi(2)
00008 x =cvar.x; y =cvar.y
00009
00010
00011 triangle = Triangle([1.3, 1.3], [2.4,2.1], [1.4, 2.8])
00012 G, x0 = geometry_mapping(triangle)
00013
00014 print G*matrix(2,1, [0,0]) + x0
00015 print G*matrix(2,1, [1,0]) + x0
00016 print G*matrix(2,1, [0,1]) + x0
00017
00018 Ginv = G.inverse()
00019 print Ginv*(matrix(2,1, [1.3,1.3]) - x0)
00020 print Ginv*(matrix(2,1, [2.4,2.1]) - x0)
00021 print Ginv*(matrix(2,1, [1.4,2.8]) - x0)
00022
00023
00024 fe = Lagrange(triangle)
00025
00026 print "global fe.N(0) ", fe.N(0)
00027
00028 reference_triangle = ReferenceTriangle()
00029 reference_fe = Lagrange(reference_triangle)
00030 xi0, xi1 = symbol("xi0"), symbol("xi1")
00031 N0 = reference_fe.N(0).subs( [x == xi0, y == xi1] )
00032
00033 p = matrix(2,1,[x, y])
00034 xi = Ginv*(p - x0)
00035 N0 = reference_fe.N(0).subs( [xi == xi[0,0], xi1 == xi[1,0]] )
00036
00037 print "fe.N(0) ", fe.N(0)
00038
00039
00040
00041 print "3D"
00042 initSyFi(3)
00043 x =cvar.x; y =cvar.y; z = cvar.z
00044
00045
00046 tetrahedron = Tetrahedron([1.3, 1.3, 1.3], [2.4, 2.1, 2.2 ], [1.4, 2.8, 1.3], [1.5, 1.3, 3.0])
00047
00048 G, x0 = geometry_mapping(tetrahedron)
00049 print "G ", G
00050 print "x0 ", x0
00051
00052 print G*matrix(3,1, [0,0,0]) + x0
00053 print G*matrix(3,1, [1,0,0]) + x0
00054 print G*matrix(3,1, [0,1,0]) + x0
00055 print G*matrix(3,1, [0,0,1]) + x0
00056
00057 Ginv = G.inverse()
00058 print Ginv*(matrix(3,1, [1.3,1.3,1.3]) - x0)
00059 print Ginv*(matrix(3,1, [2.4,2.1,2.2]) - x0)
00060 print Ginv*(matrix(3,1, [1.4,2.8,1.3]) - x0)
00061 print Ginv*(matrix(3,1, [1.5,1.3,3.0]) - x0)
00062
00063
00064 fe = Lagrange(tetrahedron)
00065
00066 print "global fe.N(0) ", fe.N(0)
00067
00068 reference_tetrahedron= ReferenceTetrahedron()
00069 reference_fe = Lagrange(reference_tetrahedron)
00070 xi0, xi1, xi2 = symbol("xi0"), symbol("xi1"), symbol("xi2")
00071 N0 = reference_fe.N(0).subs( [x == xi0, y == xi1, z == xi2] )
00072
00073 p = matrix(3,1,[x,y,z])
00074 xi = Ginv*(p - x0)
00075 N0 = reference_fe.N(0).subs( [xi == xi[0,0], xi1 == xi[1,0], xi2 == xi[2,0]] )
00076
00077 print "fe.N(0) ", fe.N(0)
00078
00079