Example Implementation of a Finite Element Computation Using the MODULEF Library

Version of 14 July 1994 for MODULEF 94

DOUGLAS N. ARNOLD and REVATHI NARASIMHAN

Introduction. This brochure describes in detail the implementation of a simple finite element computation using MODULEF. The goal is to help a beginning user of MODULEF learn to write programs calling existing modules. The computation is a very simple one, namely the solution of Laplace’s equation on a square with Dirichlet conditions on three sides and Neumann conditions on the fourth using piecewise linear elements. We have chosen a problem whose exact solution is a linear function, so that it is easy to verify that MODULEF has computed it correctly. The programs discussed here were written for illustrative purposes and are thoroughly commented, making it relatively easy for a user to modify the programs for other computations.

1. Statement of Problem. Let Ω be the unit square [0, 1] × [0, 1]. Set Γ_{N} = {1} × (0, 1), the right edge of Ω, and Γ_{D} = ∂Ω \ Γ_{N}. The boundary value problem we solve is

−Δu = 0 on Ω, u = −2 + 5x − 4y on Γ_{D}, \frac{∂u}{∂n} = 5 on Γ_{N},

summarized in figure 1.

(1,1)

u = −2 + 5x − 4y →  −Δu = 0  ← \frac{∂u}{∂n} = 5 (0,0)

Fig. 1

The weak or variational formulation of this problem is:

Find u ∈ H^{1}(Ω) satisfying u = −2 + 5x − 4y on Γ_{D} and
(1)  ∫_{Ω} ∇_{∼}u · ∇_{∼}v dx = 5 ∫_{Γ_{N}} v ds  for all v ∈ H^{1}(Ω) vanishing on Γ_{D}.

Clearly, the exact solution is u = −2 + 5x − 4y.

Typeset by AMS-TEX

We shall use MODULEF to solve this problem using the simplest finite elements, continuous piecewise linears. Thus the discrete problem we wish to solve is

Find u_{h} ∈ S_{h} satisfying u = −2 + 5x − 4y on Γ_{D} and

∫_{Ω} ∇_{∼}u_{h} · ∇_{∼}v dx = 5 ∫_{Γ_{N}} v ds for all v ∈ S_{h} vanishing on Γ_{D},

where S_{h} is the space of all continuous piecewise linear functions subordinate to some triangulation of Ω.

2. Implementation through MODULEF. The computation of the discrete problem can be broken into a series of steps for each of which there is an existing module in the MODULEF library. Thus our code will consist mainly of a sequence of calls to MODULEF subroutines.

We have chosen to break the computation into three parts, each corresponding to a separate main program:

• program MESH: mesh generation and interpolation;

• program STFLD: calculation of the stiffness matrix and load vector;

• program SOLVE: solution of the linear system. Each part calls several MODULEF modules. The division into three programs is somewhat arbitrary. For example, we could as well have computed the stiffness matrix and load vector in separate programs (although this is slightly less efficient), or have included all steps in one program.

Following the MODULEF philosophy, the communication between the three programs will be entirely through certain standardized databases which will reside on disk files. Specifically, the MESH program outputs a MAIL database and a COOR database. (The former describes the topology of the mesh and the finite element interpolation, the latter the coordinates of the mesh points.) These databases are read by STFLD, which writes a MUA and a B database (containing the assembled stiffness matrix and load vector, respectively). Finally, these are input to SOLVE which writes a second B database that stores the solution vector. This flow is summarized in figure 2. In the course of the computation, various intermediate databases are created, but not saved to permanent disk files. For example, during the course of program STFLD a TAE database is created to store the elementary stiffness matrices.

⟶^{MESH} MAIL, COOR ⟶^{ST F LD} MUA, B ⟶^{SOLV E} B

Fig. 2

Each of the programs takes some card input. MODULEF’s free format data facility allows us to comment the data files, which we have done extensively. Program STFLD also requires that we supply subprograms to describe the coefficients, forcing function, and boundary data.

MODULEF includes a library of finite elements. (It is possible to use MODULEF without using one of the elements in its library, but this involves more coding.) In particular, the present element, involving the variational formulation (1) and continuous, piecewise linear elements, is included in the element “TRIA2P1D” of the thermal library. This element will be invoked when we call the module COMAC2 from the mesh program.

Program Listings

The remainder of this brochure consists of listings of the programs and datafiles. It is highly recommended that the the reader compile, link, and run these programs. It is also advisable to run the MODULEF main program IMAGXX to print out in human readable format the contents of the MAIL, COOR, MUA, and two B databases created.

