MATLAB
MATLAB functions can be called from FEAP if you have installed the MATLAB Engine (toolbox) when you set up MATLAB on your computer. The basic steps to getting this work are gone over below for a specific use case where I would like to use a MATLAB (*.m) file to compute a material model response for a user material.
The basics steps are as follows. (1) Have MATLAB and the MATLAB Engine installed. (2) Create a user macro that will start the MATLAB engine. This permits communication between the MATLAB data space and FEAP's data space using Unix sockets. (3) Set up your user material routine to send FEAP data to MATLAB to compute the material response and to receive is back from MATLAB. (4) Edit the makefile to allow for the compilation of subroutines with special MATLAB function calls as well as linking of the main program with the MATLAB Engine.
Setup the MATLAB Engine
The code to start up the MATLAB Engine will be placed in a user macro. To initialize the MATLAB Engine we just need to call this macro command once at the start of a FEAP run. For concreteness we use umacr6. It is important to save the file using a .F extender instead of .f, so that the compiler will run the preprocessor on it before compiling; this is all tested on a Mac which should be the same as Linux too. In the routine there are two user customizable strings: (a) The path to the matlab executable, '/Applications/MATLAB_2019a.app/bin/matlab' in my case, and the path to where my MATLAB *.m files are located, '/Users/sg/Feap/ver86/main' in my case. To use the macro, just type MLENgine at the macro prompt.
!$Id:$
subroutine umacr6(lct,ctl)
! * * F E A P * * A Finite Element Analysis Program
!.... Copyright (c) 1984-2019: Regents of the University of California
! All rights reserved
!-----[--.----+----.----+----.-----------------------------------------]
! Modification log Date (dd/mm/year)
! Original version 01/11/2006
! 1. Remove 'prt' from argument list 09/07/2009
! 2. Add 'uhelpfl' comment option 19/10/2017
! 3. Add 'help' comment option - replaces 'uhelpfl' 07/01/2019
! 4. Create the ML Engine init macro 10/16/2019
!-----[--.----+----.----+----.-----------------------------------------]
! Purpose: Initialize the MATLAB Engine
! Inputs:
! lct - Command character parameters
! ctl(3) - Command numerical parameters
! Outputs:
! N.B. Users are responsible for command actions. See
! programmers manual for example.
!-----[--.----+----.----+----.-----------------------------------------]
implicit none
# include "fintrf.h"
# include "mleng.h"
include 'iofile.h'
include 'umac1.h'
mwPointer engOpen
integer :: status, engEvalString
character (len=15) :: lct
logical :: pcomp
real (kind=8) :: ctl(3)
save
! Set command word
if(pcomp(uct,'mac6',4)) then ! Usual form
uct = 'mlen' ! Specify 'name'
elseif(pcomp(lct,'help',4)) then ! Write help information
write(*,*) 'COMMAND: MLENgine inits Matlab Engine'
elseif(urest.eq.1) then ! Read restart data
elseif(urest.eq.2) then ! Write restart data
else ! Perform user operation
ep = engOpen('/Applications/MATLAB_R2019a.app/bin/matlab')
if (ep.eq.0) then
write(*,*) 'WARNING: Matlab engine failed to initialize'
call plstop
endif
status = engEvalString(ep,'addpath /Users/sg/Feap/ver86/main')
if (status.ne.0) write(*,*) 'WARNING: Matlab path not set'
endif
end subroutine umacr6
Note that this macro also includes two *.h files
# include "fintrf.h" # include "mleng.h"
The first is a MATLAB file and is located in the MATLAB application directory. The second, I have placed in my FEAP main directory. It contains the pointer to the MATLAB Engine object that is returned by the endOpen() function call.
It contents are
mwPointer ep
common /mleng/ ep
User material
First I need a user material input routine. This is standard, except I am going to be lazy and not read in any material properties (they will be hard wired in the code). I will use umati6.f and give it the user name MLMAterial.
!$Id:$
subroutine umati6(vtype,vv, d, ud, n1,n3)
! * * F E A P * * A Finite Element Analysis Program
!.... Copyright (c) 1984-2019: Regents of the University of California
! All rights reserved
!-----[--.----+----.----+----.-----------------------------------------]
! Modification log Date (dd/mm/year)
! Original version 01/11/2006
!-----[--.----+----.----+----.-----------------------------------------]
! Purpose: User input for material model 6
! Inputs:
! vtype - Name of material model
! vv(5) - Command line real data
! d(*) - Program material parameter data
! Outputs:
! ud(*) - Material parameter data for model
! n1 - Number of history items/point (time dependent)
! n3 - Number of history items/point (time independent)
!-----[--.----+----.----+----.-----------------------------------------]
implicit none
character (len=15) :: vtype
logical :: pcomp
integer :: n1,n3
real (kind=8) :: vv(5),d(*),ud(*)
! Set command name
if(pcomp(vtype,'mat6',4)) then ! Default form DO NOT CHANGE
vtype = 'mlma' ! Specify 'name'
! Input user data and save in ud(*) array
else ! Perform input for user data
endif
end subroutine umati6
The stress-strain routine
The stress strain routine will be linear elastic where the material properties will be hard wired into the MATLAB routine that computes the response based on the strain. This is done in a MATLAB umatl6.m file which I place in my FEAP main directory (since that is where I told MATLAB to find such files). The files looks like:
function [sig,dd] = umatl6(eps) % Compute the stress and moduli % Linear isotropic elastic material with hardwired % Elastic properties % 6x6 symmetric identity Isym = eye(6,6); Isym(4,4) = 0.5; Isym(5,5) = 0.5; Isym(6,6) = 0.5; % 6x6 id outer-product id ooo = zeros(6,6); ooo(1:3,1:3) = 1.0; % Lame parameters mu = 10; lambd = 10; % Elasticity moduli 6x6 dd = 2*mu*Isym + lambd*ooo; % Stress 6x1 sig = dd*eps; end
It expects a 6x1 strain vector and returns a 6x1 stress vector and a 6x6 modulus matrix. My FEAP user material subroutine will call this MATLAB function with the strain and receive back the stress and modulus. Since I am using user material 6, I place my FEAP code in umatl6.F; note the .F extender so that the preprocessor does its job correctly.
!$Id:$
subroutine umatl6(eps,theta,td,d,ud,hn,h1,nh,ii,istrt, sig,dd,isw)
! subroutine umatl6( f ,detf ,td,d,ud,hn,h1,nh,ii,istrt, sig,dd,isw)
! * * F E A P * * A Finite Element Analysis Program
!.... Copyright (c) 1984-2019: Regents of the University of California
! All rights reserved
!-----[--.----+----.----+----.-----------------------------------------]
! Modification log Date (dd/mm/year)
! Original version 01/11/2006
!-----[--.----+----.----+----.-----------------------------------------]
! Purpose: User Constitutive Model 6
! Input:
! eps(*) - Current strains at point (small deformation)
! - Deformation gradient at point (finite deformation)
! theta - Trace of strain at point
! - Determinant of deformation gradient
! td - Temperature change
! d(*) - Program material parameters (ndd)
! ud(*) - User material parameters (nud)
! hn(nh) - History terms at point: t_n
! h1(nh) - History terms at point: t_n+1
! nh - Number of history terms
! ii - Current point number
! istrt - Start state: 0 = elastic; 1 = last solution
! isw - Solution option from element
! Output:
! sig(*) - Stresses at point.
! N.B. 1-d models use only sig(1)
! dd(6,*) - Current material tangent moduli
! N.B. 1-d models use only dd(1,1) and dd(2,1)
!-----[--.----+----.----+----.-----------------------------------------]
implicit none
# include "fintrf.h"
# include "mleng.h"
mwPointer mxCreateDoubleMatrix, engGetVariable
mwPointer ml_eps, ml_sig, ml_dd
mwSize N, M
#if MX_HAS_INTERLEAVED_COMPLEX
mwPointer mxGetDoubles
#else
mwPointer mxGetPr
#endif
integer :: status, engPutVariable, engEvalString
integer :: nh,istrt,isw, ii
real (kind=8) :: td
real (kind=8) :: eps(*),theta(*),d(*),ud(*),hn(nh),h1(nh)
real (kind=8) :: sig(*),dd(6,*)
if (isw.eq.3 .or. isw.eq.6 .or. isw.eq.4 .or. isw.eq.8) then
! Compute and output stress (sig) and (moduli)
! Initialize Matlab data space (should be done just once
! instead of each time with a destroy). Important use
! mwSize typed variables for all sizes
N = 6
M = 1
ml_eps = mxCreateDoubleMatrix(N,M,0)
! Copy strains to Matlab data space
#if MX_HAS_INTERLEAVED_COMPLEX
call mxCopyReal8ToPtr(eps, mxGetDoubles(ml_eps), N)
#else
call mxCopyReal8ToPtr(eps, mxGetPr(ml_eps), N)
#endif
! Place data in the MATLAB workspace (should pass along d, ud, hn, h1)
status = engPutVariable(ep, 'eps', ml_eps)
if (status.ne.0) then
write(*,*) 'WARNING: sending strain to Matlab failed'
endif
! Call the Matlab stress-strain function
status = engEvalString(ep,'[sig,dd] = umatl6(eps);')
if (status.ne.0) then
write(*,*) 'WARNING: mymat.m eval failed'
endif
! Get the stress and tangent from Matlab
ml_sig = engGetVariable(ep, 'sig')
ml_dd = engGetVariable(ep, 'dd')
#if MX_HAS_INTERLEAVED_COMPLEX
call mxCopyPtrToReal8(mxGetDoubles(ml_sig), sig, N)
call mxCopyPtrToReal8(mxGetDoubles(ml_dd ), dd , N*N)
#else
call mxCopyPtrToReal8(mxGetPr(ml_sig), sig, N)
call mxCopyPtrToReal8(mxGetPr(ml_dd ), dd , N*N)
#endif
call mxDestroyArray(ml_eps)
call mxDestroyArray(ml_sig)
call mxDestroyArray(ml_dd )
endif
end subroutine umatl6
There is a lot to be explained here, best is to review the MATLAB Engine documentation, but here are a few main points. (1) For every variable that you will send to MATLAB you need to create a pointer to its MATLAB version. This is done with the function mxCreateDoubleMatrix() for double precision arrays. It takes 3 arguments, number of rows, number of columns, and flag for if it is complex or not. This command creates the pointer. The MATLAB pointer is hooked to your FEAP Fortran array using mxCopyReal8ToPtr() (the syntax of which depends on your MATLAB version, which the preprocessor figures out). (2) To send data from Fortran arrays to MATLAB arrays one uses the engPutVariable() function. (3) To execute MATLAB commands in the MATLAB data space one uses the engEvalString() function. In our case we evaluate the MATLAB function umatl6() which we wrote. The return values end up in sig and dd in the MATLAB data space. (3) To retrieve the values, we first set up a pointers to these MATLAB arrays using the engGetVariable() function. We then copy the data from these pointers to our FEAP Fortran arrays using the mxCopyPtrToReal8() function. (4) Lastly, I destroy all the pointers with mxDestroyArray().
Note there is a bit of waste here. I should create the pointers only once and then destroy them only once where I terminate FEAP but that requires being more sophisticated and makes it harder to understand the basic steps.
Setting the makefiles
To get this all to compile, you need to edit the makefile in FEAP's main folder. It needs to be able to find the MATLAB system include files and the location of the MATLAB libraries. What you see below, is what I needed on my Mac.
include ../makefile.in
MLROOT = /Applications/MATLAB_R2019a.app
OBJECTS = feap86.o umacr6.o umatl6.o umati6.o
#OBJECTS = feap86.o
feap: $(OBJECTS) $(ARFEAP)
ranlib $(ARFEAP)
$(FF) -Wl,-no_pie -o feap $(OBJECTS) \
$(ARFEAP) -L/$(MLROOT)/bin/maci64 \
$(LDOPTIONS) -leng -lmx
dsymutil feap
@@echo "--> FEAP executable made <--"
# $(ARCHIVELIB) \
# $(ARPACKLIB) \
# $(ARFEAP) $(LDOPTIONS)
# @@echo "--> FEAP executable made <--"
# UBUNTU and other GCC loader type machines
# Replace $(ARCHIVELIB) \ by:
# -L$(FEAPHOME8_5)/packages/arpack/archive \
# -Wl,-whole-archive -larchive -Wl,-no-whole-archive \
# Replace $(ARPACKLIB) \ by:
# -L$(FEAPHOME8_5)/packages/arpack \
# -Wl,-whole-archive -larpack -Wl,-no-whole-archive \
clean:
rm -f *.o
rm -f *genmod.mod
rm -f *genmod.f90
@@echo "--> Files cleaned <--"
fclean:
rm -f feap
rm -r -f feap*.dSYM
@@echo "--> Feap cleaned <--"
%.o: %.F
$(FF) -c $(FFOPTFLAG) -I$(FINCLUDE) -I$(MLROOT)/extern/include $< -o $@
%.o: %.f
$(FF) -c $(FFOPTFLAG) -I$(FINCLUDE) $< -o $@
%.o: %.f90
$(FF) -c $(FFOPTFLAG) -I$(FINCLUDE) $< -o $@
%.o: %.c
$(CC) -c $(CCOPTFLAG) -I$(CINCLUDE) $< -o $@