# Jump condition for elasticity on an interface

**URL:** <https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930>\
**Category:** General Discussion\
**Created:** [May 22, 2025, 2:19pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930 "2025-05-22T14:19:28Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 22, 2025, 2:19pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/1 "2025-05-22T14:19:28Z")

</div>

Dear all,

I would like to add in my variational formulation an expression as the following:

\displaystyle \frac{\beta}{h} \int\_\Gamma [u][v] dS - \frac{\beta}{h} \int\_\Gamma (\lambda {\rm div}(u\_0) I + 2 \mu \nabla^s u\_0)n [v] dS

where [\cdot] denotes the jump across the interface \Gamma, u is the solution, v the test function and u\_0 is the restriction of u in the domain 0 i.e u\_0 = u \lvert\_{\Omega\_0}.

I have two questions:

- How to restrict `intalledges` to \Gamma using the label so that I can use the `jump()` function ?
- For the second term, when integrating over \Gamma, how could I consider only the value of u in the domain 0 while considering the jump of v ?

 ![domain](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/31bf98a7439469b83b90066ef4a812d99d67fc3a.png)

Thank you in advance for your help,

Best regards,

Loïc

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 22, 2025, 6:58pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/2 "2025-05-22T18:58:36Z")

</div>

Dear Loïc,  
Instead of using `intalledges` you can use `int1d(Th,5)`, and `jump` will be available.  
For the second term, you have to take the integral “from the outside”.  
For that you need to have a clockwise orientation of \Gamma.  
About the int1d on an internal boundary of a discontinuous function (and the related orientation issue), see

> [@Normal stress in a horizontal line inside of a mesh from top side only](https://community.freefem.org/t/normal-stress-in-a-horizontal-line-inside-of-a-mesh-from-top-side-only/3717/2):
>
> Hello, Your approach is correct. Using buildmesh you can put a border that is inside the domain. Then it has a label (1 in your case) that can be used to compute int1d(Th,1)(). This “internal border” is considered as a boundary. In particular if you write int1d(Th) without mentioning a label, it will integrate on all boundaries, including the “internal boundary”. When you have a finite element function u that is discontinuous through the internal boundary (for example if u is P0 on Th), your…

But take care that in order to keep the inside region in you mesh (and not exclude it), for a clockwise orientation of \Gamma you will need to apply `buildmesh` with a negative number on border 5, like  
`mesh Th=buildmesh(b1(10)+b2(10)+b3(10)+b4(10)+b5(-20));`  
François.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 22, 2025, 7:35pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/3 "2025-05-22T19:35:33Z")

</div>

Dear @fb77 , thank you for you reply,

Concerning firstly the first term I think `jump` is not available with `int1d`. I get the following error:

> Sorry, no jump, mean, otherside in bilinear term must be in integral of type intalledges, intallVFedges or intallfaces

According to the documentation In a `problem`, `solve` or `varf` definition, the content of `int1d` must be a linear or bilinear form.

Best regards,

Loïc

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 22, 2025, 8:16pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/4 "2025-05-22T20:16:36Z")

</div>

You’re right, indeed it is available only for computing an integral, but not for defining a linear problem to solve.  
Then I would try tu use `intalledges` of the quantity you’re interested in, mutiplied by a cutoff function which is zero except on the location you want.  
I think this can be done with a cutoff function in P0edgedc, it enables to have a different value on each side of edges. But to define correctly this cutoff function is a bit of work, to find which are the involved dof.  
Maybe there is a more convenient way to do.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 23, 2025, 7:07am UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/5 "2025-05-23T07:07:25Z")

</div>

To do this I can use

> fespace Fh(Th, P0edge);  
> Fh ChiE;

and then `ChiE[][ee] = 1.;` if the edge `ee` is at the interface between the two regions.

To test this, I can make a loop on the triangles and then on the three edges, and check if the adjacent triangle belongs to the same region or not.

However, for a triangle `k` and a local numbering `e` (between 1 and 3) of the edge in the triangle, how to recover the global numbering of this edge in the whole mesh, i.e what is the mapping `ee = g(k, e)` ?

I assume that the global numbering of the edges and the numbering of the dof of `P0edge` are the same, isn’t it ?

Best regards,

Loïc,

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 23, 2025, 8:52am UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/6 "2025-05-23T08:52:26Z")

