Clear Filters
Clear Filters

I want to shorten this code but dont dont know how???

3 views (last 30 days)
function [theta, M] = assignment(fc,fy,As,xl,xu,yl,yu,nx,tol,thetaf,m)
% Input data, data can be either for under reinforced or over reinforced
function [fc,fy,As,xl,xu,yl,yu,nx,tol,thetaf,m] = data
fc = 30; % Mean compressive strength of concrete
fy = 415; % Yield strength of steel
As = [324 50 50; 324 150 50]; % [(Cross sectional area of rebar) (x-coordinate of rebar) (y-coordinate of rebar)]
xl = 0; % Lower limit of x-coordinate
xu = 200; % Upper limit of x-coordinate
yl = 0; % Lower limit of y-coordinate
yu = 500; % Upper limit of y-coordinate
nx = 200; % No of grid points in each direction
tol = 1e-3; % Tolerance for finding neutral axis
thetaf = 0.08/1000; % Final value of curvature
m = 100; % Number of increments
end
% Stress strain curve for concrete
function sigma= sigmac (e, fc)
ec2 = 2.0 / 1000 ;
ecu = 3.5 / 1000 ;
n = 2 ;
if e < 0
sigma = 0 ;
else
if e > ecu
sigma = 0 ;
elseif e > ec2
sigma = fc ;
else
sigma = fc * (1 - (1 - e / ec2)^n) ;
end
end
end
% Stress strain curve for steel
function sigma = sigmas (e, fy)
Es = 200000;
esy = fy / Es ;
esu = 10 / 1000 ;
if abs (e) > esu
sigma = 0 ;
elseif abs (e) > esy
sigma = fy * e / abs (e) ;
else
sigma = e * Es ;
end
end
% Finding the geometry of cross-section
function I = geometry(x)
xv = [0 200 200 0];
yv = [0 0 500 500];
I = inpolygon(x(1),x(2),xv,yv);
end
% Moment curvature diagram
% {Evaluating curvature considering the cross-section and eliminating other
% point which are outside the cross-section}
function [theta,M] = diagram(fc,fy,As,xl,xu,yl,yu,nx,tol,thetaf,m)
A0 = (xu-xl)*(yu-yl);
n = nx^2;
w = A0/n;
nr = size(As,1); % Number of rebars
du = 1/nx;
u = (du / 2): du: (1-du / 2) ;
k = 0 ;
for i =1: nx
for j = 1: nx
k = k + 1 ;
S(k,:) = [ u(i) u(j) ] ;
end
end
S (:, 1) = xl + S (:, 1) * (xu - xl) ;
S (:, 2) = yl + S (:, 2) * (yu - yl) ;
for i =1:n
I (i) = geometry (S(i,:)) ;
end
S(I ==0,:) = [ ] ;
ns = size (S, 1) ;
dtheta = thetaf / m;
for i =2: m
theta (i) = i * dtheta ;
M(i) = momrot (theta(i), ns, S, w, nr, As, fc, fy, yl, yu, tol) ;
end
end
% Finding the position of neutral axis
function M = momrot (theta, ns, S, w, nr, As, fc, fy, yl, yu, tol)
a = yl ;
b = yu ;
Na = resultants (theta, ns, S, w, nr, As, a, fc, fy) ;
h0 = b - a;
it = ceil (log2 (h0 / tol)) ;
for i = 1: it
c = 0.5 * (a + b) ;
Nc = resultants (theta, ns, S, w, nr, As, c, fc, fy) ;
if Na * Nc < 0
b = c ;
else
a = c ;
Na = Nc ;
end
end
yc = 0.5 * (a + b) ;
[ N, M ] = resultants (theta, ns, S, w, nr, As, yc, fc, fy) ;
end
% Finding the moment
function [ N, M ] = resultants (theta, ns, S, w, nr, As, yc, fc, fy)
N = 0 ;
M = 0 ;
for i =1: ns
y = S(i, 2) ;
e = tan(theta) * (y - yc) ;
sigma = sigmac (e, fc) ;
N = N + w * sigma ;
M = M + w * sigma * (y - yc) ;
end
for i =1: nr
y = As (i, 3) ;
e = tan(theta) * (y - yc) ;
sigma = sigmas (e, fy) ;
N = N + sigma * As (i, 1) ;
M = M + sigma * As (i, 1) * (y - yc) ;
end
end
[ fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m ]=data ;
[ theta, M ] = diagram(fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m) ;
max(M);
disp(M);
plot (theta, M);
xlabel("Curvature (1/mm)");
ylabel("Moment (Nmm)");
end
  1 Comment
Mathieu NOE
Mathieu NOE on 20 Sep 2023
why not simply store individually your many functions into a toolbox (must be declared on your path)
then your main code needs only to call (in one line) your functions
no need to have them permanently present at the end of your script

Sign in to comment.

Accepted Answer

Shreshth
Shreshth on 5 Jan 2024
Hello Kundan,
I understand that you want to shorten the given code. This can be achieved through several ways. Some of them I have mentioned below.
1.Vectorization: MATLAB is highly optimized for vector and matrix operations. Wherever possible, loops should be replaced with vectorized operations. This can significantly reduce the length of the code and often improve performance.
Original version:
for i =1: ns
y = S(i, 2);
e = tan(theta) * (y - yc);
sigma = sigmac(e, fc);
N = N + w * sigma;
M = M + w * sigma * (y - yc);
End
Vectorized version:
e = tan(theta) * (S(:, 2) - yc);
sigma = arrayfun(@(strain) sigmac(strain, fc), e);
N = sum(w * sigma);
M = sum(w * sigma .* (S(:, 2) - yc));
2.Inline Functions: MATLAB allows the creation of anonymous functions (lambdas) for simple operations. These can be used to replace standalone functions that are only used once or have a very simple body.
Original standalone function:
function sigma = sigmac(e, fc)
% ... function body ...
end
Inline version:
sigmac = @(e, fc) (e >= 0 & e <= 3.5/1000) .* fc * (1 - (1 - e / (2.0/1000)).^2);
3.Removing Redundant Code: If there are checks or initializations that are unnecessary or can be simplified, they should be removed or combined.
Original code checking for stress in concrete:
if e < 0
sigma = 0;
else
% ... more conditions ...
end
Simplified version (assuming e is always non-negative in context)
sigma = fc * (1 (1 e / ec2).^n) .* (e <= ec2) + (e > ec2 & e <= ecu) * fc;
4.Combining Functions: If functions are only called once and are specific to the problem at hand, their functionality can sometimes be combined into a single function to reduce the overhead of multiple function calls.
5.Use Built-in Functions: MATLAB has a wide range of built-in functions that can often replace custom code. For example, the ‘inpolygon’ function can be used to test if points are inside a polygon, replacing custom geometry functions.
Original custom geometry check:
function I = geometry(x)
% ... function body ...
End
Using inpolygon:
I = inpolygon(S(:,1), S(:,2), [0, 200, 200, 0], [0, 0, 500, 500]);
Using the above methods will reduce your code length significantly. Additionally, you can also try the following ways to reduce the code length further:
  • Simplifying Conditional Statements:
  • Preallocating Arrays
  • Use of ‘fzero’ for Root Finding
Hope the information is helpful.
Thank you,
Shubham Shreshth.

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!