One Stop All Solutions

Latex for Bloggers: Problems & Solutions

noreply@blogger.com (Husnain) — Sun, 12 Jun 2011 03:46:00 +0000

A few days back a visitor commented that the my equations on this blog are not showing in Google Chrome. It turned out that it was not working in any browser, let alone Google Chrome. I checked out the java script and everything seemed OK. Then I surmised that one of the servers that the script refers to might be down. I checked out and found the following to be not working in the mathtex3.js:

http://www.cyberroadie.org/cgi-bin/mathtex.cgi?.

Now one of the solutions could be the replacement of the above server with a valid working server such as the following:

http://mathcache.appspot.com/

But I chose a more simpler solution.

The Solution

The solution was quite simple. I replaced the java script location

https://sites.google.com/site/onestopallsolutions/course-materials/mathtex3.js?attredirects=0&d=1

with the working server's address

http://mathcache.appspot.com/

in the java script for latex equations. The final script looked like as shown below:

<script src="http://mathcache.appspot.com/" type="text/javascript">
</script>
<script type="text/javascript">
onload = function () {
 var els = document.getElementsByTagName("div");
 var i = els.length;
 while (i--) {
  if (els[i].className=="latexparse") {
   replaceMath(els[i]);
  }
 }
};
</script>

And now the equations are showing once again.

Setting Up Java Script for Latex in Blogger

If you are wondering how to setup a java script for latex in blogger, then you should read the post Latex for Bloggers: Watchmath not Working . It explains setting up scripts for latex equations in details.

There you go. See you next post.

Video Tutorial : How To Implement Finite Element Method 4

noreply@blogger.com (Husnain) — Mon, 23 May 2011 01:04:00 +0000

In the previous post I discussed the procedure for finding out the global stiffness matrix through local elements. In this post we are going to take this process further, develop the entire linear system, elaborate on finding out the force vector and finally solve the system by applying boundary conditions.

Linear System: Ku=F

We have already seen how to formulate the global stiffness matrix. This was one part of the entire linear system which is also called the left hand side. Traditionally, the unknowns are kept on the left hand side and the known values are put on the right hand side. In the animation which follows, the same right hand side has been discussed. This is nothing but the force vector. In all our animations, we have tried to keep things simple by assuming parameters to be constants. The same has been applied in this animation. The objective here is to build an intuition. Complex cases would be addressed later.

Let us now see this animation:

Solving the System: Boundary Conditions Application

Once the linear system is developed, the next step is to solve it (of course!). But this step also requires the application of boundary conditions. We will see in what follows that without boundary conditions, the linear system is not resolvable. The matrix is rank deficient and physically this means that the system is not constrained and so any external stimulus may cause it to undergo rigid body motion, however, it will certainly not deform.

Let us now look at the animation below to see how the boundary conditions are applied and how the linear system is resolved.

In the animation above, we also showed the stresses in the bar,

$$\sigma=E\frac{du}{dx}\ $$

Since we approximated the deformation by linear functions, therefore, the derivatives of those linear functions would be a constant. The stresses can still be recovered by joining the mid points of the constant slopes over all the elements.

Instructions for Viewing the Animations Full Screen

Here is how to play this animation full screen. Right click and save file as a gif file. If you do not already have storm codec Media Player Classic then download it using this link. Open the saved gif file with this media player and play and pause as you like. If you do not play this gif file with media player classic then you can not pause the frames and some frames might be too fast to read properly. This animation could also be played with Quicktime Player. You will find it here.

See you next post.

Video Tutorial : How To Implement Finite Element Method 3

noreply@blogger.com (Husnain) — Tue, 10 May 2011 09:06:00 +0000

From this post we saw how to formulate one row of the stiffness matrix. But this is not how you want the computer to calculate it. In this post we are going to explore, through some animations, how to calculate the entire global stiffness matrix by assembling local or elementary matrices. Please make sure to read the instructions at the end to view the animations full screen and to have play and pause control.

Think Globally Act Locally

Ultimately you would want to program your approximate solution in some computing language. The first step is to convert the global system (often a large system) into a local framework. This is exactly, think globally act locally kind of a thing. To further elucidate, let me give you another quote:

