Coverage for readers/gmsh/gmsh.py: 15%

200 statements  

« prev     ^ index     » next       coverage.py v7.7.0, created at 2025-03-20 20:51 +0100

1import logging 

2from collections import defaultdict, namedtuple 

3from functools import cached_property 

4from pathlib import Path 

5 

6import numpy as np 

7 

8from nastranio.readers.gmsh.mixins import Modificator 

9 

10Entity = namedtuple("Entity", ["entityDim", "entityTag", "bbox", "physical_tags"]) 

11PhysicalGroup = namedtuple("PhysicalGroup", ["dim", "physical_tag", "name"]) 

12Node = namedtuple("Node", ["gid", "entity", "x", "y", "z"]) 

13Element = namedtuple( 

14 "Element", 

15 [ 

16 "eid", 

17 "gids", 

18 "entity_dim", 

19 "entity_tag", 

20 "element_type", 

21 "physical_tags", 

22 "physical_groups", 

23 ], 

24) 

25 

26 

27def revert_dict(d): 

28 newd = defaultdict(set) 

29 # if `d` is kind of {k: (iterable)} 

30 if isinstance(next(iter(d.values())), (tuple, list, set)): 

31 for k, v in d.items(): 

32 for subv in v: 

33 newd[subv].add(k) 

34 else: 

35 for k, v in d.items(): 

36 newd[v] = k 

37 return dict(newd) 

38 

39 

40class GmsHParser(Modificator): 

41 def __init__(self, filepath): 

42 self.read(filepath) 

43 

44 def info(self): 

45 msg = "" 

46 msg += f"Mesh Format {self._mesh_format}\n" 

47 msg += f"nb. nodes: {len(self.nodes)}\n" 

48 msg += f"nb. elements: {len(self.elements)}\n" 

49 cnt_types = {k: len(v) for k, v in self.type2eid.items()} 

50 for typ, nb in cnt_types.items(): 

51 msg += f" Type {typ}: {nb} elements\n" 

52 coincident_gids = self.coincident_nodes() 

53 if coincident_gids: 

54 msg += f"Coincident gids: {coincident_gids}\n" 

55 return msg 

56 

57 def coincident_nodes(self, digits=0): 

58 xyzs = np.array([[n.gid, n.x, n.y, n.z] for gid, n in self.nodes.items()]) 

59 gids = xyzs[:, 0].astype(int) 

60 xyzs = xyzs[:, 1:] 

61 if digits: 

62 xyzs = np.round(xyzs, digits) 

63 dedup_values, count = np.unique(xyzs, axis=0, return_counts=True) 

64 ix2gid = dict(zip(range(len(gids)), gids)) 

65 repeated_groups = dedup_values[count > 1] 

66 coincident_gids = [] 

67 for repeated_group in repeated_groups: 

68 repeated_idx = np.argwhere(np.all(xyzs == repeated_group, axis=1)).ravel() 

69 repeated_idx = (ix2gid[ix] for ix in repeated_idx) 

70 coincident_gids.append(tuple(repeated_idx)) 

71 return coincident_gids 

72 

73 def read(self, filepath): 

74 with open(filepath) as fh: 

75 self._toc = {} 

76 self._build_toc(fh.readlines()) 

77 for section in self._toc.keys(): 

78 process_func = f"_process_{section}" 

79 if hasattr(self, process_func): 

80 getattr(self, process_func)() 

81 else: 

82 logging.warning(f"no processing function for {section}") 

83 self._make_gid_eid_type() 

84 self._make_eid_gid_physicalGroups() 

85 

86 # ========================================================================= 

87 # a few helpers 

88 # ========================================================================= 

89 # @cached_property 

90 # def type2eid(self): 

91 # eid2type = defaultdict(set) 

92 # for eid, elt in self.elements.items(): 

93 # eid2type[elt.element_type].add(eid) 

94 # return dict(eid2type) 

95 

96 @cached_property 

97 def dim2eid(self): 

98 eid2type = defaultdict(set) 

99 for eid, elt in self.elements.items(): 

100 eid2type[elt.entity_dim].add(eid) 

101 return dict(eid2type) 

102 

103 @cached_property 

104 def dim2physical_group(self): 

105 dim2pg = defaultdict(set) 

106 for (dim, _id), pg in self.physical_groups.items(): 

107 dim2pg[pg.dim].add(pg) 

