SyFi
0.3
|
00001 #!/usr/bin/env python 00002 00003 __authors__ = "Martin Sandve Alnes" 00004 __date__ = "2008-09-09 2008-09-10" 00005 00006 import unittest 00007 import os, sys, glob, shutil, commands 00008 import swiginac 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 00022 00023 from ufl import dx, ds 00024 00025 import SyFi 00026 import sfc as sfc 00027 from sfc.representation import ElementRepresentation 00028 from sfc.quadrature import find_quadrature_rule 00029 00030 import instant 00031 instant.set_logging_level("warning") 00032 00033 00034 00035 class ElementRepresentationTest(unittest.TestCase): 00036 def __init__(self, *args, **kwargs): 00037 unittest.TestCase.__init__(self, *args, **kwargs) 00038 00039 def setUp(self): 00040 SyFi.initSyFi(2) 00041 polygon = "triangle" 00042 00043 self.P1 = FiniteElement("CG", polygon, 1) 00044 self.P2 = FiniteElement("CG", polygon, 2) 00045 self.Q = FiniteElement("Q", polygon, 1) 00046 self.VQ = VectorElement("Q", polygon, 1) 00047 self.V1 = VectorElement("CG", polygon, 1) 00048 self.V2 = VectorElement("CG", polygon, 2) 00049 self.T1 = TensorElement("CG", polygon, 1, symmetry=True) 00050 self.T2 = TensorElement("CG", polygon, 2, symmetry=True) 00051 self.TH = self.V2 * self.P1 00052 self.M = MixedElement(self.T1, self.V1, self.P1) 00053 00054 self.P1rep = ElementRepresentation(self.P1) 00055 self.P2rep = ElementRepresentation(self.P2) 00056 self.Qrep = ElementRepresentation(self.Q, quad_rule = find_quadrature_rule(polygon, 2)) 00057 self.VQrep = ElementRepresentation(self.VQ, quad_rule = find_quadrature_rule(polygon, 2)) 00058 self.V1rep = ElementRepresentation(self.V1) 00059 self.V2rep = ElementRepresentation(self.V2) 00060 self.T1rep = ElementRepresentation(self.T1) 00061 self.T2rep = ElementRepresentation(self.T2) 00062 self.THrep = ElementRepresentation(self.TH) 00063 self.Mrep = ElementRepresentation(self.M) 00064 00065 self.reps = [getattr(self, d) for d in dir(self) if isinstance(getattr(self, d), ElementRepresentation)] 00066 00067 def tearDown(self): 00068 pass 00069 00070 def testSetup(self): 00071 pass 00072 00073 def testValueShapes(self): 00074 self.assertTrue( self.P1rep.value_shape == () ) 00075 self.assertTrue( self.P2rep.value_shape == () ) 00076 self.assertTrue( self.Qrep.value_shape == () ) 00077 self.assertTrue( self.VQrep.value_shape == (2,) ) 00078 self.assertTrue( self.V1rep.value_shape == (2,) ) 00079 self.assertTrue( self.V2rep.value_shape == (2,) ) 00080 self.assertTrue( self.T1rep.value_shape == (2,2) ) 00081 self.assertTrue( self.T2rep.value_shape == (2,2) ) 00082 self.assertTrue( self.THrep.value_shape == (2+1,) ) 00083 self.assertTrue( self.Mrep.value_shape == (4+2+1,) ) 00084 00085 def testGeometryDimensions(self): 00086 self.assertTrue( self.P1rep.geometric_dimension == 2 ) 00087 self.assertTrue( self.P2rep.geometric_dimension == 2 ) 00088 self.assertTrue( self.Qrep.geometric_dimension == 2 ) 00089 self.assertTrue( self.VQrep.geometric_dimension == 2 ) 00090 self.assertTrue( self.V1rep.geometric_dimension == 2 ) 00091 self.assertTrue( self.V2rep.geometric_dimension == 2 ) 00092 self.assertTrue( self.T1rep.geometric_dimension == 2 ) 00093 self.assertTrue( self.T2rep.geometric_dimension == 2 ) 00094 self.assertTrue( self.THrep.geometric_dimension == 2 ) 00095 self.assertTrue( self.Mrep.geometric_dimension == 2 ) 00096 00097 def testTopologyDimensions(self): 00098 self.assertTrue( self.P1rep.topological_dimension == 2 ) 00099 self.assertTrue( self.P2rep.topological_dimension == 2 ) 00100 self.assertTrue( self.Qrep.topological_dimension == 2 ) 00101 self.assertTrue( self.VQrep.topological_dimension == 2 ) 00102 self.assertTrue( self.V1rep.topological_dimension == 2 ) 00103 self.assertTrue( self.V2rep.topological_dimension == 2 ) 00104 self.assertTrue( self.T1rep.topological_dimension == 2 ) 00105 self.assertTrue( self.T2rep.topological_dimension == 2 ) 00106 self.assertTrue( self.THrep.topological_dimension == 2 ) 00107 self.assertTrue( self.Mrep.topological_dimension == 2 ) 00108 00109 def testArguments(self): 00110 # Tests that for each basisfunction there is at least one nonzero value component 00111 #for rep in (self.P1rep, self.V1rep, self.T1rep): 00112 #for rep in (self.V1rep, ): 00113 for rep in self.reps: 00114 if rep.ufl_element.family() == "Quadrature"\ 00115 or rep.ufl_element.family() == "Boundary Quadrature": 00116 continue 00117 #print 00118 #print rep 00119 for i in range(rep.local_dimension): 00120 zeros = 0 00121 for c in ufl.permutation.compute_indices(rep.value_shape): 00122 #print "Calling basis_function(", i, c, ") for element ", repr(rep.ufl_element) 00123 N = rep.basis_function(i, c) 00124 #print i, c, N 00125 if N == 0.0: 00126 zeros += 1 00127 self.assertTrue(zeros < rep.value_size) 00128 if rep.local_dimension > 1 and rep.value_size > 1: 00129 self.assertTrue(zeros > 0) # This should hold for compositions of scalar elements 00130 # FIXME: Test something more 00131 00132 def testDofTopology(self): 00133 def triangle_entities(): 00134 for d in range(2): 00135 for i in range((3,3,1)[d]): 00136 yield (d,i) 00137 for rep in self.reps: 00138 for (d,i) in triangle_entities(): 00139 dofs = rep.entity_dofs[d][i] 00140 self.assertTrue(len(dofs) == rep.num_entity_dofs[d]) 00141 # FIXME: Test something more 00142 00143 def testPointEvaluation(self): 00144 # Tests that the property { N_k(\xi_j) == \delta_{kj} } holds for scalar elements. 00145 for rep in (self.P1rep, self.P2rep): 00146 for i in range(rep.local_dimension): 00147 x = rep.dof_x[i] 00148 xi = rep.dof_xi[i] 00149 00150 repmap = swiginac.exmap() 00151 for j in range(2): 00152 repmap[rep.p[j]] = xi[j] 00153 00154 for k in range(rep.local_dimension): 00155 c = () # Assuming scalar element here: 00156 N = rep.basis_function(k, c) 00157 Nxi = N.subs(repmap) 00158 if i == k: 00159 self.assertTrue(Nxi == 1.0) 00160 else: 00161 self.assertTrue(Nxi == 0.0) 00162 00163 def testCoordinates(self): 00164 for rep in self.reps: 00165 #print 00166 #print rep.ufl_element 00167 for i in range(rep.local_dimension): 00168 x = rep.dof_x[i] 00169 xi = rep.dof_xi[i] 00170 #print i, xi 00171 #x2 = rep.xi_to_x(xi) 00172 #xi2 = rep.x_to_xi(x) 00173 #self.assertTrue(x2 == x) 00174 #self.assertTrue(xi2 == xi) 00175 # FIXME: Test something more 00176 00177 def testSubHierarchy(self): 00178 self.assertTrue(len(self.P1rep.sub_elements) == 0 ) 00179 self.assertTrue(len(self.P2rep.sub_elements) == 0 ) 00180 self.assertTrue(len(self.Qrep.sub_elements) == 0 ) 00181 self.assertTrue(len(self.VQrep.sub_elements) == 2 ) 00182 self.assertTrue(len(self.V1rep.sub_elements) == 2 ) 00183 self.assertTrue(len(self.V2rep.sub_elements) == 2 ) 00184 self.assertTrue(len(self.T1rep.sub_elements) == 3 ) 00185 self.assertTrue(len(self.T2rep.sub_elements) == 3 ) 00186 self.assertTrue(len(self.THrep.sub_elements) == 2 ) 00187 self.assertTrue(len(self.Mrep.sub_elements) == 3 ) 00188 00189 00190 tests = [ElementRepresentationTest] 00191 00192 if __name__ == "__main__": 00193 unittest.main() 00194