Biharmonic problem with HCT element does not converge

General Discussion
1

I’ve been trying to solve the biharmonic problem

Δu = f
u = du/dn = 0 on the border

I tried to verify it on some exact solution. I took the unit circle as the domain.
The exact solution for verification was

u(x, y) = (x^2 + y^2 - 1)^2
f = 64

The code for the testing is as follows:

load "Element_HCT"
load "Morley"

// Circle 
border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
real f = 64.;
func uexact = (x^2 + y^2 - 1.)^2;
// Testing point
real xtest = 0.0;
real ytest = 0.0;
real utest = uexact(xtest, ytest);

ofstream convergence("convergence_HCT.csv");
convergence << "ndof,L2err,xtest,ytest,utest,uexact" << endl;

for(int n = 10; n < 640; n *= 2)
{
    mesh Th = buildmesh(a(n));
    //mesh Th = square(n,n,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
    fespace Vh(Th, HCT);
    Vh [u, ux, uy], [v, vx, vy];
    solve KL1([u, ux, uy], [v, vx, vy]) = 
       int2d(Th)(dxx(u)*dxx(v) + 2.*dxy(u)*dxy(v) + dyy(u)*dyy(v))
        - int2d(Th)(f*v)
        + on(1, u=0, ux=0, uy=0);

    real L2err = sqrt(int2d(Th)((u - uexact)^2));
    real uInterp = u(xtest, ytest);
    convergence << Vh.ndof << ","  << L2err << "," << xtest << "," << ytest << "," << uInterp << "," << utest << endl;
}

so you can see that I solve the problem for finer and finer grids and track the convergence error and the divergence between the analytical and the numerical solution at a certain point of interest. The output for the script is:

ndof,L2err,xtest,ytest,utest,uexact
77,0.128422,0,0,0.87874,1
247,0.0182387,0,0,0.994747,1
935,0.0190109,0,0,1.02565,1
3535,0.0226454,0,0,1.02946,1
13925,0.0263314,0,0,1.03353,1
53521,0.0273733,0,0,1.0347,1

As you can see, the solution doesn’t seem to converge neither in L2-norm, nor in the point of interest.
On the other hand, if P2Morley is used instead of HCT, the convergence is OK:

ndof,L2err,xtest,ytest,utest,uexact
47,0.641925,0,0,1.59914,1
157,0.220036,0,0,1.216,1
609,0.0548751,0,0,1.05266,1
2329,0.0137757,0,0,1.01295,1
9229,0.00338233,0,0,1.00325,1
35573,0.000888189,0,0,1.00087,1

Is something wrong with my test script, or is there really an issue in the HCT element?
I’m using FreeFem++ v4.9

2

I know this problem with the HCT element,

  1. the geometry in FreeFEM is no the circle but a polygone,
    so the Boundary condition are wrong.
  2. the integration point are badly placed see explain in examples/plugin/bilapHCT.edp example.

the correction is

remark: this exemples they is not problem with boundary condition but it is a lucky case.

// Circle 
load "Element_HCT"
load "qf11to25" // for tripleQF function

border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
real f = 64.;
func uexact = (x^2 + y^2 - 1.)^2;
// Testing point
real xtest = 0.0;
real ytest = 0.0;
real utest = uexact(xtest, ytest);
QF2 qfHCT5 = tripleQF(qf5pT);
ofstream convergence("convergence_HCT.csv");
//convergence << "ndof,L2err,xtest,ytest,utest,uexact" << endl;
verbosity= 0; 
for(int n = 10; n < 640; n *= 2)
{
    mesh Th = buildmesh(a(n));
    //mesh Th = square(n,n,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
    fespace Vh(Th, HCT);
    Vh [u, ux, uy], [v, vx, vy];
    solve KL1([u, ux, uy], [v, vx, vy]) = 
       int2d(Th, qft=qfHCT5)(dxx(u)*dxx(v) + 2.*dxy(u)*dxy(v) + dyy(u)*dyy(v))
        - int2d(Th, qft=qfHCT5)(f*v)
        + on(1, u=0, ux=0, uy=0);

    real L2err = sqrt(int2d(Th)((u - uexact)^2));
    real uInterp = u(xtest, ytest);
    convergence << n << " " <<  Vh.ndof << " "  << L2err << " " << xtest << " " << ytest << " " << uInterp << " " << utest << endl;
}
1 Like