Plotting & Checking fiber sections in Matlab

[divamva]

What follows is a set of scripts that I am using to generate OpenSees output for fiber sections that is directly readable by a Matlab script for visualizing and checking fiber-sections. What follows is the standard "Wsection.tcl" routine modified to produce the required text-output, and then the Matlab file used to plot the output.

If for example you generate a section with tag 1 using this Wsection.tcl script, (MAKE SURE MatlabOUT==1) then load Matlab, go to the opensees analysis directory and run one of the following two commands

>> plotSection('WsectionfiberS1.out')
>> plotSection('WsectionfiberS1.out',1,1)


Wsection.tcl script:
------------------------
# Wsection.tcl: tcl procedure for creating a wide flange steel fiber section
# written: Remo M. de Souza
# date: 06/99
# modified: 08/99 (according to the new general modelbuilder)
# modified: 7/Oct/05, to include MatlabOut functionality, so
# that you can plot the section fibers in Matlab!

# +-----------------------+ + +
# | | | tf |
# +----------+ +----------+ + |
# | | | |
# y | | | |
# | | | | |
# z--+ ->| |<- tw | dw | d
# | | | |
# | | | |
# | | | |
# +----------+ +----------+ + |
# | | | tf |
# +-----------------------+ + +
# bf



# input parameters
# secID - section ID number
# matID - material ID number
# d = nominal depth
# tw = web thickness
# bf = flange width
# tf = flange thickness
# nfdw = number of fibers along web depth
# nftw = number of fibers along web thickness
# nfbf = number of fibers along flange width
# nftf = number of fibers along flange thickness

proc Wsection {secID matID d tw bf tf nfdw nftw nfbf nftf {MatlabOut 0}} {
set dw [expr $d - 2 * $tf]
if {$dw<=0} {
puts "Error: Wsection beam depth is too small"
exit
}
set y1 [expr -$d/2]
set y2 [expr -$dw/2]
set y3 [expr $dw/2]
set y4 [expr $d/2]
set z1 [expr -$bf/2]
set z2 [expr -$tw/2]
set z3 [expr $tw/2]
set z4 [expr $bf/2]
#
section fiberSec $secID {
# nfIJ nfJK yI zI yJ zJ yK zK yL zL
patch quadr $matID $nfbf $nftf $y1 $z4 $y1 $z1 $y2 $z1 $y2 $z4
patch quadr $matID $nftw $nfdw $y2 $z3 $y2 $z2 $y3 $z2 $y3 $z3
patch quadr $matID $nfbf $nftf $y3 $z4 $y3 $z1 $y4 $z1 $y4 $z4
}

if { $MatlabOut == 1 } {
set matfile [open "WsectionfiberS$secID.out" "w"]
puts $matfile "Wsection d=$d, tw=$tw, bf=$bf, tf=$tf"
puts $matfile "patch quadr $matID $nfbf $nftf $y1 $z4 $y1 $z1 $y2 $z1 $y2 $z4"
puts $matfile "patch quadr $matID $nftw $nfdw $y2 $z3 $y2 $z2 $y3 $z2 $y3 $z3"
puts $matfile "patch quadr $matID $nfbf $nftf $y3 $z4 $y3 $z1 $y4 $z1 $y4 $z4"
close $matfile
}
}




Plotsection.m (MATLAB script)
-------------------------------------------

