Coverage for pyxfem/core/level_set_mesh_cut.py: 40%
285 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Created on Fri Dec 1 16:17:21 2023
6@author: nchevaug
7"""
8from functools import partial
9from typing import Any
10import numpy as np
11import gmsh
12from pyxfem.core.config import FloatArray, IndexArray
13from pyxfem.core.mesh_simplex_one_level import DMesh
14from pyxfem.core.mesh_modifications import MeshHierarchicModification as MeshMod
15from pyxfem.core.exec_tools import Profile
16import pyxfem.core.gmsh_mesh_generator as mg
17import pyxfem.postpro.gmsh_post_pro_interface as gp
18import pyxfem.core.gmsh_mesh_model_interface as meshInterface
20''' This module deals with cutting a mesh with the iso zero of a level-set '''
21class SimpleLevelSetFunction():
22 '''class to group basic level set function for simple geometrical object'''
23 @staticmethod
24 def disk(center, radius):
25 ''' return the signed distance fonction to a circle defined by center and radius
26 (negative inside)'''
27 def phi(xy):
28 xy = np.atleast_2d(xy)
29 x, y = xy[:, 0], xy[:, 1]
30 return np.sqrt((x-center[0])**2+(y-center[1])**2) - radius
31 return phi
33 @staticmethod
34 def halfplane(point, theta):
35 '''
36 signed distance to a ligne passing thrue point makig angle theta with x axis;
37 '''
38 def phi(xy):
39 xy = np.atleast_2d(xy)
40 x, y = xy[:, 0], xy[:, 1]
41 return (x-point[0])*np.cos(theta) + (y-point[1])*np.sin(theta)
42 return phi
44 @staticmethod
45 def complement(ls):
46 ''' return the complement level set (opposite sign)'''
47 def phi(xy):
48 return -ls(xy)
49 return phi
51 @staticmethod
52 def union(lss):
53 ''' given a list of lss of levelset, return the union '''
54 def phi(xy):
55 xy = np.atleast_2d(xy)
56 allls = np.stack([ls(xy) for ls in lss], 1)
57 return np.min(allls, 1)
58 return phi
60 @staticmethod
61 def intersection(lss):
62 ''' given a list of lss of levelset, return the intersection '''
63 def phi(xy):
64 xy = np.atleast_2d(xy)
65 allls = np.stack([ls(xy) for ls in lss], 1)
66 return np.max(allls, 1)
67 return phi
70def compute_face_side(tris_ls:FloatArray):
71 ''' return on which size is a tri looking at vertices values of the ls
73 given tris_ls an 2d array such as each line i give the value of the levelset at
74 the nodes of an element
75 return an array containing the side of each element
76 (-1, 0, 1, respectivily inside, crossing o out of the iso zero of the ls)'''
77 return np.where(np.all(tris_ls >= 0., 1), 1, np.where(np.all(tris_ls <= 0., 1), -1, 0))
79def faces2cut(tris_ls:FloatArray):
80 ''' return which face need to be cut
82 given tris_ls an 2d array such as each line i give the value of
83 the levelset at the nodes of an element
84 return the indexes of the elements that needs to be cut'''
85 #np.where(np.all(np.stack((np.any(tris_ls > 0, 1), np.any(tris_ls < 0, 1)), 1), 1))[0]
86 f2c = np.all(np.stack((np.any(tris_ls > 0, 1), np.any(tris_ls < 0, 1)), 1), 1)
87 f2c = np.asarray(f2c).nonzero()[0]
88 return f2c
91def fit_to_vertices(vls:FloatArray, trisv:IndexArray, absfittol:float, relfittol:float):
92 ''' fit value of the level set so that iso-zero close to a vertex is snapped to the said vertex.
93 input :
94 - vls is an array containing the values of the levelset at each node,
95 - trisv is the element connectity (trisv[i,j] return node j of element i)
96 - absfittol : absolute fit tolerance. ls is fit to 0. where abs(ls) < absfittol
97 - relfittol : relative fit tolerance.
98 value of a level-set at a node is set to zero if one of it's edges
99 would be cut by iso-zero at parameter s < relfittol (edges are parametrized from 0 to 1)
100 '''
101 vls[np.abs(vls) < absfittol] = 0.
102 ev = trisv[:, [2, 0, 0, 1, 1, 2]].reshape((-1, 2))
103 evls = vls[ev]
104 s_index = np.where(evls[:, 0]*evls[:, 1] < 0)[0]
105 iev0ls = evls[s_index, 0]
106 iev1ls = evls[s_index, 1]
107 vls[ev[s_index[-iev0ls/(iev1ls-iev0ls) < relfittol], 0]] = 0.
110def set_levelset(vertices: Any, ls_at_vertices:FloatArray, setls):
111 ''' an internal function '''
112 for v, ls in zip(vertices, ls_at_vertices):
113 setls(v, ls)
115def fit_to_vertices_old(vertices: Any, getls, setls, absfittol:float, relfittol:float):
116 """ this is an old slow version of fit to vertices wiht too much loop.
118 still used for crackIsocut and multiLevelSetIsoCut ...
119 it work directly on mesh vertices of the dynamic mesh structure ....
120 it will desapear soon
121 """
122 for v in vertices:
123 lsv = getls(v)
124 if np.abs(lsv) < absfittol:
125 setls(v, 0.)
126 for v in vertices:
127 ls0 = getls(v)
128 if ls0 == 0.:
129 continue
130 for e in v.edges:
131 ls1 = getls(e.v1) if e.v0 == v else getls(e.v0)
132 if ls1*ls0 <= 0.:
133 s = -ls0/(ls1-ls0)
134 if s < relfittol:
135 setls(v, 0.)
136 break
139def docut(e, getls):
140 ''' return True if edge must be cut '''
141 return getls(e.v0)*getls(e.v1) < 0
144def cutpos(e, getls):
145 ''' return the position of the cut on edge e'''
146 v0, v1 = e.v0, e.v1
147 ls0, ls1 = getls(v0), getls(v1)
148 s = -ls0/(ls1-ls0)
149 xy = v0.xy*(1-s) + v1.xy*s
150 return s, xy
153def build_interpolatecallback_default(setls):
154 ''' callback that st the value of the levelset to 0 at new vertices '''
155 def interpolatecallback_default(s, vn, e):# pylint: disable=W0613
156 setls(vn, 0.)
157 return interpolatecallback_default
160def facecallback_propagate_parent(f, f0, f1):
161 ''' Call back that set parent on the new faces f0, f1 to be the same as f'''
162 parent = f.id if not hasattr(f, "parent") else f.parent
163 e = DMesh.get_common_edge(f0, f1)
164 f0.parent, f1.parent, e.parent = [parent]*3
166def lscut(m, getls, setls, fit=True, absfittol=1.e-6, relfittol=1.e-2,
167 interpolatecallback=None, edgecallback=DMesh.edge_cut_call_back_pass,
168 facecallback=DMesh.face_cut_call_back_pass):
169 ''' cut the mesh m using a level set '''
170 if interpolatecallback is None:
171 interpolatecallback = build_interpolatecallback_default(setls)
172 if fit:
173 fit_to_vertices_old(m.get_vertices(), getls,
174 setls, absfittol, relfittol)
175 edges2cut = [(e, cutpos(e, getls))
176 for e in m.get_edges() if docut(e, getls)]
177 edge2cutv = np.zeros((len(edges2cut), 2), dtype=np.dtype('int64'))
178 cutv = np.zeros((len(edges2cut), 1), dtype=np.dtype('int64'))
180 for ie, (e, (s, xy)) in enumerate(edges2cut):
181 vn = m.split_edge(e, xy, edgecallback, facecallback)
182 interpolatecallback(s, vn, e)
183 edge2cutv[ie] = np.array([e.v0.id, e.v1.id])
184 cutv[ie] = vn.id
185 return edge2cutv, cutv
187def levelset_isocut(xy:FloatArray, tris2v:IndexArray, vls:FloatArray, returnparent=True):
188 ''' improved version of levelSetIsoCut : only the part of the mesh that really need
189 to be cut are loaded in a dynamic mesh. Much faster
190 return improved data structure : meshModification
191 '''
192 modifs = MeshMod(xy, tris2v)
193 # tris_ls[i,j] contains the values of the level set at node j of element i
194 tris_ls = vls[tris2v]
195 tris2cut = faces2cut(tris_ls) # list of elements to cut
196 modifs.deleted_tris = tris2cut
197 # trisZcutv[i,j] : vertices j of triangle to cut i
198 tris2cutv = tris2v[tris2cut]
199 # list of vertices connected to triangles to cut. also used to remap from new to old numbering
200 v2cut = np.unique(tris2cutv)
201 modifs.duplicated_vertices = v2cut
202 if len(v2cut) == 0:
203 neworder = np.zeros(0, dtype=np.uint64)
204 else:
205 neworder = np.zeros(int(np.max(v2cut)+1), dtype=np.uint64)
206 neworder[v2cut] = np.arange(len(v2cut))
207 def getls(v):
208 return v.ls
209 def setls(v, ls):
210 v.ls = ls
211 m = DMesh(xy[v2cut], neworder[tris2cutv])
212 set_levelset(m.get_vertices(), vls[v2cut], setls)
213 if returnparent:
214 cutedgevertices_neworder, cut_vertices = lscut(
215 m, getls, setls, fit=False, facecallback=facecallback_propagate_parent)
216 else:
217 cutedgevertices_neworder, cut_vertices = lscut(
218 m, getls, setls, fit=False)
219 xy_cut_, tris2v_cut_ = m.get_coord_and_connect()
220 modifs.new_tris = tris2v_cut_
221 modifs.xy_new_vertices = xy_cut_
222 modifs.cut_edge_vertices = v2cut[cutedgevertices_neworder]
223 modifs.cut_vertices = cut_vertices
224 vls_cut_ = np.array([v.ls for v in m.get_vertices()])
225 if returnparent:
226 parent_id_neworder = np.array([f.id if not hasattr(
227 f, "parent") else f.parent for f in m.get_faces()], dtype=np.int64)
228 modifs.new_tris_parent = tris2cut[parent_id_neworder]
229 return modifs, vls_cut_
230 return modifs, vls_cut_
233def multi_levelset_isocut(xy, tris, list_of_ls_at_vertices, absfittol=1.e-12, relfittol=1.e-6):
234 ''' cut mesh with multiple level sets '''
235 def getls(i):
236 def getlsi(v):
237 return v.ls[i]
238 return getlsi
240 def setls(i):
241 def setlsi(v, ls):
242 v.ls[i] = ls
243 return setlsi
244 m = DMesh(xy, tris)
245 nls = len(list_of_ls_at_vertices)
246 for v in m.get_vertices():
247 v.ls = np.zeros(nls)
248 for f in m.get_faces():
249 f.side = np.zeros(nls)
250 for i in range(nls):
251 set_levelset(m.get_vertices(), list_of_ls_at_vertices[i], setls(i))
252 # for i in range(nls):
253 # def interpolatecallback(s, vn, e):
254 # vn.ls = np.zeros(nls)
255 # vn.ls[i] = 0.
256 # for k in set(range(nls)) - {i}:
257 # vn.ls[k] = e.v0.ls[k]*(1-s) + e.v1.ls[k]*s
258 # lscut(m, getls(i), setls(i), True, absfittol, relfittol,
259 # interpolatecallback=interpolatecallback)
260 def interpolatecallback_(s, vn, e, i):
261 vn.ls = np.zeros(nls)
262 vn.ls[i] = 0.
263 for k in set(range(nls)) - {i}:
264 vn.ls[k] = e.v0.ls[k]*(1-s) + e.v1.ls[k]*s
265 for i in range(nls):
266 interpolatecallback = partial(interpolatecallback_, i=i)
267 lscut(m, getls(i), setls(i), True, absfittol, relfittol,
268 interpolatecallback=interpolatecallback)
269 xy_cut, tri_cut = m.get_coord_and_connect()
270 ls_cut = np.array([v.ls for v in m.get_vertices()])
271 side_cut = np.hstack([compute_face_side(
272 ls_cut[:, i][tri_cut]).reshape((-1, 1)) for i in range(nls)])
273 return xy_cut, ls_cut, tri_cut, side_cut
275def crack_iso_cut(xy, tris, lsn, lst, fit2vertices=True, absfittol=1.e-12, relfittol=1.e-6):
276 ''' cut the mesh with a crack defined by lsn and lst '''
277 def set_lsn(v, lsn):
278 v.lsn = lsn
280 def get_lsn(v):
281 return v.lsn
283 def set_lst(v, lst):
284 v.lst = lst
286 def get_lst(v):
287 return v.lst
289 m = DMesh(xy, tris)
290 set_levelset(m.get_vertices(), lsn, set_lsn)
291 set_levelset(m.get_vertices(), lst, set_lst)
293 def interpolatecallback_lsncut(s, vn, e):
294 vn.lsn = 0.
295 vn.lst = e.v0.lst*(1-s) + e.v1.lst*s
297 def interpolatecallback_lstcut(s, vn, e):
298 vn.lst = 0.
299 vn.lsn = e.v0.lsn*(1-s) + e.v1.lsn*s
300 parent = e if not hasattr(e, "parent") else e.parent
301 vn.parent = parent
303 def edgecallback(e, e0, e1):
304 parent = e if not hasattr(e, "parent") else e.parent
305 vn = DMesh.get_common_vertex(e0, e1)
306 vn.parent, e0.parent, e1.parent = [parent]*3
308 def facecallback(f, f0, f1):
309 parent = f.id if not hasattr(f, "parent") else f.parent
310 e = DMesh.get_common_edge(f0, f1)
311 f0.parent, f1.parent, e.parent = [parent]*3
312 lscut(m, get_lsn, set_lsn, fit2vertices, absfittol, relfittol,
313 interpolatecallback_lsncut, edgecallback, facecallback)
314 lscut(m, get_lst, set_lst, fit2vertices, absfittol, relfittol,
315 interpolatecallback_lstcut, edgecallback, facecallback)
316 xy_cut, tri_cut = m.get_coord_and_connect()
317 lsn_cut = np.array([v.lsn for v in m.get_vertices()])
318 lst_cut = np.array([v.lst for v in m.get_vertices()])
319 siden_cut = compute_face_side(lsn_cut[tri_cut])
320 sidet_cut = compute_face_side(lst_cut[tri_cut])
321 parent_id = np.array([-1 if not hasattr(f, "parent")
322 else f.parent for f in m.get_faces()])
323 return xy_cut, lsn_cut, lst_cut, tri_cut, parent_id, siden_cut, sidet_cut
326def test_suite(*, element_size=0.2, profile=False, popupgmsh = False):
327 ''' function to test the tools defined in this file '''
328 Profile.doprofile = profile
329 test_fit2v = True
330 test_ls_isocut = True
331 test_multi_ls_isocut = True
332 test_crack_isocut = True
333 gmsh.initialize()
334 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
335 modelname, _ = mg.l_shape(element_size)
336 xy_appro = meshInterface.get_vertices_coord(modelname)
337 tris_appro = meshInterface.get_triangles(modelname)[0]
338 tris_appro_xy = xy_appro[tris_appro]
339 print('nb elements =', len(tris_appro))
340 print('nb vertices =', len(xy_appro))
341 if test_fit2v:
342 # tests with large value of relfittol (relative (toedges size) fit to vertex tolerance)
343 absfittol = 1.e-6
344 relfittol = 0.3
345 # very small test to check implem. result should be ('fit', [1.,2.,0.,0.])
346 tris = np.array([[0, 1, 2], [2, 1, 3]])
347 vls = np.array([1., 2., 0., 0.])
348 fit_to_vertices(vls, tris, absfittol, relfittol)
349 print('fit', vls)
350 # Next is a test on the loaded mesh. First we plot the level-set before the fit then after?
351 ls = SimpleLevelSetFunction.disk([0., 0.], .5)
352 vls = ls(xy_appro)
353 gp.listPlotFieldTri(tris_appro_xy, vls[tris_appro], P0=False, viewname='ls_before_fit',
354 IntervalsType="Filled iso-Values", NbIso=2, Range=[-1., 1.])
355 f2v = tris_appro
356 Profile.enable()
357 fit_to_vertices(vls, f2v, absfittol, relfittol)
358 Profile.disable()
359 gp.listPlotFieldTri(tris_appro_xy, vls[tris_appro], P0=False, viewname='ls_after_fit',
360 IntervalsType="Filled iso-Values", NbIso=2, Range=[-1., 1.])
361 m = DMesh(xy_appro, tris_appro)
362 vls = ls(xy_appro)
364 def getls(v):
365 return v.ls
367 def setls(v, ls):
368 v.ls = ls
370 set_levelset(m.get_vertices(), vls, setls)
371 fit_to_vertices_old(m.get_vertices(), getls,
372 setls, absfittol, relfittol)
373 xy_cut, __ = m.get_coord_and_connect()
374 # ls_cut = np.array([getls(v) for v in m.getVertices()])
375 gp.listPlotFieldTri(tris_appro_xy, vls[tris_appro], P0=False, viewname='ls_after_fit_old',
376 IntervalsType="Filled iso-Values", NbIso=2, Range=[-1., 1.])
377 if test_ls_isocut:
378 ls = SimpleLevelSetFunction.disk([0., 0.], .5)
379 vls = ls(xy_appro)
380 fit_to_vertices(vls, tris_appro, 1.e-6, 5.e-3)
381 Profile.enable()
382 # 3470091 noeud (element size = 1.e-3) -> 1.13 second 2 cut.
383 modifs, vls_cut = levelset_isocut(
384 xy_appro, tris_appro, vls, returnparent=True)
385 Profile.disable()
386 xy, tris, vls = modifs.get_new_mesh_pushdata(vls, vls_cut)
387 gp.listPlotFieldTri(xy[tris], vls[tris], P0=False, viewname="isocut_ls",
388 IntervalsType="Filled iso-Values", NbIso=10)
389 if test_multi_ls_isocut:
390 ls1 = SimpleLevelSetFunction.disk([0., 0.], .5)
391 ls2 = SimpleLevelSetFunction.disk([0., 0.], .7)
392 ls3 = SimpleLevelSetFunction.halfplane([0.3, 0.], 0.)
393 lss_xy = [ls1(xy_appro), ls2(xy_appro), ls3(xy_appro)]
394 Profile.enable()
395 xy, mphi, tris, mside = multi_levelset_isocut(
396 xy_appro, tris_appro, lss_xy, absfittol=0., relfittol=0.02)
397 Profile.disable()
398 for i in range(len(lss_xy)):
399 gp.listPlotFieldTri(xy[tris], mphi[:, i][tris],
400 P0=False, viewname=f"multi_isocut_ls_{i}")
401 gp.listPlotFieldTri(xy[tris], mside[:, i],
402 P0=True, viewname=f"multi_isocut_side_{i}")
403 if test_crack_isocut:
404 Profile.enable()
405 lsn = SimpleLevelSetFunction.halfplane([0., 0.5], np.pi/2.)
406 lst = SimpleLevelSetFunction.intersection([SimpleLevelSetFunction.halfplane(
407 [.5, 0.], 0.), SimpleLevelSetFunction.halfplane([-.5, 0.], np.pi)])
408 xy_cut, lsn_cut, lst_cut, tri_cut, parent_id, siden_cut, sidet_cut = crack_iso_cut(
409 xy_appro, tris_appro, lsn(xy_appro), lst(xy_appro),
410 absfittol=1.e-12, relfittol=1.e-6)
411 Profile.disable()
412 elem_id = np.arange(0, len(tris_appro))
413 gp.listPlotFieldTri(xy_cut[tri_cut], lsn_cut[tri_cut], P0=False,
414 viewname="crackisocut_lsn", IntervalsType="Filled iso-Values")
415 gp.listPlotFieldTri(xy_cut[tri_cut], siden_cut, P0=True,
416 viewname="crackisocut_lsn_side", IntervalsType="Filled iso-Values")
417 gp.listPlotFieldTri(xy_cut[tri_cut], lst_cut[tri_cut], P0=False,
418 viewname="crackisocut_lst", IntervalsType="Filled iso-Values")
419 gp.listPlotFieldTri(xy_cut[tri_cut], sidet_cut, P0=True,
420 viewname="crackisocut_lst_side", IntervalsType="Filled iso-Values")
421 gp.listPlotFieldTri(tris_appro_xy, elem_id, P0=True,
422 viewname="face_id", IntervalsType="Numeric values")
423 gp.listPlotFieldTri(xy_cut[tri_cut], parent_id, P0=True,
424 viewname="parentface_id", IntervalsType="Numeric values")
425 if popupgmsh:
426 gmsh.fltk.run()
427 gmsh.finalize()
428 Profile.print_stats()
430if __name__ == '__main__':
431 test_suite(popupgmsh=True)