Implementation

Implementation#

Attention

This page shows a preview of the assignment. Please fork and clone the assignment to work on it locally from GitHub

After the workshop, the solution will be added to this preview and to the GitHub-repository

In this notebook you will implement the matrix method and check it with some sanity checks.

Our matrix method implementation is now completely stored in a local package, consisting of three classes. If you need a refresher on how to code with Classes and Objects, refer to the section on Object Oriented Programming in the MUDE-book, with additionally programming assignment 1.7

import matplotlib as plt
import numpy as np
sys.path.insert(1, '/matrixmethod_solution')
import matrixmethod_solution as mm
%config InlineBackend.figure_formats = ['svg']

import numpy as np
import matplotlib as plt
import matrixmethod as mm
%config InlineBackend.figure_formats = ['svg']

1. The Node class#

This class is stored in ./matrixmethod/node.py

The purpose of this class is to store node information and keep track of the total number of DOFs of the problem. Note the automatic bookkeeping we introduce in __init__. This simple but efficient way of keeping track of which DOFs belong to which nodes will make life much easier when we need to assemble matrices from multiple elements. The Node class doesn’t need any modification.

Exercise (1.1)

To test whether you understand how the class works, create two nodes on coordinates (\(0\),\(0\)) and (\(3\),\(4\)) and print the string representation of both nodes. The clear function is called to restart the node and DOF counters. Make sure this is done whenever you start solving a new problem.

mm.Node.clear()

node1 #= mm.Node(YOUR CODE HERE

print(node1)
#YOUR CODE HERE

Hint

Your output should look like this:

This node has:
 - x coordinate=0,
 - z coordinate=0,
 - degrees of freedom=[0, 1, 2],
 - load vector=[0. 0. 0.]
This node has:
 - x coordinate=3,
 - z coordinate=4,
 - degrees of freedom=[3, 4, 5],
 - load vector=[0. 0. 0.]

2. The Element class#

This class is stored in ./matrixmethod/elements.py

This class keeps track of each element in the model, including:

  • Cross-section properties

  • Element orientation (for coordinate system transformations)

  • Which Nodes make up each element, and in turn (with help of the Node class) which DOFs belong to each element

Apart from bookkeeping element data, the other main task of this class is to provide the element stiffness matrix in the global coordinate system (for subsequent assembly) and postprocess element-level fields. For now we keep postprocessing for next week and focus only on getting the correct stiffness matrix.

Here the class describes an element combining extension and Euler-Bernoulli bending. A similar (or inherited) class could also be implemented for different element types (e.g. shear beam, Timoshenko beam, cable elements, etc). Here we also keep it simple by assuming elements are all arranged in a 2D plane.

However, the implementation is incomplete:

  • The transformation matrix is missing in __init__, which is given in Transformations. Make sure you take into account that a positive \(\Delta z\) with a positive \(\Delta x\) gives a negative angle \(\alpha\). Make use of numpy.arctan2 to return the angle between \(-\pi\) and \(\pi\), numpy.arctan returns an angle between \(-\cfrac{\pi}{2}\) and \(-\cfrac{\pi}{2}\), and therefore cannot distinguish between all four quadrants.

  • The correct stiffness matrix for this extension-bending element coordinate system is missing in stiffness. You can derive the stiffness matrix yourself using pen and paper, SymPy or Maple, or copy the given stiffness matrix from Local stiffness matrix Euler-Bernoulli element.

  • We keep the functions which add a distributed load and compute the moments / displacements untouched for this week. Next week we’ll implement those as well.

Exercise (2.1)

Add the missing pieces to the code in ./matrixmethod/elements.py, before you perform the checks below. Do you specify your stiffness matrix in the global or local coordinate system?

Whenever you make changes to your code in the ./matrixmethod/ folder, you need to reimport those. Instead of restarting the kernel, we use some magic ipython commands. Run the cell below once. Consequently, whenever you save your changes in one of the .py-files, it’s automatically reloaded.

%load_ext autoreload
%autoreload 2

Error

Note that in this online book, you cannot make changes to multiple files simultaneously. These instructions are only applicable when you’re working on this assignment locally; please fork and clone the assignment to work on it locally from GitHub.

Exercise (2.2)

First, let’s check the stiffness matrix for a beam which doesn’t require rotation. Create a horizontal element with length \(2\) and \(EI=4\) and print both the transformation matrix and the stiffness matrix.

Do the matrices match with what you’d expect?

mm.Node.clear()
mm.Element.clear()

#YOUR CODE HERE

elem #= mm.Element(#YOUR CODE HERE

section = {}
section['EI'] #= YOUR CODE HERE
elem.set_section(section)

print(elem.T)
print(elem.stiffness())

Exercise (2.3)

Now, create a vertical element with length \(2\) and \(EI=4\) and print the transformation and stiffness matrix.

Do the matrices match with what you’d expect?

#YOUR CODE HERE

Exercise (2.4)

Now, create an element rotated in \(120 ^{\circ}\) with length \(2\) and print the transformation matrix.

Do the matrices match with what you’d expect?

#YOUR CODE HERE

Exercise (2.5)