function plotSection(sectionfname,barRplotfactor,fillflag,matcolor,edgecolor)
% plotSection(sectionfname,barRplotfactor,fillflag,matcolor,edgecolor)
%_________________________________________________________________________
% INPUT
% sectionfname : filename containing section info
% barRplotfactor: determines a magnification factor for plotting fibers
% produced via the "layer straight" command. Helps visualize
% the rebars (default = 1)
% fillflag : If 1, then the fibers are plotted as filled
% polygons/circles (default = 0)
% matcolor : a cell array containing color definitions for each
% material (default = {'y','b','r','g','m','k'} )
% edgecolor : same as above for the edges of the fibers. Used only when
% fillflag==1 (default = {'k'})
%
%__________________________________________________________________________
%
% Helps with taking FiberSection output (actually generation commands)
% and recreating the section fibers for visualization
%
% Examples
%__________________________________________________________________________
% plotSection('\users\divamva\opensees\divamva\Kavala\hollowRectfiberS1.out')
% plotSection('\users\divamva\opensees\divamva\Kavala\hollowRectfiberS1.out',2,1,{[0.7,0.7,0.7],[0.5,0.5,0.5]})
% plotSection('\users\divamva\opensees\divamva\Kavala\hollowRectfiberS1.out',2,1,{[0.7,0.7,0.7]},{[0.7,0.7,0.7]})
% plotSection('\users\divamva\opensees\divamva\RCsect\BGRCfiberS2.out')
% plotSection('\users\divamva\opensees\divamva\RCsect\RCsectionfiberS3.out')
%__________________________________________________________________________
% Created by D.Vamvatsikos
%
% Upgraded 07/May/05: Now the first line in the file is not discarded
% unless it contains an invalid command (e.g. a title)
% Also "support" is provided for "fiber", "layer circ" and
% "layer patch". A check is made to ensure "patch quad"
% coordinates create a valid convex and counterclockwise
% quadrangle!
% Upgraded 07/Oct/05: Added two extra checks for colinearity which were
% missing (must check all combinations of size 3)
%_____________________________________________________________________________
% NOTES:
% 1. Opensees will run even if a clockwise quad is supplied but then weird
% things will happen! Usually axial compression becomes tension and vice-versa
% 2. use view(90,90) to rotate the axes by 90deg and maintain nice labels!
% This is useful for plotting, e.g. the G9 or G7 sections.
% 3. The commands "fiber", "layer circ" and "layer patch" are not supported
% yet. They are only reckognized as being valid, but no plotting is
% done.


%cd e:\users\divamva\opensees\divamva\RCsect
if nargin<5, edgecolor={'k'}; end
if nargin<4, matcolor={'y','b','r','g','m','k'}; end
if nargin<3, fillflag=0; end
if nargin<2
% Do not increase bar radius plotting by default
barRplotfactor=1;
end

fid=fopen(sectionfname);

if fid>=0
% Since we do not use this "name" at all, refrain from reading it.
% If the user has not supplied such a name, then the valuable first
% line will be discarded!
%name = fgetl(fid);
i=0;
while 1
tline = fgetl(fid);
% check for EOF. Fully compatible with Matlab 7!
if ~ischar(tline), break, else i=i+1; end
% read fiber generation commands
[firstcom,count,err,ni]=sscanf(tline,'%s',1);
% only the "fiber" command has a single word as declaration.
% the rest use two words (e.g. "patch circ")
if strcmp(lower(firstcom),'fiber')
% one word was enough. Then numbers follow
gencom(i).com=firstcom;
else
% two words are needed.
[gencom(i).com,count,err,ni]=sscanf(tline,'%s',2);
end
gencom(i).nums=sscanf(tline(ni:end),'%f');
end
fclose(fid);
else
error(sprintf('File "%s" not found\n',sectionfname))
end

clf;
hold on

