This is the code I used to create my tear drop region.

compute bdry info for region
function [y,dy,ldy,ny,k,sh]=reg(t);
%t is our set of sampled values in[-1 1]
%set up parameters
sh=pi/10; %amount we shift the region
t=t(:);
%set up region
r=cosh(pi*t);
dr=pi*sinh(pi*t);
ddr=pi^2*cosh(pi*t);
co=cos(pi*t+sh);%angles go from [-pi+sh,pi+sh]
si=sin(pi*t+sh);
y=[r.*co,r.*si]; %gives pts on bdry of region in cartesian coord
dy=[dr.*co-pi*r.*si,dr.*si+pi*r.*co]; %gives dy/dt
ddy=[ddr.*co-2*pi*dr.*si-pi^2*r.*co,ddr.*si+2*pi*dr.*co-pi^2*r.*si]; %gives d^2y/dt^2
ldy=sqrt(dy(:,1).^2+dy(:,2).^2); %gives lenght of dy
ny=[dy(:,2)./ldy,-dy(:,1)./ldy];%gives normal derivative at each pt on bdry
k=curv(ny,ldy,ddy); %gives 1/(4pi)*curvature at each sample pt

This next computes the dimpled region with corner at end pts.

function [y,dy,ldy,ny,k,sh]=Lbut(t);
%t is our set of sampled values in[-1 1]

%set up parameters
t=t(:);
sh=11*pi/10; %amount we shift the region
r=1+2*abs(sin(pi*t/2+pi/2));%sets up radius
dr=zeros(size(t));
ddr=zeros(size(t));
N=length(t);
for n=1:N;
if sin(pi*t(n)/2+pi/2)>0
dr(n)=pi*cos(pi*t(n)/2+pi/2);%gives first der of bdry
ddr(n)=-pi^2/2*sin(pi*t(n)/2+pi/2);%gives second der of bdry
else
dr(n)=-pi*cos(pi*t(n)/2);
ddr(n)=pi^2/2*sin(pi*t(n)/2);
end
end
co=cos(pi*t+sh);%angles go from [-pi+sh,pi+sh]
si=sin(pi*t+sh);
y=[r.*co,r.*si]; %gives pts on bdry of region in cartesian coord
dy=[dr.*co-pi*r.*si,dr.*si+pi*r.*co]; %gives dy/dt
ddy=[ddr.*co-pi*2*dr.*si-pi^2*r.*co,ddr.*si+pi*2*dr.*co-pi^2*r.*si]; %gives d^2y/dt^2
ldy=sqrt(dy(:,1).^2+dy(:,2).^2); %gives lenght of dy
ny=[dy(:,2)./ldy,-dy(:,1)./ldy];%gives normal derivative at each pt on bdry
k=curv(ny,ldy,ddy); %gives 1/(4pi)*curvature at each sample pt

The next code computes the dimpled region with corner in middle

function [y,dy,ldy,ny,k,sh]=dimp(t);
%t is our set of sampled values in[-1 1]

%set up parameters
t=t(:);
sh=pi/10; %amount we shift the region
r=1+2*abs(sin(pi*t/2));%sets up radius
dr=zeros(size(t));
ddr=zeros(size(t));
N=length(t);
for n=1:N;
if pi*t(n)/2>0
dr(n)=pi*cos(pi*t(n)/2);%gives first der of bdry
ddr(n)=-pi^2/2*sin(pi*t(n)/2);%gives second der of bdry
else
dr(n)=-pi*cos(pi*t(n)/2);
ddr(n)=pi^2/2*sin(pi*t(n)/2);
end
end
co=cos(pi*t+sh);%angles go from [-pi+sh,pi+sh]
si=sin(pi*t+sh);
y=[r.*co,r.*si]; %gives pts on bdry of region in cartesian coord
dy=[dr.*co-pi*r.*si,dr.*si+pi*r.*co]; %gives dy/dt
ddy=[ddr.*co-pi*2*dr.*si-pi^2*r.*co,ddr.*si+pi*2*dr.*co-pi^2*r.*si]; %gives d^2y/dt^2
ldy=sqrt(dy(:,1).^2+dy(:,2).^2); %gives lenght of dy
ny=[dy(:,2)./ldy,-dy(:,1)./ldy];%gives normal derivative at each pt on bdry
k=curv(ny,ldy,ddy); %gives 1/(4pi)*curvature at each sample pt

