Normal stress in a horizontal line inside of a mesh from top side only

Hi All,

I am solving two problems P1 and P2 in the same mesh. I want to create a horizontal line inside of a mesh and use normal stress (computed only from the top side of the horizontal line) from problem P1 and use it as a load in problem P2. The way I did it was I used the horizontal line to generate a mesh (Th) using buildmesh() and labeled it as 1 (My first question: Does FreeFEM treat my line as a boundary now?) . Then I passed normal stress from problem P1 in problem P2 using:

  • on int1d(Th, 1)()

I am getting good results. But I cannot confirm if the normal stress used in problem 2 are correct (I mean computed from one side only. i.e., form the top side of the horizontal line). Since the stress uses dx(), dy(), so my concern was if those values were computed using elements from the top as well as bottom of the horizontal line.

Does using + on int1d(Th, 1)() treats horizontal line as a boundary and uses only one side to compute my stress, even though the line is inside of the mesh?

I would be grateful if someone could give their thoughts on it.

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 question is “what computes int1d(Th,1)(u)?”.
The answer is that it uses the value of u only from one side: the side opposite to the direction of the normal (like if the normal were pointing outside the domain, and the value would be taken from inside the domain).
The orientation of the normal to the internal boundary [N.x,N.y] is determined by its orientation.
If you write for example
border b1(t=0.,1.){x=t ; y=1. ;label=1;}
then the orientation is from left to right, and the normal points down (N.y=-1), the value of u is taken from the up side.
If you write
border b1(t=0.,1.){x=1.-t ; y=1. ;label=1;} then the orientation is from right to left, the normal points upwards (N.y=1), the value of u is taken from the down side.

However the operators dx() and dy() are applied globally to a finite element function v.
For example if v is P1 on Th, then dx(v) is computed globally, it is P0 on Th.
If you write int1d(Th,1)(dx(v)) it will take the value of dx(v) only on the side determined by the orientation (as above).

1 Like

Thank you so much for your response and helping me out. I was concerned for dx(), but you clarified it well. Thank you again!