Skip to main content

Finite Element Method Programming - 1D bar element, C# methods, C# loops

On the previous post, we established 1D bar finite element. This time we will go more into the programming in C# and try to make some automation in our code. I will try to explain:
- loops: for
- methods in C#
- arrays 1D and 2D

In the beginning, we have to build a task for us, to not feel that we are losing time. So, let's imagine that we want to check the result for the same bar, but with several divisions. For example: Does the result will change if I will create 2 finite elements instead of 1?

The beginning of the code will be looking the same:
  • namespace FEM1d_basic_1
  • {
  •     class Program
  •     {
  •         static void Main(string[] args)
  •         {
  •             // The heading of the code
  •             Console.WriteLine("FEM 1D program for bar");

  •             // The variables
  •             double emodulus = 210000; // 1000[kN/m2]=1GPa
  •             double barlength = 10; // [m] meters
  •             double elementarea = 0.01; // [m2] meter x meter 

  •             Matrix<double> ke = DenseMatrix.OfArray(new double[,] {
  •                 {1,-1},
  •                 {-1,1}
  •             });

  •             Console.WriteLine("Element stiffness matrix:");
  •             Console.WriteLine(ke);

The element stiffness matrix will be in general the same, no matter how many elements you will have in your structure. The issue is to aggregate them properly. So let's think for a second how the model looks with 2 finite elements.
You can easily find node number 2, which is shared by element 1 and element 2. The global sitffness matrix is 3x3 because in the global model we have 3 nodes. If we start to enumerate the nodes in the local stiffness matrices according to global names, we will see that we have to add value k[1,1] from the first element and k[0,0] from the second to aggregate global stiffness matrix.


It is easy to spot, that when you will have more elements in our model, the numbers will adding only on matrix diagonal. If we will have 5 elements and 6 nodes we will have 6x6 global matrix stiffness, which looks like that:

No! Wait I will not show you that, it is too boring to draw. I will now show you how to program this thing. First of all, we have to add a variable. It will be an integer and it will be saying to program how many divisions we want to create and ! according to what I said it will be also information how big the global matrix stiffness is.
To automatically create(initialize) the global stiffness matrix full of zeros we will write a small method. The method is the procedure, which can be used in several places in the code. Simpler, it is a recipe, it needs inputs and give you result or not (void method). Methods we are adding at the end of the class in which we are scripting.
  • namespace FEM1d_basic_1
  • {
  •     class Program
  •     {
  •         static void Main(string[] args)
  •         {
  •             ... //Our main code
  •         }

  •         public static double[,] IntializeSquareArray2D(int n) //method intializing square array full
  •             {
  •             double[,] squarearray = new double[n, n]; //creating instance of the array
  •             for (int i = 0; i < n; i++)         // first loop through raws
  •             {
  •                 for (int j = 0; j < n; j++)     // second loop through columns
  •                 {
  •                     squarearray[i, j] = 0;      // each i,j element will have value 0 (zero)
  •                 }
  •             }
  •                 return squarearray;  // the method will return only one thing, the array full of zeros!
  •             }
  •     }
  • }
The method is called IntializeSquareArray2D and it is creating a 2D array of dimension n x n. Since we want to have only zeros, we need just one input: n - dimension of this array. We should write then: 
  • IntializeSquareArray2D(int n)
  • {
  •    ... // The method code
  • }
But now what should be the output? Of course an array, so we need to add at the beginning:
  • double[,] IntializeSquareArray2D(int n)
  • {
  •    ... // The method code
  • }
In the code, we have also public and static. 
public - make method readable from the main code
static - make it easier to use, because we do not have to use instance constructors (will be later, forget now)
The first line of the code inside the  method is :
  • double[,] squarearray = new double[n, n]; 

It is creating an instance of the array, we say to the debugger: Hey, we need a place for an array with n x n elements. But we do not know what values will be there so we need an loop over columns and loop over raws to fill all of the cells. 
  •             for (int i = 0; i < n; i++)         // first loop through raws
  •             {
  •                 for (int j = 0; j < n; j++)     // second loop through columns
  •                 {
  •                     squarearray[i, j] = 0;      // each i,j element will have value 0 (zero)
  •                 }
  •             }
Now all of the elements in the array have got the 0 value. The last thing is to specified what exactly method is giving you back. To do this we use return command.

How we are using the method in the main code? We just type the name of it and of course, we are giving an input (remember the number of nodes is in our case is the number of elements + 1):
  •  int nelements = 2;
  •  Matrix<double> k = DenseMatrix.OfArray(IntializeSquareArray2D(nelements+1));
So we have the global stiffness matrix, but full of zeros. Could be worst... We now have to somehow put element stiffness matrices Ke inside the global stiffness matrix K. There is a lot of ideas how to do it. Right now I will just show you very specified to this example. Let's do a new method, what we need is the K, Ke, the starting node of the element, the ending node of the element.
  • public static Matrix<double> BuildGlobalStiffnessMatrix(Matrix<double> globalk, Matrix<double> localk, int n1, int n2)
  •         {
  •             n1 = n1 - 1;
  •             n2 = n2 - 1;
  •             globalk[n1, n1] = globalk[n1, n1] + localk[0, 0];
  •             globalk[n1, n2] = globalk[n1, n2] + localk[0, 1];
  •             globalk[n2, n1] = globalk[n2, n1] + localk[1, 0];
  •             globalk[n2, n2] = globalk[n2, n2] + localk[1, 1];
  •             return globalk;
  •         }
As you can see the method is overwriting the global stiffness matrix with the new elements in the asked coordinates (nodes n1, n2). It is very simple and to make it more useful we need to know which nodes are corresponding to which elements.
  •             // connectivity array
  •             int [,] connectivityarray= new int[nelements+1,3];

  •             for (int i1 = 1; i1 < nelements+1; i1++)
  •             {
  •                 connectivityarray[i1, 0] = i1;
  •                 connectivityarray[i1, 1] = i1;
  •                 connectivityarray[i1, 2] = i1+1;
  •             }  

  •             int node1 = 0;
  •             int node2 = 0;
  •             for (int i2 = 1; i2 < nelements+1; i2++)
  •             {
  •                 node1 = connectivityarray[i2, 1];
  •                 node2 = connectivityarray[i2, 2];
  •                 k = BuildGlobalStiffnessMatrix(k, ke, node1, node2);
  •             }
The connectivity array we will explain better in next post.
Here you can find the full code:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;

namespace FEM1d_basic_1
{
    class Program
    {
        static void Main(string[] args)
        {
            // The heading of the code
            Console.WriteLine("FEM 1D program for bar");

            // The variables
            double emodulus = 210000; // 1000[kN/m2]=1GPa
            //double npoisson = 0; // no units
            double barlength = 10; // [m] meters
            double elementarea = 0.01; // [m2] meter x meter 
            int nelements = 20;  //number of finite elements

            Matrix<double> ke = DenseMatrix.OfArray(new double[,] {
                {1,-1},
                {-1,1}
            });


            Console.WriteLine("Element stiffness matrix:");
            Console.WriteLine(ke);
            
            Matrix<double> k = DenseMatrix.OfArray(IntializeSquareArray2D(nelements+1));
            Console.WriteLine("Global stiffness matrix:");
            Console.WriteLine(k);

            // connectivity array
            int [,] connectivityarray= new int[nelements+1,3];

            for (int i1 = 1; i1 < nelements+1; i1++)
            {
                connectivityarray[i1, 0] = i1;
                connectivityarray[i1, 1] = i1;
                connectivityarray[i1, 2] = i1+1;
            }  

            int node1 = 0;
            int node2 = 0;
            for (int i2 = 1; i2 < nelements+1; i2++)
            {
                node1 = connectivityarray[i2, 1];
                node2 = connectivityarray[i2, 2];
                k = BuildGlobalStiffnessMatrix(k, ke, node1, node2);
            }
            
            Console.WriteLine("Global stiffness matrix:");
            Console.WriteLine(k);

            k=ApplyBoudaryConditionToNode(k,1);

            Console.WriteLine("Element stiffness matrix(after applying BC):");
            Console.Write(k);

            k= emodulus * elementarea / (barlength/nelements) * k;

            // Inverse stiffness matrix
            Matrix<double> ik = k.Inverse();

            Console.WriteLine("Inverse element stiffness matrix:");
            Console.Write(ik);

            // The force vector
            Vector<double> f = DenseVector.OfArray(IntializeArray1D(nelements+1)); //1 kN load in node 2
            f[nelements] = 1;

            Console.WriteLine("Force vector:");
            Console.Write(f);

            // The general equation Kr=R, deformation vector re
            Vector<double> re = ik * f;

            Console.WriteLine("Deformation vector:");
            Console.WriteLine(re);

            // The elongated beam
            double newbarlength = barlength + re[nelements];

            Console.WriteLine("new bar length=" + newbarlength);

            // The finish code
            Console.ReadKey();

        }

        public static double[,] IntializeSquareArray2D(int n) //method intializing square array full of zeros
            {
            double[,] squarearray = new double[n, n]; //creating instance of the array, we say to the debugger: Hey, we need place for an array with n elements
            for (int i = 0; i < n; i++)         // first loop through raws
            {
                for (int j = 0; j < n; j++)     // second loop through columns
                {
                    squarearray[i, j] = 0;      // each i,j element will have value 0 (zero)
                }
            }
                return squarearray;  // the method will return only one thing, the array full of zeros!
            }

        public static double [] IntializeArray1D(int n)
        {
            double[] columnarray = new double[n]; //creating instance of the array, we say to the debugger: Hey, we need place for an array with n elements
            for (int j = 0; j < n; j++)     // second loop through columns
            {
                columnarray[j] = 0;      // each j element will have value 0 (zero)
            }
            return columnarray ;
        }

        public static Matrix<double> BuildGlobalStiffnessMatrix(Matrix<double> globalk, Matrix<double> localk, int n1, int n2)
        {
            n1 = n1 - 1;
            n2 = n2 - 1;
            globalk[n1, n1] = globalk[n1, n1] + localk[0, 0];
            globalk[n1, n2] = globalk[n1, n2] + localk[0, 1];
            globalk[n2, n1] = globalk[n2, n1] + localk[1, 0];
            globalk[n2, n2] = globalk[n2, n2] + localk[1, 1];
            return globalk;
        }

        public static Matrix<double> ApplyBoudaryConditionToNode(Matrix<double> globalk, int n1)
        {
            n1 = n1 - 1;

            for (int i = 0; i < globalk.ColumnCount; i++)
            {
                globalk[n1, i] = 0;
            }

            for (int j = 0; j < globalk.RowCount; j++)
            {
                globalk[j, n1] = 0;
            }

            globalk[n1, n1] = 1;

            return globalk;
        }
    }
}

Comments

Post a Comment

Popular posts from this blog

How to write C# components in Visual Studio for Grasshopper? If you have the Visual Studio community 2017, you have to also download assembler for grasshopper v6. https://marketplace.visualstudio.com/items?itemName=McNeel.GrasshopperAssemblyforv6 After downloading the file, start it. The VS community 2017 should be chosen automatically. Now you can start VS 2017 one more time, in the new project, you should find new options: If you choose Visual C# > Rhinoceros > Grasshopper Add-on for v6, a new window will appear: The missing files have to be added to references. The RhinoCommon and Rhino.exe you will probably find somewhere on C disc (in my case it is C:\Program Files\Rhinoceros 5 (64-bit) ). The Grasshopper.dll file location you should check on your disc through searching enginge ( in my case it was :  C:\Program Files\Common Files\McNeel\Rhinoceros\5.0\Plug-ins\Grasshopper (b45a29b1-4343-4035-989e-044e8580d9cf)\0.9.76.0 ). Click finish and new empty code sh
Post a Comment
Read more
The programming skills starts to be essential for nowadays engineer.  The  open source software create a great chance for developing automation in engineering, but also demand some programming knowledge from the structural engineers. Here I will try to slowly introduce engineers in to the basic programming knowledge. The main focus will be on the Finite Element Method in the parametric environment, but also the approach to some Eurocodes will be shown. The main programming language will be c# and python. The software in which I will be working will be the Grasshopper (add-on for Rhino) and Autodesk Dynamo. To start our journey with programming I have to explain the main types which you will find in codes: string - the text type, you can store letters, words or sentences in it. For example: string text1 = "architect is not engineer" ; Let's try to understand this declaration of variable in C#. First we are giving the type ( string ) then we are setti
Post a Comment
Read more