How to 1D integrate along a path inside mesh

General Discussion
1

Suppose I have a mesh Th defined in the square from x = 0 to 1, y = 0 to 1. I solve a PDE on this mesh, getting solution u. Now I want to know the line integral of u along the straight path from (0, 1/2) to (1, 1/2), that is, along a horizontal line that goes through the center of the square. What is the correct syntax for this?

I tried defining a border “bb” that runs from (0, 1/2) to (1, 1/2), then calling

   int1d(Th, bb)(u),

but this gives zero. I assume this fails because bb is not an actual border of Th.

2

I think you should use levelset keyword for the integration. something like

int1d(Th, levelset=phi)(...)

see doc here

1 Like
3

Hi. I looked at the documentation, but I didn’t really understand what “levelset” does. How would adding this flag allow me to integrate along a given path?

4

take this example.

func r = sqrt(xx +yy);
real lc = int1d(Th,levelset=r-1.)(1.) ;

here it integrates along the circle of radius one. Basically the phi in
int1d(Th, levelset=phi)(...) is the parametric equation along which you integrate.

5

Hi,

I managed to get two domains: one big domain and inside that big domain, I have a small domain. I need to get line integral around the small domain as well as area integral. I am getting some value as one d line integral but the area integral is zero. Can you suggest what is the best way of doing that? Thankyou.

example
example380×648 14.5 KB

6

see it:

7

Thanks, but in my case, I have the labels available. I have a problem in getting area integral for f(u)dxdy in the smaller domain (black mesh). Mesh Th is defined for the whole domain (green+black) and int2d(Th,1)(f(u)) doen’t give me any result. 1 is the labels of inner boundary (black).

8

I am not sure this is what you are looking for, but if you do a 1d integral, int1d(Th,BoundaryLabel), then BoundaryLabel is the label of your border. From what I see, you do a 2d integral, with a label. the 2d integrals work like this: int2d(Th,RegionLabel) and here RegionLabel is the label of the region, and in general it is different from the boundary label. I guess you have no region with number 1…

If I am correct, you should check the region number. This can be easily done by defining a P0 function equal to region, and check its value:

fespace Uh(Th,P0);
Uh Reg = region;
plot(Reg,fill=1,value=1,wait=true); //to see the value of the region flag