Now, create an element rotated in \(60 ^{\circ}\) with length \(2\) and print the transformation matrix.

Do the matrices match with what you’d expect?

#YOUR CODE HERE

Exercise (2.6)

We’re going to apply the transformation matrix to the vector \([0,0,0,\sqrt{3},1,0]\) defined in global coordinates. What would be vector in the local coordinate system of the previous element?

Check your answer using pen and paper. Tip: make a drawing instead of doing all the algebra.

print(np.matmul(elem.T,np.array([0,0,0,np.sqrt(3),1,0])))

3. The Constrainer class#

This class is stored in ./matrixmethod/constrainer.py

This small class keeps track of which DOFs have prescribed displacements and takes care of applying these constraints to the global \(\mathbf{K}\) and \(\mathbf{f}\). For now we keep it simple and assume all constraints fix the DOF values to zero. Next week we will deal with non-zero prescribed values.

However, the implementation is incomplete:

  • The constrain function is incomplete, which should mimic the process of striking rows/columns of constrained DOFs and reduce the size of the system to be solved. Remember that Constrainer stores which DOFs are constrained in self.dofs, so all the others should be free. After gathering the free DOFs in an array, you will need to select the correct blocks of \(\mathbf{K}\) and \(\mathbf{f}\). For the stiffness matrix you will need the np.ix_() helper function (check its documentation here)

  • We keep the function which calculates supports reaction untouched for this week. Next week we’ll implement that one as well.

Exercise (3.1)

Add the missing pieces to the code, before you perform the check below

Exercise (3.2)

Take the inclined element of exercise 2.5 and a bending stiffness of \(1\). What happens if you invert \(\mathbf{K}\)? Now fix all degrees of freedom of the first node. What happens when you invert your ‘constrained’ \(\mathbf{K}\)? Are the dimensions of the ‘constrained’ \(\mathbf{K}\) correct?

#YOUR CODE HERE

con = mm.Constrainer()

#YOUR CODE HERE

f = np.zeros (6) #empty load vector
Kff, Fff = con.constrain( K, F )
print(np.shape(np.linalg.inv(Kff)))

4. Full implementation extension bar#

Having made our implementations, we now check them with two simple examples that serve as sanity checks. The first is a simple bar undergoing extension:

https://raw.githubusercontent.com/ibcmrocha/public/main/extpointload.png

With \(EA = 1000\), \(F = 100\) and \(L = 1\).

Use the code blocks below to set up and solve this problem using the classes above. The steps to follow are outlined below and short explanations/hints are given. Once you have a solution for the horizontal displacement of the node at the right end of the bar, compare it to the analytical solution you obtained in the first half of the course.

mm.Node.clear()
mm.Element.clear()

Exercise (4.1)

Create two nodes here. You can store them on a list or simply create them as two separate objects (e.g. node1 and node2).

#YOUR CODE HERE

Exercise (4.2)

Here we only have a single element, so there is no need to store it in a list yet. You are also going to need a dict defining the cross-section of the element.

#YOUR CODE HERE

Exercise (4.3)

Let’s define the boundary conditions. We create an instance of the Constrainer class to deal with prescribed displacements. Take a look at its functions and inform if Node 1 is fully fixed.

You also need to pass the load \(F\) on to Node 2. Check the member functions of Node to infer how that should be done.

#YOUR CODE HERE

Exercise (4.4)

Now assemble the global stiffness matrix and force vector. Since we only have one element, there is no real assembly to be performed other than getting the stiffness matrix of the single element and storing the load at Node 2 in the correct positions of \(\mathbf{f}\).

#YOUR CODE HERE

Exercise (4.5)

Constrain the problem and solve for nodal displacements.

#YOUR CODE HERE

Exercise (4.6)

Finally, compare the displacement at the end of the bar with the one coming from the ODE solution. Note that since our element is already suitable for frames combining extension and bending, \(\mathbf{u}\) has three entries. Which one is the entry that matters to us here? Did your solutions match? If so, that is a sign your implementation is correct. Can you use the function full_disp to obtain a vector of all displacements?

#EVENTUALLY YOUR CODE HERE

5. Full implementation bending beam#

In the first example above we tested our model under extension. But that does not really guarantee it will behave correctly in bending. That is the goal of this second sanity check. Let’s solve the following problem:

https://raw.githubusercontent.com/ibcmrocha/public/main/cantilever.png

Choose appropriate values yourself

When setting up and solving your model, note that we are now interested in \(w\) displacements, our load is now vertical and the cross-section property driving our deformation is now \(EI\). Good luck!

mm.Node.clear()
mm.Element.clear()

Exercise (5.1)

Create nodes

#YOUR CODE HERE

Exercise (5.2)

Create element

#YOUR CODE HERE

Exercise (5.3)

Set boundary conditions

#YOUR CODE HERE

Exercise (5.4)

Assemble the system of equations.

#YOUR CODE HERE

Exercise (5.5)

Constrain the problem and solve for nodal displacements

#YOUR CODE HERE

Exercise (5.6)

Check with the analytical solution

Did your solutions match? If so, your implementation is correct!

#EVENTUALLY YOUR CODE HERE