function [Objectives, motorProps] = evalDesigns_training(motorProps,settingsMagnetostatic,file,maxNumberOfProcesses,designParams,msg0)
% Example of designFunction for axial flux machine designs evaluation used by
% the optimization algorithm (see section 6.3 of User Manual - Part II)
tic
designParams = round(designParams*10000)/10000;   % use only 4 digits after point
nDesigns = size(designParams,1);
motorProps=setParamafm(motorProps,'Wire size method','Fill factor');
motorProps=setParamafm(motorProps,'Coil fill factor',0.35);            % keep the same coil fill factor for all the designs
motorProps=setParamafm(motorProps,'Number Of Parallel Paths',1);           % set number of parallel paths to 1
motorProps=setParamafm(motorProps,'Number Of Turns',1);                    % set number of turns to 1 (winding losses do not depend on number of turns)
motorProps=setParamafm(motorProps,'End Winding Inductance','Auto');
motorProps=setParamafm(motorProps,'Phase Resistance','Auto');
if nargout==2
    % return only "motorProps" structure for saving to mxa-file from the Pareto plot
    motorProps = getMotorProps(motorProps,designParams,1);
    Objectives = [];
    return
end
%%% assemble all AFM designs first
cell_Geometry_ready = cell(nDesigns,1);
cell_Windings_ready = cell(nDesigns,1);
cell_Mesh_ready = cell(nDesigns,1);
cell_Materials_ready = cell(nDesigns,1);
cell_RadMesh_ready = cell(nDesigns,1);
cell_SubdomainProperty_ready = cell(nDesigns,1);
cell_Error_ready = cell(nDesigns,1);
cell_motorProps_output_ready = cell(nDesigns,1);
cell_Weight_ready = cell(nDesigns,1);
numberOfProcesses = 0; % current number of processes
i = 1;  % design number
% tic
initAssembleMXA_par();  % init AFM design assembling in parallel
msg = [msg0 ' assembling designs...)']; disp(msg);
% keep number of processes equal to maxNumberOfProcesses
while 1
    if numberOfProcesses<maxNumberOfProcesses
        motorProps = getMotorProps(motorProps,designParams,i);
        fprintf('%s assemble AFM #%i\n', msg, i);
        % assemble AFM design in parallel - add one process
        runAssembleMXA_par(file, motorProps);
        numberOfProcesses = numberOfProcesses+1;
        i = i+1;
    end
    if numberOfProcesses==maxNumberOfProcesses
        % wait when one of processes finish assembling AFM design
        while 1
            % check if any AFM design(s) are already assembled
            [cell_Geometry, cell_Windings, cell_Mesh, cell_Materials, cell_RadMesh, cell_SubdomainProperty,...
                cell_Error, cell_motorProps_output, cell_Weight, numberOfProcesses] = ...
                    getAssembleMXA_par();   % get assembled AFM design(s)
            if numberOfProcesses<maxNumberOfProcesses
                break;
            end
            pause(1); % wait for 1 second
        end
        % put assembled AFM design(s)
        for n = 1:length(cell_motorProps_output)
            if ~isempty(cell_motorProps_output{n})
                cell_Geometry_ready{n} = cell_Geometry{n};
                cell_Windings_ready{n} = cell_Windings{n};
                cell_Mesh_ready{n} = cell_Mesh{n};
                cell_Materials_ready{n} = cell_Materials{n};
                cell_RadMesh_ready{n} = cell_RadMesh{n};
                cell_SubdomainProperty_ready{n} = cell_SubdomainProperty{n};
                cell_Error_ready{n} = cell_Error{n};
                cell_motorProps_output_ready{n} = cell_motorProps_output{n};
                cell_Weight_ready{n} = cell_Weight{n};
            end
        end
    end
    if i>nDesigns
        break
    end
end

% get all (the rest) AFM designs
numberOfProcesses0 = numberOfProcesses;
while numberOfProcesses>0
    % check if any AFM design(s) are finished assembling
    [cell_Geometry, cell_Windings, cell_Mesh, cell_Materials, cell_RadMesh, cell_SubdomainProperty,...
        cell_Error, cell_motorProps_output, cell_Weight, numberOfProcesses] = ...
        	getAssembleMXA_par();   % get assembled AFM design(s)
    if numberOfProcesses~=numberOfProcesses0  % if new AFM design(s) are finished assembling (number of processes changed)
        % put assembled AFM design(s)
        for n = 1:length(cell_motorProps_output)
            if ~isempty(cell_motorProps_output{n})
                cell_Geometry_ready{n} = cell_Geometry{n};
                cell_Windings_ready{n} = cell_Windings{n};
                cell_Mesh_ready{n} = cell_Mesh{n};
                cell_Materials_ready{n} = cell_Materials{n};
                cell_RadMesh_ready{n} = cell_RadMesh{n};
                cell_SubdomainProperty_ready{n} = cell_SubdomainProperty{n};
                cell_Error_ready{n} = cell_Error{n};
                cell_motorProps_output_ready{n} = cell_motorProps_output{n};
                cell_Weight_ready{n} = cell_Weight{n};
            end
        end
        numberOfProcesses0 = numberOfProcesses;
    end
    pause(1); % wait for 1 second
end
% check if any design failed assembling and try again
for i = 1:nDesigns
    if isempty(cell_Geometry_ready{i})
        motorProps = getMotorProps(motorProps,designParams,i);
        [Geometry, Windings, Mesh, Materials, RadMesh, SubdomainProperty, Error, motorProps_output, Weight] = assembleMXA(file,motorProps);
        if isempty(Geometry)
            cell_Error_ready{i} = 'Unknown error';
        else
            cell_Geometry_ready{i} = Geometry;
            cell_Windings_ready{i} = Windings;
            cell_Mesh_ready{i} = Mesh;
            cell_Materials_ready{i} = Materials;
            cell_RadMesh_ready{i} = RadMesh;
            cell_SubdomainProperty_ready{i} = SubdomainProperty;
            cell_Error_ready{i} = Error;
            cell_motorProps_output_ready{i} = motorProps_output;
            cell_Weight_ready{i} = Weight;
        end
    end
end
% toc
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Save each step solution',0);
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Compute cogging torque',0);
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Convergence tolerance',0.01);
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Current waveform','Sinusoidal');
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Mechanical speed',1600);           % RPM
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Advance angle',0); 
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Simulation time','One electrical period');
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'RMS supply current',2000);     % initial current value
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Number of points',1);
targetTorque = 50000/(1600*2*pi/60);   % N*m  
initCurrent = 2000;
arr_Current = initCurrent*ones(nDesigns,1);
arr_Torque = zeros(nDesigns,1);
enableRun = ones(nDesigns,1);
for n = 1:nDesigns
    if ~isempty(cell_Error_ready{n})
        enableRun(n) = 0;
    end
end

% Determine current required to produce target torque for each design
msg = [msg0 ' calculating stator currents...)']; disp(msg);
iter = 1;
cell_Private_ready = cell(nDesigns,1);
cell_Coefs3D_ready = cell(nDesigns,1);
while 1
    [cell_Results, cell_Timesteppingdata, cell_Private, cell_Coefs3D] = ...
        runMStimestepping_batch(cell_Private_ready,cell_Geometry_ready,cell_Windings_ready,cell_Mesh_ready,cell_Materials_ready,cell_RadMesh_ready,cell_SubdomainProperty_ready,cell_Error_ready,...
                                [],settingsMagnetostatic,arr_Current,enableRun,maxNumberOfProcesses,msg);
    for n = 1:length(cell_Results)
        if enableRun(n)
            if isempty(cell_Results{n})
                error(['Simulation results for design #' num2str(n) ' are empty (1)']);
            end
            arr_Torque(n) = cell_Results{n}.torque;
            if 100*abs((arr_Torque(n)-targetTorque)/targetTorque)<0.5          % 0.5% torque error tolerance
                enableRun(n) = 0;
            end
            cell_Coefs3D_ready{n} = cell_Coefs3D{n};      % update the Coefs3D structure for the renewed current value
            if isempty(cell_Private_ready{n})
                cell_Private_ready{n} = cell_Private{n};
            end
        end
    end
    if any(enableRun)
        arr_Current = (targetTorque./arr_Torque).*arr_Current;
    else
        break
    end
    iter = iter+1;
    if iter>7
        break
    end
