I am encountering a problem when attempting interpolations after changing a mesh twice. I am running 2019.2.0.dev0 on a Ubuntu 22.04 box.
I am examining different fluid-structure coupling conditions, and my within time step sequence consists of the steps
- Given two identical meshes: mesh_main and mesh_alt and old solution values
- Compute provisional fluid solution using mesh_main
- Estimate necessary amount of mesh deformation
- Move mesh_alt
- Interpolate all solution functions from mesh_main to moved mesh_alt
- Move mesh_main to agree with mesh_alt
- Copy solution functions from mesh_alt to mesh_main
This approach is failing for me, and I am attaching a minimal example that fails.
In my example, I am not solving anything, I am only interpolating a polynomial function of low enough degree that all interpolations should be exact. I discover that (1) moving mesh_alt seems to affect functions defined on mesh_main; and (2) the interpolations are correctly exactly computed on the first time through the loop but incorrect on the second.
#!/usr/bin/python3
from fenics import *
import matplotlib.pyplot as plt
def mover2(msh,xi):
V=VectorFunctionSpace(msh,"Lagrange",1)
u=Function(V)
expr=Expression( ("0.","xi*x[0]*(1.-x[0])") ,degree=1, xi=xi)
u.interpolate(expr)
ALE.move(msh,u)
# generate two identical meshes
mesh_main=UnitSquareMesh(5,5)
mesh_alt=UnitSquareMesh(5,5)
# define the vector spaces
Velt=VectorElement("Lagrange",mesh_main.ufl_cell(),2)
Pelt=FiniteElement("Lagrange",mesh_main.ufl_cell(),1)
VPelt=Velt*Pelt
W=FunctionSpace(mesh_main,VPelt)
Velt_alt=VectorElement("Lagrange",mesh_alt.ufl_cell(),2)
Pelt_alt=FiniteElement("Lagrange",mesh_alt.ufl_cell(),1)
VPelt_alt=Velt_alt*Pelt_alt
W_alt=FunctionSpace(mesh_alt,VPelt_alt)
w = Function(W)
(u,p) = w.split(True)
w1 = Function(W)
(u1,p1) = w1.split(True)
w_alt= Function(W_alt)
(u_alt,p_alt) = w_alt.split(True)
w1_alt= Function(W_alt)
(u1_alt,p1_alt) = w_alt.split(True)
xivals=(.09,-.05,.03)
# u is quadratic so interpolations are exact
u_expr=Expression( ("x[0]*x[0]","x[1]*x[0]"), degree=2)
# p is linear so interpolations are exact
p_expr=Expression( "x[0]+x[1]", degree=1)
t=0
for step in range(2):
xi=xivals[step]
# set correct values into u
u.interpolate(u_expr)
p.interpolate(p_expr)
# move mesh_alt
# this will invalidate u_alt data, but u_alt is re-written later
# this SHOULD MAKE NO DIFFERENCE TO u, but it does!
u_alt.interpolate(u_expr)
print('a. step=',step,'u_alt((.2,.8))=',u_alt((.2,.8)),\
'u((.2,.8))=',u((.2,.8)))
print('a. mesh point 25 coords=',mesh_main.coordinates()[25],\
'alt coords=',mesh_alt.coordinates()[25])
mover2(mesh_alt,xi) # this statement changed u!
print('b. mesh point 25 coords=',mesh_main.coordinates()[25],\
'alt coords=',mesh_alt.coordinates()[25])
print('b. u((.2,.8)) SHOULD NOT HAVE CHANGED',u((.2,.8)))
print('b. step=',step,'after mesh_alt moved u_alt((.2,.8))=',u_alt((.2,.8)))
# this makes no difference
mesh_alt.bounding_box_tree().build(mesh_alt)
mesh_main.bounding_box_tree().build(mesh_main)
# interpolate from w to w_alt
u.set_allow_extrapolation(True)
p.set_allow_extrapolation(True)
u_alt.interpolate(u)
p_alt.interpolate(p)
# move primary mesh
mover2(mesh_main,xi)
# this makes no difference
mesh_alt.bounding_box_tree().build(mesh_alt)
mesh_main.bounding_box_tree().build(mesh_main)
# assign back to primary mesh (meshes are same now)
u.assign(u_alt)
p.assign(p_alt)
p1.interpolate(Constant(1.0))
area=assemble(p1*dx)
# compare assigned with exact functions
u1.interpolate(u_expr)
u.set_allow_extrapolation(True)
u_diff=sqrt(assemble(inner(u-u1,u-u1)*dx))/area
p1.interpolate(p_expr)
p.set_allow_extrapolation(True)
p_diff=sqrt(assemble((p-p1)*(p-p1)*dx))/area
print("step=",step,"u interpolation error=",u_diff)
print("step=",step,"p interpolation error=",p_diff)
print("\n")
Output from this code on my computer is
a. step= 0 u_alt((.2,.8))= [0.04 0.16] u((.2,.8))= [0.04 0.16]
a. mesh point 25 coords= [0.2 0.8] alt coords= [0.2 0.8]
b. mesh point 25 coords= [0.2 0.8] alt coords= [0.2 0.8144]
b. u((.2,.8)) SHOULD NOT HAVE CHANGED [0.04 0.16]
b. step= 0 after mesh_alt moved u_alt((.2,.8))= [0.04 0.15712]
Building point search tree to accelerate distance queries.
Computed bounding box tree with 99 nodes for 50 points.
Building point search tree to accelerate distance queries.
Computed bounding box tree with 99 nodes for 50 points.
step= 0 u interpolation error= 2.158268209379147e-16
step= 0 p interpolation error= 2.052474303133418e-16
a. step= 1 u_alt((.2,.8))= [0.04 0.16] u((.2,.8))= [0.04 0.16]
a. mesh point 25 coords= [0.2 0.8144] alt coords= [0.2 0.8144]
b. mesh point 25 coords= [0.2 0.8144] alt coords= [0.2 0.8064]
b. u((.2,.8)) SHOULD NOT HAVE CHANGED [0.04 0.1616]
b. step= 1 after mesh_alt moved u_alt((.2,.8))= [0.04 0.1616]
Building point search tree to accelerate distance queries.
Computed bounding box tree with 99 nodes for 50 points.
step= 1 u interpolation error= 0.004708007363347353
step= 1 p interpolation error= 0.008824209124146306
Please notice that all seems well the first time through the loop, but on the
second time through, the print statements show that
u((.2,.8)) SHOULD NOT HAVE CHANGED [0.04 0.1616]
what should not have changed has, in fact, been changed. In addition
the final two printed error values are much too large for an interpolation,
no matter what the mesh looks like.
I would appreciate any suggestions about what my mistake might be here and how to correct it. Thanks in advance.