function pbeDefinition = snyderCuboidGrowthSingle
clear all
%% Basic information
pbeDefinition.info = ['seeded dissolution of cuboids with new facets ' ...
    'evolving on edges (steppef faces) and vertices (kinked)'];

fileNameDecomposition = 'snyderCuboidDissolution_dec.matx';
if exist(fileNameDecomposition, 'file')
    dec = load(fileNameDecomposition, '-mat');
    dec = dec.dec;
else
    % Load case (contains source information)
    dec = obtainCrystalDec('snyder_cuboid');
    % generate measure information
    dec = addMeasureDataDec(dec, 'volume');
    dec = addMeasureDataDec(dec, 'surface');
    dec = addMeasureDataDec(dec, 'meanWidth');
    dec = addMeasureDataDec(dec, 'support', dec.A([1 7 11 15 19], :));
    % save file
    save(fileNameDecomposition, '-mat', 'dec');
end
% copy to pbeDef
pbeDefinition.source = dec.cdef.source;
pbeDefinition.dec = dec;

%% initial condition
% NORMAL distribution for cuboid facets. The other facets are exectly
%   disappeared, giving some dirac functions that would be multiplied.
pbeDefinition.seeding.mean = [1 2 3]*500*1e-6;
pbeDefinition.seeding.std = [0.5 1 2]'*50*1e-6;
pbeDefinition.seeding.cov = diag(pbeDefinition.seeding.std.^2);
pbeDefinition.seeding.totalValue = 1e10;
pbeDefinition.seeding.nSample = 100;

% we sample the cuboid coordinates directly..
rng(0)
disD = obtainDistributionDiscrete('normal', ...
    'nSample',pbeDefinition.seeding.nSample, ...
    'mean', pbeDefinition.seeding.mean, ...
    'cov', pbeDefinition.seeding.cov);
disD.positionMatrix = disD.positionMatrix(:,1:6);
% disD.positionMatrix = [1 2 3]'*500*1e-6;
% check standard deviation
% checkStd = std(disD.positionMatrix, [], 2)'
% ..construct the mathing additional coordinates..
additionalCoordinateMatrix = [zeros(4,3) eye(4)]*...
    mapToValidityCone([disD.positionMatrix; ...
                       100*ones(4,1)*max(disD.positionMatrix,[], 1)], ...
                       dec, 1e-12);
% ..add a small numerical correction (validity mapping is not perfect)
additionalCoordinateMatrix = bsxfun(@minus, additionalCoordinateMatrix, ...
    [1 1 1 2]'*1e-15);
% ..and finally add the information for the initial distribution:
pbeDefinition.seeding.positionMatrix = [disD.positionMatrix; ...
    additionalCoordinateMatrix];
pbeDefinition.seeding.numberVector = disD.positionMatrix(1,:)*0 + ....
    pbeDefinition.seeding.totalValue/pbeDefinition.seeding.nSample;
% ..correct validity NOW .. .. this would otherwise happen in the pbe
% solver
pbeDefinition.seeding.positionMatrix = mapToConfinementCone(...
    pbeDefinition.seeding.positionMatrix, pbeDefinition.dec, 0);

%% growth and analytical solution
% the relative growth rates are already specified in the crystal definition
%   where two versions are available. In the first organic variant 
%   (cdef.relativeGrowthOrganic), all edge and vertex facets dissolve
%   equally. In case two (cdef.relativeGrowthCaseTwo) the vertex facet 
%   dissolves faster.
pbeDefinition.growthRate = @(H, P, y, t)(...
    1./dec.cdef.relativeGrowthCaseTwo*1e-6*...
    ones(1, size(H,2))...
    );

%% pivot properties
% we need the unified partitions in P for growth rate projections
pbeDefinition.pivotProperty = @(H, up)(up);

%% time span for solution
pbeDefinition.timeSpan = 0:60:(1*60*60); %s

% options
pbeDefinition.optionStruct.distanceLimit = 1e-6; %m
pbeDefinition.optionStruct.relTol = 1e-5;
pbeDefinition.optionStruct.absTolNumber = 1e-30;
pbeDefinition.optionStruct.absTolPosition = 1e-30;
% pbeDefinition.optionStruct.odeSolver = @ode113; % 113, 45