the longest journey begins with a single step

Let me further explain the transformation with the help of the following animation:

Element Stiffness Matrix Calculation

Now that we know how to isolate a single element and what are the elemental conventions, we could calculate the element stiffness matrix. The animation below shows us how:

You must have noticed that in order to simplify the calculations, I assumed that $EA(x)=1$. The intuition behind this simplification is the fact that $EA$ is a material property and if that is equal to a constant then the stiffness matrix entries are only dependent upon the shape functions. This makes it easy to think about other problems, without knowing a lot about their physics. For example for the one dimensional electrostatic boundary value problem, the element stiffness matrix is given as:
$$k_{ij}=\int_{x_{1}^{l}}^{x_{2}^{l}} \varepsilon^{l}(x){N^{l}_{j}}^{\prime}{N^{l}_{i}}^{\prime} dx $$

Does it not look similar. It sure does. The only difference between the above expression and the one we saw in the animation is that $EA(x)$ is now replaced by the permittivity of the medium $\varepsilon (x)$, where the former is the product of the Young's modulus of elasticity $E$ and the area of the bar $A$ and the latter is the product of dielectric constant of the medium $\varepsilon_{r}$ and permittivity of the free space $\varepsilon_{o}$ ($\varepsilon=\varepsilon_{r}\varepsilon_{o}$). So for a case where $\varepsilon(x)=1$, there is no difference between the element stiffness matrices, provided the element sizes $h^{l}$ are the same.

Assembling the Element Matrices

The next step is to assemble all the local matrices in the global system. The first thing we need for this is the element connectivity table. An element connectivity table tells us how different elements are connected and which nodes connect which elements. For a one dimensional case, every element is linear and is connected by two nodes on either side. In turn each node connects two elements on either side, except for the first and the last node. To further explain this, the first node connects the first element only, the first element is connected by the first and the second node, and the second node connects the first and the second element. This means that the second node is common to both the first and the second element. Same is true for all the elements until the last element arrives. The last node only connects the last element and that is it.

Once we have the connectivity table, properly identifying the inner nodes and the boundary nodes and all the elements therein, we can move on to assemble the local matrices into a global one. What better way to do that than to show you another small animation.

Each row has three non zero entries which is coherent with what was said about the global stiffness matrix for a simple one dimensional case with linear elements. However, the first and the last rows contain only two non-zero entries. The reason is that the first and the last nodes connect to one element only.

Instructions for Viewing the Animations Full Screen

The instructions to watch the animation are repeated here from this post. Here is how to play this animation full screen. Right click and save file as a gif file. If you do not already have storm codec Media Player Classic then download it using this link. Open the saved gif file with this media player and play and pause as you like. If you do not play this gif file with media player classic then you can not pause the frames and some frames might be too fast to read properly. This animation could also be played with Quicktime Player. You will find it here.

See you next post.

Latex for Bloggers: Watchmath not Working

noreply@blogger.com (Husnain) — Sat, 07 May 2011 08:54:00 +0000

Alright a little diversion. For those of you who have seen previous posts know that I use a lot of latex equations. You can too it is simple.

Writing Latex Equations

Here are a few simple steps which you could find on a number of blogs. I am sorry if I do not remember those tens of blogs / sites to give them due credits.

1. Go to the design tab of your blogger dashboard.

2. Then click the Page Elements.

3. By default you should have your blog post on the left and some free space on the right (with some default gadgets already included). This should be the layout of your blog page.

4. On the right hand side of your page, click add a gadget. Another window will pop up. Scroll down in this window until you fine HTML/Java Script. Click on the plus sign and it will be added.

5. In the configure HTML/Java Script window copy and paste the following code:

<script src="http://www.watchmath.com/cgi-bin/mathtex3.js" type="text/javascript"></script>
<script type="text/javascript">
replaceMath( document.body );</script>

6. Now write your equations as shown below:

$ y=mx+c $
$$ y=\frac{1}{1+e^{x}}\ $$;

7. The result will be as follows:

