Strange behavior in the BEM solution of the Laplace equation

General Discussion
1

Preface: I cannot currently attach my code as I am a new user. I will provide the code when this restriction is lifted.

I am trying to use the boundary element method (BEMTool) interfaced with FreeFem to solve the 2d Laplace problem with Dirichlet boundary conditions. I have two samples of code that behave in strange ways. Note that the code (hopefully) provided below has to be run using mpi.

  1. The first file (BEM_disk_problem.edp) solves the Laplace problem on either the unit disk or an ellipse close to the unit disk. The two cases are toggled by uncommenting either line 7 or 8 of the code. The solution is plausible on the ellipse and nonsensical on the disk.

  2. The second file (BEM_oscillation_problem.edp) solves the Laplace problem on an exotic polygon. The solution seems reasonable at first glance. However, under closer examination, it becomes apparent that the solution has multiple small scale oscillations in the interior of the domain, which shouldn’t be the case for harmonic functions. In order to illustrate this behavior, the code plots a mesh adapted to the solution. Notice that is pockmarked by small regions of refined triangles. Each of those corresponds to a spurious extremum.

It is possible that my problems arise from a misunderstanding of the interface between FreeFem and BEMTool. Specifically, it is possible that I am using the wrong kernel and/or potential. I
did my best following the instuctions found at:
https://doc.freefem.org/documentation/BEM.html
That being said, I think that the code example given at the bottom of that page does not work on the current version of FreeFem (I couldn’t run it). Instead, my code is a modification of the example Helmholtz_circle_Dirichlet.edp that comes with FreeFem. There, the syntax for defining a BEM kernel is slightly different from the one found in the documentation. This seems to be the main ingredient to make the code run. However, it is no longer clear that setting k=0 in the kernel produces the kernel for the Laplace problem.

2

What is your FreeFEM version ?
Both FreeFem-sources/Helmholtz_circle_Dirichlet.edp at master · FreeFem/FreeFem-sources · GitHub and the example here should work with FreeFEM 4.7-1 or newer. Consider updating to the latest version from Releases · FreeFem/FreeFem-sources · GitHub

3

Hello Dear Professor;
How can we combine the two boundary conditions (Dirichlet and Neumann)

4

I was indeed using an old version of FreeFem. I have updated my installation to the latest version. All of the examples seem to run correctly now.

On the other hand, I cannot seem to impose boundary conditions with different functions on each boundary component. Previously, the following (severely abridged) code produced the desired outcome (up to the issues described in my original post):

meshL Th = buildmeshL(Gamma1(numpt) + Gamma2(numpt) + Gamma3(numpt) + Gamma4(numpt) );
varf vmassFirstKind(u,v) = int1d(Th,1)((BC1(x,y))*v)+ int1d(Th,2)((BC2(x,y))*v)+ int1d(Th,3)((BC3(x,y))*v)+ int1d(Th,4)((BC4(x,y))*v);
bFirstKind[] = vmassFirstKind(0,Uh);
uFirstKind[] = HFirstKind^-1*bFirstKind[];

In the above, the GammaI are boundaries, BCi are functions and HFirstKind is an Hmatrix. Currently, the above results in an all NaN solution. Changing the definition of varf to the following fixes this issue (but solves the wrong mathematical problem):

varf vmassFirstKind(u,v) = int1d(Th)((BC1(x,y))*v);

What is the correct syntax to integrate with respect to a boundary component?

(Note: I believe that the problems that I found using my original installation still remain. Once I manage to get my code to work I would like to post my code that illustrates those issues.)

5

I am guessing your bFirstKind[].linfty is zero. The trick is that a meshL is not always a boundary of another mesh, it is a mesh on its own. Thus, when integrating over part of the mesh, as in the 2D or 3D case, it is the region number that is important, not the label number. So when you write int1d(Th,1)((BC1(x,y))*v), it will integrate over all elements of region 1 ; which means that the region attribute of your mesh elements has to be defined, and not the label attribute (which is used for the boundary of a mesh).
Thus, when defining your Gamma1, Gamma2 … borders to build the meshL, you have to set the region attribute, not the label attribute. For example:

load “msh3”
border a(t=0,pi){x=cos(t);y=sin(t);region=1;}
border b(t=pi,2*pi){x=cos(t);y=sin(t);region=2;}
meshL Th = buildmeshL(a(50)+b(50));
cout << int1d(Th)(1) << " " << int1d(Th,1)(1) << " " << int1d(Th,2)(1) << endl;

