SyFi 0.3
scratch.py
Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 def boundary_callback(u, v, data):
00007     GinvT = data.GinvT()
00008     n     = data.n()
00009 
00010     Du = grad(u, GinvT)
00011     return inner(inner(n, Du), v)
00012 
00013 
00014 # example low level forms:
00015 
00016 def mass_v0(itg):
00017     """Compact version"""
00018     for i in range(itg.num_v_dofs(0)):
00019         for j in range(itg.num_v_dofs(1)):
00020             itg.A[(i,j)] = inner( itg.v_basis(0, i), itg.v_basis(1, j) )
00021 
00022 def mass_v1(itg):
00023     """More elaborate version"""
00024     for i in range(itg.num_v_dofs(0)):
00025         for j in range(itg.num_v_dofs(1)):
00026             u = itg.v_basis(0, i)
00027             v = itg.v_basis(1, j)
00028             itg.A[(i,j)] = inner(u, v)
00029 
00030 def mass_v2(itg):
00031     """'Jacobi' version"""
00032     dofs = itg.w_dofs(0)
00033     usum = itg.w_sum(0)
00034     for i in range(itg.num_v_dofs(0)):
00035         for j in range(itg.num_v_dofs(1)):
00036             u = diff(usum, dofs[i])
00037             v = itg.v_basis(1, j)
00038             itg.A[(i,j)] = inner(u, v)
00039 
00040 def stiffness_v1(itg):
00041     """symbolic-grad-version"""
00042     GinvT = itg.GinvT()
00043 
00044     for i in range(itg.num_v_dofs(0)):
00045         for j in range(itg.num_v_dofs(1)):
00046             u  = itg.v_basis(0, i)
00047             v  = itg.v_basis(1, j)
00048 
00049             Du = grad(u, GinvT)
00050             Dv = grad(v, GinvT)
00051 
00052             #Du = itg.add_token(Du) # TODO: do it like this?
00053             #Dv = itg.add_token(Dv)
00054 
00055             itg.A[(i,j)] = inner(Du, Dv)
00056 
00057 def stiffness_v2(itg):
00058     """'Manual Jacobi' version"""
00059     GinvT = itg.GinvT()
00060 
00061     dofs = itg.w_dofs(0)
00062     usum = itg.w_sum(0)
00063 
00064     Du = grad(usum, GinvT)
00065 
00066     for i in range(itg.num_v_dofs(0)):
00067         for j in range(itg.num_v_dofs(1)):
00068             v  = itg.v_basis(1, j)
00069             Dv = grad(v, GinvT)
00070 
00071             integrand = inner(Du, Dv)
00072             itg.A[(i,j)] = diff(integrand, dofs[i])
00073 
00074 
00075 
00076 # example prototype forms for alternate syntax:
00077 
00078 def stiffness_short(itg):
00079     """alternative short version"""
00080     print "Non-working prototype!"
00081     #f, g = itg.w_sums(0, 1) # TODO: do it like this?
00082     for (i, j) in itg.indices:
00083         u, v    = itg.v_basis_functions(i, j) # TODO: do it like this?
00084         Du, Dv  = itg.v_basis_function_gradients(i, j) # TODO: do it like this?
00085         itg.A[(i,j)] = inner(Du, Dv)
00086 
00087 def stiffness_v0(itg):
00088     """grad-from-itg-version"""
00089     print "Non-working prototype!"
00090     for i in range(itg.num_v_dofs(0)):
00091         for j in range(itg.num_v_dofs(1)):
00092             u  = itg.v_basis(0, i)
00093             v  = itg.v_basis(1, j)
00094 
00095             Du = itg.grad_v(0, i) # TODO: do it like this?
00096             Dv = itg.grad_v(1, j)
00097 
00098             itg.A[(i,j)] = inner(Du, Dv)
00099 
00100 
00101 # example nonlinear forms for automatic Jacobi creation:
00102 
00103 def nonlinear_boundary_F(itg):
00104     """F_i(w,f) = \int_\dOmega f(x) e^{t \cdot u(x)} (n \cdot grad u) v_i(x) ds
00105        Coefficients:
00106           w (u from last iteration)
00107           f
00108     """
00109     GinvT = itg.GinvT()
00110     n     = itg.n()
00111     t     = itg.t()
00112 
00113     w = itg.w_sum(0)
00114     f = itg.w_sum(1)
00115 
00116     Dw  = grad(w, GinvT)
00117     nDw = inner(n, Dw)
00118     wn  = inner(n, w)
00119     wt  = inner(t, w)
00120     exp_wt = exp( wt )
00121 
00122     assert isinstance(n, matrix)
00123     assert isinstance(w, matrix)
00124     assert n.nops() == w.nops()
00125     assert not isinstance(wn.evalm(), matrix)
00126     assert not isinstance(wt.evalm(), matrix)
00127 
00128     tmp = f * exp_wt * nDw
00129 
00130     for i in range(itg.num_v_dofs(0)):
00131         v = itg.v_basis(0, i)
00132         itg.A[(i,)] = inner(tmp, v)
00133 
00134 
00135 def nonlinear_F(itg):
00136     """Nonlinear test form F_i(w) = \int_\Omega w^2(x) w(x) v_i(x) dx]"""
00137     w = itg.w_sum(0)
00138     w2 = inner(w, w)
00139     for i in range(itg.num_v_dofs(0)):
00140         v = itg.v_basis(0, i)
00141         itg.A[(i,)] = w2 * inner(w, v)
00142 
00143 
00144 
00145 def stiffness_with_tokens_matrix(u, v, M, data):
00146     GinvT = data.GinvT()
00147 
00148     Du = grad(u, GinvT)
00149     Dv = grad(v, GinvT)
00150 
00151     # create manual temporary variables:
00152     Du = data.add_token("Du", Du)
00153     Dv = data.add_token("Dv", Dv)
00154 
00155     return inner(M * Du, Dv)
00156 
00157 
00158 
00159 if __name__ == "__main__":
00160 
00161     print_forms = False
00162     print_forms = True
00163 
00164     nsd = 2
00165     SyFi.initSyFi(nsd)
00166     polygon = SyFi.ReferenceTriangle()
00167     fe0 = SyFi.P0(polygon)
00168     fe1 = SyFi.Lagrange(polygon, 1)
00169     fe2 = SyFi.Lagrange(polygon, 2)
00170     vfe0 = SyFi.VectorP0(polygon)
00171     vfe1 = SyFi.VectorLagrange(polygon, 1)
00172     vfe2 = SyFi.VectorLagrange(polygon, 2)
00173     tfe0 = SyFi.TensorP0(polygon)
00174     tfe1 = SyFi.TensorLagrange(polygon, 1)
00175     tfe2 = SyFi.TensorLagrange(polygon, 2)
00176 
00177 
00178     fe_list = [fe1, fe1]
00179     #form = UserForm(rank=2, num_coefficients=0, name="mass", fe_list=fe_list, symbolic=False, quad_order=3)
00180     form = UserForm(rank=2, num_coefficients=0, name="mass", fe_list=fe_list, symbolic=True, quad_order=-1)
00181     if print_forms: print form
00182     form.sanity_check()
00183 
00184     citg = form.cell_integral()
00185     mass_v0(citg)
00186     if print_forms: print form
00187     form.sanity_check()
00188 
00189 
00190     # testing callback version of user interface:
00191 
00192     fe_list = [fe1, fe1]
00193     form = CallbackForm(name="mass", rank=2, num_coefficients=0, fe_list=fe_list, symbolic=True, quad_order=-1,
00194                       cell_integrands=[mass_callback])
00195     if print_forms: print form
00196     form.sanity_check()
00197 
00198 
00199     fe_list = [vfe1, vfe1, fe0]
00200     mf = CallbackForm(name="mass_with_c", rank=2, num_coefficients=1, fe_list=fe_list, symbolic=True, quad_order=-1,
00201                       cell_integrands=[mass_with_c_callback])
00202     if print_forms: print mf
00203 
00204     fe_list = [vfe1, vfe1, tfe0]
00205     form = CallbackForm(rank=2, fe_list=fe_list, symbolic=True, quad_order=-1,
00206                       cell_integrands=[stiffness_with_M_callback])
00207     if print_forms: print form
00208     form.sanity_check()
00209 
00210     fe_list = [fe2, fe2]
00211     form = CallbackForm(rank=2, fe_list=fe_list, symbolic=True, quad_order=-1,
00212                       exterior_facet_integrands=[boundary_callback])
00213     if print_forms: print form
00214     form.sanity_check()
00215 
00216 
00217     # testing coefficient counting
00218     fe_list = [fe1, fe0, fe2]
00219     form = UserForm(rank=0, num_coefficients=3, fe_list=fe_list)
00220     itg = form.cell_integral()
00221     if print_forms: print form
00222     form.sanity_check()
00223 
00224 
00225     # testing Jacobi version of user interface:
00226     fe_list = [vfe1, vfe1, fe0]
00227     F = UserForm(rank=1, num_coefficients=2, name="F", fe_list=fe_list, symbolic=True, quad_order=-1)
00228 
00229     itg = F.cell_integral()
00230     nonlinear_F(itg)
00231     itg.sanity_check()
00232 
00233     itg = F.exterior_facet_integral()
00234     nonlinear_boundary_F(itg)
00235     if print_forms: print F
00236     itg.sanity_check()
00237 
00238     J = Jacobi(F)
00239     if print_forms: print J
00240     J.sanity_check()
00241 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines