How to compute & display normals to one deforming boundary?

anneker · November 1, 2022, 2:36pm

Hi everyone,

I am currently working on solving a dynamic deformation problem on a 2D mesh: RectangleMesh(Point(xmin, ymin), Point(xmax, ymax), 20, 20, “right/left”). During the deformation, top boundary deforms, and my problem depends on top boundary normals.

So, I would like to compute and display the normal vectors to deforming top boundary of the rectangle mesh iteratively, at each timestep.
I have tried the below function, but when I display iteratively normals, they do not seem to be updated and stay on the initial mesh configuration.

My question:

Is my function get_new_topboundary_normals() correct (computation and display)?
How should I make sure top boundary normals are correctly updated into the problem? Recomputing them within the loop or within the residual form?
Could you maybe provide me some examples?

Many thanks,
Anne

I got inspired by code in Iterate over all facets on mesh surface to obtain average values for each facet

V = VectorFunctionSpace(mesh, “CG”, 1)
ds = Measure(“ds”, domain=mesh, subdomain_data=boundaries)

def get_new_topboundary_normals(mesh, dofs_mesh, dofs_topboundary, V, ds):
u_trial = TrialFunction(V)
v_test = TestFunction(V)
normals_topboundary = Function(V) # normals at top boundary vertices
n = FacetNormal(mesh)
a = inner(u_trial, v_test) * ds(1) # ds(1): “top” boundary
l = inner(n, v_test) * ds(1)
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
solve(A, normals_topboundary.vector(), L)
vdf.plot( normals_topboundary,
mode=‘mesh arrows’,
style=‘matplotlib’,
axes=4,
azimuth=0,
interactive=True,
viewup=[‘0.’,‘0.’,‘1.’],
scale=0.05).clear()

return normals_topboundary

dokken · November 1, 2022, 3:55pm

Please make sure that your code is appropriately formatted using markdown, i.e.

```python
def f(x):
     return x
```

Please also make sure that the code can reproduce the error message. This means that all necessary imports and definitions has to be added to the code.

anneker · November 2, 2022, 12:48pm

Hi @dokken, many thanks for the answer, please find below a formatted and reproducible code.
Could you explain me why normals to top boundary (and consequently the variational problem) are not updated at each timestep? Many thanks for your help.
Anne

mu = 1.
lambd = 1.25

def get_new_topboundary_normals(u_trial, v_test, normals_topboundary, mesh, ds):
    n = FacetNormal(mesh)
    a = inner(u_trial, v_test) * ds(1) 
    l = inner(n, v_test) * ds(1)
    A = assemble(a, keep_diagonal=True)
    L = assemble(l)
    A.ident_zeros()
    solve(A, normals_topboundary.vector(), L)
    # Plot updated normals to top boundary:
    #plot(normals_topboundary, mode="glyphs", title="Normals to top boundary", color='r') 
    #plt.show() 
    return normals_topboundary


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T) 

def sigma(u):
    return lambd*div(u)*Identity(d) + 2*mu*epsilon(u)

# Define mesh and function space
mesh = RectangleMesh(Point(0., 0.), Point(1., 0.5), 20, 20, "right/left")
V = VectorFunctionSpace(mesh, "CG", 1)
dx = Measure("dx", domain=mesh)

# Define top boundary domain
class Top(SubDomain): 
    def inside(self, x, on_boundary):
        return near(x[1], 0.5) 

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)  
top = Top()
top.mark(boundaries, 1)
ds = Measure("ds", subdomain_data=boundaries)

# Define DC condition
tol = 1E-14

def left(x, on_boundary):
    return on_boundary and x[0] < tol

def right(x, on_boundary):
    return on_boundary and x[0] > 1. - tol

bcl = DirichletBC(V, Constant((0., 0.)), left)
bcr = DirichletBC(V, Constant((0., 0.)), right)
bcs = [bcl, bcr]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
u_solution = Function(V)
normals_topboundary = Function(V)
forcing_topboundary = Function(V)
d = u_solution.geometric_dimension() 

normals_topboundary = get_new_topboundary_normals(u, v, normals_topboundary, mesh, ds) # needs to be udpated at each solve() call
forcing_topboundary = normals_topboundary * (-0.3)

a = inner(sigma(u), epsilon(v)) * dx
L = dot(forcing_topboundary, v)* ds(1) # force applied on top boundary, depending on normal vector

# Iterative problem resolution 
num_steps = 1000
dt = 5./num_steps
t = 0.0

for n in range(num_steps):
    t += dt

    # plot updated normals to top boundary
    plot(normals_topboundary, mode="glyphs", title="Normals to top boundary", color='r') 
    plt.show()

    # compute displacement
    solve(a == L, u_solution, bcs)
    plot(u_solution, mode='displacement', title='Displacement')
    plt.show()

bleyerj · November 2, 2022, 2:25pm

Hi, if you are using a small strain formulation of elasticity, computing normals on the deformed configuration will not make much sense as these are very close to the initial normals since the displacement should be very small.
In a finite-strain setting, you should compute deformed normals using Nanson’s formula. For instance, if F=Id + \nabla u is the deformation gradient and J=\det(F), then JF^{-T}N gives you a vector parallel to the new normal in the deformed configuration for a reference normal vector N. You just need to normalize it at the end.

anneker · November 2, 2022, 4:01pm

Hi @bleyerj, thank you for the idea. That’s interesting.
About your suggestion, if I compute n = J*F⁻T*N/norm(J*F⁻T*N) with N = Function(V), do you know how to make n also be a Function(V) (since J and F are not, by definition) ?

Additionally, even if I place myself in small strain elasticity, I should be able to plot normals_topboundary on the deformed mesh top boundary, is not it? The plot I get is as if normals were still in the reference configuration, and always displayed with undeformed mesh. Maybe, I miss something to update normals and display. Could you bring me some explanation, please?

dokken · November 2, 2022, 4:48pm

I would suggest using Paraview for these visualizations. Using the built-in matplotlib functionality you are limiting yourself. You could also try using Vedo: vedo/examples/other/dolfin at master · marcomusy/vedo · GitHub

anneker · November 2, 2022, 5:17pm

Thank you for your answer dokken.
Actually, I also tried with vedo and the result is the same: normals do not seem to be updated.
I aslo displayed them with Paraview and nothing moves.
Could you have another idea to update normals?

Many thanks,
Anne

dokken · November 2, 2022, 5:36pm

In your code you are not moving the mesh between time steps, so I’m not sure why you expect the mesh to deform. You would have to use the ALE.move function for this

anneker · November 3, 2022, 1:27pm

Thank you dokken, that works with ALE.move, normals move with boundary!

However, I got these updated normals (yellow) to top boundary:

Why don’t I get normals to the surface (i.e. local normals, not vertical vectors) ? Is it due to the FacetNormal function?

Many thanks

dokken · November 3, 2022, 2:14pm

As you have not included the code that you are running now, it is a bit tricky to help you.
However, my guess is that you added ALE.move inside the for-loop, but that you do not recompute the normals (calling get_new_topboundary_normals) inside the loop.

anneker · November 3, 2022, 3:09pm

Many thanks @dokken, it worked for me!

Here is the code updated with your suggestions:

from fenics import *
from vedo.dolfin import plot 
import time

mu = 1.
lambd = 1.25

def get_new_topboundary_normals(mesh, ds):
    u_trial = TrialFunction(V)
    v_test = TestFunction(V)
    normals_topboundary = Function(V)
    n = FacetNormal(mesh)
    a = inner(u_trial, v_test) * ds(1) 
    l = inner(n, v_test) * ds(1)
    A = assemble(a, keep_diagonal=True)
    L = assemble(l)
    A.ident_zeros()
    solve(A, normals_topboundary.vector(), L)
    return normals_topboundary

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T) 

def sigma(u):
    return lambd*div(u)*Identity(d) + 2*mu*epsilon(u)

# Define mesh and function space
mesh = RectangleMesh(Point(0., 0.), Point(1., 0.5), 20, 20, "right/left")
V = VectorFunctionSpace(mesh, "CG", 1)
dx = Measure("dx", domain=mesh)

# Define top boundary domain
class Top(SubDomain): 
    def inside(self, x, on_boundary):
        return near(x[1], 0.5) 

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)  
top = Top()
top.mark(boundaries, 1)
ds = Measure("ds", subdomain_data=boundaries)

# Define DC conditions
tol = 1E-14

def left(x, on_boundary):
    return on_boundary and x[0] < tol

def right(x, on_boundary):
    return on_boundary and x[0] > 1. - tol

bcl = DirichletBC(V, Constant((0., 0.)), left)
bcr = DirichletBC(V, Constant((0., 0.)), right)
bcs = [bcl, bcr]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
u_solution = Function(V)
normals_topboundary = Function(V)
forcing_topboundary = Function(V)
d = u_solution.geometric_dimension() 

# Define variational problem
forcing_topboundary = normals_topboundary * (-0.3)
a = inner(sigma(u), epsilon(v)) * dx
L = dot(forcing_topboundary, v) * ds(1) # force applied on top boundary, depending on normal vector

# Iterative problem resolution 
num_steps = 1000
dt = 5./num_steps
t = 0.0

for n in range(num_steps):
    t += dt
    
    # update normals to top boundary
    normals_topboundary.assign( get_new_topboundary_normals(mesh, ds) )
    plot(normals_topboundary, mode="mesh arrows", style="matplotlib", axes=4, azimuth=0, interactive=False, viewup=["0.","0.","1."], scale=0.05).clear() 
    time.sleep(1.)

    # get displacement
    solve(a == L, u_solution, bcs)
    """
    plot(u_solution, mode='displace', style='paraview', axes=4, azimuth=0, interactive=False, viewup=['0.','0.','1.']).clear()
    time.sleep(1.) 
    """

    # move the mesh
    ALE.move(mesh, u_solution)

Topic		Replies	Views
Obtain velocity normal to boundary	3	2110	January 6, 2023
How to plot normal unit vector of faces in a 2D mesh?	7	3918	August 10, 2020
Iterate over all facets on mesh surface to obtain average values for each facet	16	1806	September 24, 2021
How to compute curvature of a boundary	5	2155	June 5, 2020
Normals on the boundary mesh not pointing in the same direction	0	481	July 18, 2019

How to compute & display normals to one deforming boundary?

Related Topics