Continuing the recent series on using Fortran code published in Programming the Finite Element Method (5th ed. John Wiley & Sons, I.M. Smith, D.V. Griffiths and L. Margetts (2014)), this post looks at linking Excel to Fortran code compiled as a dll (dynamic link library).

It is quite possible to link Excel VBA directly to Fortran based dlls, but because Python provides a convenient link to numerous scientific and maths routines in Numpy, Scipy, and other libraries, I have chosen to create the link through Python, using xlwings.  The spreadsheet, dll files, and Python and Fortran code can be downloaded from:

PlateMC-dll.zip

To run the code requires Python and xlwings to be installed.  The Anaconda Python package includes xlwings by default, and also allows The MinGW package to be installed, which includes Gfortran.  GFortran is not required to run the application, but it is needed if you want to make any changes to the Fortran code and recompile.  The dll file may be copied to the same directory as the spreadsheet and Python file (Py64f.py), or to any directory on the system search path.

As for compiling an exe file and linking to Excel, the process is not difficult, but it is not well documented, and any deviation from the exact procedure is likely to result in errors accompanied with a collection of exceedingly unhelpful error messages.

In the previous version of this spreadsheet the data for the compiled program was written to the hard disk as text files, which were opened and read by the Fortran program, which then wrote back the results to the hard disk, to be opened and read by Excel.  This is a slow and cumbersome procedure, but it does remove all problems of differing data formats.  In contrast, when the data is transferred directly from Excel to the Fortran based dll it is essential all data and objects are in a format that is recognisable to both programs.  The procedure I have used to ensure data compatibility was:

• In Excel VBA collect the input data as a series of named worksheet ranges, converted to variant arrays.
• Transfer to Python with an xlwings UDF (user defined function).  Arrays are transferred as Python tuples.
• Convert the tuples to Numpy arrays of the appropriate data type, and if necessary convert values to the correct data type.  For this application all values were either double64 or int (equivalent to VBA double and long).
• Convert the arrays and values to c_type format, which is compatible with both C and Fortran.
• Call the required function from the dll, passing all arguments by reference.
• Return the result arrays and values to VBA using the original Python variables.
• Paste the results to named ranges in the spreadsheet.
• Update the deformed mesh plot in VBA.

The VBA code is straightforward.  The Sub P64py reads data from the ranges “Plate_def”, “Tolerance”, “srf”, and “props”, adjusting the range sizes where necessary, then calls the function xl_P64, which converts empty cells in the Plate_Def range to zeros, then calls the Python function py64 in the module Py64f.py.

The main features of the Python code are:

• Import numpy, os and xlwings.
• Import the required data types and methods from the ctypes library (byref, cdll, c_int, c_double), and import numpy.ctypeslib to handle Numpy arrays.
• Load the dll file, using cdll.LoadLibrary (see Python code below).  The code is set up to look for the dll first in the workbook directory, or if that fails look anywhere on the search path.
• The data from the “Plate_def” worksheet range is passed from VBA as a variant array with all numerical values passed as doubles.  In both Numpy and Fortran arrays must have a single data type, and the Fortran code requires some of the input data as doubles and some as integers.  For this reason the dataa array must be split into two separate arrays, containing doubles (meshdims) and integers (meshnums).
• Arrays for the output data must be sized and created in Python, and passed as ctype pointers to Fortran, together with the size of any variable dimensions.
• Having converted and created the required arrays in Numpy format, ctype pointers are created with statements such as:
  coordp = npc.as_ctypes(coorda)
  Note that it is important that the pointer has a different name to the Numpy array.
• Values are converted to ctypes using c_int or c_double, e.g:
  nsrf = c_int(nsrfi)
• The p64 function is called from p64dll, with all arguments passed “byref” (see code below)
• The data passed to the Fortran code consists of pointers to the Python variables, hence the arrays and values returned from p64 may be accessed using the original Python variables names.

The full Python code is shown below:

import numpy  as np
from ctypes import byref, cdll, c_int, c_double, c_char_p
import numpy.ctypeslib as npc
import os
import time
import xlwings as xw

