Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Medical Image Segmentation with 3D Slicer: A Beginner's Guide

3D Slicer is an open source software that is widely used for image processing, visualization, filtering, and so forth. In this article, I would like to discuss its effectiveness and application in medical image segmentation. This software is completely FREE, and there is a continuous development of many add-ins/extensions for this platform. You may download the software from 3D Slicer Webpage. After downloading, the installation procedure is straightforward, you may also follow the instructions mentioned in their webpage.

Medical image segmentation to extract the size or volume of an organ or complex airways/channels from computed tomography (CT) or micro-computed tomography (𝛍CT) is very interesting and has been playing a crucial part in biomedical engineering. For example, human nasal cavities or airways have such a complex formation that from the CT scans, we are unable to extract the volume. However, with the power of modern image processing tools, we can manually trace the nasal airways path slice by slice and finally obtain the volume of the cavities.

In this tutorial, I am going to highlight the 3D Slicer's interface, and how we can do the manual/automatic segmentation for the purpose of extracting the volume or shape information. When you launch the software, you will see the following welcome window. The window has four sectional views. The three black windows show the transverse (or, axial), sagittal (or, longitudinal), and coronal (or, frontal) views of the CT images. The other blue colored section highlights the three dimensional reconstruction of the segmented images.

3D Slicer Interface after Launching
















Medical images are often processed as .DICOM file, which can be easily opened in 3D Slicer. In the welcome window, there are options to load different format of data, but Digital Imaging and Communications in Medicine (DICOM) images are smoothly handled by the 3D Slicer. In the following, I have opened a CT scan of a human brain with a tumor, which is available in 3D Slicer's sample data for demonstration. Our task here is to find the volume of the tumor quantitatively from the CT scan images so that physicians can have more in-depth information regarding the size and growth of the tumor. They can make their next decisions confidently based on these data.

Loading the CT images in 3D Slicer
















To do so, at first, we need to open the segment editor, which is located at the top of the window where you will see an icon says 'All Module'. Click on it and then, select the segment editor as highlighted below.

Segment Editor of the 3D Slicer
















Once you select the segment editor, a new small window will pop-up just like the screenshot shown beneath. You will find some features in that new little window, such as, 'ADD', 'Show 3D' etc.

Adding a New Segment to Start the Segmentation
















You need to choose the option 'Add (+)' and then, a new segment will be created named as 'Segment_1' (see below). From the screenshot, you see that a new segment is ready to carry on for further operation.

A New Segment Created in the Property Editor
















Now, if you double click on the 'Segment_1', you will see the following options will appear to select the segmentation nature. For example, if we segment a bone, then choose 'bone' in that window; or if the segmentation is related to the brain or tooth, select the respective 'brain' and 'tooth' segmentation features. In our example, it is a tumor that we are going to segment. A tumor is a mass of abnormally growth of tissues, so we can choose the segmentation media as 'tissue'.

Selection of the Type or Nature of Segmentation
















So, we have decided how we are going to segment the tumor, now it's time to move forward using the tools to do the segmentation. There are two types of segmentations: automatic and manual. To understand the difference between the automatic and manual segmentation, let me first explain how the CT scans are formatted or processed for image analysis. There are many single images that form the CT scan of an organ, which alternatively implies that a CT scan itself an image that has subsets of many images. So, when we consider automatic segmentation, it is a process of tracing or narrowing a particular section from the CT scan that is done by the software, and manual segmentation is that you need to paint or trace that particular region by hand for all the individual images that eventually form the CT scan. As you see below, the red circled area is the tumor that we need to separate out and get the volume of it. We also see that there are some options, like 'Paint', 'Draw', 'Level Tracing' to trace the tumor. Here, I used the 'Level Tracing' option for that, traced few images or slices and then used the 'Fill between slices' option as an automatic segmentation approach.

Use of the Segmentation Tools

Once the segmentation is done, then we can estimate the volume of the tumor. For that, select again 'All Modules', and from the dropdown menu, click on the 'Quantification', and then hit 'Segment statistics'.


How to Extract Volume Information from 3D Slicer
















After that, at the bottom of your window, you will find the volume of the segmented tumor. From the screenshot below, we see the volume of the tumor is approximately 346.99 mm³. Now, remember, this is not an accurate or quantitatively perfect representation of the tumor. However, it gives a reasonable view and idea about the shape and size of the tumor. Below, we can see the looks of the tumor in a three dimensional space.

The 3D Shape Extracted from Segmentation with Volume Statistics



#ComputedTomography #3DSlicer #ImageSegmentation #ImageProcessing #3DModelling #Blog #Blogger

