ensureProperDecData

 check if unified partitions are right and correct of required

 Note that the resulting hcMatrix is either formulated in the embedding 
   proper representation (option useEmbedding

 THIS IS NO USER FUNCTION

• computePartitionIndexDec determine the partition of a given hC-vector
• mapToConfinementCone map a hC vector to the cone of confined hC-vectors
• computePositionTolVector compute tolerance vector for individual vecctors

• addNucleation add nucleation handling to solver structure
• genericDerivative compute derivative of states from structures
• newBinForNucleation remove nucleation handling to solver structure
• solvePbe solve a population balance (growth, nucleation, seeding, bulk states)

0001 function [hcMatrix, upVector, fractionCorrected] = ensureProperDecData(...
0002     hcMatrixDynamic, upVector, pbeDef)
0003 % check if unified partitions are right and correct of required
0004 %
0005 % Note that the resulting hcMatrix is either formulated in the embedding
0006 %   proper representation (option useEmbedding
0007 %
0008 % THIS IS NO USER FUNCTION
0009 
0010 % The elk-library: convex geometry applied to crystallization modeling.
0011 %   Copyright (C) 2013 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 %% shortcuts
0027 nDecCoordinate = pbeDef.nDecCoordinate;
0028 rowFilter = false(1, size(hcMatrixDynamic, 1));
0029 rowFilter(1:pbeDef.nDecCoordinate) = true;
0030 % empty hcMatrix
0031 if isempty(hcMatrixDynamic)
0032     fractionCorrected = 0;
0033     return
0034 end
0035 % incomplete decompositions
0036 if ~pbeDef.dec.isComplete
0037     fractionCorrected = 0;
0038     return
0039 end
0040 
0041 %% some preparation
0042 % tolerances
0043 tolVector = computePositionTolVector(hcMatrixDynamic(1:nDecCoordinate, :), pbeDef.optionStruct);
0044 % the correct filter to allow dynamic update of unified partition
0045 %   information in solver states (this computation here is made freshly for
0046 %   each time step.. ..the unified parition is not a state variable for the
0047 %   ODE solver)
0048 correctFilter = false(size(upVector));
0049 % if the unified partition is still OK, the vector is also still
0050 %   valid/confined.
0051 hcMatrix = hcMatrixDynamic;
0052 
0053 %% shortcut for incomplete representations
0054 if ~pbeDef.dec.isComplete
0055     % in this case, there is only one unified partition
0056     upVector = ones(1, size(hcMatrix, 2));
0057     
0058     % .. but the hcMatrix might be required in the (extended) embedding:
0059     if pbeDef.optionStruct.useEmbeddingPartitions
0060         hcMatrix(rowFilter, :) = pbeDef.dec.properEmbedding.mappingNewToProper * ...
0061             hcMatrix(rowFilter, :);
0062     end
0063 end
0064 
0065 % distinguish cases which use the embedding and which use the original
0066 if (pbeDef.dec.isProper && pbeDef.dec.isComplete) || ...
0067         ~pbeDef.optionStruct.useEmbeddingPartitions
0068     
0069     % normal case of a proper representation
0070     
0071     % if the growth rate is not adapted, validity must be ensured
0072     %   before the unified partition is checked
0073     if ~pbeDef.optionStruct.limitGrowthRate
0074         hcMatrix(rowFilter, :) = mapToConfinementCone(...
0075             hcMatrix(rowFilter, :), pbeDef.dec);
0076     end
0077         
0078         
0079     for upIndex = 1:pbeDef.dec.nPU
0080         % we handle only pivots of the current unified partition
0081         thisFilter = upVector == upIndex;
0082         % skip if nothing to do
0083         if ~any(thisFilter), continue;end;
0084         
0085         % are they in - we filter here for true, when the upVector needs
0086         %   correction
0087         thisFilter(thisFilter) = max(...
0088             pbeDef.dec.unifiedPartition{upIndex}.A * ...
0089             hcMatrix(1:pbeDef.nDecCoordinate, thisFilter), ...
0090             [], 1) > tolVector(thisFilter);
0091         
0092         % correct if required
0093         if any(thisFilter)
0094             % only if growth rate limiting is active, validity must be
0095             % corrected here. The other case is handled above.
0096             if pbeDef.optionStruct.limitGrowthRate
0097                 hcMatrix(rowFilter, thisFilter) = mapToConfinementCone(...
0098                     hcMatrix(rowFilter, thisFilter), pbeDef.dec);
0099             end    
0100             
0101             upVector(thisFilter) = computePartitionIndexDec(...
0102                 pbeDef.dec, hcMatrix(rowFilter, thisFilter), ...
0103                 tolVector(thisFilter), 1);
0104 
0105             correctFilter = correctFilter | thisFilter;
0106         end
0107         
0108     end
0109 elseif pbeDef.optionStruct.useEmbeddingPartitions
0110     % create a hcMatrix that is already in embedding but without extra
0111     %   coordinates
0112     hcMatrixLocal = pbeDef.dec.properData.mappingNewToProper * ...
0113         hcMatrix(rowFilter, :);
0114     % again, it only makes sense to check valid vectors
0115     hcMatrixLocal = mapToConfinementCone(...
0116         hcMatrixLocal, pbeDef.dec.properData);
0117     
0118     % cycle through partitions in the embedding
0119     for upIndex = 1:pbeDef.dec.properData.nPU
0120         % we handle only pivots of the current unified partition
0121         thisFilter = upVector == upIndex;
0122         % skip if nothing to do
0123         if ~any(thisFilter), continue;end;
0124         
0125         % are they in - we filter here for true, when the upVector needs
0126         %   correction
0127         thisFilter(thisFilter) = max(...
0128             pbeDef.dec.properData.unifiedPartition{upIndex}.A * ...
0129             hcMatrixLocal(:, thisFilter), ...
0130             [], 1) > tolVector(thisFilter);
0131         
0132         % correct if required
0133         if any(thisFilter)
0134             hcMatrixLocal(:, thisFilter) = mapToConfinementCone(...
0135                 hcMatrixLocal(:, thisFilter), pbeDef.dec.properData);
0136             
0137             upVector(thisFilter) = computePartitionIndexDec(...
0138                 pbeDef.dec.properData, hcMatrixLocal(:, thisFilter), ...
0139                 tolVector(thisFilter), 1);
0140             
0141             correctFilter = correctFilter | thisFilter;
0142         end
0143         
0144     end 
0145     
0146     % but we might need to add the additional coordinates again
0147     hcMatrix = [hcMatrixLocal; hcMatrix(~rowFilter, :)];
0148 else
0149     error('elk:pbeSolver:inputError', ['The given representation is neither ' ...
0150         'proper, nor improper but complete, nor incomplete. This does not ' ...
0151         'makes sense. Please check your decomposition, dude.']);
0152 end
0153 
0154 
0155 if any(correctFilter)
0156     fractionCorrected = sum(correctFilter) / length(correctFilter);
0157 else
0158     fractionCorrected = 0;
0159 end