The complete program text using Matrix.NET follows the diagram. Color coding is "standard Visual Studio".
This example is #1 in the set of 3 solved by EMOI, a program published in 1977 by John Mays, Professor of Civil and Urban Engineering, University of Colorado at Denver. Note that the dofs are numbered irregularly, indicating they were numbered ‘manually’.
#region Using directives
using
System;
using System.Collections.Generic;
using System.Text;
using EngineeringObjects.LinearAlgebra;
#endregion
namespace
SampleFEProgram {
class Program
{
static
void Main(string[] args) {
// Elements
// 1: node1 to
node2
// 2: node2 to
node3
//
// Dofs: a = angular
rotation
// 1: x at node 2
// 2: y at node 2
// 3: a at node 2
// 4: a at node
3 // note anomoly; inconsistent dof numbering
// 5: x at node 3
// 6: y at node 3
// 7: x at node 1
// 8: y at node 1
// 9: a at node 1
//
// Load:
// 5 # force in -y
direction at node 2 (dof 2)
//
// Displacement:
// 0.01 in in x
direction at node 3 (dof 5)
// coordinates
Vector3D node1 =
new Vector3D(
0, 0, 0 );
Vector3D node2 =
new Vector3D(
12 * 12, 12 * 16, 0 );
Vector3D node3 =
new Vector3D(
12 * (12 + 15), 12 * 16, 0 );
node1.Print( "node1" );
node2.Print( "node2" );
node3.Print( "node3" );
Vector e1 = node2 - node1;
// Vector in direction of element 1
Vector e2 = node3 -
node2; // Vector in direction of element 2
e1.Print(
"e1" );
e2.Print( "e2" );
double
L1 = e1.Length;
double L2 = e2.Length;
// Beam properties
double I = 100;
// Moment of intertia, in^4
double A = 2;
// Area, in^2
double E = 30000;
// Young's modulus, #/in^2
Matrix K1 = FormElementK( A,
I, E, L1 );
K1.Print( "K1" );
// direction cosines for element
1
Vector u1 =
Vector.Unit( e1 );
u1.Print( "u1" );
// element rotation matrix
Lambda
Matrix Lambda =
Rotate( u1 );
Lambda.Print( "Lambda" );
Matrix KK1 =
Matrix.TransposeMultiply( Lambda, K1 * Lambda
);
//KK1.Print( "KK1" );
Matrix K2 = FormElementK( A,
I, E, L2 );
//K2.Print( "K2" );
Matrix G = new Matrix( 9, 9 ); // global stiffness
// accumulate element 1
int[] e1d = {
// element 1 dofs
7,8,9,1,2,3
};
IntVector e1Dofs =
new IntVector(
e1d );
Accumulate( G, KK1, e1Dofs );
G.Print( "G" );
// accumulate element 2
int[] e2d = {
1,2,3,5,6,4 // <<< irregular
dof numbering
};
IntVector e2Dofs =
new IntVector(
e2d );
Accumulate( G, K2, e2Dofs );
G.Print( "G" );
// save G for later Force
recovery
Matrix H =
new Matrix( G
);
// establish right hand side
Vector rhs =
new Vector(
G.Rows );
// Assume that displacements and loads are not specified for the same degree of freedom
// Specified displacement dofs
and values
Dictionary<int,
double> displacements = GetDisplacements();
// key: dof; value: displacement
ApplyDisplacements( G, displacements, rhs );
// Specified Loads
Dictionary<int,
double> loads = GetLoads();
// key: dof; value: load
ApplyLoads( rhs, loads );
//rhs.Print("rhs");
Vector x = Solve( G, rhs );
x.Print( "Displacements" );
// recover forces
Vector F = H * x;
F.Print( "External Forces"
);
UserPause();
}
//============================================================
private
static Dictionary<int,
double> GetLoads() {
Dictionary<int,
double> loads = new
Dictionary<int,
double>();
loads.Add( 2, -5.0 );
return loads;
}
//============================================================
private
static Dictionary<int,
double> GetDisplacements() {
Dictionary<int,
double> displacements =
new Dictionary<int,
double>();
displacements.Add( 7, 0 );
displacements.Add( 8, 0 );
displacements.Add( 9, 0 );
displacements.Add( 5, 0.01 );
displacements.Add( 6, 0 );
return displacements;
}
//============================================================
private
static Matrix
FormElementK(double A,
double I, double
E, double L) {
Matrix K =
new Matrix( 6,
6 );
double twelve = 12.0;
double six = 6.0;
double four = 4.0;
double two = 2.0;
double AEonL = A * E / L;
K[1, 1] = AEonL;
K[1, 4] = K[4, 1] = -AEonL;
K[4, 4] = AEonL;
double EIonL = E * I / L;
K[2, 2] = K[5, 5] = twelve * EIonL / (L * L);
K[5, 2] = -K[2, 2];
K[2, 5] = -K[2, 2];
K[3, 2] = K[2, 3] = K[6, 2] = K[2, 6] = six * EIonL / L;
K[6, 5] = K[5, 6] = -K[3, 2];
K[5, 3] = K[3, 5] = -K[3, 2];
K[3, 3] = K[6, 6] = four * EIonL;
K[6, 3] = K[3, 6] = two * EIonL;
return K;
}
//============================================================
private
static Matrix
Rotate(Vector u1) {
Matrix L =
new Matrix( 6,
6 );
double one = 1.0;
L[1, 1] = u1[1];
L[1, 2] = u1[2];
L[2, 1] = -u1[2];
L[2, 2] = u1[1];
L[3, 3] = one;
L[4, 4] =
u1[1];
L[4, 5] = u1[2];
L[5, 4] = -u1[2];
L[5, 5] = u1[1];
L[6, 6] = one;
return L;
}
//============================================================
private
static void
Accumulate(Matrix G,
Matrix K, IntVector edofs) {
int rg, cg;
for (int
r = 1; r <= K.Rows; ++r) {
rg = edofs[r];
for (int
c = 1; c <= K.Columns; ++c) {
cg = edofs[c];
G[rg, cg] += K[r, c];
}
}
}
//============================================================
private static
void ApplyDisplacements(Matrix
G, Dictionary<int,
double> displacements,
Vector rhs) {
foreach (KeyValuePair<int,
double> kvp in
displacements) {
ApplyADisplacement( G, rhs, kvp.Key, kvp.Value );
}
}
//============================================================
private
static void
ApplyADisplacement(Matrix G,
Vector rhs, int
eqn, double disp) {
if (disp != 0.0) {
// adjust rhs
for (int
i = 1; i <= G.Rows; ++i) {
rhs[i] -= disp * G[i, eqn];
}
rhs[eqn] = disp;
}
else {
// displacement is 0
rhs[eqn] = 0.0;
}
// adjust G
G.FillColumn( eqn, 0.0 );
G.FillRow( eqn, 0.0 );
G[eqn, eqn] = 1.0;
}
//============================================================
private
static void
ApplyLoads(Vector rhs,
Dictionary<int,
double> loads) {
foreach (KeyValuePair<int,
double> load in
loads) {
rhs[load.Key] = load.Value;
}
}
//============================================================
private
static Vector
Solve(Matrix G,
Vector rhs) {
double logdet, cond;
Matrix H =
Matrix.InvertSymmetric( G,
out logdet, out
cond );
Console.WriteLine();
Console.WriteLine(
"log(determinant): {0}", logdet );
Console.WriteLine(
"condition number: {0}", cond );
Vector res = H * rhs;
return res;
}
//==================================================
private
static void
UserPause() {
Console.WriteLine(
"\nPress 'Enter' to continue..." );
Console.ReadLine();
}
}
}
Output for this program:
log(determinant): 14.5632449149509
condition number:
150618.055555556
... Displacements ...
2.066610E-002 -4.587113E-002 -1.025382E-005 3.873863E-004
1.000000E-002
0.000000E+000 0.000000E+000 0.000000E+000 0.000000E+000
... External Forces ...
0.000000E+000 -5.000000E+000 0.000000E+000 3.552714E-015
-3.555368E+000
7.363707E-002 3.555368E+000 4.926363E+000 1.351102E+001