compute function handles for edge lengths and vrep it is assumed that the edges in hrep exist and remain for inputs to the function handles. edgeLengthFunctional(h) returns a vector of edge lengths for P(H) and vrepFunctional(h) returns the corresponding vrep. THIS IS NO USER FUNCTION
0001 function [edgeLengthFunctional vrepFunctional incidenceMatrix] = ... 0002 computeEdgeLengthFunctional(hrep, vrep, zerotol) 0003 % compute function handles for edge lengths and vrep 0004 % 0005 % it is assumed that the edges in hrep exist and remain for inputs to the 0006 % function handles. edgeLengthFunctional(h) returns a vector of edge 0007 % lengths for P(H) and vrepFunctional(h) returns the corresponding vrep. 0008 % 0009 % THIS IS NO USER FUNCTION 0010 0011 % The elk-library: convex geometry applied to crystallization modeling. 0012 % Copyright (C) 2012 Alexander Reinhold 0013 % 0014 % This program is free software: you can redistribute it and/or modify it 0015 % under the terms of the GNU General Public License as published by the 0016 % Free Software Foundation, either version 3 of the License, or (at your 0017 % option) any later version. 0018 % 0019 % This program is distributed in the hope that it will be useful, but 0020 % WITHOUT ANY WARRANTY; without even the implied warranty of 0021 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0022 % General Public License for more details. 0023 % 0024 % You should have received a copy of the GNU General Public License along 0025 % with this program. If not, see <http://www.gnu.org/licenses/> 0026 0027 if isempty(hrep.A) 0028 %% empty hrep is treated differently 0029 edgeLengthFunctional = @(h)([]); 0030 vrepFunctional = @(h)(struct('V', [0 0])); 0031 incidenceMatrix = zeros([0 1]); 0032 else 0033 %% basic information 0034 % the procedure to identify vertex/edge pairs looks complicated but was 0035 % necessary because of tolerance problems. The incidence matrix 0036 % should only contain 2 entries in each row AND in each column. 0037 distanceMatrix = abs(bsxfun(@minus, hrep.A * vrep.V', hrep.h)); 0038 % select first vertex 0039 tmpMinVec = min(distanceMatrix, [], 1); 0040 incidenceMatrix = bsxfun(@eq, distanceMatrix, tmpMinVec); 0041 % this distance will not be selected again 0042 distanceMatrix(incidenceMatrix) = 2*max(distanceMatrix(:)); 0043 % try again 0044 tmpMinVec = min(distanceMatrix, [], 1); 0045 % but avoid selecting an additional vertex/edge when there are already 0046 % 2 selected 0047 tmpMinVec(sum(incidenceMatrix, 1) > 1) = 2*max(distanceMatrix(:)); 0048 incidenceMatrix = incidenceMatrix | ... 0049 bsxfun(@eq, distanceMatrix, tmpMinVec); 0050 0051 % incidenceMatrix = abs(bsxfun(@minus, hrep.A * vrep.V', hrep.h)) < zerotol; 0052 % % avoid multiple entries >> I can at least correct a single fault 0053 % if sum(sum(incidenceMatrix, 1) > 2) == 1 0054 % xIndex = find(sum(incidenceMatrix, 1) == 3); 0055 % yIndex = find(sum(incidenceMatrix, 2) == 3); 0056 % incidenceMatrix(yIndex, xIndex) = 0; 0057 % end 0058 if sum(sum(incidenceMatrix, 1) <= 2) == 1 0059 error('elk:boundary:internalError', ['Basically tolerance ' ... 0060 'problems. This should not happen']); 0061 end 0062 0063 %% compute vrep functional 0064 coefficientMatrixX = zeros(size(incidenceMatrix, 2),... 0065 size(incidenceMatrix, 1)); 0066 coefficientMatrixY = coefficientMatrixX; 0067 for iVertex = 1:size(incidenceMatrix, 2) 0068 % from the incidence matrix, we easily know the adjacent edges: 0069 edgeIndexVector = find(incidenceMatrix(:, iVertex)); 0070 % now, the vertex [x, y] is defined by: 0071 % [a1, a2]' * [x, y]' = [h1, h2]' 0072 % which becomes: 0073 % [x y] = [h1, h2] * inv([a1, a2]) 0074 % and is rewritten by: 0075 % PI = [p1, p2] = inv([a1, a2]) 0076 % x = p1' * [h1, h2]' 0077 % y = p2' * [h1, h2]' 0078 PI = inv(hrep.A(edgeIndexVector, :)'); 0079 coefficientMatrixX(iVertex, edgeIndexVector(1)) = PI(1, 1); 0080 coefficientMatrixX(iVertex, edgeIndexVector(2)) = PI(2, 1); 0081 coefficientMatrixY(iVertex, edgeIndexVector(1)) = PI(1, 2); 0082 coefficientMatrixY(iVertex, edgeIndexVector(2)) = PI(2, 2); 0083 end 0084 vrepFunctional = @(h)(struct('V', ... 0085 [coefficientMatrixX*h, coefficientMatrixY*h])); 0086 0087 %% compute edgeLengthFunctional 0088 coefficientMatrix = zeros(size(hrep.A, 1)); 0089 for iEdge = 1:size(hrep.A, 1) 0090 % identify bounding edges 0091 vertexIndexVector = find(incidenceMatrix(iEdge, :)); 0092 thisFilter = ones(1, size(hrep.A, 1)) == 1; 0093 thisFilter(iEdge) = false; 0094 thisFilter = thisFilter & (... 0095 incidenceMatrix(:,vertexIndexVector(1))' | ... 0096 incidenceMatrix(:,vertexIndexVector(2))'); 0097 edgeIndexVector = find(thisFilter); 0098 % assume, the whole polytope is moved such that the edge goes 0099 % through the current support line: 0100 % hM = h - A*(a0*h0) 0101 % now we can assign a distance of the two vertices from the origin, 0102 % when they are projected to the support line (oblique 0103 % projection): 0104 % Px = a0Perp * inv(ax'*a0Perp) * ax' 0105 % = (a0Perp*ax') / (ax'*a0Perp) 0106 % d1 = a0Perp' * P1 * a1 * hM1 0107 % d2 = a0Perp' * P2 * a2 * hM2 0108 % where aPerp is perpendicular to (a). The edge length is then 0109 % either (d1-d2) or (d2-d1). But is must be positive for lam=0. 0110 % This reduces to: 0111 % d1 = (1/(a1'*a0Perp)) * (h1 - a1'*a0*h0) 0112 % d2 = (1/(a2'*a0Perp)) * (h2 - a2'*a0*h0) 0113 % and with an arbitrary sign assumed: 0114 % d = 1/(a2'*a0Perp)*h2 ... 0115 % -1/(a1'*a0Perp)*h1 ... 0116 % +(a1'/(a1'*a0Perp) - a2'/(a2'*a0Perp))*a0*h0 0117 a0 = hrep.A(iEdge, :)'; 0118 a1 = hrep.A(edgeIndexVector(1), :)'; 0119 a2 = hrep.A(edgeIndexVector(2), :)'; 0120 a0Perp = [a0(2); -1*a0(1)]; 0121 coefficientMatrix(iEdge, iEdge) = (... 0122 a1'/(a1'*a0Perp) - a2'/(a2'*a0Perp))*a0; 0123 coefficientMatrix(iEdge, edgeIndexVector(1)) = -1/(a1'*a0Perp); 0124 coefficientMatrix(iEdge, edgeIndexVector(2)) = 1/(a2'*a0Perp); 0125 % correct sign of edge length 0126 if coefficientMatrix(iEdge, :) * hrep.h < 0 0127 coefficientMatrix(iEdge, :) = -1*coefficientMatrix(iEdge, :); 0128 end 0129 end 0130 0131 edgeLengthFunctional = @(h)((coefficientMatrix * h)'); 0132 end 0133 0134