calculate constraints to fix the position of the polytope
Syntax: constraintMatrix = fixPosition(A)
The position of polytopes is non-relevant for the topology (existence of
vertices, edges or facets). Hence, all algorithms connected to
Minkowski decomposition or the 'morphological' decomposition of the
state space require the position of each valid polytope to be fixed.
Whenever this is not the case, additional constraints are employed
without loss of generality.
This function fixes certain facets to a facet distance of 0 so that:
constraintMatrix * h = 0
ensures a fixed position of the polytope.
fixPosition(A, originalConstraintMatrix) searches for additional
constraints to fix the position of the polytope. Only the additional
constaints are returned. The condition above is then:
[constraintMatrix; originalConstraintMatrix] * h = 0
With [~, mappingMatrix] = fixPosition(A), a mappingMatrix is returned
that allows to move a given polytope (vector h) such that the new
polytope (vector hNew) fulfills the fixed position constraint. Hence,
the code for adopting a given polytope is:
[constraintMatrix, mappingMatrix] = fixPosition(hrep.A);
hrepPositionFixed = moveHrep(hrep, mappingMatrix * hrep.h);
It then holds:
constraintMatrix * hrepPositionFixed.h = 0
fixPosition(.., zerotol) sets the tolerance for underlying algorithms.
Default: elkZerotol
See also: obtainDec, decomposeMinkowski0001 function varargout = fixPosition(A, inputConstraintMatrix, zerotol) 0002 % calculate constraints to fix the position of the polytope 0003 % 0004 % Syntax: constraintMatrix = fixPosition(A) 0005 % 0006 % The position of polytopes is non-relevant for the topology (existence of 0007 % vertices, edges or facets). Hence, all algorithms connected to 0008 % Minkowski decomposition or the 'morphological' decomposition of the 0009 % state space require the position of each valid polytope to be fixed. 0010 % Whenever this is not the case, additional constraints are employed 0011 % without loss of generality. 0012 % 0013 % This function fixes certain facets to a facet distance of 0 so that: 0014 % constraintMatrix * h = 0 0015 % ensures a fixed position of the polytope. 0016 % 0017 % fixPosition(A, originalConstraintMatrix) searches for additional 0018 % constraints to fix the position of the polytope. Only the additional 0019 % constaints are returned. The condition above is then: 0020 % [constraintMatrix; originalConstraintMatrix] * h = 0 0021 % 0022 % With [~, mappingMatrix] = fixPosition(A), a mappingMatrix is returned 0023 % that allows to move a given polytope (vector h) such that the new 0024 % polytope (vector hNew) fulfills the fixed position constraint. Hence, 0025 % the code for adopting a given polytope is: 0026 % [constraintMatrix, mappingMatrix] = fixPosition(hrep.A); 0027 % hrepPositionFixed = moveHrep(hrep, mappingMatrix * hrep.h); 0028 % It then holds: 0029 % constraintMatrix * hrepPositionFixed.h = 0 0030 % 0031 % fixPosition(.., zerotol) sets the tolerance for underlying algorithms. 0032 % Default: elkZerotol 0033 % 0034 % See also: obtainDec, decomposeMinkowski 0035 0036 % The elk-library: convex geometry applied to crystallization modeling. 0037 % Copyright (C) 2012 Alexander Reinhold 0038 % 0039 % This program is free software: you can redistribute it and/or modify it 0040 % under the terms of the GNU General Public License as published by the 0041 % Free Software Foundation, either version 3 of the License, or (at your 0042 % option) any later version. 0043 % 0044 % This program is distributed in the hope that it will be useful, but 0045 % WITHOUT ANY WARRANTY; without even the implied warranty of 0046 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0047 % General Public License for more details. 0048 % 0049 % You should have received a copy of the GNU General Public License along 0050 % with this program. If not, see <http://www.gnu.org/licenses/>. 0051 0052 nFacet = size(A, 1); 0053 nDim = size(A, 2); 0054 0055 if ~exist('inputConstraintMatrix', 'var') 0056 inputConstraintMatrix = zeros(0, nFacet); 0057 end 0058 if ~exist('zerotol', 'var') 0059 zerotol = elkZerotol; 0060 end 0061 0062 %% subsequently increase the rank of Mc*A 0063 % ..by fixing selected facet distances to zero. 0064 0065 thisRank = rank(inputConstraintMatrix * A, zerotol); 0066 positionConstraintMatrix = zeros(nFacet); 0067 while thisRank < nDim 0068 0069 % singular value decomposition 0070 [U S V] = svd([inputConstraintMatrix; positionConstraintMatrix] * A); 0071 S(isZero(S(:), zerotol)) = 0; 0072 thisRank = sum(S(:) ~= 0); 0073 0074 % check rank 0075 if thisRank == nDim 0076 break; 0077 end 0078 0079 % in which direction can I move? 0080 thisDirection = V(:, end); 0081 0082 % which facet normals point into this direction? 0083 facetNormalDotProduct = A * thisDirection; 0084 if isAllZero(facetNormalDotProduct, zerotol) 0085 error('elk:decomposition:numericalProblem', ... 0086 'Failed to fix position. This should not happen.'); 0087 end 0088 0089 % take one of them as new constraint 0090 thisFacetIndex = find(... 0091 abs(facetNormalDotProduct) == max(abs(facetNormalDotProduct)), ... 0092 1, 'first'); 0093 0094 % add constraint 0095 positionConstraintMatrix(thisRank + 1, thisFacetIndex) = 1; 0096 0097 end 0098 varargout{1} = positionConstraintMatrix; 0099 0100 %% Offset generation 0101 % Now the polytope must be moved to fulfill that additional constraint. 0102 % This means we search an x0 with: Mc * (h + A*x0) = 0, which translates 0103 % to: Mc*A*x0 = -Mc*h. This can be solved by: x0 = (Mc*A)\(-Mc*h). 0104 % However, the output of this will be a matrix that transforms directly 0105 % from h to x0, which is: -1*pinv(Mc*A)*Mc. 0106 0107 if nargout == 2 0108 varargout{2} = -1 * pinv(positionConstraintMatrix * A) * ... 0109 positionConstraintMatrix; 0110 end 0111