function dec = addFacetValidityData(dec, validityConeData, zerotol, verbose)
% add data to identify invalid facets
%
% THIS FUNCTION IS NO USER FUNCTION.

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

%% data structure
%     .A       - boundary facet normals in nC-dimensional space
%     .responsibleMatrix - a matrix that indicates in each row a constraint
%                in A (first column) that must be fulfilled by a facet 
%                (second column)

%% generate facet validity data
% the responsibleMatrix indicates where the facet grows out; the facet
%   (nH+1) corresponds to the empty-polytope condidions
boundaryData.A = nan(0, dec.nC);
boundaryData.responsibleMatrix = nan(0, 2);
ineqFilterBase = validityConeData.responsibleFacetVector < 0;

if verbose,fprintf('......');end

for iFacet = 1:(dec.nH+1)
    if verbose
       fprintf('\b\b\b\b\b\b%5.1f%%', iFacet/(dec.nH+1)*100);
    end
    
    % filter required conditions (in case nH+1, ineqFilterBase is retained)
    thisIneqFilter = ineqFilterBase | ...
        (validityConeData.responsibleFacetVector == iFacet);
    
    % reduce conditions
    thisCone.A = unique(validityConeData.A(thisIneqFilter, :), 'rows');
    thisCone.h = zeros(size(thisCone.A, 1), 1);
    thisCone = reduceHrep(thisCone);
    
    % concatenate constraints
    theseIneqIdxVector = size(boundaryData.A, 1) + (1:size(thisCone.A, 1))';
    boundaryData.A = [boundaryData.A; thisCone.A];
    boundaryData.responsibleMatrix = [boundaryData.responsibleMatrix; ...
        [theseIneqIdxVector, theseIneqIdxVector*0 + iFacet]];
end
if verbose, fprintf('\b\b\b\b\b\b');end

% make constraints unique, again
[boundaryData.A dump newIndexVector] = unique(boundaryData.A, 'rows');
boundaryData.responsibleMatrix(:,1) = newIndexVector(...
    boundaryData.responsibleMatrix(:,1));
boundaryData.responsibleMatrix = unique(boundaryData.responsibleMatrix, 'rows');

% assign result
dec.facetValidityData = boundaryData;