SyFi 0.3
callback_forms.py
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines