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

% obtain decomposition
dec = obtainSynderCuboidDec;
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 = 1;
pbeDefinition.seeding.nSample = 100;

% we sample the cuboid coordinates directly..
rng(1);
disD = obtainDistributionDiscrete('normal', ...
    'nSample',pbeDefinition.seeding.nSample, ...
    'mean', pbeDefinition.seeding.mean, ...
    'cov', pbeDefinition.seeding.cov);
% 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;

%% 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.growthRateOriginal = @(H, P, y, t)(-1*...
    dec.cdef.relativeGrowthCaseTwo*2.5e-8*...
    ones(1, size(H,2))...
    );
pbeDefinition.growthRate = @(H, P, y, t)(...
    pbeDefinition.growthRateOriginal(H, P, y, t));

%% 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-10;
% pbeDefinition.optionStruct.odeSolver = @ode113; % 113, 45