Does the following script run on your machine?
It should if you have a complete FreeFEM installation, that is FreeFEM + PETSc.
It is parallel, but first, it’s best to check that it runs using the command:
FreeFem++-mpi script.edp -v 0 -wg
You should get the desired results, for example, here with n = 100 in the screenshot.
// Parameters
real sigma = 1; //value of the shift (ignore)
real a = 1./sqrt(3); //length of bottom side
string lengthName = "306090";
int n = 25; //numer of mesh points on each side
int nev = 1; //number of computed eigen values
// Mesh
border a0(t=0., .5){x=0; y=1.-t;} //top half of left side
border a1(t=0., .5){x=0; y=.5-t;} //bottom half of left side
border a2(t=0., .5){x=a*(1.-t); y=t;} //bottom half of hypotenuse
border a3(t=0., .5){x= a*(.5-t); y=.5+t;}//top half of hypotenuse
border a4(t=0., .5){x=a*t; y=0.;} //left half of bottom side
border a5(t=0., .5){x=a*(.5+t); y=0.;} //right half of bottom side
//Freefem stuff. Builds mesh and declares function space. Stuff in this section from example.
mesh Th = buildmesh(a0(n) + a1(n) + a2(n)+ a3(n)+ a4(n)+ a5(n));
plot(Th);
func Pk = P2;
fespace Vh(Th, Pk);
Vh u1, u2;
// Problem
// OP = A - sigma B ; // the shifted matrix
varf op (u1, u2)
= int2d(Th)(
dx(u1)*dx(u2)
+ dy(u1)*dy(u2)
)
+ on(a0,a1,a2,a3,a4,a5, u1=0)
;
varf b ([u1], [u2]) = int2d(Th)(u1*u2)
+ on(a0,a1,a2,a3,a4,a5, u1=0)
; //no boundary condition
load "SLEPc"
Mat OP;
include "macro_ddm.idp"
createMat(Th, OP, Pk)
OP = op(Vh, Vh, tgv = -2,sym=1); //crout solver because the matrix in not positive
matrix loc = b(Vh, Vh, tgv = -20, sym =1);
Mat B(OP, loc);
// important remark:
// the boundary condition is make with exact penalization:
// we put 1e30=tgv on the diagonal term of the lock degree of freedom.
// So take Dirichlet boundary condition just on $a$ variational form
// and not on $b$ variational form.
// because we solve $ w=OP^-1*B*v $
// Solve
real[int] ev(nev); //to store the nev eigenvalue
Vh[int] eV(nev); //to store the nev eigenvector
int k = EPSSolve(OP, B, values=ev, vectors=eV,sparams = "-st_type sinvert -eps_nev " + nev + " -eps_target " + sigma);
// Display & Plot
for (int i = 0; i < k; i++){
u1 = eV[i];
real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1));
real mm = int2d(Th)(u1*u1) ;
cout << "lambda[" << i << "] = " << ev[i] << ", err= " << int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1) << endl;
plot(eV[i], cmm="Eigen Vector "+i+" value ="+ev[i], wait=true, value=true);
}
