ALE.move class and meshes

Question

Dear all,

I am using the ALE.move class and some displacement function for moving my old mesh to a new one. Furthermore I want to reset the new mesh to my old mesh, if my new mesh doesn't fullfill some criteria.

ALE.move(mesh, displacement)

I already tried to solve my problems with the inverse direction, like

ALE.move(mesh, -1.0*displacement)

but this produces some rounding errors I want to avoid. Furthermore I tried saving the mesh to some temporary variable, like

mesh_temp = mesh
subdomains_temp = subdomains
boundaries_temp = boundaries
ALE.move(mesh, displacement)
mesh = mesh_temp
subdomains = subdomains_temp
boundaries = boundaries_temp
Vol0 = assemble(Constant(1.0)*dx(domain=mesh))

but it seems not to work the way I expect it to work. When I calculate the volume and plot the mesh, it still calculates in Vol0 the volume of the deformed mesh and still plots the deformed mesh. Is there a way to work with the old, not deformed mesh?

In another test it seems to work with saving the components to a temporary variable, thats why I think the problem might has to do with ALE.move for some reason.

Comments

Setting mesh_temp = mesh and then moving mesh using ALE.move will also move mesh_temp. Try to set mesh_temp = Mesh(mesh). For your rounding errors, could you provided a minimal working example? If you are projecting your displacement field, that might cause rounding errors. If so, you can try to interpolate in stead.

---

Using the latest version of Fenics right now, as it is 2017.01

Tested that, but it doesn't work either. I'm getting some error about

Error: Unable to create Dirichlet boundary condition
Reason: User MeshFunctions and FunctionSpace meshes are different

I'll try to provide a minimal working example asap, but I'm getting my displacement field from a solution of a PDE, however I still think moving the mesh forward and backwards should lead to the same volume and not some offset in the volume.

---

If you set mesh_temp = Mesh(mesh) then you have constructed a new mesh-object. Hence if you want to use mesh_temp in some computation you have to update boundaries and subdomains as well, using the new mesh object. For example if boundaries is a MeshFunction on the facets you can simply do

boundaries_temp = MeshFunction("size_t", mesh_temp, 2)
boundaries_temp.set_values(boundaries.array())

and likewise with subdomains.

Moving the mesh back and forth should yield the same result unless the displacement field change in between. What type of function space does your displacement belong to? If it does not belong to a linear Lagrange space then your rounding error could be explained by the projection from your displacement space to the linear Lagrange space.

---

Hi,

now this is working - atleast in my mini-example (with some minor adjustments, as my mesh is generated by gmsh and the boundaries/subdomain data is stored in some other file). My dispalcement function is in

V = VectorFunctionSpace(mesh, "Lagrange", 1)

As I said, the displacement is the solution of the linear elasticity PDE (LHS) with some shape derivative as the RHS. I'm not using

project(displacement, V)

or something.

I'd just need some way without relying on the first mesh, as there are some accepted deformations of the mesh and I don't want them to reset. Like

... # assembling a, L
solve(a == L, displacement, bcs)
... # scaling displacement
ALE.move(mesh, displacement)
... # assuming new mesh accepted, new mesh is mesh1
... # assembling a, L
solve(a == L, displacement1, bcs)
*
ALE.move(mesh1, displacment1)
...  # new mesh not accepted - restoring mesh1 without ALE.move(mesh1, -1.0*displacement)
...

If new mesh accepted, then continue from new mesh, if mesh not accepted move mesh back/reload new mesh, so at the line with the star I'd need to save mesh1 with the according subdomains and boundaries, not the ones from the "old" mesh.

Answer 1

OK, I constructed an example below which hopefully contains some answers for you.

from dolfin import *

def solve_linear_elasticity(mesh, boundaries, d):
    c = Constant(d)

    V = VectorFunctionSpace(mesh, "Lagrange", 1)
    u = TrialFunction(V)
    v = TestFunction(V)

    E, nu = 10.0, 0.3
    mu = E/(2.0*(1.0 + nu))
    lmbda = E*nu/((1.0 + nu)*(1.0 -2.0*nu))
    sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(3)
    F = inner(sigma, grad(v))*dx 
    a, L = lhs(F), rhs(F)

    bcs = [DirichletBC(V, Constant((0.0, 0.0, 0.0)), boundaries, 1),
           DirichletBC(V.sub(0), c, boundaries, 2)]

    displacement = Function(V)
    solve(a==L, displacement, bcs)
    return displacement

def update_mesh(mesh, displacement, boundaries):

    new_mesh = Mesh(mesh)
    new_boundaries = MeshFunction("size_t", new_mesh, 2)
    new_boundaries.set_values(boundaries.array())
    ALE.move(new_mesh, displacement)
    return new_mesh, new_boundaries


# Original mesh
mesh = UnitCubeMesh(8, 8, 8)
subdomain1 = CompiledSubDomain("near(x[1], 0)")
subdomain2 = CompiledSubDomain("near(x[1], 1)")
boundaries = MeshFunction("size_t", mesh, 2)
boundaries.set_all(0)
subdomain1.mark(boundaries, 1)
subdomain2.mark(boundaries, 2)
plot(mesh, title = "Original mesh")

# First iteration (accepted)
displacement1 = solve_linear_elasticity(mesh, boundaries, 0.1)
mesh1, boundaries1 = update_mesh(mesh, displacement1, boundaries)
plot(mesh1, title = "First iteration (accepted)")

# Second iteration
displacement2 = solve_linear_elasticity(mesh1, boundaries1, 0.5)
mesh2, boundaries2 = update_mesh(mesh1, displacement2, boundaries1)
plot(mesh2, title = "Second iteration (rejected)")

# Something went wrong move mesh back to previous
displacement3 = Function(displacement2.function_space())
displacement3.vector()[:] = -1.0*displacement2.vector()
mesh3, boundaries3 = update_mesh(mesh2, displacement3, boundaries2)
plot(mesh3, title = "Moved back from second iteration (same as first iteration)")

from numpy.linalg import norm as np_norm
print np_norm(mesh3.coordinates() - mesh1.coordinates())

interactive()

Comments

Hi,

thanks for your help - I think I should be able to do it.