3. Mesh generation and interpolation. This is accomplished by the program MESH:

                          program MESH
                    c*****************************************************************
                    c  Program MESH:  Generate the mesh and describe the element
                    c                 interpolation
                    c
                    c  Modules accessed:
                    c
                    c      Name     Purpose
                    c      ----     -------
                    c      INITI    initialize
                    c      APNOPO   generate geometric mesh
                    c      COMAC2   describe interpolation
                    c
                    c  Card input:  data cards for INITI and APNOPO
                    c
                    c  User supplied functions:  none
                    c
                    c  Input databases:  none
                    c
                    c  Output databases:
                    c
                    c              File
                    c      Type    name     Description
                    c      ----    ----     -----------
                    c      MAIL    mail     mesh with interpolation data from COMAC2
                    c      COOR    coor     coordinates of points of mesh from COMAC2
                    c
                    c      (if desired, the intermediate NOPO database created by
                    c      APNOPO can be saved via a "SAUV" command in the card
                    c      input for APNOPO)
                    c
                    c*****************************************************************

      parameter (lm=100000,impre=6,nnn=0,maxsbd=10)
c****
c     lm:     size of super array
c     impre: printing level during initialization
c     nnn: printing level by address utilities
c     maxsbd: max # of subdomains
c****
      dimension m(lm),icodel(4*maxsbd),nolc(1)
      double precision dm
      external ffront
      character*4  nombib
      equivalence (m(1),dm)
      common m
c
c****
c     Call INITI to initialize and read title (alternatively, one
c     could use INITIS to initialize without assigning a title)
c****
      call initi(m,lm,impre,nnn)
c
c****
c     Open files and assign units for MAIL and COOR
c       We use Modulef's trunit subroutine to assign a free
c       unit number, and Modulef's ouvrir subroutine, which
c       is a more portable replacement for the Fortran open
c       statement, to open the file
c****
      call trunit(nfmail)
      call ouvrir(nfmail,'mail','UNFORMATTED,UNKNOWN',0,iostat)
      if (iostat.ne.0) go to 5000
      nimail = 1
      call trunit(nfcoor)
      call ouvrir(nfcoor,'coor','UNFORMATTED,UNKNOWN',0,iostat)
      if (iostat.ne.0) go to 5000
      nicoor=1
c
c****
c     Call module APNOPO to generate the NOPO database (geometric
c     mesh)
c
c     INPUT:
c     ------
c     m:       super array
c     ffront:  routine to describe curved portions of the boundary
c              (not used here)
c                (Alternatively, an interpreted function can
c                be read from data cards to describe the a curved
c                portion of the boundary.)
c
c     OUTPUT:
c     -------
c     NOPO database, kept in core.
c****
      call apnopo(m,m,ffront)
c
c****
c     Call module COMAC2 to specify element interpolation
c
c     INPUT:

