Coverage for src/methodsnm/fes.py: 55%

1from abc import ABC, abstractmethod 

2import numpy as np 

3from numpy import array 

4from methodsnm.fe_1d import * 

5from methodsnm.fe_2d import * 

6 

7class FESpace(ABC): 

8 """ 

9 Abstract base class for finite element spaces. 

10 """ 

11 ndof = None 

12 mesh = None 

13 def __init__(self, mesh): 

14 pass 

15 

16 @abstractmethod 

17 def _finite_element(self, elnr): 

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

19 

20 @abstractmethod 

21 def _element_dofs(self, elnr): 

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

23 

24 def _bndry_finite_element(self, elnr): 

25 raise Exception("_bndry_finite_element not implemented - Base class should not be used") 

26 

27 def _bndry_element_dofs(self, elnr): 

28 raise Exception("_bndry_element_dofs not implemented - Base class should not be used") 

29 

30 def finite_element(self, elnr, bndry=False): 

31 """ 

32 Returns the finite element for the given element number. 

33 """ 

34 if bndry: 

35 return self._bndry_finite_element(elnr) 

36 else: 

37 return self._finite_element(elnr) 

38 

39 def element_dofs(self, elnr, bndry=False): 

40 """ 

41 Returns the dofs for the given element number. 

42 """ 

43 if bndry: 

44 return self._bndry_element_dofs(elnr) 

45 else: 

46 return self._element_dofs(elnr) 

47 

48class VertexFirstSpace_1D(FESpace): 

49 def __init__(self, mesh): 

50 self.sfe = Node_FE() 

51 self.mesh = mesh 

52 

53 def _bndry_element_dofs(self, bndry_elnr): 

54 return self.mesh.elements(bndry=True)[bndry_elnr] 

55 

56 def _bndry_finite_element(self, bndry_elnr): 

57 return self.sfe 

58 

59class P1_Segments_Space(VertexFirstSpace_1D): 

60 

61 def __init__(self, mesh, periodic=False): 

62 super().__init__(mesh) 

63 self.periodic = periodic 

64 if periodic: 

65 self.ndof = len(mesh.points) - 1 

66 else: 

67 self.ndof = len(mesh.points) 

68 self.fe = P1_Segment_FE() 

69 

70 def _finite_element(self, elnr): 

71 return self.fe 

72 

73 def _element_dofs(self, elnr): 

74 dofs = self.mesh.edges[elnr] 

75 if self.periodic and elnr == len(self.mesh.edges) - 1: 

76 return [dofs[0],0] 

77 else: 

78 return dofs 

79 

80 def _bndry_element_dofs(self, bndry_elnr): 

81 dofs = self.mesh.elements(bndry=True)[bndry_elnr] 

82 if self.periodic: 

83 return [0] 

84 else: 

85 return dofs 

86 

87class P1disc_Segments_Space(FESpace): 

88 

89 def __init__(self, mesh): 

90 self.ndof = 2*len(mesh.edges) 

91 self.mesh = mesh 

92 self.fe = P1_Segment_FE() 

93 

94 def _finite_element(self, elnr): 

95 return self.fe 

96 

97 def _element_dofs(self, elnr): 

98 return [2*elnr, 2*elnr+1] 

99 

100 def _bndry_finite_element(self, bndry_elnr): 

101 raise Exception("No boundary evaluation for Discontinuous elements") 

102 

103 def _bndry_element_dofs(self, bndry_elnr): 

104 raise Exception("No boundary evaluation for Discontinuous elements") 

105 

106class Lagrange_Segments_Space(VertexFirstSpace_1D): 

107 

108 def __init__(self, mesh, order=1): 

109 super().__init__(mesh) 

110 self.nv = len(mesh.points) 

111 self.ne = len(mesh.edges) 

112 self.order = order 

113 self.ndof = self.nv + self.ne*(order-1) 

114 self.fe=Lagrange_Segment_FE(order=self.order) 

115 

116 def _finite_element(self, elnr): 

117 return self.fe 

118 

119 def _element_dofs(self, elnr): 

120 offset = self.nv + elnr*(self.order-1) 

121 dofs = [self.mesh.edges[elnr][0]] + [offset + i for i in range(0,self.order-1)] + [self.mesh.edges[elnr][1]] 

122 return dofs 

123 

124 

125 

126class Pk_IntLeg_Segments_Space(VertexFirstSpace_1D): 

127 

128 def __init__(self, mesh, order=1): 

129 super().__init__(mesh) 

130 self.nv = len(mesh.points) 

131 self.ne = len(mesh.edges) 

132 self.order = order 

