~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 7 Feb 2005 by RCK

8

9

"""shapes.py

10

Module defining topological information for simplices in one, two, and

11

three dimensions."""

12

from math import sqrt

13

import exceptions

14

from factorial import factorial

15

16

rt2 = sqrt(2.0)

17

rt3 = sqrt(3.0)

18

rt2_inv = 1./rt2

19

rt3_inv = 1./rt3

20

21

22

class ShapeError( exceptions.Exception ):

23

"""FIAT defined exception type indicating an improper use of a shape."""

24

def __init__( self , args = None ):

25

self.args = args

26

return

27

28

# enumerated constants for shapes

29

LINE = 1

30

TRIANGLE = 2

31

TETRAHEDRON = 3

32

33

# the dimension associated with each shape

34

dims = { LINE: 1 , \

35

TRIANGLE: 2 , \

36

TETRAHEDRON: 3 }

37

38

# number of topological entities of each dimensions

39

# associated with each shape.

40

# 0 refers to the number of points that make up the simplex

41

# 1 refers to the number of edges

42

# 2 refers to the number of triangles/faces

43

# 3 refers to the number of tetrahedra

44

num_entities = { LINE: { 0: 2 , 1: 1 } , \

45

TRIANGLE: { 0: 3, 1: 3, 2:1 } , \

46

TETRAHEDRON: {0: 4, 1: 6, 2: 4, 3: 1} }

47

48

# canonical polynomial dimension associated with simplex

49

# of each dimension

50

poly_dims = { LINE: lambda n: max(0,n + 1), \

51

TRIANGLE: lambda n: max(0,(n+1)*(n+2)/2) ,\

52

TETRAHEDRON: lambda n: max(0,(n+1)*(n+2)*(n+3)/6) }

53

54

55

# dictionaries of vertices

56

vertices = { LINE : {0: (-1.0,) , 1: (1.0,)} , \

57

TRIANGLE : {0:(-1.0,-1.0),\

58

1:(1.0,-1.0),\

59

2:(-1.0,1.0)} , \

60

TETRAHEDRON :{0:(-1.,-1.,-1.), \

61

1:(1.,-1.,-1.), \

62

2:(-1.,1.,-1.), \

63

3:(-1.,-1.,1.)} }

64

65

66

# FYI -- this is the ordering scheme I'm thinking of.

67

# Should be verified against tetgen or

68

# whatever mesh generator we wind up using

69

tetrahedron_faces = {0:(1,2,3), \

70

1:(0,3,2), \

71

2:(0,1,3), \

72

3:(0,2,1)}

73

74

tetrahedron_edges = {0:(1,2), \

75

1:(2,3), \

76

2:(3,1), \

77

3:(3,0), \

78

4:(0,2), \

79

5:(0,1)}

80

81

# scaling factors mapping boundary entities of lower dimension to the reference

82

# shape of that dimension. For example, in 2d, two of the edges have length

83

# 2, which is the same length as the reference line. The third is 2 sqrt(2),

84

# so we get two scale factors of 1 and one of sqrt(2.0)

85

jac_factors = { LINE: { 0: { 0: 1.0, \

86

1: 1.0 } , \

87

1: { 0: 1.0 } } , \

88

TRIANGLE: { 0: { 0: 1.0 , \

89

1: 1.0 , \

90

2: 1.0 } , \

91

1: { 0: rt2 , \

92

1: 1.0 , \

93

2: 1.0 } , \

94

2: { 0: 1.0 } } , \

95

TETRAHEDRON : { 0: { 0: 1.0 ,

96

1: 1.0 ,

97

2: 1.0 ,

98

3: 1.0 } , \

99

1: { 0: rt2 , \

100

1: rt2 , \

101

2: rt2 , \

102

3: 1.0 , \

103

4: 1.0 , \

104

5: 1.0 } , \

105

2: { 0: rt3 , \

106

1: 1.0 , \

107

2: 1.0 , \

108

3: 1.0 } ,

109

3: { 0: 1.0 } } }

110

111

112

# dictionary of shapes -- for each shape, maps edge number to normal

113

# which is a tuple of floats.

114

normals = { LINE : { 0: (-1.0,) , 1:(1.0,) } , \

115

TRIANGLE : { 0: (rt2_inv,rt2_inv) , \

116

1: (-1.0,0) , \

117

2: (0,-1.0) } , \

118

TETRAHEDRON : { 0: (rt3_inv,rt3_inv,rt3_inv) , \

119

1: (-1.0,0.0,0.0) , \

120

2: (0.0,-1.0,0.0) , \

121

3: (0.0,0.0,-1.0) } }

122

123

tangents = { TETRAHEDRON : { 0: (-rt2_inv,rt2_inv,0.0) , \

124

1: (0.0,-rt2_inv,rt2_inv) , \

125

2: (rt2_inv,0.0,-rt2_inv) , \

126

3: (0.0,0.0,-1.0) , \

127

4: (0.0,1.0,0.0) ,\

128

5: (1.0,0.0,0.0) } }

129

130

def scale_factor( shape , d , ent_id ):

131

global jac_factors

132

return jac_factors[ shape ][ d ][ ent_id ]

133

134

def dimension( shape ):

135

"""returns the topological dimension associated with shape."""

136

global dims

137

try:

138

return dims[ shape ]

139

except:

140

raise ShapeError, "Illegal shape: shapes.dimension"

141

142

def dimension_range( shape ):

143

"""returns the list starting at zero and ending with

144

the topological dimension of shape (inclusive).

145

Hence, dimension_range( s ) is syntactic sugar for

146

range(0,dimension(shape)+1)."""

147

return range(0,dimension(shape)+1);

148

149

def entity_range( shape , dim ):

150

"""Returns the range of topological entities of dimension dim

151

associated with shape.

152

For example, entity_range( LINE , 0 ) returns the list [0,1]

153

because there are two points associated with the line."""

154

global num_entities

155

try:

156

return range(0,num_entities[ shape ][ dim ])

157

except:

158

raise ShapeError, "Illegal shape or dimension"

159

160

def polynomial_dimension( shape , n ):

161

"""Returns the number of polynomials of total degree n on the

162

shape n. This (n+1) over the line, (n+1)(n+2)/2 over the

163

triangle, and (n+1)(n+2)(n+3)/6 in three dimensions."""

164

d = dimension( shape )

165

td = 1

166

for i in xrange(0,d):

167

td = td * ( n + i + 1 )

168

return td / factorial( d )

169