I am trying to do the following. I am solving a fluid dynamics problem with FreeFEM + PETSc, and then, I would like to calculate the derivatives of the velocity field with a simple finite element projection method (mass matrix unknown derivative = varf(dx(u)v)). In general this seems to work, but near the boundaries (that correspond to walls) the derivative seems ‘noisy’, the accuracy of the derivative calculation is reduced. To enforce boundary conditions in the derivative calculation, I would like to prescribe that the tangential component of the velocity derivative is zero, and enforcing this with a Lagrange-multiplier seems to be the most straightforward method. I tried to implement the method using examples found online, but for some reason the results seem to be worse than before (see the image below). I am a bit stuck finding the issue, any help is appreciated. I attach a MWEish script, which is a modification of the FreeFem-sources/examples/hpddm/navier-stokes-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub example. Any help is appreciated!
I am not sure if I was calling the trunc() or extract()commands on an empty mesh. Maybe the problem is related to the fact that since I use a distributed mesh, some processes do not intersect with the boundary line that I want to extract so maybe.
Nevertheless, I attach a newer version of the code in which the extract()and trunc() commands should be definitely good. The error still persist. Do you have any further guesses why?
I think I figured it out. What works for me is calculating the constraints on the line mesh first, then multiplying it with the 2D-1D interpolation matrix. I attach the corredponsing code segment and the full code. Thanks @prj for the quick replies!
/* -------------------------------------------------------------- */
/* scalar FE space for the derivative calculation */
varf vMv(u1,v1) = int2d(Th)(u1 * v1);
Mv = vMv(Vh,Vh,tgv=-1);
/* line mesh extraction */
int[int] lab2Extract = [1]; // label of the surface
meshL thLG = extract(ThGlobal, label=lab2Extract);
meshL ThL;
{
fespace PhL(thLG, P0);
PhL u = chi(Th);
if(u[].max > 1.0e-2){
cout << "Rank " << mpirank << "/" << mpisize << " has thLG: " << thLG.nt << " triangles, max chi: " << u[].max << endl;
ThL = trunc(thLG,abs(u-1.0)<1.0e-2);
}
}
/* restriction matrix for line FE space */
fespace VhL(ThL,PkV);
matrix th2L = interpolate(VhL,Vh);
cout << "Rank " << mpirank << "/" << mpisize << " has th2L: " << th2L.n << "x" << th2L.m << endl;
Mat B(Mv, restriction = th2L);
Mat restrictionMatrix(B, Mv, th2L);
/* constraints for the tangential derivative */
varf vMvDer1(u1, lambda) = int1d(ThL)(lambda * u1*Tl.x);
varf vMvDer2(u1, lambda) = int1d(ThL)(lambda * u1*Tl.y);
matrix vMvDer1temp = vMvDer1(VhL,VhL);
matrix vMvDer2temp = vMvDer2(VhL,VhL);
Mat MvDer1(B, vMvDer1temp);
Mat MvDer2(B, vMvDer2temp);
/* Calculating the actual constrains with matrix - matrix multiplication */
Mat MvDer1F2L, MvDer2F2L;
MatMatMult(MvDer1, restrictionMatrix, MvDer1F2L);
MatMatMult(MvDer2, restrictionMatrix, MvDer2F2L);
Mat Ablock = [[Mv, 0, MvDer1F2L'],
[0, Mv, MvDer2F2L'],
[MvDer1F2L, MvDer2F2L, 0]];
