fiat team mailing list archive

[Jed Brown <jed@xxxxxxxx>]

Don't apply that proof-of-concept patch.  Here is a better one.  It keeps the
current choice of evenly spaced points, but puts in the plumbing.  With this
patch, you can choose point families:

  Lagrange.Lagrange(1,m,shapes.line_point_family_evenly_spaced)
  Lagrange.Lagrange(1,m,shapes.line_point_family_lgl)
  Lagrange.Lagrange(1,m,shapes.line_point_family_cgl)

If the user defines a function to evaluate interior points, they can put it in
the third argument.  If a third argument is not given, the current default
remains.  I think it would be a good idea to make LGL the default.  For
triangles and tetrahedra, we could keep the tensor product points (Legendre
would be better than the current choice) but other useful choices exist.
Especially the electrostatic points and the Fekete points.  Since these choices
would change the relationship between boundary and interior points, it's not
clear what the best abstraction is.  Therefore, the optional argument
point_family is currently ignored for triangles and tetrahedra.

Jed

diff -r bf2daae56306 FIAT/Lagrange.py
--- a/FIAT/Lagrange.py	Thu Apr 24 09:57:46 2008 -0500
+++ b/FIAT/Lagrange.py	Mon May 26 19:54:47 2008 +0200
@@ -17,7 +17,7 @@
 
 class LagrangeDual( dualbasis.DualBasis ):
     """Dual basis for the plain vanilla Lagrange finite element."""
-    def __init__( self , shape , n , U ):
+    def __init__( self , shape , n , U, point_family = 0 ):
         # Nodes are point evaluation on the lattice, which here is
         # ordered by topological dimension and entity within that
         # dimension.  For triangles we have vertex points plus
@@ -25,8 +25,9 @@
         # points on the interior of the triangle, or that many points
         # per face of tetrahedra plus max(0,(n-2)(n-1)n/6 points
         # in the interior