for i=1:length(gencom)
switch lower(gencom(i).com)
case 'patchquadr'
tmpcell=num2cell(gencom(i).nums);
[mcolor,nIJ,nJK,Iy,Iz,Jy,Jz,Ky,Kz,Ly,Lz]=deal(tmpcell{:});
% check whether they are (1)counterclockwise and (2)form a convex
% quadrilateral in the y,z coord system
%
% z L* *K
% |
% | I* *J
% |
% +-----y
%
% We will achieve this by calculating the outer products of
% the vectors IJ x IK and IK x IL. If both are positive,
% we have convex and counter-clockwise quadrangle.
% Since all vectors have no x-component, the outerproduct
% is always along the x-axis and it is calculated as
% a x b = | ya za | * (unit x-vector)
% | yb zb |
% where ya,za are the vector coords (ya=ya2-ya1, za=za2-za1)
outIJxIK=(Jy-Iy)*(Kz-Iz) - (Ky-Iy)*(Jz-Iz);
outIKxIL=(Ky-Iy)*(Lz-Iz) - (Ly-Iy)*(Kz-Iz);
% checking the above is enough to prove convexity. But the one
% below is needed to check whether points I,J,L are colinear!
outIJxIL=(Jy-Iy)*(Lz-Iz) - (Ly-Iy)*(Jz-Iz);
outJKxJL=(Ky-Jy)*(Lz-Jz) - (Ly-Jy)*(Kz-Jz);
% if they any product is less than zero, then IJKL have been defined
% clockwise or this is a non-convex quad. If any of them is zero
% then three points at least are colinear!
if outIJxIK<=0 | outIKxIL<=0 | outIJxIL <=0
warning(sprintf('Patch quad is non-convex or counter-clockwise defined or has at least 3 colinear points in line %d',i))
end
IJz=linspace(Iz,Jz,nIJ+1);
IJy=linspace(Iy,Jy,nIJ+1);
JKz=linspace(Jz,Kz,nJK+1);
JKy=linspace(Jy,Ky,nJK+1);
LKz=linspace(Lz,Kz,nIJ+1);
LKy=linspace(Ly,Ky,nIJ+1);
ILz=linspace(Iz,Lz,nJK+1);
ILy=linspace(Iy,Ly,nJK+1);
if fillflag
Z=zeros(nIJ+1,nJK+1);
Y=Z;
for j=1:(nIJ+1)
Z(j,:)=linspace(IJz(j),LKz(j),nJK+1);
Y(j,:)=linspace(IJy(j),LKy(j),nJK+1);
end
for j=1:nIJ
for k=1:nJK
patch([Z(j,k),Z(j,k+1),Z(j+1,k+1),Z(j+1,k)],...
[Y(j,k),Y(j,k+1),Y(j+1,k+1),Y(j+1,k)],...
matcolor{1+mod(mcolor,length(matcolor))},...
'edgecolor',edgecolor{1+mod(mcolor,length(edgecolor))})
end
end
else
for j=1:(nIJ+1)
plot([IJz(j) LKz(j)],[IJy(j) LKy(j)],'color',...
matcolor{1+mod(mcolor,length(matcolor))})
end
for j=1:(nJK+1)
plot([JKz(j) ILz(j)],[JKy(j) ILy(j)],'color',...
matcolor{1+mod(mcolor,length(matcolor))})
end
end
case 'layerstraight'
tmpcell=num2cell(gencom(i).nums);
[mcolor,nbars,barA,Iz,Iy,Jz,Jy]=deal(tmpcell{:});
Y=linspace(Iy,Jy,nbars);
Z=linspace(Iz,Jz,nbars);
Radius=sqrt(barA/pi);
%plot(Y,Z,'ko','markerfacecol','k')
% plot_circle is better because it captures accurately the area of
% the bar and does not rescale when we zoom in!
plot_circle(Y,Z,Radius*barRplotfactor,'','k')
case 'fiber'
warning('"fiber" is not supported yet!')
case 'layercirc'
warning('"layer circ" is not supported yet!')
case 'patchcirc'
warning('"patch circ" is not supported yet!')
end
%pause
end

% for box-girder of G9 bridge!
%axis([-15,1,-4,1])
%set(gca,'xdir','reverse')
axis equal
axis tight
set(gca,'box','on')
xlabel('z-coord (m)')
ylabel('y-coord (m)')

[divamva]

This way you will get them without the ruined formatting that I saw when posting in the forum!

http://www.ucy.ac.cy/~divamva/software/Wsection.tcl
http://www.ucy.ac.cy/~divamva/software/plotSection.m

[divamva]

Copy the stuff below to plot_circle.m, or simply download it from

http://www.ucy.ac.cy/~divamva/software/plot_circle.m
--------------------------------------------------------------------

function plot_circle(x,y,R,linestyle,facecolor,N)
% same as plot(x,y,'ko'), but plots small circles of radius R
% at each x,y point. Fully vectorized.
% Example: plot_circle([1 2 3 4], [4 5 5 6],[0.1 0.1 0.5 0.2])
%
% D.Vamvatsikos
% Updated: 8/Oct/05

if nargin<6 N=12; end
if nargin<5 facecolor=''; end
if nargin<4 linestyle='k-'; end
if nargin<3 R=repmat(1,size(x));
elseif length(R)==1 R=repmat(R,size(x));
elseif length(R)~=length(x)
error('R should be same length as x,y or scalar');
end

if length(x)~=length(y) error('x,y should be same lengths'); end
%keyboard

% use parametric form of circle
theta = linspace(0,2*pi,N)';
tx = repmat(cos(theta),1,length(x));
ty = repmat(sin(theta),1,length(x));

newx=repmat(x,N,1)+repmat(R,N,1).*tx;
newy=repmat(y,N,1)+repmat(R,N,1).*ty;

if isempty(facecolor)
plot(newx,newy,linestyle)
else
fill(newx,newy,facecolor)
end

[lightegg]

Hi! Thanks a lot.
what can we get by using plot_circle.m?

[kengawk]

thanks！
the scripts are very useful to visualize the os sections!

[linguan118]

Thank you!
Are you the author who proposed IDA?

[nahid2011]

Dear Sir

The bridge which I'm working on has a Hollow Rectangular RC section.

