SyFi
0.3
|
00001 #!/usr/bin/env python 00002 00003 __authors__ = "Martin Sandve Alnes" 00004 __date__ = "2008-09-04 -- 2012-05-16" 00005 00006 import unittest 00007 import os, sys, glob, shutil, commands 00008 00009 import ufl 00010 00011 from ufl import FiniteElement 00012 from ufl import VectorElement 00013 from ufl import TensorElement 00014 from ufl import MixedElement 00015 00016 from ufl import Argument 00017 from ufl import TestFunction 00018 from ufl import TrialFunction 00019 00020 from ufl import Coefficient 00021 from ufl import Constant, VectorConstant, TensorConstant 00022 00023 from ufl import dx, ds 00024 from ufl.common import product 00025 00026 import SyFi 00027 import sfc 00028 from sfc.representation.geometryrepresentation import GeometryRepresentation 00029 from sfc.representation.elementrepresentation import ElementRepresentation 00030 import numpy 00031 00032 from test_tempdir_base import TempDirTestBase 00033 from cell_assembly import assemble_on_cell, cell2volume, num_integrals 00034 00035 simplices = (ufl.interval, ufl.triangle, ufl.tetrahedron) 00036 quads = (ufl.quadrilateral, ufl.hexahedron) 00037 00038 #from sfc.representation import FormArgumentRepresentation 00039 from sfc.codegeneration.formargumentcg import FormArgumentCG 00040 00041 try: 00042 from uflacs.codeutils.format_code_structure import format_code_structure 00043 except: 00044 def format_code_structure(x): return str(x) 00045 00046 class ElementRepTest(unittest.TestCase): 00047 def test_foo(self): 00048 SyFi.initSyFi(2) 00049 00050 family = "CG" 00051 ufl_cell = ufl.triangle 00052 degree = 1 00053 ufl_element = FiniteElement(family, ufl_cell, degree) 00054 geomrep = GeometryRepresentation(ufl_cell.domain()) 00055 er = ElementRepresentation(ufl_element, geomrep, quad_rule=None, options=None, cache=None) 00056 print 00057 print 00058 print dir(er) 00059 print 00060 00061 ngrad = 0 00062 code = [] 00063 for dof in range(er.local_dimension): 00064 #sub_element_index = er.dof_to_sub_element_index(dof) 00065 00066 components = [()] # FIXME: permutations of value shape 00067 00068 if ngrad: 00069 assert ngrad == 1 00070 all_derivatives = [(d,) for d in range(nsd)] # FIXME: permutations of ngrad 00071 00072 for component in components: 00073 for directions in all_derivatives: 00074 expr = er.basis_function_derivative(dof, component, directions) 00075 name = "FIXME" 00076 code.append((name, expr.printc())) 00077 else: 00078 for component in components: 00079 expr = er.basis_function(dof, component) 00080 name = "FIXME" 00081 code.append((name, expr.printc())) 00082 00083 print "Element code:" 00084 for n,e in code: 00085 print "%s = %s;" % (n, e) 00086 00087 class FormArgumentCodeGenTest(object):#unittest.TestCase): 00088 def xtest_form_argument_representation_interval(self): 00089 ufl_cell = ufl.interval 00090 00091 SyFi.initSyFi(ufl_cell.d) 00092 00093 # FIXME: Must FormArgumentCG really depend on all this stuff? 00094 #formdata, element_reps, geomrep, options = ... 00095 formrep = None #FormRepresentation(formdata, element_reps, geomrep, options) 00096 farep = FormArgumentRepresentation(formrep) 00097 facg = FormArgumentCG(farep, False) 00098 decl = facg.generate_declarations() 00099 00100 self.assertTrue(decl) 00101 if 1: 00102 print 00103 print ufl_cell 00104 print format_code_structure(decl) 00105 print 00106 00107 def xtest_form_argument_representation_triangle(self): 00108 ufl_cell = ufl.triangle 00109 00110 SyFi.initSyFi(ufl_cell.d) 00111 00112 farep = FormArgumentRepresentation() 00113 facg = FormArgumentCG(farep) # FIXME 00114 decl = facg.generate_declarations() 00115 00116 self.assertTrue(decl) 00117 if 1: 00118 print 00119 print ufl_cell 00120 print format_code_structure(decl) 00121 print 00122 00123 def xtest_form_argument_representation_tetrahedron(self): 00124 ufl_cell = ufl.tetrahedron 00125 00126 SyFi.initSyFi(ufl_cell.d) 00127 00128 farep = FormArgumentRepresentation() 00129 facg = FormArgumentCG(farep) # FIXME 00130 decl = facg.generate_declarations() 00131 00132 self.assertTrue(decl) 00133 if 1: 00134 print 00135 print ufl_cell 00136 print format_code_structure(decl) 00137 print 00138 00139 def xtest_form_argument_representation_quadrilateral(self): 00140 ufl_cell = ufl.quadrilateral 00141 def xtest_form_argument_representation_hexahedron(self): 00142 ufl_cell = ufl.hexahedron 00143 00144 00145 class FormArgumentTestMixin(object): 00146 00147 def assertAlmostEqualNumpy(self, actual, expected, eps=1e-8): 00148 self.assertEqual(actual.shape, expected.shape) 00149 mdiff = expected - actual 00150 vdiff = mdiff.reshape((product(mdiff.shape),)) 00151 sdiff = numpy.dot(vdiff,vdiff) 00152 self.assertAlmostEqual(sdiff, 0.0, eps) 00153 00154 def assertAssembledAlmostEqual(self, ufl_form, ufl_cell, coeffs, expected): 00155 """Test that the integral of 1.0 over a unit cell 00156 equals the length/area/volume of the unit cell.""" 00157 ufc_form, module, formdata, prefix = sfc.jit(ufl_form, parameters=self.options) 00158 assembled = assemble_on_cell(ufc_form, ufl_cell, coeffs=coeffs) 00159 if isinstance(assembled, float): 00160 self.assertAlmostEqual(assembled, expected) 00161 else: 00162 self.assertAlmostEqualNumpy(assembled, expected) 00163 00164 00165 def _test_constant(self, cell): 00166 "Test that a single constant is handled correctly." 00167 D = cell2volume[cell] 00168 c = Constant(cell) 00169 a = c*dx 00170 coeffs = [1.2] 00171 expected = 1.2 * D 00172 self.assertAssembledAlmostEqual(a, cell, coeffs, expected) 00173 00174 def xtest_constant_interval(self): 00175 self._test_constant(ufl.interval) 00176 def xtest_constant_triangle(self): 00177 self._test_constant(ufl.triangle) 00178 def xtest_constant_tetrahedron(self): 00179 self._test_constant(ufl.tetrahedron) 00180 00181 def _test_constant_quadrilateral(self): # TODO: Not supported by dolfin yet 00182 self._test_constant(ufl.quadrilateral) 00183 def _test_constant_hexahedron(self): # TODO: Not supported by dolfin yet 00184 self._test_constant(ufl.hexahedron) 00185 00186 00187 def _test_constants(self, cell): 00188 "Test that multiple constants are handled correctly." 00189 D = cell2volume[cell] 00190 c1 = Constant(cell) 00191 c2 = Constant(cell) 00192 c3 = Constant(cell) 00193 a = c1*c2*c3*dx 00194 coeffs = [2,3,5] 00195 expected = coeffs[0]*coeffs[1]*coeffs[2] * D 00196 self.assertAssembledAlmostEqual(a, cell, coeffs, expected) 00197 00198 def xtest_constants_interval(self): 00199 self._test_constants(ufl.interval) 00200 def xtest_constants_triangle(self): 00201 self._test_constants(ufl.triangle) 00202 def xtest_constants_tetrahedron(self): 00203 self._test_constants(ufl.tetrahedron) 00204 00205 def _test_constants_quadrilateral(self): # TODO: Not supported by dolfin yet 00206 self._test_constants(ufl.quadrilateral) 00207 def _test_constants_hexahedron(self): # TODO: Not supported by dolfin yet 00208 self._test_constants(ufl.hexahedron) 00209 00210 00211 def _test_vector_constant(self, cell): 00212 "Test that a vector valued constant is handled correctly." 00213 D = cell2volume[cell] 00214 c = VectorConstant(cell) 00215 a = sum(c[i] for i in range(cell.d))*dx 00216 coeffs = [(1.2,3.4,5.6)[:cell.d]] 00217 expected = sum(coeffs[0]) * D 00218 self.assertAssembledAlmostEqual(a, cell, coeffs, expected) 00219 00220 def xtest_vector_constant_interval(self): 00221 self._test_vector_constant(ufl.interval) 00222 def xtest_vector_constant_triangle(self): 00223 self._test_vector_constant(ufl.triangle) 00224 def xtest_vector_constant_tetrahedron(self): 00225 self._test_vector_constant(ufl.tetrahedron) 00226 00227 def _test_vector_constant_quadrilateral(self): # TODO: Not supported by dolfin yet 00228 self._test_vector_constant(ufl.quadrilateral) 00229 def _test_vector_constant_hexahedron(self): # TODO: Not supported by dolfin yet 00230 self._test_vector_constant(ufl.hexahedron) 00231 00232 00233 def _test_tensor_constant(self, cell): 00234 "Test that a tensor valued constant is handled correctly." 00235 D = cell2volume[cell] 00236 c = TensorConstant(cell) 00237 sh = c.shape() 00238 a = sum(c[i,j] for i in range(sh[0]) for j in range(sh[1]))*dx 00239 rows = ( (1.1,1.2,1.3)[:sh[0]], 00240 (2.1,2.2,2.3)[:sh[0]], 00241 (3.1,3.2,3.3)[:sh[0]] )[:sh[1]] 00242 coeffs = [rows] 00243 expected = D * sum(rows[i][j] for i in range(sh[0]) for j in range(sh[1])) 00244 self.assertAssembledAlmostEqual(a, cell, coeffs, expected) 00245 00246 def xtest_tensor_constant_interval(self): 00247 self._test_tensor_constant(ufl.interval) 00248 def xtest_tensor_constant_triangle(self): 00249 self._test_tensor_constant(ufl.triangle) 00250 def xtest_tensor_constant_tetrahedron(self): 00251 self._test_tensor_constant(ufl.tetrahedron) 00252 00253 def _test_tensor_constant_quadrilateral(self): # TODO: Not supported by dolfin yet 00254 self._test_tensor_constant(ufl.quadrilateral) 00255 def _test_tensor_constant_hexahedron(self): # TODO: Not supported by dolfin yet 00256 self._test_tensor_constant(ufl.hexahedron) 00257 00258 00259 def _test_coefficients(self, cell): 00260 self._test_coefficient(cell, "DG", 0) 00261 self._test_coefficient(cell, "DG", 1) 00262 #self._test_coefficient(cell, "DG", 2) 00263 #self._test_coefficient(cell, "DG", 3) 00264 00265 #self._test_coefficient(cell, "CG", 1) 00266 #self._test_coefficient(cell, "CG", 2) 00267 #self._test_coefficient(cell, "CG", 3) 00268 00269 def _test_coefficient(self, cell, family="CG", degree=1): 00270 "Test that a single coefficient is handled correctly." 00271 D = cell2volume[cell] 00272 V = FiniteElement(family, cell, degree) 00273 c = Coefficient(V) 00274 a = c*dx 00275 coeffs = [1.2] 00276 expected = 1.2 * D 00277 self.assertAssembledAlmostEqual(a, cell, coeffs, expected) 00278 00279 # TODO: Make these work with uflacs: 00280 def xtest_coefficients_interval(self): 00281 self._test_coefficients(ufl.interval) 00282 def xtest_coefficients_triangle(self): 00283 self._test_coefficients(ufl.triangle) 00284 def xtest_coefficients_tetrahedron(self): 00285 self._test_coefficients(ufl.tetrahedron) 00286 00287 def xtest_coefficients_quadrilateral(self): # TODO: Not supported by dolfin yet 00288 self._test_coefficients(ufl.quadrilateral) 00289 def xtest_coefficients_hexahedron(self): # TODO: Not supported by dolfin yet 00290 self._test_coefficients(ufl.hexahedron) 00291 00292 00293 def test_code_snippets(self): 00294 print 00295 cell = ufl.triangle 00296 00297 # Mock up a formrep object 00298 class Mock(object): pass 00299 formrep = Mock() 00300 formrep.geomrep = Mock() 00301 formrep.geomrep.sfc_cell = Mock() 00302 formrep.geomrep.sfc_cell.nsd = cell.d 00303 00304 facg = FormArgumentCG(formrep, on_facet=False) 00305 00306 S = FiniteElement("CG", cell, 1) 00307 u = Argument(S, count=1) 00308 f = Coefficient(S, count=2) 00309 print 00310 print facg.gen_argument_name(u, (), (), None) 00311 print facg.gen_argument_name(u, (0,), (), None) 00312 print facg.gen_coefficient_name(f, (), (), None) 00313 print facg.gen_coefficient_name(f, (0,), (), None) 00314 print 00315 00316 V = VectorElement("CG", cell, 1) 00317 v = Argument(V, count=1) 00318 g = Coefficient(V, count=2) 00319 print 00320 print facg.gen_argument_name(v, (), (0,), None) 00321 print facg.gen_argument_name(v, (0,), (1,), None) 00322 print facg.gen_coefficient_name(g, (), (0,), None) 00323 print facg.gen_coefficient_name(g, (0,), (1,), None) 00324 print 00325 00326 T = TensorElement("CG", cell, 1) 00327 w = Argument(T, count=1) 00328 h = Coefficient(T, count=2) 00329 print 00330 print facg.gen_argument_name(w, (), (0,1), None) 00331 print facg.gen_argument_name(w, (0,), (1,0), None) 00332 print facg.gen_coefficient_name(h, (), (0,1), None) 00333 print facg.gen_coefficient_name(h, (0,), (1,0), None) 00334 print 00335 00336 print 00337 00338 if 1: 00339 class UflacsFormArgumentTest(TempDirTestBase, FormArgumentTestMixin): 00340 def __init__(self, *args, **kwargs): 00341 TempDirTestBase.__init__(self, *args, **kwargs) 00342 00343 def setUp(self): 00344 TempDirTestBase.setUp(self) 00345 self.options.code.integral.integration_method = "quadrature" 00346 self.set_uflacs_options() 00347 00348 def set_uflacs_options(self): 00349 try: 00350 import uflacs 00351 self.options.code.integral.use_uflacs = True 00352 except: 00353 print "UFLACS NOT AVAILABLE" 00354 00355 if 0: 00356 class SfcQuadratureFormArgumentTest(TempDirTestBase, FormArgumentTestMixin): 00357 def __init__(self, *args, **kwargs): 00358 TempDirTestBase.__init__(self, *args, **kwargs) 00359 00360 def setUp(self): 00361 TempDirTestBase.setUp(self) 00362 self.options.code.integral.integration_method = "quadrature" 00363 00364 if 0: 00365 class SfcSymbolicFormArgumentTest(TempDirTestBase, FormArgumentTestMixin): 00366 def __init__(self, *args, **kwargs): 00367 TempDirTestBase.__init__(self, *args, **kwargs) 00368 00369 def setUp(self): 00370 TempDirTestBase.setUp(self) 00371 self.options.code.integral.integration_method = "symbolic"