108 return dict(dim2pg) 

109 

110 def _make_gid_eid_type(self): 

111 """build: 

112 

113 * self.eid2gids 

114 * self.gid2eids 

115 * self.type2eid 

116 * self.eid2type 

117 

118 """ 

119 eid2gids = {} 

120 gid2eids = defaultdict(set) 

121 type2eid = defaultdict(set) 

122 eid2type = {} 

123 for eid, elt in self.elements.items(): 

124 eid2gids[eid] = set(elt.gids) 

125 type2eid[elt.element_type].add(eid) 

126 eid2type[eid] = elt.element_type 

127 for gid in elt.gids: 

128 gid2eids[gid].add(eid) 

129 self.eid2gids = eid2gids 

130 self.type2eid = dict(type2eid) 

131 self.gid2eids = dict(gid2eids) 

132 self.eid2type = eid2type 

133 

134 def _make_eid_gid_physicalGroups(self): 

135 """build: 

136 

137 * self.eid2physical_groups 

138 * self.physical_group2eids 

139 * self.physical_group2gids 

140 * self.gid2physical_groups 

141 

142 """ 

143 self.eid2physical_groups = { 

144 eid: set(elt.physical_groups) for eid, elt in self.elements.items() 

145 } 

146 physical_group2eids = defaultdict(set) 

147 gid2physical_groups = defaultdict(set) 

148 physical_group2gids = defaultdict(set) 

149 for eid, grps in self.eid2physical_groups.items(): 

150 gids = self.eid2gids[eid] 

151 for grp in grps: 

152 physical_group2eids[grp].add(eid) 

153 physical_group2gids[grp] |= gids 

154 for gid in gids: 

155 gid2physical_groups[gid] |= grps 

156 self.physical_group2eids = dict(physical_group2eids) 

157 self.physical_group2gids = dict(physical_group2gids) 

158 self.gid2physical_groups = dict(gid2physical_groups) 

159 

160 # ========================================================================= 

161 # internal stuff 

162 # ========================================================================= 

163 

164 def _build_toc(self, content): 

165 buffer = [] 

166 IN_SECTION = None 

167 for line in content: 

168 line = line.split("//")[0].strip() 

169 if line.startswith("//"): 

170 continue 

171 if line.startswith("$End"): 

172 assert line[4:] == IN_SECTION 

173 self._toc[IN_SECTION] = buffer[:] 

174 buffer = [] 

175 continue 

176 if line.startswith("$"): 

177 IN_SECTION = line[1:] 

178 continue 

179 buffer.append(line) 

180 

181 def _process_MeshFormat(self): 

182 """creates self.mesh_format""" 

183 self._mesh_format = self._toc["MeshFormat"][0].split()[0] 

184 

185 def _process_PhysicalNames(self): 

186 """ 

187 creates `self.physical_groups` mapping (dim, physical tag) to physicalGroup named tuple 

188 """ 

189 content = self._toc["PhysicalNames"] 

190 nb_physical_names = int(content.pop(0)) 

191 data = {} 

192 for line in content: 

193 dim, physical_tag, *name = line.split() 

194 dim = int(dim) 

195 name = " ".join(name) 

196 _ = PhysicalGroup(dim, int(physical_tag), name.strip('"')) 

197 data[(dim, _.physical_tag)] = _ 

198 self.physical_groups = data 

199 

200 def _process_Entities_1d2d3d(self, nb, content, dim): 

201 for _ in range(nb): 

202 itemdata = content.pop(0) 

203 if dim == 0: 

204 entityTag, minX, minY, minZ, *groups = itemdata.split() 

205 bbox = (minX, minY, minZ) 

206 else: 

207 ( 

208 entityTag, 

209 minX, 

210 minY, 

211 minZ, 

212 maxX, 

213 maxY, 

214 maxZ, 

215 *groups, 

216 ) = itemdata.split() 

217 bbox = (minX, minY, minZ, maxX, maxY, maxZ) 

218 entityTag = int(entityTag) 

219 # processing numPhysicalTags 

220 numPhysicalTags = int(groups.pop(0)) 

221 physicalTags = groups[:numPhysicalTags] 

222 assert numPhysicalTags == len(physicalTags) 

223 groups = groups[numPhysicalTags:] 

224 if dim == 0: 

225 assert len(groups) == 0 

226 else: 

227 # should be finished here with 0d 

