#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 1 16:17:21 2023
@author: nchevaug
"""
from functools import partial
from typing import Any
import numpy as np
import gmsh
from pyxfem.core.config import FloatArray, IndexArray
from pyxfem.core.mesh_simplex_one_level import DMesh
from pyxfem.core.mesh_modifications import MeshHierarchicModification as MeshMod
from pyxfem.core.exec_tools import Profile
import pyxfem.core.gmsh_mesh_generator as mg
import pyxfem.postpro.gmsh_post_pro_interface as gp
import pyxfem.core.gmsh_mesh_model_interface as meshInterface
''' This module deals with cutting a mesh with the iso zero of a level-set '''
[docs]
class SimpleLevelSetFunction():
'''class to group basic level set function for simple geometrical object'''
[docs]
@staticmethod
def disk(center, radius):
''' return the signed distance fonction to a circle defined by center and radius
(negative inside)'''
def phi(xy):
xy = np.atleast_2d(xy)
x, y = xy[:, 0], xy[:, 1]
return np.sqrt((x-center[0])**2+(y-center[1])**2) - radius
return phi
[docs]
@staticmethod
def halfplane(point, theta):
'''
signed distance to a ligne passing thrue point makig angle theta with x axis;
'''
def phi(xy):
xy = np.atleast_2d(xy)
x, y = xy[:, 0], xy[:, 1]
return (x-point[0])*np.cos(theta) + (y-point[1])*np.sin(theta)
return phi
[docs]
@staticmethod
def complement(ls):
''' return the complement level set (opposite sign)'''
def phi(xy):
return -ls(xy)
return phi
[docs]
@staticmethod
def union(lss):
''' given a list of lss of levelset, return the union '''
def phi(xy):
xy = np.atleast_2d(xy)
allls = np.stack([ls(xy) for ls in lss], 1)
return np.min(allls, 1)
return phi
[docs]
@staticmethod
def intersection(lss):
''' given a list of lss of levelset, return the intersection '''
def phi(xy):
xy = np.atleast_2d(xy)
allls = np.stack([ls(xy) for ls in lss], 1)
return np.max(allls, 1)
return phi
[docs]
def compute_face_side(tris_ls:FloatArray):
''' return on which size is a tri looking at vertices values of the ls
given tris_ls an 2d array such as each line i give the value of the levelset at
the nodes of an element
return an array containing the side of each element
(-1, 0, 1, respectivily inside, crossing o out of the iso zero of the ls)'''
return np.where(np.all(tris_ls >= 0., 1), 1, np.where(np.all(tris_ls <= 0., 1), -1, 0))
[docs]
def faces2cut(tris_ls:FloatArray):
''' return which face need to be cut
given tris_ls an 2d array such as each line i give the value of
the levelset at the nodes of an element
return the indexes of the elements that needs to be cut'''
#np.where(np.all(np.stack((np.any(tris_ls > 0, 1), np.any(tris_ls < 0, 1)), 1), 1))[0]
f2c = np.all(np.stack((np.any(tris_ls > 0, 1), np.any(tris_ls < 0, 1)), 1), 1)
f2c = np.asarray(f2c).nonzero()[0]
return f2c
[docs]
def fit_to_vertices(vls:FloatArray, trisv:IndexArray, absfittol:float, relfittol:float):
''' fit value of the level set so that iso-zero close to a vertex is snapped to the said vertex.
input :
- vls is an array containing the values of the levelset at each node,
- trisv is the element connectity (trisv[i,j] return node j of element i)
- absfittol : absolute fit tolerance. ls is fit to 0. where abs(ls) < absfittol
- relfittol : relative fit tolerance.
value of a level-set at a node is set to zero if one of it's edges
would be cut by iso-zero at parameter s < relfittol (edges are parametrized from 0 to 1)
'''
vls[np.abs(vls) < absfittol] = 0.
ev = trisv[:, [2, 0, 0, 1, 1, 2]].reshape((-1, 2))
evls = vls[ev]
s_index = np.where(evls[:, 0]*evls[:, 1] < 0)[0]
iev0ls = evls[s_index, 0]
iev1ls = evls[s_index, 1]
vls[ev[s_index[-iev0ls/(iev1ls-iev0ls) < relfittol], 0]] = 0.
[docs]
def set_levelset(vertices: Any, ls_at_vertices:FloatArray, setls):
''' an internal function '''
for v, ls in zip(vertices, ls_at_vertices):
setls(v, ls)
[docs]
def fit_to_vertices_old(vertices: Any, getls, setls, absfittol:float, relfittol:float):
""" this is an old slow version of fit to vertices wiht too much loop.
still used for crackIsocut and multiLevelSetIsoCut ...
it work directly on mesh vertices of the dynamic mesh structure ....
it will desapear soon
"""
for v in vertices:
lsv = getls(v)
if np.abs(lsv) < absfittol:
setls(v, 0.)
for v in vertices:
ls0 = getls(v)
if ls0 == 0.:
continue
for e in v.edges:
ls1 = getls(e.v1) if e.v0 == v else getls(e.v0)
if ls1*ls0 <= 0.:
s = -ls0/(ls1-ls0)
if s < relfittol:
setls(v, 0.)
break
[docs]
def docut(e, getls):
''' return True if edge must be cut '''
return getls(e.v0)*getls(e.v1) < 0
[docs]
def cutpos(e, getls):
''' return the position of the cut on edge e'''
v0, v1 = e.v0, e.v1
ls0, ls1 = getls(v0), getls(v1)
s = -ls0/(ls1-ls0)
xy = v0.xy*(1-s) + v1.xy*s
return s, xy
[docs]
def build_interpolatecallback_default(setls):
''' callback that st the value of the levelset to 0 at new vertices '''
def interpolatecallback_default(s, vn, e):# pylint: disable=W0613
setls(vn, 0.)
return interpolatecallback_default
[docs]
def facecallback_propagate_parent(f, f0, f1):
''' Call back that set parent on the new faces f0, f1 to be the same as f'''
parent = f.id if not hasattr(f, "parent") else f.parent
e = DMesh.get_common_edge(f0, f1)
f0.parent, f1.parent, e.parent = [parent]*3
[docs]
def lscut(m, getls, setls, fit=True, absfittol=1.e-6, relfittol=1.e-2,
interpolatecallback=None, edgecallback=DMesh.edge_cut_call_back_pass,
facecallback=DMesh.face_cut_call_back_pass):
''' cut the mesh m using a level set '''
if interpolatecallback is None:
interpolatecallback = build_interpolatecallback_default(setls)
if fit:
fit_to_vertices_old(m.get_vertices(), getls,
setls, absfittol, relfittol)
edges2cut = [(e, cutpos(e, getls))
for e in m.get_edges() if docut(e, getls)]
edge2cutv = np.zeros((len(edges2cut), 2), dtype=np.dtype('int64'))
cutv = np.zeros((len(edges2cut), 1), dtype=np.dtype('int64'))
for ie, (e, (s, xy)) in enumerate(edges2cut):
vn = m.split_edge(e, xy, edgecallback, facecallback)
interpolatecallback(s, vn, e)
edge2cutv[ie] = np.array([e.v0.id, e.v1.id])
cutv[ie] = vn.id
return edge2cutv, cutv
[docs]
def levelset_isocut(xy:FloatArray, tris2v:IndexArray, vls:FloatArray, returnparent=True):
''' improved version of levelSetIsoCut : only the part of the mesh that really need
to be cut are loaded in a dynamic mesh. Much faster
return improved data structure : meshModification
'''
modifs = MeshMod(xy, tris2v)
# tris_ls[i,j] contains the values of the level set at node j of element i
tris_ls = vls[tris2v]
tris2cut = faces2cut(tris_ls) # list of elements to cut
modifs.deleted_tris = tris2cut
# trisZcutv[i,j] : vertices j of triangle to cut i
tris2cutv = tris2v[tris2cut]
# list of vertices connected to triangles to cut. also used to remap from new to old numbering
v2cut = np.unique(tris2cutv)
modifs.duplicated_vertices = v2cut
if len(v2cut) == 0:
neworder = np.zeros(0, dtype=np.uint64)
else:
neworder = np.zeros(int(np.max(v2cut)+1), dtype=np.uint64)
neworder[v2cut] = np.arange(len(v2cut))
def getls(v):
return v.ls
def setls(v, ls):
v.ls = ls
m = DMesh(xy[v2cut], neworder[tris2cutv])
set_levelset(m.get_vertices(), vls[v2cut], setls)
if returnparent:
cutedgevertices_neworder, cut_vertices = lscut(
m, getls, setls, fit=False, facecallback=facecallback_propagate_parent)
else:
cutedgevertices_neworder, cut_vertices = lscut(
m, getls, setls, fit=False)
xy_cut_, tris2v_cut_ = m.get_coord_and_connect()
modifs.new_tris = tris2v_cut_
modifs.xy_new_vertices = xy_cut_
modifs.cut_edge_vertices = v2cut[cutedgevertices_neworder]
modifs.cut_vertices = cut_vertices
vls_cut_ = np.array([v.ls for v in m.get_vertices()])
if returnparent:
parent_id_neworder = np.array([f.id if not hasattr(
f, "parent") else f.parent for f in m.get_faces()], dtype=np.int64)
modifs.new_tris_parent = tris2cut[parent_id_neworder]
return modifs, vls_cut_
return modifs, vls_cut_
[docs]
def multi_levelset_isocut(xy, tris, list_of_ls_at_vertices, absfittol=1.e-12, relfittol=1.e-6):
''' cut mesh with multiple level sets '''
def getls(i):
def getlsi(v):
return v.ls[i]
return getlsi
def setls(i):
def setlsi(v, ls):
v.ls[i] = ls
return setlsi
m = DMesh(xy, tris)
nls = len(list_of_ls_at_vertices)
for v in m.get_vertices():
v.ls = np.zeros(nls)
for f in m.get_faces():
f.side = np.zeros(nls)
for i in range(nls):
set_levelset(m.get_vertices(), list_of_ls_at_vertices[i], setls(i))
# for i in range(nls):
# def interpolatecallback(s, vn, e):
# vn.ls = np.zeros(nls)
# vn.ls[i] = 0.
# for k in set(range(nls)) - {i}:
# vn.ls[k] = e.v0.ls[k]*(1-s) + e.v1.ls[k]*s
# lscut(m, getls(i), setls(i), True, absfittol, relfittol,
# interpolatecallback=interpolatecallback)
def interpolatecallback_(s, vn, e, i):
vn.ls = np.zeros(nls)
vn.ls[i] = 0.
for k in set(range(nls)) - {i}:
vn.ls[k] = e.v0.ls[k]*(1-s) + e.v1.ls[k]*s
for i in range(nls):
interpolatecallback = partial(interpolatecallback_, i=i)
lscut(m, getls(i), setls(i), True, absfittol, relfittol,
interpolatecallback=interpolatecallback)
xy_cut, tri_cut = m.get_coord_and_connect()
ls_cut = np.array([v.ls for v in m.get_vertices()])
side_cut = np.hstack([compute_face_side(
ls_cut[:, i][tri_cut]).reshape((-1, 1)) for i in range(nls)])
return xy_cut, ls_cut, tri_cut, side_cut
[docs]
def crack_iso_cut(xy, tris, lsn, lst, fit2vertices=True, absfittol=1.e-12, relfittol=1.e-6):
''' cut the mesh with a crack defined by lsn and lst '''
def set_lsn(v, lsn):
v.lsn = lsn
def get_lsn(v):
return v.lsn
def set_lst(v, lst):
v.lst = lst
def get_lst(v):
return v.lst
m = DMesh(xy, tris)
set_levelset(m.get_vertices(), lsn, set_lsn)
set_levelset(m.get_vertices(), lst, set_lst)
def interpolatecallback_lsncut(s, vn, e):
vn.lsn = 0.
vn.lst = e.v0.lst*(1-s) + e.v1.lst*s
def interpolatecallback_lstcut(s, vn, e):
vn.lst = 0.
vn.lsn = e.v0.lsn*(1-s) + e.v1.lsn*s
parent = e if not hasattr(e, "parent") else e.parent
vn.parent = parent
def edgecallback(e, e0, e1):
parent = e if not hasattr(e, "parent") else e.parent
vn = DMesh.get_common_vertex(e0, e1)
vn.parent, e0.parent, e1.parent = [parent]*3
def facecallback(f, f0, f1):
parent = f.id if not hasattr(f, "parent") else f.parent
e = DMesh.get_common_edge(f0, f1)
f0.parent, f1.parent, e.parent = [parent]*3
lscut(m, get_lsn, set_lsn, fit2vertices, absfittol, relfittol,
interpolatecallback_lsncut, edgecallback, facecallback)
lscut(m, get_lst, set_lst, fit2vertices, absfittol, relfittol,
interpolatecallback_lstcut, edgecallback, facecallback)
xy_cut, tri_cut = m.get_coord_and_connect()
lsn_cut = np.array([v.lsn for v in m.get_vertices()])
lst_cut = np.array([v.lst for v in m.get_vertices()])
siden_cut = compute_face_side(lsn_cut[tri_cut])
sidet_cut = compute_face_side(lst_cut[tri_cut])
parent_id = np.array([-1 if not hasattr(f, "parent")
else f.parent for f in m.get_faces()])
return xy_cut, lsn_cut, lst_cut, tri_cut, parent_id, siden_cut, sidet_cut
[docs]
def test_suite(*, element_size=0.2, profile=False, popupgmsh = False):
''' function to test the tools defined in this file '''
Profile.doprofile = profile
test_fit2v = True
test_ls_isocut = True
test_multi_ls_isocut = True
test_crack_isocut = True
gmsh.initialize()
gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
modelname, _ = mg.l_shape(element_size)
xy_appro = meshInterface.get_vertices_coord(modelname)
tris_appro = meshInterface.get_triangles(modelname)[0]
tris_appro_xy = xy_appro[tris_appro]
print('nb elements =', len(tris_appro))
print('nb vertices =', len(xy_appro))
if test_fit2v:
# tests with large value of relfittol (relative (toedges size) fit to vertex tolerance)
absfittol = 1.e-6
relfittol = 0.3
# very small test to check implem. result should be ('fit', [1.,2.,0.,0.])
tris = np.array([[0, 1, 2], [2, 1, 3]])
vls = np.array([1., 2., 0., 0.])
fit_to_vertices(vls, tris, absfittol, relfittol)
print('fit', vls)
# Next is a test on the loaded mesh. First we plot the level-set before the fit then after?
ls = SimpleLevelSetFunction.disk([0., 0.], .5)
vls = ls(xy_appro)
gp.listPlotFieldTri(tris_appro_xy, vls[tris_appro], P0=False, viewname='ls_before_fit',
IntervalsType="Filled iso-Values", NbIso=2, Range=[-1., 1.])
f2v = tris_appro
Profile.enable()
fit_to_vertices(vls, f2v, absfittol, relfittol)
Profile.disable()
gp.listPlotFieldTri(tris_appro_xy, vls[tris_appro], P0=False, viewname='ls_after_fit',
IntervalsType="Filled iso-Values", NbIso=2, Range=[-1., 1.])
m = DMesh(xy_appro, tris_appro)
vls = ls(xy_appro)
def getls(v):
return v.ls
def setls(v, ls):
v.ls = ls
set_levelset(m.get_vertices(), vls, setls)
fit_to_vertices_old(m.get_vertices(), getls,
setls, absfittol, relfittol)
xy_cut, __ = m.get_coord_and_connect()
# ls_cut = np.array([getls(v) for v in m.getVertices()])
gp.listPlotFieldTri(tris_appro_xy, vls[tris_appro], P0=False, viewname='ls_after_fit_old',
IntervalsType="Filled iso-Values", NbIso=2, Range=[-1., 1.])
if test_ls_isocut:
ls = SimpleLevelSetFunction.disk([0., 0.], .5)
vls = ls(xy_appro)
fit_to_vertices(vls, tris_appro, 1.e-6, 5.e-3)
Profile.enable()
# 3470091 noeud (element size = 1.e-3) -> 1.13 second 2 cut.
modifs, vls_cut = levelset_isocut(
xy_appro, tris_appro, vls, returnparent=True)
Profile.disable()
xy, tris, vls = modifs.get_new_mesh_pushdata(vls, vls_cut)
gp.listPlotFieldTri(xy[tris], vls[tris], P0=False, viewname="isocut_ls",
IntervalsType="Filled iso-Values", NbIso=10)
if test_multi_ls_isocut:
ls1 = SimpleLevelSetFunction.disk([0., 0.], .5)
ls2 = SimpleLevelSetFunction.disk([0., 0.], .7)
ls3 = SimpleLevelSetFunction.halfplane([0.3, 0.], 0.)
lss_xy = [ls1(xy_appro), ls2(xy_appro), ls3(xy_appro)]
Profile.enable()
xy, mphi, tris, mside = multi_levelset_isocut(
xy_appro, tris_appro, lss_xy, absfittol=0., relfittol=0.02)
Profile.disable()
for i in range(len(lss_xy)):
gp.listPlotFieldTri(xy[tris], mphi[:, i][tris],
P0=False, viewname=f"multi_isocut_ls_{i}")
gp.listPlotFieldTri(xy[tris], mside[:, i],
P0=True, viewname=f"multi_isocut_side_{i}")
if test_crack_isocut:
Profile.enable()
lsn = SimpleLevelSetFunction.halfplane([0., 0.5], np.pi/2.)
lst = SimpleLevelSetFunction.intersection([SimpleLevelSetFunction.halfplane(
[.5, 0.], 0.), SimpleLevelSetFunction.halfplane([-.5, 0.], np.pi)])
xy_cut, lsn_cut, lst_cut, tri_cut, parent_id, siden_cut, sidet_cut = crack_iso_cut(
xy_appro, tris_appro, lsn(xy_appro), lst(xy_appro),
absfittol=1.e-12, relfittol=1.e-6)
Profile.disable()
elem_id = np.arange(0, len(tris_appro))
gp.listPlotFieldTri(xy_cut[tri_cut], lsn_cut[tri_cut], P0=False,
viewname="crackisocut_lsn", IntervalsType="Filled iso-Values")
gp.listPlotFieldTri(xy_cut[tri_cut], siden_cut, P0=True,
viewname="crackisocut_lsn_side", IntervalsType="Filled iso-Values")
gp.listPlotFieldTri(xy_cut[tri_cut], lst_cut[tri_cut], P0=False,
viewname="crackisocut_lst", IntervalsType="Filled iso-Values")
gp.listPlotFieldTri(xy_cut[tri_cut], sidet_cut, P0=True,
viewname="crackisocut_lst_side", IntervalsType="Filled iso-Values")
gp.listPlotFieldTri(tris_appro_xy, elem_id, P0=True,
viewname="face_id", IntervalsType="Numeric values")
gp.listPlotFieldTri(xy_cut[tri_cut], parent_id, P0=True,
viewname="parentface_id", IntervalsType="Numeric values")
if popupgmsh:
gmsh.fltk.run()
gmsh.finalize()
Profile.print_stats()
if __name__ == '__main__':
test_suite(popupgmsh=True)