~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 6 April 2005

8

9

"""points.py gives locations for edge, face, tetrahedron points

10

for each simplex shape.

11

The major functions to use are

12

make_lattice( shape , n )

13

and

14

make_points( shape , dim , entity_id , order )

15

These functions are used to put down a regular lattice over a simplex

16

or to extract the points associated with particular entities of a

17

particular topological dimension in the lattice of a particular

18

order. This is very useful in determining nodal locations for finite

19

element dual bases."""

20

21

import shapes

22

from curry import curry

23

from math import sqrt

24

25

26

rt2 = sqrt(2.0)

27

28

# these three functions are called within the function

29

# make_lattice (below)

30

def make_lattice_line(n):

31

return tuple( [ ( -1.0 + (2.0 * jj) / n , )\

32

for jj in range(0,n+1) ] )

33

34

def make_lattice_triangle(n):

35

return tuple( [ ( -1.0 + (2.0 * jj) / n , -1.0 + (2.0 * ii) / n ) \

36

for ii in xrange(0,n+1) \

37

for jj in xrange(0,n-ii+1)] )

38

39

def make_lattice_tetrahedron( n ):

40

return tuple( [ (-1.+(2.*jj)/n,-1.+(2.*ii)/n,-1.+(2.*kk)/n) \

41

for ii in xrange(0,n+1) \

42

for jj in xrange(0,n-ii+1) \

43

for kk in xrange(0,n-ii-jj+1) ] )

44

45

lattice_funcs = { shapes.LINE: make_lattice_line , \

46

shapes.TRIANGLE: make_lattice_triangle, \

47

shapes.TETRAHEDRON: make_lattice_tetrahedron }

48

49

def make_lattice( shape , n ):

50

global lattice_funcs

51

try:

52

return lattice_funcs[ shape ]( n )

53

except:

54

raise shapes.ShapeError, "Illegal shape: points.make_lattice"

55

56

57

# Functions that map a single point on the line to an edge

58

# on a triangle

59

pt_to_edge_triangle = { 0:lambda x: (-x[0],x[0]), \

60

1:lambda x: (-1.0,-x[0]), \

61

2:lambda x: (x[0],-1.0) }

62

63

64

tuple_pt_to_edge_triangle = { 0:lambda x: (-x[0],x[0]), \

65

1:lambda x: (-1.0,-x[0]), \

66

2:lambda x: (x[0],-1.0) }

67

68

tuple_pt_to_edge_tetrahedron = { 0: lambda x: (-x[0],x[0],-1.), \

69

1: lambda x: (-1.,-x[0],x[0]), \

70

2: lambda x: (x[0],-1.,-x[0]), \

71

3: lambda x: (-1.,-1.,-x[0]), \

72

4: lambda x: (-1.,x[0],-1.), \

73

5: lambda x: (x[0],-1.,-1.) }

74

75

76

tuple_pt_to_edge = { shapes.TRIANGLE: tuple_pt_to_edge_triangle , \

77

shapes.TETRAHEDRON: tuple_pt_to_edge_tetrahedron }

78

79

# Functions that map a single point on the line to an edge

80

# on a tet

81

pt_to_edge_tetrahedron = { 0: lambda x: (-x[0],x[0],-1.), \

82

1: lambda x: (-1.,-x[0],x[0]), \

83

2: lambda x: (x[0],-1.,-x[0]), \

84

3: lambda x: (-1.,-1.,-x[0]), \

85

4: lambda x: (-1.,x[0],-1.), \

86

5: lambda x: (x[0],-1.,-1.) }

87

88

# functions that take a point in the plane (ref triangle) and

89

# map it to the appropriate face of the ref tet

90

pt_to_face_tetrahedron = { 0: lambda x: (-x[0]-x[1]-1.,x[0],x[1]),

91

1: lambda x: (-1.,-x[0]-x[1]-1.,x[0]),

92

2: lambda x: (x[1],-1.,-x[0]-x[1]-1.), \

93

3: lambda x: (x[0],x[1],-1.) }

94

95

# dictionaries that dispatch on shape type

96

pt_to_edge = { shapes.TRIANGLE: pt_to_edge_triangle , \

97

shapes.TETRAHEDRON: pt_to_edge_tetrahedron }

98

99

pt_to_face = { shapes.TETRAHEDRON: pt_to_face_tetrahedron }

100

101

pt_maps = { shapes.LINE : { } , \

102

shapes.TRIANGLE: { 1 : pt_to_edge_triangle } , \

103

shapes.TETRAHEDRON: { 1 : pt_to_edge_tetrahedron , \

104

2 : pt_to_face_tetrahedron } }

105

106

