Python material
Jump to navigation
Jump to search
It is possible to use Python to compute stress and modulus for FEAP. This can be helpful during early stages of material model development. Here is a very simple example for a linear elastic material model.
The idea is to create a user constitutive model (umatiXX.f and umatlXX.f). To use the material, one first opens up a SOCKET in Python (this is the server). Then FEAP, once started creates a client that sends strain and material parameters over the socket. Python computes the stress and modulus and returns it over the socket. This is a rudimentary example but it shows all the steps involved. More complex models will require the packing and sending of additional data over the socket but the ideas are the same.
Note: This has only been tested on a Mac. I imagine, it will work just fine on Linux. For Windows, you may need to include special libraries during the link, for example, -lws2_32 according to what I have found on the web. The Python version used during testing was version 3.9.
Python script
The Python script needs to be started first. It starts the socket server. Do so with python umatl3.py assuming you have called the file umatl3.py
# Simple example of a linear elastic material in Python for FEAP
# Start this file "python umatl3.py" before starting FEAP
import socket
import numpy as np
import struct
import sys
# User material model data
umat_num = 3
size_d = 2 # number of material data parameters
# Server configuration
HOST = '127.0.0.1' # Standard loopback interface address (localhost)
PORT = 65400 + umat_num # Port to listen on (non-privileged ports are > 1023) adjust by umat num
# Linear elastic model
def stress(eps,d):
# Extract material properties
mu = float(d[0]) # Ordering set by umati3.f
lam = float(d[1]) # Explicit cast to avoid Python deprecation warnings
# Build modulus tensor
dd = np.array([[2*mu+lam, lam, lam, 0., 0., 0.],
[lam, 2*mu+lam, lam, 0., 0., 0.],
[lam, lam, 2*mu+lam, 0., 0., 0.],
[0., 0., 0., mu, 0., 0.],
[0., 0., 0., 0., mu, 0.],
[0., 0., 0., 0., 0., mu]],dtype=np.float64)
sig = dd@eps # Compute stress
return sig, dd # Return stress and modulus
def recvall(sock, n):
"""
Helper function to receive exactly n bytes from the socket.
"""
data = bytearray()
while len(data) < n:
packet = sock.recv(n - len(data))
if not packet:
# Connection closed before all data was received
return None
data.extend(packet)
return data
def run_server():
"""
Sets up a TCP socket server, listens for FEAP to send in strain and properties
and sends stress and modulus back
"""
print(f"Python Server: Starting on {HOST}:{PORT}")
# Create a socket object using IPv4 and TCP
with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
# Bind the socket to the host and port
try:
s.bind((HOST, PORT))
except OSError as e:
print(f"Python Server Error: Could not bind to port {PORT}. Is it already in use?")
print(f"Details: {e}")
sys.exit(1)
# Listen for incoming connections (allow up to 1 pending connection)
s.listen(1)
print("Python Server: Listening for incoming connections...")
# Accept a connection when one arrives
conn, addr = s.accept()
with conn:
print(f"Python Server: Connected by {addr}")
# --- Receive the 6x1 strain vector from FEAP + material properties ---
# Expecting 6 double-precision floats (8 bytes each) = 48 bytes
# Plus size_d * 8 bytes
in_size = 6+size_d # input size in doubles
vector_size_bytes = in_size * struct.calcsize('d') # 'd' for double precision
# Loop to process requests
while(1) :
data = recvall(conn,vector_size_bytes)
if data is None:
print("Python Server: Connection closed.")
return
try:
# Unpack the bytes into a tuple of floats
# e.g. '6d' means 6 double-precision floats
received_floats = struct.unpack('{}d'.format(in_size), data)
data_in = np.array(received_floats, dtype=np.float64).reshape(in_size, 1)
# separate strain from material properties
eps = data_in[0:6]
d = data_in[6:6+size_d]
except struct.error as e:
print(f"Python Server Error: Failed to unpack received data. {e}")
return
# Compute stress and modulus
sig, dd = stress(eps,d)
# Pack the stress and modulus for sending
# 'd' for double-precision float
stress_bytes = struct.pack('6d', *sig.flatten())
modulus_bytes = struct.pack('36d', *dd.flatten())
# Concatenate the bytes
response_data = stress_bytes + modulus_bytes
# Send the combined results back to the Fortran client
conn.sendall(response_data)
if __name__ == "__main__":
run_server()
pysock.f90
The material model makes use of a module that contains a number of interface definitions. It importantly stores sockfd which has the handle for the socket.
! * * F E A P * * A Finite Element Analysis Program
!.... Copyright (c) 1984-2025: Regents of the University of California
! All rights reserved
!-----[--.----+----.----+----.-----------------------------------------]
! Modification log Date (dd/mm/year)
! Original version 07/12/2025
!-----[--.----+----.----+----.-----------------------------------------]
! Purpose: Module for socket connections
!-----[--.----+----.----+----.-----------------------------------------]
module pysock
use ISO_C_BINDING, ONLY : C_INT, C_SHORT, C_CHAR, C_DOUBLE, C_PTR, C_LOC, C_SIZE_T, C_NULL_CHAR !
use ISO_C_BINDING, ONLY : C_NULL_PTR
implicit none
integer(C_INT) :: sockfd ! Socket file descriptor set in umatiXX.f
! --- Constants for Socket Programming ---
integer(C_INT), parameter :: AF_INET = 2 ! IPv4 Internet protocols
integer(C_INT), parameter :: SOCK_STREAM = 1 ! TCP socket type
integer(C_INT), parameter :: IPPROTO_TCP = 6 ! TCP protocol (usually 0 is fine for SOCK_STREAM)
! --- Socket Address Structure (sockaddr_in) ---
! This derived type mirrors the C 'struct sockaddr_in' for IPv4 addresses.
! It's used to specify the server's address and port.
type, BIND(C) :: sockaddr_in_type
integer(C_SHORT) :: sin_family ! Address family (AF_INET)
integer(C_SHORT) :: sin_port ! Port number (in network byte order).
integer(C_INT) :: sin_addr ! IP address (in network byte order)
! The sin_zero padding is often omitted in Fortran, but included for strict
! C compatibility to ensure correct size and alignment.
character(kind=C_CHAR) :: sin_zero(8)
end type sockaddr_in_type
! --- C Function Interfaces (POSIX Sockets and Network Utilities) ---
! These interfaces declare how Fortran can call C socket functions.
interface
! socket(domain, type, protocol)
integer(C_INT) function socket_c(domain, ttype, protocol) BIND(C, NAME='socket')
import C_INT
integer(C_INT), value :: domain
integer(C_INT), value :: ttype
integer(C_INT), value :: protocol
end function socket_c
! connect(sockfd, addr, addrlen)
integer(C_INT) function connect_c(sockfd, addr, addrlen) BIND(C, NAME='connect')
import C_INT, C_PTR
integer(C_INT), value :: sockfd
type(C_PTR), value :: addr
integer(C_INT), value :: addrlen
end function connect_c
! send(sockfd, buf, len, flags)
integer(C_INT) function send_c(sockfd, buf, lenn, flags) BIND(C, NAME='send')
import C_INT, C_PTR, C_SIZE_T
integer(C_INT), value :: sockfd
type(C_PTR), value :: buf
integer(C_SIZE_T), value :: lenn
integer(C_INT), value :: flags
end function send_c
! recv(sockfd, buf, len, flags)
integer(C_INT) function recv_c(sockfd, buf, len, flags) BIND(C, NAME='recv')
import C_INT, C_PTR, C_SIZE_T
integer(C_INT), value :: sockfd
type(C_PTR), value :: buf
integer(C_SIZE_T), value :: len
integer(C_INT), value :: flags
end function recv_c
! close(sockfd)
integer(C_INT) function close_c(sockfd) BIND(C, NAME='close')
import C_INT
integer(C_INT), value :: sockfd
end function close_c
! htons (host to network short) - converts 16-bit integer from host to network byte order
integer(C_SHORT) function htons_c(hostint) BIND(C, NAME='htons')
import C_INT, C_SHORT
integer(C_INT), value :: hostint
end function htons_c
! inet_addr - converts IPv4 dotted-decimal string to network byte order integer
! Returns in_addr_t (typically unsigned int, so C_INT in Fortran)
integer(C_INT) function inet_addr_c(cp) BIND(C, NAME='inet_addr')
import C_CHAR, C_PTR, C_INT
character(kind=C_CHAR), dimension(*), intent(in) :: cp
end function inet_addr_c
end interface
end module pysock
umati3.f
This routine is a typical input routine for a user constitutive model in FEAP. However, it also sets up the client side of the socket for communicating with Python. Importantly, it initializes sockfd.
!$Id:$
subroutine umati3(vtype,vv, d, ud, n1,n3)
! * * F E A P * * A Finite Element Analysis Program
!.... Copyright (c) 1984-2025: 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 3
! 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)
!-----[--.----+----.----+----.-----------------------------------------]
use pysock ! sockfd (among others)
implicit none
include 'iofile.h'
! Input variables
character (len=15) :: vtype
logical :: pcomp, tinput, errck
character :: text(2)*15,termination*15
integer :: n1,n3
real (kind=8) :: vv(5),d(*),ud(*),td(2)
! Socket variables
! Define server details
integer(C_INT), parameter :: umat_num = 3 ! Set to user material number for unique port
integer(C_INT), parameter :: SERVER_PORT = 65400+umat_num ! Port number.
integer(C_INT) :: status_c ! For checking return codes of C functions
type(sockaddr_in_type), target :: server_addr ! Structure to hold server's address.
! Character buffer to hold the IP address string for inet_addr_c
character(kind=C_CHAR, len=16) :: ip_address_str ! e.g., "127.0.0.1" + null terminator
! Character buffer to zero out the entire sockaddr_in_type structure
character(kind=C_CHAR, len=sizeof(server_addr)), target :: zerobuf
! Set command name
if(pcomp(vtype,'mat3',4)) then ! Default form DO NOT CHANGE
vtype = 'pymat' ! Specify 'name'
! Input user data and save in ud(*) array
else ! Perform input for user data
text(1) = 'start'
termination = ' ' ! or for example 'my input end'
do while (.not.pcomp(text(1),termination,4))
errck = tinput(text,2,td,2) ! Read two keywords and up to 2 values
if (pcomp(text(1),'modu',4)) then
if (pcomp(text(2),'mu',2)) then
d(1) = td(1) ! Shear modulus (mu)
elseif (pcomp(text(2),'lam',3)) then
d(2) = td(1) ! Lame modulus (lambda)
endif
endif
end do
write(iow,1000) d(1),d(2) ! Output data
! Create socket
sockfd = socket_c(AF_INET, SOCK_STREAM, IPPROTO_TCP)
if (sockfd == -1) then
write(iow,*) 'Error: Failed to create socket.'
call plstop(.true.)
end if
! Zero out the entire structure first to ensure padding (sin_zero) is cleared
zerobuf = repeat(char(0, kind=C_CHAR), len(zerobuf))
server_addr = transfer(zerobuf, server_addr)
server_addr%sin_family = AF_INET
! Convert IP address string to network byte order using inet_addr_c
ip_address_str = "127.0.0.1" // C_NULL_CHAR
server_addr%sin_addr = inet_addr_c(ip_address_str)
! Convert port to network byte order using htons_c
server_addr%sin_port = htons_c(SERVER_PORT)
! Connect to the server
status_c = connect_c(sockfd, C_LOC(server_addr),
& int(sizeof(server_addr), kind=C_INT))
if (status_c == -1) then
write(iow,*) 'Failed to connect to server.'
status_c = close_c(sockfd)
call plstop(.true.)
end if
write(iow,*) 'Connected to server.'
endif
1000 format(/'Python Material Model',' : Shear modulus = ',1p,1e12.5:/
& ' Lame modulus = ',1p,1e12.5/)
end subroutine umati3
umatl3.f
This is the actual material model. It is very simple and only sends in the current strain to Python. If one had a more complex model you would have to pack more data into the send buffer and adjust umatl3.py to handle the receive. A noted defect of the programming here is that the client side socket is not elegantly closed when FEAP exits. Other bullet proofing is probably also in order.
!$Id:$
subroutine umatl3(eps,theta,td,d,ud,hn,h1,nh,ii,istrt, sig,dd,isw)
! subroutine umatl3( 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-2025: Regents of the University of California
! All rights reserved
!-----[--.----+----.----+----.-----------------------------------------]
! Modification log Date (dd/mm/year)
! Original version 01/11/2006
! 1. Dimention td(2) 12/08/2021
!-----[--.----+----.----+----.-----------------------------------------]
! Purpose: User Constitutive Model 3
! Linear elasticity evaluated using Python over a socket
! Input:
! eps(9,3)- Current strains at point (small deformation)
! f(3,3,4)- Deformation gradient at point (finite deformation)
! 1 = F_n+1; 2 = F_n; 3 = G_n+1; 4 = G_n
! theta - Trace of strain at point (small deformation)
! detf(4) - Determinant of deformation gradient (finite defm)
! 1 = detF_n+1; 2 = detF_n; 3 = detG_n+1; 4 = detG_n
! td(2) - 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)
!-----[--.----+----.----+----.-----------------------------------------]
use pysock
implicit none
include 'iofile.h'
integer :: nh,istrt,isw, ii, i
real (kind=8) :: td(2)
real (kind=8) :: eps(9,*),theta(*)
! real (kind=8) :: f(3,3,*),detf(*)
real (kind=8) :: d(*),ud(*),hn(nh),h1(nh)
! Definition need for socket communication
integer, parameter :: sentreals = 8
integer, parameter :: recvreals = 42
integer(C_SIZE_T) :: bytes_recv, bytes_expt
integer(C_INT) :: curr_recv
integer(C_INT) :: status_c
real(c_double) :: sendv(sentreals) ! send 6 strains and 2 properties
real(c_double) :: sig(6), dd(6,6) ! stress and tangent modulus
! Buffers for sending and receiving raw bytes
character(kind=C_CHAR, len=sentreals * sizeof(sendv(1))),
& target :: send_buffer
character(kind=C_CHAR, len=recvreals * sizeof(sendv(1))),
& target :: recv_buffer
! Pack up the strains and the properties
sendv(1:6) = eps(1:6,1) ! epsilon_n
sendv(7) = d(1) ! mu
sendv(8) = d(2) ! lambda1
! Send the vector to the Python server
! TRANSFER intrinsic copies the raw bytes of sendv into send_buffer.
send_buffer = transfer(sendv, send_buffer)
status_c = send_c(sockfd, C_LOC(send_buffer),
& int(len(send_buffer), kind=C_SIZE_T), 0_C_INT)
if (status_c == -1) then
write(iow,*) 'Failed to send data to Python over socket.'
status_c = close_c(sockfd)
call plstop(.true.)
end if
! Receive results from the Python server
! The Python server sends 6 doubles (sum) + 36 doubles (tensor product) = 42 doubles.
! Total bytes = 42 * 8 = 336 bytes.
bytes_recv = 0_C_SIZE_T
bytes_expt = int(len(recv_buffer), kind=C_SIZE_T)
do while(bytes_recv < bytes_expt)
curr_recv = recv_c(sockfd,
& C_LOC(recv_buffer(bytes_recv+1:)),
& bytes_expt-bytes_recv, 0_C_INT)
if (curr_recv .eq. 0) then
write(iow,*) 'Python connection closed'
status_c = close_c(sockfd)
call plstop(.true.)
else if (curr_recv .eq. -1) then
write(iow,*) 'Failed to recieve Python data'
status_c = close_c(sockfd)
call plstop(.true.)
end if
bytes_recv = bytes_recv + curr_recv
end do
if (bytes_recv .ne. bytes_expt) then
write(iow,*) 'Error: expected ',bytes_expt,
& ' received ', bytes_recv
status_c = close_c(sockfd)
call plstop(.true.)
end if
! Extract the sum vector (first 6 doubles)
sig = transfer(recv_buffer(1:6*sizeof(sendv(1))), sig)
! Extract the tensor product matrix (next 36 doubles)
dd = reshape(transfer(recv_buffer(49:), dd), shape(dd))
! Socket should be closed from the client side upon exit from FEAP
! Unclear on how to do that elegantly
! status = close_c(sockfd)
! IF (status == -1) THEN
! PRINT *, "Fortran Client Error: Failed to close socket."
! ELSE
! PRINT *, "Fortran Client: Socket closed."
! END IF
end subroutine umatl3
Ipysock
Here is a simply input file if you wish to test the coding.
feap ** Test of python socket connection for user material **
0 0 0 2 2 4
mate
solid
ucon pymat
modulus mu 1000
modulus lam 1000
block
cart 5 5
1 0 0
2 1 0
3 1 1
4 0 1
eboun
1 0 1 0
cboun
node 0 0 1 1
csurf
linear
1 1 0 2.0
2 1 1 2.0
end
inte
stop