$y=mx+c$
$$y=\frac{1}{1+e^{x}}\ $$

Watchmath not Working

Seldom you find out that the site which is hosting the javascript is down. For example the watchmath site which was hosting the latex javascript is down. Click this link and you will know. There, however, is a way around. Simply host the script at your site. Simply go to google sites and get yourself one for free. And if you do not want to do that, simply copy and paste the following code into the HTML/Java Script gagdet of your blogger and you are done.

<script src="https://sites.google.com/site/onestopallsolutions/course-materials/mathtex3.js?attredirects=0&d=1" type="text/javascript"></script>
<script type="text/javascript">
replaceMath( document.body );</script>

Known Issues

Now, if you know latex, you can easily write beautiful mathematical equations, just by enclosing them between single $ or double $$ signs. But what if you literally want to write the $ sign. The solution to this is proposed by Declan1991 at this site. Copy and paste the following code instead:

<script src="https://sites.google.com/site/onestopallsolutions/course-materials/mathtex3.js?attredirects=0&d=1" type="text/javascript">
</script>
<script type="text/javascript">
onload = function () {
 var els = document.getElementsByTagName("div");
 var i = els.length;
 while (i--) {
  if (els[i].className=="latexparse") {
   replaceMath(els[i]);
  }
 }
};
</script>

Now, whenever you want to write latex equations just enclose your argument between the $ signs as shown below:

< div class="latexparse"> $your latex argument here$ < /div>

There you go. See you next post.

Video Tutorial : How To Implement Finite Element Method 2

noreply@blogger.com (Husnain) — Fri, 06 May 2011 07:10:00 +0000

In the previous post I discussed some fundamental concepts in implementing the finite element method. In today's post we will move further. The post ends with a small animation on the calculation of the stiffness matrix.

Stiffness Matrix

For those of you who have a background in engineering, probably know Hooke's Law. It states that the restoring force in the direction opposite and equal to the direction of the applied force is directly proportional to the displacement of the material under the applied load. Mathematically this is stated as:

$f=-ku$

The negative sign simply says that the restoring force is in the opposite direction of the applied force. The constant of proportionality $k$ is called the stiffness of the material. Its magnitude is easily calculated by:
$$k=\frac{f}{u}\ $$

Now why are we discussing it here in the context of finite elements. Recall from this post that the weak formulation is given by:
$$ \int_{0}^{l}EA(x) \frac{du}{dx}\ \frac{dv}{dx}\ dx = \int_{0}^{l}f(x)v(x)dx + \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$

In the previous post we also looked at the formulation of the basis functions and how these basis functions are used in the approximation of finite element method. From the same post let us now write the following:
$$u(x)=\sum_{j=1}^{n}u_{j}\phi_{j}$$
$$v(x)=\phi_{i}$$

This substitution is proposed by Galerkin and hence this kind of finite element method is also known as Galerkin Finite Element Method. Now what is special about this Galerkin method? The special thing about this is that the trial and the test function are defined using the same basis functions. Now remember we were trying to figure out how the term stiffness is involved in finite elements. Let us now substitute the test and the trial functions in the weak form:
$$ \sum_{j=1}^{n}u_{j} \int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx= \int_{0}^{l}f(x)\phi_{i}dx $$
$$+ \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$

You may notice that we did not replace $u(x)$ in the term $$ \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$. The reason is that this term is the boundary term and is only calculated at $x=0$ and $x=l$. For an overhanging bar (fixed at one end and free at the other) this term is zero at $x=0$ and is equal to $P\phi_{j}$ at $x=l$, where $P$ is the applied load at the free end. Substituting this results in:
$$ \sum_{j=1}^{n}u_{j} \int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx= \int_{0}^{l}f(x)\phi_{i}dx + \left. P \phi_{i} \right|_{x=l}$$
Writing the above equation in matrix form gives the following:
$$K_{ij}u_{j}=F_{i}\qquad(a)$$
Where $$F_{i}=\int_{0}^{l}f(x)\phi_{i}dx + \left. P \phi_{i} \right|_{x=l}$$ is everything on the right hand side and
$$K_{ij}=\int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx$$.

