/usr/share/doc/freefem++/examples/examples++-tutorial/BEM.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 | // compute the solution of a Laplace operator in a Semi infini domain.
// with coupling of Boundary element with periodicity BC in x .
// -------------------------------------------------------------
real eps0=1;
int labup=3, labdown=1;
int nharm= 10; // Number of Harmonique
func ug= max(0.,-(x-0.5)*(x-0.75)); // boundary condition..
macro Grad(u) [dx(u),dy(u)] // eom
real Xmax=1,Ymax=0.3;
int NNx=100,NNy=NNx*Ymax;
mesh Th=square(NNx,NNy,[x*Xmax,y*Ymax]);
fespace Vh(Th,P1,periodic=[[2,y],[4,y]]);
Vh uref; // la solution de reference.
{ // calcule de la solution de reference in Huge Domaine.
mesh Th1=square(NNx,NNx*10,[x*Xmax,10*Xmax*y]); // pour la solution de reference
fespace Uh(Th1,P1,periodic=[[2,y],[4,y]]);
Uh uu,vv;
solve Pref(uu,vv)=int2d(Th1)(eps0*(Grad(uu)'*Grad(vv)))+on(labdown,uu=ug);
uref=uu;
plot(uu,wait=1,cmm=" ref sol / large Th ");
} // pour nettoyer la memoire
plot(uref,wait=1,cmm=" ref sol / Th");
varf vP(u,v)=int2d(Th)(eps0*(Grad(u)'*Grad(v)))+on(labdown,u=ug);
varf vF(u,v)=on(labdown,u=ug);
matrix<complex> A=vP(Vh,Vh); // la matrice sans BEM.
complex[int] b=vF(0,Vh);
Vh<complex> u;
{// for cleanning all local varaible at end of block.
// computation of the matrice BEM
int kdfBEM=0; // nb of DoF on border
int[int] IdfB2Vh(0:Vh.ndof-1); // for numbering IdfB2Vh[i]==i
{
int lbord = labup; // label of the BEM border
varf vbord(u1,v1) = on(lbord,u1=x-10);// negative value ..
real[int] xb=vbord(0,Vh,tgv=1);
sort(xb,IdfB2Vh); // sort of array xb and IdfB2Vh
xb = xb ? 1 : 0;// put 1 if none zero
kdfBEM = xb.sum +0.5; // number of DoF on border
IdfB2Vh.resize(kdfBEM); // IdfB2Vh[i] -> number DoF of border
}
// end of the numbering computation
// so IdfB2Vh[ibem] = iVh where ibem is a df of on bem , and iVh is a df in Vh space.
real perio=Xmax;
complex deuxpii=2*pi*1i;
int n=0;//
// Use of higher order Quadarture formular ...
varf vWn(u,w)=int1d(Th,labup,qforder=10)(exp(-deuxpii*(n)*x)*w);
//complex[int] wn=vWn(0,Vh);// with Vh numbering..
complex[int,int] ABemFull(kdfBEM,kdfBEM);// the full bem matrix in Bem numbering.
ABemFull=0;// set of 0
for ( n=-nharm;n<=nharm;++n)
{
complex[int] wwn(kdfBEM);
complex[int] wn=vWn(0,Vh);
wwn=wn(IdfB2Vh);// wwn(i) = wn(IdfB2Vh(i)) i=0 a wwn.n -1
complex Gs=+2.*pi*abs(n/perio/perio)*eps0;
ABemFull += Gs*wwn*wwn';
}
matrix<complex> ABem=ABemFull(IdfB2Vh^-1,IdfB2Vh^-1); // Build the sparse BEm matrix
// ABem(IdfB2Vh(ib),IdfB2Vh(jb)) = ABemFull(ib,jb)
A = A + ABem;
}// for cleanning all local varaible at end of block. ABem ABemFull
set(A,solver=UMFPACK);
u[]=A^-1*b;
Vh ur=real(u),ui=imag(u);
Vh err=ur-uref;
cout << " err Linty=" << err[].linfty << " / " << uref[].linfty << endl;
plot(ur,uref,wait=1,cmm="ur + uref ");
|