Raviart-Thomas implementation

larismartins · November 29, 2022, 11:12pm

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

So far I did the following

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

= 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?

aszaboa · November 30, 2022, 8:45pm

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

larismartins · November 30, 2022, 10:28pm

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

 +  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

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

aszaboa · December 1, 2022, 12:46pm

I think the error is here:

  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.

larismartins · December 1, 2022, 10:29pm

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

mCk = chk(RTm, PPm);

this is regarding

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

aszaboa · December 2, 2022, 9:23am

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

larismartins · December 2, 2022, 8:59pm

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

 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 )

aszaboa · December 4, 2022, 10:21am

Now there is no q in the varf definition …

larismartins · December 4, 2022, 2:02pm

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

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

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

aszaboa · December 4, 2022, 5:10pm

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.

larismartins · December 4, 2022, 5:30pm

it worked this way

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";
  }

Topic		Replies	Views
Execution macro error * <10LinearCombI7MGauche4C_F0E> General Discussion	2	833	May 22, 2022
Program running General Discussion	12	58	October 11, 2024
Im new to FreeFem, can someone clear up what my program is actually doing? General Discussion	2	476	May 4, 2022
Quick question about freefem output General Discussion	1	174	September 4, 2023
Toy Maxwell equation code General Discussion	0	172	February 7, 2024

Raviart-Thomas implementation

Related topics