Would you please inform me that how I can use the scripts? Already I do not have MatLab on my PC, but I can find one.

Has your script written in a way that it can be easily fixed for RC section?

Shall I run OpenSees and MatLab simultaneously?

Best Regard.

Nahid

[pixazy]

Very useful to visualize the os sections!
Thanks

[thomas.ulrich]

Thanks (9 years later...) !

A Slight addon to be able to take into consideration rectangular patches.

Thomas.

function plotSection(sectionfname,barRplotfactor,fillflag,matcolor,edgecolor)
% plotSection(sectionfname,barRplotfactor,fillflag,matcolor,edgecolor)
%_________________________________________________________________________
% INPUT
% sectionfname : filename containing section info
% barRplotfactor: determines a magnification factor for plotting fibers
% produced via the "layer straight" command. Helps visualize
% the rebars (default = 1)
% fillflag : If 1, then the fibers are plotted as filled
% polygons/circles (default = 0)
% matcolor : a cell array containing color definitions for each
% material (default = {'y','b','r','g','m','k'} )
% edgecolor : same as above for the edges of the fibers. Used only when
% fillflag==1 (default = {'k'})
%
%__________________________________________________________________________
%
% Helps with taking FiberSection output (actually generation commands)
% and recreating the section fibers for visualization
%
% Examples
%__________________________________________________________________________
% plotSection('\users\divamva\opensees\divamva\Kavala\hollowRectfiberS1.out')
% plotSection('\users\divamva\opensees\divamva\Kavala\hollowRectfiberS1.out',2,1,{[0.7,0.7,0.7],[0.5,0.5,0.5]})
% plotSection('\users\divamva\opensees\divamva\Kavala\hollowRectfiberS1.out',2,1,{[0.7,0.7,0.7]},{[0.7,0.7,0.7]})
% plotSection('\users\divamva\opensees\divamva\RCsect\BGRCfiberS2.out')
% plotSection('\users\divamva\opensees\divamva\RCsect\RCsectionfiberS3.out')
%__________________________________________________________________________
% Created by D.Vamvatsikos
%
% Upgraded 07/May/05: Now the first line in the file is not discarded
% unless it contains an invalid command (e.g. a title)
% Also "support" is provided for "fiber", "layer circ" and
% "layer patch". A check is made to ensure "patch quad"
% coordinates create a valid convex and counterclockwise
% quadrangle!
% Upgraded 07/Oct/05: Added two extra checks for colinearity which were
% missing (must check all combinations of size 3)
%_____________________________________________________________________________
% NOTES:
% 1. Opensees will run even if a clockwise quad is supplied but then weird
% things will happen! Usually axial compression becomes tension and vice-versa
% 2. use view(90,90) to rotate the axes by 90deg and maintain nice labels!
% This is useful for plotting, e.g. the G9 or G7 sections.
% 3. The commands "fiber", "layer circ" and "layer patch" are not supported
% yet. They are only reckognized as being valid, but no plotting is
% done.


%cd e:\users\divamva\opensees\divamva\RCsect
if nargin<5, edgecolor={'k'}; end
if nargin<4, matcolor={'y','b','r','g','m','k'}; end
if nargin<3, fillflag=0; end
if nargin<2
% Do not increase bar radius plotting by default
barRplotfactor=1;
end

fid=fopen(sectionfname);

if fid>=0
% Since we do not use this "name" at all, refrain from reading it.
% If the user has not supplied such a name, then the valuable first
% line will be discarded!
%name = fgetl(fid);
i=0;
while 1
tline = fgetl(fid);
% check for EOF. Fully compatible with Matlab 7!
if ~ischar(tline), break, else i=i+1; end
% read fiber generation commands
[firstcom,count,err,ni]=sscanf(tline,'%s',1);
% only the "fiber" command has a single word as declaration.
% the rest use two words (e.g. "patch circ")
if strcmp(lower(firstcom),'fiber')
% one word was enough. Then numbers follow
gencom(i).com=firstcom;
else
% two words are needed.
[gencom(i).com,count,err,ni]=sscanf(tline,'%s',2);
end
gencom(i).nums=sscanf(tline(ni:end),'%f');
end
fclose(fid);
else
error(sprintf('File "%s" not found\n',sectionfname))
end

clf;
hold on