Now since equation (a) resembles the Hooke's Law, hence the matrix $K$ is called the stiffness matrix. This is the mystery behind calling it the stiffness matrix and it remains as such even if Poisson's Equation is resolved for electric field between two conducting plates.

Stiffness Matrix Calculation

Let us now see how the stiffness matrix is calculated. This is facilitated with the help of a small animation. The directions to play the animation full screen are given just after.

Now that we have an expression of the stiffness matrix, let us see a small animation which explains how to calculate it.

Q1. Why is the second integral zero in the $i-1^{st}$ column?

A. The second integral as asked in the animation is zero because the function $\phi_{i+1}$ is zero in the range $x_{i-1}\leq x \leq x_{i}$. Pause at the frame when this question is posed in the animation. You could only pause by playing the gif file in media player classic.

Q2. Why is the first integral zero in the $i+1^{st}$ column?

A. The first integral as asked in the animation is zero because the function $\phi_{i-1}$ is zero in the range $x_{i}\leq x \leq x_{i+1}$. Pause at the frame when this question is posed in the animation. You could only pause by playing the gif file in media player classic.

This is continued until the next post.

Tutorial : How To Implement Finite Element Method

noreply@blogger.com (Husnain) — Mon, 18 Apr 2011 13:27:00 +0000

In the previous post I discussed some of the concepts in finite element method. It is important to note that although the case that we have taken up is that of longitudinal deformation of a bar, the findings and discussions are equally applicable to any problem which can be solved by finite element method. Examples of application will be presented later, this post, however deals with the general steps in implementing the method.

Meshing

The first step is to discretize the physical domain into smaller elements. For a one dimensional problem the discretization is very simple, and is shown with the help of Fig. 1 below:

Fig. 1

The discretization is done by cutting the bar into smaller linear elements. These elements are labelled as $I_{i}$ and the subscript $i$ denotes the element number. Just like the bar, these smaller elements also have end points which are called nodes of the elements. These are labelled as $x_{i}$ where the subscript $i$ denotes the node number. The nodes are shown with big dots and they mark the start and end of different elements. Note that the number of elements is one less than the number of nodes. The elements are linear, of course, because the system under consideration is one dimensional. This process of breaking up of the physical domain into smaller pieces is known as meshing.

Basis Functions

You may notice that there are coloured triangles over each of the nodes with the exception of the boundary nodes which do not have complete triangles. Moreover, these triangles are labelled as $\phi_{i}$ where $i$ corresponds to the node number. What are these triangles? These triangles are nothing but interpolation functions also called Lagrangian basis functions. Their job is to provide a linear approximation of the original solution $u(x)$. This approximation is given as: $$u_{FE}(x)= \sum_{i}u_{i}\phi_{i}(x)$$, where $u_{i}$ is the nodal amplitude.The writing of this equation implies that we are looking for a series solution for the original function $u(x)$.

Let us see through Fig. 2 below, how these basis functions which are also known as hat functions are defined.

Fig. 2

A certain basis function $\phi_{i}$ is depicted in Fig. 2. At the nodal coordinate $x_{i}$ to which this basis function belongs, its value should be 1. It can be confirmed from the figure that at $x=x_{i}$ the value of $\phi_{i}$ is 1. Basis functions are defined separately for different segments, as shown in Fig. 2. These functions (hat functions) are non-zero only in the two adjacent elements connecting the node on which they are defined. They are zero everywhere else in the domain.

The topic is continued till the next post.

Tutorial : Finite Element Method Fundamentals Cont...

noreply@blogger.com (Husnain) — Sat, 16 Apr 2011 14:29:00 +0000

In the last post we had gotten as far as the weak formulation of the problem. We had also posed a few questions. Their answers will be given in this post. The discussion will be general though and applicable to all the problems amenable to finite element approximation.

Weighted Residual!!

