"""
Created on Mon Nov 20 16:26:58 2023
@author: nchevaug
"""
from itertools import count
from typing import Literal
import numpy as np
import pyxfem.core.gmsh_mesh_model_interface as mi
from pyxfem.core.config import FloatArray
from pyxfem.core.graph import (
E2V,
T2E,
T2V,
UpwardAdjacencies,
create_edges_from_tris,
upward_adjacencies,
)
def _vertices2vertices(e2v, v2e_index, v2e):
'''
return for each vertex the vertices it is connected too with one of the edges in e2v
e2v is an array containing for each line the two vetices it connect.
'''
nv = len(v2e_index)-1
v0 = np.repeat(np.arange(nv), v2e_index[1:]-v2e_index[:-1])
tmp = e2v[v2e]
v2v = np.where(tmp[:,0] == v0, tmp[:,1], tmp[:,0])
return v2e_index, v2v
[docs]
class SMesh:
''' Class to store and interogate simplexe meshes of Triangles '''
def __init__(self, xy:FloatArray, tris:T2V):
self.xy = xy
self._t2v = tris
self.groupe2tris = {}
self.groupe2edges_v = {}
# The following are secondary attribute, that can all be computed from xy, and t2v ...
# but it can takes time . That why they are saved for fast retrievial each time
# parameter save is passed.
# Off course these need to be cleaned if xy or t2v is modified ...
self._e2v = None
self._t2e = None
self._t2eflip = None
self._e2t = None
self._v2t = None
self._v2e = None
self.vid = None
self._t2vxy = None
self.e2vxy = None
self.elen = None
self.beid = None
self.ieid = None
self.groupe2edges = {}
[docs]
def update(self, xy, t2v):
''' Modify the mesh by changing the node list and the triangle list
Old vertices id in groupe2edges_v are supposed to stay same
old triangle id in groupe2tris are supposed to stay same
'''
self.xy = xy
self._t2v= t2v
self.vid = None
self.erase_edges_data()
[docs]
def t2v(self) :
''' return the triangle to vertices connectivity '''
return self._t2v
[docs]
def e2v(self, eids = None, *, save = True):
''' return the edge 2 vertices table '''
_e2v, __, __ = self._e2v_t2e(save = save)
if eids is None:
return _e2v
return _e2v[eids]
[docs]
def t2e_flip(self, *, save = True):
''' return the tri to edge table '''
_t2e : np.ndarray[tuple[int, int], np.dtype[np.intp]]
__, _t2e, _t2eflip = self._e2v_t2e(save = save)
return _t2e, _t2eflip
[docs]
def t2e(self, *, save = True):
''' return the tri to edge table '''
__, _t2e, __ = self._e2v_t2e(save = save)
return _t2e
[docs]
def v2t(self, save = True)->UpwardAdjacencies:
''' return vertex to triangle adjacencies '''
if self._v2t is not None:
_v2t = self._v2t
else :
_v2t = upward_adjacencies(self.t2v())
if save :
self._v2t = _v2t
return _v2t
[docs]
def v2e(self, save = True)->UpwardAdjacencies:
''' return vertex to edge adjacencies '''
if self._v2e is not None:
return self._v2e
_v2e = upward_adjacencies(self.e2v(save=save))
if save :
self._v2e = _v2e
return _v2e
[docs]
def e2t(self, save = True)->UpwardAdjacencies:
''' return e to triangle upward adjacencies '''
if self._e2t is not None:
return self._e2t
__, _t2e, __ = self._e2v_t2e(save = save)
_e2t = upward_adjacencies(_t2e)
if save :
self._e2t = _e2t
return _e2t
def _e2v_t2e(self, save = True)->tuple[E2V, T2E, T2E] :
''' return the edge to vertex table and the triangle to edge table
(+ t2eflip table where t2eflip[i,j] that tell if edge j of triangle i is fliped '''
if self._e2v is not None and self._t2e is not None and self._t2eflip is not None:
e2v, t2e, t2eflip = (self._e2v, self._t2e, self._t2eflip)
else:
e2v, t2e, t2eflip = create_edges_from_tris(self._t2v)
if save :
self._e2v, self._t2e, self._t2eflip = (e2v, t2e, t2eflip)
return e2v, t2e, t2eflip
[docs]
def set_trigroupe_from_tri_id(self, groupename, tids, warn_existing = True, check = True):
''' set a group of tri in the mesh '''
if check and ((min(tids) < 0) or (max(tids) >= self._t2v.shape[0])):
print(min(tids), max(tids),self._t2v.shape[0])
raise ValueError
if warn_existing and (groupename in self.groupe2tris):
print('warning : triGroupe ', groupename, 'already exist in sMesh')
self.groupe2tris[groupename] = tids
[docs]
def set_edgegroupe_from_edge2vertex_id(self, groupename, e2v, warn_existing = True):
''' set a group of edges in the mesh '''
if warn_existing and (groupename in self.groupe2edges):
print('warning : edgeGroupe ', groupename, 'already exist in sMesh')
self.groupe2edges[groupename] = self.find_edgeids(e2v)
[docs]
def get_edgegroupe(self, groupename):
''' For groupe name, find and return the edges id
(The global id of edge, not the input order ...)
and if the "swap" status of each edge in the form of a pair of 2 arrays
'''
if not groupename in self.groupe2edges and groupename in self.groupe2edges_v:
e2v = self.groupe2edges_v[groupename]
self.groupe2edges[groupename] = self.find_edgeids(e2v)
return self.groupe2edges.get(groupename,
(np.array([], dtype = np.int64), np.array([], dtype = np.int64)))
[docs]
def get_edgegroupe_names(self):
''' return the list of created edge groupes names '''
return list(self.groupe2edges_v.keys())
[docs]
def get_trigroupe(self, groupname):
''' return an array with all the tri id contained in the group '''
return self.groupe2tris.get(groupname,
(np.array([], dtype = np.int64), np.array([], dtype = np.int64)))
[docs]
def get_trigroupe_names(self):
''' return the list of created tri groupes names '''
return list(self.groupe2tris.keys())
[docs]
def ntris(self):
''' return the number of triangles in the mesh '''
return self._t2v.shape[0]
[docs]
def nv(self, save = True):
''' return the number of vertices in the mesh '''
return self.used_v(save = save).shape[0]
[docs]
def geo_dim(self):
''' return the geometric dimension of the mesh. For Now it just check the number of columns
of the vertices array self.xy, which can only be 2 for 2 dim mesh '''
if self.xy.shape[1] != 2 :
raise ValueError("Error : code tested only for geomDim == 2 (triangle in a plane ...) ")
return 2
[docs]
def used_v(self, *, save = True):
''' return the vertices ids that are used by at least 1 triangle '''
if self.vid is not None:
return self.vid
vid = np.unique(self._t2v.flatten())
if save :
self.vid = vid
return vid
[docs]
def t2v_coord(self, fids = None, *, save = True) :
''' return a table containing for each triangle the coordinates of its 3 Nodes
ret[i, j, k] is coordinate k of node j in [0,1,2] of element i
'''
if self._t2vxy is not None:
t2vxy = self._t2vxy
else :
t2vxy = self.xy[self._t2v]
if save :
self._t2vxy = t2vxy
if fids is None:
return t2vxy
return t2vxy[fids]
[docs]
def t2v_cog(self, *, save = True):
''' return a table containing the center of gravity of each triangle '''
return np.sum(self.t2v_coord(save=save), axis = 1)/3.
[docs]
def e2v_coord(self, eids = None, *, save = True) :
''' return a table containing for each edge the coordinates of it 2 nodes
ret[i, j, k] is coordinate k of node j in [0,1] of edge i '''
if self.e2vxy is not None:
e2vxy = self.e2vxy
else :
e2vxy = self.xy[self.e2v()]
if save:
self.e2vxy = e2vxy
if eids is None :
return e2vxy
return e2vxy[eids]
[docs]
def ecog(self, eids = None):
''' return a table containing the center of gravity of each edge'''
e2xy = self.e2v_coord(eids)
cog = 0.5*np.sum(e2xy, axis = 1)
return cog
[docs]
def enormal(self, eids = None, flips = None):
''' return the normal of the edges '''
if self.geo_dim() != 2 :
raise ValueError
e2xy = self.e2v_coord(eids)
e2dir = e2xy[:,1,:] - e2xy[:,0,:]
e2dirnorm = np.linalg.norm(e2dir, axis =1)
normal = np.column_stack((e2dir[:,1]/e2dirnorm, -e2dir[:,0]/e2dirnorm))
if flips is None:
return normal
else:
return np.where( (flips == 0)[:, np.newaxis], normal, -1.*normal)
[docs]
def elength(self, eids = None, *, save = True):
''' return a table containing the lenght of each edge'''
if self.elen is not None:
elen= self.elen
else:
e2vxy = self.e2v_coord(eids = None, save = save)
elen = np.linalg.norm( e2vxy[:, 1,:] - e2vxy[:, 0, :] , axis = 1)
if save:
self.elen = elen
if eids is None:
return elen
return elen[eids]
[docs]
def vcoord(self):
''' return the table of vertices coordinates '''
return self.xy
[docs]
def nedges(self, *, save = True):
''' return the number of edge in the mesh '''
e2v, __, __ = self._e2v_t2e(save = save)
return e2v.shape[0]
[docs]
def bound_edges(self, *, save = True):
''' return a table containing the edge at the boundary of the mesh '''
if self.beid is not None :
return self.beid
e2t = self.e2t(save = save)
beid = np.where( (e2t.index[1:] - e2t.index[:-1]) == 1)[0]
if save :
self.beid = beid
return beid
[docs]
def intern_edges(self, *, save = True):
''' return a table containing the internal edge of the mesh:
edge that connect ot at least 2 trian gles '''
if self.ieid is not None :
return self.ieid
e2t = self.e2t(save)
ieid = np.where( (e2t.index[1:] - e2t.index[:-1]) == 2)[0]
if np.any ((e2t.index[1:] - e2t.index[:-1]) > 2) :
raise ValueError("expecting only 2 face per internal Edges")
if save :
self.ieid = ieid
return ieid
[docs]
def find_edgeids(self, e2v):
''' return edge ids corresponding to edge defined by there 2 nodes ids '''
e2v_sorted = np.sort(e2v, axis=1)
e2v_mesh = self.e2v()
indexes_v0_start_inmesh = np.searchsorted(e2v_mesh[:,0], e2v_sorted[:,0], side='left')
indexes_v0_end_inmesh = np.searchsorted(e2v_mesh[:,0], e2v_sorted[:,0], side='right')
eids = []
eflips = []
for i, v1 in enumerate(e2v_sorted[:,1]) :
s = indexes_v0_start_inmesh[i]
e = indexes_v0_end_inmesh[i]
eid = s+np.searchsorted(e2v_mesh[s:e,1], [v1], side='left')
if s<= eid < e :
eids.append(eid[0])
eflips.append(not np.all( e2v_mesh[eid] == e2v[i]) )
return np.array(eids), np.array(eflips, dtype =np.int64)
[docs]
def get_vertices_on_edges_groups(self, edgegroupenames):
"return all vertices id on the edge groupe"
edgesid = np.array([], dtype = np.uint64)
for gname in edgegroupenames:
gedges = self.get_edgegroupe(gname)[0].astype(np.uint64)
edgesid = np.append(edgesid, gedges)
edgesid = np.unique(edgesid)
return np.unique(self.e2v(eids=edgesid).flatten())
[docs]
def erase_edges_data(self):
''' clean up all date related to the edges '''
self.e2vxy = None
self.elen = None
self._e2v = None
self._t2e = None
self._t2eflip = None
self._e2t = None
self._e2v = None
self._v2e = None
self.beid = None
self.ieid = None
self.groupe2edges = {}
[docs]
def gmshmodel2smesh(model_name,
edge_group_names = None, tri_group_names = None, dim:Literal[2,3] = 2):
''' build a smesh from a gmsh model '''
if edge_group_names is None:
edge_group_names = []
if tri_group_names is None:
tri_group_names = []
xy = mi.get_vertices_coord(model_name, dim)
tris, gmshfacetags = mi.get_triangles(model_name)
mesh = SMesh(xy, tris)
groupe2phystag = mi.get_group_name2dim_phys_tag(model_name, edge_group_names )
for name, (gdim, phystag) in groupe2phystag.items() :
if gdim == 1:
e2v = mi.get_toplevel_elements_vertices_on_phys_tag(model_name, gdim, phystag)['edge']
mesh.groupe2edges_v[name] = e2v
#mesh.set_edgegroupe_from_edge2vertex_id(name, e2v)
groupe2phystag = mi.get_group_name2dim_phys_tag(model_name, tri_group_names )
for name, (gdim, phystag) in groupe2phystag.items() :
if dim == 2:
gmsh_faces_id = mi.get_toplevel_elements_on_phys_tag(model_name, gdim, phystag)
gmshtags2_smeshid = -np.ones(np.uint64(np.max(gmshfacetags)+1), dtype =np.int64)
gmshtags2_smeshid[gmshfacetags] = np.arange(tris.shape[0])
mesh.set_trigroupe_from_tri_id(name, gmshtags2_smeshid[gmsh_faces_id['tri']])
return mesh
[docs]
class DMesh:
''' A class to represent dynamic simplexes mesh in 2D '''
[docs]
class Vertex:
''' vertex object '''
def __init__(self, vid, xy):
self.id, self.xy, self.edges = vid, xy, set()
def __str__(self):
return 'id: '+str(self.id)+', xy :'+str(self.xy)
[docs]
class Edge:
''' edge object '''
def __init__(self, eid, vertices):
self.id, (self.v0, self.v1), self.faces = eid, vertices, set()
def __str__(self):
return 'id: '+str(self.id)+', vertices :'+str(self.v0.id) + ' ' + str(self.v1.id)
[docs]
class Face:
''' face object '''
def __init__(self, fid, edges):
self.id, (self.e0, self.e1, self.e2) = fid, edges
def __str__(self):
return 'id: '+str(self.id)+', \
edges :'+str(self.e0.id) + ' ' + str(self.e1.id) + ' ' + str(self.e2.id)
[docs]
def get_vertices_id(self):
''' return the vertices id of the current face '''
v0, v1, v2 = self.get_vertices()
return [ v0.id, v1.id, v2.id ]
[docs]
def get_vertices(self):
''' return the vertices object of the current face '''
e0, e2 = self.e0, self.e2
if e2.v1 == e0.v0 :
return [e0.v0, e0.v1, e2.v0]
if e2.v1 == e0.v1 :
return [e0.v1, e0.v0, e2.v0]
if e2.v0 == e0.v0 :
return [e0.v0, e0.v1, e2.v1]
if e2.v0 == e0.v1 :
return [e0.v1, e0.v0, e2.v1]
raise ValueError
def __init__(self, xy, tris):
e2v, t2e, __ = create_edges_from_tris(tris)
v2e = upward_adjacencies(e2v)
e2t = upward_adjacencies(t2e)
_vertices = np.fromiter( (DMesh.Vertex(vid, vxy)
for vid, vxy in enumerate(xy)),
count = len(xy) , dtype = object)
_edges = np.fromiter( (DMesh.Edge(eid, _vertices[vertices])
for eid, vertices in enumerate(e2v) ),
count = len(e2v), dtype = object )
_faces = np.fromiter( (DMesh.Face(fid, _edges[edges])
for fid, edges in enumerate (t2e)),
count = len(t2e), dtype = object )
for v, edges in zip(_vertices, (set( _edges[v2e.data[s:e]] )
for (s,e) in zip(v2e.index[:-1], v2e.index[1:]) )):
v.edges = edges
for e, faces in zip(_edges, (set( _faces[e2t.data[s:e]])
for (s,e) in zip(e2t.index[:-1], e2t.index[1:]) )):
e.faces = faces
self.vertices = dict(enumerate(_vertices))
self.edges = dict(enumerate(_edges))
self.faces = dict(enumerate(_faces))
self.vid = count(len(self.vertices))
self.eid = count(len(self.edges))
self.fid = count(len(self.faces))
[docs]
@staticmethod
def get_common_vertex(e0, e1):
''' return common vertex of eà and e0 raise if none'''
v0e0, v1e0 = e0.v0, e0.v1
v0e1, v1e1 = e1.v0, e1.v1
if v0e0 == v0e1 :
return v0e0
if v0e0 == v1e1 :
return v0e0
if v1e0 == v0e1 :
return v1e0
if v1e0 == v1e1 :
return v1e0
raise ValueError
[docs]
@staticmethod
def get_common_edge(f0, f1):
''' return common edge between face f0 and f1 or raise '''
e = {f0.e0, f0.e1, f0.e2}&{f1.e0, f1.e1, f1.e2}
if len(e) == 1:
return e.pop()
raise ValueError
[docs]
def get_coord_and_connect(self):
''' return node coordinate and triangle to vertices connectivity'''
xy = np.fromiter( (v.xy for v in self.vertices.values()),
count = len(self.vertices) , dtype = (float, 2 ))
tris = np.fromiter( (f.get_vertices_id() for f in self.faces.values()),
count = len(self.faces), dtype = (int, 3) )
return xy, tris
[docs]
def insert_vertex(self, xyz):
''' add vertex with coordinate xyz in the mesh, return the handle to the new vertex '''
vid = next(self.vid)
vn = DMesh.Vertex(vid, xyz)
self.vertices[vid] = vn
return vn
[docs]
def insert_edge(self, evertices):
''' add edge linking the 2 vertices in evertices and retuen the handle to the new edge'''
eid = next(self.eid)
en = DMesh.Edge(eid, evertices)
for v in evertices :
v.edges.add(en)
self.edges[eid] = en
return en
[docs]
def insert_face(self, fedges):
''' add face linking the 3 edges in fedges and retuen the handle to the new face'''
fid = next(self.fid)
fn = DMesh.Face(fid, fedges)
for e in fedges :
e.faces.add(fn)
self.faces[fid] = fn
return fn
def _delete_face(self, f):
for e in [f.e0, f.e1, f.e2] :
e.faces.remove(f)
del self.faces[f.id]
def _delete_edge(self, e):
if len(e.faces) > 0 :
raise ValueError
for v in [e.v0, e.v1] :
v.edges.remove(e)
del self.edges[e.id]
def _del_vertex(self, v):
if len(v.edges) > 0 :
raise ValueError
del self.vertices[v.id]
[docs]
def get_vertices(self):
''' return all vertices coordinates '''
return self.vertices.values()
[docs]
def get_edges(self):
''' return all edge to vertices '''
return self.edges.values()
[docs]
def get_faces(self):
''' return all edge 2 faces '''
return self.faces.values()
[docs]
@staticmethod
def edge_cut_call_back_pass(e, e0, e1):
''' default edge cut call back '''
[docs]
@staticmethod
def face_cut_call_back_pass(f, f0, f1):
''' default face cut call back '''
[docs]
def split_edge(self, e, xyz,
edge_cut_cb = edge_cut_call_back_pass,
face_cut_cb = face_cut_call_back_pass):
''' split edge e suppose to be part of the mesh, at position xyz '''
v0, v1 = e.v0, e.v1
vn = self.insert_vertex(xyz)
e0 = self.insert_edge([v0, vn])
e1 = self.insert_edge([v1, vn])
edge_cut_cb(e, e0, e1)
ftodel = [f for f in e.faces]
for f in e.faces:
fedges = [f.e0, f.e1, f.e2]
fie = 0 if f.e0==e else 1 if f.e1 ==e else 2
fe0, fe1, fe2 = fedges[fie], fedges[(fie+1)%3], fedges[(fie+2)%3]
fv0 = DMesh.get_common_vertex(fe2, fe0)
vc = DMesh.get_common_vertex(fe1, fe2)
fen = self.insert_edge([vn, vc])
if fv0 == v0:
f0e, f1e = [e0, fen, fe2], [e1, fe1, fen]
else :
f0e, f1e = [e1, fen, fe2], [e0, fe1, fen]
f0 = self.insert_face(f0e)
f1 = self.insert_face(f1e)
face_cut_cb(f, f0, f1)
for f in ftodel:
self._delete_face(f)
self._delete_edge(e)
return vn