end
% run magnetostatic simulation for 6 points and 1/6th of electrical period
speed = getParamafm(settingsMagnetostatic,'Mechanical speed')*2*pi/60;    % rad/s
Pmech_target = abs(speed*targetTorque);
f1 = speed*getParamafm(motorProps,'Number of pole pairs')/(2*pi);
T1 = 1/f1;

settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Simulation time','One electrical period');
nPolePairs = getParamafm(motorProps,'Number of pole pairs');
nSlots = getParamafm(motorProps,'Number of slots');
settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Number of points',setNumberOfPoints('Low accuracy', nPolePairs, nSlots));
% settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Time from',0);
% settingsMagnetostatic=setParamafm(settingsMagnetostatic,'Time to',T1/6);    % 1/6th of electrical period
enableRun = ones(nDesigns,1);
for n = 1:nDesigns
    if ~isempty(cell_Error_ready{n})
        enableRun(n) = 0;
    end
end
msg = [msg0 ' running simulations...)']; disp(msg);
[cell_Results, cell_Timesteppingdata, cell_Private, cell_Coefs3D] = ...
        runMStimestepping_batch(cell_Private_ready,cell_Geometry_ready,cell_Windings_ready,cell_Mesh_ready,cell_Materials_ready,cell_RadMesh_ready,cell_SubdomainProperty_ready,cell_Error_ready,...
                                cell_Coefs3D_ready,settingsMagnetostatic,arr_Current,enableRun,maxNumberOfProcesses,msg);
motorWeight = inf*ones(nDesigns,1);
motorLoss = inf*ones(nDesigns,1); 
motorLoss_norm = inf*ones(nDesigns,1);   % normalized to the same mechanical power
magnetLoss = inf*ones(nDesigns,1);
motorTorque = inf*ones(nDesigns,1); 
motorEfficiency = inf*ones(nDesigns,1); 
for n=1:nDesigns  
    Results = cell_Results{n};
    if enableRun(n)
        if isempty(Results)
            error(['Simulation results for design #' num2str(n) ' are empty (2)'])
        end
        Weight = cell_Weight_ready{n};
        motorWeight(n) = Weight.statoriron + Weight.rotoriron + Weight.statorwinding + Weight.magnet + Weight.conductor;
        motorLoss(n) = Results.Ps + Results.Piron_hyst + Results.Piron_eddy + Results.MagnetLoss + Results.OtherECloss;
        motorTorque(n) = Results.torque;
        motorEfficiency(n) = abs(Results.Pmech)/(abs(Results.Pmech)+motorLoss(n));
        motorLoss_norm(n) = Pmech_target*(1-motorEfficiency(n))/motorEfficiency(n);
        magnetLoss(n) = Results.MagnetLoss;
        if motorLoss_norm(n)>6000
            % punish algorithm for losses > 3000W  
            motorWeight(n) = 50;
        end
    end
end
Objectives = [1000*motorWeight motorLoss_norm magnetLoss];
toc
if any(isnan(Objectives))
    Stop=1;
end
end

%________________________________________________________________________________________________________________________


