clear;clc;
if 1

obj= readObj('./SHOW_DS/cube_4__sf.obj');%%

S.Vertices=obj.v;%%
S.Faces=obj.f.v;%%

cor=load("./SHOW_DS/cube_cor_4.txt");
elem=load("./SHOW_DS/cube_elem_4.txt");

DBC=@(X) abs(X(:,2)-0)<1e-4;
NBC=@(X) abs(X(:,2)-4)<1e-2;


fea=FEA(cor,elem);
fea.set_DBC(DBC,[1 1 1]);
order=5;
[GlobalGP,GlobalGW]=fea.set_NBC(NBC,order,[0 0 -1e-6]);
fprintf("##Total gs_w: %.2f \n",sum(GlobalGW(:)));


save NBC_GlobalGPW GlobalGP GlobalGW

fea.solve();

[Q_U,Q_S]=fea.query_US(S.Vertices);
J_FEA=fea.J;
save US.mat Q_U Q_S J_FEA
[Q_tU,Q_tS]=fea.query_corner_US(fea.cor);
Vert_tot=fea.cor;
save USt.mat Q_tU Q_tS Vert_tot
end

%
% extract FEA data & Surface data S

if ~exist("Q_U",'var')
    load US.mat Q_U Q_S J_FEA
    load USt.mat Q_tU Q_tS Vert_tot
    if ~exist("obj",'var')
        obj= readObj('./SHOW_DS/cube_4__sf.obj');%%
        S.Vertices=obj.v;%%
        S.Faces=obj.f.v;%%
    end
end

norm_U_FEA=vecnorm(Q_U,2,2);
von_S_FEA=sqrt(0.5*((Q_S(:,1)-Q_S(:,2)).^2+(Q_S(:,1)-Q_S(:,3)).^2+(Q_S(:,3)-Q_S(:,2)).^2)+...
     3*(vecnorm(Q_S(:,4:6),2,2).^2));

%---------tot U S
norm_U_total_FEA=vecnorm(Q_tU,2,2);
von_S_total_FEA=sqrt(0.5*((Q_tS(:,1)-Q_tS(:,2)).^2+(Q_tS(:,1)-Q_tS(:,3)).^2+(Q_tS(:,3)-Q_tS(:,2)).^2)+ ...
    3*(vecnorm(Q_tS(:,4:6),2,2).^2));


% extract fcm data

if ~exist('plt')
    plt=Plotter('Cube_tmp','NNew_Data');
end
plt.mp.viewxy=[90 70];
cur_c=plt.mp.recorder.intermediate_c{1};
Q_U_XV=compute_u(S.Vertices,plt.mp.recorder.U{1},plt.mp);%3x1xLen
Q_U_XV=reshape(Q_U_XV,3,[])';% lenx3
Q_S_XV=squeeze(compute_stress(S.Vertices,plt.mp.recorder.U{1},plt.mp))';% len x  6
J_XV=plt.mp.recorder.J(1);

norm_U_FCM=vecnorm(Q_U_XV,2,2);
von_S_FCM=sqrt(0.5*((Q_S_XV(:,1)-Q_S_XV(:,2)).^2+(Q_S_XV(:,1)-Q_S_XV(:,3)).^2+(Q_S_XV(:,3)-Q_S_XV(:,2)).^2)+...
     3*(vecnorm(Q_S_XV(:,4:6),2,2).^2));

%------------tot U S of FCM
Q_tU_XV=compute_u(Vert_tot,plt.mp.recorder.U{1},plt.mp);
Q_tU_XV=reshape(Q_tU_XV,3,[])';% lenx3
norm_U_total_FCM=vecnorm(Q_tU_XV,2,2);
Q_tS_XV=squeeze(compute_stress(Vert_tot,plt.mp.recorder.U{1},plt.mp))';
von_S_total_FCM=sqrt(0.5*((Q_tS_XV(:,1)-Q_tS_XV(:,2)).^2+(Q_tS_XV(:,2)-Q_tS_XV(:,3)).^2+(Q_tS_XV(:,3)-Q_tS_XV(:,1)).^2)+...
    3*(vecnorm(Q_tS_XV(:,4:6),2,2).^2));


