Unwanted growing boundary

[SimonK]

Hi,

im running the following code, which is similar to this (Precision in hyperelastic model) except that I create a XDMF file as result to view in ParaView 5.9.1.
But my result is differnt, it looks like my boundary is growing in all directions.
In the example it starts with a 100 mm circle, my circle is in the end in X- at ~55 and at X+ ~132. In the example its -50 and ~122.

Anyone having an idea where this could come from? Any help would be great. Hope this is formatted in the right way.

Im running FEniCS 2019.2.0.dev0 on Ubuntu 20.04 in Windows 10.

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

set_log_level(20)


class AllBoundaries(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]<-10

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

    bc = [DirichletBC(V, Constant((0,0)), mfunc, 1)]

    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).smooth().write('tmpmesh.xml') #vedo REMESHING + smoothing
    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)

file_results = XDMFFile("remesh_test_1.xdmf")

for i in range(N):
    mfunc = MeshFunction('size_t', mesh, 1, mesh.domains())
    mfunc.set_all(0)
    allb = AllBoundaries()
    allb.mark(mfunc, 1)

    F = np.array([2, (i-N/2)/N])

    displacement = solve_problem(mesh, mfunc, F)
    new_mesh = update(mesh, displacement)
    
    file_results.write(displacement,i)

    if do_remesh:
        mesh = remesh(new_mesh)
    else:
        mesh = new_mesh

Thank you!

[marcomusy]

Hi - interesting observation… it looks like it’s the smoothing thing, which slightly moves the boundary points, and it becomes visible as a cumulative effect. This is a revised version of the same:

import numpy as np
from dolfin import *
import vedo

set_log_level(30)


class AllBoundaries(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]<-10

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

    bc = [DirichletBC(V, Constant((0,0)), mfunc, 1)]

    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):
    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.generate_mesh(invert=vz<0).write('tmpmesh.xml') #vedo REMESHING + smoothing
    return Mesh("tmpmesh.xml")

#################################################################################
N = 40         # number of iterations of stretching
do_remesh = 0  # grab the boundary and remesh the interior at each iteration

circle = vedo.Circle(r=50)
mesh = remesh(circle)
half_circle = circle.boundaries().cut_with_plane(origin=[-10,0,0], normal='-x').z(2)
half_circle.linewidth(5).c("yellow4")

plt = vedo.Plotter(N=N, size=(2250, 1300))

meshes = [mesh]
displacements = []
for i in range(N):
    mfunc = MeshFunction('size_t', mesh, 1, mesh.domains())
    mfunc.set_all(0)
    allb = AllBoundaries()
    allb.mark(mfunc, 1)

    F = np.array([2, (i-N/2)/N])
    displacement = solve_problem(mesh, mfunc, F)
    new_mesh = update(mesh, displacement)

    if do_remesh:
        mesh = remesh(new_mesh)
    else:
        mesh = new_mesh
    meshes.append(mesh)
    displacements.append(displacement)

    arrow = vedo.Arrow2D([0,0], F*20).z(1).c("red4")
    vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()]).c("k5").lw(1)
    plt.at(i).show(f"step{i}", half_circle, vmesh, arrow)

plt.interactive().close()

no reeshing:

remeshing:

[SimonK]

Thank yout so much!
If I call the smooth() afterwards it does not happen, seems like the combination is causing it.