Coverage for src/methodsnm/fe_2d.py: 90%

1from abc import ABC, abstractmethod 

2import numpy as np 

3from numpy import array 

4from methodsnm.fe import * 

5 

6class FE_2D(FE): 

7 """ 

8 Abstract base class for finite elements in 2D. 

9 It implements a derivative evaluation using numerical differentiation. 

10 """ 

11 

12 num_diff_warned = False 

13 

14 def __init__(self): 

15 self.dim = 2 

16 

17 def _evaluate_id(self, ip): 

18 raise Exception("Not implemented - Base class should not be used") 

19 

20 def _evaluate_deriv(self, ip: np.ndarray) -> np.ndarray: 

21 ip_ad = np.empty(2, dtype=AD) 

22 x: float 

23 y: float 

24 

25 if isinstance(ip, AD): 

26 raise ValueError() 

27 elif isinstance(ip, np.ndarray): 

28 if ip.ndim == 1: 

29 x,y = ip 

30 else: 

31 x,y = ip[0] 

32 

33 ip_ad[0] = AD(x, np.array([1.,0.])) 

34 ip_ad[1] = AD(y, np.array([0.,1.])) 

35 ip_deriv = np.array(list(map(lambda f_x: f_x.deriv, self._evaluate_id(ip_ad)))).T 

36 return ip_deriv 

37 

38 def _evaluate_deriv_numerical(self, ip): 

39 # numerical differentiation - should be overwritten by subclasses 

40 # for proper accuracy and performance 

41 if not FE_2D.num_diff_warned: 

42 print("Warning: Using numerical differentiation for deriv evaluation in " + str(type(self)) + " object.") 

43 FE_2D.num_diff_warned = True 

44 eps = 1e-8 

45 left = ip.copy() - array([eps,0]) 

46 right = ip.copy() + array([eps,0]) 

47 top = ip.copy() + array([0,eps]) 

48 bottom = ip.copy() - array([0,eps]) 

49 ret = array([(self._evaluate_id(right) - self._evaluate_id(left))/(2*eps), 

50 (self._evaluate_id(top) - self._evaluate_id(bottom))/(2*eps)]) 

51 return ret 

52 

53 

54class TriangleFE(FE_2D): 

55 """ 

56 Abstract base class for finite elements on triangles. 

57 """ 

58 

59 def __init__(self): 

60 super().__init__() 

61 self.eltype = "triangle" 

62 

63 @abstractmethod 

64 def _evaluate_id(self, ip): 

65 raise Exception("Not implemented - Base class should not be used") 

66 

67class P1_Triangle_FE(TriangleFE, Lagrange_FE): 

68 """ 

69 This class represents a P1 triangle finite element. 

70 """ 

71 ndof = 3 

72 order = 1 

73 def __init__(self): 

74 super().__init__() 

75 self.nodes = [ np.array([0,0]), np.array([1,0]), np.array([0,1]) ] 

76 

77 def _evaluate_id(self, ip): 

78 """ 

79 Evaluates the P1 triangle finite element at the given integration point. 

80 

81 Parameters: 

82 ip (numpy.ndarray): The integration point at which to evaluate the finite element. 

83 

84 Returns: 

85 numpy.ndarray: The values of the P1 triangle finite element at the given integration point. 

86 """ 

87 if isinstance(ip, AD): 

88 raise ValueError() 

89 elif isinstance(ip, np.ndarray): 

90 if ip.ndim == 1: 

91 x,y = ip 

92 else: 

93 x,y = ip[0] 

94 

95 return array([1-x - y, x, y]) 

96 

97 def __str__(self): 

98 return "P1 Triangle Finite Element" 

99 

100 

101class P2_Triangle_FE(TriangleFE, Lagrange_FE): 

102 """ 

103 This class represents a P2 triangle (Lagrange) finite element. 

104 """ 

105 ndof = 6 

106 order = 2 

107 def __init__(self): 

108 super().__init__() 

109 self.nodes = [ np.array([0,0]), np.array([1,0]), np.array([0,1]), 

110 np.array([0.5,0.5]), np.array([0,0.5]), np.array([0.5,0]) ] 

111 

112 def _evaluate_id(self, ip): 

113 """ 

114 Evaluates the P2 triangle finite element at the given integration point. 

115 

116 Parameters: 

117 ip (numpy.ndarray): The integration point at which to evaluate the finite element. 

118 

119 Returns: 

120 numpy.ndarray: The values of the P2 triangle finite element at the given integration point. 

121 """ 

122 ret = np.empty(self.ndof, dtype=ip.dtype) 

123 

124 if isinstance(ip, AD): 

125 raise ValueError() 

126 elif isinstance(ip, np.ndarray): 

127 if ip.ndim == 1: 

128 x,y = ip 

129 else: 

130 x,y = ip[0] 

131 

132 lamb = [1-x-y,x,y] 

133 

134 for i in range(3): 

135 ret[i] = lamb[i]*(2*lamb[i]-1) 

136 for i in range(3): 

137 for j in range(i): 

138 ret[6-i-j] = lamb[i]*lamb[j]*4 

139 return ret 

140 

141 def __str__(self): 

142 return "P2 Triangle Finite Element\n" + super().__str__() 

143 

144 

145class P3_Triangle_FE(TriangleFE, Lagrange_FE): 

146 """ 

147 This class represents a P3 triangle (Lagrange) finite element. 

148 (first draft : Warning: we may want to change some things later on!) 

149 """ 

150 ndof = 10 

151 order = 3 

152 def __init__(self): 

153 super().__init__() 

154 self.nodes = [ np.array([0,0]), np.array([1,0]), np.array([0,1]), 

155 np.array([2/3,1/3]), np.array([1/3,2/3]), np.array([0,1/3]), np.array([0,2/3]), np.array([1/3,0]), np.array([2/3,0]), 

156 np.array([1/3,1/3]) ] 

157 

158 def _evaluate_id(self, ip): 

159 """ 

160 Evaluates the P3 triangle finite element at the given integration point. 

161 

162 Parameters: 

163 ip (numpy.ndarray): The integration point at which to evaluate the finite element. 

164 

165 Returns: 

166 numpy.ndarray: The values of the P3 triangle finite element at the given integration point. 

167 """ 

168 ret = np.empty(self.ndof, dtype = ip.dtype) 

169 

170 if isinstance(ip, AD): 

171 raise ValueError() 

172 elif isinstance(ip, np.ndarray): 

173 if ip.ndim == 1: 

174 x,y = ip 

175 else: 

176 x,y = ip[0] 

177 

178 lamb = [1-x-y,x,y] 

179 

180 #vertex dofs: 

181 for i in range(3): 

182 ret[i] = 9.0/2.0*lamb[i]*(lamb[i]-1/3)*(lamb[i]-2/3) 

183 

184 #edge dofs: 

185 for i in range(1,3): # major vertex  

186 for j in range(i): # minor vertex 

187 edge_number = 3-i-j 

188 ret[3+2*edge_number] = lamb[i]*lamb[j]*(lamb[j]-1/3)*27/2 # node next to minor vertex 

189 ret[3+2*edge_number+1] = lamb[i]*lamb[j]*(lamb[i]-1/3)*27/2 # node next to major vertex 

190 ret[9] = 27*lamb[0]*lamb[1]*lamb[2] 

191 return ret 

192 

193 def __str__(self): 

194 return "P3 Triangle Finite Element\n" + super().__str__() 

195 

196 

197 

198 

199class P1Edge_Triangle_FE(TriangleFE, Lagrange_FE): 

200 """ 

201 This class represents a P1 triangle finite element. 

202 """ 

203 ndof = 3 

204 order = 1 

205 def __init__(self): 

206 super().__init__() 

207 self.nodes = [ np.array([0.5,0.5]), np.array([0,0.5]), np.array([0.5,0]) ] 

208 

209 def _evaluate_id(self, ip): 

210 """ 

211 Evaluates the P1Edge triangle finite element at the given integration point. 

212 

213 Parameters: 

214 ip (numpy.ndarray): The integration point at which to evaluate the finite element. 

215 

216 Returns: 

217 numpy.ndarray: The values of the P1Edge triangle finite element at the given integration point. 

218 """ 

219 if isinstance(ip, AD): 

220 raise ValueError() 

221 elif isinstance(ip, np.ndarray): 

222 if ip.ndim == 1: 

223 x,y = ip 

224 else: 

225 x,y = ip[0] 

226 

227 lamb = [1-x-y,x,y] 

228 return array([1-2*lamb[i] for i in range(3)]) 

229 

230 def __str__(self): 

231 return "P1Edge Triangle Finite Element" 

232 

233