133 self.ndof = self.nv + self.ne*(order-1) 

134 self.fe=IntegratedLegendre_Segment_FE(order=self.order) 

135 

136 def _finite_element(self, elnr): 

137 return self.fe 

138 

139 def _element_dofs(self, elnr): 

140 offset = self.nv + elnr*(self.order-1) 

141 dofs = [self.mesh.edges[elnr][0]] + [self.mesh.edges[elnr][1]] + [offset + i for i in range(0,self.order-1)] 

142 return dofs 

143 

144 

145 

146class P1_Triangle_Space(FESpace): 

147 

148 def __init__(self, mesh): 

149 self.ndof = len(mesh.points) 

150 self.mesh = mesh 

151 self.fe = P1_Triangle_FE() 

152 self.sfe = P1_Segment_FE() 

153 

154 def _finite_element(self, elnr): 

155 return self.fe 

156 

157 def _bndry_finite_element(self, bndry_elnr): 

158 return self.sfe 

159 

160 def _element_dofs(self, elnr): 

161 return self.mesh.elements()[elnr] 

162 

163 def _bndry_element_dofs(self, bndry_elnr): 

164 return self.mesh.elements(bndry=True)[bndry_elnr] 

165 

166 

167class P2_Triangle_Space(FESpace): 

168 

169 def __init__(self, mesh): 

170 self.ndof = len(mesh.points) + len(mesh.edges) 

171 self.nv = len(mesh.points) 

172 self.mesh = mesh 

173 self.fe = P2_Triangle_FE() 

174 self.sfe = Lagrange_Segment_FE(order=2) 

175 

176 def _finite_element(self, elnr): 

177 return self.fe 

178 

179 def _bndry_finite_element(self, bndry_elnr): 

180 return self.sfe 

181 

182 def _element_dofs(self, elnr): 

183 return np.append(self.mesh.elements()[elnr], [self.nv + i for i in self.mesh.faces2edges[elnr]]) 

184 

185 def _bndry_element_dofs(self, bndry_elnr): 

186 verts = self.mesh.elements(bndry=True)[bndry_elnr] 

187 edge = self.mesh.bndry_edges[bndry_elnr] 

188 return [verts[0], self.nv + edge, verts[1]] 

189 

190class P3_Triangle_Space(FESpace): 

191 

192 def __init__(self, mesh): 

193 self.ndof = len(mesh.points) + 2 * len(mesh.edges) + len(mesh.faces) 

194 self.nv = len(mesh.points) 

195 self.ned = len(mesh.edges) 

196 self.nf = len(mesh.faces) 

197 self.mesh = mesh 

198 self.fe = P3_Triangle_FE() 

199 self.sfe = Lagrange_Segment_FE(order=3) 

200 

201 def _finite_element(self, elnr): 

202 return self.fe 

203 

204 def _element_dofs(self, elnr): 

205 dofs = np.empty(10, dtype=int) 

206 vnums = self.mesh.elements()[elnr] 

207 dofs[0:3] = vnums 

208 j = 3 

209 enums = self.mesh.faces2edges[elnr] 

210 ref_verts_list = [[1,2],[0,2],[0,1]] 

211 for ref_edge, edge in enumerate(enums): 

212 tverts = vnums[ref_verts_list[ref_edge]] 

213 if (tverts[0] < tverts[1]): 

214 dofs[j] = self.nv + 2*edge + 0 

215 j += 1 

216 dofs[j] = self.nv + 2*edge + 1 

217 j += 1 

218 else: 

219 dofs[j] = self.nv + 2*edge + 1 

220 j += 1 

221 dofs[j] = self.nv + 2*edge + 0 

222 j += 1 

223 dofs[9] = self.nv + 2 * self.ned + elnr 

224 return dofs 

225 

226 def _bndry_finite_element(self, bndry_elnr): 

227 return self.sfe 

228 

229 def _bndry_element_dofs(self, bndry_elnr): 

230 verts = self.mesh.elements(bndry=True)[bndry_elnr] 

231 edge = self.mesh.bndry_edges[bndry_elnr] 

232 return [verts[0], self.nv + 2 * edge, self.nv + 2 * edge + 1, verts[1]] 

233 

234 

235class P1Edge_Triangle_Space(FESpace): 

236 

237 def __init__(self, mesh): 

238 self.ndof = len(mesh.edges) 

239 self.mesh = mesh 

240 self.fe = P1Edge_Triangle_FE() 

241 

242 def _finite_element(self, elnr): 

243 return self.fe 

244 

245 def _element_dofs(self, elnr): 

246 edges = self.mesh.faces2edges[elnr] 

247 return edges