00001 00002 from swiginac import * 00003 from SyFi import * 00004 from sfc.integral_tools import integral, geometry_mapping 00005 from sfc.symbolic_utils import symbol, symbols 00006 00007 setDigits(5) 00008 00009 print "" 00010 print "--------- 1D -----------------" 00011 print "" 00012 initSyFi(1) 00013 x = symbol("x") 00014 t = ReferenceLine() 00015 00016 f = x*x*x 00017 00018 print "symbolic integral of ",f," on reference line ", integral(t,f, True) 00019 print "numerical integral of ",f," on reference line ", integral(t,f, False,3) 00020 00021 print "" 00022 print "--------- 2D -----------------" 00023 print "" 00024 initSyFi(2) 00025 x, y = symbols("xy") 00026 t = ReferenceTriangle() 00027 00028 print "symbolic integral of ",f," on reference triangle ", integral(t,f, True) 00029 print "numerical integral of ",f," on reference triangle ", integral(t,f, False,3) 00030 00031 print "" 00032 print "--------- 3D -----------------" 00033 print "" 00034 initSyFi(3) 00035 x, y, z = symbols("xyz") 00036 f = 3.7*x - 2.8*y + 7.4*z + 3.14 00037 t = ReferenceTetrahedron() 00038 00039 print "symbolic integral of ",f," on reference tetrahedron ", integral(t,f, True) 00040 print "numerical integral of ",f," on reference tetrahedron ", integral(t,f, False,3) 00041 00042 00043 00044 00045 print "" 00046 print "--------- 2D -----------------" 00047 print "" 00048 f = 3.7*x - 2.8*y + 3.14 00049 initSyFi(2) 00050 x, y = symbols("xy") 00051 t = Triangle([1.3, 1.3], [2.4,2.1], [1.4, 2.8]) 00052 G, x0 = geometry_mapping(t) 00053 detG = G.determinant() 00054 fe = Lagrange(t) 00055 integrand = fe.N(0)*f 00056 00057 print "" 00058 print "integral of ",f," on global triangle ", integral(t, integrand, True, 3) 00059 print "" 00060 00061 xi0 = symbol("xi0"); xi1 = symbol("xi1") 00062 tr = ReferenceTriangle() 00063 fe = Lagrange(tr) 00064 N0 = fe.N(0).subs( [x == xi0, y == xi1] ) 00065 integrand = N0*f 00066 print "symbolic integral of ",f," on global triangle", integral(t, integrand, True, 1) 00067 print "numeric integral of ",f," on global triangle", integral(t, integrand, False, 2) 00068 00069 00070 print "" 00071 print "--------- 3D -----------------" 00072 print "" 00073 f = 3.7*x - 2.8*y + 7.4*z + 3.14 00074 initSyFi(3) 00075 x, y, z = symbols("xyz") 00076 t = Tetrahedron([1.3, 1.3, 1.3], [2.4, 1.1, 1.2 ], [1.4, 2.8, 1.3], [1.5, 1.3, 3.0]) 00077 G, x0 = geometry_mapping(t) 00078 detG = G.determinant() 00079 fe = Lagrange(t) 00080 integrand = fe.N(0)*f 00081 00082 print "integral of ",f," on global tetrahedron ", integral(t, integrand, True, 3) 00083 00084 xi0 = symbol("xi0"); xi1 = symbol("xi1"); xi2 = symbol("xi2") 00085 tr = ReferenceTetrahedron() 00086 fe = Lagrange(tr) 00087 N0 = fe.N(0).subs( [x == xi0, y == xi1, z == xi2] ) 00088 print "symbolic integral of ",f," on global tetrahedron", integral(t, integrand, True, 1) 00089 print "numeric integral of ",f," on global tetrahedron", integral(t, integrand, False, 3) 00090 00091