Changing the shape in the surface normal directions

[Ekrem_Ekici]

Hi there,

Let’s say I wan to deform my mesh in the surface normal direction. For example, I have this geometry;

I want to update this geometry to this;

by manipualting the mesh with mesh.geometry.x[lateral_surface_indices]=....

I want to know how can I move the mesh points in the normal direction of lateral surface. How can I achieve that?

Here is the example code to generate the mesh;

import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio
L = 0.15
R = 0.2
gdim = 3
mesh_comm = MPI.COMM_WORLD
model_rank = 0

gmsh.initialize()
if mesh_comm.rank == model_rank:
    cylinder = gmsh.model.occ.add_cylinder(0,0,0,0,0,L,R)
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(2, [1], 1) # Lateral surface
    gmsh.model.addPhysicalGroup(2, [2], 2) # Top surface 
    gmsh.model.addPhysicalGroup(2, [3], 3) # Bottom surface

    gmsh.model.addPhysicalGroup(3, [1], 1) # Cylinder volume

    gmsh.option.setNumber("Mesh.MeshSizeMax", 0.03)
    gmsh.model.mesh.generate(gdim)

    gmsh.fltk.run()

mesh, subdomains, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

Thanks for your answers in advance.

[dokken]

@Ekrem_Ekici

I want to know how can I move the mesh points in the normal direction of lateral surface. How can I achieve that?

There are several things you need to consider here.
You want to move the surface (marked in yellow on your image) in a certain direction.
You want to fix the back surface (opposite of the yellow surface).
In between these two surfaces, you need some continuous (linear) extension to move all other nodes of the mesh.

In your example, you know the normal direction you want to move the nodes in, as they can be described as a linear function of one of the coordinate system directions.
Consider for instance this MWE:

import dolfinx
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl

L = 1
W = 1
NX = 25
domain = dolfinx.mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
                  [NX,6,6], cell_type=dolfinx.mesh.CellType.hexahedron)

x = domain.geometry.x

disp = 0.5
domain.geometry.x[:,0] = x[:,0] - disp * (x[:, 0]/L)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)

[Ekrem_Ekici]

Thanks for the answer @dokken.

I actually care about the green surface and I want to move it in the normal direction of it. Something like this;

[dokken]

You can Get an analytical expression for that movement, (i.e. r=(x[0]-c[0], x[1]-c[1], 0)) and apply this as a boundary condition for a vector Poisson equation to get a Smooth deformation field in the interior of your domain. Please note that you should prescribe Dirichlet conditions on the whole boundary of your mesh, otherwise you will get non-physical deformation at these boundaries.

[Ekrem_Ekici]

Could you provide MWE to explain what you mean by “vector Poisson equation”? My MWE above generates mesh with physical tags.

[dokken]

The simplest mesh deformation can be described as follows:

import dolfinx
from mpi4py import MPI
from typing import List
import numpy as np


def deform_mesh(V, bcs: List[dolfinx.fem.DirichletBCMetaClass]):
    mesh = V.mesh
    u = dolfinx.fem.Function(V)
    dolfinx.fem.petsc.set_bc(u.vector, bcs)
    deformation_array = u.x.array.reshape((-1, mesh.geometry.dim))
    mesh.geometry.x[:, :mesh.geometry.dim] += deformation_array


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD,  10, 10)

c_el = mesh.ufl_domain().ufl_coordinate_element()
V = dolfinx.fem.FunctionSpace(mesh, c_el)


def moving_boundaries(x):
    return np.isclose(x[1], 1) | np.isclose(x[0], 1)


mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
all_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
top_facets = dolfinx.mesh.locate_entities_boundary(
    mesh, mesh.topology.dim-1, moving_boundaries)

bc_fixed_facets = []
for facet in all_facets:
    if facet not in top_facets:
        bc_fixed_facets.append(facet)

fixed_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim-1, bc_fixed_facets)
c = dolfinx.fem.Constant(mesh, (0., 0.))
bc_fixed = dolfinx.fem.dirichletbc(c, fixed_dofs, V)


def radial(x):
    return (x[0], x[1])


u_radial = dolfinx.fem.Function(V)
u_radial.interpolate(radial)
top_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim-1, top_facets)
bc = dolfinx.fem.dirichletbc(u_radial, top_dofs)

bcs = [bc_fixed, bc]
deform_mesh(V, bcs)

with dolfinx.io.XDMFFile(mesh.comm, "mesh_deformed.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

as you can see here, only the boundary nodes are moved:

Something more advanced would be a vector Poisson equation, i.e. use:

def deform_mesh(V, bcs: List[dolfinx.fem.DirichletBCMetaClass]):
    mesh = V.mesh
    uh = dolfinx.fem.Function(V)
    dolfinx.fem.petsc.set_bc(uh.vector, bcs)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx
    L = ufl.inner(dolfinx.fem.Constant(mesh, (0., 0.)), v)*ufl.dx
    problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs, uh)
    problem.solve()
    deformation_array = uh.x.array.reshape((-1, mesh.geometry.dim))
    mesh.geometry.x[:, :mesh.geometry.dim] += deformation_array

which would give you

You could also use the equations of Linear elasticity instead of the vector Poisson equation

[sermarc]

Thanks @Ekrem_Ekici @dokken for the great discussion. I’m fighting with the same problem but with a small twist: I’m interested in sensitivity analysis so the perturbations are infinitesimal. My impression is that in this case, it is not necessary to smooth out the entire mesh, and for first-order meshes indeed just perturbing the boundary would be enough.

However, for higher-order meshes, my impression is that perturbing only the boundary nodes will break some ‘property’ of the mesh where mid-vertex nodes are in the middle of the vertex. It is also not entirely clear to me that Poisson smoothing will (completely) fix that, as the resulting map might be nonlinear. Is this line of thought correct?

[dokken]

The node never needs to be in the middle of the facet for a higher order mesh. There is no requirement for this.