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