Errors in determine circumference through levelset function

General Discussion
1

Hi, can anyone help with my following question.

I’m trying to get the circumference of an ellipse. I have the following mesh which have two domains. The correct circumference for my ellipse (a = 1.5, b = 1) should be c ~= 7.9, this could be obtained through my following code if I set mesh density factor n = 1. However, if I change the mesh density factor n to 2 (OR 3,4,…), this c will change to c ~= 3.9. Strangely, if I change it to a circle with a = 1 and b = 1, then regardless of the value of n, it always gives me right answer with c ~= 6.28.

I realize that by using a different mesh, say only define the 4 sides and directly mesh the whole domain, I can get the correct ellipse circumference with n >= 2. But how to explain such difference by using different mesh discretization algorithms?

Codes are as follows. Thanks! circumference.edp (648 Bytes)

/* circumference of an ellipse or circle */

/* mesh and FE space /
border bottom(t=-2.0,2.0){x=t;y=-2.0;}
border right(t=-2.0,2.0){x=2.0;y=t;}
border top(t=-2.0,2.0){x=-t;y=2.0;}
border left(t=-2.0,2.0){x=-2.0;y=-t;}
real a = 1.5;
real b = 1.;
real n = 2.;
border ellipse(t=0.,2.pi){x=acos(t);y=bsin(t);}
mesh Th1 = buildmesh(bottom(n20)+right(n20)+top(n20)+left(n20)+ellipse(-n40));
mesh Th2 = buildmesh(ellipse(n40));
mesh Th = Th1 + Th2;
fespace Vh(Th,P1);

/* levelset */
func phio = (x^2./a^2.+y^2./b^2.-1.);

/* circumference */
real c = int1d(Th,levelset=phio)(1.);
cout << "n: " << n << ", c: " << c << endl;

2

Using the above code, I further did a parametric study with parameter “a” changed from 0.6 to 1.7, while all other parameters kept. I compared the circumference provided by FreeFEM++ and theoretic formula, and below are the results.

1
11475×175 30.4 KB

It shows that when the ellipse is close to a circle, then the results by FreeFEM++ in this case is reliable. However, if it is too elliptic, or the difference between a and b is large, then the calculations by FreeFEM++ in this case are not reliable. Is there any idea on this phenomenon? Thanks!

3

Have you solved this problem? I have a similar problem recently, and I hope you can give me some help. My current solution is to give him a small offset, for example, func phio = (x^2./a^2.+y^2./b^2.-1.0001); however, it introduce some errors.

4

Thank to all, you find a bug, I will try correct this, but it is difficult because it is du to round-off.