← Back to branch summary

~chaffra/fiat/main

« back to all changes in this revision

Viewing changes to FIAT/Nedelec.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:
FIAT/Nedelec.py
 
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 modified 2 May 2005
 
8
 
 
9
import dualbasis, polynomial, functionalset, functional, shapes, points, \
 
10
       quadrature, Numeric, RaviartThomas
 
11
 
 
12
# (P_k)^d \circplus { p \in (P_k^H)^d : p(x)\cdot x = 0 }
 
13
# Only defined on tetrahedra
 
14
# Arnold decomposes as curl^t of RT plus grad P_k
 
15
# indexint starts at zero
 
16
 
 
17
def NedelecSpace( k ):
 
18
    shape = shapes.TETRAHEDRON
 
19
    d = shapes.dimension( shape )
 
20
    Vh = RaviartThomas.RTSpace( shape , k )
 
21
 
 
22
    Wh = polynomial.OrthogonalPolynomialSet( shape , k + 1 )
 
23
 
 
24
    curl_trans_rts = [ polynomial.curl_transpose( v ) \
 
25
                       for v in Vh ]
 
26
    
 
27
    curl_trans_rts_coeffs = Numeric.array( [ ctv.dof \
 
28
                                             for ctv in curl_trans_rts ] )
 
29
    curlTVh = polynomial.PolynomialSet( Vh.base , curl_trans_rts_coeffs )
 
30
 
 
31
    grad_whs = [ polynomial.gradient( w ) for w in Wh ]
 
32
    grad_whs_coeffs = Numeric.array( [ gw.dof for gw in grad_whs ] )
 
33
    gradWh = polynomial.PolynomialSet( Wh.base , grad_whs_coeffs )
 
34
 
 
35
    return polynomial.poly_set_union( curlTVh , gradWh )
 
36
 
 
37
class NedelecDual( dualbasis.DualBasis ):
 
38
    def __init__( self , U , k ):
 
39
        shape = shapes.TETRAHEDRON
 
40
        # tangent at k+1 points on each edge
 
41
        
 
42
        edge_pts = [ points.make_points( shape , \
 
43
                                         1 , i , k+2 ) \
 
44
                     for i in shapes.entity_range( shape , \
 
45
                                                   1 ) ]
 
46
 
 
47
        mdcb = functional.MakeDirectionalComponentBatch
 
48
        ls_per_edge = [ mdcb( edge_pts[i] , \
 
49
                              shapes.tangents[shape][i] ) \
 
50
                        for i in shape.entity_range( shape , 1 ) ]
 
51
 
 
52
        edge_ls = reduce( lambda a,b:a+b , ls_per_edge )
 
53
 
 
54
        # cross with normal at dim(P_{k-1}) points per face
 
55
        face_pts = [ points.make_points( shape , \
 
56
                                         2 , i , k+1 ) \
 
57
                     for i in shapes.entity_range( shape , \
 
58
                                                   2 ) ]
 
59
 
 
60
        # internal moments of dim( P_{k-2}) points
 
61
        internal_pts = points.make_points( shape , \
 
62
                                           3 , i , k+2 ) 
 
63
 
 
64
        entity_ids = {}
 
65
        entity_ids[0] = {}
 
66
        entity_ids[1] = {}
 
67
        cur = 0
 
68
        for i in shapes.entity_range(shape,1):
 
69
            entity_ids[1][i] = []
 
70
            for j in range(k+1):
 
71
                entity_ids[1][i].append(cur)
 
72
                cur += 1
 
73
        entity_ids[2] = {}
 
74
        for i in shapes.entity_range(shape,2):
 
75
            entity_ids[2][i] = []
 
76
            for j in range(len(face_pts[0])):
 
77
                entity_ids[2][i].append[cur]
 
78
                cur += 1
 
79
 
 
80
        entity_ids[3] = {}
 
81
        entity_ids[3][0] = []
 
82
        for i in len( face_pts ):
 
83
            entity_ids[3][0].append( cur )
 
84
            cur += 1
 
85
 
 
86
 
 
87
class Nedelec( polynomial.FiniteElement ):
 
88
    def __init__( self , k ):
 
89
        U = NedelecSpace( k )
 
90
        Udual = NedelecDual( k , U )
 
91
        polynomial.FiniteElement.__init__( self , Udual , U )
 
92