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) ) //
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];
This script is missing a few definitions… Can you add the necessary things to make it run?
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];
}
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"
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];
}
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