A Taylor Series Expansion to Approximate a Function near a Point

The Taylor series expansion up to seventh terms for approximation of near point x =1 is,
Taylor series expansion up to seventh terms





Letting starting point xₒ is at 0 and the step size (xi - xₒ) is 1. Now, the value of eˣ at x = 1 is to be calculated by Taylor series expansion on the basis of the starting point and its derivatives at x = 0 with step size 1. A MATLAB program is written to evaluate this value from zero order approximation to sixth order approximation and show the improvement of accuracy and minimization of truncation error as well.
The true value of is 2.718281828459046. The truncation error is the difference between true and approximate values. The following MATLAB program describes each steps of the calculation. 

Taylor Series Implementation by a MATLAB Program

% Taylor Series expansion to approximate exp(x) near point x=1
 
function e = taylorseries(x0,x1,h)
 
e0=exp(x0);              % calculates exp(x) value when x=0
 
derivative1=1;           % all the six derivatives calculated at x=0
derivative2=1;
derivative3=1;
derivative4=1;
derivative5=1;
derivative6=1;
 
% zero order approximation
 
        e=e0;   
        error=exp(x1)-e;      
        disp (e);        % display of exp(x) and error at command prompt          
        disp(error);
        
% first order approximation        
            
        e=e0+(derivative1*h);       
        error=exp(x1)-e;       
        disp (e);       
        disp(error);
             
