How to Apply Surface Loads

Original Post - 12 Jul 2025 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

Applying surface loads (pressure loads) to solid and shell elements in OpenSees is difficult. The typical approach is to use equivalent nodal loads, but that’s intractable for anything beyond simple meshes. Do you want to calculate equivalent nodal loads (in your script, after the model is defined) for a triangulated mesh? Yeah, me neither.

It would be better to have OpenSees compute equivalent nodal loads for you and your model. Such capability has been in OpenSees for some time, but it’s not well known and the implementation has been limited.

Chris McGann implemented the SurfaceLoad class as a subclass of Element. As described on its wiki page, the SurfaceLoad element applies “energetically-conjugate forces” on four nodes, e.g., normal to one face of a brick element or normal to a quadrilateral shell element. “Energetically-conjugate forces” is another way of saying the element uses bilinear shape functions to compute consistent nodal loads from uniform pressure.

A copy-paste-modify implementation, TriSurfaceLoad, is suited to one face of a tetrahedral element or a triangular shell element. For this element, the energetically-conjugate forces are based on the shape functions of a constant strain triangle.

Pressure and nodal loads on quad and triangular elements
Pressure loading shown here is negative, opposite to the direction produced by the right hand rule around nodes I-J-K-L and I-J-K

After I write this post, I will make a two node LineSurfaceLoad element suitable for applying pressure loads on the edges of quad and triangle elements in 2D analysis as well as on the edges of shells.

Anyway, at first glance, this surface-load-as-an-element approach may appear bloated because the element only provides unbalanced nodal loads and does not return stiffness, mass, or damping. However, SurfaceLoad and TriSurfaceLoad objects act like wrappers that are agnostic to the elements on which the surface load is applied, making the implementations quite general.

Minimal Examples

You can use SurfaceLoad or TriSurfaceLoad without any underlying elements at all. Below is an overly minimal example that makes the point. The pressure is defined as positive according to the normal vector implied by the right hand rule going around the horn of I-J-K-L nodes for SurfaceLoad and the I-J-K nodes for TriSurfaceLoad.

import openseespy.opensees as ops
from math import isclose

ops.wipe()
ops.model('basic','-ndm',3,'-ndf',3)

ops.node(1,0,0,0); ops.fix(1,1,1,1)
ops.node(2,1,0,0); ops.fix(2,1,1,1)
ops.node(3,1,1,0); ops.fix(3,1,1,1)
ops.node(4,0,1,0); ops.fix(4,1,1,1)

meshType = 'quad'
# meshType = 'tri'

if meshType == 'quad':
    #                        tag I J K L   p
    ops.element('SurfaceLoad',1, 1,2,3,4,-1.0)
if meshType == 'tri':
    #                           tag I J K   p
    ops.element('TriSurfaceLoad',1, 1,2,3,-1.0)
    ops.element('TriSurfaceLoad',2, 1,3,4,-1.0)

ops.system('UmfPack')
ops.analysis('Static','-noWarnings')
ops.analyze(1)
ops.reactions()

if meshType == 'quad':
    for nd in [1,2,3,4]:
        assert isclose(0.25,ops.nodeReaction(nd,3))
if meshType == 'tri':
    for nd in [1,3]:
        assert isclose(1/3,ops.nodeReaction(nd,3))
    for nd in [2,4]:
        assert isclose(1/6,ops.nodeReaction(nd,3))

The assertion shows the pressure load becomes reactions on the fixed nodes. For the quadrilateral elements, the pressure is evenly distributed at the four nodes. With triangular elements, more load is applied to the nodes shared between the two elements. The assertions verify these conditions are satisfied.

If the nodes were not in a unit square, the reactions would adjust according to the nodal locations based on the shape functions in the surface load elements.

When the surface loads are defined outside a load pattern, the load factor applied to the surface loads is always 1.0. To get time varying surface loads, inside a time series/load pattern combination you can associate an element load to surface load elements. In the minimal example below, a single brick element is defined with a surface load on its top face. An element load of type surfaceLoad is associated with a load pattern/time series combo.

import openseespy.opensees as ops
from math import isclose

E = 29000
v = 0.3

ops.wipe()
ops.model('basic','-ndm',3,'-ndf',3)

ops.node(1,0,0,0); ops.fix(1,1,1,1)
ops.node(2,1,0,0); ops.fix(2,1,1,1)
ops.node(3,1,1,0); ops.fix(3,1,1,1)
ops.node(4,0,1,0); ops.fix(4,1,1,1)
ops.node(5,0,0,1)
ops.node(6,1,0,1)
ops.node(7,1,1,1)
ops.node(8,0,1,1)

ops.nDMaterial('ElasticIsotropic',1,E,v)
ops.element('stdBrick',1, 1,2,3,4,5,6,7,8, 1)

