I am trying to compute incompressible steady Stokes equation by using discontinuous Galerkin SIPG method. But in FreeFem++ the operators jump and mean giving error message for using vector valued finite element function (for example jump([ux,uy]) ). How can I use jump or mean operator for vector valued function?
Here is the code:
macro dn(u) (N.xdx(u)+N.ydy(u) ) // def the normal derivative
macro gradx(u) (dx(u))//
macro grady(u) (dy(u))//
real[int] L2error(5); //an array of two values
int i=1;
for (i=1;i<6;i++)
{
mesh Th = square(2^i,2^i); // unite square
plot (Th);
fespace Vh(Th,P2dc); // Discontinous P2 finite element
Vh ux, uy;
Vh vx, vy;
Vh uex1, uex2, f1,f2;
fespace Ph(Th, P1dc);
Ph p=0, pp;
Ph pex;
real pena=10*(2^i); // a paramater to add penalisation
macro grad(z) [dx(z), dy(z)] //
macro Grad(u) [grad(u#x), grad(u#y)] //
macro div(ux,uy) (dx(ux) + dy(uy)) //
varf Ans([ux, uy, p], [vx, vy, pp])=
int2d(Th)(dx(ux)*dx(vx)
+ dy(ux)*dy(vx)
+ dx(uy)dx(vy)
+ dy(uy)dy(vy)
- div(ux,uy) * pp
- div(vx,vy) * p)
+intalledges(Th)(
// loop on all edge of all triangle
// the edge are see nTonEdge times so we / nTonEdge
// remark: nTonEdge =1 on border edge and =2 on internal
// we are in a triange th normal is the exterior normal
// def: jump = external - internal value; on border exter value =0
// mean = (external + internal value)/2, on border just internal value
(mean(Grad(u))[N.x,N.y])’
([jump(vx),jump(vy)])+ (mean(Grad(v))[N.x,N.y])’([jump(ux),jump(uy)])
- pena*([jump(ux),jump(uy)])([jump(vx),jump(vy)]) +mean§([jump(vx),jump(vy)])’
[N.x,N.y]+mean(pp)([jump(ux),jump(uy)])’[N.x,N.y]) / nTonEdge
);
uex1= - 2.0 * x ^ 2 * ( x - 1.0 ) ^ 2 * y * ( y - 1.0 ) * ( 2.0 * y - 1.0 );
uex2= 2.0 * x * ( x - 1.0 ) * ( 2.0 * x - 1.0 ) y ^ 2 * ( y - 1.0 ) ^ 2;
f1 = 2.0
* ( 12.0 * x ^ 2 - 12.0 * x + 2.0 )
* ( 2.0 * y ^ 3 - 3.0 * y ^ 2 + y )
+ 2.0
* ( x ^ 4 - 2.0 * x ^ 3 + x ^ 2 )
* ( 12.0 * y - 6.0 );
f2 = - 2.0
* ( 12.0 * x - 6.0 )
* ( y ^ 4 - 2.0 * y ^ 3 + y ^ 2 )
- 2.0
* ( 2.0 * x ^ 3 - 3.0 * x ^ 2 + x )
* ( 12.0 * y ^ 2 - 12.0 * y + 2.0 );
func g=0;
Vh u,v;
problem stokes([ux,uy,p],[vx,vy,pp] ,solver=UMFPACK) = Ans
- int2d(Th)([f1,f2]’*[vx,vy]);
stokes; // solve DG
}
// Plot
plot([ux, uy], p, wait=1);