add data to identify vanished facets and map to validity cone THIS FUNCTION IS NO USER FUNCTION.
0001 function dec = addValidityMappingData(dec, zerotol) 0002 % add data to identify vanished facets and map to validity cone 0003 % 0004 % THIS FUNCTION IS NO USER FUNCTION. 0005 0006 % The elk-library: convex geometry applied to crystallization modeling. 0007 % Copyright (C) 2012 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 %% data structure 0023 % Draft of the gathered information in boundaryData: 0024 % .A - boundary facet normals in nC-dimensional space 0025 % .nBoundary - total number of rows in A 0026 % .responsibleVector - a vector that indicates the facet (group) that 0027 % disappears given corresponding to the constraint in A. 0028 % .directionMatrix - the collection of directions (column vectors) in 0029 % hC-space which is used to correct a specific boundary 0030 % (column). It contains nan, when mapping is not possible 0031 % without changing a valid facet and inf when mapping 0032 % to the origin would be appropriate. In the latter case, 0033 % the boundary indicates the disappearance of all facets at 0034 % once. 0035 0036 %% prepare setup 0037 % copy constraints 0038 dec.confinementMappingData.A = dec.facetValidityData.A; 0039 % predefine assigned directions (per constraint) 0040 dec.confinementMappingData.responsibleVector = ... 0041 nan(1, size(dec.confinementMappingData.A, 1)); 0042 % predefine possible projection sets (per facet) 0043 dec.confinementMappingData.projectionData = struct([]); 0044 % predefine the directions 0045 dec.confinementMappingData.directionMatrix = nan([dec.nC dec.nC+1]); 0046 % local copy 0047 responsibleMatrix = dec.facetValidityData.responsibleMatrix; 0048 0049 %% generate zero-projection data 0050 dec.confinementMappingData.responsibleVector(... 0051 responsibleMatrix(responsibleMatrix(:, 2) == (dec.nH + 1), 1)) = dec.nC+1; 0052 dec.confinementMappingData.projectionData(dec.nC+1).A = zeros(0, dec.nC); 0053 dec.confinementMappingData.projectionData(dec.nC+1).responsibleMatrix = ... 0054 zeros(0, 3); 0055 dec.confinementMappingData.projectionData(dec.nC+1).projectionList = ... 0056 {zeros(dec.nC)}; 0057 % delete all entries for these constraints, they appear for all facets 0058 for iConstraint = find(dec.confinementMappingData.responsibleVector == (dec.nC+1)) 0059 responsibleMatrix(responsibleMatrix(:, 1) == iConstraint, :) = []; 0060 end 0061 0062 %% generate facet group projections 0063 % construct unique non-zero rows 0064 [directionMatrix facetIndexVector, directionIndexVector] = ... 0065 unique(dec.mappingReducedToFull, 'rows'); 0066 % Zero rows in the direction matrix occure when the position was fixed 0067 % by the decomposition algorithm. These direction must be considered 0068 % with some precaution. They must be treated separately. Therefore, we 0069 % identify the original zero rows and multiply them appropriately. 0070 zeroRowFilter = all(dec.mappingReducedToFull < zerotol, 2); 0071 zeroRowIndexFilter = find(zeroRowFilter); 0072 % if there is only one zero row then the above is already correct. If there 0073 % are more zero rows, we must append some zero-row directions. 0074 if sum(zeroRowFilter) > 1 0075 % find the current zero-row index 0076 thisZeroIndex = find(all(directionMatrix < zerotol, 2)); 0077 % add zero rows 0078 directionMatrix = [directionMatrix zeros(sum(zeroRowFilter)-1, dec.nC)]; 0079 % substitute direction indices 0080 zeroDirIndexVector = find(directionIndexVector == thisZeroIndex); 0081 directionIndexVector(zeroDirIndexVector) = zeroRowIndexFilter; 0082 end 0083 0084 nDir = size(directionMatrix, 1); 0085 0086 for iDir = 1:nDir 0087 0088 % --- 0089 % identify facet(s) for direction 0090 thisFacetIndex = find(directionIndexVector == iDir, 1, 'first'); 0091 isZeroDir = any(zeroRowIndexFilter == thisFacetIndex); 0092 0093 % --- 0094 % identify mapping direction (null space) 0095 thisDirectionFilter = true(1, nDir); 0096 thisDirectionFilter(iDir) = false; 0097 % Things are different when iDir indicates a zero row, therefore the 0098 % separate treatment: 0099 if ~isZeroDir 0100 % The projection direction equals the null space of the oblique 0101 % projection that is to be constructed. While this direction is 0102 % not important for the handling of the proper representation, it 0103 % becomes valuable as soon as an improper representation is 0104 % constructed. 0105 dec.confinementMappingData.directionMatrix(:, iDir) = null(... 0106 directionMatrix(thisDirectionFilter, :)); 0107 0108 % For the construction of the obligue projection, we require the 0109 % orthogonal complement of the null space (see wikipedia). 0110 %! DEV: this line was used beforehand.. ..and can't be right in 0111 % general. It is only a proper result if the selected 0112 % direction is already orthogonal to all other directions. 0113 % thisOrthNull = orth(directionMatrix(thisDirectionFilter, :)'); 0114 % we use this new line: 0115 thisOrthNull = null(dec.confinementMappingData.directionMatrix(:, iDir)'); 0116 else 0117 % Things are different for a zero row that was introduced to fix the 0118 % position. The disappearance of this facet also implies that the 0119 % position is not fixed anymore. The nullspace is, hence, computed 0120 % from this degree of freedom. Assuming this facet has 0121 % disappeared, we temporarly ignore the facet 0122 tmpA = dec.A; tmpA(thisFacetIndex, :) = []; 0123 tmpMc = dec.constraintMatrix; tmpMc(:, thisFacetIndex) = []; 0124 tmpMapToH = dec.mappingReducedToFull; tmpMapToH(thisFacetIndex, :) = []; 0125 tmpMapFromH = pinv(tmpMapToH); 0126 % obtain free direction via singular value decomposition 0127 [U S V] = svd(tmpMc * tmpA); 0128 % The direction in which we can move freely is: 0129 % x0 = V(:, end) 0130 % which is applied according to: 0131 % h' = h + A*x0. 0132 % Scince h fulfills tmpMc*h=0 and h' is computed so that it 0133 % fulfills tmpMc*h', its difference also fulfills this 0134 % constraint so that hC is: 0135 thisDirHc = tmpMapFromH*tmpA*V(:, end); 0136 % This is the null space and we need the orthogonal complement. 0137 thisOrthNull = null(thisDirHc'); 0138 end 0139 0140 % --- 0141 % determine range and construct projection(s) 0142 0143 % identify assigned constraints (range) 0144 dec.confinementMappingData.responsibleVector(responsibleMatrix(... 0145 responsibleMatrix(:, 2) == thisFacetIndex, 1)) = iDir; 0146 thisConstraintFilter = dec.confinementMappingData.responsibleVector == iDir; 0147 0148 % fill projection data (one for each constraint) 0149 thisProjectionData = appendProjectionData([], dec.nC); 0150 for iConstraint = find(thisConstraintFilter) 0151 % the facet is spanned by the following base vectors that are 0152 % constructed from the orthogonal null space, spanned by the 0153 % facet normal. Hence, we take the null space of the null space 0154 % for the range of the projection onto the facet. 0155 thisRange = null(dec.confinementMappingData.A(iConstraint, :)); 0156 0157 % the obligue projection is now: 0158 thisProjection = thisRange * inv(... 0159 thisOrthNull' * thisRange) * thisOrthNull'; 0160 % it must map to the group validity cone: 0161 % Av * P * hC <= 0 0162 thisConstraint = ... 0163 dec.confinementMappingData.A(thisConstraintFilter, :) * ... 0164 thisProjection; 0165 % .. and - the projection must actually be valid. This is redundant 0166 % to the constraints that are checked in confinementMappingData.A 0167 % but this check is required for completeness when dealing with 0168 % projections of proper decompositions. There, the whole validity 0169 % of this projection is required. 0170 thisConstraint = [thisConstraint; ... 0171 -1*dec.confinementMappingData.A(iConstraint, :)]; 0172 % assign results 0173 thisProjectionData = appendProjectionData(thisProjectionData, ... 0174 thisConstraint, thisProjection, zerotol); 0175 end 0176 dec.confinementMappingData.projectionData(iDir) = thisProjectionData; 0177 end 0178 0179 % add something that is used only for improper representations 0180 dec.confinementMappingData.ignoreFadeOut = [];