# Compute Hdiv norm for RT elements

Hi there,

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

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

https://doc.freefem.org/documentation/finite-element.html

it does come up in some examples AFAICT,

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

``````

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

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

Based on what I can gather from this documentation, it seems like you are doing it right?

``````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);
``````

These norms are kicking my assâ€¦ when I do it like this:

``````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)

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

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

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

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

Can you share a MWE of your problem?

Sure! This is the variational formulation

In particular, I can compute the first equation strongly, such as

For the RT0 case, I only have equation (14), then:

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

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];
``````

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

``````load "Element_P3"

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) ) //

//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];

}

``````

Hi, this still doesnâ€™t work. Can you provide a MWE?

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

Sorry about before! Now its working

``````load "Element_P3"

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;

//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];

}

``````

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

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;`

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

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.

1 Like