# Compute Hdiv norm for RT elements

**URL:** <https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247>\
**Category:** General Discussion\
**Created:** [January 25, 2023, 5:50pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247 "2023-01-25T17:50:20Z")\
**Posts on this page:** 17\
**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:** [January 25, 2023, 5:50pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/1 "2023-01-25T17:50:20Z")

</div>

Hi there,

I couldnt find anything regarding the Hdiv norm for a RT element in the documentation. Can someone help?

---

<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:** [January 26, 2023, 12:45am UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/2 "2023-01-26T00:45:50Z")

</div>

You can search this for Raviart or did you find this?

[https://doc.freefem.org/documentation/finite-element.html](https://doc.freefem.org/documentation/finite-element.html)

it does come up in some examples AFAICT,

```auto
 find .. -type f | xargs -d \\n grep -i Raviart
../plugin/splitmesh6.edp:fespace RT6(Th6,RT0); // Raviart Thomas ordre 0
../plugin/test-suite.log: 11 : fespace RT6(Th6,RT0); // Raviart Thomas ordre 0
../plugin/splitmesh6.edp.log: 11 : fespace RT6(Th6,RT0); // Raviart Thomas ordre 0

I guess the best way to answer specific questions however is make simple test cases 
or read the source code :)

```

---

<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:** [February 1, 2023, 11:40pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/3 "2023-02-01T23:40:19Z")

</div>

My main difficultie is on how to compute the divergent of a vector valued function  
This is my RTm space

```auto
macro div(u,v) ( dx(u)+dy(v) ) //

func bsigExatoX = -1*kappa*dxu;
func bsigExatoY = -1*kappa*dyu;

func RT = RT0;
fespace RTm(Th[K], RT);
RTm [sigmaHx, sigmaHy];

RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
   
 real rsigmaHx, rsigmaHy;
 for(int i=0; i<sigmaKHx.n; i++){
 rsigmaHx = rsigmaHx + errorRTx[][i]^2;
 rsigmaHy = rsigmaHy + errorRTy[][i]^2;
 }
     
      
L2errorRT[K] = int2d(Th[K])(rsigmaHx + rsigmaHy);
L2diverrorRT[K] = int2d(Th[K])((div(errorRTx, errorRTy))^2);
HdiverrorRT[K] = L2errorRT[K] + L2diverrorRT[K]; 

```

This way of using the macro div, it doesn’t seems right

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [February 3, 2023, 10:26pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/4 "2023-02-03T22:26:53Z")

</div>

Based on what I can gather from [this](https://doc.freefem.org/documentation/finite-element.html#raviart-thomas-element) documentation, it seems like you are doing it right?

```auto
VhRT(Th, RT0);
VhRT [f1h, f2h] = [f1(x, y), f2(x, y)];
// divf = dx(f1h) + dy(f2h);
VhP2(Th, [P2, P2]);
VhP2 [g1h, g2h] = [g1(x, y), g2(x, y)];
// divg = dx(g1h) + dy(g2h);

```

---

<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:** [February 13, 2023, 11:22pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/5 "2023-02-13T23:22:31Z")

</div>

These norms are kicking my ass… when I do it like this:

```auto
RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
   
   real rsigmaHx, rsigmaHy;
   for(int i=0; i<ndofPPm; i++){
    rsigmaHx = rsigmaHx + errorRTx[][i];
    rsigmaHy = rsigmaHy + errorRTy[][i];
   }
        
   L2errorRT[K] = int2d(Th[K])(rsigmaHx^2 + rsigmaHy^2);

```

I have convergence for the L2 norm, but when I do it like this (that seems the right way)

```auto
RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
        
L2errorRT[K] = int2d(Th[K])(errorRTx^2 + errorRTy^2);

```

the L2 norm does not converges.  
Do you know which one is the right way? I can’t find anything in the documentation regarting norms for vector spaces such as Raviar Thomas

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [February 14, 2023, 6:01pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/6 "2023-02-14T18:01:32Z")

</div>

What type of convergence you are referring to?

In the first case, you are integrating a constant over the domain (`rsigmaHx` and `rsigmaHy` are just `real` numbers). I agree that this is not what you want to do.

Maybe this is what you are looking for?  
`real L2norm = errorRTx[].l2;`  
[https://doc.freefem.org/examples/developers.html#array](https://doc.freefem.org/examples/developers.html#array)

---

<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:** [February 14, 2023, 10:28pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/7 "2023-02-14T22:28:04Z")

</div>

L2 convergence for sigma and the L2 convergence for the divergent of sigma  
sigma is a function of the space RT  
for the L2 norm of the divergente I`m doing like this

```auto
L2diverrorRT[K] = int2d(Th[K])((divbsigExato-div(sigmaHx,sigmaHy))^2);

```

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [February 15, 2023, 6:51am UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/8 "2023-02-15T06:51:35Z")

</div>

Can you share a MWE of your problem?

---

<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:** [February 15, 2023, 2:23pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/9 "2023-02-15T14:23:52Z")

</div>

Sure! This is the variational formulation  
 ![posprocessamento](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/0e154bc8e56ed6fb814b00d03fa89b38e39ad81c.png)  
In particular, I can compute the first equation strongly, such as  
 ![posprocessamento2](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/519b735560de52129db278cd967d44c11f94dbec.png)  
For the RT0 case, I only have equation (14), then:

```auto
//-------------Global Mesh----------------
int globalH = 2^coarseMesh;
mesh calP = square(globalH,globalH);

//exact solution
func u = sin(2*pi*x)*sin(2*pi*y);
func dxu = 2*pi*cos(2*pi*x)*sin(2*pi*y);
func dyu = 2*pi*sin(2*pi*x)*cos(2*pi*y);

macro div(u,v) ( dx(u)+dy(v) ) //

load "Element_Mixte"

real[int] L2errorRT(alloc), L2diverrorRT(alloc), HdiverrorRT(alloc);
real[int] L2errormsl(alloc), L2diverrormsl(alloc), Hdiverrormsl(alloc);
real globalSumL2RT = 0, globalSumHdivRT = 0, globalSumL2divRT = 0;
real globalSumHdivMSL = 0;

func bsigExatoX = -1*kappa*dxu;
func bsigExatoY = -1*kappa*dyu;

func divbsigExato= -1*kappa*(-8pi*pi*sin(2*pi*x)*sin(2*pi*y));

real[int] sigmaKHx(nDoFl0K),sigmaKHy(nDoFl0K);

func RT = RT0;
func PPm = P0edge;

for(int K = 0; K < NbTriangles; K++){

  //size of edge
  fespace Ph(calP, PPm);
  varf vlenedge(u,v) = intalledges(calP)(1*v/nTonEdge);
  Ph le;
  le[]= vlenedge(0,Ph);

  //space RTm
  fespace RTm(Th[K], RT);

  RTm [sigmaHx, sigmaHy];

//for para construir o vetor de coordenadas com os valores corretas. cada base associada a uma face, lambda da face0*tamface0
  for(int i=0;i<nDoFl0K;i++){
    int ii = locationMatrix(K,i);

    real bLH = lambdaH[][ii];

    sigmaKHx(i) = bLH*le[][Ph(K,i)];
    sigmaKHy(i) = bLH*le[][Ph(K,i)];

 /* cout << bLH << "\n";
    cout << varphix[][i] << "\n";
    cout << le[][Ph(K,i)] << "\n";*/
  }
  

   sigmaHx[] = sigmaKHx;
   sigmaHy[] = sigmaKHy;

RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
     
      
   L2errorRT[K] = errorRTx[].l2 + errorRTy[].l2;
   L2diverrorRT[K] = int2d(Th[K])((divbsigExato-div(sigmaHx,sigmaHy))^2);
   HdiverrorRT[K] = L2errorRT[K] + L2diverrorRT[K]; 
   

   globalSumL2RT = globalSumL2RT + L2errorRT[K];
   globalSumL2divRT = globalSumH1semiRT + L2diverrorRT[K];
   globalSumHdivRT = globalSumHdivRT + HdiverrorRT[K];

```

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [February 15, 2023, 7:27pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/10 "2023-02-15T19:27:31Z")

</div>

This script is missing a few definitions… Can you add the necessary things to make it run?

---

<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:** [February 23, 2023, 11:02pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/11 "2023-02-23T23:02:28Z")

</div>

```auto
load "Element_P3"
load "Element_PkEdge"

func Pm = P0edge; //polinomial order edge
func Pk = P2; //polinomial order volume

//------------------Global Mesh---------------------
int globalH = 2^coarseMesh;
mesh calP = square(globalH,globalH);

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

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]);
}

//----------------DoF of V0 space-------------------
fespace V0(calP,P0);
V0 u0;
V0 u0x = x, u0y = y;

int nDoFu0 = V0.ndof; //number of DoF
int nDoFu0K = V0.ndofK; //number of DoF per elem

/*cout << "DoF V0 = " << nDoFu0 << "\n";
cout << "Number of DoF per elem V0 = " << nDoFu0K
     << "\n"; */

//function to identify global elements
for(int i = 0; i < NbTriangles; i++){
 u0[][i] = i;
}

//------------DoF of LambdaH space-----------------
fespace LambdaH(calP, Pm);
LambdaH l0;
LambdaH lX = x, lY = y;

int nDoFl0 = LambdaH.ndof;
int nDoFl0K = LambdaH.ndofK;

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^subMesh);

 int NbTrianglesLocal = Th[K].nt;
}

