Heterogeneous list

#1

It’s been mentioned in another thread: What is a pure function? - #23 by everythingfunctional that heterogeneous lists would be helpful. Here is an example from @everythingfunctional: Brad Richardson / heterogeneous_list · GitLab.

Say you have a finite element mesh with different kinds of elements (hex, tet, prism, pyramid), and you need to store them somehow. One way to do it is to just have homegenous lists (arrays) for each element separately, so 4 different arrays.

@gnikit you also mentioned “An extremely useful feature to have if you want to experiment with different finite element discretisations in unstructured geometries”. Can you give an example?

4 Likes
What is a pure function?
#2

I believe this is an example of that concept:

1 Like
#3

If your PDE has tensor products, like in radiation transport (3D in space X 2D in angle in its simplest form) you have to discretise both the spatial and angular dimensions. You can think of the angular dimensions as the direction of propagation of particles at each point in space using a spherical coordinate system (azimuthal and polar angles being the 2D in angle).

If you choose to discretise the angular domain with multiple methods (there are good reasons why you would want that) e.g. finite elements and spherical harmonics and haar wavelets, storing the flux and boundary condition data in heterogeneous lists would be more convenient IMO than having an array for each discretisation method. The problem gets exacerbated if you choose to adapt the angular meshes which would result in the Degrees of Freedom vastly varying from spatial mesh node to spatial mesh node.

1 Like
#4

A heterogeneous list could be used to model DataFrames. A bit unrelated, yes, but still on topic I suppose.

5 Likes
#5

Is there a silver bullet? At some point you arrive at indirect addresing using arrays of integers.
If you don’t want to manage the complexity yourself you write/use a finite-element form compiler (e.g. FreeFEM, FEniCSx, FireDrake, libCEED).

The classic approach you see in Fortran codes is to have essentially a linked-list, with the pointers stored in a separate offset array:

Example 
program main
implicit none

integer, parameter :: TRI03 = 5 !> Linear triangle
integer, parameter :: QUA04 = 9 !> Linear quad

integer :: npoin, nelem
integer, allocatable :: lconn(:), lcoff(:), ltype(:)
real, allocatable :: x(:), y(:)

npoin = 6
nelem = 3

! 2-----3-----5
! | v   |     |
! |   v |     |
! 0-----1-----4  

x = [0.0, 1.0, 0.0, 1.0, 2.0, 2.0]
y = [0.0, 0.0, 1.0, 1.0, 0.0, 1.0]

! Connectivity
lconn = [0, 1, 2, 1, 3, 2, 1, 4, 5, 3]  !> Connectivity
lcoff = [3, 6, 10]                      !> Offset
ltype = [TRI03, TRI03, QUA04]           !> Element type

call dumpvtk(npoin,x,y,nelem,lconn,lcoff,ltype)

contains

subroutine dumpvtk(npoin,x,y,nelem,lconn,lcoff,ltype)
   use, intrinsic :: iso_fortran_env, only: int8
   integer, intent(in) :: npoin
   real, intent(in) :: x(:), y(:)
   integer, intent(in) :: nelem
   integer, intent(in) :: lconn(:), lcoff(:), ltype(:)

   character(len=512) :: buffer
   character(len=1), parameter :: nl = new_line('a')
   integer :: ivtu, i

   open(newunit=ivtu,file="dump.vtu")

   write(ivtu,'(A)') '<?xml version="1.0"?>'
   write(ivtu,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'

   write(ivtu,'(A)') '<UnstructuredGrid>'

   write(buffer,"('<Piece NumberOfPoints=',A,' NumberOfCells=',A,'>')") dq(npoin), dq(nelem)
   write(ivtu,'(A)') trim(buffer)
   
   write(ivtu,'(A)') '<Points>'
   write(ivtu,'(A)') '<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
   write(ivtu,'(*(E12.5,1X))') [(x(i),y(i),0.0, i=1,npoin)]
   write(ivtu,'(A)') '</DataArray>'
   write(ivtu,'(A)') '</Points>'

   write(ivtu,'(A)') '<Cells>'
   write(ivtu,'(A)') '<DataArray type="Int32" Name="connectivity" format="ascii">'
   write(ivtu,'(*(I0,1X))') lconn
   write(ivtu,'(A)') '</DataArray>'
   write(ivtu,'(A)') '<DataArray type="Int32" Name="offsets" format="ascii">'
   write(ivtu,'(*(I0,1X))') lcoff
   write(ivtu,'(A)') '</DataArray>'
   write(ivtu,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
   write(ivtu,'(*(I0,1X))') int(ltype,int8)
   write(ivtu,'(A)') '</DataArray>'
   write(ivtu,'(A)') '</Cells>'

   write(ivtu,'(A)') '</Piece>'
   write(ivtu,'(A)') '</UnstructuredGrid>'
   write(ivtu,'(A)') '</VTKFile>'
   
   close(ivtu)

end subroutine

function dq(i)
   integer, intent(in) :: i
   character(len=:), allocatable :: dq
   dq = '"'//itoa(i)//'"'
end function

function itoa(i)
   integer, intent(in) :: i
   character(len=:), allocatable :: itoa
   character(len=32) :: tmp
   write(tmp,'(I0)') i
   itoa = trim(tmp)
end function

end program
2 Likes
#6

It’s useful for neural networks as well. For example, in neural-fortran, the user instantiates a network from an array of layer instances:

Though the array of layers is not exactly heterogeneous, heterogeneity is emulated with each layer having concrete implementation as a polymorphic component.

5 Likes
#7

The question is not whether heterogeneous lists are helpful, for the Fortran language via its standard has long supported them (albeit with verbose syntax and a long, long, long waiting period for compiler implementations).

The question is how such a need translates to a compelling use case for PURE subroutines with polymorphic received arguments with the INTENT(OUT) attribute:

module m
   type :: t
   end type
contains
   subroutine sub1( a ) !<-- Ok
      class(t), intent(out) :: a
   end subroutine 
   pure subroutine sub2( a )  !<-- Ok
      type(t), intent(out) :: a
   end subroutine 
   pure subroutine sub3( a ) !<-- not allowed by the standard
      class(t), intent(out) :: a
   end subroutine 
end module 

C:\temp>gfortran -c m.f90
m.f90:8:26:

    8 |    pure subroutine sub3( a )
      |                          1
Error: INTENT(OUT) argument 'a' of pure procedure 'sub3' at (1) may not be polymorphic
2 Likes
#8

Thanks everybody for the examples.

I think Fortran already supports a homogeneous list of the base polymorphic class, so that allows to emulate a heterogeneous list. Is there anything else that needs to be done?

Regarding the pure issue, I think it’s because of the finalizer, correct?

#9

@gnikit, as I have been trying to convey, it is accepted the usefulness and benefits of having heterogenous lists in Fortran code, particularly in “setting up” computations (the so-called preprocessing layers) and consuming and presenting the results (post-processing) from all the number-crunching. Staring Fortran 2003 thru’ the 2008, 2018, 2023 revisions, the language does support such lists. The question is different as noted upthread.

1 Like