Finite Element Method Programming in C#  CodeProject
: 19
Table of Content
Introduction
This article is about BriefFiniteElement.NET, an open source library for analyzing solids and structures with FEM in C# hosted on Codeplex. Development of library and documentation (including this article on Codeproject) is ongoing.
You can check the library features from Features of BriefFiniteElement.NET
Download
As every several days the source code is changing (because of adding new features or fixing bugs) so I suggest you to download the source code and simply build them to be sure that you have latest binaries of BriefFiniteElement.NET.
Check the Downloading Latest Source Code for BriefFiniteElement.NET for instructions on simply downloading the latest source code of project.
Before starting the examples, you have to add a reference to BriefFiniteElement.NET
project or dll file.
Example 1: A Simple Truss
In this example, I want to analyse a simple truss with 4 members as shown in picture.
All members sections are the same, a square steel section with dimension of 3 cm. So the properties of members will be:
E = 210 GPa = 210*10^9 Pa = 210*10^9 N/M^2
A = 0.03m * 0.03m = 9 * 10^4 m^2
We should do 5 steps before we solve the model:
 Create Model, Members and Nodes.
 Add the Nodes and Elements to Model
 Assign geometrical and mechanical properties to Elements
 Assign Constraints to Nodes (fix the DoF s)
 Assign Load to Node
And finally solve model with Model.Solve()
method and then extract analysis results like support reactions or member internal forces or nodal deflections.
Creating Model, Members and Nodes
Creating Model
We should create a Finite Element model first and then add members and nodes to it:
// Initiating Model, Nodes and Members
var model = new Model();
Creating Nodes
We should create nodes like this. In BriefFiniteElement.NET
, every node and element have a property of type string
named Label
and another one named Tag
both of which are inherited from BriefFiniteElementNet.StructurePart
. In every Model
, Label
of every member should be unique among all members (both Nodes and Elements) unless the Label
be equal to null
which is by default. In the below code, we are creating 5 nodes of truss and assigning a unique Label to each one.
var n1 = new Node(1, 1, 0);
n1.Label = "n1";//Set a unique label for node
var n2 = new Node(1, 1, 0) {Label = "n2"};//using object initializer for assigning Label
var n3 = new Node(1, 1, 0) {Label = "n3"};
var n4 = new Node(1, 1, 0) {Label = "n4"};
var n5 = new Node(0, 0, 1) {Label = "n5"};
Creating Elements
Then we have to create the elements. In BriefFiniteElement.NET
, the TrussElement
class represents a truss element in 3D.
var e1 = new TrussElement2Node(n1, n5) {Label = "e1"};
var e2 = new TrussElement2Node(n2, n5) {Label = "e2"};
var e3 = new TrussElement2Node(n3, n5) {Label = "e3"};
var e4 = new TrussElement2Node(n4, n5) {Label = "e4"};
2 Adding Nodes and Elements to Model
You can simply add the elements and nodes we created into the Model
. Model
has two members Elements and Nodes which both represents an IList<T>
of nodes and members, plus an AddRange
method.
model.Nodes.Add(n1, n2, n3, n4, n5);
model.Elements.Add(e1, e2, e3, e4);
Please note that if Node
or Element
’s Label
property is something else than null
, then it should be unique among all nodes and elements.
3 Assigning geometrical and mechanical properties to Elements
As elastic module for all members equals to 210 GPa and area of all members equals to 0.0009 m^2 we can set the element properties like this:
e1.A = e2.A = e3.A = e4.A = 9e4;
e1.E = e2.E = e3.E = e4.E = 210e9;
4 Assigning Constraint to Nodes to fix the DoFs
Now, we should make some DoFs of structure fix in order to make analysis logically possible.
In BriefFiniteElement.NET
, every node has 6 degree of freedom: X, Y, and Z rotations and X, Y, and Z translations. For a every truss model, we have to fix rotational DoFs for each Node (X,Y and Z rotation). Also the nodes 1 to 4 are also movement fixed, then nodes 1 to 4 should be totally fixed and node 5 should be rotation fixed. In BriefFiniteElement.NET
, a struct
named Constraint
represents a constraint that is applicable to a 6 DoF node, it have Dx, Dy, Dz, Rx, Ry and Rz properties of type DofConstraint
which is an enum
and have two possible values 0 (Released) and 1 (Fixed). For making work easier, the Constraint struct
has some predefined Constraint
s in its static
properties for example Constraint.Fixed
or Constraint.Free
. Here is more detailed information:
Property name 
Description 

