00001
00002
00003 from swiginac import *
00004 import sfc
00005 from sfc import *
00006 from sfc.symbolic_utils import *
00007 import sfc.common.options
00008 from sfc.common.utilities import get_callable_name
00009
00010
00011
00012 def constant_scalar(itg):
00013 """a(;) = \int 1 dx"""
00014 return numeric(1)
00015
00016 def L2_scalar(w, itg):
00017 """a(;w) = \int w^2 dx"""
00018 return inner(w, w)
00019
00020 def H1_semi_scalar(w, itg):
00021 """a(;w) = \int grad(w)^2 dx"""
00022 GinvT = itg.GinvT()
00023 Dw = grad(w, GinvT)
00024 return inner(Dw, Dw)
00025
00026 def H1_scalar(w, itg):
00027 """a(;w) = \int w^2 + grad(w)^2 dx"""
00028 GinvT = itg.GinvT()
00029 Dw = grad(w, GinvT)
00030 return inner(w, w) + inner(Dw, Dw)
00031
00032 scalar_forms = [constant_scalar, L2_scalar, H1_semi_scalar, H1_scalar]
00033
00034
00035
00036
00037
00038 def constant_vector(v, itg):
00039 """a(v;) = \int 1 dx"""
00040 return numeric(1)
00041
00042 def constant_source_vector(v, itg):
00043 """a(v;) = \int v dx"""
00044 return v
00045
00046 def source_vector(v, f, itg):
00047 """a(v; f) = \int f . v dx"""
00048 return inner(f, v)
00049
00050
00051 vector_forms = [constant_vector, constant_source_vector, source_vector]
00052
00053
00054
00055
00056
00057 def load_vector(v, t, itg):
00058 """a(v; t) = \int t . v dx"""
00059 return inner(t, v)
00060
00061
00062
00063
00064 def constant_matrix(v, u, itg):
00065 """a(v, u; ) = \int 1 dx"""
00066 return numeric(1)
00067
00068 def mass_matrix(v, u, itg):
00069 """a(v, u; ) = \int u . v dx"""
00070 return inner(v, u)
00071
00072 def mass_with_c_matrix(v, u, c, itg):
00073 """a(v, u; c) = \int c (u . v) dx"""
00074 return c * inner(v, u)
00075
00076 def stiffness_matrix(v, u, itg):
00077 GinvT = itg.GinvT()
00078 Du = grad(u, GinvT)
00079 Dv = grad(v, GinvT)
00080 return inner(Du, Dv)
00081
00082 def stiffness_with_M_matrix(v, u, M, itg):
00083 GinvT = itg.GinvT()
00084 Du = grad(u, GinvT)
00085 Dv = grad(v, GinvT)
00086 return inner(M * Du, Dv)
00087
00088 matrix_forms = [constant_matrix, mass_matrix, mass_with_c_matrix, stiffness_matrix, stiffness_with_M_matrix]
00089
00090
00091
00092
00093
00094 def mass_boundary_matrix(v, u, itg):
00095 """a(v, u; ) = \int u . v ds"""
00096 return inner(v, u)
00097
00098
00099
00100
00101
00102
00103 if __name__ == "__main__":
00104 import sys
00105
00106 sfc.common.options.add_debug_code = False
00107
00108 sfc.common.options.print_options()
00109
00110 args = set(sys.argv[1:])
00111
00112 print_forms = True if "p" in args else False
00113 generate = True if "g" in args else False
00114 compile = True if "c" in args else False
00115
00116 def check(form):
00117 form.sanity_check()
00118 if print_forms:
00119 print form
00120 if generate or compile:
00121 if compile:
00122 compiled_form = compile_form(form)
00123 print "Successfully compiled form:"
00124 print compiled_form
00125 print dir(compiled_form)
00126 else:
00127 res = write_ufc_code(form)
00128 print "Successfully generated form code:"
00129 print res
00130
00131 quad_order = 3
00132
00133 formcount = 0
00134 def form_name(callback):
00135 global formcount
00136 name = "form_%s_%d" % (get_callable_name(callback), formcount)
00137 formcount += 1
00138 return name
00139
00140 for nsd in [2, 3]:
00141 polygon = { 2: "triangle", 3: "tetrahedron" }[nsd]
00142 print "Using polygon = ", polygon
00143
00144 fe0 = FiniteElement("P0", polygon, 0)
00145 fe1 = FiniteElement("Lagrange", polygon, 1)
00146 fe2 = FiniteElement("Lagrange", polygon, 2)
00147 scalar_elements = [fe0, fe1, fe2]
00148
00149 vfe0 = VectorElement("P0", polygon, 0)
00150 vfe1 = VectorElement("Lagrange", polygon, 1)
00151 vfe2 = VectorElement("Lagrange", polygon, 2)
00152 vector_elements = [vfe0, vfe1, vfe2]
00153
00154 tfe0 = TensorElement("P0", polygon, 0)
00155 tfe1 = TensorElement("Lagrange", polygon, 1)
00156 tfe2 = TensorElement("Lagrange", polygon, 2)
00157 tensor_elements = [tfe0, tfe1, tfe2]
00158
00159
00160 scalar_elements = [fe1]
00161 vector_elements = [vfe1]
00162 tensor_elements = [tfe1]
00163
00164 all_elements = scalar_elements + vector_elements + tensor_elements
00165
00166
00167 for symbolic in [False, True]:
00168 print "Using symbolic = ", symbolic
00169
00170 options = { "symbolic": symbolic, "quad_order": quad_order }
00171
00172
00173
00174 print "Testing scalar forms"
00175
00176 for fe in all_elements:
00177
00178 callback = L2_scalar
00179 basisfunctions = []
00180 coefficients = [Function(fe)]
00181 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00182 check(form)
00183
00184 for fe in scalar_elements + vector_elements:
00185
00186 callback = H1_semi_scalar
00187 basisfunctions = []
00188 coefficients = [Function(fe)]
00189 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00190 check(form)
00191
00192 callback = H1_scalar
00193 basisfunctions = []
00194 coefficients = [Function(fe)]
00195 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00196 check(form)
00197
00198
00199
00200
00201 print "Testing vector forms"
00202
00203 for fe in all_elements:
00204 callback = constant_vector
00205 basisfunctions = [TestFunction(fe)]
00206 coefficients = []
00207 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00208 check(form)
00209
00210 for fe in scalar_elements:
00211 callback = constant_source_vector
00212 basisfunctions = [TestFunction(fe)]
00213 coefficients = []
00214 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00215 check(form)
00216
00217 for fe in scalar_elements:
00218 for f_fe in scalar_elements:
00219 callback = source_vector
00220 basisfunctions = [TestFunction(fe)]
00221 coefficients = [Function(f_fe)]
00222 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00223 check(form)
00224
00225 callback = load_vector
00226 basisfunctions = [TestFunction(fe)]
00227 coefficients = [Function(f_fe)]
00228 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, exterior_facet_integrands=[callback], options=options)
00229 check(form)
00230
00231 for fe in vector_elements:
00232 for f_fe in vector_elements:
00233 callback = source_vector
00234 basisfunctions = [TestFunction(fe)]
00235 coefficients = [Function(f_fe)]
00236 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00237 check(form)
00238
00239 callback = load_vector
00240 basisfunctions = [TestFunction(fe)]
00241 coefficients = [Function(f_fe)]
00242 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, exterior_facet_integrands=[callback], options=options)
00243 check(form)
00244
00245 for fe in tensor_elements:
00246 for f_fe in tensor_elements:
00247 callback = source_vector
00248 basisfunctions = [TestFunction(fe)]
00249 coefficients = [Function(f_fe)]
00250 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00251 check(form)
00252
00253 callback = load_vector
00254 basisfunctions = [TestFunction(fe)]
00255 coefficients = [Function(f_fe)]
00256 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, exterior_facet_integrands=[callback], options=options)
00257 check(form)
00258
00259
00260
00261
00262 print "Testing matrix forms"
00263
00264 for fe in all_elements:
00265
00266 callback = constant_matrix
00267 basisfunctions = [TestFunction(fe), TrialFunction(fe)]
00268 coefficients = []
00269 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00270 check(form)
00271
00272 callback = mass_matrix
00273 basisfunctions = [TestFunction(fe), TrialFunction(fe)]
00274 coefficients = []
00275 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00276 check(form)
00277
00278 callback = mass_boundary_matrix
00279 basisfunctions = [TestFunction(fe), TrialFunction(fe)]
00280 coefficients = []
00281 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, exterior_facet_integrands=[callback], options=options)
00282 check(form)
00283
00284 for c_fe in scalar_elements:
00285 callback = mass_with_c_matrix
00286 basisfunctions = [TestFunction(fe), TrialFunction(fe)]
00287 coefficients = [Function(c_fe)]
00288 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00289 check(form)
00290
00291 for fe in scalar_elements + vector_elements:
00292 callback = stiffness_matrix
00293 basisfunctions = [TestFunction(fe), TrialFunction(fe)]
00294 coefficients = []
00295 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00296 check(form)
00297
00298 for M_fe in scalar_elements + tensor_elements:
00299 callback = stiffness_with_M_matrix
00300 basisfunctions = [TestFunction(fe), TrialFunction(fe)]
00301 coefficients = [Function(M_fe)]
00302 form = Form(name=form_name(callback), basisfunctions=basisfunctions, coefficients=coefficients, cell_integrands=[callback], options=options)
00303 check(form)
00304