Mixed DG method for Biharmonic equation using FreeFem

mixed_dg_freefem.cpp

/* Program to solve the Biharmonic equation 
  \Delta^2 u = f(x,y) in \Omega
   with u = \Grad u.n = 0 on boundary
*/
include "getARGV.idp"

load "Element_P3dc"
load "Element_P4dc"
load "Element_P4"
verbosity = 0.;
bool debug = true;

macro grad(u) [dx(u), dy(u)]//macro for gradient
macro dn(u)(N.x*dx(u)+N.y*dy(u))//macro for the normal derivative

int n, M, r=4;
int[int] ps = [4,8,16,32];

real[int] error1(r), order(r-1), h1(r), h1err(r), h1order(r-1);


cout<<"P4 Discontinuous Galerkin Method"<<endl;
for(int p=0; p<r; p++){

    n = ps[p];

    real pena = 2*(n)^1; //Change to 2*(n)^3 for superpenalization
    real pena1 = 0;
    mesh Th = square(n,n,flags=0);

    fespace Vh(Th,P4dc); //Change to P1dc/ P2dc etc
    fespace Wh(Th,P4dc);

    Vh uh,f,ue,phi;
    Wh vh,psi;

    f = 24000*x^4*y^4*(x - 1)^4 + 24000*x^4*y^4*(y - 1)^4 + 24000*x^4*(x - 1)^4*(y - 1)^4 + 24000*y^4*(x - 1)^4*(y - 1)^4 + 384000*x*y^4*(x - 1)^3*(y - 1)^4 + 384000*x^4*y*(x - 1)^4*(y - 1)^3 + 192000*x^3*y^4*(2*x - 2)*(y - 1)^4 + 192000*x^4*y^3*(2*y - 2)*(x - 1)^4 + 288000*x^2*y^2*(x - 1)^4*(y - 1)^4 + 768000*x^2*y^3*(x - 1)^4*(y - 1)^3 + 864000*x^2*y^4*(x - 1)^2*(y - 1)^4 + 288000*x^2*y^4*(x - 1)^4*(y - 1)^2 + 768000*x^3*y^2*(x - 1)^3*(y - 1)^4 + 2048000*x^3*y^3*(x - 1)^3*(y - 1)^3 + 768000*x^3*y^4*(x - 1)^3*(y - 1)^2 + 288000*x^4*y^2*(x - 1)^2*(y - 1)^4 + 864000*x^4*y^2*(x - 1)^4*(y - 1)^2 + 768000*x^4*y^3*(x - 1)^2*(y - 1)^3 + 288000*x^4*y^4*(x - 1)^2*(y - 1)^2;

    problem BH([uh,vh],[phi,psi]) = int2d(Th)(grad(uh)'*grad(phi)) + intalledges(Th)(mean(dn(phi))*jump(uh)/nTonEdge + mean(dn(uh))*jump(phi)*(nTonEdge-1)/nTonEdge) - int2d(Th)(vh*phi) + int2d(Th)(grad(vh)'*grad(psi)) + intalledges(Th)(mean(dn(vh))*jump(psi)/nTonEdge + mean(dn(psi))*jump(vh)*(nTonEdge-1)/nTonEdge + pena*jump(uh)*jump(psi)/nTonEdge) - int2d(Th)(f*psi);

    BH;

    plot(uh,fill=debug,wait=!debug,value=debug);
    ue = 1000*(x^4*y^4*(x - 1)^4*(y - 1)^4);
    Vh ve = 1000*(- 12*x^2*y^4*(x - 1)^4*(y - 1)^4 - 32*x^3*y^4*(x - 1)^3*(y - 1)^4 - 12*x^4*y^2*(x - 1)^4*(y - 1)^4 - 32*x^4*y^3*(x - 1)^4*(y - 1)^3 - 12*x^4*y^4*(x - 1)^2*(y - 1)^4 - 12*x^4*y^4*(x - 1)^4*(y - 1)^2);

    Vh h = hTriangle;
    h1(p) = h[].max;

    Vh err = ue - uh;
    error1(p) = sqrt(int2d(Th)((err)^2));
    h1err(p) = sqrt( int2d(Th)((dx(err))^2 + (dy(err))^2 ) + intalledges(Th)(pena*(jump(err))^2) );
    //cout<<"At iteration p = "<<p<<"\nError = "<<error1(p)<<"\t"<<"h = "<<h1(p)<<endl;
    cout<<"L2 error = "<<error1[p]<<endl;

    //n = n*2;
 }

for(int q=0; q<r-1; q++)
 {
     h1order(q) = log(h1err(q)/h1err(q+1))/log(h1(q)/h1(q+1));
     order(q) = log(error1(q)/error1(q+1))/log(h1(q)/h1(q+1));
     //cout<<order(q)<<"\t\t"<<h1order(q)<<endl;
     //cout<<order(q)<<endl;
 }
cout<<order<<endl;
cout<<"\n";