All 6 DoFs are fixed 

All 6 DoFs are released 

3 translation DoFs are fixed and 3 rotation DoFs are released 

3 translation DoFs are released and 3 rotation DoFs are fixed 
We can fix DoFs of nodes 1 to 4 like this:
n1.Constraints = n2.Constraints = n3.Constraints = n4.Constraints = new Constraint(dx:DofConstraint.Fixed, dy:DofConstraint.Fixed, dz:DofConstraint.Fixed, rx:DofConstraint.Fixed, ry:DofConstraint.Fixed, rz:DofConstraint.Fixed);
or:
n1.Constraints = n2.Constraints = n3.Constraints = n4.Constraints = Constraint.Fixed
and should fix the rotational DoFs of node 5:
n5.Constraints = Constraint.RotationFixed
5 Assigning Loads to Nodes
In BriefFiniteElement.NET
, there is a struct
named Force
which represent a concentrated force in 3D space which contains of 3 force components in X, Y and Z directions and three moment components in X, Y and Z directions. It have 6 double properties named Fx, Fy, Fz, Mx, My and Mz that are representing the load components. There are also two properties of type Vector
for this struct
named Forces
and Moments
. On setting or getting, they will use the Fx, Fy, Fz, Mx, My and Mz to perform operations:
/// <summary>
/// Gets or sets the forces.
/// </summary>
/// <value>
/// The forces as a <see cref="Vector"/>.
/// </value>
public Vector Forces
{
get
{
return new Vector(fx,fy,fz);
}
set
{
this.fx = value.X;
this.fy = value.Y;
this.fz = value.Z;
}
}
Same is with Moments
property. The Forces
and Moments
property do not actually store values in something other than 6 obvious properties.
As LoadCase
and LoadCombination
concepts are supported in BriefFiniteElement.NET
, every Load
should have a LoadCase
. A LoadCase
is simply a struct
that has two properties: CaseName
with string
type and LoadType
with LoadType
type which is an enum
and has some possible values:
public enum LoadType
{
Default = 0,
Dead,
Live,
Snow,
Wind,
Quake,
Crane,
Other
}
The LoadType.Default
is a load type that is created for built in usage in library and it do not meant to have meaning like Dead
, Live
, etc. The LoadCase struct
has a static
property named LoadCase.DefaultLoadCase
:
/// <summary>
/// Gets the default load case.
/// </summary>
/// <value>
/// The default load case.
/// </value>
/// <remarks>
/// Gets a LoadCase with <see cref="LoadType"/> of <see cref="BriefFiniteElementNet.LoadType.Default"/> and empty <see cref="CaseName"/></remarks>
public static LoadCase DefaultLoadCase
{
get { return new LoadCase(); }
}
Which represents a LoadCase
with LoadType
of Default
and CaseName
of null
. We will call such a LoadCase
as DefaultLoadCase
. For simplicity of usage in BriefFiniteElement.NET
everywhere that you’ll prompt for a LoadCase
, if you do not provide a LoadCase
then the LoadCase
is assumed DefualtLoadCase
by the library. For example, when you want to assign a load to a node, you should provide a LoadCase
for it, like this:
var load = new NodalLoad(new Force(0, 0, 1000, 0, 0, 0), new LoadCase("Case1",LoadType.Dead));
but if you do not provide the LoadCase
in the above code like this:
var load = new NodalLoad(new Force(0, 0, 1000, 0, 0, 0));
then the load case will be assumed DefaultLoadCase
by the library.
Ok, next we have to add 1KN load to node 5 like this, will do it with DefaultLoadCase
:
var force = new Force(0, 0, 1000, 0, 0, 0);
n5.Loads.Add(new NodalLoad(force));//adds a load with LoadCase of DefaultLoadCase to node loads
And finally solve the model with model.Solve()
method. Actually solving the model is done in two stages:
 First stage is creating stiffness matrix and factorizing stiffness matrix which will take majority of time for analyzing
 Second phase is analyzing structure against each load case which takes much less time against first stage (say for example 13 sec for first stage and 0.5 sec for second stage).
