%% Script for remeshing of example 4, run 3
% We remesh the distorted part and then repeat the analysis
clear
load("ex4_run3_remesh.mat");
nodes_deleted = [109 110 111 126 127 143 144 145 160 176];
%% Renumbering of nodes
nnodes_remesh = nnodes - length(nodes_deleted);
nodes_remesh = zeros(nnodes_remesh,3);
shift = 0;
for i = 1:nnodes
    if ismember(i,nodes_deleted)
        shift = shift + 1;
    else
        nodes_remesh(i-shift,1:2) = nodes(i,1:2);
        nodes_remesh(i-shift,3) = i; % keep old node number
    end
end
%% Identify edges to delete
edges_deleted(1:nedges) = 0;
k = 0;
for i = 1:nedges
    if ismember(edges(i,1),nodes_deleted)
        k = k + 1;
        edges_deleted(k) = i;
    elseif ismember(edges(i,2),nodes_deleted)
        k = k + 1;
        edges_deleted(k) = i;
    end
end
edges_deleted(k+1:end) = [];
nedges_deleted = length(edges_deleted);
%% Create new edges array
nedges_remesh = nedges - nedges_deleted;
edges_remesh = zeros(nedges_remesh,3);
k = 0;
for i = 1:nedges
    if  ~ismember(i,edges_deleted)
        k = k + 1;
        node1_old = edges(i,1);
        [node1_new,~] = find(nodes_remesh(:,3)==node1_old);
        node2_old = edges(i,2);
        [node2_new,~] = find(nodes_remesh(:,3)==node2_old);
        edges_remesh(k,1:2) = [node1_new,node2_new];
        edges_remesh(k,3) = i;
    end
end
%% Add new edges
% first with old numbering
edges_new = ...
    [175 191;
    175 177;
    159 177;
    142 177;
    142 161;
    125 161;
    125 162;
    108 162;
    108 146;
    108 128;
    91 128;
    92 128;
    92 129;
    93 129;
    93 112];
for i = 1:15
    k = k + 1;
    node1_old = edges_new(i,1);
    [node1_new,~] = find(nodes_remesh(:,3)==node1_old);
    node2_old = edges_new(i,2);
    [node2_new,~] = find(nodes_remesh(:,3)==node2_old);
    edges_remesh(k,1:2) = [node1_new,node2_new];
end
nedges_remesh_copied = nedges_remesh;
nedges_remesh = size(edges_remesh,1);
%% Copy beam data
b_list_remesh = zeros(nedges_remesh,1);
for i = 1:nedges_remesh_copied
    b_list_remesh(i,1) = b_list_p(edges_remesh(i,3),1);
end
b_list_remesh(i+1:end,1) = b_min;
%% Identify plates to delete
plates_deleted(1:nplates) = 0;
k = 0;
for i = 1:nplates
    if ismember(plates(i,1),nodes_deleted)
        k = k + 1;
        plates_deleted(k) = i;
    elseif ismember(plates(i,2),nodes_deleted)
        k = k + 1;
        plates_deleted(k) = i;
    elseif ismember(plates(i,3),nodes_deleted)
        k = k + 1;
        plates_deleted(k) = i;
    end
end
plates_deleted(k+1:end) = [];
nplates_deleted = length(plates_deleted);
%% Create new plates array
nplates_remesh = nplates - nplates_deleted;
plates_remesh = zeros(nplates_remesh,3);
k = 0;
for i = 1:nplates
    if  ~ismember(i,plates_deleted)
        k = k + 1;
        node1_old = plates(i,1);
        [node1_new,~] = find(nodes_remesh(:,3)==node1_old);
        node2_old = plates(i,2);
        [node2_new,~] = find(nodes_remesh(:,3)==node2_old);
        node3_old = plates(i,3);
        [node3_new,~] = find(nodes_remesh(:,3)==node3_old);
        plates_remesh(k,1:3) = [node1_new,node2_new,node3_new];
    end
end
%% Add new plates
% first with old numbering
plates_new = ...
    [175 191 190;
    175 177 191;
    159 177 175;
    142 177 159;
    142 161 177;
    125 161 142;
    125 162 161;
    108 162 125;
    108 146 162;
    108 128 146;
    91 128 108;
    92 128 91;
    92 129 128;
    93 129 92;
    93 112 129;
    93 94 112];
