# fortran_shapefuncs > 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: | Constant | Element | Reference domain | Nodes, order p | Orders | |---|---|---|---|---| | `CUB_LIN` | line | [-1, 1] | p+1 | >= 0 | | `CUB_QUA` | quadrilateral | [-1, 1]^2 | (p+1)^2 | >= 0 | | `CUB_HEX` | hexahedron | [-1, 1]^3 | (p+1)^3 | >= 0 | | `CUB_TRI` | triangle | x, y >= 0, x + y <= 1 | (p+1)(p+2)/2 | >= 1 | | `CUB_TET` | tetrahedron | x, y, z >= 0, x + y + z <= 1 | (p+1)(p+2)(p+3)/6 | >= 1 | | `CUB_WED` | wedge (prism) | triangle x [-1, 1] | (p+1)^2(p+2)/2 | >= 1 | | `CUB_QUA`, serendipity | quadrilateral | [-1, 1]^2 | 4p | 1 to 3 | | `CUB_HEX`, serendipity | hexahedron | [-1, 1]^3 | 8 + 12(p-1) | 1 to 3 | 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 | Call | Result | |---|---| | `type(shapefunc)` | an unset object: `elm = 0`, arrays not allocated, `is_valid()` false | | `shapefunc(q, order)` | Lagrange shape functions of order `order` on element `q%elm` at the points of the cubature `q` | | `shapefunc(q, order, infin)` | with infinitude `infin(:)` per direction: size 1 (every direction) or `q%dim` | | `shapefunc(q, order, family=SHP_SERENDIPITY)` | serendipity quadrilateral or hexahedron, order 1 to 3 | | `shapefunc(q, order, nodes=perm)` | caller's node numbering: node `k` is node `perm(k)` of the default numbering | 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 | Constant | Meaning | |---|---| | `SHP_FIN` | finite direction: Lagrange polynomials on [-1, 1] (default) | | `SHP_INF` | infinite direction: mapping functions, to interpolate coordinates | | `SHP_CHP` | infinite direction: chopped functions, to interpolate a solution that vanishes at infinity | | `SHP_LAGRANGE` | family: full Lagrange elements (default) | | `SHP_SERENDIPITY` | family: serendipity quadrilaterals and hexahedra, orders 1 to 3 | ## Components All public; read them directly and set them only through the constructor or `set`. | Component | Meaning | |---|---| | `s%elm` | element type, a `CUB_*` constant; 0 if unset | | `s%dim` | spatial dimension, 1 to 3 | | `s%order` | polynomial order p | | `s%family` | `SHP_LAGRANGE` or `SHP_SERENDIPITY` | | `s%nnodes` | number of nodes | | `s%npoints` | number of cubature points | | `s%infin(3)` | infinitude per direction, `SHP_*`; unused entries are `SHP_FIN` | | `s%lattice(:, nnodes)` | integer lattice indices of each node (see "Node numbering") | | `s%coords(dim, nnodes)` | `real64` nodal coordinates in the reference element | | `s%func(nnodes, npoints)` | `real64` values N_i at each point | | `s%derv(dim, nnodes, npoints)` | `real64` first derivatives dN_i/dxi_j | | `s%curv(dim, dim, nnodes, npoints)` | `real64` second derivatives d2N_i/dxi_j dxi_k, symmetric in j, k | ## Type-bound procedures | Call | Result | |---|---| | `call s%set(q, order [, infin] [, family] [, nodes])` | builds the shape functions in place, like the constructor, discarding previous contents | | `call s%eval(xi, func, derv, curv)` | evaluates at any point `xi(dim)` of the reference element into `func(nnodes)`, `derv(dim, nnodes)`, `curv(dim, dim, nnodes)`; the object must be set; infinite directions are singular at xi = 1 | | `s%is_valid()` | logical: the object is set and its array shapes are consistent | | `call s%summary([unit])` | writes element, dimension, order, family, infinitude, numbers of nodes and points (default `output_unit`) | | `call s%show([unit])` | writes the summary, then the value and first derivatives for every point and node, and per point the sums over the nodes (1 and 0 unless a direction is chopped) | | `call s%numbering([unit])` | writes the summary, then the number and reference coordinates of every node | | `call s%destroy()` | deallocates and resets to the unset state | ## 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: | Element | Edges | |---|---| | quadrilateral | 1-2, 2-3, 3-4, 4-1 | | hexahedron | 1-2, 2-3, 3-4, 4-1, 5-6, 6-7, 7-8, 8-5, 1-5, 2-6, 3-7, 4-8 | | triangle | 1-2, 2-3, 3-1 | | tetrahedron | 1-2, 2-3, 3-1, 1-4, 2-4, 3-4 | | wedge | 1-2, 2-3, 3-1, 4-5, 5-6, 6-4, 1-4, 2-5, 3-6 | - 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. ```fortran 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 ``` ```fortran 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 ``` ```fortran 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 ``` ```fortran 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 ``` ```fortran 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 ``` ```fortran 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 ```