#Load p64s.dll
try:
    wkbpath = os.path.dirname(os.path.abspath(__file__))
    p64dll = cdll.LoadLibrary(wkbpath + '/p64s.dll')
except:
    p64dll = cdll.LoadLibrary('./p64s.dll')
    
def py64(dataa, tola, srfa, propa):
# Convert dataa to Numpy array, and split into double and integer arrays
    dataa = np.array(dataa)
    meshdims = dataa[0:2,:]
    meshnums = dataa[2:4,0:2].astype(int)

#Calculate size of output arrays    
    np_types = int(dataa[4,0])
    nx1 = meshnums[0,0]
    nx2 = meshnums[0,1]
    ny1 = meshnums[1,0]
    ny2 = meshnums[1,1]
    nye=ny1+ny2
    nels=nx1*nye+ny2*nx2
    nn=(3*nye+2)*nx1+2*nye+1+(3*ny2+2)*nx2
#Create Numpy arrays and convert to ctypes pointers
#Convert integers to c_int
    coorda = np.zeros((nn,2))
    coordp = npc.as_ctypes(coorda)
    g_num = np.zeros((nels,8), dtype=np.int)
    g_nump = npc.as_ctypes(g_num)
    srfa = np.array(srfa)[:,0]
    nsrfi = srfa.shape[0]
    nsrf = c_int(nsrfi)
    nodedisa = np.zeros((nn,nsrfi*3))
    nodedisp = npc.as_ctypes(nodedisa)
    nnp = c_int(nn)
    meshdimp = npc.as_ctypes(meshdims)
    meshnump = npc.as_ctypes(meshnums)
    tolp = npc.as_ctypes(np.array(tola)[0])
    srfp = npc.as_ctypes(srfa)
    propa = np.array(propa)
    propp = npc.as_ctypes(propa)
    resa = np.zeros((5,nsrfi))
    resp = npc.as_ctypes(resa)
    np_types = c_int(np_types)
    nelsp = c_int(nels)
    
 #Call p64 with all ctype arguments passed by reference   
    iErr = p64dll.p64(byref(resp), byref(meshdimp), byref(meshnump), byref(tolp), byref(srfp), byref(nsrf), 
    byref(propp), byref(np_types), byref(coordp), byref(nnp), byref(nodedisp), byref(g_nump), byref(nelsp))
 
 #Return results using original Python variables
    return [resa, coorda, nodedisa, g_num, nn,  nels] 

The Fortran code requires an interface to read input data and return results. Note that dynamic arrays must have their size defined by integer arguments passed from the calling program. Also note that Fortran arrays are defined in the opposite direction to Python and VBA arrays (i.e. columns, rows):

module main
implicit none
integer, parameter:: dp=SELECTED_REAL_KIND(15)

contains

...  Functions called by p64

END module main

integer function p64(resa, meshdims, meshnums, tola, srfa, nsrf, prop, &
   np_types, g_coord, nn, NodeDisp,g_num, nels)
!-------------------------------------------------------------------------
! Program 6.4 Plane strain slope stability analysis of an elastic-plastic
!             (Mohr-Coulomb) material using 8-node rectangular
!             quadrilaterals. Viscoplastic strain method.
!-------------------------------------------------------------------------
  use main
  implicit none 
 integer, intent(in):: meshnums(2,2), nsrf, np_types, nn, nels
 integer, intent(out):: g_num(8,nels)
 real(dp),INTENT(IN)::meshdims(3,2), tola(3), srfa(nsrf), prop(7,np_types)
 real(dp),INTENT(OUT)::resa(nsrf,5), g_coord(2,nn), NodeDisp(nsrf*3,nn)

...
  

The Fortran code is compiled to a dll as shown below:

• gfortran -shared -fno-underscoring -o p64.dll P64dll.f03

As for the compilation of an exe file, if the program is to run on a computer without MinGW it must be compiled with the -static option:

• gfortran -shared -static -fno-underscoring -o p64s.dll P64dll.f03

The resulting output is shown below:

The VBA plotting routine has been re-written to allow changes in the number of elements to be plotted.

See the next post for further examples and discussion of results with varying element size.