SyFi 0.3
|
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 BasisFunction 00014 from ufl import TestFunction 00015 from ufl import TrialFunction 00016 00017 from ufl import Function 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 = Function(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 = Function(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 = Function(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 = Function(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