SyFi
0.3
|
00001 #!/usr/bin/env python 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 # ----------------------------------------------- Scalar forms: 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 # ----------------------------------------------- Vector forms: 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 # ----------------------------------------------- Vector boundary forms: 00056 00057 def load_vector(v, t, itg): 00058 """a(v; t) = \int t . v dx""" 00059 return inner(t, v) 00060 00061 00062 # ----------------------------------------------- Matrix forms: 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 # ----------------------------------------------- Boundary matrix forms: 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 # ----------------------------------------------- Testing: 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 # quicker, for debugging: 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 # creating scalar callback forms: 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 # creating vector callback forms: 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 # creating matrix callback forms: 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