UMACR1.F - - METIS Example

The umacr[0-9].f files in the user directory allow one to create user macro commands. To set up a user macro command one needs to edit the file umacr1.f (for example) to include the code one wishes to have executed when the user mesh command is called from inside FEAP. The interface can be programed in fortran or it can be programed in C. In the example that follows a fortran interface is utilzed. The example shows how one can interface the popular graph partitioning library metis with FEAP.

The main steps in creating a user macro command whether in C or Fortran are

Set the name of the macro command to whatever 4-letter name you wish it to have when FEAP calls it.
Set up the source code you wish to have executed when the command is called from FEAP

The macro command name is set in the common variable uct in the common block umac1. When uct is equal to 'mac1' (for umacr1.f or .c) it should be set to whatever you wish for the name of your mesh command.

A example umacr1.f is given below. The name of the user macro is set to graph. The macro command itself is used to create element and node graphs from FEAP's mesh data and to partition the graphs using calls to the METIS library which can be downloaded from the METIS site. The command takes as options the keywords elem or node and the number of domains desired.

As these routines us the user memory allocation system. The needed ualloc file is also contained below.

Source file:

      subroutine umacr1(lct,ctl,prt)
c
c      * * F E A P * * A Finite Element Analysis Program
c
c.... Copyright (c) 1984-1994: Robert L. Taylor, Sanjay Govindjee
c
c     User defined routine to add command to FEAP
c
c     lct    - character field 'xxx' on 'ctrl,xxx,ctl(1),ctl(2),ctl(3)'
c     ctl(3) - real field on input command "   "    "      "      "
c     prt    - true if printing is active
c
c
c     Create element and node graphs and partition them
c      with the METIS libraries
c
c      Node    Graph partitioning not implemented
c      Element Graph partitioning     implemented
c
c------------------------------------------------------------------

implicit none

      include 'cdata.h'
      include 'sdata.h'
      include 'umac1.h'

      include 'pointer.h'
      include 'upointer.h'
      include 'comblk.h'

logical pcomp,prt, setval,ualloc, palloc
logical egraph, ngraph

      integer     mnsize,nodglen,elmglen,ixinvlen,edgecut, i
      integer     domains
      integer     options(5), moptions(3)

character*4 lct

real*8 ctl(3)

save

if(pcomp(uct,'mac1',4)) then

c     Set name to user defined
        uct = 'grap'
        egraph = .true.
        ngraph = .true.

setval = ualloc(1,'VAR1',numnp+1,1)
call size_ixinv(mr(np(33)),mr(up(1)),ixinvlen)

        setval = ualloc(2,'VAR2',ixinvlen,1)
        setval = palloc(111,'TEMP1',numnp,1)
        call create_ixinv(mr(np(33)),mr(up(1)),mr(np(111)),mr(up(2)))
        setval = palloc(111,'TEMP1',0,1)
        elmglen = 0

c Create a node graph option
elseif(pcomp(lct,'node',4)) then

        if(ngraph) then
          setval = ualloc(5,'XADJN',numnp+1,1)
          setval = ualloc(6,'NODG',1,1)

c                              IX         V         IXINV
c                              XADJN      NODG      LEN
          call make_node_graph(mr(np(33)),mr(up(1)),mr(up(2)),
     &                         mr(up(5)), mr(up(6)),nodglen)

          setval = ualloc(6,'NODG',nodglen,1)
          ngraph = .false.
        endif

c Create an element graph option
elseif(pcomp(lct,'elem',4)) then

        if(egraph) then
          setval = ualloc(3,'XADJE',numel+1,1)
          setval = ualloc(4,'ELMG',1,1)
c                              IX         V         IXINV
c                              XADJE      ELMG      LEN
          call make_elem_graph(mr(np(33)),mr(up(1)),mr(up(2)),
     &                         mr(up(3)), mr(up(4)),elmglen)

setval = ualloc(4,'ELMG',elmglen,1)

          options(1) = 0
          egraph = .false.
        endif

        if(up(7).eq.0) then
          setval = ualloc(7,'PARTE',numel,1)
        else
          call pzeroi(mr(up(7)),numel)
        endif

call pzeroi(moptions,3)

domains = nint(ctl(1))

c      Partition the element graph
        if(domains.le.8) then
          call Metis_PartGraphRecursive(numel,mr(up(3)),mr(up(4)),
     &                           moptions(1),moptions(2),moptions(3),1,
     &                           domains,options,edgecut,mr(up(7)))
        else
          call Metis_PartGraphKway(numel,mr(up(3)),mr(up(4)),
     &                           moptions(1),moptions(2),moptions(3),1,
     &                           domains,options,edgecut,mr(up(7)))
        endif

c
c       Stuff partition number into region number slot
c       for later use
c
        do i = 1, numel
         mr(np(33)+i*nen1-2) = mr(up(7)+i-1)
        end do

      else
        write(*,*) 'ENTER: GRAPh,NODE,#; or GRAPh,ELEMent,#'
      end if

end

c    Routine needed for determining the size of the inverse ix list and the needed pointers
      subroutine size_ixinv(ix,v,ixinvlen)
c
c     Input: ix       array of nodes on elements
c     Output: v        array of pointers to start of records for ixinv
c             ixinvlen length of ixinv array
c
c
      implicit none

