Interpolate a function from a 2D mesh to a 3D mesh

General
1

Hi,

I am currently working on a composite problem and here are my plans on solving this:

  • First calculate the residual stress field in the inclusion with a 2D axisymmetric formulation;

  • Then interpolate this stress field to the 3D mesh, which is resulting from the revolution of the 2D mesh.

I know if one wants to interpolate between nonmatching meshes with same dimension, one can use nmm_interpolation_data to specify this interpolation data. But can I interpolate a function from 2D to 3D?

Thank you for any help!

First

2

Yes, you can, see for instance: dolfinx/python/test/unit/fem/test_interpolation.py at main · FEniCS/dolfinx · GitHub

3

Hi dokken,

Thank you for your prompt response! I guess I haven’t made myself clear. If I understand the function create_nonmatching_meshes_interpolation_data correctly, it interpolates a function that defines on a 2D mesh into a surface in a 3D mesh, where these two surfaces shares the same coordinates.

However, what I want to do will require this function defined on 2D mesh to be equaly interpolated onto all cross-sections of the 3D mesh. Will it be possible to do so?

4

Extending 2D data in 3D is not trivial.
What I would do is:
map 2D data to surface in 3D, then extend the data.
Of course this depends on how you define your 3D grid, a minimal example of what kind of grids you want to use would be helpful.

Another possibility would be to “move” the slice along the 3D grid and interpolate it at various positions (this would be expensive).

5

Hi dokken, sorry for the late reply and here is a minimal example.

import dolfinx.io.gmshio as gmshio
from mpi4py import MPI
from dolfinx import mesh, fem, io
import pygmsh
import numpy as np

# Create 2d mesh of a unit square
domain_2d = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
# Use pygmsh to generate 1/8 of a cylinder with tetrahedron mesh
with pygmsh.geo.Geometry() as geom:
    p1 = geom.add_point([0.0, 0.0, 0.0], 0.1)
    p2 = geom.add_point([1.0, 0.0, 0.0], 0.1)
    p3 = geom.add_point([1.0, 1.0, 0.0], 0.1)
    p4 = geom.add_point([0.0, 1.0, 0.0], 0.1)

    l1 = geom.add_line(p1, p2)
    l2 = geom.add_line(p2, p3)
    l3 = geom.add_line(p3, p4)
    l4 = geom.add_line(p4, p1)

    cl = geom.add_curve_loop([l1, l2, l3, l4])
    plane = geom.add_plane_surface(cl)

    left, cylinder, lateral = geom.revolve(plane, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0], -0.5 * np.pi)
    geom.generate_mesh()
    pygmsh.write("cylinder.msh")

# Load mesh and meshtags
domain_3d, cell_markers, facet_markers = gmshio.read_from_msh("cylinder.msh", MPI.COMM_WORLD, gdim=3)

# Define function spaces
V_2d = fem.functionspace(domain_2d, ("CG", 1))
V_3d = fem.functionspace(domain_3d, ("CG", 1))

# Define a function in 2d
u_2d = fem.Function(V_2d)
u_expr = lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2
u_2d.interpolate(u_expr)

# Define the interpolation data between non-matching meshes
interpolation_data = fem.create_nonmatching_meshes_interpolation_data(domain_3d._cpp_object, V_3d.element, domain_2d._cpp_object)

# Create the function in 3d function space and interpolate
u_3d = fem.Function(V_3d)
u_3d.interpolate(u_2d, nmm_interpolation_data = interpolation_data)

# Write to .xdmf file and view in paraview:

if not os.path.exists("mve"):
    os.makedirs("mve")

with io.XDMFFile(MPI.COMM_WORLD, "mve/u_2d.xdmf", "w") as file:
    file.write_mesh(domain_2d)
    file.write_function(u_2d)

with io.XDMFFile(MPI.COMM_WORLD, "mve/u_3d.xdmf", "w") as file:
    file.write_mesh(domain_3d)
    file.write_function(u_3d)

The result is as follows

image
image1412×706 40.4 KB

As one can see that this interpolation works perfectly for the plane z=0 (front face). What I would like to do is to interpolate this to all cross-sections of this cylinder, i.e. the function u_3d is axisymmetric with respect to y-axis.