''' module for some basic mesh modification of SMESH '''
import numpy as np
#from numpy.typing import NDArray
from pyxfem.core.mesh_simplex_one_level import SMesh
from pyxfem.core.graph import upward_adjacencies, UpwardAdjacencies, CompressedDataIndex, \
build_undirected_connection, coloring_undirected_connection, IndexType
[docs]
def duplicate_edges(mesh:SMesh, array_edges_id):
""" duplicate the edges in array_edges in the mesh and modify connectivity
in the mesh accordingly return the modifyed mesh and an array dupedge containing
the id of the new edges
so that dupedge[i] is the id of the reversed edge array_edges_id[i]
"""
e2v = mesh.e2v()
t2e, t2eflip = mesh.t2e_flip()
e2t = mesh.e2t()
mesh.erase_edges_data()
nbfaces = e2t.index[array_edges_id+1] - e2t.index[array_edges_id]
if not np.all(nbfaces == 2):
print("error : Trying to duplicate boundary edges ...")
raise ValueError()
nedge0 = e2v.shape[0]
e2v_new = np.array(e2v[array_edges_id])[:, [1,0]]
dupedges = nedge0 + np.arange(array_edges_id.shape[0], dtype = IndexType)
e2v = np.vstack((e2v, e2v_new))
#print(e2tindex.dtype)
e2tris_id = np.column_stack( (e2t.data[e2t.index[array_edges_id]],
e2t.data[e2t.index[array_edges_id]+1 ]))
for (ie, e_id) in enumerate(array_edges_id) :
#print("Edge Eid TO dup", e_id)
it0 = e2tris_id[ie,0]
t0_2e = t2e[it0]
ie_in_t = 0 if t0_2e[0] == e_id else 1 if t0_2e[1] == e_id else 2
if t2eflip[it0][ie_in_t] :
t2e[it0, ie_in_t] = ie + nedge0
t2eflip[it0, ie_in_t] = False
else :
it1 = e2tris_id[ie,1]
t1_2e = t2e[it1]
ie_in_t = 0 if t1_2e[0] == e_id else 1 if t1_2e[1] == e_id else 2
t2e[it1, ie_in_t] = ie + nedge0
mesh._e2v, mesh._t2e, mesh._t2eflip = e2v, t2e, t2eflip
return mesh, dupedges
[docs]
def merge_edges(mesh:SMesh, duplicated_edges, tomerge) : #-> tuple[SMesh, any]|None:
''' merging the edge '''
edge2edge = duplicated_edges[tomerge]
e2v = mesh.e2v()
t2e0, t2eflip0 = mesh.t2e_flip()
if edge2edge.shape[0] > 0 :
e2t = mesh.e2t()
t2e = t2e0.copy()
t2eflip = t2eflip0.copy()
e1 = edge2edge[:,1]
t2v = mesh.t2v()
for (e0, e1) in edge2edge:
for t in e2t[e1]:
for i in np.arange(3):
if t2e[t,i] == e1 :
t2eflip[t,i] = False if t2v[t, i] == e2v[e0,0] else True
nedgetomerge = edge2edge.shape[0]
nbedges = e2v.shape[0]
oldtonew_e = np.arange(nbedges, dtype = IndexType)
for i in np.arange(edge2edge.shape[0], dtype = IndexType):
e0, e1 = edge2edge[i]
oldtonew_e[e1] = e0
idshift_start= e1+1
idshift_end = nbedges if i == nedgetomerge-1 else edge2edge[i+1,1]
oldtonew_e[idshift_start:idshift_end] -= IndexType(i+1)
duplicated_edges = np.delete(duplicated_edges, tomerge, 0).reshape((-1,2))
duplicated_edges = oldtonew_e[duplicated_edges]
t2e = oldtonew_e[t2e]
nbedges = e2v.shape[0]
e2v = np.delete(e2v, edge2edge[:,1], axis = 0)
mesh.erase_edges_data()
mesh._e2v, mesh._t2e, mesh._t2eflip = e2v, t2e, t2eflip
return mesh, duplicated_edges
[docs]
def duplicate_nodes(mesh:SMesh, cut_eid):
''' duplicate node and update the connectivity so that field is discontinuous across
the cut edge : cut_eid '''
if np.all(mesh.e2t().ndata_perentry()[cut_eid] != 2):
# with the curent algorithm, we can not work with boundary edge ....
raise ValueError
cut_e2v = mesh.e2v()[cut_eid]
nodes_to_check = np.unique(cut_e2v)
cut_v2e = upward_adjacencies(cut_e2v)
cut_v2e.data = cut_eid[cut_v2e.data]
def connect_patch(ndata, connectors:CompressedDataIndex):
patchconn = [[] for i in range(ndata)]
for i in range(len(connectors.index)-1):
t = connectors[i]
nt = len(t)
for i in range(nt):
patchconn[t[i]]+=list(t[:i])+ list(t[i+1:])
return build_undirected_connection(patchconn)
class Patch():
''' representing patch information around vertex.
it store the vertex id (vid)
the ids of the triangles connected to vertex vid (tris)
the colors of each triangle (one color for each connected part of the triangles patch)
(2 triangle are connected if they share an edge)
'''
def __init__(self, vid=-1,
tris=np.array([], dtype = IndexType),
colors=np.array([], dtype = IndexType)):
self.vid = vid
self.tris = tris
self.colors = colors
self.ncolors = np.max(colors, initial=0)+1
def __str__(self):
return f"v{self.vid}, tris{self.tris}, colors{self.colors}"
def __repr__(self):
return str(self)
patches = []
for v in nodes_to_check:
v_tris_id = mesh.v2t()[v]
v_edges_id = np.setdiff1d(mesh.v2e()[v], cut_v2e[v])
e2t = mesh.e2t()
v_e2t_index = np.zeros(len(v_edges_id)+1, dtype = IndexType)
ss = e2t.index[v_edges_id]
es = e2t.index[v_edges_id+1]
np.cumsum(es - ss, out = v_e2t_index[1:])
v_e2t_data = np.concatenate( sum([[e2t.data[s:e]] for e, s in zip(es, ss)],[]) )
len_v_e2t_data = len(v_e2t_data)
v_e2t_data_tmp = np.concatenate([v_e2t_data, v_tris_id])
v_tris_id, v_e2t_data_renum = np.unique(v_e2t_data_tmp, return_inverse=True)
v_e2t_data_renum = v_e2t_data_renum[:len_v_e2t_data]
v_e2t = UpwardAdjacencies(v_e2t_index, v_e2t_data_renum)
v_t_connect = connect_patch(len(v_tris_id), v_e2t)
v_colors = coloring_undirected_connection(v_t_connect)
if np.max(v_colors, initial=0):
v_patch = Patch(v, v_tris_id, v_colors)
patches.append(v_patch)
vid = len(mesh.xy)
new_vxy= []
mesht2v = mesh.t2v()
for patch in patches :
patch:Patch= patch
vxy = mesh.xy[patch.vid]
ncopy = patch.ncolors-1
new_vxy += [vxy]*ncopy
patch_v_copies_id = [patch.vid] + [i for i in range(vid, vid+ncopy)]
for t, c in zip(patch.tris, patch.colors) :
#t2v = mesh.t2v()[t]
t2v = mesht2v[t]
iv = np.argwhere(t2v == patch.vid)
t2v[iv] = patch_v_copies_id[c]
#mesh._t2v[t] = t2v
mesht2v[t] = t2v
vid+=ncopy
#mesh.xy = np.vstack([mesh.xy]+ new_vxy)
meshxy = np.vstack([mesh.xy]+ new_vxy)
mesh.update(meshxy, mesht2v)
[docs]
class MeshHierarchicModification():
''' A class to store and do basic numbering for hierarchic modification;
currently used by levelsetmeshcut
'''
def __init__(self, xy, tris):
self.dim = 2
self.xy = xy
self.tris = tris
self.coord_dtype = self.xy.dtype
self.deleted_tris = np.zeros(0, dtype = IndexType)
self.new_tris = np.zeros((0,3), dtype = IndexType)
self.xy_new_vertices = np.zeros((0,self.dim), dtype = self.coord_dtype)
self.new_tris_parent = np.zeros(0, dtype = IndexType)
self.duplicated_vertices = np.zeros(0, dtype = IndexType)
self.split_edge_dtype = np.dtype( [('v0', IndexType), ('v1', IndexType),
('s', self.coord_dtype), ('vid', IndexType)])
self.split_edges = np.zeros(0, dtype = self.split_edge_dtype)
self.cut_edge_vertices = np.zeros((0,2), dtype = IndexType)
self.cut_vertices = np.zeros((0,1), dtype = IndexType)
[docs]
def get_new_mesh(self):
''' return coordinates of vertices and triangles '''
xy = np.vstack((self.xy, self.xy_new_vertices[len(self.duplicated_vertices):]))
nbv = len(self.xy)
nbvnew = nbv + len(self.xy_new_vertices)- len(self.duplicated_vertices)
v2cut_added = np.hstack((self.duplicated_vertices,
np.arange(nbv, nbvnew, dtype =IndexType)),
dtype = IndexType)
tris = np.vstack((np.delete(self.tris, self.deleted_tris, axis=0),
v2cut_added[self.new_tris]),
dtype = IndexType)
return xy, tris
[docs]
def get_new_mesh_pushdata(self, verticesdata_oldmesh, verticesdata_modifyedmesh):
''' return coordinates of vertices and triangles and data for the new mesh obtained
after modifications '''
xy, tris = self.get_new_mesh()
verticesdata = np.vstack((verticesdata_oldmesh[:, np.newaxis],
verticesdata_modifyedmesh[len(self.duplicated_vertices):, np.newaxis]))
return xy, tris, verticesdata.squeeze()
[docs]
def parents2sons(self):
''' return parent to son data '''
son2parent = self.new_tris_parent
sons = np.argsort(son2parent)
index_parent2sons = np.cumsum(np.bincount(son2parent[sons]+1, minlength=len(self.tris)+1))
return index_parent2sons, sons
if __name__ == "__main__":
def test_dup_node():
''' testing node duplication on a simple mesh '''
vertices = np.array([[0.,0.],[-1.,0.], [-0.5, 0.5], [0.5, 0.5], [0.5, -0.5],[-0.5, -0.5]])
tris = np.array([[0,1,2],[0,2,3], [0,3,4], [0,4,5], [0,5,1]])
mesh = SMesh(vertices, tris)
cut_eid, __ = mesh.find_edgeids([[0,1],[0,3]])
print('edge_cut', mesh.e2v()[cut_eid])
duplicate_nodes(mesh, cut_eid)
print(mesh.xy)
print(mesh.t2v)
test_dup_node()