SyFi
0.3
|
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