//------------------Flux Recovery---------------------------

macro div(u,v) ( dx(u)+dy(v) ) //

load "Element_Mixte"

//parameters
real kappa = 1.0;

real[int] L2errorRT(alloc), L2diverrorRT(alloc), HdiverrorRT(alloc);
real[int] L2errormsl(alloc), L2diverrormsl(alloc), Hdiverrormsl(alloc);
real globalSumL2RT = 0, globalSumHdivRT = 0, globalSumL2divRT = 0;

func bsigExatoX = -1*kappa*dxu;
func bsigExatoY = -1*kappa*dyu;

func divbsigExato = -1*kappa*(-8*pi*pi*sin(2*pi*x)*sin(2*pi*y));

func RT = RT0;
func PPme = P0edge;

for(int K = 0; K < NbTriangles; K++){

  real[int] edgeLabel = [11,22,33];
  //size of edge
  fespace PPm(calP, PPme);
  varf vlenedge(u,v) = intalledges(calP)(1*v/nTonEdge);
  PPm le;
  le[]= vlenedge(0,PPm);
   

  //space RTm
  fespace RTm(Th[K], RT);
  RTm [sigmaHx, sigmaHy];
  RTm [varphix, varphiy];
  
  int ndofRTm = RTm.ndof;
  int ndofPPm = PPm.ndofK;
  
  real[int] sigmaKHx(ndofPPm),sigmaKHy(ndofPPm);
  for(int i=0;i<ndofPPm;i++){
      
    real bLH = 1.;

    sigmaKHx(i) = bLH*le[][PPm(K,i)];
    sigmaKHy(i) = bLH*le[][PPm(K,i)];
    
  }
  

   sigmaHx[] = sigmaKHx;
   sigmaHy[] = sigmaKHy; 
         
    cout << sigmaHx[] << "\n";
    //cout << sigmaHy[] << "\n";
    
    cout << divbsigExato << "\n";
    cout << dx(sigmaHx) << "\n";
   
   //salto da componente normal tem q ser 0
   //cout << "Jump elem " << K << ": " << jump(sigmaHx*N.x+sigmaHy*N.y) << "\n";

   RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
   
         
   L2errorRT[K] = errorRTx[].l2 + errorRTy[].l2;
   //L2diverrorRT[K] = int2d(Th[K])(div(errorRTx, errorRTy));
   L2diverrorRT[K] = int2d(Th[K])((divbsigExato-div(sigmaHx,sigmaHy))^2);
   HdiverrorRT[K] = L2errorRT[K] + L2diverrorRT[K]; 
   
   /*cout << divbsigExato << "\n";
   cout << div(sigmaHx,sigmaHy) << "\n";*/
   
   //cout << L2errorRT[K] << "\n";
      
   //plot([sigmaHx, sigmaHy], wait=1, fill=1, value=1);
   

   globalSumL2RT = globalSumL2RT + L2errorRT[K];
   globalSumL2divRT = globalSumL2divRT + L2diverrorRT[K];
   globalSumHdivRT = globalSumHdivRT + HdiverrorRT[K];  

}

