OOP for an Extensible FEM code

Aug 18 2016

7 minute read

fem , python

Introduction

After a solid 2 days of research, multiple articles and thesis, I figured that the whole OOP paradigm might be the best alternative to an extensible fem code.

The thing is: I need something simple enough so people can learn how the code works only with basic knowledge of the language and the theory.

The current framework

Top level description

User input as the Mesh file and the Problem parameters. Then, we call the solver that return an array with the unknown variable. Lastly we process the output by plotting.

agent Mesh
agent Problem
agent Solver
agent PostProcessor

Mesh -r-> Problem
Problem -r-> Solver
Solver -r-> PostProcessor

img/oopfem.png

Inside the solver is where the analysis encounter the model description.

Solver description

This module consists of assembling the global matrices and solve the linear system formed. The solver was created based on a procedural approach. Which means that the data is attached to a particular flow.

As an example on how the data is attached to the structure:

 1import numpy as np
 2
 3
 4class Model(object):
 5    """Builds a Model object"""
 6    def __init__(self):
 7        self.ne = 4             # number of elements
 8        self.nn = 9             # number of nodes
 9        self.ndof = 18         # number of degrees of freedom
10
11        self.XYZ = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [.5, 0], 
12                             [1, .5], [.5, 1], [0, .5], [.5, .5]])
13
14        self.CONN = np.array([[0, 4, 8, 7], [4, 1, 5, 8], [8, 5, 2, 6], 
15                              [7, 8, 6, 3]])
16
17        self.DOF = np.array([[0, 1, 8, 9, 16, 17, 14, 15],
18                             [8, 9, 2, 3, 10, 11, 16, 17],
19                             [16, 17, 10, 11, 4, 5, 12, 13],
20                             [14, 15, 16, 17, 12, 13, 6, 7]])
21        self.surf_of_ele = [0, 0, 0, 0]  # surface 0 for all elements
22
23def K_matrix(model, material):
24    """Build the GLOBAL stiffnees (K) matrix"""
25    K = np.zeros((model.ndof, model.ndof))
26
27    for e, conn in enumerate(model.CONN):
28        xyz = model.XYZ[conn]   # coordinates of nodes in element e
29        surf = model.surf_of_ele[e]
30        dof = model.DOF[e]
31
32        mat = material[surf]          # material from the surf tag
33
34        k = k_matrix(model, xyz, mat)  # This is where the ANALYSIS encounter the MODEL
35
36        id = np.ix_(dof, dof)
37
38        K[id] += k              # Assemble
39    return K
40
41def k_matrix(model, xyz, mat):
42    """Build the ELEMENT stiffness (k) matrix"""
43    k = np.zeros((8, 8))
44    return k
45
46model = Model()                 # instanciate the Model class
47material = {0: 10}              # material property for each surface key
48K = K_matrix(model, material)

We can see that if we try to implement a new element type, the k_matrix would have to change. For the assembler part it doesn't matter how the element stiffness is built. It only requires the element matrix and the degrees of freedom it affects. Therefore it makes sense to build a element object that has methods to build the stiffness matrix.

Objective

I need a architecture that allows extensibility of the code.
1. For example, I need a way to implement a new type of element, 8 nodes, without affecting the current solver.
2. This can be achieved by building the element matrices independently of the assembling (which is the current procedure).
3. In order to assemble the global matrix the only requirement is the element dof (degrees of freedom) and its element matrix.
4. Therefore, I need an array that stores the element type by its index, loop over it, call the specific constructor for this element type and store it at the global matrix.
5. The object Model seems to be the one responsible for caring the elements type.

New approach

Using the same example as before but now with a more robust framework.

object `Model`

First a class for the Model. This class will take a file as argument and extract data from it.

 1import numpy as np
 2
 3
 4class Model(object):
 5    """Builds a Model object"""
 6    def __init__(self):
 7        self.ne = 4             # number of elements
 8        self.nn = 9             # number of nodes
 9        self.ndof = 18         # number of degrees of freedom
10
11        self.XYZ = np.array([[0, 0],
12                             [1, 0],
13                             [1, 1],
14                             [0, 1],
15                             [.5, 0], 
16                             [1, .5],
17                             [.5, 1],
18                             [0, .5],
19                             [.5, .5]])
20
21        self.CONN = np.array([[0, 4, 8, 7], 
22                              [4, 1, 5, 8],
23                              [8, 5, 2, 6], 
24                              [7, 8, 6, 3]])
25
26        self.DOF = np.array([[0, 1, 8, 9, 16, 17, 14, 15],
27                             [8, 9, 2, 3, 10, 11, 16, 17],
28                             [16, 17, 10, 11, 4, 5, 12, 13],
29                             [14, 15, 16, 17, 12, 13, 6, 7]])
30
31        self.surf_of_ele = [0, 0, 0, 0]  # surface 0 for all elements
32
33        self.physical_surf = [0]         # physical surfaces tag
34
35        self.TYPE = [3, 3, 3, 3]
36
37model = Model()                 # Instanciate a Model object

