User manual

Fortran module shapefuncs: the derived type shapefunc, which holds the shape functions of a finite element and their first and second derivatives at the points of a cubature from fortran_cubatures. Lines, triangles, quadrilaterals, tetrahedra, hexahedra and wedges of any order, serendipity quadrilaterals and hexahedra, mapped infinite elements, and a node numbering that the caller can permute. Fortran 2018, fpm or make. This file is the complete API reference for tools and coding agents; it is checked against the source in CI.

Use: use shapefuncs (everything) or use shapefuncs, only: shapefunc, SHP_INF, ..., together with use cubatures for cubature and the element constants CUB_LIN ... CUB_WED. fpm dependency: fortran_shapefuncs = { git = "https://github.com/willklausler/fortran_shapefuncs", tag = "v0.2.0" } (it pulls in fortran_cubatures v0.2.0). Kinds are not exported: values and coordinates are real64 from iso_fortran_env.

Conventions (read first)

• A shapefunc is built from a cubature q and an order p. The element is q%elm, and the functions are evaluated at the points q%abscissae. Point g and node i index the arrays as s%func(i,g) = N_i(xi_g), s%derv(j,i,g) = dN_i/dxi_j, s%curv(j,k,i,g) = d2N_i/dxi_j dxi_k (symmetric in j, k).
• Reference elements are those of fortran_cubatures:

