-- FreeFem++ v4.12 (Sat 06 May 2023 07:56:14 PM CEST - git v4.12-72-gb0015903)
   file : ./LapDG3d.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //    Discontinous Galerlin Method
    2 : //   based on paper from 
    3 : // Riviere, Beatrice; Wheeler, Mary F.; Girault, Vivette
    4 : // title: 
    5 : // A priori error estimates for finite element 
    6 : // methods based on discontinuous approximation spaces
    7 : //  for elliptic problems.
    8 : //  SIAM J. Numer. Anal. 39 (2001), no. 3, 902--931 (electronic).
    9 : //  ---------------------------------
   10 : //  Formulation given by Vivette Girault
   11 : //  ------ 
   12 : // Author: F. Hecht , december 2003
   13 : // -------------------------------
   14 : //   nonsymetric bilinear form
   15 : //   ------------------------
   16 : //  solve $ -\Delta u = f$ on $\Omega$ and $u= g$ on $\Gamma$
   17 : load "msh3"
   18 : macro grad(u) [dx(u),dy(u),dz(u)]  )  //
   19 : macro dn(u) (N'*grad(u) )  )  //  def the normal derivative 
   20 : int nn=5;
   21 : mesh3 Th = cube(nn,nn,nn); // unite square 
   22 : //savemesh(Th,"cube10.mesh");
   23 : //mesh3 Th("cube10.mesh");
   24 : int[int] labs=labels(Th);
   25 : fespace Vh(Th,P2dc3d);     // Discontinous P2 finite element
   26 : fespace Xh(Th,P2);
   27 : //  if param = 0 => Vh must be P2 otherwise we need some penalisation  
   28 : real pena=1e2; // a paramater to add penalisation 
   29 : func f=1;
   30 : func g=0;
   31 : Vh u,v;
   32 : Xh uu,vv;
   33 : problem A(u,v,solver=sparsesolver) = 
   34 :    int3d(Th)( grad(u)   [dx(u),dy(u),dz(u)] '*grad(v)   [dx(v),dy(v),dz(v)] )
   35 :  + intallfaces(Th)(//  loop on all  edge of all triangle 
   36 :            ( jump(v)*mean(dn(u)    (N'*grad(u)   [dx(u),dy(u),dz(u)]  ) ) -  jump(u)*mean(dn(v)    (N'*grad(v)   [dx(v),dy(v),dz(v)]  ) ) 
   37 :           + pena*jump(u)*jump(v) ) / nElementonB 
   38 : )
   39 :  
   40 : - int3d(Th)(f*v) 
   41 : - int2d(Th)(g*dn(v)    (N'*grad(v)   [dx(v),dy(v),dz(v)]  )   + pena*g*v) 
   42 : 
   43 : ;
   44 : problem A1(uu,vv,solver=sparsesolver) 
   45 : = 
   46 :  int3d(Th)(grad(uu)   [dx(uu),dy(uu),dz(uu)] '*grad(vv)   [dx(vv),dy(vv),dz(vv)] ) - int3d(Th)(f*vv) + on(labs,uu=g);
   47 :  
   48 :  A; // solve  DG
   49 :  A1; // solve continuous
   50 : 
   51 :  real err= sqrt(int3d(Th) ( sqr(u-uu)));
   52 :  cout << " err= "<< err<< " " << u[].linfty << " " << uu[].linfty << endl; 
   53 : plot(u,uu,cmm="Discontinue Galerkin",wait=1,value=1);
   54 : plot(u,cmm="Discontinue Galerkin",wait=1,value=1,fill=1);
   55 : assert(err< 0.01);
   56 :  sizestack + 1024 =5152  ( 4128 )

  -- Cube  nv=216 nt=750 nbe=300 kind= 6
  -- Build Nodes/DF on mesh :   n.v. 216, n. elmt. 750, n b. elmt. 300
     nb of Nodes 750    nb of DoF   7500  DFon=00010
  -- FESpace: Nb of Nodes 750 Nb of DoF 7500
  -- Build Nodes/DF on mesh :   n.v. 216, n. elmt. 750, n b. elmt. 300
     nb of Nodes 1331    nb of DoF   1331  DFon=1100
  -- FESpace: Nb of Nodes 1331 Nb of DoF 1331