Fortran_cubatures and fortran_shapefuncs: FE quadrature and shape functions of any order

Announcements
1

Greetings,

I’d like to share two companion fpm packages for writing finite element codes:

  • fortran_cubatures: Gaussian integration points and weights on the six reference elements (line, triangle, quadrilateral, tetrahedron, hexahedron, wedge) for any degree of exactness. No dependencies.
  • fortran_shapefuncs: shape functions and their first and second derivatives at those points.
use cubatures
use shapefuncs
type(cubature)  :: q
type(shapefunc) :: s

q = cubature(CUB_HEX, 4)                     ! integration points
s = shapefunc(q, 2)                          ! 27-node hexahedron at those points
print *, s%nnodes, sum(s%func(:,1))          ! 27  1.0
s = shapefunc(q, 2, family=SHP_SERENDIPITY)  ! 20 nodes

Packages’ contents

  • Any degree or order, with separate degrees per direction on tensor-product elements. Wedges included.
  • All weights positive, all points inside the element. Tests check that each rule integrates every monomial up to its degree exactly.
  • Lagrange elements of any order, serendipity elements, and mapped infinite elements for unbounded domains. Node numbering is the classical one or set at run time.
  • Pure procedures, no module state.
  • CI on Linux with gfortran, ifx, LLVM flang and nvfortran, plus gfortran on macOS and gfortran and ifx on Windows.

As far as I could find, nothing else in Fortran covers this: stdlib and quadpack are one-dimensional, triasymq is FORTRAN 77 without line or wedge rules, and solvers like Elmer keep their shape functions internal. Corrections are welcome if I missed something.

How they were made
I originally wrote these modules for use in FEAP for my research. For this release, I used Claude Code extensively to generalize both modules to arbitrary element order, add serendipity elements, expand the tests, and polish them for GitHub release. I’ve started verifying the packages with an FE modal analysis project, also on my page. The code hasn’t been independently reviewed, so I’d especially welcome eyes on the higher order elements.

Each repo also has an AGENTS.md with instructions for coding agents working on the package, and an llms.txt with the whole API in one plain-text file (it’s also the source of the user manual).

Install (shapefuncs pulls in cubatures):

[dependencies]
fortran_shapefuncs = { git = "https://github.com/willklausler/fortran_shapefuncs", tag = "v0.2.0" }

Docs: cubatures · shapefuncs
DOI: 10.5281/zenodo.23300827 · 10.5281/zenodo.23301496

Question: would hierarchical bases, or node numberings that match Gmsh or VTK, be useful to you?

5 Likes
2 
userland@localhost:~$ gfortran features.f -ffree-form
userland@localhost:~$ ./a.out
Quadrilateral shape functions and cubature
            x            y       det(J)       sum(N)
     0.422650     0.211325     0.500000     1.000000
     1.577350     0.211325     0.500000     1.000000
     0.422650     0.788675     0.500000     1.000000
     1.577350     0.788675     0.500000     1.000000

Computed area:             2.000000000
Exact area:                2.000000000
Integral of x^2 + y^2:     3.333333333
Exact integral:            3.333333333
Maximum sum(N) error:     1.1102E-16

userland@localhost:~$ cat features.f

program cubature_demo
  use iso_fortran_env, only: real64
  implicit none

  integer, parameter :: dp = real64

  real(dp) :: xn(4), yn(4), gp(2)
  real(dp) :: n(4), dxi(4), deta(4)
  real(dp) :: xi, eta, x, y, detj
  real(dp) :: xxi, xeta, yxi, yeta
  real(dp) :: area, integral, error_sum
  integer  :: i, j

  ! Counterclockwise nodes of a 2-by-1 rectangle.
  xn = [0.0_dp, 2.0_dp, 2.0_dp, 0.0_dp]
  yn = [0.0_dp, 0.0_dp, 1.0_dp, 1.0_dp]

  ! Two-point Gauss rule: both weights are 1.
  gp = [-1.0_dp, 1.0_dp] / sqrt(3.0_dp)

  area      = 0.0_dp
  integral  = 0.0_dp
  error_sum = 0.0_dp

  print '(a)', 'Quadrilateral shape functions and cubature'
  print '(a)', '            x            y       det(J)       sum(N)'

  do j = 1, 2
    eta = gp(j)
    do i = 1, 2
      xi = gp(i)

      call shapefunc(xi, eta, n, dxi, deta)

      ! Map reference coordinates to physical coordinates.
      x = dot_product(n, xn)
      y = dot_product(n, yn)

      ! Jacobian of the coordinate transformation.
      xxi  = dot_product(dxi,  xn)
      xeta = dot_product(deta, xn)
      yxi  = dot_product(dxi,  yn)
      yeta = dot_product(deta, yn)
      detj = xxi*yeta - xeta*yxi

      if (detj <= 0.0_dp) error stop 'Invalid element Jacobian'

      area      = area + detj
      integral  = integral + (x*x + y*y)*detj
      error_sum = max(error_sum, abs(sum(n) - 1.0_dp))

      print '(4f13.6)', x, y, detj, sum(n)
    end do
  end do

  print *
  print '(a,f14.9)',  'Computed area:          ', area
  print '(a,f14.9)',  'Exact area:             ', 2.0_dp
  print '(a,f14.9)',  'Integral of x^2 + y^2:  ', integral
  print '(a,f14.9)',  'Exact integral:         ', 10.0_dp/3.0_dp
  print '(a,es12.4)', 'Maximum sum(N) error:   ', error_sum

contains

  subroutine shapefunc(xi, eta, n, dxi, deta)
    real(dp), intent(in)  :: xi, eta
    real(dp), intent(out) :: n(4), dxi(4), deta(4)
    real(dp), parameter :: sx(4) = [-1.0_dp,  1.0_dp, 1.0_dp, -1.0_dp]
    real(dp), parameter :: sy(4) = [-1.0_dp, -1.0_dp, 1.0_dp,  1.0_dp]

    n    = 0.25_dp*(1.0_dp + sx*xi)*(1.0_dp + sy*eta)
    dxi  = 0.25_dp*sx*(1.0_dp + sy*eta)
    deta = 0.25_dp*sy*(1.0_dp + sx*xi)
  end subroutine shapefunc

end program cubature_demo