# Biharmonic problem with HCT element does not converge

**URL:** <https://community.freefem.org/t/biharmonic-problem-with-hct-element-does-not-converge/1314>\
**Category:** General Discussion\
**Created:** [November 8, 2021, 2:31pm UTC](https://community.freefem.org/t/biharmonic-problem-with-hct-element-does-not-converge/1314 "2021-11-08T14:31:11Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![viviaxenov](https://avatars.discourse-cdn.com/v4/letter/v/439d5e/32.png) [@viviaxenov](https://community.freefem.org/u/viviaxenov)\
**Post date:** [November 8, 2021, 2:31pm UTC](https://community.freefem.org/t/biharmonic-problem-with-hct-element-does-not-converge/1314/1 "2021-11-08T14:31:11Z")

</div>

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:

```auto
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:

```auto
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:

```auto
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

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [November 9, 2021, 9:53am UTC](https://community.freefem.org/t/biharmonic-problem-with-hct-element-does-not-converge/1314/2 "2021-11-09T09:53:46Z")

</div>

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.

```auto
// 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;
}

```