SyFi 0.3
element_representations.py
Go to the documentation of this file.
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 BasisFunction
00017 from ufl import TestFunction
00018 from ufl import TrialFunction
00019 
00020 from ufl import Function
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 testBasisFunctions(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 basisfunction(", i, c, ") for element ", repr(rep.ufl_element)
00123                     N = rep.basisfunction(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.basisfunction(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( self.P1rep.num_sub_elements == 0 )
00179         self.assertTrue( self.P2rep.num_sub_elements == 0 )
00180         self.assertTrue( self.Qrep.num_sub_elements  == 0 )
00181         self.assertTrue( self.VQrep.num_sub_elements == 2 )
00182         self.assertTrue( self.V1rep.num_sub_elements == 2 )
00183         self.assertTrue( self.V2rep.num_sub_elements == 2 )
00184         self.assertTrue( self.T1rep.num_sub_elements == 3 )
00185         self.assertTrue( self.T2rep.num_sub_elements == 3 )
00186         self.assertTrue( self.THrep.num_sub_elements == 2 )
00187         self.assertTrue( self.Mrep.num_sub_elements  == 3 )
00188 
00189 
00190 tests = [ElementRepresentationTest]
00191 
00192 if __name__ == "__main__":
00193     unittest.main()
00194 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines