Several questions about ALE.move

1

Dear all,

I am working on viscoelasticity equation with point forces, which are achieved by PointSource function and I need to solve the elasticity equation at each timestep. Thus, after each timestep, I used ALE.move(mesh,u_mechan) where u_mechan is the solution of the elasticity equation. Here are my questions about ALE.move.

  1. Do I need to update anything like subdomains, boundaries, function space etc I defined from the original mesh? Or these are all updated automatically and I don’t need to write down any other codes?

  2. After moving the mesh I also plot out the mesh, and sometimes when I apply PointSource function or evaluate a function at some point, I always got the error saying “The point is out of the domain”. However, according to the plot of new mesh, the point is definitely inside the domain and additionally I always use the point in the central region of the computational domain. Any idea or explanations of this scenario?

  3. If Q2 is due to the source code, can I use the change of variable to rewrite the weak form in the original mesh? Then how should I rewrite sigma(u) and epsilon(u) where

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

## sigma depends on the value of Young's modulus
def sigma_elas(u):
    return Es/(1+nu)*(epsilon(u)+tr(epsilon(u))*Identity(dim_mechan)*nu/(1-2*nu))

def sigma_visco(u):
    return Es/(2*(1+nu))*(epsilon(u)-epsilon(u_mechan_n))/dt+Es/(3*(1-2*nu))*(div(u)-div(u_mechan_n))*Identity(dim_mechan)/dt
​

and the weak form code is

F_mechan = inner(sigma_elas(u_mechan) + sigma_visco(u_mechan), nabla_grad(v_mechan)) * dx + kappa_mechan * inner(u_mechan, v_mechan) * ds
​

with some point sources?

Any help is appreciated. Thanks in advance.

Best,
Alice

2

Hi Alice, I will give a partial answer, which I hope it is useful. Keep in mind that I have never used PointSource, so I will miss many things surely.

  1. Do I need to update anything like subdomains, boundaries, function space etc I defined from the original mesh? Or these are all updated automatically and I don’t need to write down any other codes?

This depends on how your conditions are implemented. If they are calculated at every solve, then you might be missing points, i.e, if they are set with something like near(x[1], 1) and then your boundary moves away from that line, then the points won’t be marked. If otherwise your BC are set with markers, then they only depend on the label of the points, and wherever the points are located, the conditions will be set correctly. This is the approach I have used, and it works fine. More details here:

https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html

  1. After moving the mesh I also plot out the mesh, and sometimes when I apply PointSource function or evaluate a function at some point, I always got the error saying “The point is out of the domain”. However, according to the plot of new mesh, the point is definitely inside the domain and additionally I always use the point in the central region of the computational domain. Any idea or explanations of this scenario?

Maybe the pointsource gets fixed at instantiation and then it doesn’t read the mesh anymore, so it will not see the new coordinates. Honestly, no idea.

  1. If Q2 is due to the source code, can I use the change of variable to rewrite the weak form in the original mesh? Then how should I rewrite sigma(u) and epsilon(u) where…

Yes you can and it is what I would try, simply use the chain rule to change your differential operators: If X is the old configuration and x is the new one, you get

eps(u) = 0.5 (du/dX + du/dX^T) = 0.5 ( du/dx dx/dX + (du/dx dx/dX)^T),

this added to changing your measure with the Jacobian should work.

Best regards,
Nicolas

1 Like
3

Dear Nicolas,

Thanks a lot for your kind help and sorry for the late reply.

  1. I defined my boundary by doing something similar as near(x[1], 1) etc. but I then marked them so I guess it is going to work? Or since I only need the boundary of the computational domain, so I don’t clarify the boundary. Will that be a problem for moving mesh?

  2. I did some experiments and it seemed that after moving the mesh, Fenics does have a problem determining whether a point is inside the mesh or not. Thus, I don’t think I can figure it out by myself.

  3. Yes I managed to rewrite it wrt the domain integral but wrt the boundary integral, it is not hard to write based on single line segment but I can’t write a general form for the weak form.

Kind regard,
Alice

4

Hi Alice,

I have not tested these answeres by implementing a minimal working example, but they are rather based on experience. It would be easier to solve the problems if you post a minimal working example which represents each problem.

  1. If you have marked the boundaries of the mesh using a MeshFunction you do not need to update anything when moving the mesh. The boundary condition will be applied to the correct vertices.

  2. If you do not do this already, the answer to add mesh.bounding_box_tree().build(mesh) after you perform ALE.move(mesh, u_mechan).

Good luck!

Best,
Aslak

5

Hi Aslak,

Thanks for your kind help and the first one I think it works well now in my code but for the second one, there are still some problems, so I decided to solve the PDE in original mesh by change of variables. Anyway, thanks again for your kind reply.

Kind regards,
Alice