% second order approximation
 
        e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% third order approximation
                     e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% fourth order approximation        
                  e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)));       
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% fifth order approximation                      
        e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)))+(derivative5*(h^5/factorial(5)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% sixth order approximation
        e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)))+(derivative5*(h^5/factorial(5)))+(derivative6*(h^6/factorial(6)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
end


Table 1: Outputs from the Above Program to Approximate at x=1

Approximation of an exponential function by Taylor series










The above table tells us as the order of the Taylor series increases the truncation error also decreases. However, after sixth and higher order approximations, improvement of the accuracy is not much significant. So, we can stop the approximation to first seven terms of the Taylor series and the approximate value of ex near point x = 1 is 2.718.

Finding Roots of a Equation by Newton-Raphson Method

In this tutorial, we are interested to the roots of the following function.

f(x) = sin x + 2 - 0.1 x

To find the value of x or roots of the equation, Newton-Raphson method is applied. For the proper initial guess for x, graphical method is applied to visualize the characteristics of this nonlinear function. This graphical technique also helps to find multiple roots for the unknown function. The following MATLAB commands generate the plot of the function, f(x). This function is plotted with respect to the values of x from -50 to +50.

>> X=-100:100;
>> plot(X,(sin(X)+2-(0.1*X)))

From the following figure 1, it is clear that the equation has seven roots and plotting the function helps to get the roots. The function changes its sign from, x = 10 - 11, 11 -12, 16 -17, 18 - 19, 21 - 22, 25 - 26, and 27 - 28. And, the roots of the equation must lie between these intervals. Now, to get good approximations of the roots, a Newton-Raphson algorithm is written in MATLAB.

Plot of the function f(x) to identify the intervals where roots of the function lie















Figure 1: Plot of the function f(x) to identify the intervals where roots of the function lie.

The program is developed by MATLAB user defined function called “sinusoidal”. The initial assumptions for approximating seven roots of x are 10.5, 11.5, 16.5, 18.5, 21.5, 25.5 and 27.5 respectively. The codes are given below. It is straightforward to run the program seven times as we have seven initial guesses for the corresponding roots.

Newton-Raphson Algorithm by MATLAB

% Definition of a function "sinusoidal" to solve equation by Newton-Raphson
function [x,ea] = sinusoidal(X, etol)
format long;
 
% Input: X=21.5, etol=Tolerance definition for error
% Output: x=Solution of equation, ea=Calculated error in loop
 
% Program Initialization
 
solution = 21.5;
 
% Iterative Calculations
 
while (1)
    solutionprevious = solution;

% Implementation of Newton-Raphson Method 
           
    solution = X - (sin(X)+2-(0.1*X))/(cos(X)-0.1);               
    X=solution;       
           
if solution~=0                         % Approximate percent relative error
                                 
    ea=abs((solution - solutionprevious)/solution)*100;      
end
 
if ea<=etol,                  % Condition to meet specified error tolerance
   
    break,                
end
 
x = solution;
 
disp(x);
 
end

% Display of output parameters
 
disp(x);                      
disp(ea);
 
end

Program Outputs: 

>> sinusoidal (10.5,1e-10)
ans =
  10.636782424281707
>> sinusoidal(11.5, 1e-10)
ans =
  11.562055865585652
>> sinusoidal (16.5,1e-10)
ans =
  16.107752970689450
>> sinusoidal (18.5,1e-10)
ans =
  18.721338781142521
>> sinusoidal (21.5,1e-10)
ans =
  21.809224297055675
>> sinusoidal (25.5,1e-10)
ans =
  25.744697403517950
>> sinusoidal (27.5,1e-10)
ans =
  27.435908953388580

So, the seven roots from the program are 10.637, 11.562, 16.108, 18.721, 21.809, 25.745 and 27.436. In all the cases, the error tolerance limit is kept extremely small (0.0000000001) to get more accurate and precise results.


#Blog #Blogger #NewtonRaphson #Matlab #Sinusoidal

Explicit Euler Method to Solve System of ODEs in MATLAB

In this tutorial, I am going to show a simple way to solve system of first order ordinary differential equations (ODE) by using explicit Euler method. Let' say, we have following three first order ODEs.

First Order Ordinary Differential Equations









Here, we have 3 ODEs, 3 dependent variables (x, y, z), and 1 independent variable, t. The initial conditions when time is zero are, x(0) = y(0) = z(0) =1. Our goal is to solve these differential equations with explicit Euler approach and plot the solutions afterwards.

Step 1: Define the Equations
The first step is to define all the differential equations in MATLAB. I did this by using MATLAB function handle, which is shown below.

Defining differential equations in MATLAB












Step 2: Choose a Numerical Approach 
The next step is to select a numerical method to solve the differential equations. In this example, we will use explicit Euler method. I have created a function to implement the algorithm. The following image shows the application of the explicit Euler method.

Application of the explicit Euler method in MATLAB
















Step 3: Call the Function
This is the final step where we need to call the function. In the previous step, we have created a function, which we are going to call in this step to solve the equations.

Final stage where we call the function



 











So, that's it. We are done with the process, and we can now visualize the solution. I have also attached the MATLAB codes for this problem at the end.

Plot showing the solutions of differential equations














MATLAB Program:

close all;

clc;
format long;

% Defining the three differential equations of the problem
f = @(t,y) [-5*y(1)+5*y(2); 14*y(1)-2*y(2)-y(1)*y(3); -3*y(3)+y(1)*y(2)];
[x,y] = explicit_euler(f,[0,5],[1;1;1],0.001); % Calling Euler function
plot(x,y);
title('When Time Step is 0.001');
legend('x(t)', 'y(t)', 'z(t)', 'Location', 'NorthEast')
xlabel('t')
ylabel('Solutions')


function [x, y] = explicit_euler( f, xRange, y_initial, h )
% This function uses Euler’s explicit method to solve the ODE
% dv/dt=f(t,v); x refers to independent and y refers to dependent variables
% f defines the differential equation of the problem
% xRange = [x1, x2] where the solution is sought on
% y_initial = column vector of initial values for y at x1
% numSteps = number of equally-sized steps to take from x1 to x2
% x = row vector of values of x
% y = matrix whose k-th column is the approximate solution at x(k)
x(1) = xRange(1);
numSteps = ( xRange(2) - xRange(1) ) /h ;
y(:,1) = y_initial(:);
for k = 1 : numSteps
x(k + 1) = x(k) + h; 
y(:,k+1) = y(:,k) + h * f( x(k), y(:,k) );
end



#ExplicitEulerMethod #Matlab #NumericalMethod #Blog #Blogger

Discussion on Numerical Integration Approaches with MATLAB

In this tutorial, we are going to implement four numerical integration schemes in MATLAB. The methods are:

i)   Midpoint Rectangle Rule
ii)  Trapezoidal Rule
iii) Simpson’s 1/3 Rule
iv) Simpson’s 3/8 Rule

Let's say, we would like to integrate the following error function where the upper interval, a is 1.5. At first, we will evaluate the true integral of the error function, and then, we will find the integral numerically with the above methods (i) -(iv).

Error function integral




Finally, we will compare each method with the true value and make some comments based on the deviations if we have any. For each method, we will find out the following parameters:

Table showing integration parameters




Here, n refers to the numbers, 1,2,4,8,16, ...,128. I_num is the numerical value of integration, h is the step size, and E is the error, which is the absolute difference between true and numerical results. First, we are going to determine the true value of the error function while integrating from 0 to 1.5. Then, we will apply the corresponding methods and compare the absolute error.

MATLAB Code for Midpoint Rectangle Rule

close all;
clear;
clc;
 
% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);
 
fprintf('\n     n          h         I_NUM          E')
 
for j = 1:L
    h(j)=(b-a)/n(j);
   
