Matlab椭圆网格划分

clc;clear all;close all;
% Example: (Ellipse)
    fd=@(p) p(:,1).^2/2^2+p(:,2).^2/1^2-1;
    [p,t]=distmesh2d(fd,@huniform,0.2,[-2,-1;2,1],[]);
   figure(2);[X,Y]=meshgrid(p(:,1),p(:,2));
   plot(p(:,1),p(:,2),'ro');  

function d=drectangle(p,x1,x2,y1,y2)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

d=-min(min(min(-y1+p(:,2),y2-p(:,2)),-x1+p(:,1)),x2-p(:,1));
end
    
function d=dcircle(p,xc,yc,r)
%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.
d=sqrt((p(:,1)-xc).^2+(p(:,2)-yc).^2)-r;
end

function d=ddiff(d1,d2), d=max(d1,-d2);
%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.
end
%cc=1;

function d=dpoly(p,pv)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

np=size(p,1);
nvs=size(pv,1)-1;

ds=dsegment(p,pv);
%ds=zeros(np,nvs);
%for iv=1:nvs
%  ds(:,iv)=donesegment(p,pv(iv:iv+1,:));
%end
d=min(ds,[],2);
d=(-1).^(inpolygon(p(:,1),p(:,2),pv(:,1),pv(:,2))).*d;
% MEXED
%function ds=donesegment(p,pv)
%
%e=ones(size(p,1),1);
%
%v=diff(pv,1);
%w=p-e*pv(1,:);
%
%c1=sum(w.*v(e,:),2);
%c2=sum(v(e,:).^2,2);
%
%ds=0*e;
%
%ix=c1<=0;
%ds(ix)=sqrt(sum((p(ix,:)-pv(1*ones(sum(ix),1),:)).^2,2));
%
%ix=c1>=c2;
%ds(ix)=sqrt(sum((p(ix,:)-pv(2*ones(sum(ix),1),:)).^2,2));
%
%ix=c1>0 & c2>c1;
%nix=sum(ix);
%if nix>0
%  Pb=ones(nix,1)*pv(1,:)+c1(ix)./c2(ix)*v;
%  ds(ix)=sqrt(sum((p(ix,:)-Pb).^2,2));
%end
end

    function tri=surftri(p,t)
%SURFTRI Find surface triangles from tetrahedra mesh
%   TRI=SURFTRI(P,T)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

% Form all faces, non-duplicates are surface triangles
faces=[t(:,[1,2,3]);
       t(:,[1,2,4]);
       t(:,[1,3,4]);
       t(:,[2,3,4])];
node4=[t(:,4);t(:,3);t(:,2);t(:,1)];
faces=sort(faces,2);
[foo,ix,jx]=unique(faces,'rows');
vec=histc(jx,1:max(jx));
qx=find(vec==1);
tri=faces(ix(qx),:);
node4=node4(ix(qx));

% Orientation
v1=p(tri(:,2),:)-p(tri(:,1),:);
v2=p(tri(:,3),:)-p(tri(:,1),:);
v3=p(node4,:)-p(tri(:,1),:);
ix=find(dot(cross(v1,v2,2),v3,2)>0);
tri(ix,[2,3])=tri(ix,[3,2]);
end

function v=simpvol(p,t)
%SIMPVOL Simplex volume.
%   V=SIMPVOL(P,T)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

switch size(p,2)
 case 1
  d12=p(t(:,2),:)-p(t(:,1),:);
  v=d12;
 case 2
  d12=p(t(:,2),:)-p(t(:,1),:);
  d13=p(t(:,3),:)-p(t(:,1),:);
  v=(d12(:,1).*d13(:,2)-d12(:,2).*d13(:,1))/2;
 case 3
  d12=p(t(:,2),:)-p(t(:,1),:);
  d13=p(t(:,3),:)-p(t(:,1),:);
  d14=p(t(:,4),:)-p(t(:,1),:);
  v=dot(cross(d12,d13,2),d14,2)/6;
 otherwise
  v=zeros(size(t,1),1);
  for ii=1:size(t,1)
    A=zeros(size(p,2)+1);
    A(:,1)=1;
    for jj=1:size(p,2)+1
      A(jj,2:end)=p(t(ii,jj),:);
    end
    v(ii)=det(A);
  end
  v=v/factorial(size(p,2));
end

end

function h=huniform(p,varargin)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

h=ones(size(p,1),1);
end

function [p,t,pix]=fixmesh(p,t,ptol)
%FIXMESH  Remove duplicated/unused nodes and fix element orientation.
%   [P,T]=FIXMESH(P,T)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

if nargin<3, ptol=1024*eps; end
if nargin>=2 & (isempty(p) | isempty(t)), pix=1:size(p,1); return; end

snap=max(max(p,[],1)-min(p,[],1),[],2)*ptol;
[foo,ix,jx]=unique(round(p/snap)*snap,'rows');
p=p(ix,:);

