% This file generates a reference case 
% 10-by-10 on fixed columns, orthogonal grid
close all;
clear;
clc;
commandwindow; 
%% Generate the triangular mesh: rectangular geometry
Xmin = 0; Xmax = 10;
Ymin = 0; Ymax = 10;
geometry.xmin = Xmin; geometry.xmax = Xmax;
geometry.ymin = Ymin; geometry.ymax = Ymax;
% Define two meshgrids
resol = 16; 
dx = Xmax/resol; dy = Ymax/resol;
[Xgrid1,Ygrid1] = meshgrid(0:dx:Xmax,0:dy:Ymax);
[Xgrid2,Ygrid2] = meshgrid(dx/2:dx:Xmax-dx/2,dy/2:dy:Ymax-dy/2);
% Node list
nodes1 = [Xgrid1(:),Ygrid1(:)];
nodes2 = [Xgrid2(:),Ygrid2(:)];
nnodes1 = size(nodes1,1);
nnodes2 = size(nodes2,1);
nodes = [nodes1;nodes2];
nnodes = nnodes1 + nnodes2;
% Loop over blocks
ntri = 4*resol^2;
plates = zeros(ntri,3);
k = 1;
for i = 1:resol
    for j = 1:resol
        n1 = (i-1)*(resol+1)+j+1;
        n2 = n1 - 1;
        n3 = nnodes1 + (i-1)*(resol)+j;
        plates(k,1:3) = [n1,n2,n3];
        n1 = (i-1)*(resol+1)+j;
        n2 = n1 + resol + 1;
        n3 = nnodes1 + (i-1)*(resol)+j;
        plates(k+1,1:3) = [n1,n2,n3];
        n1 = (i)*(resol+1)+j;
        n2 = n1 + 1;
        n3 = nnodes1 + (i-1)*(resol)+j;
        plates(k+2,1:3) = [n1,n2,n3];
        n1 = (i)*(resol+1)+j+1;
        n2 = n1 - (resol+1);
        n3 = nnodes1 + (i-1)*(resol)+j;
        plates(k+3,1:3) = [n1,n2,n3];
        k = k + 4;
    end
end

% Create a pde mesh and view
model = createpde();
geometryFromMesh(model,nodes',plates');
% figure(1);
% pdeplot(model,NodeLabels='on');
% axis equal tight
% keyboard;

% Find connectivities
[conn,connnum,count] = neighborelem(plates,max(max(plates)));

% Compute areas of triangles
[~,plates] = areaplate(nodes,plates);

%% Material, load, setup
E = 30e6; nu = 0.2; massden = 25; 
loadmag = 5; 
defallow = Xmax/750;
type = '4corners_fixed'; %'simply_supp';% 'simply_supp_SYM'; %'4corners_fixed'; %'4corners_fixed_SYM';% 'SSSYM'; ; 'SSSYM'
%% Generate data for beams
% Create list of edges
nplates = size(plates,1);
nedges = nplates*3; % Number of non-unique edges
edges = zeros(nedges,2);
for i = 1:nplates
    edges(3*i-2,1:2) = [plates(i,1),plates(i,2)];
    edges(3*i-1,1:2) = [plates(i,2),plates(i,3)];
    edges(3*i,1:2) = [plates(i,3),plates(i,1)];
end
% Organize such that n1<n2
for i = 1:nedges
    if edges(i,1) > edges(i,2)
        tmp = edges(i,1);
        edges(i,1) = edges(i,2);
        edges(i,2) = tmp;
    end
end
edges_srt = sortrows(edges);
edges = unique(edges_srt,'rows');
nedges = size(edges,1);
nbeams = nedges;
%% Apply BC
ndof = 6*nnodes;
supdofs = applyBC(nodes,type,geometry);
freedofs = setdiff(1:ndof,supdofs);
%% Initialization
ndv1 = nbeams + 2; % b, h, th
ndv2 = nnodes*2; % X and Y
ndv = ndv1 + ndv2;
% Initial design
b_peri = 0.1; % Beams on perimeter
b_intrn = 0.1; % Internal beams
b_list = zeros(nbeams,1);
lastnode = nnodes1;
for i = 1:nedges
    n1 = edges(i,1);
    n2 = edges(i,2);
    % Skip edges with nodes outside the ortho grid 
    if max(n1,n2)>lastnode; continue; end % Skip
    if (n2-n1)==1 % y-dir
        if n1<(resol+1) % left edge
            b_list(i,1) = b_peri;
        elseif n1>(resol*(resol+1)) % right edge
            b_list(i,1) = b_peri;
        else 
            b_list(i,1) = b_intrn;
        end
    end
    if (n2-n1)==(resol+1) % x-dir
        if mod(n1,resol+1) == 1 % bottom edge
            b_list(i,1) = b_peri;
        elseif mod(n1,resol+1) == 0 % bottom edge
            b_list(i,1) = b_peri;
        else 
            b_list(i,1) = b_intrn;
        end
    end
end
h_total = 0.53; % Total thickness
th = 0.08; % Plate thickness 
h_beam = h_total-th; % Beam height 
%% Computations
% Volume of beams
V = 0;
dv1 = zeros(ndv1,1);
dv2 = zeros(2*nnodes,1);
[V,~,~] = volbeam(V,dv1,dv2,nodes,edges,b_list,h_beam);
% Volume of plates
[totalArea,plates] = areaplate(nodes,plates);
V = V + totalArea*th;
% Assemble external load: distributed in Z
Fext = zeros(ndof,1);
for i = 1:nplates
    n1 = plates(i,1); n2 = plates(i,2); n3 = plates(i,3);
    loaddof = [6*n1-3 6*n2-3 6*n3-3];
    Fext(loaddof,1) = Fext(loaddof,1) - loadmag*plates(i,4)/3;
end
% Assemble self weight load
Fsw = 0*Fext;
[Fsw,~] = assmFplate(Fsw,nodes,plates,th,massden);
[Fsw,~] = assmFbeam(Fsw,nodes,edges,b_list,h_beam,massden);
% Assemble ribbed plate stiffness
K = zeros(ndof);
K = assmKplate(K,nodes,plates,E,nu,th);
K = assmKbeam(K,nodes,edges,E,nu,b_list,h_beam,th);
% Solve
F = Fext + Fsw;
U = 0*F;
U(freedofs,1) = K(freedofs,freedofs) \ F(freedofs,1);
comp = F'*U;
% Compute maximum deflection
Uver = U(3:6:end,1);
maxUver = max(-Uver);
pN = 16;
apprxmaxUver = (sum(Uver.^pN))^(1/pN);
% Draw
msh = model.Mesh;
pdeplot(msh,XYData=Uver,Mesh="off");
axis equal tight off
view(2)
