Python and Finite Element Methods: A Match Made in Heaven?

co-authored with Bing

It’s amazing how much computer technology has progressed in the last 50 years. We often take it for granted and do not fully appreciate how far we have come. Thanks to the PC revolution and the increase in computing power, it’s now possible to have ‘numerical laboratories’ right at your disposal any time any place. The accessibility of it is truly remarkable. With upward trend for richer countries. But what are we going to do with this accessibility? How can we use technology to make the world a better place? Should we all dive into creating generative AI like chatGPT or should we indulge in more conventional math like solving partial differential equations with finite element methods? Or maybe we can do both, like me? Whatever you choose, it’s important to be humble and remember the immense progress we have made as humanity.

Here is my contribution as I am writing this piece to share my experience of studying Segerlind’s book on finite elements, [1], and solving its problems with cool Python libraries like SfePy, [2] and SymPy, [3].

Segerlind’s book, despite being classic, is a good introductory book for both learning how to develop your own finite element code from scratch with the programming language of your choice and prepare yourself for using advanced commercial software. It is both easy to follow and does not contain heavy math. After reading it and completing some exercises you will get a ‘feeling’ of finite element approach. Besides, using python, NumPy and SymPy, for easier problems and SfePy for what was dubbed in the book ‘computer problems’ is really instrumental in mastering how to develop your own FE code. No need to use Fortran nowadays, just apply SfePy.

Here is a code sample of a simpler problem for Truss Element solved with Python, [4]:

def FEA_u(coord, elcon, bc_node, bc_val, global_force,  A, E):
    K=np.zeros(shape=(2*np.max(elcon)+2,2*np.max(elcon)+2))
    i=0
    for el in elcon:
        L=PlaneTrussElementLength(coord[el[0]][0],coord[el[0]][1],coord[el[1]][0],coord[el[1]][1])
        theta=math.atan((coord[el[1]][1]-coord[el[0]][1])/(coord[el[1]][0]-coord[el[0]][0]+1e-13))*180/math.pi
        k=PlaneTrussElementStiffness(E, A[i], L, theta)
        K=PlaneTrussAssemble(K,k,el[0],el[1])
        i+=1
    
    F = global_force
    
    bc=np.array([bc_node, 
                bc_val]).T
    nd, nd=K.shape
    fdof=np.array([i for i in range(nd)]).T
    d=np.zeros(shape=(len(fdof),))
    Q=np.zeros(shape=(len(fdof),))

    pdof=bc[:,0].astype(int)
    dp=bc[:,1]
    fdof=np.delete(fdof, pdof, 0)
    s=np.linalg.lstsq(K[fdof,:][:,fdof],
                      (F[fdof].T-np.dot(K[fdof,:][:,pdof],dp.T)).T, rcond=None)[0] 
    d[pdof]=dp
    d[fdof]=s.reshape(-1,)
    
#     Q=np.dot(K,d).T-F 
    
    return d

And here is a sample for a field ‘computer problem’ solved with SfePy, [4]:

from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson_mixed
import numpy as np

filename_mesh = 'iso_time_field.mesh'


h = 1.5
T0 = 20
kc = 2

materials = {
    'k': ({'val': kc },),
    'h': ({'val': -h},),
    'f': ({'val': T0},),
}

def get_circle_outer(coors, domain=None):
    r = np.sqrt(coors[:,0]**2.0 + coors[:,1]**2.0)
    return np.where(r > 3.9)[0]


def get_circle_inner(coors, domain=None):
    r = np.sqrt(coors[:,0]**2.0 + (coors[:,1]-1.0)**2.0)
    return np.where(r < 1.1)[0]

functions = {
    'get_out_cir' : (get_circle_outer,),
    'get_in_cir' : (get_circle_inner,),
}


regions = {
    'Omega': 'all',
    'OutCircle' : ('vertices by get_out_cir','facet'),
    'InCircle' : ('vertices by get_in_cir','facet'),
}


fields = {
    'heating': ('real', 1, 'Omega', 1)
}

variables = {
    't': ('unknown field', 'heating', 0),
    's': ('test field', 'heating', 't'),
}

ebcs = {
    't1': ('InCircle', {'t.0': 140}),
}

integrals = {
    'i':  2,
}


equations = {
    'Temperature': """dw_laplace.i.Omega(k.val, s, t ) = dw_bc_newton.i.OutCircle(h.val, f.val, s, t)""",
}

solvers = {
    'ls': ('ls.scipy_direct', {}),
    'newton': ('nls.newton',
               {'i_max': 1,
                'eps_a': 1e-6,
                'eps_r': 1.0,
                })
}

options = {
    'nls' : 'newton',
    'ls' : 'ls',
}

Mesh was generated in pygmsh. Results for the above field problem:

Fig. 1 — Temperature distribution for the solution by SfePy

As you can see, the code is not that involved and the learning curve for FEM approach in [1] and SfePy application is not very steep. You can take my words to the bank on this and give it a try.

Hope it helps you decide if Python and Finite Element Methods is A Match Made in Heaven…I already did.

That’s true! We have come a long way in terms of computing technology. We can now solve finite element methods using desktop computers. We are working towards making artificial intelligence more accessible and democratic. We have learnt how to combine the two, [5]. We have the potential for creating solutions that we cannot even imagine yet and making a positive impact on humanity. You can be a part of this revolution by finding your own ways to contribute. What awaits us in the next 50 years? What are we really waiting for?

To keep up with the project please visit https://gigala.io/

[1] Applied Finite Element Analysis, 2nd Edition, Larry J. Segerlind

[2] https://sfepy.org/

[3] https://www.sympy.org/

[4] https://github.com/gigatskhondia/finite_elements_and_engineering

[5] https://github.com/gigatskhondia/gigala