computeAdditionalFacetMatrix

 search for additional facets in pairs of SEs

 THIS IS NO USER FUNCTION

• addMinkowskiHrep calculate the Minkowski sum of two polytopes in H-representation
• computeVolumeHrep compute the volume of a polytope in H-representation
• convertHrepToVrep convert the H-representation of a polytope into V-representation

• computeUnifyingDec create embedding unifying representation from decomposition

0001 function newFacetNormalMatrix = computeAdditionalFacetMatrix(dec, zerotol)
0002 % search for additional facets in pairs of SEs
0003 %
0004 % THIS IS NO USER FUNCTION
0005 
0006 % The elk-library: convex geometry applied to crystallization modeling.
0007 %   Copyright (C) 2014 Alexander Reinhold
0008 %
0009 % This program is free software: you can redistribute it and/or modify it
0010 %   under the terms of the GNU General Public License as published by the
0011 %   Free Software Foundation, either version 3 of the License, or (at your
0012 %   option) any later version.
0013 %
0014 % This program is distributed in the hope that it will be useful, but
0015 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0016 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0017 %   General Public License for more details.
0018 %
0019 % You should have received a copy of the GNU General Public License along
0020 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0021 
0022 % the following loop goes through all pairwise sums of structuring elements
0023 %   that are not adjacent (are not added to represent crystals / do not
0024 %   have a unified partition in common).
0025 % These pairs are expected to deliver a new set of facets. If no extra
0026 %   facets are created, an error is thrown. The facets are listed in
0027 %   newFacetNormalMatrix and a logical matrix tracks in which summations
0028 %   which facets were created. These sets are expected to be distinct,
0029 %   which is also explicitely checked (only pairwise).
0030 % From this facet creation matrix the new facet-groups are identified.
0031 [rowIndexVector, colIndexVector] = find(...
0032     (dec.seIncidenceMatrix + tril(ones(dec.nS) - eye(dec.nS))) == 0);
0033 newFacetNormalMatrix = [];
0034 % facetCreationMatrix = [];
0035 
0036 for iPair = 1:length(rowIndexVector)
0037     % The two summands
0038     hrepSummand1 = dec.se(rowIndexVector(iPair), :) * dec.mappingReducedToFull';
0039     hrepSummand1 = struct('A', dec.A, 'h', hrepSummand1');
0040     hrepSummand2 = dec.se(colIndexVector(iPair), :) * dec.mappingReducedToFull';
0041     hrepSummand2 = struct('A', dec.A, 'h', hrepSummand2');
0042     
0043     % Minkowski addition of pair
0044     hrepSum = addMinkowskiHrep(hrepSummand1, hrepSummand2);
0045     
0046     % it might happen that the resulting shape is a flat object. In this
0047     %   case we have (exactly) two opposing facet normals with distance 0
0048     %   and the other facet normals have absolutely no meaning. The facet
0049     %   distance might also be set off by translation.
0050     if computeVolumeHrep(hrepSum, zerotol) < zerotol
0051         vrepSum = convertHrepToVrep(hrepSum);
0052         % create a matrix with each row corresponding to a vertex and each
0053         %   column to a facet.
0054         facetFilter = (hrepSum.A * vrepSum.V')';
0055         % The (two) coumns where all values are equal correspond to the two
0056         %   facet normals we search for. The first is sufficient.
0057         for iFacet = 1:size(facetFilter, 2)
0058             if all(abs(facetFilter(:,iFacet) - facetFilter(1,iFacet)) < zerotol)
0059                 break;
0060             end
0061         end
0062         if iFacet == size(facetFilter, 2)
0063             error('elk:decomposition:numericalError', ...
0064                   'Could not find the facets of a flat object. This should not happen');
0065         end
0066         % obtain facet normals
0067         thisFacetNormalMatrix = [hrepSum.A(iFacet, :); ...
0068                               -1*hrepSum.A(iFacet, :)];
0069         
0070         % Filter original facets
0071         thisFacetNormalMatrix = eliminateRows(thisFacetNormalMatrix, ...
0072                                               dec.A, zerotol);
0073     else
0074         % Filter original facets
0075         thisFacetNormalMatrix = eliminateRows(hrepSum.A, dec.A, zerotol);
0076     end
0077  
0078     % throw error, if there are no additional facets
0079     if isempty(thisFacetNormalMatrix)
0080         error('elk:decomposition:numericalError', ...
0081             ['Structuring elements ' num2str(rowIndexVector(iPair)) ' and ' ...
0082              num2str(colIndexVector(iPair)) ' do not share a unified partition ' ...
0083              ' but also do not create additional facets by Minkowski addition.']);
0084     end
0085     
0086     % Filter already known facets
0087     [thisFacetNormalMatrix, thisCreatedFacetIndexVector] = eliminateRows(...
0088         thisFacetNormalMatrix, newFacetNormalMatrix, zerotol);
0089     
0090     % Append new facets
0091     thisCreatedFacetIndexVector = [thisCreatedFacetIndexVector, ...
0092           (1:size(thisFacetNormalMatrix, 1)) + size(newFacetNormalMatrix, 1)];
0093     newFacetNormalMatrix = [newFacetNormalMatrix; thisFacetNormalMatrix];
0094     
0095     % Enter information in facet creation matrix
0096 %     facetCreationMatrix(thisCreatedFacetIndexVector, iPair) = 1;
0097 end