-        pts = [ [ shapes.make_points(shape,d,e,n) \
-                      for e in shapes.entity_range(shape,d) ] \
+        
+        pts = [ [ shapes.make_points(shape,d,e,n,point_family)
+                  for e in shapes.entity_range(shape,d) ] \
                     for d in shapes.dimension_range(shape) ]
 
         pts_flat = reduce( lambda a,b:a+b , \
@@ -57,14 +58,14 @@
         return
 
 class Lagrange( polynomial.FiniteElement ):
-    def __init__( self , shape , n ):
+    def __init__( self , shape , n, point_family = 0 ):
         if n < 1:
             raise RuntimeError, \
                   "Lagrange elements are only defined for n >= 1"
         self.shape = shape
         self.order = n
         U = polynomial.OrthogonalPolynomialSet( shape , n )
-        Udual = LagrangeDual( shape , n , U )
+        Udual = LagrangeDual( shape , n , U, point_family )
         polynomial.FiniteElement.__init__( self , \
                                            Udual , U )
 
diff -r bf2daae56306 FIAT/polynomial.py
--- a/FIAT/polynomial.py	Thu Apr 24 09:57:46 2008 -0500
+++ b/FIAT/polynomial.py	Mon May 26 19:54:47 2008 +0200
@@ -24,7 +24,7 @@
 class OrthonormalPolynomialBase( AbstractPolynomialBase ):
 	"""Defines an object that can tabulate the orthonormal polynomials
 	of a particular degree on a particular shape """
-	def __init__( self , shape , n ):
+	def __init__( self , shape , n , point_family = 0 ):
 		"""shape is a shape code defined in shapes.py	
 		n is the polynomial degree"""
 		self.bs = expansions.make_expansion( shape , n )
@@ -35,7 +35,7 @@
 		else:
 			dtildes = [ numpy.zeros( (len(self.bs),len(self.bs)) , "d" ) \
                       for i in range(0,d) ]
-			pts = shapes.make_points( shape , d , 0 , n + d + 1 )
+			pts = shapes.make_points( shape , d , 0 , n + d + 1 , point_family )
 			v = numpy.transpose( expansions.tabulators[shape](n,pts) )
 			vinv = numpy.linalg.inv( v )
 			dtildes = [numpy.transpose(a) \
diff -r bf2daae56306 FIAT/shapes.py
--- a/FIAT/shapes.py	Thu Apr 24 09:57:46 2008 -0500
+++ b/FIAT/shapes.py	Mon May 26 19:54:47 2008 +0200
@@ -17,6 +17,7 @@
 from factorial import factorial
 from math import sqrt
 import numbering, reference # UFC vs. FIAT modules.
+import jacobi
 
 
 def strike_col( A , j ):
@@ -216,16 +217,34 @@
         td = td * ( n + i + 1 )
     return td / factorial( d )
 
-def make_lattice_line( n , interior = 0 ):
+def line_point_family_evenly_spaced( n ):
+    return [ -1 + 2.0 * jj / n for jj in xrange (1, n) ]
+
+def line_point_family_lgl( n ):
+    return jacobi.compute_gauss_jacobi_points(1,1,n-1)
+
+def line_point_family_cgl( n ):
+    return [ numpy.cos( (n-jj) * numpy.pi / n ) for jj in xrange(1,n) ]
+
+def make_lattice_line( n , interior = 0 , point_family = 0 ):
+    if not point_family:
+        point_family = line_point_family_evenly_spaced
     if n == 0:
         return tuple( )
     v0 = vertices[ LINE ][0][0]
     v1 = vertices[ LINE ][1][0]
-    h = abs( v1 - v0 )
-    return tuple( [ ( v0 + ( h * jj ) / n , ) \
-                    for jj in xrange( interior , n + 1 - interior ) ] )
+    h = ( v1 - v0 ) / 2
+    if n < 2:
+        xi = []
+    else:
+        xi = point_family( n )
+    if interior:
+        xs = xi
+    else:
+        xs = [-1] + xi + [1]
+    return tuple( [ (v0 + h * (x + 1),) for x in xs ] )
 
-def make_lattice_triangle( n , interior = 0 ):
+def make_lattice_triangle( n , interior = 0, point_family = 0 ):
     if n == 0: return tuple( )
     vs = [ numpy.array( vertices[ TRIANGLE ][ i ] ) \
              for i in vertices[ TRIANGLE ].iterkeys() ]
@@ -234,7 +253,7 @@
                     for ii in xrange( interior , n + 1 - interior ) \
                         for jj in xrange( interior , n - ii + 1 - interior ) ] )
 
-def make_lattice_tetrahedron( n , interior = 0 ):
+def make_lattice_tetrahedron( n , interior = 0, point_family = 0 ):
     if n == 0: return tuple( )
     vs = [ numpy.array( vertices[ TETRAHEDRON ][ i ] ) \
              for i in vertices[ TETRAHEDRON ].iterkeys() ]
@@ -253,8 +272,8 @@
                   TRIANGLE : make_lattice_triangle , \
                   TETRAHEDRON : make_lattice_tetrahedron }
 
-def make_lattice( shp , n , interior = 0 ):
-    return lattice_funcs[ shp ]( n , interior )    
+def make_lattice( shp , n , interior = 0, point_family = 0 ):
+    return lattice_funcs[ shp ]( n , interior, point_family ) 
 
 def make_pt_to_edge( shp , verts ):
     """verts is the tuple of vertex ids on the reference shape.
@@ -311,12 +330,12 @@
                             2 : pt_to_face[ TETRAHEDRON ] , \
                             3 : lambda y: lambda x: x } }
                             
-def make_points( shp , dim , entity_id , order ):
+def make_points( shp , dim , entity_id , order, point_family = 0 ):
     if dim == 0:
         return ( vertices[ shp ][ entity_id ] , )
     if dim == shp:
         if entity_id == 0:
-            return make_lattice( shp , order , 1 )
+            return make_lattice( shp , order , 1, point_family )
         else:
             raise ShapeError, "Can't make those points"
     else:

Attachment: pgpV7AmPxzsal.pgp
Description: PGP signature