# How to impose Robin boundary conditions for a eigenproblem

**URL:** <https://community.freefem.org/t/how-to-impose-robin-boundary-conditions-for-a-eigenproblem/2240>\
**Category:** General Discussion\
**Created:** [January 22, 2023, 8:10pm UTC](https://community.freefem.org/t/how-to-impose-robin-boundary-conditions-for-a-eigenproblem/2240 "2023-01-22T20:10:21Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Christian](https://avatars.discourse-cdn.com/v4/letter/c/b2d939/32.png) [@Christian](https://community.freefem.org/u/Christian)\
**Post date:** [January 22, 2023, 8:10pm UTC](https://community.freefem.org/t/how-to-impose-robin-boundary-conditions-for-a-eigenproblem/2240/1 "2023-01-22T20:10:21Z")

</div>

Hi everyone!  
I need to solve the diffusion neutron equaiton in a 2D geometry; In particular i want to solve the eigenproblem. Here below the mathematical description of the problem is shown:

![Screenshot_20230122_203515](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/8/87a7f51acf4ad55c161ce0daf91f1992897c6a9c.png)

 ![Screenshot_20230122_203525](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/350197511acddbc8124ff00c4d618b72a45f4166.png)

As can be seen from the formulation, the energy is discretized in “G” groups, and the parameter “K\_eff” is the eigenvalue.

At this point i would like to impose an “incoming current” = 0 and the icoming current is defined as:

![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/cfcc76d71552d594a4046ca049c5bfbef4ec55ca.png)

And the vector “J\_g” is defined as:  
 ![Screenshot_20230122_205303](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/bfac6dba9f2f182b18f98906714263beab81eb09.png)

so the final boundary condition should be:

![Screenshot_20230122_210202](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/eba940b8176d0956ebb9776ef46dd2ad80fa341c.png)

How can i impose this kind of boundary condition for solving the eigenproblem by using the “varf” formulation?

Here is reported the variational form of the problem in my script (NOTICE: here zero FLUX condition is imposed, so a Dirichlet BC, but it is not correct)

```auto

func Pk = [P1,P1,P1]; 
fespace Vh(coremsh,P1);
fespace VhV(coremsh,Pk);
Vh phi1, phi2, phi3, v1, v2, v3, fisspow, phiTOT;
macro div(D,phi,v) D*(dx(phi)*dx(v) + dy(phi)*dy(v)) // end of macro

varf righthandside([phi1,phi2,phi3],[v1,v2,v3])
 = int2d(coremsh) (CHI[0]*NSF[0]*phi1*v1+CHI[0]*NSF[1]*phi2*v1+CHI[0]*NSF[2]*phi3*v1+CHI[1]*NSF[0]*phi1*v2
   +CHI[1]*NSF[1]*phi2*v2+CHI[1]*NSF[2]*phi3*v2+CHI[2]*NSF[0]*phi1*v3+CHI[2]*NSF[1]*phi2*v3+CHI[2]*NSF[2]*phi3*v3);

varf lefthandside([phi1,phi2,phi3],[v1,v2,v3]) 
= int2d(coremsh) (div(DFC[0],phi1,v1)+RXS[0]*phi1*v1+div(DFC[1],phi2,v2)+RXS[1]*phi2*v2+div(DFC[2],phi3,v3)+RXS[2]*phi3*v3 
- (S1[0]*v1*phi1+S2[0]*v1*phi2+S3[0]*v1*phi3 + S1[1]*v2*phi1+S2[1]*v2*phi2+S3[1]*v1*phi3 + S1[2]*v3*phi1+S2[2]*v3*phi2+S3[2]*v3*phi3))
     + on(0, phi1=0,phi2=0,phi3=0);

real sigma=0;
int nev=1;
real[int] ev(nev); // to store nev eigenvalues
//Vh[int] Vev(nev);
VhV[int] [e1,e2,e3](nev); // to store nev eigenvectors

matrix A = lefthandside(VhV,VhV, tgv = tgv,solver=GMRES);
matrix B = righthandside(VhV,VhV, tgv = tgv,solver=GMRES);

int kk = EigenValue(A,B, sym=false, sigma=0, value=ev, vector=e1, maxit=0, ncv=0,tol=1e-20);

```

---

<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:** [January 25, 2023, 10:19am UTC](https://community.freefem.org/t/how-to-impose-robin-boundary-conditions-for-a-eigenproblem/2240/2 "2023-01-25T10:19:00Z")

</div>

I think you need to use the `int1d` command as in the thermal conduction example in the documentation.

---

<div class="post-metadata">

**Author:** ![Christian](https://avatars.discourse-cdn.com/v4/letter/c/b2d939/32.png) [@Christian](https://community.freefem.org/u/Christian)\
**Post date:** [February 2, 2023, 8:34am UTC](https://community.freefem.org/t/how-to-impose-robin-boundary-conditions-for-a-eigenproblem/2240/3 "2023-02-02T08:34:04Z")

</div>

Thanks for your answer!  
Yes, it is the right command to use for Robin’s BC; It finally turned out that the problem i was facing was about the syntax.

The BC’s expression was something like that :

```auto
varf problemexpression = int2d(Th) (...) +int1d(Th,label_i)(1./2.*phi*w);

```

But it did not work. I finally undestood that, for some reasons, FF seems to reject the usage of “1./2.” directilly inside the BC’s expression, so it is necessary to previously define a real variable, like:

```auto
real coeff=1./2.;
varf problemexpression = int2d(Th) (...) +int1d(Th,label_i)(coeff*phi*w);

```

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [February 9, 2023, 8:57am UTC](https://community.freefem.org/t/how-to-impose-robin-boundary-conditions-for-a-eigenproblem/2240/4 "2023-02-09T08:57:25Z")

</div>

remark, the eigen value problem must be defined on vector space so your Robin condition must homogeneous.

un exemple:

```auto
verbosity=1;
mesh Th=square(20,20,[pi*y,pi*x]);
fespace Vh(Th,P2);
Vh u1,u2;

real sigma = 00; // value of the shift 

varf a(u1,u2)= int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma* u1*u2 )
          + int1d(Th,1)(u1*u2)
                    + on(2,3,4,u1=0) ; // Boundary condition
                   
varf b([u1],[u2]) = int2d(Th)( u1*u2 ) ; // no Boundary condition

matrix A= a(Vh,Vh,solver=sparsesolver); 
matrix B= b(Vh,Vh,solver=CG,eps=1e-20); 

// important remark:
// the boundary condition is make with exact penalisation:
// we put 1e30=tgv on the diagonal term of the lock degre of freedom.
// So take dirichlet boundary condition just on $a$ variationnal form
// and not on $b$ variationnanl form.
// because we solve
// $$ w=A^-1*B*v $$

int nev=20; // number of computed eigen valeu close to sigma

real[int] ev(nev); // to store nev eigein value
Vh[int] eV(nev); // to store nev eigen vector

int k=EigenValue(A,B,sym=true,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=0,ncv=0,which="LM"); // which="LM" default 
// tol= the tolerace
// maxit= the maximal iteration see arpack doc.
// ncv see arpack doc.
// the return value is number of converged eigen value.

```