function [p,pneg,ppos] = distpremesh(fd,fh,bpm,varargin)

geps = bpm.mesh_tol * bpm.hmin; %.001*bpm.hmin; 

% 1. Create initial distribution in bounding box (equilateral triangles)
[x,y] = meshgrid(bpm.bbox(1,1):bpm.hmin:bpm.bbox(2,1), bpm.bbox(1,2):bpm.hmin*sqrt(3)/2:bpm.bbox(2,2));
x(2:2:end,:) = x(2:2:end,:) + bpm.hmin/2;  % Shift even rows
p = [x(:),y(:)];                           % List of node coordinates

%pneg = p(feval(fd,p)<geps, :);  % keep these points
ppos = p( feval(fd,p) >= geps, :); % neglect these points

% 2. Remove points outside the region, apply the rejection method
p = p(feval(fd,p,varargin{:})<geps,:);                 % Keep only d<0 points
r0 = 1./feval(fh,p,varargin{:}).^2;                    % Probability to keep point: randomly get rid of points whose metric is larger than normalized metric of sizing function 
p = p(rand(size(p,1),1) < r0./max(r0), :);             % Rejection method

% fh
% maxr0 = max(r0)
% minr0 = min(r0)
% figure(11), clf, 
% hold on
% plot(r0,'b')
% plot(p(:,1),'r')
% plot(p(:,2),'g')

% rand(size(p,1),1)

if ~isempty(bpm.pfix)                   % Remove duplicated nodes
    p = setdiff(p,bpm.pfix,'rows'); 
end 

pfix = unique(bpm.pfix,'rows'); 
p = [pfix; p];                                         % Prepend fix points

pneg = p;   % keep these points
