% Script for 1st LPT and DA Challenge output writing, DA cases

% Mesh extent and number of points in each direction
% Example of a domain covering [-50,50]x[-25,25]x[0.01,30.01] (X Y Z), with 251x126x76 points

Xlin = linspace(-50,50,251);
Ylin = linspace(-25,25,126);
Zlin = linspace(0.01,30.01,76);

% 3D mesh: first index describes variation in X, second in Y, third in Z

[Y,X,Z] = meshgrid(Ylin,Xlin,Zlin);

VX = zeros(size(X));
VY = zeros(size(X));
VZ = zeros(size(X));
dVXdX = zeros(size(X));
dVXdY = zeros(size(X));
dVXdZ = zeros(size(X));
dVYdX = zeros(size(X));
dVYdY = zeros(size(X));
dVYdZ = zeros(size(X));
dVZdX = zeros(size(X));
dVZdY = zeros(size(X));
dVZdZ = zeros(size(X));
p = zeros(size(X));
T = zeros(size(X));

%If not in that format already, write/interpolate your results into the 3D numpy arrays for velocity, velocity gradient and pressure
%X = ...
%...
%VX = ...
%...
%p = ...
%T = ...

% Pressure normalization: with respect to a specified point, in the present example point (X,Y,Z)=(0,0.2,0.01)

p = p-p(126,64,1);

% Output file naming
Team = 'ONERA';
% In case you upload two results obtained with two different algorithms, then use different team names, e.g. ONERA_XXXX and ONERA_YYYY
Case = 1;
ImageStrategy = "TR";
Ppp = 0.16;
InstantNumber = 25;

if Case == 1
    Frequency = 1 # only for case 2
    InertialPart = 0 # only for case 2
    resultFileTemplate = '%s_DA_CASE%02d_%s_ppp_0_%03d_FR_%d_I%d_PartFile_%04d.h5'
    resultFile = sprintf(resultFileTemplate,Team,Case,ImageStrategy,round(Ppp*1000),Frequency,InertialPart,InstantNumber)
else
    resultFileTemplate = '%s_DA_CASE%02d_%s_ppp_0_%03d_PartFile_%04d.h5'
    resultFile = sprintf(resultFileTemplate,Team,Case,ImageStrategy,round(Ppp*1000),Frequency,InertialPart,InstantNumber)
end

% At file writing, convert the 3D arrays to a 1D linear array with ordering as requested in the instructions.
% 
% Note that VXout = VX(:) is equivalent to e.g.
% 
% ii=1
% for k=1:length(Zlin)
%     for j=1:length(Ylin)
%         for i=1:length(Xlin)
%             VXout(ii) = VX(i,j,k)
%             ii = ii+1
%          end
%     end
% end

h5create(resultFile,'/X',length(X(:)));
h5write(resultFile,'/X',X(:));
h5create(resultFile,'/Y',length(Y(:)));
h5write(resultFile,'/Y',Y(:));
h5create(resultFile,'/Z',length(Z(:)));
h5write(resultFile,'/Z',Z(:));

h5create(resultFile,'/VX',length(VX(:)));
h5write(resultFile,'/VX',VX(:));
h5create(resultFile,'/VY',length(VY(:)));
h5write(resultFile,'/VY',VY(:));
h5create(resultFile,'/VZ',length(VZ(:)));
h5write(resultFile,'/VZ',VZ(:));

h5create(resultFile,'/dVXdX',length(dVXdX(:)));
h5write(resultFile,'/dVXdX',dVXdX(:));
h5create(resultFile,'/dVYdX',length(dVYdX(:)));
h5write(resultFile,'/dVYdX',dVYdX(:));
h5create(resultFile,'/dVZdX',length(dVZdX(:)));
h5write(resultFile,'/dVZdX',dVZdX(:));

h5create(resultFile,'/dVXdY',length(dVXdY(:)));
h5write(resultFile,'/dVXdY',dVXdY(:));
h5create(resultFile,'/dVYdY',length(dVYdY(:)));
h5write(resultFile,'/dVYdY',dVYdY(:));
h5create(resultFile,'/dVZdY',length(dVZdY(:)));
h5write(resultFile,'/dVZdY',dVZdY(:));

h5create(resultFile,'/dVXdZ',length(dVXdZ(:)));
h5write(resultFile,'/dVXdZ',dVXdZ(:));
h5create(resultFile,'/dVYdZ',length(dVYdZ(:)));
h5write(resultFile,'/dVYdZ',dVYdZ(:));
h5create(resultFile,'/dVZdZ',length(dVZdZ(:)));
h5write(resultFile,'/dVZdZ',dVZdZ(:));

h5create(resultFile,'/p',length(p(:)));
h5write(resultFile,'/p',p(:));

h5create(resultFile,'/T',length(T(:)));
h5write(resultFile,'/T',T(:));

% display contents of created file
h5disp(resultFile);

