I simulate ALE FSI problem in 3D with 2nd order meshes. In my framework, a remeshing process is needed to interpolate data from the old mesh to the new mesh in every several steps. Here a strange interpolation result is presented and seek for help on solving this issue.
To demonstrate this issue, two test sets of meshes (mesh1_order[1|2].bp and mesh2_order[1|2].bp) are provided – https://drive.google.com/file/d/1YjC6tCtiJaxQ7Sg7FVakZ8or23QNNMAc/view?usp=sharing. The MWE below interpolates function f1 on mesh1 to f2 on mesh2. The problem occurs on the interpolated result f2.
In the mesh.bp files, mesh and cell tags are stored. The mesh is a cubic domain ranging from [-1, -1, -1] to [1, 1, 1]. A sphere with radius 0.5 is centered at [0, 0, 0]. The cell tags for the sphere and the rest part are 0 and 1 to represent solid and fluid tags, respectively.
The result can be found in the attached figure. I tested it on DOLFINx v0.9.0. The original function f1 is all 0 around and with [0, 1] value in the spherical solid part. To compare the interpolation performance, 1st order and 2nd order meshes are used. I can sense that the wiggling result in 1st order (figure g,h) is due to numerical error with lower resolution on the sphere. But it seems a mistake is in 2nd order result (figure c,d) as [-0.12, 0] values on f2 is computed in the halo region of the sphere. I have no clue why there is a miscalculation on the second order meshes interpolation.
MWE:
# interpolation.py
# Run with `mpirun -n 3 python interpolation.py`
# Test on DOLFINx v0.9.0
import os
os.environ['OMP_NUM_THREADS'] = '1'
from mpi4py import MPI
from dolfinx import mesh, io, fem
import adios4dolfinx
import numpy as np
def assign_xyz(submesh, xyz) -> None:
xlims = []
xmins = []
for d in range(ndims):
geom = submesh.geometry.x[:, d]
xmin = geom.min() if len(geom) > 0 else 1e6
xmax = geom.max() if len(geom) > 0 else -1
xmin = comm.allreduce(xmin, op = MPI.MIN)
xlim = comm.allreduce(xmax, op = MPI.MAX) - xmin
xmins.append(xmin)
xlims.append(xlim)
def exp_xyz(x):
r = []
for d in range(ndims):
r.append((x[d]-xmins[d])/xlims[d])
return r
xyz.interpolate(exp_xyz)
def nonmatching_meshes_interpolation(dest: fem.Function, src: fem.Function):
V_dest = dest.function_space
V_src = src.function_space
tdim = V_src.mesh.topology.dim
dest_mesh_cell_map = V_dest.mesh.topology.index_map(tdim)
num_cells_on_proc = dest_mesh_cell_map.size_local + dest_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype = np.int32)
interpolation_data = fem.create_interpolation_data(V_dest, V_src, cells, padding=1e-8)
dest.interpolate_nonmatching(src, cells, interpolation_data=interpolation_data)
dest.x.scatter_forward()
comm = MPI.COMM_WORLD
ndims = 3
order = 2 # 1 or 2
file1 = f"mesh1_order{order}.bp"
file2 = f"mesh2_order{order}.bp"
mesh1 = adios4dolfinx.read_mesh(file1, comm, engine="BP4")
mesh2 = adios4dolfinx.read_mesh(file2, comm, engine="BP4")
celltags1 = adios4dolfinx.read_meshtags(file1, mesh1, meshtag_name = "celltags", engine = "BP4")
celltags2 = adios4dolfinx.read_meshtags(file2, mesh2, meshtag_name = "celltags", engine = "BP4")
V1 = fem.functionspace(mesh1, ("Lagrange", order, (ndims,)))
V2 = fem.functionspace(mesh2, ("Lagrange", order, (ndims,)))
f1 = fem.Function(V1)
f2 = fem.Function(V2)
submesh1, cellmap1, _, _ = mesh.create_submesh(mesh1, ndims, celltags1.find(0))
submesh2, cellmap2, _, _ = mesh.create_submesh(mesh2, ndims, celltags2.find(0))
xyz1 = fem.Function(fem.functionspace(submesh1, ("Lagrange", order, (ndims,))))
xyz2 = fem.Function(fem.functionspace(submesh2, ("Lagrange", order, (ndims,))))
# make functions on submesh1 and mesh1
assign_xyz(submesh1, xyz1)
n = len(cellmap1)
f1.interpolate(xyz1, cells0=np.arange(n), cells1=cellmap1)
f1.x.scatter_forward()
# interpolate into f2 with nonmatching meshes interpolation
nonmatching_meshes_interpolation(f2, f1)
# make function on submesh2 with identical mesh interpolation
n = len(cellmap2)
xyz2.interpolate(f2, cells0=cellmap2, cells1=np.arange(n))
xyz2.x.scatter_forward()
# visualization
with io.XDMFFile(comm, f"f1_order{order}.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh1)
xdmf.write_function(f1)
with io.XDMFFile(comm, f"f2_order{order}.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh2)
xdmf.write_function(f2)
# assistant visualization
with io.XDMFFile(comm, f"xyz1_order{order}.xdmf", "w") as xdmf:
xdmf.write_mesh(submesh1)
xdmf.write_function(xyz1)
with io.XDMFFile(comm, f"xyz2_order{order}.xdmf", "w") as xdmf:
xdmf.write_mesh(submesh2)
xdmf.write_function(xyz2)
