function fh = pdesample
fh = localfunctions;
end

function out = pdesample_input(in)
% pde model calculation
%
% Copyright 2016 The MathWorks, Inc.

material = in(1);
x1 = in(1+1);
y1 = in(2+1);
rx = in(3+1);
ry = in(4+1);
dh = in(5+1);
theta = in(6+1);

try
% Create a rectangle and two adjoining circles
% Rectangle
rect1 = [3; 4; 0; 10; 10; 0; 0; 0; 5; 5];
% Geometry description:
rect2 = [3; 4; 10 - 2; 10; 10; 10 - 2; 5; 5; 5 + dh; 5+dh];

% Ellipse
C1 = [4; x1; y1; rx; ry; theta];

% Append extra zeros to the circles
C1 = [C1;zeros(length(rect1) - length(C1),1)];
% Combine the shapes into one matrix
gd = [rect1,rect2,C1];

%% 
% Give names for the three shapes
ns = char('rect1','rect2','C2');
ns = ns';

% Specify the union of the rectangle and C1, and subtract C2
sf = '(rect1+rect2)-C2';

[dl,bt] = decsg(gd,sf,ns); % combines the basic shapes using the set formula
[dl2,~] = csgdel(dl,bt); % removes subdomain boundaries
model = createpde(2);
geometryFromEdges(model,dl2);
% pdegplot(model,'EdgeLabels','on');
% axis equal

%%
clampedBC = applyBoundaryCondition(model,'Edge',6,'u',[0 0]);
sigma = 1e3;
presBC = applyBoundaryCondition(model,'Edge',1,'g',[0 -sigma]);
% setInitialConditions(model,0,0);

if material == 1
E = 200e9; % Young's modulus of the material, lbs/in^2
gnu = .3; % Poisson's ratio of the material
rho = .3/386; % Density of the material
elseif material == 2 % alm
E = 73e9; % Young's modulus of the material, lbs/in^2
gnu = .34; % Poisson's ratio of the material
rho = 2.70e-04; % Density of the material
end

G = E/(2.*(1+gnu));
mu = 2*G*gnu/(1-gnu);
c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
f = [0 0]'; % No body forces
specifyCoefficients(model, 'm', 0, 'd', 0, 'c', c, 'a', 0, 'f', f);

generateMesh(model,'Hmax',0.2);
results = solvepde(model);
u = results.NodalSolution;
color = sqrt(u(:,1).^2 + u(:,2).^2);
[p,e,t] = meshToPet(model.Mesh);
p2 = p + u'*1000000; % scale up
pdeplot(p2,e,t,'xydata',color,'mesh','off','contour','on');
colorbar('off')
xlim([-1 11])
axis equal
idx = ismember(p',[10 0],'rows');
out = abs(color(idx)'*1000000); % scale up
if out == 0
    out = 100;
end

catch
    out = 100;
end
end
%% 
% pdesample_input_doe_sn (S/N calculation)

function [sn resp] = pdesample_input_doe_sn(in)
% SN calculation
%
% Copyright 2016 The MathWorks, Inc.

%% in

L = orthogonalarray('L8');
idx = L == 1;
alp = 0.01;
L8 = ones(8,7) - alp*idx + alp*(~idx);
data = [ones(8,1) L8].*in;
resp = ones(8,1);
%%
for k = 1:8
    in1 = data(k,:);
    
    material = in1(1);
    x1 = in1(1+1);
    y1 = in1(2+1);
    rx = in1(3+1);
    ry = in1(4+1);
    dh = in1(5+1);
    theta = in1(6+1);
    
    
    %%
    
    % Create a rectangle and two adjoining circles
    % Rectangle
    rect1 = [3; 4; 0; 10; 10; 0; 0; 0; 5; 5];
    % Geometry description:
    rect2 = [3; 4; 10 - 2; 10; 10; 10 - 2; 5; 5; 5 + dh; 5+dh];
    
    % Ellipse
    % pdeellip(7,2,1,1,0,'E1');
    C1 = [4; x1; y1; rx; ry; theta];
    
    % Append extra zeros to the circles
    C1 = [C1;zeros(length(rect1) - length(C1),1)];
    % Combine the shapes into one matrix
    gd = [rect1,rect2,C1];
    
    %%
    % Give names for the three shapes
    ns = char('rect1','rect2','C2');
    ns = ns';
    
    % Specify the union of the rectangle and C1, and subtract C2
    sf = '(rect1+rect2)-C2';
    
    [dl,bt] = decsg(gd,sf,ns); % combines the basic shapes using the set formula
    [dl2,~] = csgdel(dl,bt); % removes subdomain boundaries
    model = createpde(2);
    geometryFromEdges(model,dl2);
    pdegplot(model,'EdgeLabels','on');
    axis equal
    clampedBC = applyBoundaryCondition(model,'Edge',6,'u',[0 0]);
    sigma = 1e3;
    presBC = applyBoundaryCondition(model,'Edge',1,'g',[0 -sigma]);
    % setInitialConditions(model,0,0);
    
    if material == 1
        E = 200e9; % Young's modulus of the material, lbs/in^2
        gnu = .3; % Poisson's ratio of the material
        rho = .3/386; % Density of the material
        
    elseif material == 2
        E = 73e9; % Young's modulus of the material, lbs/in^2
        gnu = .34; % Poisson's ratio of the material
        rho = 2.70e-04; % Density of the material
    end
    
    G = E/(2.*(1+gnu));
    mu = 2*G*gnu/(1-gnu);
    c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
    f = [0 0]'; % No body forces
    specifyCoefficients(model, 'm', 0, 'd', 0, 'c', c, 'a', 0, 'f', f);
    
    generateMesh(model,'Hmax',0.2);
    results = solvepde(model);
    u = results.NodalSolution;
    color = sqrt(u(:,1).^2 + u(:,2).^2);
    [p,e,t] = meshToPet(model.Mesh);
    p2 = p + u'*1000000;
    pdeplot(p2,e,t,'xydata',color,'mesh','off','contour','on');
    colorbar('off')
    xlim([-1 11])
    axis equal
    idx = ismember(p',[10 0],'rows');
    out = abs(color(idx)'*1000000);
    resp(k,1) = out;
end
% SN
sn = 10*log((sum(mean(resp).^2))/sum((resp - mean(resp)).^2)/7);
end


%% COST Function(Single Optimization)

function out = myobjfun(in)

out = pdesample_input([1 in(1) in(2) 0.5 0.2 in(3) in(4)]);

end
%% COST Function (Multi Optimizaion)

function out = myobjfun2(in)

out(1) = pdesample_input([1 in 0.2 0*pi/12]);
out(2) = -pi*in(3)*in(4);

end

%% nonlinear constraint

function [c,ceq] = mynoncon(x)
%
% Copyright 2016 The MathWorks, Inc.

c(1) = x(1) + abs(0.5*cos(x(4))) - 9.8 ;
c(3) =  -x(2) +  abs(0.5*sin(x(4))) + 0.2;
if 0.5*cos(x(4)) > 0.2 
c(2) = x(2) +  abs(0.5*sin(x(4))) - 4.8;
else
c(2) = x(2) +  0.2 - 4.9;
end
    
ceq = [];

end

