%19.01.2011
%input and output of the numerical models
close all
clear all
clc
format bank
%input files
fe_model_input_file_name = 'sample_mesh.key';
fe_model_output_file_name = 'sample_mesh_trias_.key';
fe_model_diagnostic_file_name = 'sample_mesh_diagnostic_.key';
%
fe_file_in = fopen(fe_model_input_file_name);
fe_file_out = fopen(fe_model_output_file_name,'w');
diag_file_out = fopen(fe_model_diagnostic_file_name,'w');
%input from file
id_element=1; id_node=1;
el_count = 1; new_el_ID = 1000000;
method_=2; %1 by nodes; 2 by min diagonal
min_AR = 0.80; %minimum Aspect Ratio
while true
text_line_in=fgetl(fe_file_in);
%look for EOF
if ~ischar(text_line_in)
break
end %if
%check keywords
keyword_line=(text_line_in(1:1));
if keyword_line == '*'
is_node = strcmp(text_line_in, '*NODE');
is_element = strcmp(text_line_in, '*ELEMENT_SHELL');
is_section = strcmp(text_line_in, '*SECTION_SHELL');
is_curve = strcmp(text_line_in, '*DEFINE_CURVE_TITLE');
is_material = strcmp(text_line_in, '*MAT_MODIFIED_PIECEWISE_LINEAR_PLASTICITY_TITLE');
is_end = strcmp(text_line_in, '*END');
end
%input nodes
if is_node == 1
if length(text_line_in) > 15
node_(id_node,1)=id_node; %ID
node_(id_node,2)=str2num(text_line_in(1:8)); %NODE ID
node_(id_node,3)=str2num(text_line_in(9:25)); %coord X
node_(id_node,4)=str2num(text_line_in(26:41)); %coord Y
node_(id_node,5)=str2num(text_line_in(42:57)); %coord Z
id_node=id_node+1;
end
end
%input elements
if is_element == 1
if length(text_line_in) > 15
%sample data
%12345678901234567890123456789012345678901234567890
%-------1-------2-------3-------4-------5-------6
% 122 9000 109451 4333 108459 108459
element_(id_element,1)=str2num(text_line_in(1:8)); %ELEMENT ID
element_(id_element,2)=str2num(text_line_in(9:16)); %PID
element_(id_element,3)=str2num(text_line_in(17:24)); %node 1
element_(id_element,4)=str2num(text_line_in(25:32)); %node 2
element_(id_element,5)=str2num(text_line_in(33:40)); %node 3
element_(id_element,6)=str2num(text_line_in(41:48)); %node 4
element_(id_element,7)=0;
%check for trias node 3 == node 4
is_tria = false;
if element_(id_element,5)==element_(id_element,6)
is_tria = true;
end
if is_tria == true
new_el_ID=100000+el_count;
new_element_(el_count,1)=new_el_ID;
new_element_(el_count,2)=element_(id_element,2);
new_element_(el_count,3)=element_(id_element,3);
new_element_(el_count,4)=element_(id_element,4);
new_element_(el_count,5)=element_(id_element,5);
new_element_(el_count,6)=element_(id_element,6);
corespondence_(id_element,1)=element_(id_element,1);
corespondence_(id_element,2)=new_el_ID;
n1_=[node_(n_id(1),3) node_(n_id(1),4) node_(n_id(1),5)];
n2_=[node_(n_id(2),3) node_(n_id(2),4) node_(n_id(2),5)];
n3_=[node_(n_id(3),3) node_(n_id(3),4) node_(n_id(3),5)];
t1_=[1 1 1; n1_(1) n2_(1) n3_(1); n1_(2) n2_(2) n3_(2)];
t2_=[1 1 1; n1_(2) n2_(2) n3_(2); n1_(3) n2_(3) n3_(3)];
t3_=[1 1 1; n1_(3) n2_(3) n3_(3); n1_(1) n2_(1) n3_(1)];
area_(1)=1/2*sqrt((det(t1_))^2+(det(t2_))^2+(det(t3_))^2);
l_(1)=sqrt((n1_- n2_)*(n1_- n2_)');
l_(2)=sqrt((n2_- n3_)*(n2_- n3_)');
l_(3)=sqrt((n3_- n1_)*(n3_- n1_)');
AR_(1)=4*sqrt(3)*area_(1)/(l_(1)^2+l_(2)^2+l_(3)^2);
if AR_(1) < min_AR
element_(id_element, 7)= element_(id_element,2)+1;
new_element_(el_count,2)=element_(id_element,7);
end
sh_(1)=2*area_(1)/l_(1)/l_(1);
sh_(2)=2*area_(1)/l_(2)/l_(2);
sh_(3)=2*area_(1)/l_(3)/l_(3);
corespondence_(id_element,4)=area_(1);
corespondence_(id_element,6)=AR_(1);
corespondence_(id_element,8)=1/min(sh_);
%
el_count=el_count+1;
end
if is_tria == false
%compute diagonals
n_id(1)=find(node_(:,2)==element_(id_element,3));
n_id(2)=find(node_(:,2)==element_(id_element,4));
n_id(3)=find(node_(:,2)==element_(id_element,5));
n_id(4)=find(node_(:,2)==element_(id_element,6));
n1_=[node_(n_id(1),3) node_(n_id(1),4) node_(n_id(1),5)];
n2_=[node_(n_id(2),3) node_(n_id(2),4) node_(n_id(2),5)];
n3_=[node_(n_id(3),3) node_(n_id(3),4) node_(n_id(3),5)];
n4_=[node_(n_id(4),3) node_(n_id(4),4) node_(n_id(4),5)];
d_n1n3 = sqrt((n1_-n3_)*(n1_-n3_)');
d_n2n4 = sqrt((n2_-n4_)*(n2_-n4_)');
[min_d, min_d_id] = min([d_n1n3 d_n2n4]);
if method_ == 1 min_d_id = 1; end;
if method_ == 2 min_d_id = min_d_id; end;
if min_d_id == 1
new_el_ID=100000+el_count;
new_element_(el_count,1)=new_el_ID;
new_element_(el_count,2)=element_(id_element,2);
new_element_(el_count,3)=element_(id_element,3);
new_element_(el_count,4)=element_(id_element,4);
new_element_(el_count,5)=element_(id_element,5);
new_element_(el_count,6)=element_(id_element,5);
corespondence_(id_element,1)=element_(id_element,1);
corespondence_(id_element,1)=element_(id_element,1);
corespondence_(id_element,2)=new_el_ID;
t1_=[1 1 1; n1_(1) n2_(1) n3_(1); n1_(2) n2_(2) n3_(2)];
t2_=[1 1 1; n1_(2) n2_(2) n3_(2); n1_(3) n2_(3) n3_(3)];
t3_=[1 1 1; n1_(3) n2_(3) n3_(3); n1_(1) n2_(1) n3_(1)];
area_(1)=1/2*sqrt((det(t1_))^2+(det(t2_))^2+(det(t3_))^2);
l_(1)=sqrt((n1_- n2_)*(n1_- n2_)');
l_(2)=sqrt((n2_- n3_)*(n2_- n3_)');
l_(3)=sqrt((n3_- n1_)*(n3_- n1_)');
AR_(1)=4*sqrt(3)*area_(1)/(l_(1)^2+l_(2)^2+l_(3)^2);
sh_(1)=2*area_(1)/l_(1)/l_(1);
sh_(2)=2*area_(1)/l_(2)/l_(2);
sh_(3)=2*area_(1)/l_(3)/l_(3);
corespondence_(id_element,4)=area_(1);
corespondence_(id_element,6)=AR_(1);
if AR_(1) < min_AR
element_(id_element, 7)= element_(id_element,2)+1;
new_element_(el_count,2)=element_(id_element,7);
end
corespondence_(id_element,8)=1/min(sh_);
el_count=el_count+1;
%
new_el_ID=100000+el_count;
new_element_(el_count,1)=new_el_ID;
new_element_(el_count,2)=element_(id_element,2);
new_element_(el_count,3)=element_(id_element,5);
new_element_(el_count,4)=element_(id_element,6);
new_element_(el_count,5)=element_(id_element,3);
new_element_(el_count,6)=element_(id_element,3);
corespondence_(id_element,3)=new_el_ID;
%compute areas
t1_=[1 1 1; n1_(1) n3_(1) n4_(1); n1_(2) n3_(2) n4_(2)];
t2_=[1 1 1; n1_(2) n3_(2) n4_(2); n1_(3) n3_(3) n4_(3)];
t3_=[1 1 1; n1_(3) n3_(3) n4_(3); n1_(1) n3_(1) n4_(1)];
area_(2)=1/2*sqrt((det(t1_))^2+(det(t2_))^2+(det(t3_))^2);
l_(1)=sqrt((n1_- n3_)*(n1_- n3_)');
l_(2)=sqrt((n3_- n4_)*(n3_- n4_)');
l_(3)=sqrt((n4_- n1_)*(n4_- n1_)');
AR_(2)=4*sqrt(3)*area_(2)/(l_(1)^2+l_(2)^2+l_(3)^2);
sh_(1)=2*area_(2)/l_(1)/l_(1);
sh_(2)=2*area_(2)/l_(2)/l_(2);
sh_(3)=2*area_(2)/l_(3)/l_(3);
corespondence_(id_element,5)=area_(2);
corespondence_(id_element,7)=AR_(2);
if AR_(2) < min_AR
element_(id_element, 7)= element_(id_element,2)+1;
new_element_(el_count,2)=element_(id_element,7);
end
corespondence_(id_element,9)=1/min(sh_);
el_count=el_count+1;
end
if min_d_id == 2
new_el_ID=100000+el_count;
new_element_(el_count,1)=new_el_ID;
new_element_(el_count,2)=element_(id_element,2);
new_element_(el_count,3)=element_(id_element,3);
new_element_(el_count,4)=element_(id_element,4);
new_element_(el_count,5)=element_(id_element,6);
new_element_(el_count,6)=element_(id_element,6);
corespondence_(id_element,1)=element_(id_element,1);
corespondence_(id_element,1)=element_(id_element,1);
corespondence_(id_element,2)=new_el_ID;
%compute areas
t1_=[1 1 1; n1_(1) n2_(1) n4_(1); n1_(2) n2_(2) n4_(2)];
t2_=[1 1 1; n1_(2) n2_(2) n4_(2); n1_(3) n2_(3) n4_(3)];
t3_=[1 1 1; n1_(3) n2_(3) n4_(3); n1_(1) n2_(1) n4_(1)];
area_(1)=1/2*sqrt((det(t1_))^2+(det(t2_))^2+(det(t3_))^2);
l_(1)=sqrt((n1_- n2_)*(n1_- n2_)');
l_(2)=sqrt((n2_- n4_)*(n2_- n4_)');
l_(3)=sqrt((n4_- n1_)*(n4_- n1_)');
AR_(1)=4*sqrt(3)*area_(1)/(l_(1)^2+l_(2)^2+l_(3)^2);
sh_(1)=2*area_(1)/l_(1)/l_(1);
sh_(2)=2*area_(1)/l_(2)/l_(2);
sh_(3)=2*area_(1)/l_(3)/l_(3);
corespondence_(id_element,4)=area_(1);
corespondence_(id_element,6)=AR_(1);
corespondence_(id_element,8)=1/min(sh_);
if AR_(1) < min_AR
element_(id_element, 7)= element_(id_element,2)+1;
new_element_(el_count,2)=element_(id_element,7);
end
corespondence_(id_element,9)=1/min(sh_);
el_count=el_count+1;
%
new_el_ID=100000+el_count;
new_element_(el_count,1)=new_el_ID;
new_element_(el_count,2)=element_(id_element,2);
new_element_(el_count,3)=element_(id_element,4);
new_element_(el_count,4)=element_(id_element,5);
new_element_(el_count,5)=element_(id_element,6);
new_element_(el_count,6)=element_(id_element,6);
corespondence_(id_element,3)=new_el_ID;
t1_=[1 1 1; n2_(1) n3_(1) n4_(1); n2_(2) n3_(2) n4_(2)];
t2_=[1 1 1; n2_(2) n3_(2) n4_(2); n2_(3) n3_(3) n4_(3)];
t3_=[1 1 1; n2_(3) n3_(3) n4_(3); n2_(1) n3_(1) n4_(1)];
area_(2)=1/2*sqrt((det(t1_))^2+(det(t2_))^2+(det(t3_))^2);
l_(1)=sqrt((n2_- n3_)*(n2_- n3_)');
l_(2)=sqrt((n3_- n4_)*(n3_- n4_)');
l_(3)=sqrt((n4_- n2_)*(n4_- n2_)');
AR_(2)=4*sqrt(3)*area_(2)/(l_(1)^2+l_(2)^2+l_(3)^2);
sh_(1)=2*area_(2)/l_(1)/l_(1);
sh_(2)=2*area_(2)/l_(2)/l_(2);
sh_(3)=2*area_(2)/l_(3)/l_(3);
corespondence_(id_element,5)=area_(2);
corespondence_(id_element,7)=AR_(2);
corespondence_(id_element,9)=1/min(sh_);
if AR_(2) < min_AR
element_(id_element, 7)= element_(id_element,2)+1;
new_element_(el_count,2)=element_(id_element,7);
end
corespondence_(id_element,9)=1/min(sh_);
el_count=el_count+1;
end
end
id_element = id_element + 1;
end
end
end
save('corespondence_table.dat', 'corespondence_', '-ascii', '-tabs');
save('model_elements.dat' , 'element_' , '-ascii' , '-tabs');
fprintf(fe_file_out,'*ELEMENT_SHELL'); fprintf(fe_file_out,'\r'); fprintf(fe_file_out,'\n');
for ii_=1:el_count-1
fprintf(fe_file_out,'%8.0f',new_element_(ii_,1));
fprintf(fe_file_out,'%8.0f',new_element_(ii_,2));
fprintf(fe_file_out,'%8.0f',new_element_(ii_,3));
fprintf(fe_file_out,'%8.0f',new_element_(ii_,4));
fprintf(fe_file_out,'%8.0f',new_element_(ii_,5));
fprintf(fe_file_out,'%8.0f',new_element_(ii_,6));
fprintf(fe_file_out,'\r'); fprintf(fe_file_out,'\n');
end
fclose(fe_file_in);
fclose(fe_file_out);
%diagnostic
fprintf(diag_file_out,'*ELEMENT_SHELL'); fprintf(diag_file_out,'\r'); fprintf(diag_file_out,'\n');
for ii_=1:length(element_(:,1))
fprintf(diag_file_out,'%8.0f',element_(ii_,1));
if element_(ii_,7) ~= 0
fprintf(diag_file_out,'%8.0f',element_(ii_,7));
else
fprintf(diag_file_out,'%8.0f',element_(ii_,2));
end
fprintf(diag_file_out,'%8.0f',element_(ii_,3));
fprintf(diag_file_out,'%8.0f',element_(ii_,4));
fprintf(diag_file_out,'%8.0f',element_(ii_,5));
fprintf(diag_file_out,'%8.0f',element_(ii_,6));
fprintf(diag_file_out,'\r'); fprintf(diag_file_out,'\n');
end
fclose(diag_file_out)
%statistics
stat_=[0 0 0 0 0 0 0 0 0 0];
for ii_=1:length(corespondence_(:,1))
if corespondence_(ii_,6)<0.1 && corespondence_(ii_,6)>0.0 stat_(1) =stat_(1)+1; end
if corespondence_(ii_,6)<0.2 && corespondence_(ii_,6)>0.1 stat_(2) =stat_(2)+1; end
if corespondence_(ii_,6)<0.3 && corespondence_(ii_,6)>0.2 stat_(3) =stat_(3)+1; end
if corespondence_(ii_,6)<0.4 && corespondence_(ii_,6)>0.3 stat_(4) =stat_(4)+1; end
if corespondence_(ii_,6)<0.5 && corespondence_(ii_,6)>0.4 stat_(5) =stat_(5)+1; end
if corespondence_(ii_,6)<0.6 && corespondence_(ii_,6)>0.5 stat_(6) =stat_(6)+1; end
if corespondence_(ii_,6)<0.7 && corespondence_(ii_,6)>0.6 stat_(7) =stat_(7)+1; end
if corespondence_(ii_,6)<0.8 && corespondence_(ii_,6)>0.7 stat_(8) =stat_(8)+1; end
if corespondence_(ii_,6)<0.9 && corespondence_(ii_,6)>0.8 stat_(9) =stat_(9)+1; end
if corespondence_(ii_,6)<1.0 && corespondence_(ii_,6)>0.9 stat_(10)=stat_(10)+1; end
%
if corespondence_(ii_,7)<0.1 && corespondence_(ii_,7)>0.0 stat_(1) =stat_(1)+1; end
if corespondence_(ii_,7)<0.2 && corespondence_(ii_,7)>0.1 stat_(2) =stat_(2)+1; end
if corespondence_(ii_,7)<0.3 && corespondence_(ii_,7)>0.2 stat_(3) =stat_(3)+1; end
if corespondence_(ii_,7)<0.4 && corespondence_(ii_,7)>0.3 stat_(4) =stat_(4)+1; end
if corespondence_(ii_,7)<0.5 && corespondence_(ii_,7)>0.4 stat_(5) =stat_(5)+1; end
if corespondence_(ii_,7)<0.6 && corespondence_(ii_,7)>0.5 stat_(6) =stat_(6)+1; end
if corespondence_(ii_,7)<0.7 && corespondence_(ii_,7)>0.6 stat_(7) =stat_(7)+1; end
if corespondence_(ii_,7)<0.8 && corespondence_(ii_,7)>0.7 stat_(8) =stat_(8)+1; end
if corespondence_(ii_,7)<0.9 && corespondence_(ii_,7)>0.8 stat_(9) =stat_(9)+1; end
if corespondence_(ii_,7)<1.0 && corespondence_(ii_,7)>0.9 stat_(10)=stat_(10)+1; end
end
close all
clc
interval_=[0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9];
bar_=bar(interval_+.05, stat_);
ch = get(bar_,'Children'); fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData');
upper_lim = 0.05 * sum(stat_);
axis([0.1 1.0 0 upper_lim])
for ii_=1:length(stat_)
if stat_(ii_)~=0
txt_=text(0.05+interval_(ii_), .8*stat_(ii_)/sum(stat_)*upper_lim, num2str(stat_(ii_)));
set(txt_, 'FontSize', 12);
set(txt_, 'FontWeight', 'Bold'); set(txt_, 'Rotation', 90);
set(txt_, 'BackgroundColor',[1 1 1]); set(txt_, 'Margin',4);
end
if interval_(ii_) > 0.3 && interval_(ii_) < 0.6
fvcd(fvd(ii_,1)) = 3;
set(ch,'FaceVertexCData',fvcd);
end
if interval_(ii_) > 0.5 && interval_(ii_) < 0.7
fvcd(fvd(ii_,1)) = 1;
set(ch,'FaceVertexCData',fvcd);
end
if interval_(ii_) > 0.7
fvcd(fvd(ii_,1)) = 2;
set(ch,'FaceVertexCData',fvcd);
end
end
xlabel('Aspect Ratio', 'FontSize', 12, 'FontWeight', 'Bold');
ylabel('Number of elements', 'FontSize', 12, 'FontWeight', 'Bold');
title('Normalized aspect ratio. Model Statistics', 'FontSize', 14, 'FontWeight', 'Bold')
colormap(lines(3));
