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
« 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
6import numpy as np
8from nastranio.readers.gmsh.mixins import Modificator
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)
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)
40class GmsHParser(Modificator):
41 def __init__(self, filepath):
42 self.read(filepath)
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
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
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()
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)
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)
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)
110 def _make_gid_eid_type(self):
111 """build:
113 * self.eid2gids
114 * self.gid2eids
115 * self.type2eid
116 * self.eid2type
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
134 def _make_eid_gid_physicalGroups(self):
135 """build:
137 * self.eid2physical_groups
138 * self.physical_group2eids
139 * self.physical_group2gids
140 * self.gid2physical_groups
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)
160 # =========================================================================
161 # internal stuff
162 # =========================================================================
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)
181 def _process_MeshFormat(self):
182 """creates self.mesh_format"""
183 self._mesh_format = self._toc["MeshFormat"][0].split()[0]
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
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
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()}
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 = {}
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
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 )
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
331 if __name__ == "__main__":
332 import doctest
334 doctest.testmod(optionflags=doctest.ELLIPSIS | doctest.NORMALIZE_WHITESPACE)