Integral of a function and the projection of this function are not the same

General Discussion
1

Hi everyone!

I’m trying to prove the following result

Screenshot from 2023-09-22 16-50-20
Screenshot from 2023-09-22 16-50-20779×98 9.05 KB

but FreeFem++ is giving me different values for each DoF of the space P1(K).

This is my code:

load "qf11to25"
//------------------Global Mesh---------------------
int globalH = 2^1;
mesh calP = square(globalH,globalH);

int NbTriangles = calP.nt; //number of triangles
int NbVertices = calP.nv; //number of vertices

fespace V0(calP,P0);
V0 u0;
for(int i = 0; i < NbTriangles; i++){
    u0[][i] = i;
}

//-------------------Mesh info---------------------
//-------------------------vector: nodes in elem K
real[int,int] elemNodesGlobal(NbTriangles,3);
for (int i = 0; i < NbTriangles; i++){
  for(int j = 0; j < 3; j++)
    elemNodesGlobal(i,j) = int(calP[i][j]);
}
//cout << elemNodesGlobal <<"\n";
  

func f = 8*pi*pi*sin(2*pi*x)*sin(2*pi*y);
//f integration order
int forder = 25;

mesh[int] Th(NbTriangles);
for(int K = 0; K < NbTriangles; K++){

  real coordX0, coordX1, coordX2;
  real coordY0, coordY1, coordY2;

  for(int i = 0; i< NbVertices; i++){
    if(elemNodesGlobal(K,0) == i){
      coordX0 = calP(i).x;
      coordY0 = calP(i).y;
    }
    if(elemNodesGlobal(K,1) == i){
      coordX1 = calP(i).x;
      coordY1 = calP(i).y;
     }
     if(elemNodesGlobal(K,2) == i){
       coordX2 = calP(i).x;
       coordY2 = calP(i).y;
     }
  }


  real[int, int] nodes = [[coordX0, coordY0], [coordX1,coordY1],
                           [coordX2, coordY2]];

 //   cout << nodes << "\n";

  border f0(t=0,1){P.x=nodes(2,0)*t + nodes(1,0)*(1-t);P.y=nodes(2,1)*t + nodes(1,1)*(1-t); label=11;};
  border f1(t=0,1){P.x=nodes(0,0)*t + nodes(2,0)*(1-t);P.y=nodes(0,1)*t + nodes(2,1)*(1-t); label=22;};
  border f2(t=0,1){P.x=nodes(1,0)*t + nodes(0,0)*(1-t);P.y=nodes(1,1)*t + nodes(0,1)*(1-t); label=33;};

  Th[K] = buildmesh(f0(1) + f1(1) + f2(1));
  Th[K] = trunc(Th[K], abs(u0 -K) < 1e-5, split=2^1);

  //projection space
    fespace Qh(Th[K], P1);
    Qh qk;

     varf localMassF(f, qk)
     = int2d(Th[K], qforder=forder)(-1*(8*pi*pi*sin(2*pi*x)*sin(2*pi*y))*qk);
    real[int] veclocalMassF = localMassF(0,Qh);

    cout << "\int_K f*qk, qk \in P1(K): " << "\n";
    cout << veclocalMassF << "\n";

    Qh projF = -1*f;

    cout << "\int_K Pi(f)*qk, qk \in P1(K), Pi(.) projecao em P1(K): " << "\n";
    varf localMassprojF(f, qk)
     = int2d(Th[K], qforder=forder)(projF*qk);
    real[int] veclocalMassprojF = localMassprojF(0,Qh);

    cout << veclocalMassprojF << "\n";
 

}

This is the result for one element K:

int_K f*qk, qk in P1(K): 
6
        -0.2267604553   -0.2267604553   -0.0901406829   -1.636619772    -0.9098593171
        -0.9098593171