% plot surface

%--------------------------------------------------------------------plot deformation
% close all;
plot_surf(S,norm_U_FEA,'norm deformation of FEA');% plot norm_U_FEA
plot_surf(S,norm_U_FCM,'norm deformation of FCM');% plot norm_U_FCM
%--------------------------------------------------------------------plot stress
plot_surf(S,von_S_FEA,'Von Mises stress of FEA')% plot von_S_FEA
plot_surf(S,von_S_FCM,'Von Mises stress of FCM')% plot von_S_FEA

%--------------------------------------------------------hist of U
figure;
histogram(norm_U_FEA);
hold on;
histogram(norm_U_FCM);
xlabel('deformation distribution of FEA and FCM');
legend({'FEA','FCM'},'Location',"best");
%--------------------------------------------------------hist of S
figure;
histogram(von_S_FEA);
hold on;
histogram(von_S_FCM);
xlabel('von Mises stress distribution of FEA and FCM');
legend({'FEA','FCM'},'Location',"best");

%
%% Surface Statistic

%-----------------------------------------compliance
compliance_FEA=J_FEA
compliance_FCM=J_XV
abs(compliance_FEA-compliance_FCM)./compliance_FCM
%----------------------------
max_U_FEA=max(norm_U_FEA(:))
max_U_FCM=max(norm_U_FCM(:))
abs(max_U_FEA-max_U_FCM)/max_U_FCM
max_S_FEA=max(von_S_FEA(:))
max_S_FCM=max(von_S_FCM(:))
abs(max_S_FEA-max_S_FCM)/max_S_FCM
%------------------------------------------relative error
threshold=0.2;
%---------------------------------for U
sorted_norm_U_FEA=sort(norm_U_FEA);
thresh_val=sorted_norm_U_FEA(round(size(sorted_norm_U_FEA,1)*threshold));
mini_U_idx=(norm_U_FEA<thresh_val);
% mini_U_idx=(norm_U_FEA<0);%% wait to del
rela_error_U=abs(norm_U_FEA-norm_U_FCM)./norm_U_FEA;
rela_error_U(mini_U_idx)=0;
rela_error_U(rela_error_U==Inf)=0;
mean(rela_error_U)
max(rela_error_U)
%----------------------------------for S
sorted_von_S_FEA=sort(von_S_FEA);
thresh_val=sorted_von_S_FEA(round(size(sorted_von_S_FEA,1)*threshold));
mini_S_idx=von_S_FEA<thresh_val;
mini_S_idx=(mini_S_idx<0);%% wait to del
rela_error_S=abs(von_S_FEA-von_S_FCM)./von_S_FEA;
rela_error_S(mini_S_idx)=0;
rela_error_S(rela_error_S==Inf)=0;
mean(rela_error_S)
max(rela_error_S)


%-------------------------------------------------------plot rela error U
plot_surf(S,rela_error_U,'relative error of deformation between FEA and FCM');
plot_surf(S,rela_error_S,'relative error of Von Mises stress between FEA and FCM');
%-------------------------------------------------------hist relative error
figure;
histogram(rela_error_U);
xlabel('hist of relative error of deformation');

figure;
histogram(rela_error_S);
xlabel('hist of relative error of Von Mises stress');
% legend({'FEA','FCM'},'Location',"best");

%% total vertex analysis

