Coverage for pyxfem/core/mesh_modifications.py: 17%
166 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1''' module for some basic mesh modification of SMESH '''
2import numpy as np
3#from numpy.typing import NDArray
4from pyxfem.core.mesh_simplex_one_level import SMesh
5from pyxfem.core.graph import upward_adjacencies, UpwardAdjacencies, CompressedDataIndex, \
6 build_undirected_connection, coloring_undirected_connection, IndexType
10def duplicate_edges(mesh:SMesh, array_edges_id):
11 """ duplicate the edges in array_edges in the mesh and modify connectivity
12 in the mesh accordingly return the modifyed mesh and an array dupedge containing
13 the id of the new edges
14 so that dupedge[i] is the id of the reversed edge array_edges_id[i]
15 """
16 e2v = mesh.e2v()
17 t2e, t2eflip = mesh.t2e_flip()
18 e2t = mesh.e2t()
19 mesh.erase_edges_data()
20 nbfaces = e2t.index[array_edges_id+1] - e2t.index[array_edges_id]
21 if not np.all(nbfaces == 2):
22 print("error : Trying to duplicate boundary edges ...")
23 raise ValueError()
24 nedge0 = e2v.shape[0]
25 e2v_new = np.array(e2v[array_edges_id])[:, [1,0]]
26 dupedges = nedge0 + np.arange(array_edges_id.shape[0], dtype = IndexType)
27 e2v = np.vstack((e2v, e2v_new))
28 #print(e2tindex.dtype)
29 e2tris_id = np.column_stack( (e2t.data[e2t.index[array_edges_id]],
30 e2t.data[e2t.index[array_edges_id]+1 ]))
31 for (ie, e_id) in enumerate(array_edges_id) :
32 #print("Edge Eid TO dup", e_id)
33 it0 = e2tris_id[ie,0]
34 t0_2e = t2e[it0]
35 ie_in_t = 0 if t0_2e[0] == e_id else 1 if t0_2e[1] == e_id else 2
36 if t2eflip[it0][ie_in_t] :
37 t2e[it0, ie_in_t] = ie + nedge0
38 t2eflip[it0, ie_in_t] = False
39 else :
40 it1 = e2tris_id[ie,1]
41 t1_2e = t2e[it1]
42 ie_in_t = 0 if t1_2e[0] == e_id else 1 if t1_2e[1] == e_id else 2
43 t2e[it1, ie_in_t] = ie + nedge0
44 mesh._e2v, mesh._t2e, mesh._t2eflip = e2v, t2e, t2eflip
45 return mesh, dupedges
47def merge_edges(mesh:SMesh, duplicated_edges, tomerge) : #-> tuple[SMesh, any]|None:
48 ''' merging the edge '''
49 edge2edge = duplicated_edges[tomerge]
50 e2v = mesh.e2v()
51 t2e0, t2eflip0 = mesh.t2e_flip()
52 if edge2edge.shape[0] > 0 :
53 e2t = mesh.e2t()
54 t2e = t2e0.copy()
55 t2eflip = t2eflip0.copy()
56 e1 = edge2edge[:,1]
57 t2v = mesh.t2v()
58 for (e0, e1) in edge2edge:
59 for t in e2t[e1]:
60 for i in np.arange(3):
61 if t2e[t,i] == e1 :
62 t2eflip[t,i] = False if t2v[t, i] == e2v[e0,0] else True
63 nedgetomerge = edge2edge.shape[0]
64 nbedges = e2v.shape[0]
65 oldtonew_e = np.arange(nbedges, dtype = IndexType)
66 for i in np.arange(edge2edge.shape[0], dtype = IndexType):
67 e0, e1 = edge2edge[i]
68 oldtonew_e[e1] = e0
69 idshift_start= e1+1
70 idshift_end = nbedges if i == nedgetomerge-1 else edge2edge[i+1,1]
71 oldtonew_e[idshift_start:idshift_end] -= IndexType(i+1)
72 duplicated_edges = np.delete(duplicated_edges, tomerge, 0).reshape((-1,2))
73 duplicated_edges = oldtonew_e[duplicated_edges]
74 t2e = oldtonew_e[t2e]
75 nbedges = e2v.shape[0]
76 e2v = np.delete(e2v, edge2edge[:,1], axis = 0)
77 mesh.erase_edges_data()
78 mesh._e2v, mesh._t2e, mesh._t2eflip = e2v, t2e, t2eflip
79 return mesh, duplicated_edges
81def duplicate_nodes(mesh:SMesh, cut_eid):
82 ''' duplicate node and update the connectivity so that field is discontinuous across
83 the cut edge : cut_eid '''
84 if np.all(mesh.e2t().ndata_perentry()[cut_eid] != 2):
85 # with the curent algorithm, we can not work with boundary edge ....
86 raise ValueError
87 cut_e2v = mesh.e2v()[cut_eid]
88 nodes_to_check = np.unique(cut_e2v)
89 cut_v2e = upward_adjacencies(cut_e2v)
90 cut_v2e.data = cut_eid[cut_v2e.data]
91 def connect_patch(ndata, connectors:CompressedDataIndex):
92 patchconn = [[] for i in range(ndata)]
93 for i in range(len(connectors.index)-1):
94 t = connectors[i]
95 nt = len(t)
96 for i in range(nt):
97 patchconn[t[i]]+=list(t[:i])+ list(t[i+1:])
98 return build_undirected_connection(patchconn)
100 class Patch():
101 ''' representing patch information around vertex.
102 it store the vertex id (vid)
103 the ids of the triangles connected to vertex vid (tris)
104 the colors of each triangle (one color for each connected part of the triangles patch)
105 (2 triangle are connected if they share an edge)
106 '''
107 def __init__(self, vid=-1,
108 tris=np.array([], dtype = IndexType),
109 colors=np.array([], dtype = IndexType)):
110 self.vid = vid
111 self.tris = tris
112 self.colors = colors
113 self.ncolors = np.max(colors, initial=0)+1
114 def __str__(self):
115 return f"v{self.vid}, tris{self.tris}, colors{self.colors}"
116 def __repr__(self):
117 return str(self)
119 patches = []
120 for v in nodes_to_check:
121 v_tris_id = mesh.v2t()[v]
122 v_edges_id = np.setdiff1d(mesh.v2e()[v], cut_v2e[v])
123 e2t = mesh.e2t()
124 v_e2t_index = np.zeros(len(v_edges_id)+1, dtype = IndexType)
125 ss = e2t.index[v_edges_id]
126 es = e2t.index[v_edges_id+1]
127 np.cumsum(es - ss, out = v_e2t_index[1:])
128 v_e2t_data = np.concatenate( sum([[e2t.data[s:e]] for e, s in zip(es, ss)],[]) )
129 len_v_e2t_data = len(v_e2t_data)
130 v_e2t_data_tmp = np.concatenate([v_e2t_data, v_tris_id])
131 v_tris_id, v_e2t_data_renum = np.unique(v_e2t_data_tmp, return_inverse=True)
132 v_e2t_data_renum = v_e2t_data_renum[:len_v_e2t_data]
133 v_e2t = UpwardAdjacencies(v_e2t_index, v_e2t_data_renum)
134 v_t_connect = connect_patch(len(v_tris_id), v_e2t)
135 v_colors = coloring_undirected_connection(v_t_connect)
136 if np.max(v_colors, initial=0):
137 v_patch = Patch(v, v_tris_id, v_colors)
138 patches.append(v_patch)
140 vid = len(mesh.xy)
141 new_vxy= []
142 mesht2v = mesh.t2v()
143 for patch in patches :
144 patch:Patch= patch
145 vxy = mesh.xy[patch.vid]
146 ncopy = patch.ncolors-1
147 new_vxy += [vxy]*ncopy
148 patch_v_copies_id = [patch.vid] + [i for i in range(vid, vid+ncopy)]
149 for t, c in zip(patch.tris, patch.colors) :
150 #t2v = mesh.t2v()[t]
151 t2v = mesht2v[t]
152 iv = np.argwhere(t2v == patch.vid)
153 t2v[iv] = patch_v_copies_id[c]
154 #mesh._t2v[t] = t2v
155 mesht2v[t] = t2v
156 vid+=ncopy
157 #mesh.xy = np.vstack([mesh.xy]+ new_vxy)
158 meshxy = np.vstack([mesh.xy]+ new_vxy)
159 mesh.update(meshxy, mesht2v)
163class MeshHierarchicModification():
164 ''' A class to store and do basic numbering for hierarchic modification;
165 currently used by levelsetmeshcut
166 '''
167 def __init__(self, xy, tris):
168 self.dim = 2
169 self.xy = xy
170 self.tris = tris
171 self.coord_dtype = self.xy.dtype
172 self.deleted_tris = np.zeros(0, dtype = IndexType)
173 self.new_tris = np.zeros((0,3), dtype = IndexType)
174 self.xy_new_vertices = np.zeros((0,self.dim), dtype = self.coord_dtype)
175 self.new_tris_parent = np.zeros(0, dtype = IndexType)
176 self.duplicated_vertices = np.zeros(0, dtype = IndexType)
177 self.split_edge_dtype = np.dtype( [('v0', IndexType), ('v1', IndexType),
178 ('s', self.coord_dtype), ('vid', IndexType)])
179 self.split_edges = np.zeros(0, dtype = self.split_edge_dtype)
181 self.cut_edge_vertices = np.zeros((0,2), dtype = IndexType)
182 self.cut_vertices = np.zeros((0,1), dtype = IndexType)
184 def get_new_mesh(self):
185 ''' return coordinates of vertices and triangles '''
186 xy = np.vstack((self.xy, self.xy_new_vertices[len(self.duplicated_vertices):]))
187 nbv = len(self.xy)
188 nbvnew = nbv + len(self.xy_new_vertices)- len(self.duplicated_vertices)
189 v2cut_added = np.hstack((self.duplicated_vertices,
190 np.arange(nbv, nbvnew, dtype =IndexType)),
191 dtype = IndexType)
192 tris = np.vstack((np.delete(self.tris, self.deleted_tris, axis=0),
193 v2cut_added[self.new_tris]),
194 dtype = IndexType)
195 return xy, tris
197 def get_new_mesh_pushdata(self, verticesdata_oldmesh, verticesdata_modifyedmesh):
198 ''' return coordinates of vertices and triangles and data for the new mesh obtained
199 after modifications '''
200 xy, tris = self.get_new_mesh()
201 verticesdata = np.vstack((verticesdata_oldmesh[:, np.newaxis],
202 verticesdata_modifyedmesh[len(self.duplicated_vertices):, np.newaxis]))
203 return xy, tris, verticesdata.squeeze()
205 def parents2sons(self):
206 ''' return parent to son data '''
207 son2parent = self.new_tris_parent
208 sons = np.argsort(son2parent)
209 index_parent2sons = np.cumsum(np.bincount(son2parent[sons]+1, minlength=len(self.tris)+1))
210 return index_parent2sons, sons
212if __name__ == "__main__":
213 def test_dup_node():
214 ''' testing node duplication on a simple mesh '''
215 vertices = np.array([[0.,0.],[-1.,0.], [-0.5, 0.5], [0.5, 0.5], [0.5, -0.5],[-0.5, -0.5]])
216 tris = np.array([[0,1,2],[0,2,3], [0,3,4], [0,4,5], [0,5,1]])
217 mesh = SMesh(vertices, tris)
218 cut_eid, __ = mesh.find_edgeids([[0,1],[0,3]])
219 print('edge_cut', mesh.e2v()[cut_eid])
220 duplicate_nodes(mesh, cut_eid)
221 print(mesh.xy)
222 print(mesh.t2v)
224 test_dup_node()