vtkDelaunay2D, Edge not recovered, polygon fill suspect

[SimonK]

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[1]<0
        
class Disp(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1]>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"][0][2] # 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!

[marcomusy]

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[1] < 0


class Disp(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1] > 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()[0]
    vz = vmesh.celldata["Normals"][0][2]  # 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

[SimonK]

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