prepareSeeding

 ensure proper seeding structure

 THIS IS NO USER FUNCTION

• solvePbe solve a population balance (growth, nucleation, seeding, bulk states)

0001 function [seeding seedingTimeVector] = prepareSeeding(seeding, pbeDef, ...
0002     optionStruct)
0003 % ensure proper seeding structure
0004 %
0005 % THIS IS NO USER FUNCTION
0006 
0007 % The elk-library: convex geometry applied to crystallization modeling.
0008 %   Copyright (C) 2012 Alexander Reinhold
0009 %
0010 % This program is free software: you can redistribute it and/or modify it
0011 %   under the terms of the GNU General Public License as published by the
0012 %   Free Software Foundation, either version 3 of the License, or (at your
0013 %   option) any later version.
0014 %
0015 % This program is distributed in the hope that it will be useful, but
0016 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0017 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0018 %   General Public License for more details.
0019 %
0020 % You should have received a copy of the GNU General Public License along
0021 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0022 
0023 seedingTimeVector = nan(1, length(seeding));
0024 flagDensity = pbeDef.optionStruct.density;
0025 dec = pbeDef.dec;
0026 
0027 %! Dev: Multiple seed populations cannot be supported since the
0028 %    distribution density might be wrong after a second seeding event. This
0029 %    does also affect the error estimates for integral properties, though,
0030 %    the integrals itself should be properly calculated.
0031 % It should be straight forward to allow multiple seeding in this
0032 %   implementation, but it is not required in this work, so far.
0033 % If you attempt to add seeding, please consider that multiple seeding
0034 %   might imply having multiple crystal types in the system that might use
0035 %   different decomposition data. This change is straight forward, but much
0036 %   more complicated since it affects most of the functions in this folder.
0037 %/
0038 if length(seeding) > 1
0039     error('elk:pbeSolver:wrongInput', ['Handling of multiple seed ' ...
0040         'populations is currently not properly supported. View notes in ' ...
0041         'this function for details']);
0042 end
0043 %/
0044 
0045 for iSeed = 1:length(seeding)
0046     %% seeding time
0047     if isfield(seeding, 'time') && ~isempty(seeding(iSeed).time)
0048         seedingTimeVector(iSeed) = seeding(iSeed).time;
0049     elseif iSeed == 1
0050         seedingTimeVector(iSeed) = 0;
0051     else
0052         error('elk:pbeSolver:wrongInput', ['You need to provide the seeding ' ...
0053             'time with pbeDefinition.seeding(iSeed).time.']);
0054     end
0055     
0056     % check case when user provides disC or fhandle and disD
0057     if isfield(seeding, 'disD') && ...
0058             (isfield(seeding, 'disC') || isfield(seeding, 'fhandle'))
0059         seeding(iSeed).positionMatrix = seeding(iSeed).disD.positionMatrix;
0060         seeding(iSeed).probabilityVector = seeding(iSeed).disD.probabilityVector;
0061         if isfield(seeding, 'fhandle')
0062             thisHandle = seeding(iSeed).fhandle;
0063         else
0064             thisHandle = seeding(iSeed).disC.fhandle;
0065         end
0066         seeding(iSeed).numberVector = ...
0067             1/length(seeding(iSeed).probabilityVector) * ...
0068             thisHandle(seeding(iSeed).positionMatrix) ./ ...
0069             seeding(iSeed).probabilityVector;
0070     end
0071     
0072     % check probability vector
0073     if ~isfield(seeding, 'probabilityVector') || ...
0074             isempty(seeding(iSeed).probabilityVector')       
0075         % throw error if density should have been provided
0076         if flagDensity
0077             error('elk:pbeSolver:wrongInput', ['To provide ' ...
0078              'distribution density information, you need to provide the ' ...
0079              'probability density of selecting the pivots in seeding ' ...
0080              'population ' num2str(iSeed)]);
0081         end;
0082         % Dev: we change the probability vector to nan to ensure it cannot
0083         %   be used
0084         seeding(iSeed).probabilityVector = ...
0085             seeding(iSeed).numberVector*0 + nan;
0086     end
0087     
0088     %% validate positionMatrix, numberVector and probabilityVector
0089     % enforce proper sizing
0090     seeding(iSeed).numberVector = seeding(iSeed).numberVector(:)';
0091     seeding(iSeed).probabilityVector = seeding(iSeed).probabilityVector(:)';
0092     % check dimension
0093     if size(seeding(iSeed).positionMatrix, 1) - pbeDef.nExtraCoordinate ...
0094             ~= pbeDef.dec.nC
0095         error('elk:populationBalance:wrongInput', ['Number of rows in ' ...
0096             ' positionMatrix of seeding(' num2str(iSeed) ') does not ' ...
0097             'match the dimension of the problem.']);
0098     end
0099 
0100     nSample = size(seeding(iSeed).positionMatrix, 2);
0101     if length(seeding(iSeed).numberVector) ~= nSample || ...
0102             length(seeding(iSeed).probabilityVector) ~= nSample
0103         error('elk:populationBalance:wrongInput', ['Elements in in ' ...
0104             ' numberVector or probabilityVector of seeding(' ...
0105              num2str(iSeed) ') does not match the number of sample points']);
0106     end
0107            
0108     %% ensure positions are valid (ignored for incompletes)
0109     if dec.isComplete
0110     toleranceValue = optionStruct.distanceLimit;
0111     if (dec.isOrdinary && ...
0112         any(any(dec.confinementCone.A * seeding(iSeed).positionMatrix(1:dec.nC, :) > toleranceValue)) ...
0113        ) || (~dec.isOrdinary && ...
0114         any(any(dec.confinementConeHull.A * seeding(iSeed).positionMatrix(1:dec.nC, :) > toleranceValue)) ...
0115        )
0116         error('elk:populationBalance:wrongInput', ...
0117             ['Not all pivots of seeding population ' num2str(iSeed) ...
0118              ' reside in the confinement cone']);
0119     end
0120     end
0121     
0122 end