228 numBounding = int(groups[0]) 

229 pointTags = list(map(int, groups[1:])) 

230 assert len(pointTags) == numBounding 

231 entity = Entity( 

232 entityDim=dim, 

233 entityTag=entityTag, 

234 bbox=tuple(map(float, bbox)), 

235 physical_tags=tuple(map(int, physicalTags)), 

236 ) 

237 self.entities[dim][entityTag] = entity 

238 return content 

239 

240 def _process_Entities(self): 

241 """ 

242 creates a `self.entities` dictionnary, mapping 

243 {0: {entityTag: {"bbox": (0, 0, 0), "physicalTags": []}}} 

244 """ 

245 content = self._toc["Entities"] 

246 self.entities = {k: defaultdict(dict) for k in range(4)} 

247 nb = content.pop(0).split() 

248 # ---------------------------------------------------------------- 

249 # processing points 

250 numPoints, numCurves, numSurfaces, numVolumes = map(int, nb) 

251 content = self._process_Entities_1d2d3d(numPoints, content, 0) 

252 content = self._process_Entities_1d2d3d(numCurves, content, 1) 

253 content = self._process_Entities_1d2d3d(numSurfaces, content, 2) 

254 content = self._process_Entities_1d2d3d(numVolumes, content, 3) 

255 self.entities = {dim: dict(data) for dim, data in self.entities.items()} 

256 

257 def _process_Nodes(self): 

258 content = self._toc["Nodes"] 

259 numEntityBlocks, numNodes, minNodeTag, maxNodeTag = map( 

260 int, content.pop(0).split() 

261 ) 

262 self.nodes = {} 

263 

264 for blockno in range(numEntityBlocks): 

265 entityDim, entityTag, parametric, numNodesInBlock = map( 

266 int, content.pop(0).split() 

267 ) 

268 if numNodesInBlock == 0: 

269 continue 

270 node_ids = tuple(map(int, ",".join(content[:numNodesInBlock]).split(","))) 

271 xyzs = "\n".join(content[numNodesInBlock : numNodesInBlock * 2]) 

272 content = content[2 * numNodesInBlock :] 

273 xyzs = np.fromstring(xyzs, dtype=float, sep=" ") 

274 xyzs = xyzs.reshape(int(len(xyzs) / 3), 3) 

275 for i, gid in enumerate(node_ids): 

276 xyz = xyzs[i] 

277 node = Node( 

278 gid=gid, 

279 entity=self.entities[entityDim][entityTag], 

280 x=xyz[0], 

281 y=xyz[1], 

282 z=xyz[2], 

283 ) 

284 self.nodes[gid] = node 

285 

286 def _process_Elements(self): 

287 self.elements = {} 

288 # --------------------------------------------------------------------- 

289 # parse content 

290 content = self._toc["Elements"] 

291 numEntityBlocks, numElements, minElementTag, maxElementTag = map( 

292 int, content.pop(0).split() 

293 ) 

294 

295 for blockno in range(numEntityBlocks): 

296 entityDim, entityTag, elementType, numElementsInBlock = map( 

297 int, content.pop(0).split() 

298 ) 

299 physicalTags = tuple(self.entities[entityDim][entityTag].physical_tags) 

300 try: 

301 physicalNames = tuple( 

302 ( 

303 self.physical_groups[(entityDim, id_)].name 

304 for id_ in physicalTags 

305 ) 

306 ) 

307 except KeyError: 

308 # GMSH V4.12 API changes??? 

309 physicalNames = tuple( 

310 ( 

311 self.physical_groups[(entityDim, -id_)].name 

312 for id_ in physicalTags 

313 ) 

314 ) 

315 buffer = content[:numElementsInBlock] 

316 content = content[numElementsInBlock:] 

317 for row in buffer: 

318 eid_and_gids = map(int, row.split()) 

319 eid, *gids = eid_and_gids 

320 element = Element( 

321 eid=eid, 

322 gids=tuple(gids), 

323 entity_dim=entityDim, 

324 entity_tag=entityTag, 

325 element_type=elementType, 

326 physical_tags=physicalTags, 

327 physical_groups=physicalNames, 

328 ) 

329 self.elements[eid] = element 

330 

331 if __name__ == "__main__": 

332 import doctest 

333 

334 doctest.testmod(optionflags=doctest.ELLIPSIS | doctest.NORMALIZE_WHITESPACE)