Common node counts: order 1: line 2, quad 4, hex 8, triangle 3, tet 4, wedge 6; order 2: 3, 9, 27, 6, 10, 18 (serendipity quad 8, hex 20); order 3: 4, 16, 64, 10, 20, 40 (serendipity quad 12, hex 32). - Nodes are equispaced. There is no upper limit on the order, but interpolation becomes ill-conditioned at high orders (Runge's phenomenon). The tests cover orders up to 8 for lines and triangles, 6 for quadrilaterals and tetrahedra, 5 for wedges and 4 for hexahedra. - Node numbering is hierarchical: vertices, then edge nodes, then face nodes, then interior nodes. Orders 1 to 3 give the classical numbering of the 2-, 3-, 4-node line, 4-, 9-, 16-node quadrilateral, 8-, 27-node hexahedron, 3-, 6-, 10-node triangle, 4-, 10-, 20-node tetrahedron and 6-, 18-node wedge. See "Node numbering" for the rules. The nodes argument of set replaces it with any other numbering. - Infinitude is chosen per direction with SHP_FIN (finite, the default), SHP_INF (mapping functions, for the geometry) or SHP_CHP (chopped functions, for the solution). An infinite direction extends to infinity at xi = 1. See "Infinite elements". - Error policy: invalid arguments to the constructor or set stop with error stop and a message that starts with shapefunc%set:. Errors cannot be caught. set stops when - q is not set; - order is negative, or below 1 for triangles, tetrahedra, wedges or any infinite direction; - infin has a size other than 1 or q%dim, or a value other than SHP_FIN, SHP_INF, SHP_CHP; - a triangle or tetrahedron has an infinite direction, or a wedge has one other than direction 3; - family is not SHP_LAGRANGE or SHP_SERENDIPITY, or a serendipity element is not a finite quadrilateral or hexahedron of order 1 to 3; - nodes does not have nnodes entries or is not a permutation of 1:nnodes. - No module state. The constructor, set, eval, is_valid and destroy are pure; only summary, show and numbering (output) are not. A set shapefunc can be shared read-only between threads.

Types and constructors

The optional arguments infin, family and nodes combine. Assignment between shapefunc objects is intrinsic (deep copy). The shape functions depend only on the reference element, so build one object per element type and reuse it for every element of the mesh.

Constants

Components

All public; read them directly and set them only through the constructor or set.

Type-bound procedures

Shape functions

• Lines, quadrilaterals, hexahedra: tensor products of the Lagrange polynomials of degree p on the equispaced nodes -1 + 2a/p, a = 0..p. Order 0 is one constant function with its node at the centre.
• Triangles, tetrahedra: the complete polynomials of total degree p. With barycentric coordinates lambda (for a triangle lambda = (xi_1, xi_2, 1 - xi_1 - xi_2)), the node at lambda = a/p has N = prod_k P_{a_k}(lambda_k), where P_a(l) = prod_{j=0}^{a-1} (p l - j)/(j + 1) (Silvester's polynomials).
• Wedges: the triangle functions in (xi_1, xi_2) times the line functions in xi_3.
• Serendipity quadrilaterals and hexahedra (family=SHP_SERENDIPITY, orders 1 to 3): nodes on the vertices and edges only. They span the monomials of superlinear degree <= p, the degree that counts only the variables of exponent 2 or more (Arnold and Awanou). For p = 2 these are the classical 8-node quadrilateral, with corner and midside functions N = (1 + x x_i)(1 + y y_i)(x x_i + y y_i - 1)/4 and N = (1 - x^2)(1 + y y_i)/2, and the 20-node hexahedron. set computes each serendipity function once as a fixed combination of the tensor-product Lagrange functions of the same order.
• Every shape function is a product of one-dimensional functions of factor coordinates that are affine in xi (the coordinates themselves, or the barycentric coordinates), so the derivatives are exact: they follow from the product rule and a constant chain-rule matrix.

Node numbering

• Vertices. Quadrilateral: counter-clockwise from (-1, -1). Hexahedron: the bottom face xi_3 = -1 counter-clockwise from (-1, -1, -1), then the top face in the same order. Triangle and tetrahedron: vertex k is where barycentric coordinate k is 1, so vertex 1 is (1, 0) or (1, 0, 0) and the last vertex is the origin. Wedge: the triangle vertices at xi_3 = -1, then at xi_3 = 1.
• Edges, with their nodes running from the first vertex of the edge to the second:

• Faces. Hexahedron: xi_3 = -1, xi_3 = 1, xi_1 = -1, xi_1 = 1, xi_2 = -1, xi_2 = 1; a face in (xi_j, xi_k), j < k, is numbered as a quadrilateral in those coordinates. Tetrahedron: 1-2-3, 1-2-4, 2-3-4, 1-3-4, each numbered as a triangle with its vertices in that order. Wedge: the bottom and top triangles, then the quadrilaterals on edges 1-2, 2-3, 3-1, each numbered in (position along the edge, xi_3).
• Interiors of quadrilaterals, hexahedra, triangles and tetrahedra (and of their faces) are numbered recursively as the same element of order p-2 (quadrilateral, hexahedron), p-3 (triangle) or p-4 (tetrahedron). Wedge interior nodes are triangle interiors, layer by layer from xi_3 = -1.
• Serendipity elements keep the vertex and edge nodes of the Lagrange element in the same order: the 8-node quadrilateral is the 9-node one without its centre, the 20-node hexahedron is the 27-node one without its face and centre nodes.
• The 16-node quadrilateral (p = 3, xi_1 to the right, xi_2 up) and the 15-node triangle (p = 4, vertex 3 at the origin):

text 4 --10 -- 9 -- 3 2 | | | \ 11 16 15 8 7 6 | | | \ 12 13 14 7 8 14 5 | | | \ 1 -- 5 -- 6 -- 2 9 15 13 4 | \ 3 --10--11--12-- 1

• Lattice indices s%lattice(:, i) give the integer position of node i: one index 0..p per direction for lines, quadrilaterals and hexahedra; the barycentric indices a_k, which sum to p, for triangles and tetrahedra; three barycentric indices then the axial index 0..p for wedges. They locate nodes without comparing floating-point coordinates: for example a triangle node is on the boundary where one of its indices is 0.
• Caller-defined numbering: mesh formats number higher-order nodes differently. With nodes=perm, the caller's node k is node perm(k) of the default numbering, and every output (lattice, coords, func, derv, curv, eval) follows the caller's numbering at no extra cost per evaluation. Build perm by matching reference coordinates: perm(k) = minloc(sum(abs(s%coords - spread(xmine(:,k), 2, s%nnodes)), dim=1), 1) for a default s and the caller's node coordinates xmine(:,k) (see demo_serendipity).

Infinite elements

Mapped infinite elements model unbounded domains (soil under a foundation, a far field, a body in free space): one layer of them surrounds the finite mesh. They are the elements of Bettess in the form of Zienkiewicz, Emson and Bettess, generalised to order p as by Marques and Owen.

• Any direction of a line, quadrilateral or hexahedron, and direction 3 of a wedge, can be infinite. Triangles and tetrahedra cannot. An infinite direction of order p >= 1 has p+1 nodes at t_a = -1 + 2a/(p+1), a = 0..p: the grid of order p+1 without the node t_{p+1} = 1 at infinity (linear: -1, 0; quadratic: -1, -1/3, 1/3).
• An infinite element uses two shapefunc objects at the same cubature points and with the same nodes: SHP_INF for the geometry and SHP_CHP for the field.
• Mapping functions (SHP_INF): M_a(xi) = l_a(xi) (1 - t_a)/(1 - xi), with l_a the Lagrange polynomials on t_0..t_p. They are 1 at their own node, sum to 1, reproduce q(xi)/(1 - xi) for q of degree <= p, and are singular at xi = 1. Linear: M_0 = -2 xi/(1 - xi), M_1 = (1 + xi)/(1 - xi).
• The map x(xi) = sum_a M_a(xi) x_a sends xi -> 1 to infinity. In one dimension with nodes x_0, x_1 it is x = x_P + 2(x_1 - x_0)/(1 - xi) with the pole x_P = 2 x_0 - x_1: place the nodes on rays from a pole inside the finite domain, with x_1 - x_0 = x_0 - x_P. Then 1 - xi is proportional to 1/r, r the distance from the pole. For p = 2 the nodes sit at distances from the pole in the ratio 2 : 3 : 6.
• Chopped functions (SHP_CHP): the Lagrange polynomials of order p+1 on t_0..t_{p+1} = 1 without the function of the node at infinity. They vanish at xi = 1 and do not sum to 1. Because 1 - xi ~ 1/r, they represent u(r) = a_1/r + a_2/r^2 + ... + a_{p+1}/r^{p+1}. For a field that tends to a nonzero constant, subtract the constant first.
• Integration: the integrands are rational, so no rule is exact in general. Gauss-Legendre rules keep their points away from xi = 1; raise the degree until the result converges. With u ~ 1/r, the 1/(1 - xi)^2 of the Jacobian and the (1 - xi)^2 of u^2 often cancel to a polynomial (see demo_infinite, which is exact).

Choosing the cubature

For an affine element of order p, the mass matrix needs a cubature of degree 2p and the stiffness matrix 2p - 2. Curved and distorted elements need more, and infinite elements are not polynomial (see above). In an element loop with nodal coordinates xn(dim, nnodes), the Jacobian at point g is matmul(xn, transpose(s%derv(:,:,g))) and the physical gradients are J^(-T) times s%derv(:,:,g) (see demo_stiffness).

Compiler notes

Tested in CI with gfortran 14 (Linux, macOS, Windows), gfortran 13 (Linux), ifx 2026 (Linux, Windows), LLVM flang 21 and nvfortran 26 (Linux). No compiler-specific code and no known compiler limitations.

Differs from

• FEM_BASIS (John Burkardt, Fortran 90): Lagrange basis functions of any degree on simplices of any dimension, evaluated point by point; no quadrilaterals, hexahedra, wedges, serendipity or infinite elements, no node-numbering convention for meshes, no fpm package.
• FEM2D_PACK (John Burkardt): low-order 2D triangles and quadrilaterals only.
• FELIPE (University of Bath, Fortran 77): a teaching FE package whose programs include 1D and 2D mapped infinite elements, built into the solvers rather than offered as a library.
• Solvers written in Fortran, such as Elmer and OpenRadioss, have shape functions inside the solver that cannot be used on their own.
• Basix (FEniCS, C++ and Python) tabulates many element families on reference cells, which are [0, 1]-based, with a different numbering; it has no infinite elements. Map the coordinates and permute the nodes (nodes=) if you mix them.

Module layout

One module, shapefuncs, in src/shapefuncs.f90. It uses the module cubatures from fortran_cubatures; the make build compiles both into libfortran_shapefuncs.a. Nothing is generated.

Examples

These are the programs in example/, compiled and run in CI.

program demo_stiffness
  !! The element loop of a finite element code: the stiffness matrix of the
  !! Laplacian on a bilinear quadrilateral, from the isoparametric map and
  !! the physical gradients dN/dx = J**(-T) dN/dxi
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cubatures
  use shapefuncs
  implicit none
  real(dp), parameter :: xn(2,4) = reshape([0._dp, 0._dp, 1._dp, 0._dp, &
                                          & 1._dp, 1._dp, 0._dp, 1._dp], [2, 4])  ! unit square
  type(cubature) :: q
  type(shapefunc) :: s
  real(dp) :: k(4,4), jac(2,2), jinv(2,2), detj, dndx(2,4)
  integer :: g

  q = cubature(CUB_QUA, 2)                     ! degree 2p for order p = 1
  s = shapefunc(q, 1)
  k = 0
  do g = 1, q%npoints
    jac  = matmul(xn, transpose(s%derv(:,:,g)))
    detj = jac(1,1)*jac(2,2) - jac(1,2)*jac(2,1)
    jinv = reshape([jac(2,2), -jac(2,1), -jac(1,2), jac(1,1)], [2, 2])/detj
    dndx = matmul(transpose(jinv), s%derv(:,:,g))
    k = k + matmul(transpose(dndx), dndx)*detj*q%weights(g)
  end do

  print '(4f9.5)', k(1,:)                      ! 0.66667 -0.16667 -0.33333 -0.16667
  print '(a,es8.1)', 'max row sum:', maxval(abs(sum(k, dim=2)))  ! 0 up to rounding
end program demo_stiffness

program demo_isoparametric
  !! Area of the quarter annulus 1 <= r <= 2, 0 <= theta <= pi/2, meshed by
  !! one 9-node quadrilateral whose nodes lie on the exact geometry
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cubatures
  use shapefuncs
  implicit none
  real(dp), parameter :: pi = acos(-1._dp)
  type(cubature) :: q
  type(shapefunc) :: s
  real(dp) :: xn(2,9), jac(2,2), area, r, t
  integer :: g, n

  q = cubature(CUB_QUA, 6)
  s = shapefunc(q, 2)

  do n = 1, s%nnodes                           ! xi_1 -> radius, xi_2 -> angle
    r = 1.5_dp + 0.5_dp*s%coords(1,n)
    t = pi/4*(1 + s%coords(2,n))
    xn(:,n) = r*[cos(t), sin(t)]
  end do

  area = 0                                     ! sum_g w_g det(dx/dxi)
  do g = 1, q%npoints
    jac = matmul(xn, transpose(s%derv(:,:,g))) ! dx/dxi = sum_n x_n dN_n/dxi
    area = area + (jac(1,1)*jac(2,2) - jac(1,2)*jac(2,1))*q%weights(g)
  end do

  print '(a,f12.8)', 'isoparametric:', area    ! 2.32842712
  print '(a,f12.8)', 'exact 3*pi/4: ', 3*pi/4  ! 2.35619449
end program demo_isoparametric

program demo_table
  !! Keep one set of quadratic shape functions per element type, with a
  !! matching cubature, in a table indexed by the element constants
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cubatures
  use shapefuncs
  implicit none
  type(shapefunc) :: table(6)
  integer :: elm

  do elm = CUB_LIN, CUB_WED                    ! CUB_* are 1 to 6
    table(elm) = shapefunc(cubature(elm, 4), 2)
  end do

  call table(CUB_WED)%summary()                ! WED, order 2, 18 nodes, 18 points
  print '(a,6i4)', 'nodes: ', table%nnodes     ! 3   6   9  10  27  18
  print '(a,6i4)', 'points:', table%npoints    ! 3   6   9  27  27  18
end program demo_table

program demo_numbering
  !! Print the node numbering of the quadratic tetrahedron and check that
  !! each shape function is 1 at its own node and 0 at the others
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cubatures
  use shapefuncs
  implicit none
  type(shapefunc) :: s
  real(dp) :: func(10), derv(3,10), curv(3,3,10), err
  integer :: i

  s = shapefunc(cubature(CUB_TET, 2), 2)       ! 10-node tetrahedron
  call s%numbering()                           ! vertices 1-4, then edge midpoints 5-10

  err = 0
  do i = 1, s%nnodes
    call s%eval(s%coords(:,i), func, derv, curv)
    func(i) = func(i) - 1
    err = max(err, maxval(abs(func)))
  end do
  print '(a,es8.1)', 'max |N_j(x_i) - delta_ij| =', err  ! 0.0E+00
end program demo_numbering

program demo_serendipity
  !! Use the 20-node serendipity hexahedron in the node ordering of a mesh
  !! format that lists the edge nodes by their first vertex, as Gmsh does.
  !! The permutation for `nodes` comes from matching reference coordinates.
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cubatures
  use shapefuncs
  implicit none
  integer, parameter :: edges(2,12) = reshape([1,2, 1,4, 1,5, 2,3, 2,6, 3,4, &
                                             & 3,7, 4,8, 5,6, 5,8, 6,7, 7,8], [2, 12])
  type(cubature) :: q
  type(shapefunc) :: s
  real(dp) :: xmesh(3,20)
  integer :: perm(20), k

  q = cubature(CUB_HEX, 4)
  s = shapefunc(q, 2, family=SHP_SERENDIPITY)

  xmesh(:,1:8) = s%coords(:,1:8)               ! same vertices, then the mesh's edges
  do k = 1, 12
    xmesh(:,8+k) = (s%coords(:,edges(1,k)) + s%coords(:,edges(2,k)))/2
  end do

  do k = 1, 20                                 ! mesh node k is node perm(k) of s
    perm(k) = minloc(sum(abs(s%coords - spread(xmesh(:,k), 2, 20)), dim=1), 1)
  end do
  print '(a,20i3)', 'nodes =', perm            ! 1 ... 8  9 12 17 10 18 11 19 20 13 16 14 15

  s = shapefunc(q, 2, family=SHP_SERENDIPITY, nodes=perm)
  print '(a,l2)', 'same order:', all(abs(s%coords - xmesh) < 1e-12_dp)  ! T
end program demo_serendipity

program demo_infinite
  !! Integral of u**2 from 1 to infinity for u = 1/x, which is 1, on one
  !! infinite line element with nodes at x = 1 and x = 2. Linear mapping
  !! functions give x(xi) = 2/(1 - xi), so xi -> 1 is x -> infinity; linear
  !! chopped functions interpolate u and represent 1/x exactly.
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cubatures
  use shapefuncs
  implicit none
  real(dp), parameter :: xn(2) = [1._dp, 2._dp]   ! nodal coordinates
  type(cubature) :: q
  type(shapefunc) :: map, field
  real(dp) :: un(2), u, dxdxi, total
  integer :: g

  q = cubature(CUB_LIN, 2)
  map = shapefunc(q, 1, [SHP_INF])                ! geometry
  field = shapefunc(q, 1, [SHP_CHP])              ! solution, vanishes at infinity
  un = 1/xn

  total = 0
  do g = 1, q%npoints
    dxdxi = dot_product(map%derv(1,:,g), xn)
    u = dot_product(field%func(:,g), un)
    total = total + u**2*dxdxi*q%weights(g)
  end do

  print '(a,f12.8)', 'infinite element:', total   ! 1.00000000
end program demo_infinite