The next code is due to [Tref 3]1
and is used to calculate gauss quad points and weights

%gauss nodes x (Legendre points) and weights w
% for gausss quadrature

function [x,w]=gauss(N)
beta=.5./sqrt(1-(2*(1:N-1)).^(-2));
T=diag(beta,1)+diag(beta,-1);
[V,D]=eig(T);
x=diag(D); [x,i]=sort(x);
w=2*V(1,i).^2;

The next code is also due to [Tref 3]1
and is used to calculate Clenshaw-Curtis quad points and weights

%clencurt nodes x (Chebyshev points) and weights w
%For clenshaw curtis quadrature

function [x,w]=clencurt(N)
theta=pi*(0:N)'/N; x=cos(theta);
w=zeros(1,N+1); ii=2:N; v=ones(N-1,1);
if mod(N,2)==0
w(1)=1/(N^2-1);w(N+1)=w(1);
for k=1:N/2-1,v=v-2*cos(2*k*theta(ii))/(4*k^2-1); end
v=v-cos(N*theta(ii))/(N^2-1);;
else
w(1)=1/N^2; w(N+1)=w(1);
for k=1:(N-1)/2, v=v-2*cos(2*k*theta(ii))/(4*k^2-1); end
end
w(ii)=2*v/N;

The next code was used to calculate the intergral rep using the Nystrom method
given a set of quad points and regions. You can change quad points by changing what is commented out and you
can change regions by changing the code that is called

function [tau,y,ny,W]=nystint(N)

[x,w]=gauss(N);%calculates Legendre pts and wghts for gauss quadrature
%x=[-1:2/N:1-2/N]';%calculates even wght pts
%w=2/N*ones(size(x'));
%[x,w]=clencurt(N);%calc clen curt pts and wghts
%x=x([1:N]);w=w([1:N]);
[y,dy,ldy,ny,k,sh]=Lbut(x);%sets up info for region at pts x can change region here
theta=pi*x+sh;%gives angles associated to x
%figure
%plot(y(:,1),y(:,2),'.','markersize',16)
%line(y(:,1),y(:,2))
%title('region')

%set up boundary function
har=y(:,1).*y(:,2);
%set up nyst method for interior problem
K=zeros(N,N); %gives size of nyst ker
for j=1:N
K(:,j)=dipole(y(:,1),y(:,2),y(j,1),y(j,2),ny(j,1),ny(j,2));%gives dipole kernel at given pts
end
for j=1:N
K(j,j)=k(j); %gives diag of dipole
end
W=w.*ldy';
L=repmat(W,N,1);%sets up weights
Ker=L.*K; %gives nyst kern
tau=(Ker-.5*eye(N))\har;%gives representation of h=x1x2

The next function will plot the solution to interior dirchlet problem for the region and quaduture specified in nystint.m

function dirintpl(N); %plots the internal dirchlet problem
[tau,y,ny,W]=nystint(N);%gives int rep tau, interp points y, normal derivatives ny and weights W

%est integral
g1=[-3:.1:2]';%setting up grid pts should change them to suit the region you are working with
g2=[-3:.1:3]';
nu1=length(g1);nu2=length(g2);
[x1,x2]=meshgrid(g1,g2);
X1=reshape(x1,nu1*nu2,1);
X2=reshape(x2,nu1*nu2,1);
Z=estu(X1,X2,y(:,1),y(:,2),ny(:,1),ny(:,2),tau,W'); %estimates h in region

%plot

Z=reshape(Z,nu2,nu1);
figure
mesh(x1,x2,Z)
view(90,90)

figure
contour(x1,x2,Z,[-4:.2:4])

This function will compute the value of the est solution for the internal dirchlet problem for the region and quaduture specified in nystint.m and point (pi,p2)

function Z=convdir(N,m) %lists convergence at point p1, p2 of inter dir up to N sample pts at intervals m
%should have m|N
p1=-1%gives x coord of test point
p2=-1%gives y coord of test point
Z=zeros(N/m-1,1);
for n=m:m:N;
[tau,y,ny,W]=nystint(n);%gives int rep tau, interp points y, normal derivatives ny and weights W
Z(n/m)=estu(p1,p2,y(:,1),y(:,2),ny(:,1),ny(:,2),tau,W');
end