include 'cdata.h'
include 'sdata.h'

integer e,j,n,ixinvlen
integer ix(nen1,numel),v(*)

      do e = 1, numel
        do j=1,nen
          n = ix(j,e)
          if (n.ne.0) v(n+1) = v(n+1) + 1
        end do
      end do

      ixinvlen = 0
      v(1) = 1
      do j = 2, numnp
        ixinvlen = ixinvlen + v(j)
        v(j) = v(j-1) + v(j)
      end do
      ixinvlen = ixinvlen + v(numnp+1)
      v(numnp+1) = ixinvlen + 1

end

subroutine create_ixinv(ix,v,skip,ixinv)

c Create ix-inverse in compressed form with a pointer list
c in v()

c Input: ix, v array
c Output: ixinv

implicit none

include 'cdata.h'
include 'sdata.h'

integer e,j,n
integer ix(nen1,numel),v(*),skip(*),ixinv(*)

      do e=1, numel
        do j=1, nen
          n = ix(j,e)
          if (n.ne.0) then
            ixinv(v(n) + skip(n)) = e
            skip(n) = skip(n) + 1
          endif
        end do ! j
      end do ! e

end

c Routine for creating an element graph
subroutine make_elem_graph(ix,v,ixinv,point,elemg,len)

c     Inputs: ix
c              v
c              ixinv
c
c     Outputs: elemg    -- element graph
c              point    -- pointer array for element graph
c              len      -- length of element graph array

implicit none

include 'cdata.h'
include 'sdata.h'

      logical insert
      integer jj,e,n,i,ce, len,skip,t
      integer ix(nen1,numel),v(*),ixinv(*),point(*),elemg(*)

point(1) = 1
c len = 0

do e = 1, numel
skip = 0

do n = 1, nen

jj = ix(n,e)

        if(jj.ne.0) then
         do i = v(jj), v(jj+1)-1
          ce = ixinv(i)
          insert = .true.

          do t = 0, skip
           if(ce.eq.elemg(point(e)+t) .or.
     &        ce.eq.e) insert=.false.
          end do ! t

          if(insert) then
            elemg(point(e)+skip) = ce
            skip = skip + 1
          endif

         end do ! i
        endif
       end do ! n

point(e+1) = point(e) + skip

end do ! e

len = point(numel+1) - 1

end

c Routine for creatining a node graph
subroutine make_node_graph(ix,v,ixinv,point,nodeg,len)

c     Inputs: ix
c              v
c              ixinv
c
c     Outputs: nodeg    -- node graph
c              point    -- pointer array for node graph
c              len      -- length of node graph array

implicit none

include 'cdata.h'
include 'sdata.h'

      logical insert
      integer jj,e,n,i,cn,k, len,skip,t
      integer ix(nen1,numel),v(*),ixinv(*),point(*),nodeg(*)

      point(1) = 1
      do n = 1, numnp
       skip = 0

do k = v(n), v(n+1)-1

e = ixinv(k)

        do i = 1, nen
         insert = .true.
         cn = ix(i,e)

          do t = 0, skip
           if(cn.eq.nodeg(point(n)+t) .or.
     &        cn.eq.n .or. cn.eq.0) insert=.false.
          end do ! t

          if(insert) then
            nodeg(point(n)+skip) = cn
            skip = skip + 1
          endif

end do ! i
end do ! k

point(n+1) = point(n) + skip

end do ! n

len = point(numnp+1) - 1

end

c User allocation routine needed for the METIS umacro interface

c$Id: ualloc.f,v 1.1 1997/07/29 19:15:14 rlt Exp $
logical function ualloc(num,name,length,precis)

c * * F E A P * * A Finite Element Analysis Program

c-----[--+---------+---------+---------+---------+---------+---------+-]
c      Purpose: Define, delete, or resize a user dictionary entry.
c               Pointer defined for integer (single) and real
c               (double precision) arrays.

c N.B. Currently limited to 200 names by dimension of
c common blocks 'allotd','allotn','pointer','upointer'

c Inputs:

c         num        - Entry number for array (see below)
c         name       - Name of array          (see below)
c         length     - Length of array defined: =0 for delete
c         precis     - Precision of array: 1 = integers; 2 = reals

c Output:

c up(num) - Pointer to first word of array in blank common
c-----[--+---------+---------+---------+---------+---------+---------+-]

implicit none

      logical          usetmem
      character*(*)    name
      integer          i, num,length,precis

include 'allotd.h'

c Storage definitions for UALLOC variables

integer list
parameter (list = 7)

character*5 names(list)

save

c Define and store names

data (names(i),i=1,list)/

& 'VAR1','VAR2','XADJE','ELMG','XADJN','NODG','PARTE'/

c Short description of variables

c              'VAR1',     !     1: Start here with user defined names
c              'VAR2',     !     2: Start here with user defined names
c              'XADJE',    !     3: Start here with user defined names
c              'ELMG',     !     4: Start here with user defined names
c              'XADJN',    !     5: Start here with user defined names
c              'NODG',     !     6: Start here with user defined names
c              'PARTE',    !     7: Start here with user defined names

c Do memory management operations

ualloc = usetmem(list,names,num,name,length,precis)

end