Transfer data on a common interface

General
1

Hello everyone, I am working on FeniCSX 0.6 and I have a question regarding of using a fem function (from mesh01) on mesh02. You can find more details as following:

I have defined two mesh: mesh01 = (0, 1) (0, 1) and mesh02 = (0, 1) (-1, 1). There is a common interface between these two meshes y = 0.

Now, I have got a FEM function u01 (on mesh01) after solving a PDE system, I want to set u01 as the Dirichlet boundary condition on the interface y = 0 for another FEM function u02 (on mesh02).

Is there any convenient way to reach my goal? I know we can do it by setting u01.set_allow_extrapolation(True) on old version. But I failed on fenicsx.

import ufl
import numpy as np
import dolfinx 
from petsc4py import PETSc
from mpi4py import MPI
import dolfinx.fem.petsc

def f_real(x):
    return np.sin(np.pi*x[0]*(1 - x[1])) + np.cos(np.pi*x[1]*(1 - x[0]))

cell_type01 = dolfinx.mesh.CellType.triangle
mesh01 = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,
                                       [[0.0, 0.0], [1.0, 1.0]],
                                       [10, 10],
                                       cell_type01)
x01  = mesh01.geometry.x
tdim01 = mesh01.topology.dim
fdim01 = tdim01 - 1
def bottom01(x): 
    return np.isclose(x[1], 0.0)
def right01(x):  
    return np.isclose(x[0], 1.0)
def top01(x):    
    return np.isclose(x[1], 1.0)
def left01(x):   
    return np.isclose(x[0], 0.0)

bottom_facet01 = dolfinx.mesh.locate_entities_boundary(mesh01, fdim01, bottom01)
right_facet01  = dolfinx.mesh.locate_entities_boundary(mesh01, fdim01, right01)
top_facet01     = dolfinx.mesh.locate_entities_boundary(mesh01, fdim01, top01)
left_facet01    = dolfinx.mesh.locate_entities_boundary(mesh01, fdim01, left01)

facets01 = np.concatenate([bottom_facet01, right_facet01, top_facet01, left_facet01])
values01 = np.concatenate([
    np.full_like(bottom_facet01, 1, dtype=np.int32),
    np.full_like(right_facet01,  2, dtype=np.int32),
    np.full_like(top_facet01,    3, dtype=np.int32),
    np.full_like(left_facet01,   4, dtype=np.int32)
])

sorted_facets01 = np.argsort(facets01)
facet_tag01 = dolfinx.mesh.meshtags(mesh01, fdim01, facets01[sorted_facets01], 
                                    values01[sorted_facets01])

ds01 = ufl.Measure("ds", domain=mesh01, subdomain_data=facet_tag01)


Vh01_scal = dolfinx.fem.FunctionSpace(mesh01, ("CG", 2))
u01 = dolfinx.fem.Function(Vh01_scal)
u01.interpolate(f_real)


mesh02 = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,
                                       [[0.0, -1.0], [1.0, 0.0]],
                                       [50, 50],
                                       cell_type01)
x02  = mesh02.geometry.x
tdim02 = mesh02.topology.dim
fdim02 = tdim02 - 1
def bottom02(x): 
    return np.isclose(x[1], -1.0)
def right02(x):  
    return np.isclose(x[0], 1.0)
def top02(x):    
    return np.isclose(x[1], 0.0)
def left02(x):   
    return np.isclose(x[0], 0.0)

bottom_facet02 = dolfinx.mesh.locate_entities_boundary(mesh02, fdim02, bottom02)
right_facet02  = dolfinx.mesh.locate_entities_boundary(mesh02, fdim02, right02)
top_facet02     = dolfinx.mesh.locate_entities_boundary(mesh02, fdim02, top02)
left_facet02    = dolfinx.mesh.locate_entities_boundary(mesh02, fdim02, left02)

facets02 = np.concatenate([bottom_facet02, right_facet02, top_facet02, left_facet02])
values02 = np.concatenate([
    np.full_like(bottom_facet02, 1, dtype=np.int32),
    np.full_like(right_facet02,  2, dtype=np.int32),
    np.full_like(top_facet02,    3, dtype=np.int32),
    np.full_like(left_facet02,   4, dtype=np.int32)
])

sorted_facets02 = np.argsort(facets02)
facet_tag02 = dolfinx.mesh.meshtags(mesh02, fdim02, facets02[sorted_facets02], 
                                    values02[sorted_facets02])

ds02 = ufl.Measure("ds", domain=mesh02, subdomain_data=facet_tag02)

Vh02_scal  = dolfinx.fem.FunctionSpace(mesh02, ("CG", 2))
u02   = dolfinx.fem.Function(Vh02_scal)

facet_top02 = facet_tag02.find(3)
dofs_top02  = dolfinx.fem.locate_dofs_topological(Vh02_scal, fdim02, facet_top02)


bc_top02 = dolfinx.fem.dirichletbc(u01, dofs_top02, Vh02_scal)
dolfinx.fem.set_bc(u02.x.array, [bc_top02])

a01 = dolfinx.fem.assemble_scalar(dolfinx.fem.form(u02*ds02(3)))
print(f"------{a01}----------")
2

I would:

  1. Update from fenicsx 0.6 as many improvements (especially on non-matching interpolation) has been done in the newer versions.

Then you can use the following:

from mpi4py import MPI
import dolfinx
import numpy as np
import basix.ufl
import ufl

# Define data on mesh 1
mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 20, 20, cell_type=dolfinx.mesh.CellType.triangle
)
el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
V = dolfinx.fem.functionspace(mesh, el)

k = dolfinx.fem.Function(V, name="k")
k.interpolate(lambda x: 2 + np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]))

# Define mesh 2
mesh2 = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD,
    [[-1, 0], [0, 1]],
    [10, 10],
    cell_type=dolfinx.mesh.CellType.quadrilateral,
)
el2 = basix.ufl.element("Lagrange", mesh2.basix_cell(), 1)
V2 = dolfinx.fem.functionspace(mesh2, el2)
k2 = dolfinx.fem.Function(V2, name="k_moved")


# Cells of mesh2 that we want to transfer data to
def interface(x):
    return np.isclose(x[0], 0.0)


facets2 = dolfinx.mesh.locate_entities(mesh2, mesh2.topology.dim - 1, interface)
mesh2.topology.create_connectivity(mesh2.topology.dim - 1, mesh2.topology.dim)
cells2 = dolfinx.mesh.compute_incident_entities(
    mesh2.topology, facets2, mesh2.topology.dim - 1, mesh2.topology.dim
)


interpolation_data = dolfinx.fem.create_interpolation_data(V2, V, cells2, padding=1e-6)
k2.interpolate_nonmatching(k, cells2, interpolation_data)

with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "initial_function.bp", [k]) as writer:
    writer.write(0.0)

with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "transferred_function.bp", [k2]) as writer:
    writer.write(0.0)

yielding

image
image1920×1075 120 KB

1 Like
3

Thank you for your reply, and I kindly wonder which version should I update? Fenics 0.9?

4

We just released 0.10 (yesterday). I would advice to go for that release.