First stage is done in model.Solve()
method and second stage will done if they’ll be need to.
There are loads with different LoadCase
s that are applied to the Nodes and Elements. So the Node.GetSupportReaction()
method have an overload which gets a LoadCombination
and returns the support reactions based on the load combination. LoadCombination
has a static
property named LoadCombination.DefaultLoadCombination
which has only one LoadCase
in it (the DefaultLoadCase
) with factor of 1.0. also everywhere that you should provide a LoadCombination
, if you do not provide any, then DefaultLoadCombination
will be considered by library. I’ve used DefaultLoadCase
and DefaultLoadCombination
in library to make working with library easier for people who are not familiar with load case and load combination stuff.
For getting the support reaction for the truss, we can simply call Node.GetSupportReaction()
to get support reaction for every node:
Force r1 = n1.GetSupportReaction();
Force r2 = n2.GetSupportReaction();
Force r3 = n3.GetSupportReaction();
Force r4 = n4.GetSupportReaction();
The plus operator is overloaded for Force struct
, so we can check the sum of support reactions:
var rt = r1 + r2 + r3 + r4;//shows the Fz=1000 and Fx=Fy=Mx=My=Mz=0.0
The forces (Fx, Fy and Fz) amount should be equal to sum of external loads and direction should be opposite to external loads to satisfy the structure static equilibrium equations.
All Codes Together
This is all codes above for truss example.
Please note that these codes are available in BriefFiniteElementNet.CodeProjectExamples project in library solution.
private static void Example1()
{
Console.WriteLine("Example 1: Simple 3D truss with four members");
// Initiating Model, Nodes and Members
var model = new Model();
var n1 = new Node(1, 1, 0);
n1.Label = "n1";//Set a unique label for node
var n2 = new Node(1, 1, 0) {Label = "n2"};//using object initializer for assigning Label
var n3 = new Node(1, 1, 0) {Label = "n3"};
var n4 = new Node(1, 1, 0) {Label = "n4"};
var n5 = new Node(0, 0, 1) {Label = "n5"};
var e1 = new TrussElement2Node(n1, n5) {Label = "e1"};
var e2 = new TrussElement2Node(n2, n5) {Label = "e2"};
var e3 = new TrussElement2Node(n3, n5) {Label = "e3"};
var e4 = new TrussElement2Node(n4, n5) {Label = "e4"};
//Note: labels for all members should be unique, else you will receive InvalidLabelException when adding it to model
e1.A = e2.A = e3.A = e4.A = 9e4;
e1.E = e2.E = e3.E = e4.E = 210e9;
model.Nodes.Add(n1, n2, n3, n4, n5);
model.Elements.Add(e1, e2, e3, e4);
//Applying restrains
n1.Constraints = n2.Constraints = n3.Constraints = n4.Constraints = Constraint.Fixed;
n5.Constraints = Constraint.RotationFixed;
//Applying load
var force = new Force(0, 1000, 1000, 0, 0, 0);
n5.Loads.Add(new NodalLoad(force));//adds a load with LoadCase of DefaultLoadCase to node loads
//Adds a NodalLoad with Default LoadCase
model.Solve();
var r1 = n1.GetSupportReaction();
var r2 = n2.GetSupportReaction();
var r3 = n3.GetSupportReaction();
var r4 = n4.GetSupportReaction();
var rt = r1 + r2 + r3 + r4;//shows the Fz=1000 and Fx=Fy=Mx=My=Mz=0.0
Console.WriteLine("Total reactions SUM :" + rt.ToString());
}
Example 2: Distributed Loads
In this example I’m going to analyze a single sloped steel frame which its elements are under uniform distributed loads as shown in image.
The purpose is to find the internal forces of members based on shown loads.
Creating members and nodes
First thing we are going to do is to create the members and assign the geometrical and mechanical properties to them
*Note: in next code, two columns are ‘e1’ and ‘e4’ and two beams are ‘e2’ and ‘e3’.
var model = new Model();
var n1 = new Node(10, 0, 0);
var n2 = new Node(10, 0, 6);
var n3 = new Node(0, 0, 8);
var n4 = new Node(10, 0, 6);
var n5 = new Node(10, 0, 0);
model.Nodes.Add(n1, n2, n3, n4, n5);
var secAA = SectionGenerator.GetISetion(0.24, 0.67, 0.01, 0.006);
var secBB = SectionGenerator.GetISetion(0.24, 0.52, 0.01, 0.006);
var e1 = new FrameElement2Node(n1, n2);
var e2 = new FrameElement2Node(n2, n3);
var e3 = new FrameElement2Node(n3, n4);
var e4 = new FrameElement2Node(n4, n5);
e1.Geometry = e4.Geometry = secAA;
e2.Geometry = e3.Geometry = secBB;
e1.E = e2.E = e3.E = e4.E = 210e9;
e1.UseOverridedProperties = e2.UseOverridedProperties = e3.UseOverridedProperties = e4.UseOverridedProperties = false;
In briefFiniteElement.NET, there is two ways to define section for a frame member. First is to define each of geometrical parameters (like area and second moments of area) by settings the FrameElement2Node.A or FrameElement2Node.Iy or FrameElement2Node.Iz, another way to do this is to define a polygon as section for frame element and library will automatically calculate the gemetrical parameters of section like area and second moments of area. By default library assumes that you are defining geometrical parameters by settings FrameElement2Node.A, Iy or Iz. if you want to provide a polygon as section, you should set the FrameElement2Node.Geometry property (which its type is PolygonYz) and right after that you should set the FrameElement2Node.UseOverridedProperties to false in order to say the library to use Geometry instead of A, Iy, Iz etc. There is also a useful static SectionGenerator class which provides some predefined sections like I sections or rectangular sections.
After creating elements and nodes, we can visualize the model in debug mode like image below.
(for more information see Add Debugger Visualizers for Visual Studio 2013, 2012 and 2010)
*Note: elements with UseOverridedProperties = true are shown with square sections (dimentsion of section automatically tunes for better visualization of elements) but elements with UseOverridedProperties = false will be shown whith their real section with real dimestion.
Next we have to set boundary conditions or supports. As this problem is 2D (in xz plane) we have to fix all movement in Y direction and also fix all rotations in x and z directions. As n1 and n5 are hinged supports we also have to fix the movement of n1 and n5 in both x and z directions too.
n1.Constraints = n2.Constraints = n3.Constraints = n4.Constraints = n5.Constraints = Constraint.FixedDy & Constraint.FixedRx & Constraint.FixedRz;//Dy fixed and Rx fixed and Rz fixed
n1.Constraints = n1.Constraints & Constraint.MovementFixed;
n5.Constraints = n5.Constraints & Constraint.MovementFixed;
Just have to note that the ‘&’ operator between two Constraint objects work like this: if both or one of Dx s are fixed then resut Dx is fixed, else Dx is released. It’s like Boolean operation while 'fixed' means '1' and 'free' means '0' !
Assigning Loads to Model
There are two distributed loads which are applying to model. Both are 1 Tonf/m which is equals to 10000 N/m but there is a difference as is illustrated before. The load which is applying to e2 element is in vertical direction but the load is applying to e3 element is perpendicular to e3 element. So it can be said the load on element e2 is a uniform load in GLOBAL z direction with amount of 10000 N/m and load which is applying to e3 is in element LOCAL z direction of coordination system and amount of 10000 N/m.
var ll = new UniformLoad1D(10000, LoadDirection.Z, CoordinationSystem.Global);
var lr = new UniformLoad1D(10000, LoadDirection.Z, CoordinationSystem.Local);
e2.Loads.Add(ll);
e3.Loads.Add(lr);
Checking for errors in model
In this example if we do not set e1.Geometry property, will get an error when equation system is solving because of structure is not stable. in this case finding the error maybe seems simple but its not always like this. for this reason there is a Trace property for Model class which will let user to see status of model while its working. For example on calling Model.Solve() for first time, several separated task are being done. first is calculating stiffness matrix of each element, next they should be assembeled into global stiffness matrix KG, then kff should extract from KG and finally equation system that contains kff should be solved. User do not see the information but information about these datas are possible to achieve from Model.Trace. It like System.Diagnostics.Trace functionality. The Trace member have TraceListeners property which is a list of ITraceListener objects. once an ITraceListener object is added to this list, it will receive trace info of model. There is a ready ITraceListener class named WpfTraceListener in BriefFiniteElementNet.Controls project that can be used for seeing the information. here is example:
var wnd = WpfTraceListener.CreateModelTrace(model);
new ModelWarningChecker().CheckModel(model);
wnd.ShowDialog();
the first line creates a WpfTraceListener object and add to Model.Trace.TraceListeners list. model.CheckForErrors () will check model for usual (obvious) errors and if any found, will be written into Model.Trace so we can see it our window. and finally wnd.ShowDialog() will shows the window. if for example we do not set Geometry property for e2 element, the will see this:
As you can see there is a table and a record inside of it. First field is TimeStamp which is showing the time, next is a brief message about the record, next is Error ID (in case that it be error or warning) next is Level. Level have 3 possible values:
* Error
* Warning
* Info
For example of Level, look at next pictures. Next is TargetField whici is showing the target of record, in our case it is e2 element and it have a link for more information (in this case it redirect to error list page) ad you can see description about error number MA10010 in the link.
If model do not contains error and we call wnd.ShowDialog() after call to model.Solve() like this:
var wnd = WpfTraceListener.CreateModelTrace(model);
new ModelWarningChecker().CheckModel(model);
wnd.ShowDialog();
model.Solve();
Will see something like this:
That contains information about nested tasks in solve process.
For more information see: Cheking model for errors
All together (Please note that these codes are available in BriefFiniteElementNet.CodeProjectExamples project in library solution.)
private static void Example2()
{
Console.WriteLine("Example 1: Simple 3D Frame with distributed loads");
var model = new Model();
var n1 = new Node(10, 0, 0);
var n2 = new Node(10, 0, 6);
var n3 = new Node(0, 0, 8);
var n4 = new Node(10, 0, 6);
var n5 = new Node(10, 0, 0);
model.Nodes.Add(n1, n2, n3, n4, n5);
var secAA = SectionGenerator.GetISetion(0.24, 0.67, 0.01, 0.006);
var secBB = SectionGenerator.GetISetion(0.24, 0.52, 0.01, 0.006);
var e1 = new FrameElement2Node(n1, n2);
e1.Label = "e1";
var e2 = new FrameElement2Node(n2, n3);
e2.Label = "e2";
var e3 = new FrameElement2Node(n3, n4);
e3.Label = "e3";
var e4 = new FrameElement2Node(n4, n5);
e4.Label = "e4";
e1.Geometry = e4.Geometry = secAA;
e2.Geometry = e3.Geometry = secBB;
e1.E = e2.E = e3.E = e4.E = 210e9;
e1.G = e2.G = e3.G = e4.G = 210e9/(2*(1 + 0.3));//G = E / (2*(1+no))
e1.UseOverridedProperties =
e2.UseOverridedProperties = e3.UseOverridedProperties = e4.UseOverridedProperties = false;
model.Elements.Add(e1, e2, e3, e4);
n1.Constraints =
n2.Constraints =
n3.Constraints =
n4.Constraints =
n5.Constraints =
Constraint.FixedDY & Constraint.FixedRX & Constraint.FixedRZ;//DY fixed and RX fixed and RZ fixed
n1.Constraints = n1.Constraints & Constraint.MovementFixed;
n5.Constraints = n5.Constraints & Constraint.MovementFixed;
var ll = new UniformLoad1D(10000, LoadDirection.Z, CoordinationSystem.Global);
var lr = new UniformLoad1D(10000, LoadDirection.Z, CoordinationSystem.Local);
e2.Loads.Add(ll);
e3.Loads.Add(lr);
var wnd = WpfTraceListener.CreateModelTrace(model);
new ModelWarningChecker().CheckModel(model);
wnd.ShowDialog();
model.Solve();
}
Analyzing the model and getting internal forces
For analyzing the model we should simply call the Model.Solve() Method. After model is solved we can get the internal force of members in any position of member's length with calling GetInternalForce method of FrameElement2Node class.
Following is the moment diagram based on output of FrameElement2Node.GetInternalForce() for this model.
Performance of BriefFiniteElement.NET
Look at Performance of BriefFiniteElement.NET
Points of Interest
At this time, the project have some features that are listed on Features of BriefFiniteElement.NET. There are also a set of features I planned to include in this library which I am working.
Support
If you have questions about the library or you've found a bug or have any idea/feature request, please let me know via DISCUSSION or ISSUEs on Codeplex, or mail me.
Contribute
If you have good knowledge in linear fem, look at HERE and help us solve our problems.
TODO
These examples will be added to this article:
 Example 1: Simple Truss (Done)
 Example 2: Element Distributed Load (Done)
 Example 3: Load Combination
 Example 4: Analyzing grid and checking result for statically equilibrium.
History
 9 July 2014  First version
 26 July 2014  Add Example #2 (Element Distributed Load)