Although I agree that the label attribute of the border could be translated into the region number when building a meshL. We will see about that.

6

Sorry, it is currently not possible to easily impose mixed boundary conditions for BEM (Dirichlet on a part of the boundary, Neumann on another part) ; I am guessing that is what you want to do ?
We will work on making mixed boundary conditions possible and add an example in the distribution when it is done, hopefully in the near future.

7

Thank you for your help. My code now executes as previously, but on the new version of FreeFem. Thus, we can now discuss the issues I initially wanted to bring to your attention. In order to do this, I need the permission to post attachments. What is the procedure for this?

8

Unfortunately there is no easy solution to grant you permission to post attachments. The easiest way would be that you use Preformatted text to embed your code in the post (put your code between ```
```), or host it on any hosting website and link it here if it really is too long. Sorry for this inconvenience.

9

Here is the first program. In order to see the the strange behavior, toggle between the lines defining stretchFactor.

load "bem"
load "msh3"

// CHANGE THESE SETTINGS TO SEE THE PROBLEM (uncomment the relevent line)

real stretchFactor = 1; // This produces a disk. THIS PRODUCES A STRANGE BEHAVIOR.
// real stretchFactor = 1.01; // This produces an ellipse. THIS BEHAVES CORRECTLY.

// SETTINGS END HERE.

int k = 0;

func real heaviside(real input){ //Boundary condition.
	
	if(input<0){
		return 0;}
	else{
	return 1;}
};

func bc = heaviside(x);

//  Mesh
int n = 1000;
border circle(t=0, 2*pi){x=stretchFactor*cos(t); y=sin(t);} // <---------- THIS IS WHERE THE SETTING stretchFactor IS USED.
meshL Th = buildmeshL(circle(n));

fespace Uh(Th,P1);

BemKernel ker1("SL",k=k);
varf vk1(u,v)=int1dx1d(Th)(Th)(BEM(ker1,u,v)) ;  
HMatrix<complex> HFirstKind = vk1(Uh,Uh);

Uh<complex> uFirstKind, bFirstKind;
varf vmassFirstKind(u,v) = int1d(Th)((bc)*v);
bFirstKind[] = vmassFirstKind(0,Uh);
uFirstKind[] = HFirstKind^-1*bFirstKind[];

mesh ThOut = buildmesh(circle(200));
fespace UhOut(ThOut,P1);

BemPotential Pot("SL",k=k);
varf vp(u,v)=int1d(Th)(POT(Pot,u,v)) ;  

HMatrix<complex> HPot = vp(Uh,UhOut,eta=10,eps=1e-3,minclustersize=10,maxblocksize=1000000);

UhOut<complex> vFirstKind;
vFirstKind[] = HPot*uFirstKind[];
UhOut vFirstKindReal = real(vFirstKind);

plot(vFirstKindReal, dim=1, fill=1,value=1, nbiso=50, wait=true);
10

Here is the second program. In this case there is nothing to uncomment.

load "bem";
load "msh3";
 
int k = 0;
real epsilon = 1e-10;

func real findt(real xx, real yy, real X1, real Y1, real X2,real Y2){
real t;
if(abs(X1-X2) < epsilon){
t = (yy-Y1)/(Y2-Y1);
return t;
}
else{
t = (xx-X1)/(X2-X1);
return t;
}
};