x_1 = a:h(j):b;
% Using MATLAB built-in command for repeated arrays
x_2 = repelem(x_1(2:end-1),2); 
xh = [a x_2 b];
 
% Numerical integration by mid-point rectangular rule
I_NUM=0;
for i = 1:length(x_1)-1
    I_NUM = I_NUM + f(a  +(i-1/2).*h).*h;
end
 
% Defining the error
E(j) = abs(I-I_NUM(j));
 
fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
 
end

The following table shows the parameters determined by this approach:
Show comparison with true integration value












MATLAB Code for Trapezoidal Rule

close all;
clear;
clc;
 
% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);
 
fprintf('\n     n          h         I_NUM          E')
 
for j = 1:L
    h(j)=(b-a)/n(j);

% Numerical integration using Trapezoidal rule
I_NUM=0;
for i = 1:n(j)
    I_NUM = I_NUM+(f(a+(i-1).*h) + f(a+i.*h)).*h/2;
end
 
% Defining the error
E(j) = I - I_NUM(j);
 
fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
end

The following table shows the parameters determined by this approach:
Show comparison with true integration value












MATLAB Code for Simpson’s 1/3 Rule

close all;
clear;
clc;
 
% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);
 
fprintf('\n     n          h         I_NUM          E')

for j = 1:L    
    h(j) = (b-a)/(2*n(j));  
 
% Numerical integration using Simpson’s 1/3 rule
 
I_NUM = 0;
 
for i = 1:n(j)
    I_NUM = I_NUM+(f(a+(2*i-2).*h) + 4*f(a+(2*i-1).*h) + f(a+(2*i).*h)).*h./3;
end
E(j) = abs(I - I_NUM(j));
 
fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
end

The following table shows the parameters determined by this approach:
Showing comparison with true integration value












MATLAB Code for Simpson’s 3/8 Rule

close all;
clear;
clc;
 
% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);
 
fprintf('\n     n          h         I_NUM          E')

for j = 1:L  
h(j) = (b-a)/(3*n(j)); 
 
% Numerical integration using Simpson’s 3/8 rule
I_NUM = 0;
for i=1:n(j)
    I_NUM = I_NUM + (f(a+(3*i-3).*h) + 3*f(a+(3*i-2).*h) + 3*f(a+(3*i-1).*h) + f(a+(3*i).*h))*3.*h/8;
end
 
% Defining the error
E(j) = abs(I - I_NUM(j));
 
fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
 
end

The following table shows the parameters determined by this approach:
Showing comparison with true integration value













We see from the above comparisons that both Simpson's rules give us a good approximation of the true value, especially when n is greater than 4. For the rectangular approach, the value reaches close to the true value when n is 128. The same is also evident for the trapezoidal rule. So, we can conclude that both the Simpson's rules are more efficient and accurate approaches compared to the rest.



#TrapezoidalRule #SimpsonsRule #Matlab #NumericalIntegration #Blog #Blogger

Stress and Strain Analysis of Linear Finite Elements

In this tutorial, we will use principle of virtual work method to determine the stress and displacement of the nodes of linear one dimensional finite elements. Following figure 1 represents a cantilever beam which is stressed axially. Let's assume the parameters L, A, E, c have the same value of 1 unit for simplicity. Our goal is to compare the finite element results with the exact solutions of a typical cantilever problem.
For this problem, we need to consider two quadratic one dimensional element. So, the displacement is interpolated among the three nodes and has the following form:
displacement


Cantilever Beam
Fig.1: Showing one dimensional axially load bar given in the assignment problem.

















piecewise linear interpolation
Fig.2: Showing one dimensional piecewise linear interpolation.













Mathematica Program for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]]

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

(*Shape Functions*)
Nn = {N1, N2, N3}

(*Applied Axial Load*)
p = c*x;

(*Nodal Force Vectors*)
fe = Integrate[Nn*p, {x, 0, L/2}];
c = 1;
L = 1;
fe // MatrixForm

Program Outputs

Displacement Vector for Element 1:

Displacement Vector




Shape Functions for Element 1:

Shape Functions




Nodal Force Vector for Element 1:

Nodal Force Vector






Mathematica Program for Element 2

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]]

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

(*Shape Functions*)
Nn = {N1, N2, N3}

(*Applied Axial Load*)
p = c*(x + L/2);

(*Nodal Force Vectors*)
fe = Integrate[Nn*p, {x, 0, L/2}];
c = 1;
L = 1;
fe // MatrixForm

Program Outputs

Displacement Vector for Element 2

Displacement Vector




Shape Functions for Element 2

Shape Functions




Nodal Force Vector for Element 2

