~chaffra/fiat/main

« back to all changes in this revision

Viewing changes to FIAT/PhiK.py

  • Committer: kirby
  • Date: 2005-05-31 15:03:57 UTC
  • Revision ID: devnull@localhost-20050531150357-n7ka5987yo6ylpaz
[project @ 2005-05-31 15:03:57 by kirby]
Initial revision

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
# Written by Robert C. Kirby
 
2
# Copyright 2005 by The University of Chicago
 
3
# Distributed under the LGPL license
 
4
# This work is partially supported by the US Department of Energy
 
5
# under award number DE-FG02-04ER25650
 
6
 
 
7
# last edited 9 May 2005
 
8
 
 
9
import shapes, functional, polynomial, points, functionalset
 
10
 
 
11
def PhiK( shape , n , U ):
 
12
    """The set PhiK defined by Brezzi & Fortin of vector-valued
 
13
    polynomials of degree n with vanishing divergence and
 
14
    vanishing normal components around the boundary."""
 
15
    U = polynomial.OrthogonalPolynomialArraySet( shape , n )
 
16
    # need Phi and U to share a common base, so I have to
 
17
    # take degree n and then throw away some to get degree n-1
 
18
    Phi = polynomial.OrthogonalPolynomialSet( shape , n )
 
19
    Phi_lower = Phi[:shapes.polynomial_dimension( shape , n-1 )]    
 
20
    DCPE = functional.DirectionalComponentPointEvaluation
 
21
    d = shapes.dimension( shape )
 
22
    pts_per_edge = [ [ x \
 
23
                       for x in points.make_points( shape , \
 
24
                                                    d-1 , \
 
25
                                                    i , \
 
26
                                                    n+d ) ] \
 
27
                    for i in shapes.entity_range( shape , d-1 ) ]
 
28
    pts_flat = reduce( lambda a,b:a+b , pts_per_edge )
 
29
 
 
30
    ls = [ functional.IntegralMomentOfDivergence( U , phi ) \
 
31
           for phi in Phi_lower ]
 
32
 
 
33
    for i in range(d+1):
 
34
        nrml = shapes.normals[ shape ][ i ]
 
35
        ls_cur = [ DCPE(U,nrml,pt) for pt in pts_per_edge[i] ]
 
36
        ls.extend(ls_cur)
 
37
    fset = functionalset.FunctionalSet(U,ls)
 
38
 
 
39
    return polynomial.ConstrainedPolynomialSet( fset )