Questions regarding about function extrapolation

General
1

Hello, everyone:
The purpose of this post is to discuss the following problem: When solving multiphysics field equations, suppose we have two domains - Domain F is a rectangle [0.0, 1.0]×[0.0, 1.0] and Domain S is [0.0, 1.0]×[-1.0, 0.0] (with non-conforming mesh), sharing a common interface [0.0, 1.0]×[0.0, 0.0].
I need to use a function u_s (the solution of a PDE solved on Domain S, a finite element function) in the equations for Domain F. I’m considering two approaches:
Method A: Define a function u_f1 on Domain F by extrapolating u_s to u_f1.
Method B: Define u_f2 on Domain F by assigning u_s as Dirichlet boundary condition for u_f2.

My question is: Are there any essential differences between these two approaches?

from fenics import *
import mshr
import numpy as np

domain_f = mshr.Rectangle(Point(0.0, 0.0), Point(1.0, 1.0))
domain_s = mshr.Rectangle(Point(0.0, -1.0), Point(1.0, 0.0))
mesh_f = mshr.generate_mesh(domain_f, 30)
mesh_s = mshr.generate_mesh(domain_s, 50)
# mark boundary
class BottomBoundary_f(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)

class RightBoundary_f(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

class TopBoundary_f(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1.0)

class LeftBoundary_f(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

bottom_f = BottomBoundary_f()
right_f = RightBoundary_f()
top_f = TopBoundary_f()
left_f = LeftBoundary_f()
bmf_f = MeshFunction("size_t", mesh_f, mesh_f.topology().dim() - 1)
bmf_f.set_all(0)
bottom_f.mark(bmf_f, 1) 
right_f.mark(bmf_f, 2)   
top_f.mark(bmf_f, 3)   
left_f.mark(bmf_f, 4)   
ds_f = Measure("ds", domain=mesh_f, subdomain_data=bmf_f)

class BottomBoundary_s(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], -1.0)

class RightBoundary_s(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

class TopBoundary_s(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)

class LeftBoundary_s(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

bottom_s = BottomBoundary_s()
right_s = RightBoundary_s()
top_s = TopBoundary_s()
left_s = LeftBoundary_s()
bmf_s = MeshFunction("size_t", mesh_s, mesh_s.topology().dim() - 1)
bmf_s.set_all(0)
bottom_s.mark(bmf_s, 1) 
right_s.mark(bmf_s, 2)   
top_s.mark(bmf_s, 3)   
left_s.mark(bmf_s, 4)   
ds_s = Measure("ds", domain=mesh_s, subdomain_data=bmf_s)
# Function Expression
func01 = Expression(('100*sin(x[0]) + 100*sin(x[1]) + 3.0', '100*x[0]*x[0] + 100*x[1]*x[1] + 3.0'), degree = 5)

V_f = VectorFunctionSpace(mesh_f, "CG", 2)
u_f01 = Function(V_f)
u_f02 = Function(V_f)

V_s = VectorFunctionSpace(mesh_s, "CG", 2)
u_s01 = Function(V_s)
u_s02 = Function(V_s)

u_s01.interpolate(func01)
u_s02.interpolate(func01)
u_s01.set_allow_extrapolation(True)
u_s02.set_allow_extrapolation(True)
u_f01.set_allow_extrapolation(True)
u_f02.set_allow_extrapolation(True)

# Mathod A
u_f01.interpolate(u_s01)
# Mathod B
B_f02 = DirichletBC(V_f, u_s02, bmf_f, 1)
B_f02.apply(u_f02.vector())

# int on ds(1)
print("u_f01*ds(1):", assemble(dot(u_f01, u_f01)*ds_f(1)))
print("u_f02*ds(1):", assemble(dot(u_f02, u_f02)*ds_f(1)))

# vector on ds(1)
v_f = TestFunction(V_f)
v_s = TestFunction(V_s)
form_f01 = assemble(dot(u_s01, v_f)*ds_f(1)).get_local()
form_f02 = assemble(dot(u_f02, v_f)*ds_f(1)).get_local()

bool01 = form_f01 == form_f02
for a01 in bool01:
    if a01:
        continue
    else:
        print(a01)

error = form_f01 - form_f02
print("error:", np.linalg.norm(error, ord = 2))
print("----done---"

The motivation behind my question stems from solving a complex fluid-structure interaction problem, where I observed that when using Method A, the solution fails to converge, whereas Method B yields convergent results.

Additionally, this raises a broader question regarding multiphysics simulations: What are the recommended approaches when one needs to use functions defined on domain B (with non-matching grids) in computations for domain A?

I would be truly grateful if anyone could provide professional advice or suggestions on this matter.

2

As far as I am aware, those two methods should be more or less equivalent in legacy FEniCS.
However, are you experiencing the differences when running in serial or parallel?

Secondly, you can play around with the various methods you can send into “DirichletBC”, i.e. topological", “geometric”, “pointwise” and see if either of them give the same result as interpolation with extrapolation.

Secondly, as you have code for each method, how does the two solutions look when you visualize then, i.e in your script above, do you see clear differences between u_f01 and u_f02 if you save them to file and visualize them with for instance Paraview?