Coverage for src/methodsnm/fe_1d.py: 85%

1from abc import ABC, abstractmethod 

2import numpy as np 

3from numpy import array 

4from methodsnm.fe import * 

5from methodsnm.diff import AD 

6 

7class FE_1D(FE): 

8 """ 

9 Abstract base class for finite elements in 1D. 

10 It implements a derivative evaluation using numerical differentiation. 

11 """ 

12 num_diff_warned = False 

13 

14 def __init__(self): 

15 self.eltype = "segment" 

16 self.dim = 1 

17 

18 @abstractmethod 

19 def _evaluate_id(self, ip): 

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

21 

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

23 ip_ad = np.empty(1, dtype=AD) 

24 ip_ad[0] = AD(ip[0], np.array([1])) 

25 ip_deriv = np.array(list(map(lambda f_x: f_x.deriv, self._evaluate_id(ip_ad))), dtype=np.float64) 

26 return ip_deriv.reshape((1,self.ndof)) 

27 

28 def _evaluate_deriv_numerical(self, ip): 

29 # numerical differentiation - should be overwritten by subclasses 

30 # for proper accuracy and performance 

31 if not FE_1D.num_diff_warned: 

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

33 FE_1D.num_diff_warned = True 

34 eps = 1e-8 

35 left = ip.copy() - eps 

36 right = ip.copy() + eps 

37 return ((self._evaluate_id(right) - self._evaluate_id(left))/(2*eps)).reshape((1,self.ndof)) 

38 

39class P1_Segment_FE(FE_1D, Lagrange_FE): 

40 """ 

41 This class represents a P1 segment finite element. 

42 """ 

43 ndof = 2 

44 order = 1 

45 def __init__(self): 

46 super().__init__() 

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

48 

49 def _evaluate_id(self, ip): 

50 """ 

51 Evaluates the P1 segment finite element at the given integration point. 

52 

53 Parameters: 

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

55 

56 Returns: 

57 numpy.ndarray: The values of the P1 segment finite element at the given integration point. 

58 """ 

59 return array([1-ip[0], ip[0]]) 

60 

61 def __str__(self): 

62 return "P1 Segment Finite Element\n" + super().__str__() 

63 

64 #def _evaluate_deriv(self, ip): 

65 # return np.full(shape=ip.shape + (2,), fill_value=[-1,1]) 

66 

67class P2_Segment_FE(FE_1D): 

68 """ 

69 This class represents a P1 segment finite element. 

70 """ 

71 ndof = 3 

72 order = 2 

73 def __init__(self): 

74 super().__init__() 

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

76 

77 def _evaluate_id(self, ip): 

78 """ 

79 Evaluates the P2 segment 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 P2 segment finite element at the given integration point. 

86 """ 

87 return array([1-ip[0], ip[0], 4*ip[0]*(1-ip[0])]) 

88 

89 def __str__(self): 

90 return "P2 Segment Finite Element\n" + super().__str__() 

91 

92 #def _evaluate_deriv(self, ip): 

93 # deriv = np.zeros(ip.shape + (3,)) 

94 # deriv[:, 0] = -1 

95 # deriv[:, 1] = 1 

96 # deriv[:, 2] = 4 - 8 * ip 

97 # return deriv 

98 

99class Lagrange_Segment_FE(Lagrange_FE, FE_1D): 

100 """ 

101 This class represents a Lagrange finite element on [0,1]. 

102 """ 

103 def __init__(self, order, nodes=None): 

104 super().__init__() 

105 self.order = order 

106 self.ndof = order+1 

107 if nodes is not None: 

108 if len(nodes) != self.ndof: 

109 raise Exception("Invalid number of nodes") 

110 self.nodes = nodes 

111 else: 

112 self.nodes = [ np.array(x) for x in np.linspace(0, 1, self.ndof) ] 

113 self.barycentric_weights = np.ones(self.ndof) 

114 for i in range(self.ndof): 

115 for j in range(self.ndof): 

116 if i != j: 

117 self.barycentric_weights[i] /= (self.nodes[i] - self.nodes[j]) 

118 

119 def _evaluate_id(self, ip): 

