Coverage for pyxfem/core/mesh_simplex_one_level.py: 75%
374 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1"""
2Created on Mon Nov 20 16:26:58 2023
4@author: nchevaug
5"""
7from itertools import count
8from typing import Literal
10import numpy as np
12import pyxfem.core.gmsh_mesh_model_interface as mi
13from pyxfem.core.config import FloatArray
14from pyxfem.core.graph import (
15 E2V,
16 T2E,
17 T2V,
18 UpwardAdjacencies,
19 create_edges_from_tris,
20 upward_adjacencies,
21)
24def _vertices2vertices(e2v, v2e_index, v2e):
25 '''
26 return for each vertex the vertices it is connected too with one of the edges in e2v
27 e2v is an array containing for each line the two vetices it connect.
28 '''
29 nv = len(v2e_index)-1
30 v0 = np.repeat(np.arange(nv), v2e_index[1:]-v2e_index[:-1])
31 tmp = e2v[v2e]
32 v2v = np.where(tmp[:,0] == v0, tmp[:,1], tmp[:,0])
33 return v2e_index, v2v
35class SMesh:
36 ''' Class to store and interogate simplexe meshes of Triangles '''
37 def __init__(self, xy:FloatArray, tris:T2V):
38 self.xy = xy
39 self._t2v = tris
40 self.groupe2tris = {}
41 self.groupe2edges_v = {}
42 # The following are secondary attribute, that can all be computed from xy, and t2v ...
43 # but it can takes time . That why they are saved for fast retrievial each time
44 # parameter save is passed.
45 # Off course these need to be cleaned if xy or t2v is modified ...
46 self._e2v = None
47 self._t2e = None
48 self._t2eflip = None
49 self._e2t = None
50 self._v2t = None
51 self._v2e = None
52 self.vid = None
53 self._t2vxy = None
54 self.e2vxy = None
55 self.elen = None
56 self.beid = None
57 self.ieid = None
58 self.groupe2edges = {}
60 def update(self, xy, t2v):
61 ''' Modify the mesh by changing the node list and the triangle list
62 Old vertices id in groupe2edges_v are supposed to stay same
63 old triangle id in groupe2tris are supposed to stay same
64 '''
65 self.xy = xy
66 self._t2v= t2v
67 self.vid = None
68 self.erase_edges_data()
70 def t2v(self) :
71 ''' return the triangle to vertices connectivity '''
72 return self._t2v
74 def e2v(self, eids = None, *, save = True):
75 ''' return the edge 2 vertices table '''
76 _e2v, __, __ = self._e2v_t2e(save = save)
77 if eids is None:
78 return _e2v
79 return _e2v[eids]
81 def t2e_flip(self, *, save = True):
82 ''' return the tri to edge table '''
83 _t2e : np.ndarray[tuple[int, int], np.dtype[np.intp]]
84 __, _t2e, _t2eflip = self._e2v_t2e(save = save)
85 return _t2e, _t2eflip
87 def t2e(self, *, save = True):
88 ''' return the tri to edge table '''
89 __, _t2e, __ = self._e2v_t2e(save = save)
90 return _t2e
92 def v2t(self, save = True)->UpwardAdjacencies:
93 ''' return vertex to triangle adjacencies '''
94 if self._v2t is not None:
95 _v2t = self._v2t
96 else :
97 _v2t = upward_adjacencies(self.t2v())
98 if save :
99 self._v2t = _v2t
100 return _v2t
102 def v2e(self, save = True)->UpwardAdjacencies:
103 ''' return vertex to edge adjacencies '''
104 if self._v2e is not None:
105 return self._v2e
106 _v2e = upward_adjacencies(self.e2v(save=save))
107 if save :
108 self._v2e = _v2e
109 return _v2e
111 def e2t(self, save = True)->UpwardAdjacencies:
112 ''' return e to triangle upward adjacencies '''
113 if self._e2t is not None:
114 return self._e2t
115 __, _t2e, __ = self._e2v_t2e(save = save)
116 _e2t = upward_adjacencies(_t2e)
117 if save :
118 self._e2t = _e2t
119 return _e2t
121 def _e2v_t2e(self, save = True)->tuple[E2V, T2E, T2E] :
122 ''' return the edge to vertex table and the triangle to edge table
123 (+ t2eflip table where t2eflip[i,j] that tell if edge j of triangle i is fliped '''
125 if self._e2v is not None and self._t2e is not None and self._t2eflip is not None:
126 e2v, t2e, t2eflip = (self._e2v, self._t2e, self._t2eflip)
127 else:
128 e2v, t2e, t2eflip = create_edges_from_tris(self._t2v)
129 if save :
130 self._e2v, self._t2e, self._t2eflip = (e2v, t2e, t2eflip)
131 return e2v, t2e, t2eflip
133 def set_trigroupe_from_tri_id(self, groupename, tids, warn_existing = True, check = True):
134 ''' set a group of tri in the mesh '''
135 if check and ((min(tids) < 0) or (max(tids) >= self._t2v.shape[0])):
136 print(min(tids), max(tids),self._t2v.shape[0])
137 raise ValueError
138 if warn_existing and (groupename in self.groupe2tris):
139 print('warning : triGroupe ', groupename, 'already exist in sMesh')
140 self.groupe2tris[groupename] = tids
142 def set_edgegroupe_from_edge2vertex_id(self, groupename, e2v, warn_existing = True):
143 ''' set a group of edges in the mesh '''
144 if warn_existing and (groupename in self.groupe2edges):
145 print('warning : edgeGroupe ', groupename, 'already exist in sMesh')
146 self.groupe2edges[groupename] = self.find_edgeids(e2v)
148 def get_edgegroupe(self, groupename):
149 ''' For groupe name, find and return the edges id
150 (The global id of edge, not the input order ...)
151 and if the "swap" status of each edge in the form of a pair of 2 arrays
152 '''
153 if not groupename in self.groupe2edges and groupename in self.groupe2edges_v:
154 e2v = self.groupe2edges_v[groupename]
155 self.groupe2edges[groupename] = self.find_edgeids(e2v)
156 return self.groupe2edges.get(groupename,
157 (np.array([], dtype = np.int64), np.array([], dtype = np.int64)))
159 def get_edgegroupe_names(self):
160 ''' return the list of created edge groupes names '''
161 return list(self.groupe2edges_v.keys())
163 def get_trigroupe(self, groupname):
164 ''' return an array with all the tri id contained in the group '''
165 return self.groupe2tris.get(groupname,
166 (np.array([], dtype = np.int64), np.array([], dtype = np.int64)))
168 def get_trigroupe_names(self):
169 ''' return the list of created tri groupes names '''
170 return list(self.groupe2tris.keys())
172 def ntris(self):
173 ''' return the number of triangles in the mesh '''
174 return self._t2v.shape[0]
176 def nv(self, save = True):
177 ''' return the number of vertices in the mesh '''
178 return self.used_v(save = save).shape[0]
180 def geo_dim(self):
181 ''' return the geometric dimension of the mesh. For Now it just check the number of columns
182 of the vertices array self.xy, which can only be 2 for 2 dim mesh '''
183 if self.xy.shape[1] != 2 :
184 raise ValueError("Error : code tested only for geomDim == 2 (triangle in a plane ...) ")
185 return 2
187 def used_v(self, *, save = True):
188 ''' return the vertices ids that are used by at least 1 triangle '''
189 if self.vid is not None:
190 return self.vid
191 vid = np.unique(self._t2v.flatten())
192 if save :
193 self.vid = vid
194 return vid
196 def t2v_coord(self, fids = None, *, save = True) :
197 ''' return a table containing for each triangle the coordinates of its 3 Nodes
198 ret[i, j, k] is coordinate k of node j in [0,1,2] of element i
199 '''
200 if self._t2vxy is not None:
201 t2vxy = self._t2vxy
202 else :
203 t2vxy = self.xy[self._t2v]
204 if save :
205 self._t2vxy = t2vxy
206 if fids is None:
207 return t2vxy
208 return t2vxy[fids]
210 def t2v_cog(self, *, save = True):
211 ''' return a table containing the center of gravity of each triangle '''
212 return np.sum(self.t2v_coord(save=save), axis = 1)/3.
214 def e2v_coord(self, eids = None, *, save = True) :
215 ''' return a table containing for each edge the coordinates of it 2 nodes
216 ret[i, j, k] is coordinate k of node j in [0,1] of edge i '''
217 if self.e2vxy is not None:
218 e2vxy = self.e2vxy
219 else :
220 e2vxy = self.xy[self.e2v()]
221 if save:
222 self.e2vxy = e2vxy
223 if eids is None :
224 return e2vxy
225 return e2vxy[eids]
227 def ecog(self, eids = None):
228 ''' return a table containing the center of gravity of each edge'''
229 e2xy = self.e2v_coord(eids)
230 cog = 0.5*np.sum(e2xy, axis = 1)
231 return cog
233 def enormal(self, eids = None, flips = None):
234 ''' return the normal of the edges '''
235 if self.geo_dim() != 2 :
236 raise ValueError
237 e2xy = self.e2v_coord(eids)
238 e2dir = e2xy[:,1,:] - e2xy[:,0,:]
239 e2dirnorm = np.linalg.norm(e2dir, axis =1)
240 normal = np.column_stack((e2dir[:,1]/e2dirnorm, -e2dir[:,0]/e2dirnorm))
241 if flips is None:
242 return normal
243 else:
244 return np.where( (flips == 0)[:, np.newaxis], normal, -1.*normal)
246 def elength(self, eids = None, *, save = True):
247 ''' return a table containing the lenght of each edge'''
248 if self.elen is not None:
249 elen= self.elen
250 else:
251 e2vxy = self.e2v_coord(eids = None, save = save)
252 elen = np.linalg.norm( e2vxy[:, 1,:] - e2vxy[:, 0, :] , axis = 1)
253 if save:
254 self.elen = elen
255 if eids is None:
256 return elen
257 return elen[eids]
259 def vcoord(self):
260 ''' return the table of vertices coordinates '''
261 return self.xy
263 def nedges(self, *, save = True):
264 ''' return the number of edge in the mesh '''
265 e2v, __, __ = self._e2v_t2e(save = save)
266 return e2v.shape[0]
268 def bound_edges(self, *, save = True):
269 ''' return a table containing the edge at the boundary of the mesh '''
270 if self.beid is not None :
271 return self.beid
272 e2t = self.e2t(save = save)
273 beid = np.where( (e2t.index[1:] - e2t.index[:-1]) == 1)[0]
274 if save :
275 self.beid = beid
276 return beid
278 def intern_edges(self, *, save = True):
279 ''' return a table containing the internal edge of the mesh:
280 edge that connect ot at least 2 trian gles '''
281 if self.ieid is not None :
282 return self.ieid
283 e2t = self.e2t(save)
284 ieid = np.where( (e2t.index[1:] - e2t.index[:-1]) == 2)[0]
285 if np.any ((e2t.index[1:] - e2t.index[:-1]) > 2) :
286 raise ValueError("expecting only 2 face per internal Edges")
287 if save :
288 self.ieid = ieid
289 return ieid
291 def find_edgeids(self, e2v):
292 ''' return edge ids corresponding to edge defined by there 2 nodes ids '''
293 e2v_sorted = np.sort(e2v, axis=1)
294 e2v_mesh = self.e2v()
295 indexes_v0_start_inmesh = np.searchsorted(e2v_mesh[:,0], e2v_sorted[:,0], side='left')
296 indexes_v0_end_inmesh = np.searchsorted(e2v_mesh[:,0], e2v_sorted[:,0], side='right')
297 eids = []
298 eflips = []
299 for i, v1 in enumerate(e2v_sorted[:,1]) :
300 s = indexes_v0_start_inmesh[i]
301 e = indexes_v0_end_inmesh[i]
302 eid = s+np.searchsorted(e2v_mesh[s:e,1], [v1], side='left')
303 if s<= eid < e :
304 eids.append(eid[0])
305 eflips.append(not np.all( e2v_mesh[eid] == e2v[i]) )
306 return np.array(eids), np.array(eflips, dtype =np.int64)
308 def get_vertices_on_edges_groups(self, edgegroupenames):
309 "return all vertices id on the edge groupe"
310 edgesid = np.array([], dtype = np.uint64)
311 for gname in edgegroupenames:
312 gedges = self.get_edgegroupe(gname)[0].astype(np.uint64)
313 edgesid = np.append(edgesid, gedges)
314 edgesid = np.unique(edgesid)
315 return np.unique(self.e2v(eids=edgesid).flatten())
317 def erase_edges_data(self):
318 ''' clean up all date related to the edges '''
319 self.e2vxy = None
320 self.elen = None
321 self._e2v = None
322 self._t2e = None
323 self._t2eflip = None
324 self._e2t = None
325 self._e2v = None
326 self._v2e = None
327 self.beid = None
328 self.ieid = None
329 self.groupe2edges = {}
331def gmshmodel2smesh(model_name,
332 edge_group_names = None, tri_group_names = None, dim:Literal[2,3] = 2):
333 ''' build a smesh from a gmsh model '''
334 if edge_group_names is None:
335 edge_group_names = []
336 if tri_group_names is None:
337 tri_group_names = []
338 xy = mi.get_vertices_coord(model_name, dim)
339 tris, gmshfacetags = mi.get_triangles(model_name)
340 mesh = SMesh(xy, tris)
341 groupe2phystag = mi.get_group_name2dim_phys_tag(model_name, edge_group_names )
342 for name, (gdim, phystag) in groupe2phystag.items() :
343 if gdim == 1:
344 e2v = mi.get_toplevel_elements_vertices_on_phys_tag(model_name, gdim, phystag)['edge']
345 mesh.groupe2edges_v[name] = e2v
346 #mesh.set_edgegroupe_from_edge2vertex_id(name, e2v)
348 groupe2phystag = mi.get_group_name2dim_phys_tag(model_name, tri_group_names )
349 for name, (gdim, phystag) in groupe2phystag.items() :
350 if dim == 2:
351 gmsh_faces_id = mi.get_toplevel_elements_on_phys_tag(model_name, gdim, phystag)
352 gmshtags2_smeshid = -np.ones(np.uint64(np.max(gmshfacetags)+1), dtype =np.int64)
353 gmshtags2_smeshid[gmshfacetags] = np.arange(tris.shape[0])
354 mesh.set_trigroupe_from_tri_id(name, gmshtags2_smeshid[gmsh_faces_id['tri']])
355 return mesh
357class DMesh:
358 ''' A class to represent dynamic simplexes mesh in 2D '''
359 class Vertex:
360 ''' vertex object '''
361 def __init__(self, vid, xy):
362 self.id, self.xy, self.edges = vid, xy, set()
363 def __str__(self):
364 return 'id: '+str(self.id)+', xy :'+str(self.xy)
365 class Edge:
366 ''' edge object '''
367 def __init__(self, eid, vertices):
368 self.id, (self.v0, self.v1), self.faces = eid, vertices, set()
369 def __str__(self):
370 return 'id: '+str(self.id)+', vertices :'+str(self.v0.id) + ' ' + str(self.v1.id)
371 class Face:
372 ''' face object '''
373 def __init__(self, fid, edges):
374 self.id, (self.e0, self.e1, self.e2) = fid, edges
375 def __str__(self):
376 return 'id: '+str(self.id)+', \
377 edges :'+str(self.e0.id) + ' ' + str(self.e1.id) + ' ' + str(self.e2.id)
378 def get_vertices_id(self):
379 ''' return the vertices id of the current face '''
380 v0, v1, v2 = self.get_vertices()
381 return [ v0.id, v1.id, v2.id ]
382 def get_vertices(self):
383 ''' return the vertices object of the current face '''
384 e0, e2 = self.e0, self.e2
385 if e2.v1 == e0.v0 :
386 return [e0.v0, e0.v1, e2.v0]
387 if e2.v1 == e0.v1 :
388 return [e0.v1, e0.v0, e2.v0]
389 if e2.v0 == e0.v0 :
390 return [e0.v0, e0.v1, e2.v1]
391 if e2.v0 == e0.v1 :
392 return [e0.v1, e0.v0, e2.v1]
393 raise ValueError
395 def __init__(self, xy, tris):
396 e2v, t2e, __ = create_edges_from_tris(tris)
397 v2e = upward_adjacencies(e2v)
398 e2t = upward_adjacencies(t2e)
399 _vertices = np.fromiter( (DMesh.Vertex(vid, vxy)
400 for vid, vxy in enumerate(xy)),
401 count = len(xy) , dtype = object)
402 _edges = np.fromiter( (DMesh.Edge(eid, _vertices[vertices])
403 for eid, vertices in enumerate(e2v) ),
404 count = len(e2v), dtype = object )
405 _faces = np.fromiter( (DMesh.Face(fid, _edges[edges])
406 for fid, edges in enumerate (t2e)),
407 count = len(t2e), dtype = object )
408 for v, edges in zip(_vertices, (set( _edges[v2e.data[s:e]] )
409 for (s,e) in zip(v2e.index[:-1], v2e.index[1:]) )):
410 v.edges = edges
411 for e, faces in zip(_edges, (set( _faces[e2t.data[s:e]])
412 for (s,e) in zip(e2t.index[:-1], e2t.index[1:]) )):
413 e.faces = faces
414 self.vertices = dict(enumerate(_vertices))
415 self.edges = dict(enumerate(_edges))
416 self.faces = dict(enumerate(_faces))
417 self.vid = count(len(self.vertices))
418 self.eid = count(len(self.edges))
419 self.fid = count(len(self.faces))
421 @staticmethod
422 def get_common_vertex(e0, e1):
423 ''' return common vertex of eà and e0 raise if none'''
424 v0e0, v1e0 = e0.v0, e0.v1
425 v0e1, v1e1 = e1.v0, e1.v1
426 if v0e0 == v0e1 :
427 return v0e0
428 if v0e0 == v1e1 :
429 return v0e0
430 if v1e0 == v0e1 :
431 return v1e0
432 if v1e0 == v1e1 :
433 return v1e0
434 raise ValueError
436 @staticmethod
437 def get_common_edge(f0, f1):
438 ''' return common edge between face f0 and f1 or raise '''
439 e = {f0.e0, f0.e1, f0.e2}&{f1.e0, f1.e1, f1.e2}
440 if len(e) == 1:
441 return e.pop()
442 raise ValueError
444 def get_coord_and_connect(self):
445 ''' return node coordinate and triangle to vertices connectivity'''
446 xy = np.fromiter( (v.xy for v in self.vertices.values()),
447 count = len(self.vertices) , dtype = (float, 2 ))
448 tris = np.fromiter( (f.get_vertices_id() for f in self.faces.values()),
449 count = len(self.faces), dtype = (int, 3) )
450 return xy, tris
451 def insert_vertex(self, xyz):
452 ''' add vertex with coordinate xyz in the mesh, return the handle to the new vertex '''
453 vid = next(self.vid)
454 vn = DMesh.Vertex(vid, xyz)
455 self.vertices[vid] = vn
456 return vn
457 def insert_edge(self, evertices):
458 ''' add edge linking the 2 vertices in evertices and retuen the handle to the new edge'''
459 eid = next(self.eid)
460 en = DMesh.Edge(eid, evertices)
461 for v in evertices :
462 v.edges.add(en)
463 self.edges[eid] = en
464 return en
465 def insert_face(self, fedges):
466 ''' add face linking the 3 edges in fedges and retuen the handle to the new face'''
467 fid = next(self.fid)
468 fn = DMesh.Face(fid, fedges)
469 for e in fedges :
470 e.faces.add(fn)
471 self.faces[fid] = fn
472 return fn
473 def _delete_face(self, f):
474 for e in [f.e0, f.e1, f.e2] :
475 e.faces.remove(f)
476 del self.faces[f.id]
477 def _delete_edge(self, e):
478 if len(e.faces) > 0 :
479 raise ValueError
480 for v in [e.v0, e.v1] :
481 v.edges.remove(e)
482 del self.edges[e.id]
483 def _del_vertex(self, v):
484 if len(v.edges) > 0 :
485 raise ValueError
486 del self.vertices[v.id]
487 def get_vertices(self):
488 ''' return all vertices coordinates '''
489 return self.vertices.values()
490 def get_edges(self):
491 ''' return all edge to vertices '''
492 return self.edges.values()
493 def get_faces(self):
494 ''' return all edge 2 faces '''
495 return self.faces.values()
497 @staticmethod
498 def edge_cut_call_back_pass(e, e0, e1):
499 ''' default edge cut call back '''
501 @staticmethod
502 def face_cut_call_back_pass(f, f0, f1):
503 ''' default face cut call back '''
505 def split_edge(self, e, xyz,
506 edge_cut_cb = edge_cut_call_back_pass,
507 face_cut_cb = face_cut_call_back_pass):
508 ''' split edge e suppose to be part of the mesh, at position xyz '''
509 v0, v1 = e.v0, e.v1
510 vn = self.insert_vertex(xyz)
511 e0 = self.insert_edge([v0, vn])
512 e1 = self.insert_edge([v1, vn])
513 edge_cut_cb(e, e0, e1)
514 ftodel = [f for f in e.faces]
515 for f in e.faces:
516 fedges = [f.e0, f.e1, f.e2]
517 fie = 0 if f.e0==e else 1 if f.e1 ==e else 2
518 fe0, fe1, fe2 = fedges[fie], fedges[(fie+1)%3], fedges[(fie+2)%3]
519 fv0 = DMesh.get_common_vertex(fe2, fe0)
520 vc = DMesh.get_common_vertex(fe1, fe2)
521 fen = self.insert_edge([vn, vc])
522 if fv0 == v0:
523 f0e, f1e = [e0, fen, fe2], [e1, fe1, fen]
524 else :
525 f0e, f1e = [e1, fen, fe2], [e0, fe1, fen]
526 f0 = self.insert_face(f0e)
527 f1 = self.insert_face(f1e)
528 face_cut_cb(f, f0, f1)
529 for f in ftodel:
530 self._delete_face(f)
531 self._delete_edge(e)
532 return vn