107

108

# functions that get points of a particular order

109

# on a mesh component of particular dimension

110

# these are not needed in the public interface,

111

# which is simply the function make_points.

112

113

def make_vertex_points_line( vid , order ):

114

try:

115

return ( shapes.vertices[shapes.LINE][vid] , )

116

except:

117

raise RuntimeError, "Illegal vertex number"

118

119

def make_vertex_points_triangle( vid , order ):

120

global vertices

121

try:

122

return ( shapes.vertices[shapes.TRIANGLE][vid], )

123

except:

124

raise RuntimeError, "Illegal vertex number"

125

126

def make_vertex_points_tetrahedron( vid , order ):

127

global vertices

128

try:

129

return ( shapes.vertices[shapes.TETRAHEDRON][vid], )

130

except:

131

raise RuntimeError, "Illegal vertex number"

132

133

def make_edge_points_line( eid , order ):

134

if eid != 0:

135

raise RuntimeError, "Illegal edge number"

136

else:

137

return tuple( map(lambda x: (x,) ,line_range(order)) )

138

139

def make_edge_points_triangle( eid , order ):

140

global pt_to_edge_triangle

141

try:

142

return tuple( [ pt_to_edge_triangle[eid]( x ) \

143

for x in line_range( order ) ] )

144

except:

145

raise RuntimeError, "Illegal edge number."

146

147

def make_edge_points_tetrahedron( eid , order ):

148

global pt_to_edge_tetrahedron

149

try:

150

return tuple( [ pt_to_edge_tetrahedron[eid]( x ) \

151

for x in line_range( order ) ] )

152

except:

153

raise RuntimeError, "Illegal edge number"

154

155

def make_face_points_tetrahedron( fid , order ):

156

try:

157

return tuple( [ pt_to_face_tetrahedron[fid]( x ) \

158

for x in \

159

make_interior_points_triangle( 0 , order ) ] )

160

except:

161

raise RuntimeError, "Can't make triangle points"

162

163

164

def make_interior_points_triangle( iid , order ):

165

if iid != 0:

166

raise RuntimeError, "Illegal id"

167

return tuple( [ ( -1.0 + (2.0 * jj) / order , \

168

-1.0 + (2.0 * ii) / order ) \

169

for ii in range(1,order-1) \

170

for jj in range(1,order-ii)] )

171

172

def make_interior_points_tetrahedron( iid , n ):

173

if iid != 0:

174

raise RuntimeError, "Illegal id"

175

return tuple( [ (-1.+(2.*jj)/n,-1.+(2.*ii)/n,-1.+(2.*kk)/n) \

176

for ii in xrange(1,n) \

177

for jj in xrange(1,n-ii) \

178

for kk in xrange(1,n-ii-jj) ] )

179

180

181

def line_range( n ):

182

f = curry( lambda m,i : (-1.0 + (2.0*i)/m,) , n )

183

return tuple( map( f , range(1,n) ) )

184

185

186

#def line_range(n):

187

# f = curry( lambda m,i : -1.0 + (2.0*i)/m , n )

188

# return tuple( map( f , range(1,n) ) )

189

190

# dictionaries mapping dimension of mesh entities to

191

# functions for generating points of a particular order

192

# on those mesh entities

193

194

line_pts = { 0: make_vertex_points_line , \

195

1: make_edge_points_line }

196

197

tri_pts = { 0: make_vertex_points_triangle , \

198

1: make_edge_points_triangle , \

199

2: make_interior_points_triangle }

200

201

tet_pts = { 0: make_vertex_points_tetrahedron , \

202

1: make_edge_points_tetrahedron , \

203

2: make_face_points_tetrahedron , \

204

3: make_interior_points_tetrahedron }

205

206

pt_funcs = { shapes.LINE: line_pts , \

207

shapes.TRIANGLE: tri_pts , \

208

shapes.TETRAHEDRON: tet_pts }

209

210

def make_points( shape , dim , entity_id , order ):

211

"""Main public interface for getting points associated with mesh

212

components on reference domain. You give it a shape

213

( defined in shapes.py) the topological dimension of the

214

entity (0 for vertices, 1 for lines, etc), the id of the mesh

215

entity of that dimension, and the order of points, and it gives

216

you the appropriate tuple."""

217

global pt_funcs

218

try:

219

return pt_funcs[ shape ][ dim ]( entity_id , order )

220

except:

221

raise RuntimeError, "Can't make those points"

222

223

224

def barycentric_to_ref_tri( pt_bary ):

225

(L1,L2,L3) = pt_bary

226

x = -L1 + L2 - L3

227

y = -L1 - L2 + L3

228

return (x,y)

229