for i=1:length(gencom)
%lower(gencom(i).com)
switch lower(gencom(i).com)
case {'patchquadr','patchrect'}
tmpcell=num2cell(gencom(i).nums);
if strcmp(lower(gencom(i).com),'patchrect')
[mcolor,nIJ,nJK,Iy,Iz,Ky,Kz]=deal(tmpcell{:});
Jy=Ky;
Jz=Iz;
Ly=Iy;
Lz=Kz;
else
[mcolor,nIJ,nJK,Iy,Iz,Jy,Jz,Ky,Kz,Ly,Lz]=deal(tmpcell{:});
end
% check whether they are (1)counterclockwise and (2)form a convex
% quadrilateral in the y,z coord system
%
% z L* *K
% |
% | I* *J
% |
% +-----y
%
% We will achieve this by calculating the outer products of
% the vectors IJ x IK and IK x IL. If both are positive,
% we have convex and counter-clockwise quadrangle.
% Since all vectors have no x-component, the outerproduct
% is always along the x-axis and it is calculated as
% a x b = | ya za | * (unit x-vector)
% | yb zb |
% where ya,za are the vector coords (ya=ya2-ya1, za=za2-za1)
outIJxIK=(Jy-Iy)*(Kz-Iz) - (Ky-Iy)*(Jz-Iz);
outIKxIL=(Ky-Iy)*(Lz-Iz) - (Ly-Iy)*(Kz-Iz);
% checking the above is enough to prove convexity. But the one
% below is needed to check whether points I,J,L are colinear!
outIJxIL=(Jy-Iy)*(Lz-Iz) - (Ly-Iy)*(Jz-Iz);
outJKxJL=(Ky-Jy)*(Lz-Jz) - (Ly-Jy)*(Kz-Jz);
% if they any product is less than zero, then IJKL have been defined
% clockwise or this is a non-convex quad. If any of them is zero
% then three points at least are colinear!
if outIJxIK<=0 | outIKxIL<=0 | outIJxIL <=0
warning(sprintf('Patch quad is non-convex or counter-clockwise defined or has at least 3 colinear points in line %d',i))
end
IJz=linspace(Iz,Jz,nIJ+1);
IJy=linspace(Iy,Jy,nIJ+1);
JKz=linspace(Jz,Kz,nJK+1);
JKy=linspace(Jy,Ky,nJK+1);
LKz=linspace(Lz,Kz,nIJ+1);
LKy=linspace(Ly,Ky,nIJ+1);
ILz=linspace(Iz,Lz,nJK+1);
ILy=linspace(Iy,Ly,nJK+1);
if fillflag
Z=zeros(nIJ+1,nJK+1);
Y=Z;
for j=1:(nIJ+1)
Z(j,:)=linspace(IJz(j),LKz(j),nJK+1);
Y(j,:)=linspace(IJy(j),LKy(j),nJK+1);
end
for j=1:nIJ
for k=1:nJK
patch([Z(j,k),Z(j,k+1),Z(j+1,k+1),Z(j+1,k)],...
[Y(j,k),Y(j,k+1),Y(j+1,k+1),Y(j+1,k)],...
matcolor{1+mod(mcolor,length(matcolor))},...
'edgecolor',edgecolor{1+mod(mcolor,length(edgecolor))})
end
end
else
for j=1:(nIJ+1)
plot([IJz(j) LKz(j)],[IJy(j) LKy(j)],'color',...
matcolor{1+mod(mcolor,length(matcolor))})
end
for j=1:(nJK+1)
plot([JKz(j) ILz(j)],[JKy(j) ILy(j)],'color',...
matcolor{1+mod(mcolor,length(matcolor))})
end
end
case 'layerstraight'
tmpcell=num2cell(gencom(i).nums);
[mcolor,nbars,barA,Iz,Iy,Jz,Jy]=deal(tmpcell{:});
Y=linspace(Iy,Jy,nbars);
Z=linspace(Iz,Jz,nbars);
Radius=sqrt(barA/pi);
%plot(Y,Z,'ko','markerfacecol','k')
% plot_circle is better because it captures accurately the area of
% the bar and does not rescale when we zoom in!
plot_circle(Y,Z,Radius*barRplotfactor,'','k')
case 'fiber'
warning('"fiber" is not supported yet!')
case 'layercirc'
warning('"layer circ" is not supported yet!')
case 'patchcirc'
warning('"patch circ" is not supported yet!')
end
%pause
end

% for box-girder of G9 bridge!
%axis([-15,1,-4,1])
%set(gca,'xdir','reverse')
axis equal
axis tight
set(gca,'box','on')
xlabel('z-coord (m)')
ylabel('y-coord (m)')