SyFi  0.3
test_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 -- 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 )
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines