function d = dexpr ( p, fin, nit, alpha ) %*****************************************************************************80 % %% DEXPR computes the signed distance for a general implicit expression FIN. % % Discussion: % % The assumption is that the boundary is defined as a level set function, % f(x,y) = 0, with f(x,y) < 0 for points inside the region, and 0 < f(x,y) % for points outside the region. % % To find the point (X,Y) on the level set that is nearest to a given point % (X0,Y0), we require that f(X,Y) be zero and that the vector from (X,Y) to % (X0,Y0) be parallel to the gradient (fx,fy) at (X,Y), which defines a pair % of nonlinear equations L(X,Y): % % f(X,Y) = 0 % ( X - X0 ) * fy - ( Y - Y0 ) * fx = 0 % % These equations are solved by Newton's method. % % The use of this function requires access to the Symbolic Toolbox, although % it would be possible to rewrite it to accept derivatives of FIN as input % as well. This routine is quite slow relative to a C version. % % Licensing: % % (C) 2004 Per-Olof Persson. % See COPYRIGHT.TXT for details. % % Reference: % % Per-Olof Persson and Gilbert Strang, % A Simple Mesh Generator in MATLAB, % SIAM Review, % Volume 46, Number 2, June 2004, pages 329-345. % % Parameters: % % Input, real P(NP,2), the point coordinates. % % Input, string FIN, the expression that defines F(X,Y). % % Input, integer NIT, the number of iterations. The default value of NIT is 20. % % Input, real ALPHA, the damping parameter for the Newton method. % The default value of ALPHA is 0.1. % % Output, real D(NP), the signed distance of each point to the boundary. % if ( nargin < 3 ) nit = 20; end if ( nargin < 4 ) alpha = 0.1; end fx = inline ( vectorize(maple(['diff(',fin,',x)'])), 'x', 'y' ); fy = inline ( vectorize(maple(['diff(',fin,',y)'])), 'x', 'y' ); fxx = inline ( vectorize(maple(['diff(',fin,',x\$2)'])), 'x', 'y' ); fyy = inline ( vectorize(maple(['diff(',fin,',y\$2)'])), 'x', 'y' ); fxy = inline ( vectorize(maple(['diff(diff(',fin,',x),y)'])), 'x', 'y' ); f = inline ( vectorize(fin), 'x', 'y' ); x0 = p(:,1); y0 = p(:,2); x = x0; y = y0; for it = 1 : nit cf = feval(f,x,y); cfx = feval(fx,x,y); cfy = feval(fy,x,y); cfxx = feval(fxx,x,y); cfxy = feval(fxy,x,y); cfyy = feval(fyy,x,y); % % Evaluate L: % F1 = cf; F2 = (x-x0) .* cfy - (y-y0) .* cfx; % % Compute the jacobian of L: % % fx fy % fy + ( x - x0 ) * fxy - ( y - y0 ) * fxx -fx - ( y - y0 ) * fxy + ( x - x0 ) * fyy % J11 = cfx; J12 = cfy; J21 = cfy + (x-x0) .* cfxy - (y-y0) .* cfxx; J22 = -cfx - (y-y0) .* cfxy + (x-x0) .* cfyy; detJ = J11 .* J22 - J12 .* J21; detJ(detJ==0) = inf; % % Apply the Newton correction to (X,Y). % x = x - alpha * ( J22 .* F1 - J21 .* F2 ) ./ detJ; y = y - alpha * ( -J12 .* F1 + J11 .* F2 ) ./ detJ; end % % Having found the nearest point (X,Y) on the level set, we can compute the signed distance. % d = sqrt((x-x0).^2+(y-y0).^2) .* sign(feval(f,x0,y0)); return end