% examples of mesh generation. IMPORTANT: there are multiple solutions! You are welcome to try other ways of
% solving the exercises.

%% hole
%-------------------------
clear
x_c = 0.6;
y_c = 0.4;
rad_c = 0.1;

% way 1: element-wise multiplication of d in each region: d1<0 inside the square, d2<0 inside the circle, so
% inside the circle d1*d2>0 -> no mesh as desired
square = @(P) max(abs( P(:,1)-0.5), abs(P(:,2) - 0.5)) -0.5 ;
circle = @(P) sqrt((P(:,1)-x_c).^2 + (P(:,2)-y_c).^2) - rad_c ;
hole_d = @(P)  circle(P) .* square(P);
[P,T]=distmesh2d(hole_d,@huniform,0.05,[-1 -1; 2 2;],[0 0; 1 1; 0 1; 1 0]);

% fixing just two points out of four will make the convergence of distmesh more difficult ("the mesh shakes").
% Try the following command, 
% [P,T]=distmesh2d(hole_d,@huniform,0.05,[-1 -1; 2 2;],[0 0; 1 1;]);


% way 2: use max: d1<0 inside the square, d2<0 inside the circle, so -d2<0 *outside* the circle. Then 
% -- inside  the circle max(d1,-d2) = -d2 >0 
% -- outside the circle max(d1,-d2) = max(two negative numbers) < 0


hole_d2 = @(P) max(square(P),-circle(P));
figure
[P,T]=distmesh2d(hole_d2,@huniform,0.05,[-1 -1; 2 2;],[0 0; 1 1; 0 1; 1 0]);

% way 3: use distmesh commands, such as dcircle, drectangle, ddiff
hole_d3 = @(P) ddiff(drectangle(P,0,1,0,1),dcircle(P,x_c,y_c,rad_c));
figure
[P,T]=distmesh2d(hole_d3,@huniform,0.05,[-1 -1; 2 2;],[0 0; 1 1; 0 1; 1 0]);



%% half circle
%-----------------------
clear
x_c = 0.0;
y_c = 0.0;
rad_c = 1;
halfc_d = @(P) (-P(:,1)) .* ( sqrt((P(:,1)-x_c).^2 + (P(:,2)-y_c).^2) - rad_c  );

[P,T]=distmesh2d(halfc_d,@huniform,0.05,[-1 -1; 0 1;],[0 -1; 0 1]);


%% moon
%-----------------------
clear
x_c = 0.0;
y_c = 0.0;
rad_c = 1;
x_c2 = 1;
y_c2 = 0.0;
rad_c2 = sqrt(rad_c^2 + x_c2^2);

moon_d = @(P)  ( sqrt((P(:,1)-x_c).^2 + (P(:,2)-y_c).^2) - rad_c  ) .* (-rad_c2^2 + ((P(:,1)-x_c2).^2 + (P(:,2)-y_c2).^2) );

% this one won't work (although formally correct), three constraints are difficult to handle
%moon_d = @(P)  (-P(:,1)).*( sqrt((P(:,1)-x_c).^2 + (P(:,2)-y_c).^2) - rad_c  ) .* (-rad_c2^2 + ((P(:,1)-x_c2).^2 + (P(:,2)-y_c2).^2) );

figure
[P,T]=distmesh2d(moon_d,@(P) ones(size(P,1),1) ,0.05,[-1 -1; 0 1],[0 -1; 0 1]);

%% non uniform meshes: 
% ---------------------------

% the second function gives a number proportional to hsize, so we want it to be larger and larger when x gets
% far from 0
figure
[P,T]=distmesh2d(moon_d,@(P) exp(-1.4*P(:,1)) ,0.02,[-1 -1; 0 1],[0 -1; 0 1]);

figure
[P,T]=distmesh2d(moon_d,@(P) 1./(cos(P(:,1)).^2) ,0.02,[-1 -1; 0 1],[0 -1; 0 1]);

% the first one is more effective in refining "only at x=0"

%% recover border nodes

[P,T]=distmesh2d(moon_d,@(P) ones(size(P,1),1) ,0.05,[-1 -1; 0 1],[0 -1; 0 1]);

%
val = moon_d(P);

% this won't work, due to machine precision issues. Also, distmesh proceeds iteratively, and will stop 
% as soon as it has found the position of the border d==0 "up to a certain tolerance"

B = P(val==0,:); 

B = P(abs(val)<1e-4,:); 

hold on
plot(B(:,1),B(:,2),'or','MarkerFaceColor','r')

%% horn --> problem with the tip
%-----------------------

% if the two circles have equal tangent in x=-1, then "no triangle fits in that corner" (cuspid) and the
% mesh generator cannot catch it properly. To see this, run the code below, then change y_c2=0, xmin_c2=-1 and
% run again.

clear
x_c = 0.0;
y_c = 0.0;
rad_c = 1;
x_c2 = -0.6; % the x coord of the center of a second circle. Removing this one from the previous one we get the "horn" shape
y_c2 = -0.2; % the y coord of the center of a second circle. Removing this one from the previous one we get the "horn" shape
xmin_c2 = -0.8; % the leftmost point where the second circle instersects the horizontal axis. Note that x_c2-xmin_c2 is
                % in general *not* equal to the radius of the second circle. It is so only if y_c2=0
rad_c2 =sqrt((x_c2-xmin_c2)^2 + y_c2^2);

twoc = @(P) ( sqrt((P(:,1)-x_c).^2 + (P(:,2)-y_c).^2) - rad_c  ) .* ( sqrt((P(:,1)-x_c2).^2 + (P(:,2)-y_c2).^2) - rad_c2 );

% again, three constraints are tough, so I use a combination of * and max here
horn_d =@(P) max( twoc(P),-P(:,2) ); 

[P,T]=distmesh2d(horn_d,@huniform,0.04,[-1 -1; 1 1;],[-1 0; 1 0; 0 1; xmin_c2 0; x_c2+rad_c2 0]);