```

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [February 24, 2023, 5:14pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/12 "2023-02-24T17:14:57Z")

</div>

Hi, this still doesn’t work. Can you provide a MWE?

```auto
 The Identifier coarseMesh does not exist 

 Error line number 9, in file test.edp, before token coarseMesh

  current line = 9
Compile error : 
	line number :9, coarseMesh
error Compile error : 
	line number :9, coarseMesh
 code = 1 mpirank: 0

```

---

<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:** [February 26, 2023, 12:40pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/13 "2023-02-26T12:40:31Z")

</div>

Sorry about before! Now its working

```auto
load "Element_P3"
load "Element_PkEdge"

func Pm = P0edge; //polinomial order edge
func Pk = P2; //polinomial order volume

int coarseMesh = 0;
int subMesh = 0;

//------------------Global Mesh---------------------
int globalH = 2^coarseMesh;
mesh calP = square(globalH,globalH);

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

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]);
}

//----------------DoF of V0 space-------------------
fespace V0(calP,P0);
V0 u0;
V0 u0x = x, u0y = y;

int nDoFu0 = V0.ndof; //number of DoF
int nDoFu0K = V0.ndofK; //number of DoF per elem

/*cout << "DoF V0 = " << nDoFu0 << "\n";
cout << "Number of DoF per elem V0 = " << nDoFu0K
     << "\n"; */