What is meant by weighted residual? Remember from the previous post, the residual of the system were all the terms moved to the left hand side of the differential equation. The right hand side then had no terms left and it was obviously equal to zero. In other words, the residual was equal to zero. If we multiply this residual by some weighing function, the resultant is still zero (no surprises!). Now the idea of a weighted residual formulation is to multiply the residual by some weighing function and then to integrate this product over the entire domain of the problem. For a one dimensional problem, this is equivalent to integrating over the length of the domain. The integral is then put equal to zero. Of course, because the residual is zero.

$$ \int_{\Omega}r(x)w(x)dx=0$$

Here $\Omega$ is nothing but the physical domain of the problem. For a one dimensional problem and in particular for the problem of longitudinal deformation of a bar, this is the length of the bar.

A very delicate point in all this process is that since we are approximating a function, the residual is not strictly equal to zero; owing to approximation errors. To further elaborate, let us reconsider the expression of the residual function:
$$r(x)= -\frac{d}{dx}\left ( EA\frac{du}{dx} \right) - f(x)$$
Since we know that the original function $u(x)$ would be replaced by some finite element trial function $\hat{u}(x)$, which only approximates the original function as close as possible; therefore, there will be a difference in the two solutions and hence the residual is not strictly equal to zero. This brings us to another very important point.

"Since the residual is not strictly equal to zero, the whole idea of a finite element solution is to find so good an approximation that the residual nearly comes out to be zero and so that the weighted residual integral becomes valid if it is put equal to zero."

Test Function vs Weighing Function

In the previous post I wrote that the weighing function is deliberately omitted and that the test function $v(x)$ is used instead. This point is not a big deal. The weak form of the problem remains the same no matter which of the two are used. The only point worth mentioning here is that somewhere in the history we started using the following convention:

$u$ for the trial function
$v$ for the test function
$w$ for the weight function

One way of obtaining finite element formulation is to derive the weak form using the weight function $w(x)$ and then substitute for it the test function $v(x)$. When this step is done, it is said that $v(x)$ is an arbitrary admissible function. One key property of this admissible function is that when the original function is supposed to be zero, the test function is also zero. This property comes in handy when applying boundary conditions. See the previous post for boundary condition application.

Why Call it Weak Form?

The next step is to apply integration by parts to the weighted residual formulation or the formulation with the test function, and render the weak form of the integral. Why call it 'weak'? It is called a weak form because in the original integral the unknown function is differentiated twice, and the weighing function or the test function is without any differential operator. However, once integration by parts is applied, both the trial and the test functions have the same number of differential operators.
$$\int_{\Omega} EA(x) \frac{du}{dx} \frac{dv}{dx} dx = \int_{\Omega}f(x)v(x)dx + \left. EA \left( \frac{du}{dx} \right) v \right |_{0}^{l}$$
This process of taking off one differential operator of the original trial function reduces the minimum smoothness criterion. In other words, for the above integral to be finite only the first derivative of the trial function $u(x)$ needs to be defined. If on the other hand we had tried solving the original integral then we would have needed the second derivative of $u(x)$ to be defined also. This is viewed as weakening of the requirement of some higher order derivative to be defined in the integral and hence, this form is called the weak form of the problem. To see the complete steps, refer to my previous post.

In the next post we will discuss the fundamental steps in realization of finite element approximation.

Tutorial : Finite Element Method Fundamentals

noreply@blogger.com (Husnain) — Fri, 15 Apr 2011 21:24:00 +0000

hello again,

In the previous post I left you guys with a simple code which gives a one dimensional finite element approximation to longitudinal deformation of a bar. Let us learn about the problem in more details.

Problem Definition and its Differential Equation

The problem is explained with the help of the following diagram

Fig. 1

As shown in Fig. 1, a bar fixed at one end is subjected to uniform axial force $f_{o}$ and an end load $P$. At the fixed end the deformation $u$ is zero. The objective is to obtain an expression for the deformation $u$ as a function of $x$. To elaborate the implementation of FEM on this simple problem, let us consider its differential equation.

This equation in terms of the unknown function $u(x)$ is valid in the domain $0 < x< l$ .

Boundary Conditions Application

The boundary conditions given at the two extremities of $x$ are:

Since the bar is fixed at $x=0$, it is very easy to see that there is no deformation and consequently $u(x)$ is zero at this point. This is called as the "Dirichlet" type of boundary condition. At the other end, no deformation is imposed and consequently the bar is free to deform under the applied load. It can be observed from Fig. 1 that there is an end load present. This end load which is a point force can be expressed in terms of the first derivative of the unknown deformation function $u(x)$. This expression is demonstrated by the second boundary condition and it also goes by the name "Neumann" type of boundary condtion.

Obtaining the Weak Form

The next step is to find the weak form of the differential equation. Let $r(x)$ be the residual of the differential equation given by:

$$r(x)=-\frac{d}{dx} \left( EA(x)\frac{du}{dx} \right ) - f(x)$$

Let us now introduce a test function $v(x)$ which is an admissible function and behaves more or less like the original trial function $u(x)$. Note that we are intentionally skipping the part where the admissible function $v(x)$ is termed as weight function $w(x)$ and multiplied by the residue $r(x)$. This omitted step has the name "weighted residual". However, we are going to follow almost identical steps as with weighted residual without having to use it explicitly. Since, the residual is zero we may write;

$$\int_{0}^{l}\left \{ -\frac{d}{dx} \left ( EA (x)\frac{du}{dx} \right ) - f(x) \right \} v(x) dx = 0$$

The most important step in FEM formulation is the application of integration by parts. Let us apply this technique first and then see why was it necessary.

$$\int_{0}^{l}EA(x)\frac{du}{dx} \frac{dv}{dx} dx = \int_{0}^{l}f(x)v(x) + \left. EA(x)\frac{du}{dx} v(x) \right |_{0}^{l}$$

The second term on the right hand side is the boudary term. At $x=0$, $u=0$ and therefore $v$ is also zero. This is one of the assumptions for the test function to be admissible as mentioned earlier. At $x=l$:

$$\left. EA(x)\frac{du}{dx} v(x) \right |_{x=l}=\left. Pv(x)\right |_{x=l}$$

Simplifying and substituting the boundary conditions, the weak form is finally given as:

$$\int_{0}^{l}EA(x)\frac{du}{dx} \frac{dv}{dx} dx = \int_{0}^{l}f(x)v(x) +\left. Pv(x)\right |_{x=l}$$

This topic is continued until the next post. Till then take care.

Tutorial : 1D Finite Element Method Matlab Code

noreply@blogger.com (Husnain) — Mon, 27 Dec 2010 22:39:00 +0000

Hello Everybody,

This is my first post. What better way to do this than to elaborate the title of my blog. I must confess, I had a tough time coming up with a suitable title, primarily because my interests are very diverse and there were not too many words which could reflect what this blog is all about. I finally managed to come up with this title ‘One Stop All Solutions’. This implies that this would be a single place where you would get almost everything regarding engineering, data mining, acquisition and IT. I know this is exaggerated, but this is what it literally means.

On a more practical note, this blog is about my professional experience in the industry, my passion for IT and my recent affair with data mining.

I have done my PHD in vibrations and acoustics. I taught too while I was doing my PHD and I love sharing knowledge in general, so, through this blog, I will write articles about numerical and computational techniques, give you short MATLAB codes (occasionally Python or C / C++ too if possible) on diverse fields that I have been involved with (machine learning, Bayesian inference, Genetic Algorithms, Data Acquisition, FEM etc.) and not too often I will bore you with anecdotes of my life too, before they are eternally buried in the sands of time.

To end my first post, I leave you with a Matlab code for a very simple 1D FEM problem. I will elaborate on this in the next post. Till then, take care.

function fem_1D
% This is a simple 1D FEM program. The FEM
% solution is based on linear elements also called hat functions.The
% problem addressed is the extension of a bar under the action of applied
% forces. Use mesh parameters under the heading mesh of this code to change
% values. For example change the number of nodes to 2 to really see the
% difference between the exact and the FEM solution. In general, change
% constants the way you like.
%
% The second plot of stresses in the bar suggests that for each of the
% finite elements in the bar the solution (that is the slope of the
% extension) is a constant. This is illustrated with the help of horizontal
% lines. Each line also represents one element. It would be shown in later
% codes that the slope of the approximate solution (stresses) coincides
% with the slope of the exact solution exactly at the mid point of the
% element.

