# vtkDelaunay2D, Edge not recovered, polygon fill suspect

Hi,

im running the following code, which is based on this (Precision in hyperelastic model ) except that I increase the deformation to force my problem to occur. Its may not a total MWE but I tried to shorten my problem on this example.

The problem occurs after some iterations that the remeshing seems not to be able extracting the boundary right for remeshing. And this warning comes up:

``````vtkDelaunay2D.cxx:1419  WARN| vtkDelaunay2D (0x3286090): Edge not recovered, polygon fill suspect
``````

Anyone having an idea for a workaround or is this just a limit reached for this remeshing process?

``````from dolfin import *
from mshr import *
import numpy as np
import vedo

set_log_level(20)

class Fixed(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x<0

class Disp(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x>5

displ_max = 50
tend = 40

displ = Expression(('displ_max*time/tend','0.0'), displ_max = displ_max, time = 0, tend = tend, degree = 2)

def solve_problem(mesh, mfunc, force, displ, tend):
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
displacement = Function(V)

bc = [DirichletBC(V, Constant((0.,0.)), Fixed()),
DirichletBC(V, displ, Disp())]

F = Constant(force)
E = Constant(5000)
nu = Constant(0.3)
mu = E / (2.0*(1+nu))
lmbda = E * nu / ((1.0+nu) * (1-2*nu))
sigma = 2.0 * mu * sym(grad(u)) + lmbda * tr( sym(grad(u)) ) * Identity(2)
solve(inner(sigma, grad(v)) * dx == inner(F, v) * dx, displacement, bc)
displacement.set_allow_extrapolation(True)
return displacement

def update(mesh, displacement):
new_mesh = Mesh(mesh)
ALE.move(new_mesh, displacement)
return new_mesh

def remesh(mesh, res=50):
if isinstance(mesh, vedo.Mesh):
vmesh = mesh
else:
vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()])
bpts = vmesh.compute_normals(cells=True).boundaries().join(reset=1) #extract boundary
vz = vmesh.celldata["Normals"] # check z component of normals at first point
bpts.tomesh(invert=vz<0).write('tmpmesh.xml') #vedo REMESHING
return Mesh("tmpmesh.xml")

N = 40
res = 15
do_remesh = True
vedo.settings.use_parallel_projection = True

circle = Circle(Point(0, 0), 50)
mesh = generate_mesh(circle, res)

meshes = [mesh]
displacements = []

file_results = XDMFFile("remesh_test.xdmf")

for i in range(N):
displ.time = i

mfunc = MeshFunction('size_t', mesh, 1, mesh.domains())
mfunc.set_all(0)
allb = Fixed()
allb.mark(mfunc, 1)

F = np.array([0, 0])

displacement = solve_problem(mesh, mfunc, F, displ, tend)
new_mesh = update(mesh, displacement)

file_results.write(displacement,i)

if do_remesh:
mesh = remesh(new_mesh)
mesh.smooth(1)
else:
mesh = new_mesh
``````

Thanks for any idea!

I think the issue is that the line connecting the sharp boundaries is interpolated and might become invalid - infact already after a few iteration the mesh triangle look invalid:

you can reduce (but not remove completely) the problem by increasing the boundary and mesh resolution or provide a custom contour to generate the mesh

``````import numpy as np
from dolfin import *
import vedo
from vedo.dolfin import plot

set_log_level(20)

class Fixed(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x < 0

class Disp(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x > 5

def solve_problem(mesh, mfunc, force, displ, tend):
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
displacement = Function(V)
bc = [DirichletBC(V, Constant((0.0, 0.0)), Fixed()), DirichletBC(V, displ, Disp())]
F = Constant(force)
E = Constant(5000)
nu = Constant(0.3)
mu = E / (2.0 * (1 + nu))
lmbda = E * nu / ((1.0 + nu) * (1 - 2 * nu))
sigma = 2.0 * mu * sym(grad(u)) + lmbda * tr(sym(grad(u))) * Identity(2)
solve(inner(sigma, grad(v)) * dx == inner(F, v) * dx, displacement, bc)
displacement.set_allow_extrapolation(True)
return displacement

def update(mesh, displacement):
new_mesh = Mesh(mesh)
ALE.move(new_mesh, displacement)
return new_mesh

def remesh(mesh, res=50):
if isinstance(mesh, vedo.Mesh):
vmesh = mesh
else:
vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()])
bpts = vmesh.compute_normals(cells=True).boundaries().join_segments()
vz = vmesh.celldata["Normals"]  # check z component of normals at first point
m = bpts.generate_mesh(invert=vz < 0, mesh_resolution=res, line_resolution=res*3)
m.write("tmpmesh.xml")
return Mesh("tmpmesh.xml")

displ_max = 50
tend = 40
N = 40
res = 100
do_remesh = True
vedo.settings.use_parallel_projection = True

circle = vedo.Circle(r=50)
mesh = remesh(circle, res)

displ = Expression(
("displ_max*time/tend", "0.0"),
displ_max=displ_max, time=0, tend=tend, degree=2
)
file_results = XDMFFile("remesh_test.xdmf")

for i in range(N):
displ.time = i

mfunc = MeshFunction("size_t", mesh, 1, mesh.domains())
mfunc.set_all(0)
allb = Fixed()
allb.mark(mfunc, 1)

F = np.array([0, 0])

displacement = solve_problem(mesh, mfunc, F, displ, tend)
new_mesh = update(mesh, displacement)

file_results.write(displacement, i)
plot(displacement).clear()

if do_remesh:
mesh = remesh(new_mesh, res)
mesh.smooth(1)
else:
mesh = new_mesh
``````

2 Likes

Thank you for the fast response, looks like the right direction!