Subject: [FIAT-dev] Lagrange basis

From: Jed Brown <[email-masked]>
Date: Mon, 26 May 2008 00:00:11 +0200

I was playing with FIAT today and became somewhat confused with the choice of
Lagrange basis.  I'm just working in one dimension because I'm building high
order tensor products bases.

It appears that the Lagrange basis is based on evenly spaced points, but that's
terrible if the polynomial order grows.  Indeed, when I plotted the basis for
modest order, the Runge effect quickly becomes a problem.  So why is the
Lagrange basis not imeplemented on Legendre-Gauss-Lobatto nodes?  It's not
really harder to do since

  jacobi.compute_gauss_jacobi_points(1,1,p-1)

gives the interior points of an order p basis).  They are pretty much the same
at order 3, but by order 4 there is already an advantage to using the LGL
points.  Is there any advantage (ever) to using equally spaced nodes?  If you
want the points to nest, the Chebyshev points would be fine.  Am I missing
something?  I think it would be nice to be able to choose LGL/CGL/equal spacing,
but this requires some (maybe not much) refactoring.  Is it sufficient to have
make_lattice_xxx() take an extra parameter and add corresponding plumbing up
to the front-ends?

On another note, in K&S, the `standard modal' basis is chosen so that the
element mass and stiffness (for linear Laplacian) matrices are banded.  This
might be a nice basis to use for p-type refinement, but unfortunately it's
sparseness doesn't hold for nonlinear problems.  Does anyone (potential FIAT
users) care about hierarchical bases and sparser matrices for high order
Laplacians?

Jed

---

From: Jed Brown <[email-masked]>
Date: Mon, 26 May 2008 10:12:32 +0200

Here is a patch which adds LGL and CGL to the possible Lagrange LINE nodes.
Someone else can extend it to simplices if desired.  I didn't add plumbing up to
the user interface, but it's a one line change to choose a different set of
nodes.

Jed
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 10:11:41 2008 +0200
@@ -17,7 +17,7 @@
 from factorial import factorial
 from math import sqrt
 import numbering, reference # UFC vs. FIAT modules.
-
+import jacobi
 
 def strike_col( A , j ):
     m,n = A.shape
@@ -216,7 +216,7 @@
         td = td * ( n + i + 1 )
     return td / factorial( d )
 
-def make_lattice_line( n , interior = 0 ):
+def make_lattice_line_evenly_spaced( n , interior = 0 ):
     if n == 0:
         return tuple( )
     v0 = vertices[ LINE ][0][0]
@@ -224,6 +224,29 @@
     h = abs( v1 - v0 )
     return tuple( [ ( v0 + ( h * jj ) / n , ) \
                     for jj in xrange( interior , n + 1 - interior ) ] )
+
+def compute_lgl_points( n, interior = 0 ):
+    if n == 0:
+        return []
+    xi = jacobi.compute_gauss_jacobi_points(1,1,n-2)
+    return interior and xi or [-1] + xi + [1]
+
+def compute_cgl_points( n, interior = 0 ):
+    if n == 0:
+        return []
+    xi = n == 1 and [] or [ numpy.cos(i * numpy.pi / (n-1)) for i in xrange(1,n-1) ]
+    return interior and xi or [-1] + xi + [1]
+
+def make_lattice_line_gen( n , interior = 0 , pointgen = compute_lgl_points ):
+    if n == 0:
+        return tuple( )
+    v0 = vertices[ LINE ][0][0]
+    v1 = vertices[ LINE ][1][0]
+    h = ( v1 - v0 ) / 2
+    return tuple( [ (v0 + h * (x + 1),) for x in pointgen(n+1, interior) ] )
+
+def make_lattice_line( n , interior = 0 ):
+    return make_lattice_line_gen( n, interior, compute_lgl_points )
 
 def make_lattice_triangle( n , interior = 0 ):
     if n == 0: return tuple( )

---

From: Jed Brown <[email-masked]>
Date: Mon, 26 May 2008 19:55:15 +0200

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 )
 
-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:

---

From: "Robert Kirby" <[email-masked]>
Date: Mon, 26 May 2008 21:37:55 -0500

Dear Jed,Thanks for your interest in FIAT and your comments (I see the
patches that came in later emails).
I am aware about several of the issues you raise.  Let me respond:
1.) I agree that using "better" points have definite advantages in terms of
conditioning.  It's mainly a matter of FIAT being little used in 1d and
focusing on enabling hard elements on simplices -- when I wrote the initial
version, you had to go through quite a bit to obtain the simplicial (near)
Fekete points.  I just didn't do it.  If FIAT's 1d capability were used to
form tensor products, this would be easier to do.  Since I wrote FIAT, I
believe Tim Warburton has come up with a much better way to explicitly get
good points.  The make_lattice functions have the benefit of being very
simple to code in complete generality, both in 1d and higher-dimensional
simplices.
2.) The main use of FIAT is currently within FFC, which is better set up
algorithmically to handle low to moderate orders (very low by the standards
of the spectral world).  We aren't really hitting the point where equispaced
points hurt.
3.) Finding structure, both in the resulting mass/stiffness matrices as well
as in the generalized Vandermonde matrix, is a very interesting issue, but
one I haven't tackled.  I'm trying to be as general as possible at this
point.  I would like to see how to automatically infer these kinds of
structure within FIAT.  The disadvantage of the K&S approach is that you
work very hard to get one particular basis, but I want "all" the elements.
 Being able to compute with all of the elements has a certain value, even if
we don't get the full benefits of the specialized spectral element approach.

2008/5/25 Jed Brown <[email-masked]>:

> _______________________________________________
> FIAT-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/fiat-dev
>
>