Click here to Skip to main content
Click here to Skip to main content

Tagged as

Finite Element Method Programming in C#

, 3 Nov 2014 GPL3
Rate this:
Please Sign up or sign in to vote.
How to use BriefFiniteElement.NET to analyse solids and structures

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.

Examples

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:

  1. Create Model, Members and Nodes.
  2. Add the Nodes and Elements to Model
  3. Assign geometrical and mechanical properties to Elements
  4. Assign Constraints to Nodes (fix the DoF s)
  5. 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 TrussElement(n1, n5) {Label = "e1"};
var e2 = new TrussElement(n2, n5) {Label = "e2"};
var e3 = new TrussElement(n3, n5) {Label = "e3"};
var e4 = new TrussElement(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.AddRange(n1, n2, n3, n4, n5);
model.Elements.AddRange(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 = 9e-4;
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 Constraints in its static properties for example Constraint.Fixed or Constraint.Free. Here is more detailed information:

Property name

Description

Fixed

All 6 DoFs are fixed

Free

All 6 DoFs are released

MovementFixed

3 translation DoFs are fixed and 3 rotation DoFs are released

RotationFixed

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:

  1. First stage is creating stiffness matrix and factorizing stiffness matrix which will take majority of time for analyzing
  2. 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 LoadCases 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.

        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 TrussElement(n1, n5) {Label = "e1"};
            var e2 = new TrussElement(n2, n5) {Label = "e2"};
            var e3 = new TrussElement(n3, n5) {Label = "e3"};
            var e4 = new TrussElement(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 = 9e-4;
            e1.E = e2.E = e3.E = e4.E = 210e9;

            model.Nodes.AddRange(n1, n2, n3, n4, n5);
            model.Elements.AddRange(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, 0, -1000, 0, 0, 0);
            n5.Loads.Add(new NodalLoad(force));//adds a load whith 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          
        }

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 x-z 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);
model.CheckForErrors();
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);
model.CheckForErrors();
model.Solve();
wnd.ShowDialog();

Will see something like this:

That contains information about nested tasks in solve process.

For more information see: Cheking model for errors

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.

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)

License

This article, along with any associated source code and files, is licensed under The GNU General Public License (GPLv3)

Share

About the Author

Ehsan.MA
Software Developer
Iran (Islamic Republic Of) Iran (Islamic Republic Of)
I'm a civil engineer with several years of experiences in C# and WPF.

Comments and Discussions

 
QuestionMeaning of bitwise & operator PinmemberHiresdev3-Nov-14 19:28 
AnswerRe: Meaning of bitwise & operator PinmemberEhsan.MA3-Nov-14 21:38 
GeneralRe: Meaning of bitwise & operator PinmemberHiresdev5-Nov-14 1:47 
GeneralRe: Meaning of bitwise & operator PinmemberEhsan.MA5-Nov-14 1:51 
GeneralRe: Meaning of bitwise & operator PinmemberHiresdev5-Nov-14 2:46 
GeneralRe: Meaning of bitwise & operator PinmemberEhsan.MA6-Nov-14 7:03 
GeneralRe: Meaning of bitwise & operator PinmemberHiresdev6-Nov-14 18:58 
GeneralRe: Meaning of bitwise & operator PinmemberHiresdev7-Nov-14 2:48 
GeneralRe: Meaning of bitwise & operator PinmemberEhsan.MA7-Nov-14 5:40 
QuestionBoundary Conditions for Frame Example PinmemberHiresdev2-Nov-14 3:21 
AnswerRe: Boundary Conditions for Frame Example PinmemberEhsan.MA2-Nov-14 3:34 
QuestionI have a problem with constraints... PinmemberDeep Cloud21-Sep-14 10:33 
AnswerRe: I have a problem with constraints... PinmemberEhsan.MA21-Sep-14 10:54 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 11:06 
GeneralRe: I have a problem with constraints... PinmemberEhsan.MA21-Sep-14 11:10 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 11:33 
GeneralRe: I have a problem with constraints... PinmemberEhsan.MA21-Sep-14 12:08 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 11:46 
GeneralRe: I have a problem with constraints... PinmemberEhsan.MA21-Sep-14 11:53 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 12:04 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 12:08 
GeneralRe: I have a problem with constraints... PinmemberEhsan.MA21-Sep-14 12:16 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 12:37 
GeneralRe: I have a problem with constraints... PinmemberDeep Cloud21-Sep-14 12:15 
QuestionCan't find model.Nodes.AddRange PinmemberMember 1103117223-Aug-14 3:41 
AnswerRe: Can't find model.Nodes.AddRange PinmemberEhsan.MA23-Aug-14 7:50 
QuestionCongratulations. PinmemberMember 1055260618-Aug-14 2:55 
AnswerRe: Congratulations. PinmemberEhsan.MA19-Aug-14 9:43 
QuestionMy vote of 5 PinmemberloPetS15-Jul-14 5:41 
AnswerRe: My vote of 5 PinmemberEhsan.MA15-Jul-14 9:53 
GeneralMy vote of 4 PinprofessionalCatchExAs14-Jul-14 12:35 
GeneralRe: My vote of 4 PinmemberEhsan.MA15-Jul-14 0:57 
GeneralExamples.sln not found PinmemberMember 1094496114-Jul-14 4:29 
GeneralRe: Examples.sln not found PinmemberEhsan.MA15-Jul-14 0:52 
QuestionFinite Element Method PinmemberMember 990918810-Jul-14 10:09 
AnswerRe: Finite Element Method PinmemberEhsan.MA10-Jul-14 11:52 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.1411022.1 | Last Updated 3 Nov 2014
Article Copyright 2014 by Ehsan.MA
Everything else Copyright © CodeProject, 1999-2014
Layout: fixed | fluid