border Gamma1(t=0., 1){x=0.38209*(1-t) + 0.42502*(t); y=0.23174*(1-t) + 0.20594*(t); region = 1;}
border Gamma2(t=0., 1){x=0.42502*(1-t) + 0.43293*(t); y=0.20594*(1-t) + 0.22413*(t); region = 2;}
border Gamma3(t=0., 1){x=0.43293*(1-t) + 0.43395*(t); y=0.22413*(1-t) + 0.23477*(t); region = 3;}
border Gamma4(t=0., 1){x=0.43395*(1-t) + 0.46091*(t); y=0.23477*(1-t) + 0.24775*(t); region = 4;}
border Gamma5(t=0., 1){x=0.46091*(1-t) + 0.45655*(t); y=0.24775*(1-t) + 0.26457*(t); region = 5;}
border Gamma6(t=0., 1){x=0.45655*(1-t) + 0.4702*(t); y=0.26457*(1-t) + 0.25033*(t); region = 6;}
border Gamma7(t=0., 1){x=0.4702*(1-t) + 0.48683*(t); y=0.25033*(1-t) + 0.22076*(t); region = 7;}
border Gamma8(t=0., 1){x=0.48683*(1-t) + 0.52613*(t); y=0.22076*(1-t) + 0.25234*(t); region = 8;}
border Gamma9(t=0., 1){x=0.52613*(1-t) + 0.48608*(t); y=0.25234*(1-t) + 0.25757*(t); region = 9;}
border Gamma10(t=0., 1){x=0.48608*(1-t) + 0.49534*(t); y=0.25757*(1-t) + 0.27633*(t); region = 10;}
border Gamma11(t=0., 1){x=0.49534*(1-t) + 0.49881*(t); y=0.27633*(1-t) + 0.28588*(t); region = 11;}
border Gamma12(t=0., 1){x=0.49881*(1-t) + 0.53623*(t); y=0.28588*(1-t) + 0.30591*(t); region = 12;}
border Gamma13(t=0., 1){x=0.53623*(1-t) + 0.49897*(t); y=0.30591*(1-t) + 0.29396*(t); region = 13;}
border Gamma14(t=0., 1){x=0.49897*(1-t) + 0.47359*(t); y=0.29396*(1-t) + 0.27985*(t); region = 14;}
border Gamma15(t=0., 1){x=0.47359*(1-t) + 0.4471*(t); y=0.27985*(1-t) + 0.28364*(t); region = 15;}
border Gamma16(t=0., 1){x=0.4471*(1-t) + 0.44953*(t); y=0.28364*(1-t) + 0.27544*(t); region = 16;}
border Gamma17(t=0., 1){x=0.44953*(1-t) + 0.41959*(t); y=0.27544*(1-t) + 0.28921*(t); region = 17;}
border Gamma18(t=0., 1){x=0.41959*(1-t) + 0.39141*(t); y=0.28921*(1-t) + 0.28922*(t); region = 18;}
border Gamma19(t=0., 1){x=0.39141*(1-t) + 0.38209*(t); y=0.28922*(1-t) + 0.23174*(t); region = 19;}

