Hi there,
I am trying to interpolate a solution from the boundary of one domain to the boundary of an adjoining domain, as seen in the pictures.
Domain B is the first mesh and domain A is the second, and for illustration they have been separated: domains A and B have a shared boundary but nonmatching mesh nodes. Basically, I want to solve a PDE in domain A, then use that solution of the PDE as a dirichlet boundary condition for domain B.
My question is: is there a way to do this in FEniCSx?
Note: I cannot (as far as I can tell) use subdomains for this applications. My problem requires that domain B has a very steep anisotropy towards its boundary, and there seems to be no way to create this one sided anisotropy in GMSH (domain A is more complicated than this example and cannot use transfinite meshes.
I attach here a minimum working example
import gmsh
from dolfinx import fem, mesh, io
from mpi4py import MPI
import meshio
gmsh.initialize()
# Mesh for domain A
lc = 0.5
gmsh.model.add("domain_A")
gmsh.model.geo.addPoint(0, -4, 0, lc, 1)
gmsh.model.geo.addPoint(3, -4, 0, lc, 2)
gmsh.model.geo.addPoint(3, 4, 0, lc, 3)
gmsh.model.geo.addPoint(0, 4, 0, lc, 4)
gmsh.model.geo.addPoint(0, 2, 0, lc, 5)
gmsh.model.geo.addPoint(1, 2, 0, lc, 6)
gmsh.model.geo.addPoint(1, -2, 0, lc, 7)
gmsh.model.geo.addPoint(0, -2, 0, lc, 8)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 5, 4)
gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 1, 8)
gmsh.model.geo.addCurveLoop([1, 2, 3, 4,5,6,7,8], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
AB_boundary = 1
right_boundary = 2
surface = 1
gmsh.model.addPhysicalGroup(1, [5, 6, 7], AB_boundary)
gmsh.model.addPhysicalGroup(1, [8, 1, 2,3,4], right_boundary)
gmsh.model.addPhysicalGroup(2, [1], surface)
gmsh.model.mesh.generate(2)
gmsh.write("domain_A.msh")
gmsh.open("domain_A.msh")
gmsh.model.setCurrent("domain_A")
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
# Read in mesh
msh_A = meshio.read("domain_A.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh_A = create_mesh(msh_A, "triangle", prune_z=True)
line_mesh_A = create_mesh(msh_A, "line", prune_z=True)
meshio.write("mesh_A.xdmf", triangle_mesh_A)
meshio.write("mt_A.xdmf", line_mesh_A)
# Create an XDMF File
with io.XDMFFile(MPI.COMM_WORLD, "mesh_A.xdmf", "r") as xdmf:
mymesh_A = xdmf.read_mesh(name="Grid")
ct_A = xdmf.read_meshtags(mymesh_A, name="Grid")
mymesh_A.topology.create_connectivity(mymesh_A.topology.dim, mymesh_A.topology.dim-1)
with io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
ft_A = xdmf.read_meshtags(mymesh_A, name="Grid")
V_A = fem.FunctionSpace(mymesh_A, ("CG", 1))
# Function defined on domain A (usually this is the solution to a PDE)
uh = fem.Function(V_A)
uh.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# Mesh for domain B
gmsh.model.add("domain_B")
progr = 0.8 # coefficient of geometric sequence
nx = 10 # number of evenly spaced nodes on boundary
gmsh.model.geo.addPoint(0, -2, 0, lc, 1)
gmsh.model.geo.addPoint(1, -2, 0, lc, 2)
gmsh.model.geo.addPoint(1, 2, 0, lc, 3)
gmsh.model.geo.addPoint(0, 2, 0, lc, 4)
gmsh.model.geo.addPoint(0, 0, 0, lc, 5)
gmsh.model.geo.addPoint(1, 0, 0, lc, 6)
gmsh.model.geo.addLine(1,2, 1)
gmsh.model.geo.addLine(2,6, 2)
gmsh.model.geo.addLine(6,5, 3)
gmsh.model.geo.addLine(5,1, 4)
gmsh.model.geo.addLine(6,3, 5)
gmsh.model.geo.addLine(3,4, 6)
gmsh.model.geo.addLine(4,5, 7)
# Curved loops
gmsh.model.geo.addCurveLoop([1,2, 3, 4], 1)
gmsh.model.geo.addCurveLoop([5,6,7,-3], 2)
BA_boundary = gmsh.model.geo.addCurveLoop([1,2,5,6],3)
# Plane surfaces
surface1 = gmsh.model.geo.addPlaneSurface([1], 1)
surface2 = gmsh.model.geo.addPlaneSurface([2], 2)
# Set transfinite curve
gmsh.model.geo.mesh.setTransfiniteCurve(1, nx, "Progression", progr)
gmsh.model.geo.mesh.setTransfiniteCurve(2, nx, "Progression", -progr)
gmsh.model.geo.mesh.setTransfiniteCurve(3, nx, "Progression", -progr)
gmsh.model.geo.mesh.setTransfiniteCurve(4, nx, "Progression", progr)
gmsh.model.geo.mesh.setTransfiniteCurve(5, nx, "Progression", progr)
gmsh.model.geo.mesh.setTransfiniteCurve(6, nx, "Progression", -progr)
gmsh.model.geo.mesh.setTransfiniteCurve(7, nx, "Progression", -progr)
gmsh.model.geo.mesh.setTransfiniteSurface(1, "AlternateLeft", [1, 2, 6, 5])
gmsh.model.geo.mesh.setTransfiniteSurface(2, "AlternateRight", [5,6,3,4])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(2, [surface1, surface2], 1)
gmsh.model.addPhysicalGroup(1, [BA_boundary], 1)
gmsh.model.mesh.generate(2)
gmsh.write("domain_B.msh")
gmsh.open("domain_B.msh")
gmsh.model.setCurrent("domain_B")
# Read in mesh
msh_B = meshio.read("domain_B.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh_B = create_mesh(msh_B, "triangle", prune_z=True)
line_mesh_B = create_mesh(msh_B, "line", prune_z=True)
meshio.write("mesh_B.xdmf", triangle_mesh_B)
meshio.write("mt_B.xdmf", line_mesh_B)
# Create an XDMF File
with io.XDMFFile(MPI.COMM_WORLD, "mesh_B.xdmf", "r") as xdmf:
mymesh_B = xdmf.read_mesh(name="Grid")
ct_B = xdmf.read_meshtags(mymesh_B, name="Grid")
mymesh_B.topology.create_connectivity(mymesh_B.topology.dim, mymesh_B.topology.dim-1)
with io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
ft_B = xdmf.read_meshtags(mymesh_B, name="Grid")
V_B = fem.FunctionSpace(mymesh_B, ("CG", 1))
for some reason I didn’t think that would work, but it did!