Nodal Force Vector






Local Stiffness Matrix for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L, EA]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]];

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

B = {{D[N1, x], D[N2, x], D[N3, x]}};

(*Local Stiffness Matrix*)
K = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

L = 1;
EA = 1;
K // MatrixForm

Program Outputs

Program Outputs






Comparison between the Exact and Solution using Four Elements

We know that the exact solution for the axially loaded bar shown in Fig. 1 may be written as,

exact solution




Where, for the consecutive four elements, we may write the displacement for each element as,

each element displacement











From the assignment problem, the parameter values are,

parameter values




From the virtual work method, we may write,

virtual work method





The internal virtual work (IVW) may be expressed as,

internal virtual work




And, the external virtual work may be presented as,

external virtual work




MATHEMATICA PROGRAM

The following program is written to determine the nodal displacement using the finite element method considering the principle of virtual work. 


Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]];

(*Shape Functions*)
N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];
Nn = {N1, N2, N3};

B = {{D[N1, x], D[N2, x], D[N3, x]}};

c = 1;
L = 1;
EA = 1;

(*Axial Force on Element 1*)
p1 = c*x;
fe1 = Integrate[Nn*p1, {x, 0, L/2}];

(*Axial Force on Element 2*)
p2 = c*(x + L/2);
fe2 = Integrate[Nn*p2, {x, 0, L/2}];

K1 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

K2 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

Ndof = 5;
Ne = 2;
ID1 = {1, 2, 3};
ID2 = {3, 4, 5};

K = ConstantArray[0, {Ndof, Ndof}];
K[[ID1, ID1]] = K1 + K[[ID1, ID1]];
K[[ID2, ID2]] = K2 + K[[ID2, ID2]];

fe = ConstantArray[0, {Ndof, 1}];
fe[[ID1]] = fe1 + fe[[ID1]];
fe[[ID2]] = fe2 + fe[[ID2]];

IDE = {1};
IDN = {2, 3, 4, 5};

KR = Delete[K, IDE];
KR = Delete[Transpose[KR], IDE];
KR = Transpose[KR];

feR = Delete[fe, IDE];

(*Calculation of Nodal Displacements*)
DisR = LinearSolve[KR, feR];

Dis = ConstantArray[0, {Ndof, 1}];
Dis[[IDN]] = DisR + Dis[[IDN]];
Dis // MatrixForm

Program Outputs

Program Outputs









So, we see that the exact and finite element solutions are the same.

From the following two plots, we notice that although the nodal displacement is continuous from the FEA result, it is not true for the stress. The stress distribution is indeed discontinuous.

nodal displacement











stress distribution


MATHEMATICA PROGRAM FOR COMPARISON

(*Exact Solution*)
uexact = c*L^2/2/EA*x - c/6/EA*x^3;

(*Finite Element Approximation*)
N1 = Piecewise[{{4 x/L, 0 <= x <= L/4}, {4/L (L/2 - x), 
     L/4 <= x <= L/2}}, 0];
N2 = Piecewise[{{4/L (x - L/4), L/4 <= x <= L/2}, {4/L (3 L/4 - x), 
     L/2 <= x <= 3 L/4}}, 0];
N3 = Piecewise[{{4/L (x - L/2), L/2 <= x <= 3 L/4}, {4/L (L - x), 
     3 L/4 <= x <= L}}, 0];
N4 = Piecewise[{{4/L (x - 3 L/4), 3 L/4 <= x <= L}}, 0];

u1 = c*L^3/EA*(47/384);
u2 = c*L^3/EA*(11/48);
u3 = c*L^3/EA*(39/128);
u4 = c*L^3/EA*(1/3);

(*Calculating Nodal Displacement*)
u = u1*N1 + u2*N2 + u3*N3 + u4*N4;
uexactp = uexact /. {c -> 1, EA -> 1, L -> 1}

sigmaexact = D[uexactp, x]
up = u /. {c -> 1, EA -> 1, L -> 1}
sigmap = D[up, x]

(*Plotting to Show the Comparison between Exact and FEA Solution*)
Plot[{uexactp, up}, {x, 0, 1}, 
 AxesLabel -> {"x", "Nodal Displacement"}, 
 PlotLegends -> {"uexact", "u_Finite Element"}]
Plot[{sigmaexact, sigmap}, {x, 0, 1}, AxesLabel -> {"x", "sigma_11"}, 
 PlotLegends -> {"sigma_exact", "sigma_Finite Element"}]


For further information on finite element analysis and Mathematica, readers are encouraged to visit this page.


#FiniteElement #StressStrain #Mathematica #Abaqus #FEA #Blog #Blogger