# Some problems about solving coupling problem

**URL:** https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686
**Category:** Legacy DOLFIN
**Created:** [February 23, 2024, 6:54am UTC](https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686 "2024-02-23T06:54:46Z")
**Posts on this page:** 4
**Page:** 1

<div class="post-metadata">

### Author: ![xiaohe](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/xiaohe/32/7777_2.png) [@xiaohe](https://fenicsproject.discourse.group/u/xiaohe)
#### Post date: [February 23, 2024, 6:54am UTC](https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686/1 "2024-02-23T06:54:46Z")

</div>

Hello, everyone.  
I am using fenics 2019.1 to solve a coupling problem, and i have some questions as follows:  
I have two mesh named mesh01 and mesh02 ,there is a public interface **Gamma(marked 999)** between the two mesh. I have solved a pde on the mesh01 and got the solution u1. But I need calculate the integral of u1 on the interface when i solve anthor pde on the mesh02. I wonder if the following approach is correct

```auto
u1.set_allow_extrapolation(True)
# function defined on mesh02
u_1_to_2 = Function(V02) 
u_1_to_2.interpolate(u_1)
or
u_1_to_2 = project(u_1, V02)
# 999 : the public interface
# phi2 : testfunction on FunctionSpace V02 defined mesh02
rhs = u_1_to_2*phi2*ds(999)

```

In addition , is rhere any difference between project and interpolate in above code

---

<div class="post-metadata">

### Author: ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [February 23, 2024, 7:03am UTC](https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686/2 "2024-02-23T07:03:30Z")

</div>

If your meshes are separate objects, you do not really have a common interface, Even if the facets align.  
See for instance

> [@How to represent non-uniform first and second type boundary conditions?](https://fenicsproject.discourse.group/t/how-to-represent-non-uniform-first-and-second-type-boundary-conditions/13654/6):
>
> First of all, your coordinates make no sense, as will never match with Please pay special attention to details before posting code. Your code is also way to complicated, introducing boundary meshes, sub meshes of boundary meshes etc. Here is a minimal working example of what you want to achieve: from dolfin import \* import numpy as np mesh = UnitSquareMesh(4, 4) VT = FunctionSpace(mesh, "CG", 1) T\_data = np.array([0, 50, 100, 150, 0]) T\_coordinates = np.array([[0, 0], [0., 0.25], [0., 0.…

and other related posts, and please add a minimal reproducible example for any further posts

---

<div class="post-metadata">

### Author: ![xiaohe](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/xiaohe/32/7777_2.png) [@xiaohe](https://fenicsproject.discourse.group/u/xiaohe)
#### Post date: [February 23, 2024, 8:04am UTC](https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686/3 "2024-02-23T08:04:48Z")

</div>

Thank you for your response!  
There is a minimal example

```auto
from fenics import *
from mshr import *

domain01 = Rectangle(Point(0.0, 0.0), Point(1.0, 1.0))
domain02 = Rectangle(Point(1.0, 0.0), Point(2.0, 1.0))
mesh01 = generate_mesh(domain01, 0.05)
mesh02 = generate_mesh(domain02, 0.01)
# mark boundary
class Lift02(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

tb01 = Lift02()
bm01 = MeshFunction("size_t", mesh01, mesh01.topology().dim()-1, 0)
bm01.set_all(0)
tb01.mark(bm01, 1)
ds01 = Measure("ds", domain= mesh01, subdomain_data= bm01)

bm02 = MeshFunction("size_t", mesh02, mesh02.topology().dim()-1, 0)
bm02.set_all(0)
tb01.mark(bm02, 2)
ds02 = Measure("ds", domain= mesh02, subdomain_data= bm02)

c01 = Constant(1.0)
V01 = FunctionSpace(mesh01, "P", 2)
u01 = Function(V01)
u01.interpolate(c01)
u01.set_allow_extrapolation(True)
print("mesh01_domain_int", assemble(u01*dx))
print("mesh01_boundary", assemble(u01*ds01))
print("mesh01_boundary1", assemble(u01*ds01(1)))

V02 = FunctionSpace(mesh02, "P", 2)
u02 = Function(V02)
# u02.interpolate(u01)
u02 = project(u01, V02)
print("mesh02_domain_int", assemble(u02*dx))
print("mesh02_boundary", assemble(u02*ds02))
print("mesh02_boundary2", assemble(u02*ds02(2)))

```

In the example, we get function u01 by interpolating Constant(1.0). In my work, we get function u01 from a pde system defined mesh01.  
I have the same question :how to transfer data between the two mesh.

---

<div class="post-metadata">

### Author: ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [February 23, 2024, 10:53am UTC](https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686/4 "2024-02-23T10:53:45Z")

</div>

That is already explained in the referenced post, as you just have to extract the data from the first mesh (and associated coordinates), and use the method described in that points. This can easily be done with dolfin.Function.vectot().get\_local() and dolfin.FunctionSpace.tabulate\_dof\_coordinates()
