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

**URL:** <https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721>\
**Category:** General Discussion\
**Created:** [September 22, 2023, 3:52pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721 "2023-09-22T15:52:44Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![larismartins](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/larismartins/32/2712_2.png) [@larismartins](https://community.freefem.org/u/larismartins)\
**Post date:** [September 22, 2023, 3:52pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/1 "2023-09-22T15:52:44Z")

</div>

Hi everyone!

I’m trying to prove the following result

 ![Screenshot from 2023-09-22 16-50-20](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/e4f74454d64f5ec2ecd13d2f07baa7e53a01430e.png)  
but FreeFem++ is giving me different values for each DoF of the space P1(K).

This is my code:

```auto
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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![larismartins](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/larismartins/32/2712_2.png) [@larismartins](https://community.freefem.org/u/larismartins)\
**Post date:** [September 26, 2023, 9:31am UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/2 "2023-09-26T09:31:57Z")

</div>

Has someone bumped into the same issue as me?

---

<div class="post-metadata">

**Author:** ![mathmerizing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mathmerizing/32/1941_2.png) [@mathmerizing](https://community.freefem.org/u/mathmerizing)\
**Post date:** [September 26, 2023, 11:44am UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/3 "2023-09-26T11:44:29Z")

</div>

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](https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#ac6b404bf03cb2a742b290421cc2789fe). 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.

---

<div class="post-metadata">

**Author:** ![larismartins](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/larismartins/32/2712_2.png) [@larismartins](https://community.freefem.org/u/larismartins)\
**Post date:** [September 26, 2023, 12:06pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/4 "2023-09-26T12:06:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 26, 2023, 12:32pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/5 "2023-09-26T12:32:17Z")

</div>

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](https://github.com/mmarchywka/mjmdatascope)

---

<div class="post-metadata">

**Author:** ![larismartins](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/larismartins/32/2712_2.png) [@larismartins](https://community.freefem.org/u/larismartins)\
**Post date:** [September 26, 2023, 12:44pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/6 "2023-09-26T12:44:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 26, 2023, 2:59pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/7 "2023-09-26T14:59:12Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![larismartins](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/larismartins/32/2712_2.png) [@larismartins](https://community.freefem.org/u/larismartins)\
**Post date:** [September 26, 2023, 5:38pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/8 "2023-09-26T17:38:54Z")

</div>

Yes, but based on the definition:

 ![Screenshot from 2023-09-26 18-38-20](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/1/159451091790d6133539bf7f7dcf69d4e0a2dec1.png)  
I should get the same values for the integrals, right?

---

<div class="post-metadata">

**Author:** ![mathmerizing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mathmerizing/32/1941_2.png) [@mathmerizing](https://community.freefem.org/u/mathmerizing)\
**Post date:** [September 26, 2023, 7:49pm UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/9 "2023-09-26T19:49:25Z")

</div>

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:

```auto
// 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.)

---

<div class="post-metadata">

**Author:** ![larismartins](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/larismartins/32/2712_2.png) [@larismartins](https://community.freefem.org/u/larismartins)\
**Post date:** [September 27, 2023, 11:35am UTC](https://community.freefem.org/t/integral-of-a-function-and-the-projection-of-this-function-are-not-the-same/2721/10 "2023-09-27T11:35:11Z")

</div>

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!