Code covered by the BSD License

# Automated Failure Boundary Mapping

### Stuart Kozola (view profile)

18 Aug 2009 (Updated )

Demo files from July 21, 2009 webinar

boundarySearch(X0,lb,ub,limitFunction, varargin)
```function [xb,P, a] = boundarySearch(X0,lb,ub,limitFunction, varargin)
%BOUNDARYSEARCH searches a 2-D space to identify the failure region(s)
%  BOUNDARYSEARCH(X, LIMITFUNCTION) searches over a 2D space to identify
%  regions where the limit function returns (True) 1 if failure is present.
%  X0 (nx2) is initial search starting points of (x,y) pairs. LIMITFUNCTION
%  is a handle to the limit state function that that takes X as an input.
%
%  Optional input are:
%    'showplot'     show plot update of algorithm (voronoi cells and
%                   next points to update)
%    'ptol'         Perimeter change tolerance, followed by numeric value.
%
%  SYNTAX
%     [xb,P]   = boundarySearch(X0, limitFunction)
%     [xb,P,A] = boundarySearch(X0, limitFunction,'showplot','ptol',val)
%
%  Example
%    [x,y,z] = peaks(10);
%    X0 = [x(:) y(:)];
%    ub = [3 3];
%    lb = [-3 -3];
%    [xb P a] = boundarySearch(X0,lb,ub,@peaksLimitFunction,'showplot','ptol',0.1)
%
%

%   Copyright 2009 The MathWorks, Inc.

%% Input Parsing
debug = false; showplot = false;
ptol = 1e-3; vin = 1; maxregion = 5;
if nargin > 2
while vin <= length(varargin)
switch lower(varargin{vin})
case 'debug'
debug = true;
case 'showplot'
showplot = true;
case 'ptol'
vin = vin+1;
ptol = varargin{vin};
case 'maxregion'
vin = vin+1;
maxregion = varargin{vin};
otherwise
warning('BOUNDARYSEARCH:InvalidInputs',...
end
vin = vin+1;
end
end

%% Initialization
X = [];
failed = logical([]);
nextX = X0;
P = [];
Plast = 0;
deltaP = 100;
iter = 1;
while deltaP > ptol && iter < 1000
X = [X; nextX];
%% Call LimitFunction with New Inputs
nextFail = limitFunction(nextX);
failed = [failed; nextFail];
nextX = [];

%% Use Delaunay Triangulation
% Later we'll use the convex hull and voronoi diagrams to update search
% region
DT = DelaunayTri(X(:,1),X(:,2));
[V,R] = voronoiDiagram(DT);

if debug || showplot
voronoi(X(:,1),X(:,2)), drawnow
axis([lb(1) ub(1) lb(2) ub(2)])
end

%% Test for failure point detected
if sum(failed) > 0 % Found failure point(s)

%% Use clustering to identify number of differnt regions
zones = clusterdata([X,failed],maxregion);
if debug || showplot
%unique(zones)
hold on
scatter(X(:,1),X(:,2),50,zones,'filled')
axis square equal
axis([lb(1) ub(1) lb(2) ub(2)])
hold off
drawnow
end

%% Determine Number of Failure Regions
zoneNumbers = unique(zones.*failed);
zoneNumbers = zoneNumbers( zoneNumbers>0);

%% Find Failure Boundary
col = 'rgbmcy';
for zn = 1:length(zoneNumbers)
newDT = DelaunayTri(X(:,1),X(:,2));
[V,R] = voronoiDiagram(newDT);
% Remove other zones from newDT
indx = zones == zoneNumbers(zn);
if sum(indx) <= 2
% if only one/two points, grab voronoi cell points and continue
pts = find(indx > 0);
for i = 1:length(pts)
nextX = [nextX; V(R{pts(i)},:)];
end
continue
else
indx = find(~indx); % convert to integer index
newDT.X(indx,:) = [];
end
%% Find Each Failure Region's Convex Hull
try
[ch a(zn)] = convexHull(newDT);
catch eDT
% Check for empty triangulation (collinear points likely)
if strcmpi(eDT.identifier,'MATLAB:DelaunayTri:TriIsEmptyErrId')
pts = find(indx > 0);
for i = 1:length(pts)
nextX = [nextX; V(R{pts(i)},:)];
end
continue
else
rethrow(eDT)
end
end

% Calculate the perimeter of the convex hull
xb{zn} = newDT.X(ch,:);
dX = diff(newDT.X(ch,:));
P(zn)  = sum( sqrt( sum(dX.^2,2) ) );

% Make note of any points that are in the convex hull or on it
in = inpolygon(X(:,1),X(:,2),newDT.X(ch,1),newDT.X(ch,2));

% Use the point in the convex hull to find the voronoi cells
vc = R(in);

% Use the voronoi cell points as the next points to evaluate
idxNewPts = [];
for j = 1:length(vc)
idxNewPts = [idxNewPts, vc{j}];
end
idxNewPts = unique(idxNewPts);
newX = V(idxNewPts,:);

% Remove any points in/on the convex hull
in = inpolygon(newX(:,1),newX(:,2),newDT.X(ch,1),newDT.X(ch,2));
newX = newX(~in,:);

% Update next points array
nextX = [nextX;newX];

if showplot || debug
hold on
plot(newDT.X(ch,1),newDT.X(ch,2),col(zn))
plot(nextX(:,1),nextX(:,2),'r+')
axis([lb(1) ub(1) lb(2) ub(2)])
drawnow
if debug, pause, end
hold off
end
iter = iter+1;
end %zn loop
end
%% Remove INF and repeat (x,y) pairs
if isempty(nextX)      % No failure points found, use voronoi diagram points to update
nextX = V;
deltaP = deltaP - 1;
else
if isempty(P)
deltaP = deltaP -1;
else
deltaP = sum(P) - Plast;
Plast = sum(P);
end
end
nextX = unique(nextX,'rows');
remove = (nextX(:,1) == Inf | nextX(:,2) == Inf);
nextX(remove,:) = [];

%% Remove any Points Not In Bounds
% Remove x's
nextX = nextX(nextX(:,1) <= ub(1),:);
nextX = nextX(nextX(:,1) >= lb(1),:);
nextX = nextX(nextX(:,2) <= ub(2),:);
nextX = nextX(nextX(:,2) >= lb(2),:);

end %while
%repeat```