int_K Pi(f)*qk, qk in P1(K), Pi(.) projecao em P1(K): 
6
        -0.2056167584   -0.2056167584   -2.51807905e-17 -1.23370055     -0.4112335167
        -0.4112335167

Can someone help me? Maybe I’m doing something wrong.

2

Has someone bumped into the same issue as me?

3

I don’t fully understand how you are performing your L2-projection by creating a different FE space. To the best of my knowledge, you need to solve a linear system consisting of a mass matrix and \int f \phi as right hand side, cf. deal.II’s L2-projection. I think that your code might be interpolating the function f instead on the FE mesh and interpolation in general leads to different results than projection.

4

The fespace is the same for each K. So in each element K of the mesh, I should have the same result, right?

5

I hate to hijack the thread but I noticed that you are citing a different FEM code
and curious if you want to start another thread on FEM topics for either 1) visualization of
intermediate results in unobtrusive async way or 2) transform domain
such as FT or laplace on arbitrary mesh. FF supports FFTW but AFAICT
is limited to uniform grid. Given that these linear equations tend to have forms
that more or less conform to polynomials in Laplace domain that seems like
there is a lot of potential there.

For visualization without a bunch of files ( beyond a fifo although UDP is partially
supported ) and a possible ability to assemble distributed data independent of
the computing cluster, I’m working on this ( AFAICT FF doesn’t have a lot for 1D data
display in any case),

https://github.com/mmarchywka/mjmdatascope

6

Hi mike, do you have any ideia of what is happening in my case?

7

I think I downloaded your code but have not looked at it yet. Generally you are projecting
continuous functions onto a finite set of dof’s and want to do it in such a way to minimize
the difference between the original and your reconstruction. Projections aren’t usually
invertible so there is some error or ambiguity you can expect.

8

Yes, but based on the definition:

Screenshot from 2023-09-26 18-38-20
Screenshot from 2023-09-26 18-38-20942×274 22 KB

I should get the same values for the integrals, right?

9

The integrals should be the same, since this is the definition of the L^2 and H^1 projections.
I decided to re-implement the L^2 projection as well in FreeFEM:

// Parameters
int nn = 4;	//mesh quality
real L = 1.; //geometry length
real H = 1.;	//geometry height

// Mesh
mesh Th = square(L*nn, H*nn, [L*x, H*y]);

// Fespace
fespace Uh(Th, P1);
Uh u;
Uh uh;
Uh F = 8*pi*pi*sin(2*pi*x)*sin(2*pi*y);
func f = 8*pi*pi*sin(2*pi*x)*sin(2*pi*y);

// Problem
problem Projection (u, uh)
	= int2d(Th)(
		  u * uh
	)
	- int2d(Th)(
		  f * uh
	)
	//+ on(1, 2, 3, 4, u=0)
	;

// Solve
Projection;

// Check that projection works correctly
for (int i=0; i < Uh.ndof; i++)
{
  uh[][i] = 1.; // choose i.th basis function as test function
  cout << "LHS of L2-projection @ i=" << i << ": " << int2d(Th)(u * uh) << endl;
  cout << "RHS of L2-projection @ i=" << i << ": " << int2d(Th)(f * uh) << endl;
  uh[][i] = 0.;
}
  
// f interpolated on FE Mesh
plot(F, fill=true, value=true);

// f projected on FE Mesh (L2-projection)
plot(u, fill=true, value=true);

// difference between interpolated and projected right hand side function f
Uh diff = F-u;
plot(diff, fill=true, value=true);

If you run this code, then you see that the left hand side (LHS) of the projection problem (\Pi(u),v_h) equals the right hand side (RHS) of the projection problem (u,v_h) for each test function up to machine precision. (As additional information, I added a comparison of the L^2 projection and the interpolation of the function at the DoF values of the FE space.)

1 Like
10

Yes! That works!
But then again, maybe there’s something wrong with the definition on FreeFem++? Or maybe I’m doing something wrong.
Either way, thanksfor your help!