for i = 1:16
    k = k + 1;
    node1_old = plates_new(i,1);
    [node1_new,~] = find(nodes_remesh(:,3)==node1_old);
    node2_old = plates_new(i,2);
    [node2_new,~] = find(nodes_remesh(:,3)==node2_old);
    node3_old = plates_new(i,3);
    [node3_new,~] = find(nodes_remesh(:,3)==node3_old);
    plates_remesh(k,1:3) = [node1_new,node2_new,node3_new];
end
nplates_remesh_copied = nplates_remesh;
nplates_remesh = size(plates_remesh,1);

%% Re-analyze
% Apply BC
ndof = 6*nnodes_remesh;
% Renumbering of nodes for BC
innerBound_remesh = 0*innerBound;
for i = 1:size(innerBound,1)
    node_old = innerBound(i,1);
    [node_new,~] = find(nodes_remesh(:,3)==node_old);
    innerBound_remesh(i,1) = node_new;
end
outerBound_remesh = 0*outerBound;
for i = 1:size(outerBound,1)
    node_old = outerBound(i,1);
    [node_new,~] = find(nodes_remesh(:,3)==node_old);
    outerBound_remesh(i,1) = node_new;
end
edge1Bound_remesh = 0*edge1Bound;
for i = 1:size(edge1Bound,1)
    node_old = edge1Bound(i,1);
    [node_new,~] = find(nodes_remesh(:,3)==node_old);
    edge1Bound_remesh(i,1) = node_new;
end
edge2Bound_remesh = 0*edge2Bound;
for i = 1:size(edge2Bound,1)
    node_old = edge2Bound(i,1);
    [node_new,~] = find(nodes_remesh(:,3)==node_old);
    edge2Bound_remesh(i,1) = node_new;
end
supdofs = applyBCtri(innerBound_remesh,outerBound_remesh,...
    edge1Bound_remesh,edge2Bound_remesh);
freedofs = setdiff(1:ndof,supdofs);
% % Volume of beams
% [V_remesh,~,~] = volbeam(V,dv1,dv2,nodes,edges,b_list,h_beam);
%     dv2(2*activeNodeX-1,1) = LapTX(activeNodeX,activeNodeX)*dv2(2*activeNodeX-1,1);
%     dv2(2*activeNodeY,1) = LapTY(activeNodeY,activeNodeY)*dv2(2*activeNodeY,1);
%     dv2(activeDOFs==0) = [];
% Volume of plates
[totalArea,plates] = areaplate(nodes_remesh,plates_remesh);
% V = V + totalArea*th;
% dv1(ndv1,1) = totalArea;
% Assemble external load: distributed in Z
Fext = zeros(ndof,1);
for i = 1:nplates_remesh
    n1 = plates_remesh(i,1); n2 = plates_remesh(i,2); n3 = plates_remesh(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_remesh,plates_remesh,th,massden);
[Fsw,~] = assmFbeam(Fsw,nodes_remesh,edges_remesh,b_list_remesh,h_beam,massden);
% Assemble plate-beam stiffness
K = zeros(ndof);
K = assmKplate(K,nodes_remesh,plates_remesh,E,nu,th);
K = assmKbeam(K,nodes_remesh,edges_remesh,E,nu,b_list_remesh,h_beam,th);
% Solve
F = Fext + Fsw;
U_remesh = 0*F;
U_remesh(freedofs,1) = K(freedofs,freedofs) \ F(freedofs,1);
comp = F'*U_remesh;
% Compute maximum deflection
Uver = U_remesh(3:6:end,1);
maxUver = max(-Uver);
pN = 16;
apprxmaxUver = (sum(Uver.^pN))^(1/pN);
% Create a pde mesh and view
model = createpde();
geometryFromMesh(model,nodes_remesh(:,1:2)',plates(:,1:3)');
figure(1);
pdeplot(model,NodeLabels='on');
axis equal tight
figure(2)
msh = model.Mesh;
pdeplot(msh,XYData=Uver,ZData=Uver,Mesh="off");
axis equal tight off
view(2)