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
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
00053
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
00077
00078 def stiffness_short(itg):
00079 """alternative short version"""
00080 print "Non-working prototype!"
00081
00082 for (i, j) in itg.indices:
00083 u, v = itg.v_basis_functions(i, j)
00084 Du, Dv = itg.v_basis_function_gradients(i, j)
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)
00096 Dv = itg.grad_v(1, j)
00097
00098 itg.A[(i,j)] = inner(Du, Dv)
00099
00100
00101
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
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
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
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
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
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