SyFi
0.3
|
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 #xi = matrix(2,1,[xi0, xi1]) 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 #xi = matrix(2,1,[xi0, xi1]) 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