This object contains all parameters that describe the physical entity analyzed.

object `Material`

This object contains material parameters. The class instance doesn't require any extra arguments to initialize. The materials are set in the most general way, the key is the material property and the value is a dictionary whose key is the surface and the value is the physical value of the property. The __dict__ attribute contains the instance attributes as dictionary and the update() method, an dictionary method, adds a new entry to the dictionary.

 1class Material(object):
 2    """Builds a Material object"""
 3    def __init__(self, **kw):
 4        self.__dict__.update(kw)
 5
 6material = Material(E={}, nu={0: 0.2})
 7
 8def E_time_dependent(t=1):
 9    return 1e5*t
10
11material.E[model.physical_surf[0]] = E_time_dependent
12print(material.E, material.nu, material.__dict__)

{0: <function E_time_dependent at 0x0000021DD38717B8>} {0: 0.2} {'nu': {0: 0.2}, 'E': {0: <function E_time_dependent at 0x0000021DD38717B8>}}

This implementation is good because:

I used to hard code the material properties, like material.E, as a attribute.
Any materials properties can be set. Therefore, if the problem is static, there is no need for inertial parameters, and if its a thermoelastic problem we need thermal properties.
Even material properties that vary through time.

Cautions:

Because the user is setting any material property we need to define a set of variable conventions that need to be checked when the parameter is going to be used. So, for example, when calculating the stiffness matrix we need to check if material.E and material.nu were defined, otherwise raise an error saying that those properties need to be assigned.

object `Element`

The element object contains the constructor of element matrices and vectors. It is created only when the solver is assembling the global matrices.

 1class Element(object):
 2    """Build an Element object with a constructor for a specific element type"""
 3
 4    def __init__(self, eid, model):
 5        Element.type = model.TYPE[eid]
 6        Element.conn = model.CONN[eid]
 7        Element.xyz = model.XYZ[self.conn]
 8        Element.dof = model.DOF[eid]  # dof must be know a priori
 9        Element.surf = model.surf_of_ele[eid]
10
11        if Element.type == 3:
12            self.constructor = Quad4()
13        elif Element.type == 4:
14            self.constructor = Quad9()
15        else:
16            raise Exception('Element not implemented yet!')
17
18
19class Quad4(object):
20    """Constructor of a 4-node quadrangle (TYPE 3) element"""
21
22    @staticmethod
23    def stiffness(material):
24        try:
25            E = material.E[Element.surf]
26        except:
27            raise Exception('Material not assigned for surface!')
28        
29        k = np.zeros((8, 8))
30        return k
31        
32
33class Quad9(object):
34    """Constructur of a 9-node quadrangle (TYPE 10) element"""
35    @staticmethod
36    def stiffness(material):
37        pass
38
39
40element_test_obj = Element(1, model)
41print(element_test_obj.type)
42print(element_test_obj.constructor.stiffness(material))

3
[[ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.]]

Attributes:

~~At first I don't see any use of attributes~~ Attributes derived from the model will be the same for all element types. This is done using class attributes instead of instance attribute (self)
Element.type seems useful. And it is, it's going to be used to decide which constructor object to instantiate.
The thing is: I don't want to repeat the model attributes for each element type, so I assigned them as class attributes. I'm feeling that there must be a better way, but I've already spent too much time.

Methods:

Static methods seems useful. Indeed, they are. I will use static methods to build the element matrices using the constructor object.

Observations

According to Archer (1996), there is a Element base class and the each element type inherit from it.

module `Solver`

The solver module contain the procedure to solve the problem. If it is a statics problem, then we call the statics.solver(). At first, this procedural approach for the solver seems reasonable since the methodology consists of a step-by-step framework. There is no need to create an class for the solver.

Inside the solver module there are the assembler functions (in different modules), for instance, the stiffness matrix assembler.

 1def K_matrix(model, material):
 2    """Build the GLOBAL stiffness matrix"""
 3    K = np.zeros((model.ndof, model.ndof))
 4    
 5    for eid in range(model.ne):
 6        element = Element(eid, model)  # Build the element based on type
 7        k = element.constructor.stiffness(material)  # constructor for each element type
 8        id = np.ix_(element.dof, element.dof)
 9        K[id] += k
10    return K

1def solver(model, material):
2    """Solves the problem"""
3    K = K_matrix(model, material)
4
5solver(model, material)

Introduction

The current framework

Top level description

Solver description

Objective

New approach

object Model

object Material

object Element

Observations

module Solver

See Also

object `Model`

object `Material`

object `Element`

module `Solver`