~chaffra/fiat/main : revision 2

1

# Written by Robert C. Kirby

2

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

11

12

13

def RT0Space( shape ):

14

d = shapes.dimension( shape )

15

vec_P1 = polynomial.OrthogonalPolynomialArraySet( shape , 1 )

16

dimP1 = shapes.polynomial_dimension( shape , 1 )

17

dimP0 = shapes.polynomial_dimension( shape , 0 )

18

19

vec_P0 = vec_P1.take( reduce( lambda a,b:a+b , \

20

[ range(i*dimP1,i*dimP1+dimP0) \

21

for i in range(d) ] ) )

22

P1 = polynomial.OrthogonalPolynomialSet( shape , 1 )

23

P0H = P1[:dimP0]

24

25

Q = quadrature.make_quadrature( shape , 2 )

26

27

P0Hxcoeffs = Numeric.array( [ [ polynomial.projection( P1 , \

28

lambda x:x[i]*p(x), \

29

Q ).dof \

30

for i in range(d) ] \

31

for p in P0H ] )

32

33

P0Hx = polynomial.VectorPolynomialSet( P1.base , P0Hxcoeffs )

34

35

return polynomial.poly_set_union( vec_P0 , P0Hx )

36

37

38

# (P_k)^2 + x (P_k)

39

def RTSpace( shape , k ):

40

if k == 0:

41

return RT0Space( shape )

42

d = shapes.dimension( shape )

43

vec_Pkp1 = polynomial.OrthogonalPolynomialArraySet( shape , k+1 )

44

dimPkp1 = shapes.polynomial_dimension( shape , k+1 )

45

dimPk = shapes.polynomial_dimension( shape , k )

46

dimPkm1 = shapes.polynomial_dimension( shape , k-1 )

47

vec_Pk = vec_Pkp1.take( reduce( lambda a,b:a+b , \

48

[ range(i*dimPkp1,i*dimPkp1+dimPk) \

49

for i in range(d) ] ) )

50

Pkp1 = polynomial.OrthogonalPolynomialSet( shape , k + 1 )

51

PkH = Pkp1[dimPkm1:dimPk]

52

53

Q = quadrature.make_quadrature( shape , 2 * k )

54

55

PkHxcoeffs = Numeric.array( [ [ polynomial.projection( Pkp1 , \

56

lambda x:x[i]*p(x), \

57

Q ).dof \

58

for i in range(d) ] \

59

for p in PkH ] )

60

61

PkHx = polynomial.VectorPolynomialSet( Pkp1.base , PkHxcoeffs )

62

63

return polynomial.poly_set_union( vec_Pk , PkHx )

64

65

66

class RTDual( dualbasis.DualBasis ):

67

def __init__( self , shape , k , U ):

68

mdcb = functional.make_directional_component_batch

69

d = shapes.dimension( shape )

70

pts_per_edge = [ [ x \

71

for x in points.make_points( shape , \

72

d-1 , \

73

i , \

74

d+k ) ] \

75

for i in shapes.entity_range( shape , d-1 ) ]

76

nrmls = shapes.normals[shape]

77

ls = reduce( lambda a,b:a+b , \

78

[ mdcb(U,nrmls[i],pts_per_edge[i]) \

79

for i in shapes.entity_range(shape,d-1) ] )

80

81

if k > 0:

82

Pkp1 = polynomial.OrthogonalPolynomialArraySet( shape , k+1 )

83

dim_Pkp1 = shapes.polynomial_dimension( shape , k+1 )

84

dim_Pkm1 = shapes.polynomial_dimension( shape , k-1 )

85

86

Pkm1 = Pkp1.take( reduce( lambda a,b:a+b , \

87

[ range(i*dim_Pkp1,i*dim_Pkp1+dim_Pkm1) \

88

for i in range(d) ] ) )

89

90

91

interior_moments = [ functional.IntegralMoment( U , p ) \

92

for p in Pkm1 ]

93

94

ls.extend( interior_moments )

95

else:

96

interior_moments = []

97

98

99

entity_ids = {}

100

for i in range(d-1):

101

entity_ids[i] = {}

102

for j in shapes.entity_range(shape,i):

103

entity_ids[i][j] = []

104

pts_per_bdry = len(pts_per_edge[0])

105

entity_ids[d-1] = {}

106

node_cur = 0

107

for j in shapes.entity_range(shape,d-1):

108

for k in range(pts_per_bdry):

109

entity_ids[d-1][j] = node_cur

110

node_cur += 1

111

entity_ids[d] = range(node_cur,\

112

node_cur+len(interior_moments))

113

114

115

dualbasis.DualBasis.__init__( self , \

116

functionalset.FunctionalSet( U , ls ) , \

117

entity_ids )

118

119

class RT0( polynomial.FiniteElement ):

120

def __init__( self , shape ):

121

U = RT0Space( shape )

122

Udual = RTDual( shape , 0 , U )

123

polynomial.FiniteElement.__init__( self , Udual , U )

124

125

126

class RaviartThomas( polynomial.FiniteElement ):

127

def __init__( self , shape , n ):

128

U = RTSpace( shape , n )

129

Udual = RTDual( shape , n , U )

130

polynomial.FiniteElement.__init__( self , Udual , U )

131