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 ?
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
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.
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.
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 ?
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 ?
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.
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.
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
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.
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.