</div>

Maybe just with `ee = Fh(k, e)`. I will try.

Loïc

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 23, 2025, 1:28pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/7 "2025-05-23T13:28:34Z")

</div>

I have succeeded to localize the interface with:

> fespace Eh(Th, P0edge);  
> Eh ChiE;
> 
> int NbTriangles = Th.nt;  
> for (int k = 0; k \< NbTriangles; k++){  
> for (int e = 0; e \< 3; e++){  
> int ee = e;  
> int adjacent = Th[k].adj(ee);  
> if (Th[k].region != Th[adjacent].region){  
> ChiE[Eh(k,e)] = 1; }  
> }}

And I have add the jump at the interface in my variational formulation for linear elasticity.  
I control the intensity of the jump in the normal and tangential direction with the parameters `KN` and `KT`.

However adding the jump **change absolutely nothing in the results**. To test this, I compute the integral of the displacement in the circle, and the value of this integral is the same with or without the jump.

As I have not to much experience, with linear elasticity, could someone check if what I have done is right or make sense ?

I put my code in attachment.

Thank you in advance for your help,

Loïc

[MWE\_elasticity\_jump\_Loic.edp](https://community.freefem.org/uploads/short-url/f4W8kboNvZPOGOKIJtaU8u9dYGd.edp) (3.3 KB)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 23, 2025, 8:20pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/8 "2025-05-23T20:20:17Z")

</div>

Your definition of ChiE to localize the circle is correct.

However your intalledges do nothing because your functions are in P2 on Th, hence they are continuous through the circle, and jump()=0.  
I guess that what you want to do is to have an unknown which is P2 on each region, and discontinuous through the circle.  
For this you need to have two meshes, one for the inside region and one for the outside region, and have separate unknowns u0 and u1, which are P2 on the respective meshes.  
Then you have to solve the coupled problem, for example using a composite space.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 24, 2025, 9:42am UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/9 "2025-05-24T09:42:38Z")

</div>

I get it now. I was thinking that since the properties are discontinuous trough the interface, `jump` would be able to catch this discontinuity.

Now I have two meshes: a square with a hole with label 1, 2, 3, 4, 5 and a full circle with label 5.

I have a setting with composite space that works, however now in my variational formulation I don’t know how to prescribe the jump between the two meshes. I think that I have to do manually the difference between the unknowns u0 and u1, but since it is a two side coupling, I don’t know on which mesh I have to compute the jump.

I would be grateful, if you have a minimal example, to show me how does it work.

Thank you in advance,

Best regards,

Loïc

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 24, 2025, 11:03am UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/10 "2025-05-24T11:03:37Z")

</div>

It is not clear to me what transmission conditions you want to set on the interface between the two regions.  
If it is u\_0=u\_1 and \sigma\_0 N=0 (\sigma\_0 the stress from the region 0),  
it means that you can first solve u\_0 on region 0 with Neumann BC, then solve u\_1 on region 1 with nonhomogeneous Dirichlet condition u\_1=u\_0.

If the conditions are more complicate and really coupled, an example is

> [@Problem in implementation of DG code for Stokes-Darcy interface problem](https://community.freefem.org/t/problem-in-implementation-of-dg-code-for-stokes-darcy-interface-problem/3616/6):
>
> Sorry for the DG implementation the code was wrong. A correct version is [composite-space-DGStokesDarcy.edp](https://community.freefem.org/uploads/short-url/nJQLZPcMuzN3uLBzurdFy3y66jO.edp) (4.4 KB)

The description of the problem and corresponding coupling interface conditions is in the pdf file in the beginning of that discussion.

The interface is considered either as a boundary of the domain above (boundary label 5 of Th1 `int1d(Th1,5)`), either as a boundary of the domain below (boundary label 3 of Th2 `int1d(Th2,3)`), depending if the test function corresponds to the unknown in the domain above or below.  
These interface integrals can involve unknowns from both domains.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [May 25, 2025, 3:28pm UTC](https://community.freefem.org/t/jump-condition-for-elasticity-on-an-interface/3930/11 "2025-05-25T15:28:18Z")

</div>

Yes, the conditions are really coupled (maybe I I didn’t explain the problem properly). I had a look on the document, now, I think I have understand how to solve my problem.

Thanks for you help,

Best regards,

Loïc,
