addConfinementPartitionData

 add confined partitions

 THIS IS NO USER FUNCTION

• appendPartitionDataForEmbedding append embedding (unified) partition data
• appendUnifiedPartition add unified partition, calculate simplex partitions and perform book
• closeCone compute the closure of the given cone
• openCone remove the prior closure of the cone
• scaleToClose scale input matrix of hC-vectors to the closing condition of the cone
• scaleToReal scale input matrix of hC-vectors to some real geometry
• computeVolumeVrep compute the volume of a polytope in V-representation
• convertHrepToVrep convert the H-representation of a polytope into V-representation
• convertVrepToHrep convert a V-representation into a H-representation

• projectCompleteDec project a complete proper representation

0001 function dec = addConfinementPartitionData(dec, partitionData, zerotol)
0002 % add confined partitions
0003 %
0004 % THIS IS NO USER FUNCTION
0005 
0006 % The elk-library: convex geometry applied to crystallization modeling.
0007 %   Copyright (C) 2013 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 %% construct and accumulate SEs
0023 % each partition contains only the open hrep, we need to construct the
0024 %   vrep, extract the ray elements and append them as SEs
0025 dec.seScaledToClose = zeros(0, dec.nC);
0026 for iPart = 1:length(partitionData)
0027 %     iPart
0028     % skip unconfined partitions
0029     if ~partitionData(iPart).isConfined,continue;end;
0030     
0031     % get current structuring elements
0032     thisHrepClosed = closeCone(partitionData(iPart).hrep, dec.centerRay);
0033     thisVrepClosed = convertHrepToVrep(thisHrepClosed);
0034     thisVrep = openCone(thisVrepClosed, zerotol);
0035 %     thisSeMatrix = scaleToReal(thisVrep.V, dec, zerotol)
0036         
0037     % filter known SEs
0038     [thisSeMatrix thisSeIndexVector] = eliminateRows(thisVrep.V, ...
0039         dec.seScaledToClose, zerotol);
0040     % append new SEs
0041     thisSeIndexVector = [thisSeIndexVector, size(dec.seScaledToClose, 1) + ...
0042         (1:size(thisSeMatrix, 1))];
0043     dec.seScaledToClose = [dec.seScaledToClose; thisSeMatrix];
0044     % store result
0045     partitionData(iPart).seIndexVector = sort(thisSeIndexVector);
0046 end
0047 % scale all structuring elements to physical
0048 dec.se = scaleToReal(dec.seScaledToClose, dec, zerotol);
0049 
0050 %% recalculate center ray
0051 dec.centerRay = mean(dec.se, 1)';
0052 dec.centerRay = dec.centerRay / norm(dec.centerRay);
0053 % scale SEs
0054 dec.seScaledToClose = scaleToClose(dec.se, dec);
0055 % volume of validity cone is, hence, invalid - instead of recalculating it,
0056 %   it is deleted
0057 dec.validityCone = rmfield(dec.validityCone, 'volume');
0058 % se indices
0059 dec.validityCone.seIndexVector = find(any(abs(dec.validityCone.A * ...
0060     dec.seScaledToClose') < zerotol, 1));
0061 
0062 %% construct confinement cone
0063 % all structuring elements span the confinement cone, even though it is not
0064 %   a proper positive basis.
0065 vrepClosedCone.V = closeCone(dec.seScaledToClose);
0066 hrepClosedCone = convertVrepToHrep(vrepClosedCone, zerotol);
0067 dec.confinementCone = openCone(hrepClosedCone);
0068 % add SE indices
0069 dec.confinementCone.seIndexVector = find(any(abs(dec.confinementCone.A * ...
0070     dec.seScaledToClose') < zerotol, 1));
0071 % add cone volume
0072 dec.confinementCone.volume = computeVolumeVrep(vrepClosedCone);
0073 
0074 %% add unified & simplex partitions
0075 % we will track simultaneaously the total volume of all unified partitions.
0076 %   Which must sum up to the confinement cone volume for non-odd
0077 %   representations
0078 % initialize
0079 dec = appendUnifiedPartition(dec);
0080 unifiedVolumeSum = 0;
0081 % execute
0082 for iPart = 1:length(partitionData)
0083     % skip if not confined
0084     if ~partitionData(iPart).isConfined
0085         continue
0086     end
0087     
0088     % add partition as usual
0089     dec = appendUnifiedPartition(dec, partitionData(iPart).seIndexVector, ...
0090         zerotol, partitionData(iPart).isValid, 1);
0091     
0092     dec = appendPartitionDataForEmbedding(dec, partitionData(iPart), zerotol);
0093     
0094     % accumulate unified partition volume
0095     unifiedVolumeSum = unifiedVolumeSum + ...
0096         dec.unifiedPartition{end}.volumeScaled;
0097     
0098 end
0099 dec.nPU = length(dec.unifiedPartition);
0100 dec.nPS = length(dec.simplexPartition);
0101 dec.nS = size(dec.se, 1);
0102 
0103 %% check oddity
0104 % the confinement cone for odd representations is not convex. Here, the
0105 %   scaled volumes of the unified partitions don't sum up to one.
0106 if abs(unifiedVolumeSum - 1) <= zerotol
0107     dec.isOrdinary = 1;
0108     return
0109 else
0110     dec.isOrdinary = 0;
0111     % if the result is singular, we must add some more partitions
0112 end
0113 
0114 %% add unconfined partitions
0115 for iPart = 1:length(partitionData)
0116     % skip if confined
0117     if partitionData(iPart).isConfined, continue;end;
0118     
0119     % compute SEs (physically scaled)
0120     thisHrepClosed = closeCone(struct(...
0121         'A', [dec.confinementCone.A; partitionData(iPart).hrep.A], ...
0122         'h', [dec.confinementCone.h; partitionData(iPart).hrep.h]), ...
0123         dec.centerRay);
0124     thisVrepClosed = convertHrepToVrep(thisHrepClosed);
0125     thisVrep = openCone(thisVrepClosed, zerotol);
0126     thisSeMatrix = scaleToReal(thisVrep.V, dec, zerotol);
0127     
0128     % as we intersect a partition with the confinement cone the result can
0129     %   be a lower dimensional entity at the boundary of the confinement
0130     %   cone.
0131     if computeVolumeVrep(thisVrepClosed)/dec.confinementCone.volume < zerotol
0132         continue
0133     end
0134     
0135     % We expect that a confined partition is either fully contained in the
0136     %   confinement cone or it is fully outside of the confinement cone.
0137     %   Hence, we check if the extreme rays reside fully outside or fully
0138     %   inside:
0139     distanceMatrix = dec.confinementCone.A * thisSeMatrix';
0140     if all(distanceMatrix(:) > -zerotol)
0141         % partition completely outside
0142         continue
0143     elseif all(distanceMatrix(:) < zerotol)
0144         % partition completely inside, continue below
0145     else
0146         error('elk:decomposition:numericalError', ['I''m confused since a ' ...
0147             'partition should either fully reside inside the confinement ' ...
0148             'cone or fully outside. Most likely, the numerics fail, here. ' ...
0149             'This should not happen.']);
0150     end
0151         
0152     % filter known SEs
0153     [thisSeMatrix thisSeIndexVector] = eliminateRows(thisSeMatrix, ...
0154         dec.se, zerotol);
0155     % append new SEs
0156     thisSeIndexVector = [thisSeIndexVector, size(dec.se, 1) + ...
0157         (1:size(thisSeMatrix, 1))];
0158     dec.se = [dec.se; thisSeMatrix];
0159     dec.seScaledToClose = [dec.seScaledToClose; ...
0160                            scaleToClose(thisSeMatrix, dec)];
0161     
0162     % add partition as usual
0163     dec = appendUnifiedPartition(dec, thisSeIndexVector, zerotol, ...
0164         partitionData(iPart).isValid, 0);
0165     dec = appendPartitionDataForEmbedding(dec, partitionData(iPart), zerotol);
0166       
0167     % accumulate unified partition volume
0168     unifiedVolumeSum = unifiedVolumeSum + ...
0169         dec.unifiedPartition{end}.volumeScaled;
0170       
0171 end
0172 
0173 % recompute some values:
0174 dec.nPU = length(dec.unifiedPartition);
0175 dec.nPS = length(dec.simplexPartition);
0176 dec.nS = size(dec.se, 1);
0177 
0178 % if abs(unifiedVolumeSum - 1) > zerotol
0179 %     error('elk:decomposition:numericalError', ['I''ve considered any ' ...
0180 %         'partition I know and still don''t get the volume of the confinement ' ...
0181 %         'cone. :/ This should not happen.'])
0182 % end