function [cell_Results_ready, cell_Timesteppingdata_ready, cell_Private_output, cell_Coefs3D_output] = ...
    runMStimestepping_batch(cell_Private_input,cell_Geometry,cell_Windings,cell_Mesh,cell_Materials,cell_RadMesh,cell_SubdomainProperty,cell_Error,cell_Coefs3D_input,...
                            settingsMagnetostatic,arr_Current,enableRun,maxNumberOfProcesses,msg)
%%% run magnetostatic FE simulations for "nDesigns" designs in parallel
nDesigns = length(arr_Current);
cell_Results_ready = cell(nDesigns,1);
cell_Timesteppingdata_ready = cell(nDesigns,1);
cell_Private_output = cell(nDesigns,1);
cell_Coefs3D_output = cell(nDesigns,1);
numberOfProcesses = 0; % current number of processes
i = 1;
initMStimesteppingafm_par(); % init running magnetostatic FE simulations in parallel
% keep number of processes equal to maxNumberOfProcesses
while 1
    if numberOfProcesses<maxNumberOfProcesses
        Geometry = cell_Geometry{i};
        Windings = cell_Windings{i};
        Mesh = cell_Mesh{i};
        Materials = cell_Materials{i};
        RadMesh = cell_RadMesh{i};
        SubdomainProperty = cell_SubdomainProperty{i};
        Error = cell_Error{i};
        if ~isempty(cell_Private_input)
            Private = cell_Private_input{i};
        else
            Private = [];
        end
        if ~isempty(cell_Coefs3D_input)
            Coefs3D = cell_Coefs3D_input{i};
        else
            Coefs3D = [];
        end
        if isempty(Error) && enableRun(i) % && cmpstruct(motorProps,motorProps_output)
            settingsMagnetostatic=setParamafm(settingsMagnetostatic,'RMS supply current',arr_Current(i));
            fprintf('%s run simulation #%i\n', msg, i);
            % run magnetostatic FE simulation in parallel - add one process
            runMStimesteppingafm_par(Private,Geometry,Mesh,Windings,settingsMagnetostatic,RadMesh,Materials,SubdomainProperty,[],Coefs3D);
            numberOfProcesses = numberOfProcesses+1;
        end
        i = i+1;
    end
    if numberOfProcesses == maxNumberOfProcesses
        % wait when one of processes finish calculation
        while 1
            % check if any results are ready
            [cell_Results, cell_Timesteppingdata, cell_Coefs3D, ~, cell_Private, numberOfProcesses] = getMStimesteppingafm_par();   % get magnetostatic analysis results
            if numberOfProcesses<maxNumberOfProcesses
                break;
            end
            pause(1); % wait for 1 second
        end
        % put results
        k = 1;
        for n = 1:length(cell_Results)
            while enableRun(k) == 0
                k = k+1;
            end
            if ~isempty(cell_Results{n})
                cell_Results_ready{k} = cell_Results{n};
                cell_Timesteppingdata_ready{k} = cell_Timesteppingdata{n};
                cell_Private_output{k} = cell_Private{n};
                cell_Coefs3D_output{k} = cell_Coefs3D{n};
            end
            k = k+1;
        end
    end
    if i>nDesigns
        break
    end
end

% get all (the rest) results
numberOfProcesses0 = numberOfProcesses;
while numberOfProcesses>0
    % check if any results are ready
    [cell_Results, cell_Timesteppingdata, cell_Coefs3D, ~, cell_Private, numberOfProcesses] = getMStimesteppingafm_par();
    if numberOfProcesses ~= numberOfProcesses0  % if new results are ready (number of processes changed)
        % put results
        k = 1;
        for n = 1:length(cell_Results)
            while enableRun(k) == 0
                k = k+1;
            end
            if ~isempty(cell_Results{n})
                cell_Results_ready{k} = cell_Results{n};
                cell_Timesteppingdata_ready{k} = cell_Timesteppingdata{n};
                cell_Private_output{k} = cell_Private{n};
                cell_Coefs3D_output{k} = cell_Coefs3D{n};
            end
            k = k+1;
        end
        numberOfProcesses0 = numberOfProcesses;
    end
    pause(1); % wait for 1 second
