SyFi  0.3
tests/sfc_jit/test.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 
00003 import unittest
00004 import os, sys, glob, shutil, commands
00005 
00006 import ufl
00007 
00008 from ufl import FiniteElement
00009 from ufl import VectorElement
00010 from ufl import TensorElement
00011 from ufl import MixedElement
00012 
00013 from ufl import Argument
00014 from ufl import TestFunction
00015 from ufl import TrialFunction
00016 
00017 from ufl import Coefficient
00018 from ufl import Constant
00019 
00020 from ufl import dx, ds
00021 
00022 import SyFi
00023 import newsfc as sfc
00024 
00025 from dolfin import Mesh, MeshEditor, assemble
00026 
00027 #import instant
00028 #instant.set_logging_level("error")
00029 #instant.set_logging_level("warning")
00030 #instant.set_logging_level("info")
00031 #instant.set_logging_level("debug")
00032 
00033 
00034 def num_integrals(form):
00035     return (form.num_cell_integrals(), form.num_exterior_facet_integrals(), form.num_interior_facet_integrals())
00036 
00037 cell2dim = { "interval": 1, "triangle": 2, "tetrahedron": 3, "quadrilateral": 2, "hexahedron": 3 }
00038 
00039 cell2volume = { "interval": 1.0, "triangle": 0.5, "tetrahedron": 1.0/6.0, "quadrilateral": 1.0, "hexahedron": 1.0 }
00040         
00041 
00042 def UnitCell(celltype):
00043     tdim = cell2dim[celltype]
00044     gdim = tdim
00045     mesh = Mesh()
00046     editor = MeshEditor()
00047     editor.open(mesh, celltype, tdim, gdim)
00048     if celltype == "interval":
00049         vertices = [(0.0,),
00050                     (1.0,)]
00051     if celltype == "triangle":
00052         vertices = [(0.0, 0.0),
00053                     (1.0, 0.0),
00054                     (0.0, 1.0)]
00055     if celltype == "tetrahedron":
00056         vertices = [(0.0, 0.0, 0.0),
00057                     (1.0, 0.0, 0.0),
00058                     (0.0, 1.0, 0.0),
00059                     (0.0, 0.0, 1.0)]
00060     if celltype == "quadrilateral":
00061         vertices = [(0.0, 0.0),
00062                     (1.0, 0.0),
00063                     (1.0, 1.0),
00064                     (0.0, 1.0)]
00065     if celltype == "hexahedron":
00066         vertices = [(0.0, 0.0, 0.0),
00067                     (1.0, 0.0, 0.0),
00068                     (1.0, 1.0, 0.0),
00069                     (0.0, 1.0, 0.0),
00070                     (0.0, 0.0, 1.0),
00071                     (1.0, 0.0, 1.0),
00072                     (1.0, 1.0, 1.0),
00073                     (0.0, 1.0, 1.0)]
00074     editor.initVertices(len(vertices))
00075     editor.initCells(1)
00076     for i, p in enumerate(vertices):
00077         editor.addVertex(i, *p)
00078     editor.addCell(0, *range(len(vertices)))
00079     editor.close()
00080     return mesh
00081 
00082 
00083 def assemble_on_cell(form, celltype, coeffs):
00084     "Assemble UFC form on a unit cell mesh and return the result as a float or numpy array."
00085     mesh = UnitCell(celltype)
00086     A = assemble(form, mesh, coeffs)
00087     if isinstance(A, float):
00088         return A
00089     return A.array()
00090 
00091 _test_temp_dir = "temp_dir"
00092 _done_test_temp_dir = "done_temp_dir"
00093 class SFCJitTest(unittest.TestCase):
00094     def setUp(self):
00095         #print "Running sfc jit test in testdir"
00096         #print "Imported SyFi from location", SyFi.__file__
00097         #print "Imported sfc  from location", sfc.__file__
00098         # Generate code in a clean directory:
00099         shutil.rmtree(_test_temp_dir, ignore_errors=True)
00100         os.mkdir(_test_temp_dir)
00101         os.chdir(_test_temp_dir)
00102     
00103     def tearDown(self):
00104         dirs = glob.glob("*")
00105         os.chdir("..")
00106         for d in dirs:
00107             os.rename(os.path.join(_test_temp_dir, d), os.path.join(_done_test_temp_dir, d))
00108     
00109     def testSetup(self):
00110         pass
00111     
00112     def _testJitVolume(self, polygon):
00113         "Test that the integral of 1.0 over a unit cell equals the length/area/volume of the unit cell."
00114         c = Constant(polygon)
00115         a = c*dx
00116         form = sfc.jit(a)
00117         self.assertTrue(form.rank() == 0)
00118         self.assertTrue(form.num_coefficients() == 1)
00119         self.assertTrue(num_integrals(form) == (1,0,0))
00120         A = assemble_on_cell(form, polygon, coeffs=[1.0])
00121         self.assertAlmostEqual(A, cell2volume[polygon])
00122     
00123     def testJitVolumeInterval(self):
00124         self._testJitVolume("interval")
00125 
00126     def testJitVolumeTriangle(self):
00127         self._testJitVolume("triangle")
00128     
00129     def testJitVolumeTetrahedron(self):
00130         self._testJitVolume("tetrahedron")
00131     
00132     def _testJitVolumeQuadrilateral(self): # Not supported by dolfin yet
00133         self._testJitVolume("quadrilateral")
00134     
00135     def _testJitVolumeHexahedron(self): # Not supported by dolfin yet
00136         self._testJitVolume("hexahedron")
00137 
00138     def _testJitConstant(self, polygon, degree):
00139         """Test that the integral of a constant coefficient over a unit 
00140         cell mesh equals the constant times the volume of the unit cell."""
00141         element = FiniteElement("CG", polygon, degree)
00142         f = Coefficient(element)
00143         a = f*dx
00144         form = sfc.jit(a)
00145         self.assertTrue(form.rank() == 0)
00146         self.assertTrue(form.num_coefficients() == 1)
00147         self.assertTrue(num_integrals(form) == (1,0,0))
00148         const = 1.23
00149         A = assemble_on_cell(form, polygon, coeffs=[const])
00150         self.assertAlmostEqual(A, const*cell2volume[polygon])
00151     
00152     def testJitConstantInterval(self):
00153         polygon = "interval"
00154         self._testJitConstant(polygon, 1)
00155         self._testJitConstant(polygon, 2)
00156 
00157     def testJitConstantTriangle(self):
00158         polygon = "triangle"
00159         self._testJitConstant(polygon, 1)
00160         self._testJitConstant(polygon, 2)
00161     
00162     def testJitConstantTetrahedron(self):
00163         polygon = "tetrahedron"
00164         self._testJitConstant(polygon, 1)
00165         self._testJitConstant(polygon, 2)
00166     
00167     def _testJitConstantQuadrilateral(self): # Not supported by dolfin yet
00168         polygon = "quadrilateral"
00169         self._testJitConstant(polygon, 1)
00170         self._testJitConstant(polygon, 2)
00171 
00172     def _testJitConstantHexahedron(self): # Not supported by dolfin yet
00173         polygon = "hexahedron"
00174         self._testJitConstant(polygon, 1)
00175         self._testJitConstant(polygon, 2)
00176     
00177     def testJitSource(self):
00178         "Test the source vector."
00179         element = FiniteElement("CG", "triangle", 1)
00180         v = TestFunction(element)
00181         f = Coefficient(element)
00182         a = f*v*dx
00183         form = sfc.jit(a)
00184         self.assertTrue(form.rank() == 1)
00185         self.assertTrue(form.num_coefficients() == 1)
00186         self.assertTrue(num_integrals(form) == (1,0,0))
00187         A = assemble_on_cell(form, "triangle", coeffs=[3.14])
00188         # TODO: Assert correct result
00189     
00190     def testJitMass(self):
00191         "Test the mass matrix."
00192         element = FiniteElement("CG", "triangle", 1)
00193         v = TestFunction(element)
00194         u = TrialFunction(element)
00195         f = Coefficient(element)
00196         a = f*u*v*dx
00197         form = sfc.jit(a)
00198         self.assertTrue(form.rank() == 2)
00199         self.assertTrue(form.num_coefficients() == 1)
00200         self.assertTrue(num_integrals(form) == (1,0,0))
00201         A = assemble_on_cell(form, "triangle", coeffs=[5.43])
00202         # TODO: Assert correct result
00203 
00204     def testJitSplitTerms(self):
00205         "Test a form split over two foo*dx terms, using the mass matrix."
00206         element = FiniteElement("CG", "triangle", 1)
00207         v = TestFunction(element)
00208         u = TrialFunction(element)
00209         f = Coefficient(element)
00210         a = u*v*dx + f*u*v*dx
00211         form = sfc.jit(a)
00212         self.assertTrue(form.rank() == 2)
00213         self.assertTrue(form.num_coefficients() == 1)
00214         self.assertTrue(num_integrals(form) == (1,0,0))
00215         A = assemble_on_cell(form, "triangle", coeffs=[4.43])
00216         # TODO: Assert correct result
00217 
00218 
00219 def test(verbosity=0):
00220     shutil.rmtree(_done_test_temp_dir, ignore_errors=True)
00221     os.mkdir(_done_test_temp_dir)
00222     
00223     classes = [SFCJitTest]
00224     suites = [unittest.makeSuite(c) for c in classes]
00225     testsuites = unittest.TestSuite(suites)
00226     unittest.TextTestRunner(verbosity=verbosity).run(testsuites)
00227 
00228 if __name__ == "__main__":
00229     test()
00230 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines