# Raviart-Thomas implementation

**URL:** <https://community.freefem.org/t/raviart-thomas-implementation/2149>\
**Category:** General Discussion\
**Created:** [November 29, 2022, 11:12pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149 "2022-11-29T23:12:52Z")\
**Posts on this page:** 11\
**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:** [November 29, 2022, 11:12pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/1 "2022-11-29T23:12:52Z")

</div>

Hi there, I’m trying to implement the following set of equations on FreeFem++

 ![Screenshot from 2022-11-29 20-10-23](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b3a486e612c9296ef2a80cd0f9a1096c8520dc8a.png)

So far I did the following

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

load "Element_Mixte"

func RT = RT1;
func PP = P0;
func PPm = P1edge;

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

  //vectorial space
  fespace Pmd(Th[K], [PP, PP]);
  Pmd [taux, tauy];
  Pmd [phix, phiy];
  
  //space RTm
  fespace RTm(Th[K], RT);
  RTm [sigmaHx, sigmaHy];
  RTm [varphix, varphiy];
  
  //edge space
  fespace PPm(Th[K], PPm);
  PPm p, q;
  
  real[int] edgeLabel = [11,22,33];
  
  matrix mCk;
  for(int e=0; e<3;e++){
    varf chk([sigmaHx, sigmaHy, p], [varphix, varphiy, q])
      = int1d(Th[K],edgeLabel[e])((varphix*N.x + varphiy*N.y)*q);
      
    mCk = chk(RTm, PPm);
  }
  
  real[int] bLH(nDoFl0K);
  for(int i=0;i<nDoFl0K;i++){
    int ii = locationMatrix(K,i);

    bLH(i) = lambdaH[][ii];
  }
  
  matrix mbl;
  for(int e=0; e<3;e++){
    varf lhK(muHk, q)
      = int1d(Th[K],edgeLabel[e])(bLH(e)*q);  
    mbl = lhK(LambdaH,PPm);
    
    //cout<< mbl <<"\n";
    
  }
  
  matrix mDk;
  varf dhK([sigmaHx, sigmaHy], [taux, tauy])
    = int2d(Th[K])((varphix*phix) + (varphiy*phiy));
  mDk = dhK(RTm, Pmd);
  
  fespace VhKaux(Th[K],Pk);
  VhKaux auxuHhK;
  auxuHhK[] = uHhg(K,:);
  real[int] sigmaHmsl = kappa*[dx(auxuHhK), dy(auxuHhK)];
  matrix mh;
  varf hhK(auxuHhK, [taux, tauy])
    =int2d(Th[K])(-1*kappa*((dx(auxuHhK)*phix)+(dy(auxuHhK)*phiy)));
  mh = hhK(Vh, Pmd);
  
      
   
   
   cout << "SigmaH x elem " << K << "\n";
   cout << sigmaHx[] << "\n"; 
   cout << "SigmaH y elem " << K << "\n";
   cout << sigmaHy[] << "\n";
   
   cout << "SigmaH msl " << "\n";
   cout << sigmaHmsl << "\n";

}

```

But I’m getting the following error

```auto
= int1d(Th[K],edgeLabel[e])((varphix*N.x + varphiy*N.y)*q) error operator * <10LinearCombI6MDroit4C_F0E>, <10LinearCombI6MDroit4C_F0E> 

```

Could someone help me find what I did wrong?

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [November 30, 2022, 8:45pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/2 "2022-11-30T20:45:08Z")

</div>

I think in the following line all the terms are from the test function space:  
`varf chk([sigmaHx, sigmaHy, p], [varphix, varphiy, q]) = int1d(Th[K],edgeLabel[e])((varphix*N.x + varphiy*N.y)*q);`

Maybe it should be (?)  
`varf chk([sigmaHx, sigmaHy, p], [varphix, varphiy, q]) = int1d(Th[K],edgeLabel[e])((sigmaHx*N.x + sigmaHy*N.y)*q);`

---

<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:** [November 30, 2022, 10:28pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/3 "2022-11-30T22:28:06Z")

</div>

Thank you! Now it seems I have a different problem.

```auto
 + 0 0*2 0
 + 1 0*2 0
  current line = 652
Exec error : Check BilinearOperator N M
   -- number :1
Exec error : Check BilinearOperator N M
   -- number :1
 err code 8 , mpirank 0

```

the code is now as follows

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

load "Element_Mixte"

real[int] L2errorRT(alloc), errorMSL(alloc);
real[int] H1semierrorRT(alloc);
real[int] H1errorRT(alloc);
real globalSumL2RT = 0, globalSumHdivRT = 0, globalSumH1semiRT = 0;
real globalSumMSL = 0;

real[int] sigmaExato = -1*kappa*[dxu, dyu];

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

func RT = RT1;
func PP = P0;
func PPm = P1edge;

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

  //vectorial space
  fespace Pmd(Th[K], [PP, PP]);
  Pmd [taux, tauy];
  Pmd [phix, phiy];
  
  //space RTm
  fespace RTm(Th[K], RT);
  RTm [sigmaHx, sigmaHy];
  RTm [varphix, varphiy];
  
  //edge space
  fespace PPm(Th[K], PPm);
  PPm p, q;
  
  real[int] edgeLabel = [11,22,33];
  
  matrix mCk;
  for(int e=0; e<3;e++){
    varf chk([sigmaHx, sigmaHy, p], [varphix, varphiy, q])
      = int1d(Th[K],edgeLabel[e])((sigmaHx*N.x + sigmaHy*N.y)*q);
      
    mCk = chk(RTm, PPm);
  }
  
  real[int] bLH(nDoFl0K);
  for(int i=0;i<nDoFl0K;i++){
    int ii = locationMatrix(K,i);

    bLH(i) = lambdaH[][ii];
  }
  
  matrix mbl;
  for(int e=0; e<3;e++){
    varf lhK(muHk, q)
      = int1d(Th[K],edgeLabel[e])(bLH(e)*q);  
    mbl = lhK(LambdaH,PPm);
    
    //cout<< mbl <<"\n";
    
  }
  
  matrix mDk;
  varf dhK([sigmaHx, sigmaHy], [taux, tauy])
    = int2d(Th[K])((sigmaHx*taux) + (sigmaHy*tauy));
  mDk = dhK(RTm, Pmd);
  
  fespace VhKaux(Th[K],Pk);
  VhKaux auxuHhK;
  auxuHhK[] = uHhg(K,:);
  real[int] sigmaHmsl = kappa*[dx(auxuHhK), dy(auxuHhK)];
  
  fespace vVhK(Th[K], [Pk,Pk]);
  vVhK [vsolx, vsoly] = [dx(auxuHhK), dy(auxuHhK)];
  matrix mh;
  varf hhK([vsolx, vsoly], [taux, tauy])
    =int2d(Th[K])(-1*kappa*(vsolx*taux + vsoly*tauy));
  mh = hhK(vVhK, Pmd);
  
      
   
   
   cout << "SigmaH x elem " << K << "\n";
   cout << sigmaHx[] << "\n"; 
   cout << "SigmaH y elem " << K << "\n";
   cout << sigmaHy[] << "\n";
   
   cout << "SigmaH msl " << "\n";
   cout << sigmaHmsl << "\n";

}

cout << "Sigma exato" << "\n";
cout << sigmaExato << "\n";

```

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [December 1, 2022, 12:46pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/4 "2022-12-01T12:46:37Z")

</div>

I think the error is here:

```auto
  matrix mbl;
  for(int e=0; e<3;e++){
    varf lhK(muHk, q)
      = int1d(Th[K],edgeLabel[e])(bLH(e)*q);  
    mbl = lhK(LambdaH,PPm);
    
    //cout<< mbl <<"\n";
    
  }

```

You want a matrix as an input, but the varf is missing `muHk` I think.

---

<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:** [December 1, 2022, 10:29pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/5 "2022-12-01T22:29:44Z")

</div>

According to the code line that its accusing the error, is right here

```auto
mCk = chk(RTm, PPm);

```

this is regarding

```auto
  real[int] edgeLabel = [11,22,33];
  
  matrix mCk;
  for(int e=0; e<3;e++){
    varf chk([sigmaHx, sigmaHy, p], [varphix, varphiy, q])
      = int1d(Th[K],edgeLabel[e])((sigmaHx*N.x + sigmaHy*N.y)*q);
      
    mCk = chk(RTm, PPm);
  }
  

```

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [December 2, 2022, 9:23am UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/6 "2022-12-02T09:23:05Z")

</div>

I see. I think the problem is the dimensions of the FE spaces, e.g, `RT1` is a two-variable vectorial FE space, not a 3 variable one (`RTm [sigmaHx, sigmaHy];` vs. ` varf chk([sigmaHx, sigmaHy, p], [varphix, varphiy, q])`

---

<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:** [December 2, 2022, 8:59pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/7 "2022-12-02T20:59:43Z")

</div>

Yes, but when I put only two variables, I get the following error

```auto
 647 : matrix mCk;
  648 : for(int e=0; e<3;e++){
  649 : varf chk([sigmaHx, sigmaHy], [varphix, varphiy])
  650 : = int1d(Th[K],edgeLabel[e])((sigmaHx*N.x + sigmaHy*N.y)*q) error operator <20CDomainOfIntegration>, <10LinearCombI7MGauche4C_F0E> 
 List of choices 
	 ( <12FormBilinear> : <20CDomainOfIntegration>, <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> )
	 ( <St7complexIdE> : <20CDomainOfIntegration>, <St7complexIdE> )
	 ( <d> : <20CDomainOfIntegration>, <d> )
	 ( <10FormLinear> : <20CDomainOfIntegration>, <10LinearCombI6MDroit4C_F0E> )

 Error line number 650, in file mhmV1pp_corrigido.edp, before token )

```

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [December 4, 2022, 10:21am UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/8 "2022-12-04T10:21:25Z")

</div>

> [@larismartins](#):
>
> `varf chk([sigmaHx, sigmaHy], [varphix, varphiy])`

Now there is no `q` in the varf definition …

---

<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:** [December 4, 2022, 2:02pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/9 "2022-12-04T14:02:09Z")

</div>

For me, the correct way for the bilinear form would be

```auto
matrix mCk;
  for(int e=0; e<3;e++){
    varf chk([sigmaHx, sigmaHy], q)
      = int1d(Th[K],edgeLabel[e])((sigmaHx*N.x + sigmaHy*N.y)*q);
      
    mCk = chk(RTm, PPm);
    cout << mCk << "\n";
  }

```

but it does not work

```auto
647 : matrix mCk;
  648 : for(int e=0; e<3;e++){
  649 : varf chk([sigmaHx, sigmaHy], q)
 Error line number 649, in file mhmV1pp_corrigido.edp, before token )
 Sorry you mixte formulation with and without array 
  current line = 649
Compile error : Sorry you mixte formulation with and without array 
	line number :649, )

```

I’m not find much to go on in the documentation

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [December 4, 2022, 5:10pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/10 "2022-12-04T17:10:19Z")

</div>

Maybe you should also define a vectorial FE space for q, then specify [q1,q2] in the varf, but you do not need to use both q1 and q2.

---

<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:** [December 4, 2022, 5:30pm UTC](https://community.freefem.org/t/raviart-thomas-implementation/2149/11 "2022-12-04T17:30:50Z")

</div>

it worked this way

```auto
matrix mCk;
  for(int e=0; e<3;e++){
    varf chk([sigmaHx, sigmaHy], [q])
      = int1d(Th[K],edgeLabel[e])(q*(sigmaHx*N.x + sigmaHy*N.y));
      
    mCk = chk(RTm, PPm);
    cout << mCk << "\n";
  }

```