end
% check if any simulation failed and try again
for i = 1:nDesigns  
    Error = cell_Error{i};
    if isempty(Error) && enableRun(i)
        if isempty(cell_Results_ready{i})
            Geometry = cell_Geometry{i};
            Windings = cell_Windings{i};
            Mesh = cell_Mesh{i};
            Materials = cell_Materials{i};
            RadMesh = cell_RadMesh{i};
            SubdomainProperty = cell_SubdomainProperty{i};
            settingsMagnetostatic=setParamafm(settingsMagnetostatic,'RMS supply current',arr_Current(i));
            try
                [Results, Timesteppingdata] = runMStimesteppingafm([],Geometry,Mesh,Windings,settingsMagnetostatic,RadMesh,Materials,SubdomainProperty,0); % run magnetostatic analysis
            catch ME
                rethrow(ME)
            end
            cell_Results_ready{i} = Results;
            cell_Timesteppingdata_ready{i} = Timesteppingdata;
        end
    end
end
        
end

%________________________________________________________________________________________________________________________


function motorProps = getMotorProps(motorProps,designParams,i)
innerDiameter = designParams(i,1);                     % motor inner diameter
slotWidth = designParams(i,2);                         % slot width
statorHeight = designParams(i,3);                      % stator height
magnetSpacing = designParams(i,4);                     % magnet spacing
yokeHeight = designParams(i,5);                        % yoke height
magnetHeight = designParams(i,6);                      % magnet height
airGap = designParams(i,7);                            % air gap
outerDiameter = 350;                                   % motor outer diameter
innerDiameter = checkInnerDiameter(innerDiameter,slotWidth,getParamafm(motorProps,'Number of slots'));
innerDiameter = round(innerDiameter*10000)/10000;      % use only 4 digits after point
motorProps=setParamafm(motorProps,'Stator inner diameter',innerDiameter);            % stator inner diameter
motorProps=setParamafm(motorProps,'Stator outer diameter',outerDiameter);            % stator outer diameter
motorProps=setParamafm(motorProps,'Stator tooth inner diameter',innerDiameter);      % tooth inner diameter
motorProps=setParamafm(motorProps,'Stator tooth outer diameter',outerDiameter);      % tooth outer diameter
motorProps=setParamafm(motorProps,'Slot width',slotWidth);                           % slot width
motorProps=setParamafm(motorProps,'Stator height',statorHeight);                     % stator height
motorProps=setParamafm(motorProps,'Rotor inner diameter',innerDiameter);             % rotor inner diameter
motorProps=setParamafm(motorProps,'Rotor outer diameter',outerDiameter);             % rotor outer diameter
motorProps=setParamafm(motorProps,'Magnet inner diameter',innerDiameter);            % magnet inner diameter
motorProps=setParamafm(motorProps,'Magnet outer diameter',outerDiameter);            % magnet outer diameter
motorProps=setParamafm(motorProps,'Magnet spacing',magnetSpacing);                   % magnet spacing
motorProps=setParamafm(motorProps,'Rotor yoke height',yokeHeight);                   % yoke height
motorProps=setParamafm(motorProps,'Rotor height',magnetHeight+yokeHeight);           % rotor height
motorProps=setParamafm(motorProps,'Air gap',airGap);                                 % air gap
end  

%________________________________________________________________________________________________________________________
 

function innerDiameter = checkInnerDiameter(innerDiameter,slotWidth,nSlots)
% check if slot width is consitent with inner diameter
N = nSlots;                  % number of polygon angles
S = slotWidth;               % polygon side
alfa = pi/N;
Rmin1 = S/2/sin(alfa);       % minimum radius
Rmin = 1.03*Rmin1;            % 3% margin
if innerDiameter<2*Rmin
    innerDiameter = 2*Rmin;  % increase inner diameter
end

end
