Hello,
I am working on the Stokes equation now. I noticed FreeFEM introduces an extra small L2 term on the pressure to stabilize the scheme. However, I am trying to impose the pressure with a point value to gain the well-posedness of the system. I found the d.o.f of pressure via the following approach (sorry, the forum website does not allow me to upload the file, so I will paste my code here):
//************************************************************************************
fespace Xh(Th,[P2,P2,P1]);
Xh [uind,vind,pind] = [1,2,3];
for(int irow=0;irow<Ahat.n;irow++){//Xh.ndof
//int irow = 2;
if(uind[irow]==3){
/* Change to p = pvalue!
irowInd = irow;
//cout<<"Ahat = "<<Ahat(irow,irow)<<endl;
for(int icol=0;icol<Ahat.n;icol++){
if(abs(Ahat(irow,icol))>1e-8){
Ahat(irow,icol) = 0.0;
}
}
for(int icol=0;icol<Ahat.n;icol++){
if(abs(Ahat(icol,irow))>1e-8){
bhat(icol) = bhat(icol) - Ahat(icol,irow)ppiont(nu, t,ux[][irow],uy[][irow]);
}
}
Ahat(irow,irow) = 1.0;
bhat(irow) = ppiont(nu, t,ux[][irow],uy[][irow]);//tgv
/
// Penalty approach
Ahat(irow,irow) = 1e30;
bhat(irow) = 1e30ppiont(nu, t,ux[irow],uy[irow]);
break;
}
}//Xh.ndof
solhat = Ahat^-1bhat;
uhat[] = solhat;
//***********************************************************************************
Unfortunately, it does not give the exact solution (ppiont(nu, t,ux[irow],uy[irow])) at uhat[irowInd]. As you can see, I also try the other way to impose such boundary condition, but it also fails. Can anyone help me out?