func real BC1(real xx, real yy){
real T;
T = findt(xx, yy,0.38209,0.23174,0.42502,0.20594);
real ans;
ans = -0.61135 + -0.27391*sin(( pi/2 +  pi*(1- 1))*T) + 0.24543*sin(( pi/2 +  pi*(2- 1))*T) + 0.96603*sin(( pi/2 +  pi*(3- 1))*T) + -0.63572*sin(( pi/2 +  pi*(4- 1))*T) + -0.67363*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC2(real xx, real yy){
real T;
T = findt(xx, yy,0.42502,0.20594,0.43293,0.22413);
real ans;
ans = -0.20258 + -0.075077*sin(( pi/2 +  pi*(1- 1))*T) + -0.83038*sin(( pi/2 +  pi*(2- 1))*T) + 0.80321*sin(( pi/2 +  pi*(3- 1))*T) + 0.56767*sin(( pi/2 +  pi*(4- 1))*T) + -0.94681*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC3(real xx, real yy){
real T;
T = findt(xx, yy,0.43293,0.22413,0.43395,0.23477);
real ans;
ans = -0.15854 + -0.37852*sin(( pi/2 +  pi*(1- 1))*T) + -0.071237*sin(( pi/2 +  pi*(2- 1))*T) + 0.23449*sin(( pi/2 +  pi*(3- 1))*T) + -0.66526*sin(( pi/2 +  pi*(4- 1))*T) + -0.97925*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC4(real xx, real yy){
real T;
T = findt(xx, yy,0.43395,0.23477,0.46091,0.24775);
real ans;
ans = -0.54532 + -0.46983*sin(( pi/2 +  pi*(1- 1))*T) + -0.098644*sin(( pi/2 +  pi*(2- 1))*T) + -0.87349*sin(( pi/2 +  pi*(3- 1))*T) + -0.27098*sin(( pi/2 +  pi*(4- 1))*T) + 0.39175*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC5(real xx, real yy){
real T;
T = findt(xx, yy,0.46091,0.24775,0.45655,0.26457);
real ans;
ans = -1.1273 + -0.635*sin(( pi/2 +  pi*(1- 1))*T) + -0.61341*sin(( pi/2 +  pi*(2- 1))*T) + -0.95555*sin(( pi/2 +  pi*(3- 1))*T) + 0.17258*sin(( pi/2 +  pi*(4- 1))*T) + -0.72678*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC6(real xx, real yy){
real T;
T = findt(xx, yy,0.45655,0.26457,0.4702,0.25033);
real ans;
ans = -3.0038 + 0.28593*sin(( pi/2 +  pi*(1- 1))*T) + -0.47272*sin(( pi/2 +  pi*(2- 1))*T) + 0.69709*sin(( pi/2 +  pi*(3- 1))*T) + 0.070972*sin(( pi/2 +  pi*(4- 1))*T) + -0.95075*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC7(real xx, real yy){
real T;
T = findt(xx, yy,0.4702,0.25033,0.48683,0.22076);
real ans;
ans = -2.5697 + 0.26014*sin(( pi/2 +  pi*(1- 1))*T) + 0.48737*sin(( pi/2 +  pi*(2- 1))*T) + 0.52573*sin(( pi/2 +  pi*(3- 1))*T) + -0.30955*sin(( pi/2 +  pi*(4- 1))*T) + -0.13706*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC8(real xx, real yy){
real T;
T = findt(xx, yy,0.48683,0.22076,0.52613,0.25234);
real ans;
ans = -2.0987 + 0.16838*sin(( pi/2 +  pi*(1- 1))*T) + -0.54246*sin(( pi/2 +  pi*(2- 1))*T) + 0.24401*sin(( pi/2 +  pi*(3- 1))*T) + -0.93799*sin(( pi/2 +  pi*(4- 1))*T) + 0.27936*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC9(real xx, real yy){
real T;
T = findt(xx, yy,0.52613,0.25234,0.48608,0.25757);
real ans;
ans = 0.073463 + 0.54455*sin(( pi/2 +  pi*(1- 1))*T) + -0.12756*sin(( pi/2 +  pi*(2- 1))*T) + 0.34777*sin(( pi/2 +  pi*(3- 1))*T) + -0.51289*sin(( pi/2 +  pi*(4- 1))*T) + -0.24397*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC10(real xx, real yy){
real T;
T = findt(xx, yy,0.48608,0.25757,0.49534,0.27633);
real ans;
ans = 1.3623 + -0.54388*sin(( pi/2 +  pi*(1- 1))*T) + 0.40629*sin(( pi/2 +  pi*(2- 1))*T) + 0.56381*sin(( pi/2 +  pi*(3- 1))*T) + 0.23853*sin(( pi/2 +  pi*(4- 1))*T) + -0.75657*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC11(real xx, real yy){
real T;
T = findt(xx, yy,0.49534,0.27633,0.49881,0.28588);
real ans;
ans = -0.019211 + 0.75053*sin(( pi/2 +  pi*(1- 1))*T) + -0.92659*sin(( pi/2 +  pi*(2- 1))*T) + 0.12716*sin(( pi/2 +  pi*(3- 1))*T) + -0.4567*sin(( pi/2 +  pi*(4- 1))*T) + 0.97077*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC12(real xx, real yy){
real T;
T = findt(xx, yy,0.49881,0.28588,0.53623,0.30591);
real ans;
ans = 3.2125 + 0.94017*sin(( pi/2 +  pi*(1- 1))*T) + 0.93232*sin(( pi/2 +  pi*(2- 1))*T) + 0.70643*sin(( pi/2 +  pi*(3- 1))*T) + -0.47362*sin(( pi/2 +  pi*(4- 1))*T) + -0.98286*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC13(real xx, real yy){
real T;
T = findt(xx, yy,0.53623,0.30591,0.49897,0.29396);
real ans;
ans = 3.4176 + -0.64805*sin(( pi/2 +  pi*(1- 1))*T) + 0.59052*sin(( pi/2 +  pi*(2- 1))*T) + 0.72157*sin(( pi/2 +  pi*(3- 1))*T) + 0.77674*sin(( pi/2 +  pi*(4- 1))*T) + -0.12991*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC14(real xx, real yy){
real T;
T = findt(xx, yy,0.49897,0.29396,0.47359,0.27985);
real ans;
ans = 1.9939 + 0.87496*sin(( pi/2 +  pi*(1- 1))*T) + 0.87264*sin(( pi/2 +  pi*(2- 1))*T) + 0.88164*sin(( pi/2 +  pi*(3- 1))*T) + -0.58886*sin(( pi/2 +  pi*(4- 1))*T) + 0.79885*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC15(real xx, real yy){
real T;
T = findt(xx, yy,0.47359,0.27985,0.4471,0.28364);
real ans;
ans = 4.2656 + -0.5964*sin(( pi/2 +  pi*(1- 1))*T) + 0.45966*sin(( pi/2 +  pi*(2- 1))*T) + 0.33778*sin(( pi/2 +  pi*(3- 1))*T) + 0.4788*sin(( pi/2 +  pi*(4- 1))*T) + -0.73609*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC16(real xx, real yy){
real T;
T = findt(xx, yy,0.4471,0.28364,0.44953,0.27544);
real ans;
ans = 2.3324 + 0.45231*sin(( pi/2 +  pi*(1- 1))*T) + -0.92794*sin(( pi/2 +  pi*(2- 1))*T) + -0.53618*sin(( pi/2 +  pi*(3- 1))*T) + 0.22167*sin(( pi/2 +  pi*(4- 1))*T) + 0.74808*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC17(real xx, real yy){
real T;
T = findt(xx, yy,0.44953,0.27544,0.41959,0.28921);
real ans;
ans = 3.7029 + -0.85721*sin(( pi/2 +  pi*(1- 1))*T) + -0.84537*sin(( pi/2 +  pi*(2- 1))*T) + -0.41625*sin(( pi/2 +  pi*(3- 1))*T) + 0.55012*sin(( pi/2 +  pi*(4- 1))*T) + -0.22203*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC18(real xx, real yy){
real T;
T = findt(xx, yy,0.41959,0.28921,0.39141,0.28922);
real ans;
ans = 2.5027 + 0.47231*sin(( pi/2 +  pi*(1- 1))*T) + -0.98573*sin(( pi/2 +  pi*(2- 1))*T) + 0.40533*sin(( pi/2 +  pi*(3- 1))*T) + -0.29559*sin(( pi/2 +  pi*(4- 1))*T) + -0.30102*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

func real BC19(real xx, real yy){
real T;
T = findt(xx, yy,0.39141,0.28922,0.38209,0.23174);
real ans;
ans = 4.3606 + -1.8784*sin(( pi/2 +  pi*(1- 1))*T) + 0.84928*sin(( pi/2 +  pi*(2- 1))*T) + -0.50152*sin(( pi/2 +  pi*(3- 1))*T) + 0.93544*sin(( pi/2 +  pi*(4- 1))*T) + -0.80731*sin(( pi/2 +  pi*(5- 1))*T);
return ans;
 };

int numpt = 60;

meshL Th = buildmeshL(Gamma1(numpt) + Gamma2(numpt) + Gamma3(numpt) + Gamma4(numpt) + Gamma5(numpt) + Gamma6(numpt) + Gamma7(numpt) + Gamma8(numpt) + Gamma9(numpt) + Gamma10(numpt) + Gamma11(numpt) + Gamma12(numpt) + Gamma13(numpt) + Gamma14(numpt) + Gamma15(numpt) + Gamma16(numpt) + Gamma17(numpt) + Gamma18(numpt) + Gamma19(numpt));

mesh ThOut = buildmesh(Gamma1(numpt*2) + Gamma2(numpt*2) + Gamma3(numpt*2) + Gamma4(numpt*2) + Gamma5(numpt*2) + Gamma6(numpt*2) + Gamma7(numpt*2) + Gamma8(numpt*2) + Gamma9(numpt*2) + Gamma10(numpt*2) + Gamma11(numpt*2) + Gamma12(numpt*2) + Gamma13(numpt*2) + Gamma14(numpt*2) + Gamma15(numpt*2) + Gamma16(numpt*2) + Gamma17(numpt*2) + Gamma18(numpt*2) + Gamma19(numpt*2));

fespace Uh(Th,P1);
fespace UhOut(ThOut,P1);

BemKernel ker1("SL",k=k);
varf vk1(u,v)=int1dx1d(Th)(Th)(BEM(ker1,u,v)) ;  
HMatrix<complex> HFirstKind = vk1(Uh,Uh,eta=10,eps=1e-3,minclustersize=10,maxblocksize=1000000);

Uh<complex> uFirstKind, bFirstKind;
varf vmassFirstKind(u,v) = int1d(Th,1)((BC1(x,y))*v)+ int1d(Th,2)((BC2(x,y))*v)+ int1d(Th,3)((BC3(x,y))*v)+ int1d(Th,4)((BC4(x,y))*v)+ int1d(Th,5)((BC5(x,y))*v)+ int1d(Th,6)((BC6(x,y))*v)+ int1d(Th,7)((BC7(x,y))*v)+ int1d(Th,8)((BC8(x,y))*v)+ int1d(Th,9)((BC9(x,y))*v)+ int1d(Th,10)((BC10(x,y))*v)+ int1d(Th,11)((BC11(x,y))*v)+ int1d(Th,12)((BC12(x,y))*v)+ int1d(Th,13)((BC13(x,y))*v)+ int1d(Th,14)((BC14(x,y))*v)+ int1d(Th,15)((BC15(x,y))*v)+ int1d(Th,16)((BC16(x,y))*v)+ int1d(Th,17)((BC17(x,y))*v)+ int1d(Th,18)((BC18(x,y))*v)+ int1d(Th,19)((BC19(x,y))*v);
bFirstKind[] = vmassFirstKind(0,Uh);

uFirstKind[] = HFirstKind^-1*bFirstKind[];
 
varf vp(u,v) = int1d(Th)(POT(BemPotential("SL",k=k),u,v));

HMatrix<complex> HPot = vp(Uh,UhOut);
UhOut<complex> vFirstKind;
vFirstKind[] = HPot*uFirstKind[];
UhOut vFirstKindReal = real(vFirstKind);
real maxsol = vFirstKindReal[].max;
real minsol = vFirstKindReal[].min;

real[int] viso(20);
for (int i = 0; i < viso.n; i++)
 viso[i] = minsol + i*(maxsol-minsol)/viso.n;

real[int] colorhsv=[
0, 0 , 0,
1, 0. , 1
];

mesh OptMesh = adaptmesh(ThOut,vFirstKindReal);
fespace OptUh(OptMesh,P1);
matrix I = interpolate(OptUh, UhOut);
OptUh OptSol;
OptSol[] = I*vFirstKindReal[];

plot(OptMesh, wait=true, cmm="Optimized Mesh: Notice the pockmarks.");
plot(OptSol, dim=1, viso=viso(0:viso.n), fill=true, wait=true, hsv=colorhsv, cmm="Activate the mesh to see the location of spurious extrema (dense meshing away from boundary).");
11

It appears that in both cases the problem comes from the compression error of the low-rank compression algorithm when assembling the H-Matrix for the potential. If you set eps (the target compression error) to be tighter, for example eps=1e-5 on line 45 of the first code

HMatrix<complex> HPot = vp(Uh,UhOut,eta=10,eps=1e-5,minclustersize=10,maxblocksize=1000000);

and same on line 212 of the second code

HMatrix<complex> HPot = vp(Uh,UhOut,eps=1e-5);

it appears to solve both problems. Also keep in mind that for validation/debugging purposes, you can set the geometric admissibility criterion eta to a negative value, i.e. eta = -1, to disable low-rank compression (and thus remove any compression error that might pollute your solution).

Another remark is that if the BEM boundary is included in the output mesh used for the reconstruction of the solution, it may lead to singularities or near-singularities in the Potential, altering the quality of your reconstructed solution. For example, you will notice that you will get a more regular solution near the boundary if you define your output mesh as

border circle2(t=0, 2*pi){x=0.99*stretchFactor*cos(t); y=0.99*sin(t);} 
mesh ThOut = buildmesh(circle2(200));
12

Setting eps=1e-5 and eta=-1 improves the results, but does not solve the whole issue. I still get qualitatively different results for the two different values of stretchFactor in the first program. However, the strange-looking solution is now smooth.

Regarding your comments on the functions near the boundary, how would you proceed to evaluate and illustrate the solutions, especially in the case of complicated domains such as the one used in my second program? I was hoping that a projection onto a fine finite element mesh would do the trick, but it seems that it is not so.

13

For the Laplace equation in 2D, the well-posedness of the problem depends on the diameter of the geometry at hand: from

@Book{Steinbach2008NAM,
author = {Steinbach, Olaf},
publisher = {Springer, New York},
title = {Numerical approximation methods for elliptic boundary value problems},
year = {2008},
doi = {10.1007/978-0-387-68805-3}
}

Thm

Indeed, the problem disappears if you consider a circle of radius 0.9 instead of 1.
For a more general solution, I am afraid that for now you would have to change the value of the relevant constant in the C++ code of the Laplace Potential in BEMTool and recompile the plugin, which is not super convenient.

For your second comment, I am afraid there is no easy solution for now. Using any non-matching background grid as an output mesh usually does the trick, but if you are unlucky some vertex could be too close to a boundary vertex which could generate instabilities in the reconstructed solution near the boundary.