//function to identify global elements
for(int i = 0; i < NbTriangles; i++){
 u0[][i] = i;
}

//------------DoF of LambdaH space-----------------
fespace LambdaH(calP, Pm);
LambdaH l0;
LambdaH lX = x, lY = y;

int nDoFl0 = LambdaH.ndof;
int nDoFl0K = LambdaH.ndofK;

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^subMesh);

 int NbTrianglesLocal = Th[K].nt;
}

//------------------Flux Recovery---------------------------

macro div(u,v) ( dx(u)+dy(v) ) //

int alloc = NbTriangles;

load "Element_Mixte"

//parameters
real kappa = 1.0;

real[int] L2errorRT(alloc), L2diverrorRT(alloc), HdiverrorRT(alloc);
real[int] L2errormsl(alloc), L2diverrormsl(alloc), Hdiverrormsl(alloc);
real globalSumL2RT = 0, globalSumHdivRT = 0, globalSumL2divRT = 0;

func dxu = 2*pi*cos(2*pi*x)*sin(2*pi*y);
func dyu = 2*pi*sin(2*pi*x)*cos(2*pi*y);

func bsigExatoX = -1*kappa*dxu;
func bsigExatoY = -1*kappa*dyu;

func divbsigExato = -1*kappa*(-8*pi*pi*sin(2*pi*x)*sin(2*pi*y));

func RT = RT0;
func PPme = P0edge;

for(int K = 0; K < NbTriangles; K++){

  real[int] edgeLabel = [11,22,33];
  //size of edge
  fespace PPm(calP, PPme);
  varf vlenedge(u,v) = intalledges(calP)(1*v/nTonEdge);
  PPm le;
  le[]= vlenedge(0,PPm);
   

  //space RTm
  fespace RTm(Th[K], RT);
  RTm [sigmaHx, sigmaHy];
  RTm [varphix, varphiy];
  
  int ndofRTm = RTm.ndof;
  int ndofPPm = PPm.ndofK;
  
  real[int] sigmaKHx(ndofPPm),sigmaKHy(ndofPPm);
  for(int i=0;i<ndofPPm;i++){
      
    real bLH = 1.;

    sigmaKHx(i) = bLH*le[][PPm(K,i)];
    sigmaKHy(i) = bLH*le[][PPm(K,i)];
    
  }
  

   sigmaHx[] = sigmaKHx;
   sigmaHy[] = sigmaKHy; 
         
    cout << sigmaHx[] << "\n";
    //cout << sigmaHy[] << "\n";
    
    cout << divbsigExato << "\n";
    cout << dx(sigmaHx) << "\n";
   
   //salto da componente normal tem q ser 0
   //cout << "Jump elem " << K << ": " << jump(sigmaHx*N.x+sigmaHy*N.y) << "\n";

   RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
   
         
   L2errorRT[K] = errorRTx[].l2 + errorRTy[].l2;
   //L2diverrorRT[K] = int2d(Th[K])(div(errorRTx, errorRTy));
   L2diverrorRT[K] = int2d(Th[K])((divbsigExato-div(sigmaHx,sigmaHy))^2);
   HdiverrorRT[K] = L2errorRT[K] + L2diverrorRT[K]; 
   
   /*cout << divbsigExato << "\n";
   cout << div(sigmaHx,sigmaHy) << "\n";*/
   
   //cout << L2errorRT[K] << "\n";
      
   //plot([sigmaHx, sigmaHy], wait=1, fill=1, value=1);
   

   globalSumL2RT = globalSumL2RT + L2errorRT[K];
   globalSumL2divRT = globalSumL2divRT + L2diverrorRT[K];
   globalSumHdivRT = globalSumHdivRT + HdiverrorRT[K];  

}