c     ------
c     m:      super array
c     ndsd:   number of subdomains
c     nblc:   number of curved segments
c     nolc:   array of reference #'s of curved segments (not
c             used if nblc = 0 )
c     nombib: name of the Modulef library containing the element
c             to use
c     icodel: array describing the element to use (for each subdomain
c             this array contains the number of different types of
c             straight elements and for each its name, followed
c             by the number of different types of curved elements
c             and for each its name
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c     ntXXXX: number of user supplied arrays to associate
c             with the XXXX database (usually none)
c
c     OUTPUT:
c     -------
c     MAIL and COOR databases.
c****
      ndsd = 1
      nblc = 0
      nombib = 'THER'
      do 20 i = 1,ndsd
         icodel(4*(i-1)+1) = 1
         icodel(4*(i-1)+2) = ichar4('TRIA')
         icodel(4*(i-1)+3) = ichar4('2P1D')
         icodel(4*(i-1)+4) = 0
20    continue
      nfnopo = 0
      ninopo = 2
c****
c     Level 2 nopo is set in the data cards to APNOPO
c****
      ntmail = 0
      ntcoor = 0
      call comac2(m,ndsd,nblc,nolc,nombib,
     &            icodel,nfnopo,ninopo,nfmail,nimail,
     &            ntmail,nfcoor,nicoor,ntcoor)
c
      stop
c
5000  print*,' error in opening file: iostat =',iostat
      stop
c
      end

The data for the program MESH is as follows.

$------------------------------------------------------------
$ Data card for INITI: title
$------------------------------------------------------------
'EXAMPLE, step 1: program MESH'
$------------------------------------------------------------
$ Data cards for APNOPO
$------------------------------------------------------------

$  [Remark:  for the unit square it is simpler to triangulate
$  using the keyword OBJE, but in the interest of generality
$  we use define the domain via points and lines and use TRIA.]
$  Enter points
'POIN'
$  printing   # of points
$  IMPRE        NPOINT
    2             4
$  For each point, give
$  # of pt.  ref #     x     y
$  NOP(I)  NREF(I)    X(I)  Y(I)
    1        1         0.    0.
    2        1         1.    0.
    3        1         1.    1.
    4        1         0.    1.
$  Enter lines
'LIGN'
$  printing   # of lines
    2            4
$  For each line, give
$     line #  # nodes   #end 1   #end 2   ref. # line opt.  ratio
$      NOLIG   NOELIG   NEXTR1   NEXTR2   NOREFL  NFFRON   RAISON
         1        3        1        2        1       0       1.0
         2        3        2        3        2       0       1.0
         3        3        3        4        1       0       1.0
         4        3        4        1        1       0       1.0
$   Triangulate
'TRIA'
$ printing   level   subdomain#   # of lines  sense
$ IMPRE     NIVEAU     NUDSD        NBRELI     NPROPA
    2          1         1            4         1
$ Lines describing (sub)domain, in order
    1     2     3   4
$  components of bdry   internal pts to add    generation option
   1                    0                      1
$ number of lines in each component of boundary
     4
$ Renumber the nodes and elements to decrease matrix band width
'RENE'
$ printing   level in   level out
$ IMPRE        NIVO1      NIVO2
    2            1          2
$ Done
'FIN'

4. Calculation of the stiffness matrix. This is accomplished by the program STFLD:

      program  STFLD
c************************************************************
c  Program STFLD:  Construct the stiffness matrix and the
c                  load vector
c
c  Modules accessed:
c

c      Name     Purpose
c      ----     -------
c      INITI    initialize
c      COMILI   describe storage of coefficient information
c      COFORC   describe storage of forcing function and
c               Neumann data
c      THENEW   create elementary stiffness matrix and
c               elementary load vector
c      CORNOE   compute coordinates of nodes
c      COBDC1   describe Dirichlet boundary conditions
c      PREPAC   prepare pointers for assembly
c      ASSMUA   assemble stiffness matrix
c      ASEMBV   assemble load vector
c      CLIMPC   incorporate Dirichlet boundary conditions into
c               the stiffness matrix and load vector
c
c  Card input: data cards for INITI, COMILI, and COFORC
c
c  User supplied functions:
c
c      Name     Purpose
c      ----     ------
c      MILIEV   describe coefficients
c      FORCEV   describe forcing function and Neumann data
c      VALCLR   describe Dirichlet data
c
c  Input databases:
c
c              File
c      Type    name     Description
c      ----    ----     -----------
c      MAIL    mail     mesh with interpolation data from MESH
c      COOR    coor     coordinates of points of mesh from MESH
c
c  Output databases:
c
c              File
c      Type    name     Description
c      ----    ----     -----------
c      MUA     mua      assembled stiffness matrix from CLIMPC
c      B        b       assembled load vector from CLIMPC
c
c**************************************************************
      parameter (lm=100000,iprint=6,ipradr=0)
      dimension  m(lm),iopt(4),nofr(1,3),val(1),ita5(1,3),ta6(1,2)
      double precision dm,d1(1),d2(1)
      character*4 nomcou
      logical fonint
      external valclr, valcld, tnmass, tnraid, tnsecm, tncntr,
     &         forcev, miliev
      equivalence (m(1),dm)
      common m
c****
c     Call INITI to initialize
c****
      call initi(m,lm,iprint,ipradr)
c
c****
c     Assign units and open files for input MAIL and COOR databases
c     and output MUA and B databases.  Also open a scratch file for

c     TAE database.
c****
      call trunit(nfmail)
      call ouvrir(nfmail,'mail','UNFORMATTED,OLD',0,iostat)
      if (iostat.ne.0) go to 5000
      nimail = 1
      call trunit(nfcoor)
      call ouvrir(nfcoor,'coor','UNFORMATTED,OLD',0,iostat)
      if (iostat.ne.0) go to 5000
      nicoor = 1
      call trunit(nfmua2)
      call ouvrir(nfmua2,'mua','UNFORMATTED,UNKNOWN',0,iostat)
      if (iostat.ne.0) go to 5000
      nimua2 = 1
      call trunit(nfb2)
      call ouvrir(nfb2,'b','UNFORMATTED,UNKNOWN',0,iostat)
      if (iostat.ne.0) go to 5000
      nib2 = 1
      call trunit(nftae)
      call ouvrir(nftae,' ','UNFORMATTED,SCRATCH',0,iostat)
      if (iostat.ne.0) go to 5000
      nitae = 1
c
c****
c     Call module COMILI to create database MILI, which
c     will be used along with subroutine MILIEV to
c     describe the coefficient function for module
c     THENEW.
c
c     INPUT:
c     ------
c     m     : super array
c     nfmili: unit number of the mili database (0=in core)
c     nimili: level number of the mili database
c     ntmili: number of user supplied arrays to the database
c
c     OUTPUT:
c     -------
c     the database MILI
c****
      nfmili = 0
      nimili = 1
      ntmili = 0
      call comili(m,nfmili,nimili,ntmili)
c
c****
c     Call module COFORC to create database FORC, which
c     will be used along with subroutine FORCEV to
c     describe the forcing function for module
c     THENEW.
c
c     INPUT:
c     ------
c     m     : super array
c     nfforc: unit number of the forc database (0=in core)
c     niforc: level number of the forc database
c     ntforc: number of user supplied arrays to the database
c
c     OUTPUT:
c     -------

c     the database FORC
c****
      nfforc = 0
      niforc = 1
      ntforc = 0
      call coforc(m,nfforc,niforc,ntforc)
c
c****
c     Call module THENEW to create the TAE database containing the
c     elementary stiffness matrix and elementary load vectors
c
c     INPUT:
c     ------
c     m:      super array
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c     ieXXXX: set to 1 if the database XXXX is used; 0 if not
c     iopt:   array containing  options for calculating the mass
c             matrix, stiffness matrix, load vector, and fluxes
c     nl:     number of linear combinations of the mass and stiffness
c             matrices
c     d1,d2:  arrays containing the the coefficients of the i'th
c             linear combination; used only when nl is nonzero
c     nntae:  name of array to be associated with TAE.  Used when
c             nl is nonzero
c     nprov:  indicates type of problem: 1 for Thermal, 2 for

c             Elasticity, ...

c     nfpoba: direct access file containing precomputed information
c             for certain quadrature rules (not used here)
c     nthela: option for thermo-elastic problems
c     nfbo,nibo:     databases containing information on initial
c                    deformations (not used here)
c     nfndl1,nindl1: database used if all the nodes do not

c                    have the same degrees of freedom

c     lvect:  the number of elements to treat in this call (lvect > 1
c             is used for vectorization)
c     nomcou: the name of an array to store the element "colors",
c             used only if lvect > 1
c     thmass: name of subprogram to call the elementary mass matrix routine
c     thraid: name of subprogram to call the elementary stiffness matrix
c             routine
c     thload: name of subprogram to call the elementary load vector routine
c     thcntr: name of subprogram to call the elementary stress tensor routine
c     forcev:  name of subroutine giving loading data
c     miliev: name of subroutine giving material coefficients
c
c     OUTPUT:
c     -------
c     the database TAE
c****
      nttae = 0
      iemili = 1
      ieforc = 1
      iopt(1) = 0
      iopt(2) = 1
      iopt(3) = 1
      iopt(4) = 0
      nl = 0
      nprov = 1

      nfpoba = 0
      nthela = 0
      nfbo = 0
      nibo = 0
      nfndl1 = 0
      nindl1 = 0
      lvect = 1
      nomcou = 'NONE'
      call thenew(m,m,m,nfmail,nimail,nfcoor,nicoor,
     &            iemili,nfmili,nimili,ieforc,nfforc,
     &            niforc,nftae,nitae,nttae,iopt,nl,d1,d2,
     &            nntae,nprov,nfpoba,nthela,nfbo,nibo,
     &            nfndl1,nindl1,lvect,nomcou,tnmass,tnraid,
     &            tnsecm,tncntr,forcev,miliev)
c
c****
c     The next step is to create a BDCL database describing the
c     Dirichlet boundary conditions using the modulef COBDC1.
c     First we call module CORNOE to make a second COOR database, this
c     one containing the coordinates of the nodes (rather than the points).
c     Such a database is required by modulef COBDC1.  Actually this step
c     is unnecessary for P1 elements, since nodes and points coincide,
c     but we include it for completeness.
c****
      nfcoor2 = 0
      nicoor2 = 2
      ntcoor2 = 0
      call cornoe(m,nfmail,nimail,nfcoor,nicoor,nfcoor2,nicoor2,
     &            ntcoor2)
c
c****
c     Call module COBDC1 which creates the BDCL database containing
c     information on the Dirichlet data
c
c
c     INPUT:
c     ------
c     m:      super array
c     iconst: indicates if the b.v. of the d.o.f. is given
c             by an array of constants (iconst=1) or by a
c             function (iconst=2)
c     nbfr:   number of references numbers for which degrees
c             of freedom are to be imposed
c     nofr:   array of dimension (nbfr,3) describing the
c             degrees of freedom to impose:
c             for i = 1 to nbfr
c                 nofr(i,1) = reference number
c                 nofr(i,2) = number of the unknown
c                 nofr(i,3) = mnemonic of the d.o.f.
c             NOTE: nofr must be declared as nofr(nbfr,3)
c     nty:    type of the values to impose (2: real)
c     val:    array of type nty used to give the values to be
c                imposed if iconst=1
c     nclrl:  number of linear combination boundary conditions
c                to impose
c     nmtrl:  number of terms in the linear combination (irrelevant
c                if nclrl=0)
c     ita5,ta6: arrays which describe the linear combinations
c               (irrelevant here since nclrl=0)

c     ieXXXX: set to 1 if the database XXXX already exists; 0 if not
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c             Note:  the TYNO and MPLR databases will be created
c             since they don't exist (ietyno=iemplr=0).  The
c             NDL1 database is not needed since all nodes have
c             the same degrees of freedom.
c     nopfi : option number of storing the database MAIL
c                0: all in central memory
c                1: only one element is in central memory
c     fonint: b.c. are to specified by an "interpreted function",
c             i.e., an expression to be read as data and parsed
c     valclr: function subprogram giving boundary conditions if
c             fonint = .false., single precision (nty = 2)
c     valcld: function subprogram giving boundary conditions if
c             fonint = .false., double precision
c
c     OUTPUT:
c     -------
c     the BDCL database
c****
      iconst = 2
      nbfr = 1
      nofr(1,1) = 1
      nofr(1,2) = 1
      nofr(1,3) = ichar4('VN  ')
      nty  = 2
      nclrl = 0
      nmtrl = 0
      iemail = 1
      nftyno = 0
      nityno = 1
      ietyno = 0
      nfmplr = 0
      nimplr = 0
      iemplr = 0
      nfbdcl = 0
      nibdcl = 1
      ntbdcl = 0
      nopfi = 0
      fonint = .false.
      call cobdc1(m,iconst,nbfr,nofr,nty,val,nclrl,nmtrl,
     &            ita5,ta6,iemail,nfmail,nimail,ietyno,
     &            nftyno,nityno,iemplr,nfmplr,nimplr,nfndl1,
     &            nindl1,nfcoor2,nicoor2,nfbdcl,nibdcl,ntbdcl,
     &            nopfi,fonint,valclr,valcld)
c
c****
c     Call module PREPAC which prepares for the assembly of the
c     global stiffness matrix and global load vector
c
c     INPUT:
c     ------
c     m:      super array
c     nnomax: maximum number of nodes in an element
c     nd:     the number of d.o.f. at each node if same throughout.
c             0, if not; (see CONDL1 for the latter case)
c     ndsm:   number of right hand sides
c     nclrl:  number of linear combination boundary conditions

c             to impose
c     ncodst: code for the storage of the matrix A: -1, non-symmetric;

c             0, symmetric; 1, symmetric

c     nty:    type of the arrays A and B; 2, real single precision;
c             5, real double precision;
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c
c     OUTPUT:
c     -------
c     lbdp:   bandwidth of MUA
c     bdpdl:  largest separation between d.o.f. of
c             the same element
c     mcndl1: words of memory occupied by
c             the NDL1 database
c     mcmua:  words of memory occupied by
c             the MUA database
c     mcb:    words of memory occupied by
c             the B database
c****
      nnomax = 3
      nd = 1
      ndsm = 1
      nclrl = 0
      ncodst = 1
      nty = 2
      nfmua = 0
      nimua = 1
      ntmua = 0
c
      call prepac(m,nnomax,nd,ndsm,nclrl,ncodst,nty,
     &            nfmail,nimail,nfndl1,nindl1,nfbdcl,nibdcl,
     &            nfmua,nimua,ntmua,
     &            lbdp,lbdpdl,mcndl1,mcmua,mcb)
c
c****
c     Call module ASSMUA to assemble the global stiffness matrix
c
c     INPUT:
c     ------
c     not:    number of the array in the TAE to assemble
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c
c     OUTPUT:
c     -------
c     the assembled MUA database
c
c     NOTE:  the same MUA database is used both as input and output.
c     Consequently it is overwritten.
c****
      not = 1
      call assmua(m,not,nftae,nitae,nfndl1,nindl1,nfmua,nimua,
     &            nfmua,nimua)
c
c****
c     Call module ASEMBV to assemble the global load vector
c
c     INPUT:
c     ------

c     m:      super array
c     not:    number of the array in the TAE to assemble
c     ndsm:   number of load vectors
c     ncodsb: always = -1
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c
c     OUTPUT:
c     -------
c     the assembled B database
c****
      not = 2
      ndsm = 1
      ncodsb = -1
      nfb = 0
      nib = 1
      ntb = 0
      call asembv(m,not,ndsm,ncodsb,nftae,nitae,nfndl1,
     &            nindl1,nfb,nib,ntb)
c
c****
c     Call module CLIMPC to incorporate the Dirichlet boundary
c     conditions into the stiffness matrix and load vector
c
c     INPUT:
c     ------
c     m:      super array
c     nivo:   option for how to process
c     vtg:    parameter for approximate imposition of b.c.s
c     nfXXXX: file number of the XXXX database (0 if kept in core)
c     niXXXX: level of the XXXX database
c
c     OUTPUT:
c     -------
c     the MUA and B databases after imposition of boundary conditions.
c****
      nivo = 3
      vtg = 1.0
      call climpc(m,nivo,vtg,nfndl1,nindl1,nfmua,nimua,nfb,nib,
     &            nfbdcl,nibdcl,nfmua2,nimua2,nfb2,nib2)
      stop
c
5000  print*,'error in opening file: iostat =',iostat
      stop
c
      end

In addition we need to provide three subprograms:

      subroutine miliev(m,lopt,x,ndim,npo,tar,ltar,lvect,iadr,
     &                  i1,nare,ia)
c**************************************************************
c
c     This subroutine supplies the coefficient information
c     for the THENEW module for the computation of the
c     elementary stiffness matrix.
c
c     Input Parameters:
c     -----------------

c     M : the super array
c     LOPT : parameter relating the number of the subdomain
c            with its assigned information in the data cards
c            for COMILI
c     X : array of dimension NPO by NDIM, containing the
c         coordinates of the vertices of the current
c         triangle
c     LTAR : length of the array TAR (in words) needed for the
c            current problem; calculated by THENEW
c     LVECT : number of elements to treat on this call
c     IADR : beginning index in TAR; given in the datacards
c            for COMILI
c     I1 : this parameter is zero for standard uses; it
c          is a translation parameter for some internal arrays
c          in THENEW
c     NARE : number of the current face (3-D) or side (2-D)
c            under consideration
c     IA : address-1 in the D.B. MAIL which can be used to
c          retrieve information about the current element.
c
c     Output Parameter:
c     -----------------
c
c     TAR : contains coefficient information for the computation
c           of the elementary stiffness matrix.  For the
c           thermal problem, isotropic case, TAR(1) contains
c           the value of the coefficient function at the
c           barycenter.  TAR(2)=TAR(3)=0 since we're considering
c           the isotropic case.
c
c**********************************************************************
      dimension  m(*),tar(*),x(npo,ndim)
98    format ('*** FATAL ERROR in subroutine MILIEV ',
     &        I2,' is an invalid value for LOPT.')
c
      if (lopt .eq. 1) then
         tar(1) = 1.0
         tar(2) = 0.0
         tar(3) = 0.0
      else
         write(6,98) lopt
         call tilt
      endif
c
      return
      end
      subroutine forcev(m,lopt,x,ndim,npo,tar,ltar,lvect,iadr,
     &                 il,nare,ia)
c
c***********************************************************
c
c     This subroutine supplies information about the forcing
c     function and the Neumann boundary data, which is used
c     in the computation of the elementary load vector.
c
c     Input Parameters:
c     -----------------
c     M : the super array
c     LOPT : parameter relating the number of the subdomain

c            or the reference number of part of boundary with its
c            assigned information in the data cards for COFORC
c     X : array of dimension NPO by NDIM, containing the
c         coordinates of the vertices of the current
c         triangle
c     LTAR : length of the array TAR (in words) needed for the
c            current problem; calculated by THENEW
c     LVECT : number of elements to treat on this call
c     IADR : beginning index in TAR; given in the datacards
c            for COFORC
c     I1 : this parameter is zero for standard uses; it
c          is a translation parameter for some internal arrays
c          in THENEW
c     NARE : number of the current face (3-D) or side (2-D)
c            under consideration
c     IA : address -1 in the D.B. MAIL which can be used to
c          retrieve information about the current element.
c
c     Output Parameters:
c     ------------------
c
c     The output must be consistent with the integration rule
c     for the TRIA 2P1D element of the MODULEF thermal library.
c
c     TAR: array which passes back information about the forcing
c          function or Neumann boundary data on the current element
c          as follows:
c
c          TAR(I)  :  value of forcing function at node I; I=1,2,3,
c                          --OR--
c                     value of Neumann data at endpoint I of the
c                     side; I=1,2.
c
c******************************************************************
      dimension m(*),tar(*),x(npo,ndim)
c
98    format ('*** FATAL ERROR in subroutine FORCEV ',
     &        I2,' is an invalid value for LOPT ')
99    format ('*** FATAL ERROR in subroutine FORCEV ',
     &        I2,' is an invalid value for NARE ')
c
      if (lopt .eq. 1) then
c**** Forcing function is identically zero
      do 10 j=1,3
        tar(j)=0.0
10    continue
      else
         if (lopt .eq. 2) then
c****
c        Neumann data, in this case identically = 5
c
c        The following computed GOTO statement is unnecessary
c        for this case as g is constant.  However, this structure
c        allows different formulas for the different sides of the
c        triangle.
c****
            goto(11,12,13) nare
            write(6,99) nare
            call tilt

11          tar(1) = 5.
            tar(2) = 5.
            goto 100
12          tar(1) = 5.
            tar(2) = 5.
            goto 100
13          tar(1) = 5.
            tar(2) = 5.
         else
            write(6,98) lopt
            call tilt
         endif
      endif
100   return
      end
      real function valclr(i,x,y,z)
c****
c    Function describing the Dirichlet boundary conditions
c    for the sample problem. This routine is needed
c    by the subroutine COBDC1.
c
c     Input parameters:
c     -----------------
c     I     : ranges from 1 to the number of reference numbers where
c             Dirichlet boundary conditions are imposed
c             NOTE: I is NOT necessarily the actual reference number.
c     X,Y,Z : coordinates of the point under consideration.
c
c     The output value is the value to impose.
c***
c
98    format ('*** FATAL ERROR in function VALCLR ',
     &        I2,' is an invalid value for I.')
c
      goto (10) i
      write(6,98) i
      call tilt
10    valclr = -2. + 5.*x - 4.*y
      return
      end

The data for the program STFLD is as follows.

$------------------------------------------------------------
$ Data card for INITI: title
$------------------------------------------------------------
'EXAMPLE, step 2: program STFLD'
$------------------------------------------------------------
$ Data cards for the subroutine COMILI
$------------------------------------------------------------
$ mil2(i); i=1,8
$  The MILI database may contain or not each of the eight
$  tables, MIL3, MIL4, ..., MILA.  The first four are used
$  to give volumetric, superficial, linear, and punctual
$  coefficient information by reference number/subdomain,
$  respectively; the next four are used to give the
$  the same four kinds of information by element number.
$  The eight values mil2(i), i=1,8, flag the existence

$  of these eight tables.  A zero value indicates
$  non-existence, while a positive value gives how many
$  reference or subdomain numbers or elements are described.
$  Here we give only 2-dimensional information by subdomain
$  number.
  0 1 0 0 0 0 0 0
$
$ Loop 1 to mil2(2), giving the number of the subdomain
$ and other information pertinent to  the description
$ of the coefficients
$
$ NUDSD: number of the subdomain
$ NFRO: only for 3-D elements; assigned a value
$       of 0 for 2-D elements
$ NOPT: 1 for isotropic; cf. documentation on MILI D.S.
$ ITRAIT: option of treatment of the coefficient function:
$         1: by array
$         2: by subroutine MILIEV
$         3: by user's own subroutine
$ LOPT: parameter to pass to subr. MILIEV corresponding
$       to the subdomain under consideration
$ IADR: parameter to pass to subr. MILIEV; corresponds to
$       to the first index of the TAR array in that
$       subprogram (usually set to 1)
$
  1 0 1 2       $ NUDSD   NFRO  NOPT  ITRAIT
  1 1           $ LOPT    IADR
$-------------------------------------------------------------
$ Data cards for COFORC
$-------------------------------------------------------------
  1    $ NDSM : number of load vectors
$
$ for2(i); i=1,8
$  Describes which of the eight tables FOR3, ..., FORA is
$  present and how many reference numbers, subdomain
$  numbers, or elements they describe, as for MILI above.
$  Here we describe 2D information for 1 subdomain, and
$  1D information for 1 reference number.
  0 1 1 0 0 0 0 0
$
$ Loop 1 to for2(2), giving the ref. # of the subdomain
$ and other information pertinent to the forcing function:
$
$ NUDSD: number of the subdomain
$ NFRO: only for 3-D elements; assigned a value
$       of 0 for 2-D elements
$ NOPT: 1 for isotropic
$ ITRAIT: option of treatment of the forcing function:
$         1: by array
$         2: by subroutine FORCEV
$         3: by user's own subroutine
$ LOPT: parameter to pass to subr. FORCEV corresponding
$       to the subdomain under consideration
$ IADR: parameter to pass to subr. FORCEV; corresponds to
$       to the first index of the TAR array in that
$       subprogram (usually set to 1)
$

  1 0 1 2       $ NUDSD  NFRO  NOPT  ITRAIT
  1 1           $ LOPT   IADR
$
$ for2(3): number of references numbers used for Neumann
$          data
$
$ Loop 1 to for2(3), giving the ref. # of the side
$ and other information pertinent to Neumann data:
$
$ NOSD : number of the side.
$ LOPT : parameter to pass to subr. FORCEV corresponding
$        to the side under consideration
$
$ The other parameters carry the same definitions as above.
$
  2 0 1 2       $ NOSD  NFRO  NOPT  ITRAIT
  2 1           $ LOPT  IADR

5. Solution of the linear system. This is accomplished by the program SOLVE:

      program SOLVE
c*************************************************************
c  Program SOLVE : Solve the resulting linear system
c
c  Modules accessed:
c
c      Name     Purpose
c      ----     -------
c      CHOLPC   compute Cholesky decomposition of the matrix
c      DRCHPC   solve linear system using forward and
c               backward substitution
c
c  Card Input: Data cards for INITI
c
c  User Supplied Functions: none
c
c  Input databases:
c
c              File
c      Type    name     Description
c      ----    ----     -----------
c      MUA     mua      assembled stiffness matrix from STFLD
c      B        b       assembled load vector from STFLD
c
c  Output database:
c
c              File
c      Type    name     Description
c      ----    ----     -----------
c      B       soln     the solution vector
c
c*************************************************************
c
      parameter(lm=100000,iprint=6,ipradr=0)
      double precision dm
      equivalence (m(1),dm)
      common m(lm)
c
c****

c     Initialize and read title
c****
      call initi(m,lm,iprint,ipradr)
c
c****
c     Assign units and open files for input MUA and B databases,
c     and output B database for the solution
c****
      call trunit(nfmua)
      call ouvrir(nfmua,'mua','UNFORMATTED,OLD',0,iostat)
      if (iostat.ne.0) go to 5000
      nimua = 1
      call trunit(nfb)
      call ouvrir(nfb,'b','UNFORMATTED,OLD',0,iostat)
      if (iostat.ne.0) go to 5000
      nib = 1
      call trunit(nfsoln)
      call ouvrir(nfsoln,'soln','UNFORMATTED,UNKNOWN',0,iostat)
      if (iostat.ne.0) go to 5000
      nisoln = 1
c
c****
c     Call CHOLPC to compute the Cholesky decomposition
c     of the global stiffness matrix
c
c     INPUT:
c     ------
c     eps     : machine precision for single precision
c                  numbers (our example is in single precision)
c     nentre  : parameter set to zero if the subroutine is
c               to return to its called location at the presence
c               of L  such that L   <= eps* A
c                   ii           ii          ii
c     nfxxxx  : file number of the xxxx database
c     nixxxx  : level of the xxxx database
c               (nfndl1=0 since there is no ndl1 database here)
c
c     OUTPUT:
c     -------
c     the factored global matrix, MUA2: resides in core ;
c     the parameter nretou: =1 if a pivot less than eps*A
c                                                        ii
c****
      eps = 1.e-6
      nentre = 0
      nfndl1 = 0
      nindl1 = 0
      nfmua2  =  0
      nimua2  =  1
      call cholpc(m,eps,nentre,nfndl1,nindl1,nfmua,nimua,
     &            nfmua2,nimua2,nretou)
      if (nretou .ne. 0) goto 6000
c
c****
c     Call DRCHPC to solve the system using fwd. and
c     bwd. substitution
c
c     INPUT:
c     ------
c     The factored matrix (MUA2) and the assembled right hand side (B)

c     m       : super array
c     niveau  : flag, set to 2 to perform complete solution of
c
c                             t
c                      [L] [L] {X} = {B}
c
c     nclrl   : number of linear relations in the Dirichlet
c               B.C.'s (none in this example)
c     nfxxxx  : file number of the xxxx database
c     nixxxx  : file number of the xxxx database
c               (nfndl1=0 since there is no ndl1 database here and
c                nfbdcl doesn't apply since nclrl=0)
c
c     OUTPUT:
c     -------
c     the solution vector SOLN
c****
      niveau=2
      nclrl=0
      call drchpc(m,niveau,nclrl,nfndl1,nindl1,nfmua2,nimua2,
     &            nfb,nib,nfbdcl,nibdcl,nfsoln,nisoln)
c
      stop
c
5000  print*,'error in opening file, iostat =',iostat
      stop
c
6000  print*,'*** FATAL ERROR in solve:  found nearly vanishing pivot'
      stop
c
      end

The data for the program SOLVE is as follows.

$------------------------------------------------------------
$ Data card for INITI: title
$------------------------------------------------------------
'EXAMPLE, step 3: program SOLVE'