00001
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
00028
00029
00030
00031
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
00096
00097
00098
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):
00133 self._testJitVolume("quadrilateral")
00134
00135 def _testJitVolumeHexahedron(self):
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):
00168 polygon = "quadrilateral"
00169 self._testJitConstant(polygon, 1)
00170 self._testJitConstant(polygon, 2)
00171
00172 def _testJitConstantHexahedron(self):
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
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
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
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