```

---

<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:** [March 27, 2023, 10:50pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/14 "2023-03-27T22:50:29Z")

</div>

Hi there! Did you get the chance to look at it?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [March 29, 2023, 10:15am UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/15 "2023-03-29T10:15:26Z")

</div>

Yes I see that the code is working, but I still do not exactly understand your question. Could you be more specific? One thing I notice is that you are adding the errors incorrectly. The combined l2 norm is `L2errorRT[K] = sqrt(errorRTx[].l2^2 + errorRTy[].l2^2);`

EDIT: since you have a multi component FE array, the combined norm is actually given directly by `L2errorRT[K] = errorRTx[].l2;`

> [@larismartins](#):
>
> ```auto
> load "Element_P3"
> load "Element_PkEdge"
> 
> func Pm = P0edge; //polinomial order edge
> func Pk = P2; //polinomial order volume
> 
> int coarseMesh = 0;
> int subMesh = 0;
> 
> //------------------Global Mesh---------------------
> int globalH = 2^coarseMesh;
> mesh calP = square(globalH,globalH);
> 
> int NbTriangles = calP.nt; //number of triangles
> int NbVertices = calP.nv; //number of vertices
> 
> 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]);
> }
> 
> //----------------DoF of V0 space-------------------
> fespace V0(calP,P0);
> V0 u0;
> V0 u0x = x, u0y = y;
> 
> int nDoFu0 = V0.ndof; //number of DoF
> int nDoFu0K = V0.ndofK; //number of DoF per elem
> 
> /*cout << "DoF V0 = " << nDoFu0 << "\n";
> cout << "Number of DoF per elem V0 = " << nDoFu0K
> << "\n"; */
> 
> //function to identify global elements
> for(int i = 0; i < NbTriangles; i++){
> u0[][i] = i;
> }
> 
> //------------DoF of LambdaH space-----------------
> fespace LambdaH(calP, Pm);
> LambdaH l0;
> LambdaH lX = x, lY = y;
> 
> int nDoFl0 = LambdaH.ndof;
> int nDoFl0K = LambdaH.ndofK;
> 
> 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^subMesh);
> 
> int NbTrianglesLocal = Th[K].nt;
> }
> 
> //------------------Flux Recovery---------------------------
> 
> macro div(u,v) ( dx(u)+dy(v) ) //
> 
> int alloc = NbTriangles;
> 
> load "Element_Mixte"
> 
> //parameters
> real kappa = 1.0;
> 
> real[int] L2errorRT(alloc), L2diverrorRT(alloc), HdiverrorRT(alloc);
> real[int] L2errormsl(alloc), L2diverrormsl(alloc), Hdiverrormsl(alloc);
> real globalSumL2RT = 0, globalSumHdivRT = 0, globalSumL2divRT = 0;
> 
> func dxu = 2*pi*cos(2*pi*x)*sin(2*pi*y);
> func dyu = 2*pi*sin(2*pi*x)*cos(2*pi*y);
> 
> func bsigExatoX = -1*kappa*dxu;
> func bsigExatoY = -1*kappa*dyu;
> 
> func divbsigExato = -1*kappa*(-8*pi*pi*sin(2*pi*x)*sin(2*pi*y));
> 
> func RT = RT0;
> func PPme = P0edge;
> 
> for(int K = 0; K < NbTriangles; K++){
> 
> real[int] edgeLabel = [11,22,33];
> //size of edge
> fespace PPm(calP, PPme);
> varf vlenedge(u,v) = intalledges(calP)(1*v/nTonEdge);
> PPm le;
> le[]= vlenedge(0,PPm);
>    
> 
> //space RTm
> fespace RTm(Th[K], RT);
> RTm [sigmaHx, sigmaHy];
> RTm [varphix, varphiy];
>   
> int ndofRTm = RTm.ndof;
> int ndofPPm = PPm.ndofK;
>   
> real[int] sigmaKHx(ndofPPm),sigmaKHy(ndofPPm);
> for(int i=0;i<ndofPPm;i++){
>       
> real bLH = 1.;
> 
> sigmaKHx(i) = bLH*le[][PPm(K,i)];
> sigmaKHy(i) = bLH*le[][PPm(K,i)];
>     
> }
>   
> 
> sigmaHx[] = sigmaKHx;
> sigmaHy[] = sigmaKHy; 
>          
> cout << sigmaHx[] << "\n";
> //cout << sigmaHy[] << "\n";
>     
> cout << divbsigExato << "\n";
> cout << dx(sigmaHx) << "\n";
>    
> //salto da componente normal tem q ser 0
> //cout << "Jump elem " << K << ": " << jump(sigmaHx*N.x+sigmaHy*N.y) << "\n";
> 
> RTm [errorRTx, errorRTy] = [bsigExatoX, bsigExatoY] - [sigmaHx, sigmaHy];
>    
>          
> L2errorRT[K] = errorRTx[].l2 + errorRTy[].l2;
> //L2diverrorRT[K] = int2d(Th[K])(div(errorRTx, errorRTy));
> L2diverrorRT[K] = int2d(Th[K])((divbsigExato-div(sigmaHx,sigmaHy))^2);
> HdiverrorRT[K] = L2errorRT[K] + L2diverrorRT[K]; 
>    
> /*cout << divbsigExato << "\n";
> cout << div(sigmaHx,sigmaHy) << "\n";*/
>    
> //cout << L2errorRT[K] << "\n";
>       
> //plot([sigmaHx, sigmaHy], wait=1, fill=1, value=1);
>    
> 
> globalSumL2RT = globalSumL2RT + L2errorRT[K];
> globalSumL2divRT = globalSumL2divRT + L2diverrorRT[K];
> globalSumHdivRT = globalSumHdivRT + HdiverrorRT[K];  
> 
> }
> 
> ```

---

<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:** [April 3, 2023, 11:31pm UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/16 "2023-04-03T23:31:07Z")

</div>

But what about the L2 norm of the divergent? Is it correct?  
Because my problem its that, in theory, the norm converges, but my diverges. I dont know if the problem is in the way I’m computing or if it is something else

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [April 4, 2023, 9:20am UTC](https://community.freefem.org/t/compute-hdiv-norm-for-rt-elements/2247/17 "2023-04-04T09:20:02Z")

</div>

I’m sorry, but I am not familiar enough with RT elements to help with that. Perhaps someone else on the forum can help.

One thing I can say for certain is that your definition of `L2diverror` is missing a square root.