% This code is inspired from a manual written by Jack Chessa on FEM
% implementation in MATLAB.

% Written 26-Dec-2010
% Author: Husnain Inayat Hussain (husnain.inayat@gmail.com)

% FEM_1D

% 1. PreProcessor
% 2. Processor
% 3. PostProcessor

% ---------------- PREPROCESSOR -------------------

% Preprocessor consists of defining the problem, material, elemental
% discretization, geometry, applied conditions such as forces and
% displacements.

% PROBLEM

% MATERIAL
E = 1; % Young's Modulus
A = 1; % Area of the bar
L = 1; % Length of the bar
f_o = 1;
u_o =0;
P = 10;

% MESH
num_nodes = 3;
num_elems = num_nodes - 1;
nodes = linspace (0,L,num_nodes);
elm_cnc = [1:1:num_elems ; 2:1:num_nodes]';
h_e = nodes (2) - nodes (1);

% ----------------   PROCESSOR   -------------------

% Stiffness Matrix
K_e = E*A/h_e * [1 -1;
                -1 1];
K = zeros (num_nodes, num_nodes);

F_e = f_o*h_e/2;
F = zeros (num_nodes, 1);

for i = 1 : size(elm_cnc,1)
    asm = elm_cnc (i,:);
    K(asm , asm) = K(asm , asm) + K_e;
    F(asm) = F(asm) + F_e;
end

F = F - K(:,1) * u_o;
F (1) = u_o;
F (end) = F (end) + P;

K(1,:) = 0;
K(:,1) = 0;
K(1,1) = 1;

u = K\F;

% ----------------   POSTPROCESSOR   -------------------

% Plot the displacements and the stressed of the analytical and FEM
% solutions.

% FEM Solution & Plot
plot(nodes, u,'Marker','o','LineWidth',1,'LineStyle','--',...
    'Color',[1 0 0],...
    'DisplayName','FEM');

% Analytical Solution & Plot
x = linspace(0,L,100);
u1 = - f_o/(2*E*A)*x.^2 + (P + f_o*L)*x;
hold
plot(x, u1,'LineWidth',1,'DisplayName','Analytical');
h = legend('FEM','Analytical',2);
set(h,'Interpreter','none')
grid on

% Create title
title({'Extension of a Longitudianl Bar','Analytical vs.FEM Solution'},...
    'FontWeight','normal',...
    'FontSize',12);

% Create xlabel
xlabel({'Lenght (m)'},'FontWeight','normal','FontSize',11);

% Create ylabel
ylabel({'Extension (m)'},'FontWeight','normal','FontSize',11);

% Plot the stresses of the analytical and FEM solutions

% FEM Solution & Plot
figure
du = zeros (5,1);
plot(nan,nan)
hold

h1 = zeros (size(elm_cnc,1),1);

for i = 1 : size(elm_cnc,1)
    du(i) = (u(i+1) - u (i))/ (nodes(i+1) - nodes(i));
    asm = elm_cnc (i,:);
    h1(i) = plot (nodes(asm), repmat(du(i),1,length(nodes(asm))),...
        'Marker','o','LineWidth',1,'LineStyle','--','Color',[1 0 0],...
        'DisplayName','FEM');
end

% Analytical Solution & Plot
du1 = - f_o/(E*A)*x + (P + f_o*L);
h2 = plot (x, du1, 'LineWidth',1,'LineStyle','-',...
    'Color',[0 0 1], 'DisplayName','Analytical');
h = legend([h1(end) h2], 'FEM','Analytical', 1);
set(h,'Interpreter','none')
grid on

% Create title
title({'Stresses in a Longitudianl Bar','Analytical vs.FEM Solution'},...
    'FontWeight','normal',...
    'FontSize',12);

% Create xlabel
xlabel({'Lenght (m)'},'FontWeight','normal','FontSize',11);

% Create ylabel
ylabel({'Stress (N/m^2)'},'FontWeight','normal','FontSize',11);