convert the V-representation of a polytope into a H-representation Syntax: hrep = convertVrepToHrepByConvhulln(vrep, zerotol) The input should already be scaled by calling function. THIS FUNCTION IS NO USER FUNCTION.
0001 function hrep = convertVrepToHrepByConvhulln(vrep, zerotol) 0002 % convert the V-representation of a polytope into a H-representation 0003 % 0004 % Syntax: hrep = convertVrepToHrepByConvhulln(vrep, zerotol) 0005 % 0006 % The input should already be scaled by calling function. 0007 % 0008 % THIS FUNCTION IS NO USER FUNCTION. 0009 0010 % The elk-library: convex geometry applied to crystallization modeling. 0011 % Copyright (C) 2012 Alexander Reinhold 0012 % 0013 % This program is free software: you can redistribute it and/or modify it 0014 % under the terms of the GNU General Public License as published by the 0015 % Free Software Foundation, either version 3 of the License, or (at your 0016 % option) any later version. 0017 % 0018 % This program is distributed in the hope that it will be useful, but 0019 % WITHOUT ANY WARRANTY; without even the implied warranty of 0020 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0021 % General Public License for more details. 0022 % 0023 % You should have received a copy of the GNU General Public License along 0024 % with this program. If not, see <http://www.gnu.org/licenses/>. 0025 0026 %% handle low-dim cases 0027 nDim = size(vrep.V, 2); 0028 if nDim < 1 0029 hrep.A = []; 0030 hrep.h = []; 0031 return 0032 elseif nDim == 1 0033 hrep.A = [1; -1]; 0034 hrep.h = [max(vrep.V(:)); -1*min(vrep.V(:))]; 0035 return 0036 end 0037 0038 %% call convhulln 0039 patchIndexMatrix = convhulln(vrep.V, {'Qt', 'Qs'}); 0040 nPatch = size(patchIndexMatrix, 1); 0041 0042 %% identify facets 0043 facetNormalMatrix = nan(nPatch, nDim); 0044 facetDistanceVector = nan(nPatch, 1); 0045 for iPatch = 1:nPatch 0046 % generate matrix of vertices 0047 thisVertexMatrix = vrep.V(patchIndexMatrix(iPatch, :), :); 0048 % move so that one vertex equals the origin 0049 thisVertexMatrix = thisVertexMatrix - ones(nDim, 1) * thisVertexMatrix(1,:); 0050 0051 % an additional check of dimensionality could be implemented here it 0052 % should hold: 0053 % rank(thisVertexMatrix) == nDim - 1 0054 0055 % the normal vector v will fulfill the following equation: 0056 % thisV * v = 0 0057 % Hence, it equals the eigenvector that corresponds to the eigenvalue 0058 % of zero. While this is true, there are frequently multiple 0059 % eigenvalues and the precission of the corresponsing eigenvectors is 0060 % extremely bad. 0061 % To overcome this accuracy problem, the singular value decomposition 0062 % is applied. 0063 [~, thisS, thisV] = svd(thisVertexMatrix); 0064 thisS = diag(thisS); 0065 thisVecIndex = find(abs(thisS) < zerotol); 0066 if numel(thisVecIndex) ~= 1 0067 % error('elk:polytope:numericalError', ... 0068 % 'This patch is empty. This should not happen.'); 0069 facetNormalMatrix(iPatch, :) = zeros(1, nDim); 0070 facetDistanceVector(iPatch) = 0; 0071 continue; 0072 else 0073 facetNormalMatrix(iPatch, :) = thisV(:, thisVecIndex)'; 0074 end 0075 vertexDistanceVector = facetNormalMatrix(iPatch, :) * vrep.V'; 0076 facetDistanceVector(iPatch) = ... 0077 vertexDistanceVector(patchIndexMatrix(iPatch, 1)); 0078 0079 % The normal vector should be an outer normal. This cannot be decided 0080 % on the facet vertices itself. A vertex that belongs to the polytope 0081 % aside the facet is required. This vertes must lie in the negative 0082 % direction of the normal vector. 0083 vertexDistanceVector = vertexDistanceVector - facetDistanceVector(iPatch); 0084 maxDistance = max(vertexDistanceVector); 0085 minDistance = min(vertexDistanceVector); 0086 if maxDistance > zerotol 0087 % change sign, if an additional vertex lies in the direction of the 0088 % normal vector 0089 facetDistanceVector(iPatch) = -1 * facetDistanceVector(iPatch); 0090 facetNormalMatrix(iPatch, :) = -1 * facetNormalMatrix(iPatch, :); 0091 % in this case minDistance should be around zero 0092 if ~(abs(minDistance) < zerotol) 0093 error('elk:polytope:numericalError', ... 0094 'The polytope has vertices inside and ouside of an identified halfspace'); 0095 end 0096 elseif ~(minDistance < -zerotol) 0097 error('elk:polytope:numericalError', ... 0098 'The polytope has vertices inside and ouside of an identified halfspace'); 0099 end 0100 end 0101 0102 %% make unique normal vectors 0103 % This happens, because the facets are triangulated and the normal vectors 0104 % for these simplices are identified. 0105 [~, idx] = eliminateRows(facetNormalMatrix, ... 0106 [zeros(1, nDim); ... 0107 facetNormalMatrix], zerotol); 0108 idx(idx == 1) = []; 0109 idx = idx - 1; 0110 idx = unique(idx); 0111 0112 %% write output 0113 hrep.A = facetNormalMatrix(idx, :); 0114 hrep.h = facetDistanceVector(idx);