if nargin>=2
    t=reshape(jx(t),size(t));
    
    [pix,ix1,jx1]=unique(t);
    t=reshape(jx1,size(t));
    p=p(pix,:);
    pix=ix(pix);
    
    if size(t,2)==size(p,2)+1
        flip=simpvol(p,t)<0;
        t(flip,[1,2])=t(flip,[2,1]);
    end
end
end

function [p,t]=distmesh2d(fd,fh,h0,bbox,pfix,varargin)
%DISTMESH2D 2-D Mesh Generator using Distance Functions.
%   [P,T]=DISTMESH2D(FD,FH,H0,BBOX,PFIX,FPARAMS)
%
%      P:         Node positions (Nx2)
%      T:         Triangle indices (NTx3)
%      FD:        Distance function d(x,y)
%      FH:        Scaled edge length function h(x,y)
%      H0:        Initial edge length
%      BBOX:      Bounding box [xmin,ymin; xmax,ymax]
%      PFIX:      Fixed node positions (NFIXx2)
%      FPARAMS:   Additional parameters passed to FD and FH
%
%   Example: (Uniform Mesh on Unit Circle)
%      fd=@(p) sqrt(sum(p.^2,2))-1;
%      [p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);
%
%   Example: (Rectangle with circular hole, refined at circle boundary)
%      fd=@(p) ddiff(drectangle(p,-1,1,-1,1),dcircle(p,0,0,0.5));
%      fh=@(p) 0.05+0.3*dcircle(p,0,0,0.5);
%      [p,t]=distmesh2d(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);
%
%   Example: (Polygon)
%      pv=[-0.4 -0.5;0.4 -0.2;0.4 -0.7;1.5 -0.4;0.9 0.1;
%          1.6 0.8;0.5 0.5;0.2 1;0.1 0.4;-0.7 0.7;-0.4 -0.5];
%      [p,t]=distmesh2d(@dpoly,@huniform,0.1,[-1,-1; 2,1],pv,pv);
%
%   Example: (Ellipse)
%      fd=@(p) p(:,1).^2/2^2+p(:,2).^2/1^2-1;
%      [p,t]=distmesh2d(fd,@huniform,0.2,[-2,-1;2,1],[]);
%
%   Example: (Square, with size function point and line sources)
%      fd=@(p) drectangle(p,0,1,0,1);
%      fh=@(p) min(min(0.01+0.3*abs(dcircle(p,0,0,0)), ...
%                   0.025+0.3*abs(dpoly(p,[0.3,0.7; 0.7,0.5]))),0.15);
%      [p,t]=distmesh2d(fd,fh,0.01,[0,0;1,1],[0,0;1,0;0,1;1,1]);
%
%   Example: (NACA0012 airfoil)
%      hlead=0.01; htrail=0.04; hmax=2; circx=2; circr=4;
%      a=.12/.2*[0.2969,-0.1260,-0.3516,0.2843,-0.1036];
%
%      fd=@(p) ddiff(dcircle(p,circx,0,circr),(abs(p(:,2))-polyval([a(5:-1:2),0],p(:,1))).^2-a(1)^2*p(:,1));
%      fh=@(p) min(min(hlead+0.3*dcircle(p,0,0,0),htrail+0.3*dcircle(p,1,0,0)),hmax);
%
%      fixx=1-htrail*cumsum(1.3.^(0:4)');
%      fixy=a(1)*sqrt(fixx)+polyval([a(5:-1:2),0],fixx);
%      fix=[[circx+[-1,1,0,0]*circr; 0,0,circr*[-1,1]]'; 0,0; 1,0; fixx,fixy; fixx,-fixy];
%      box=[circx-circr,-circr; circx+circr,circr];
%      h0=min([hlead,htrail,hmax]);
%
%      [p,t]=distmesh2d(fd,fh,h0,box,fix);

%
%   See also: MESHDEMO2D, DISTMESHND, DELAUNAYN, TRIMESH.

%   distmesh2d.m v1.1
%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

dptol=.001; ttol=.1; Fscale=1.2; deltat=.2; geps=.001*h0; deps=sqrt(eps)*h0;
densityctrlfreq=30;

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

% 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
p=p(rand(size(p,1),1)<r0./max(r0),:);                % Rejection method
if ~isempty(pfix), p=setdiff(p,pfix,'rows'); end     % Remove duplicated nodes
pfix=unique(pfix,'rows'); nfix=size(pfix,1);
p=[pfix; p];                                         % Prepend fix points
N=size(p,1);                                         % Number of points N

count=0;
pold=inf;                                            % For first iteration
clf,view(2),axis equal,axis off
while 1
  count=count+1;
  % 3. Retriangulation by the Delaunay algorithm
  if max(sqrt(sum((p-pold).^2,2))/h0)>ttol           % Any large movement?
    pold=p;                                          % Save current positions
    t=delaunayn(p);                                  % List of triangles
    pmid=(p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;    % Compute centroids
    t=t(feval(fd,pmid,varargin{:})<-geps,:);         % Keep interior triangles
    % 4. Describe each bar by a unique pair of nodes
    bars=[t(:,[1,2]);t(:,[1,3]);t(:,[2,3])];         % Interior bars duplicated
    bars=unique(sort(bars,2),'rows');                % Bars as node pairs
    % 5. Graphical output of the current mesh
    cla,patch('vertices',p,'faces',t,'edgecol','k','facecol',[.8,.9,1]);
    drawnow
  end

  % 6. Move mesh points based on bar lengths L and forces F
  barvec=p(bars(:,1),:)-p(bars(:,2),:);              % List of bar vectors
  L=sqrt(sum(barvec.^2,2));                          % L = Bar lengths
  hbars=feval(fh,(p(bars(:,1),:)+p(bars(:,2),:))/2,varargin{:});
  L0=hbars*Fscale*sqrt(sum(L.^2)/sum(hbars.^2));     % L0 = Desired lengths
  
  % Density control - remove points that are too close
  if mod(count,densityctrlfreq)==0 & any(L0>2*L)
      p(setdiff(reshape(bars(L0>2*L,:),[],1),1:nfix),:)=[];
      N=size(p,1); pold=inf;
      continue;
  end
  
  F=max(L0-L,0);                                     % Bar forces (scalars)
  Fvec=F./L*[1,1].*barvec;                           % Bar forces (x,y components)
  Ftot=full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2));
  Ftot(1:size(pfix,1),:)=0;                          % Force = 0 at fixed points
  p=p+deltat*Ftot;                                   % Update node positions

  % 7. Bring outside points back to the boundary
  d=feval(fd,p,varargin{:}); ix=d>0;                 % Find points outside (d>0)
  dgradx=(feval(fd,[p(ix,1)+deps,p(ix,2)],varargin{:})-d(ix))/deps; % Numerical
  dgrady=(feval(fd,[p(ix,1),p(ix,2)+deps],varargin{:})-d(ix))/deps; %    gradient
  dgrad2=dgradx.^2+dgrady.^2;
  p(ix,:)=p(ix,:)-[d(ix).*dgradx./dgrad2,d(ix).*dgrady./dgrad2];    % Project

  % 8. Termination criterion: All interior nodes move less than dptol (scaled)
  if max(sqrt(sum(deltat*Ftot(d<-geps,:).^2,2))/h0)<dptol, break; end
end

% Clean up and plot final mesh
[p,t]=fixmesh(p,t);
simpplot(p,t)
end

function simpplot(p,t,expr,bcol,icol,nodes,tris)

%   Copyright (C) 2004-2012 Per-Olof Persson. See COPYRIGHT.TXT for details.

dim=size(p,2);
switch dim
 case 2
  if nargin<4 | isempty(bcol), bcol=[.8,.9,1]; end
  if nargin<5 | isempty(icol), icol=[0,0,0]; end
  if nargin<6, nodes=0; end
  if nargin<7, tris=0; end
  
  trimesh(t,p(:,1),p(:,2),0*p(:,1),'facecolor',bcol,'edgecolor','k');
  if nodes==1
    line(p(:,1),p(:,2),'linest','none','marker','.','col',icol,'markers',24);
  elseif nodes==2
    for ip=1:size(p,1)
      txtpars={'fontname','times','fontsize',12};
      text(p(ip,1),p(ip,2),num2str(ip),txtpars{:});
    end
  end
  if tris==2
    for it=1:size(t,1)
      pmid=mean(p(t(it,:),:),1);
      txtpars={'fontname','times','fontsize',12,'horizontala','center'};
      text(pmid(1),pmid(2),num2str(it),txtpars{:});
    end
  end
  view(2)
  axis equal
  axis off
  ax=axis;axis(ax*1.001);
 case 3
  if nargin<4 | isempty(bcol), bcol=[.8,.9,1]; end
  if nargin<5 | isempty(icol), icol=[.9,.8,1]; end
  
  if size(t,2)==4
    tri1=surftri(p,t);
    if nargin>2 & ~isempty(expr)
      incl=find(eval(expr));
      t=t(any(ismember(t,incl),2),:);
      tri1=tri1(any(ismember(tri1,incl),2),:);
      tri2=surftri(p,t);
      tri2=setdiff(tri2,tri1,'rows');
      h=trimesh(tri2,p(:,1),p(:,2),p(:,3));
      set(h,'facecolor',icol,'edgecolor','k');
      hold on
    end
  else
    tri1=t;
    if nargin>2 & ~isempty(expr)
      incl=find(eval(expr));
      tri1=tri1(any(ismember(tri1,incl),2),:);
    end
  end
  h=trimesh(tri1,p(:,1),p(:,2),p(:,3));
  hold off
  set(h,'facecolor',bcol,'edgecolor','k');
  axis equal
  cameramenu
 otherwise
  error('Unimplemented dimension.');
end
end