%-----------------------------------------compliance
compliance_total_FEA=J_FEA
compliance_total_FCM=J_XV
abs(compliance_total_FEA-compliance_total_FCM)./compliance_total_FCM
%--------------------------------------max
max_U_total_FEA=max(norm_U_total_FEA(:))
max_U_total_FCM=max(norm_U_total_FCM(:))
abs(max_U_total_FEA-max_U_total_FCM)/max_U_total_FCM
max_S_total_FEA=max(von_S_total_FEA(:))
max_S_total_FCM=max(von_S_total_FCM(:))
abs(max_S_total_FEA-max_S_total_FCM)/max_S_total_FCM

%------------------------------------------relative error
threshold=0.2;
%---------------------------------for U
sorted_norm_U_FEA=sort(norm_U_total_FEA);
thresh_val=sorted_norm_U_FEA(round(size(sorted_norm_U_FEA,1)*threshold));
mini_U_idx=(norm_U_total_FEA<thresh_val);
% mini_U_idx=(norm_U_FEA<0);%% wait to del
rela_error_U=abs(norm_U_total_FEA-norm_U_total_FCM)./norm_U_total_FEA;
rela_error_U(mini_U_idx)=0;
rela_error_U(rela_error_U==Inf)=0;
mean(rela_error_U)
max(rela_error_U)
%----------------------------------for S
sorted_von_S_FEA=sort(von_S_total_FEA);
thresh_val=sorted_von_S_FEA(round(size(sorted_von_S_FEA,1)*threshold));
mini_S_idx=von_S_total_FEA<thresh_val;
% mini_S_idx=(mini_S_idx<0);%% wait to del
rela_error_S=abs(von_S_total_FEA-von_S_total_FCM)./von_S_total_FEA;
rela_error_S(mini_S_idx)=0;
rela_error_S(rela_error_S==Inf)=0;
mean(rela_error_S)
max(rela_error_S)

%
%% fcmlab contrast

load /home/cflin/MatlabCode/fcmlab/drivers/FCM_3D_AnalysisDriver/PlasticAnalysis/FCMlab_X_U_u_stress_C.mat X u C U stress
norm_U_fcmlab=vecnorm(u,2,2);
von_S_fcmlab=cal_von_Mises_Stress(stress);
U_fcmlab=U;
if ~exist('plt')
    plt=Plotter('Cube_tmp','NNew_Data');
end
plt.mp.viewxy=[90 70];
cur_c=plt.mp.recorder.intermediate_c{1};
Q_U_FCM=compute_u(X,plt.mp.recorder.U{1},plt.mp);%3x1xLen
Q_U_FCM=reshape(Q_U_FCM,3,[])';% lenx3
Q_S_FCM=squeeze(compute_stress(X,plt.mp.recorder.U{1},plt.mp))';% len x  6
J_XV=plt.mp.recorder.J(1);

norm_U_FCM=vecnorm(Q_U_FCM,2,2);
von_S_FCM=cal_von_Mises_Stress(Q_S_FCM);

% analysis
compliance_fcmlab=C
compliance_FCM=J_XV

abs(compliance_fcmlab-compliance_FCM)./compliance_fcmlab
relative_error=zeros(size(norm_U_FCM));
zero_U_idx=norm_U_fcmlab==0|norm_U_FCM==0;
relative_error(~zero_U_idx)=abs(norm_U_FCM(~zero_U_idx)-norm_U_fcmlab(~zero_U_idx))./norm_U_fcmlab(~zero_U_idx);
% relative_error=abs(norm_U_FCM-norm_U_fcmlab)./norm_U_fcmlab;
% relative_error(relative_error==Inf)=0;
mean(relative_error)

max(relative_error)

relative_error=abs(von_S_FCM-von_S_fcmlab)./von_S_fcmlab;
relative_error(relative_error==Inf)=0;
mean(relative_error)

max(relative_error)
% U_FCM=plt.mp.recorder.U{1};
% figure
% histogram(U_FCM);
% hold on;
% histogram(U_fcmlab);
% legend({'FCM','fcmlab'},'Location',"best");
% 
% sU_FCM=sort(U_FCM);
% sU_fcmlab=sort(U_fcmlab);