120 """ 

121 Evaluates the Lagrange segment finite element at the given integration point. 

122 

123 Uses the barycentric form of the Lagrange polynomials,  

124 see https://en.wikipedia.org/wiki/Lagrange_polynomial#Barycentric_form 

125 

126 l_j(x) = prod_{i!=j} (x-x_i)/(x_j-x_i)  

127 = w_j * prod_{i!=j} (x-x_i) with w_j = prod_i (1/(x_j-x_i)) 

128 = w_j / (x-x_j) * l(x) with l(x) = prod_i (x-x_i) 

129 With further 

130 1 = sum_i l_i(x) = sum_i (w_i / (x-x_i)) * l(x)  

131 = l(x) * sum_i (w_i / (x-x_i)) 

132 we have 

133 l_j(x) = w_j / (x-x_j) * sum_i (w_i / (x-x_i)) 

134 where the last sum is a does not depend on j. 

135 

136 Evaluation costs are hence O(ndof) instead of O(ndof^2) for the naive approach. 

137 

138 Parameters: 

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

140 

141 Returns: 

142 numpy.ndarray: The values of the Lagrange segment finite element at the given integration point. 

143 """ 

144 if ip[0] in self.nodes: 

145 ret = np.full(self.ndof, ip[0] - ip[0], dtype=ip.dtype) 

146 ret[self.nodes.index(ip[0])] = ip[0]**0 

147 else: 

148 denom = sum([self.barycentric_weights[i]/(ip[0]-self.nodes[i]) for i in range(self.ndof)]) 

149 ret = self.barycentric_weights.copy() 

150 if isinstance(ip[0], AD): 

151 ret = np.array(list(map(lambda x: AD(x, np.array([0.])), ret))) 

152 for i in range(self.ndof): 

153 ret[i] /= (ip[0]-self.nodes[i]) * denom 

154 return ret 

155 

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

157 """ 

158 Fall back to numerical differentiation, 

159 because the AD method is not suitable for the definition of 

160 the Lagrange polynomials. 

161 """ 

162 return super()._evaluate_deriv_numerical(ip) 

163 

164 def __str__(self): 

165 return f"Lagrange Segment Finite Element(order={self.order})\n" + super().__str__() 

166 

167from methodsnm.recpol import * 

168class RecPol_Segment_FE(FE_1D): 

169 """ 

170 This class represents a Recursive Polynomial finite element on [0,1]. 

171 """ 

172 def __init__(self, order, recpol): 

173 super().__init__() 

174 self.order = order 

175 self.ndof = order+1 

176 self.recpol = recpol 

177 

178 def _evaluate_id(self, ip): 

179 return self.recpol.evaluate_all(2*ip-1, self.order).reshape(-1) 

180 

181 def _evaluate_id_array(self, ip): 

182 return self.recpol.evaluate_all(2*ip[:,0]-1, self.order) 

183 

184 def __str__(self): 

185 return f"RecPol Segment Finite Element(recpol={self.recpol}, order={self.order})\n" + super().__str__() 

186 

187def Legendre_Segment_FE(order): 

188 return RecPol_Segment_FE(order, LegendrePolynomials()) 

189 

190def Jacobi_Segment_FE(order, alpha, beta): 

191 return RecPol_Segment_FE(order, JacobiPolynomials(alpha, beta)) 

192 

193 

194class IntegratedLegendre_Segment_FE(FE_1D): 

195 """ 

196 This class represents a finite element on [0,1] 

197 that combines the lowest order P1 element with 

198 integrated Legendre polynomials of higher order. 

199 """ 

200 def __init__(self, order): 

201 super().__init__() 

202 self.order = order 

203 self.ndof = order+1 

204 self.recpol = IntegratedLegendrePolynomials() 

205 

206 def _evaluate_id(self, ip): 

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

208 ret[0:2] = np.array([1-ip[0],ip[0]]) 

209 ret[2::] = self.recpol.evaluate_all(2*ip-1, self.order-1)[0,1::] 

210 return ret 

211 

212 def _evaluate_id_array(self, ip): 

213 ret = np.empty((len(ip),self.ndof)) 

214 ret[:,0] = 1-ip[:,0] 

215 ret[:,1] = ip[:,0] 

216 ret[:,2::] = self.recpol.evaluate_all(2*ip[:,0]-1, self.order-1)[:,1::] 

217 return ret 

218 

219 def __str__(self): 

220 return f"RecPol Segment Finite Element(recpol={self.recpol}, order={self.order})\n" + super().__str__()