// This FreeFEM script computes eigenvalues of the Zaremba Laplacian in the domain Omega, with Dirichlet conditions eveywhere except the top curved boundary
//
//---- A. DECLARATIONS ----
// IMPORTANT - every command ends with the semicolon!
// IMPORTANT - all variables must be declared before(or during) first use
//
real L=1.;	// declare variable L to be real and assign value 1.0 to it; decimal point indicates it is not integer
int npoints=30; // declare variable npoints to be integer and assign value 30 to it
int N=50; // declare variable N  to be integer and assign value 50 to it 
real[int] Evalues(N); // an array of N real numbers parametrised by integers; in FreeFem, indices of an array of length N run from 0 to N-1
real t; // real uninitialised variable used later
//---- END OF DECLARATIONS so far ----
// some more to come later
//
//----  B. BOUNDARY DEFINITIONS ----
// just use parametric curve definitions going in COUNTERCLOCKWISE direction; (t=t0, t1) in lines below means that t changes from t0 to t1 for this piece
//
border d1Omega(t=0, L) { x = pi*t; y = 0; label=1; };
border d2Omega(t=0, 1) { x = pi*L; y = pi*t; label=1; };
border d3Omega(t=pi*L, 0) { x = t; y = pi+t*(pi-t)/pi; label=2; };
border d4Omega(t=pi, 0) { x = 0; y = t; label=1; };
border d5Omega(t=0, 2*pi) { x = pi/3 + (pi/4)*cos(t); y = pi/2 - (pi/4)*sin(t); label=1; }
// "label=...;" parts are essential here
//---- END OF BOUNDARY DEFINITIONS ----
//
//---- C. MESH CREATION USING buildmesh COMMAND ----
// format of the argument: piece_defined_by_border(number_of_mesh_points)+...
//
mesh Th=buildmesh(d1Omega(npoints*L*pi)+d2Omega(npoints*pi)+d3Omega(1.2*npoints*L*pi)+d4Omega(npoints*pi)+d5Omega(2*pi*pi/4*npoints)); //the number of mesh points per side need not be proportional to side length and may not be integer (it's rounded down)
//---- END OF MESH DEFINITIONS ----
//
//---- VISUALISE THE MESH, may be commented out
//
plot(Th,wait=1);
//
//---- D. DECLARE THE FEM SPACE AND FEM VARIABLES ----
//
fespace Vh(Th,P2);
Vh u,v;
//---- END OF FEM SPACE DEFINITIONS ----
//
//---- E. DEFINE THE QUADRATIC FORMS ----
varf q(u,v)=int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v) )+on(1, u=0);
varf b(u,v)=int2d(Th)(u*v);
//---- END OF QUADRATIC FORMS DEFINITIONS ----
//
//----  F. CREATE THE MATRICES ----
matrix S=q(Vh,Vh);
matrix M=b(Vh,Vh);
//---- END OF MATRIX CREATION ----
//
//----  DECLARE THE ARRAY TO HOLD EIGENFUNCTIONS ----
Vh[int] Efunctions(N);
//
//---- G. SOLVE THE PROBLEM ----
int k=EigenValue(S,M,sym=true,value=Evalues,vector=Efunctions);
//
//---- END OF SOLVER ----
//---- H. PRINT THE EIGENVALUES ----
cout << "We asked for " << N << " eigenvalues and computed " << k << " eigenvalues:\n" << Evalues;
//
//---- PLOT THE 6th EIGENFUNCTION ----
plot(Efunctions[5]);
// press '?' on the image to see options for graphics
//
//---- END ----