ops.timeSeries('Linear',1)
ops.pattern('Plain',1,1)
ops.element('SurfaceLoad',2, 5,6,7,8, -1.0)
ops.eleLoad('-ele',2,'-type','-surfaceLoad')

ops.system('UmfPack')
ops.analysis('Static','-noWarnings')

ops.analyze(10)
ops.reactions()

lam = ops.getLoadFactor(1)
for nd in [1,2,3,4]:
    assert isclose(0.25*lam,ops.nodeReaction(nd,3))

Remember that the eleLoad command maps to the surface load element, not to the element whose nodes will be loaded.

Surface Loads on Solid Elements

Whether using SurfaceLoad elements to apply pressure loading on a mesh of brick elements or TriSurfaceLoad elements to apply pressure loading on a tetrahedral mesh, using surface load elements can be difficult for two reasons:

You have to figure out which solid elements have a face on the loaded surface and then the three or four element nodes on that face.
The I-J-K-L or I-J-K node ordering on the face of each element of the loaded surface can give a normal vector that is either “in” or “out” of the surface.

You may be able to circumvent these issues if you use a third party mesh generator, but I’ll focus on the mesh commands available in OpenSees.

First, define a function that returns a list of element nodes that are on a surface.

def nodesOnSurface(eleTag, surfaceNodes):
'''
Return a list of element nodes that appear in the list
of surface nodes
'''
    eleNodes = ops.eleNodes(eleTag)

    nodes = []
    for nd in eleNodes:
        if nd in surfaceNodes:
            nodes.append(nd)

    return nodes

Yeah, a list comprehension would be clever, but an if statement embedded in a for loop gets the job done.

Also define a function that computes the normal vector based on three nodes.

def normal(ndI, ndJ, ndK):
'''
Compute the normal vector (cross product) formed by
position vectors I-J and I-K
'''
    xyzI = ops.nodeCoord(ndI)
    xyzJ = ops.nodeCoord(ndJ)
    xyzK = ops.nodeCoord(ndK)

    e1 = np.subtract(xyzJ,xyzI) # e1 = J - I
    e2 = np.subtract(xyzK,xyzI) # e2 = K - I

    return np.cross(e1,e2) # e1 x e2

Consider the mesh of tetrahedral elements shown in this post. I won’t repeat the model definition here, but let’s apply a pressure loading to the top surface of the mesh.

As shown in the previous post, the mesh of nodes on the top surface has tag 15 while the mesh of tetrahedral elements has tag 20.

p = 1.0 # reference pressure
ops.timeSeries('Linear',1)
ops.pattern('Plain',1,1)

eleTag = max(ops.getEleTags())
topNodes = ops.getNodeTags('-mesh',15)
for ele in ops.getEleTags('-mesh',20):
    eleTag += 1
    nodesToLoad = nodesOnSurface(ele, topNodes)
    c = len(nodesToLoad)
    if c != 3: # This element does not have three nodes on surface
        continue
    n = normal(*nodesToLoad)

    ops.element('TriSurfaceLoad',eleTag,*nodesToLoad,p*np.sign(n[2]))
    ops.eleLoad('-ele',eleTag,'-type','-surfaceLoad')

Minor modifications to the above script are necessary in order to apply surface loads to brick elements having four nodes on the loaded surface.

Surface Loads on Shell Elements

Surface loads on shell elements work only after PR #1634, and will be available in versions of OpenSees.exe after 3.7.1 and after OpenSeesPy 3.7.1.2.

Using the surface load elements is a little more straightforward with shell models because the element meshes are planar. As a result, we don’t need to determine element nodes on a surface, but we still need the normal vectors to ensure the pressure loads are applied in the correct direction.

Consider the shell mesh shown in this post. Again, I won’t repeat the model definition here, but let’s assume we want to apply pressure loading to the top surface of the mesh.

As shown in the previous post, the shell elements are defined in two meshes with tags 10 and 11. I used two meshes to make a different point in the previous post, so don’t take the following code too literally.

p = 1.0 # reference pressure
ops.timeSeries('Linear',1)
ops.pattern('Plain',1,1)

eleTag = max(ops.getEleTags())
for meshTag in [10,11]:
    for ele in ops.getEleTags('-mesh',meshTag):
        eleTag += 1
        nodesToLoad = ops.eleNodes(ele)
        n = normal(*nodesToLoad)

        ops.element('TriSurfaceLoad',eleTag,*nodesToLoad,p*np.sign(n[2]))
        ops.eleLoad('-ele',eleTag,'-type','-surfaceLoad')

Again, only minor modifications are required in order to use this script for meshes of four noded shell elements.

Yeah, I agree–using surface load elements to apply pressure loading to OpenSees models is not so straightforward. But the approach shown in this press sure beats defining consistent nodal loads yourself.

OpenSees Cloud

OpenSees AMI

How to Apply Surface Loads

Minimal Examples

Surface Loads on Solid Elements

Surface Loads on Shell Elements