Example Integration of a New Element into the MODULEF Library

Version of 14 July 1994 for MODULEF 94

DOUGLAS N. ARNOLD and REVATHI NARASIMHAN

Introduction. This brochure illustrates in detail the integration of a new element into the MODULEF library. The goal is to help the user test new elements with a minimal amount of programming. We have chosen a basic piecewise linear element, so that the steps involved can be clearly illustrated in as simple a possible setting. We assume that the reader is familiar with MODULEF at the level covered in “Example Implementation of a Finite Element Computation Using the MODULEF Library,” which we shall henceforth refer to as “Example Implementation.”

1. Statement of the Problem. Let Ω be a polygonal domain in ℝ^{2}. We consider the following boundary value problem:

Find u : Ω → ℝ such that

(1)

−∇ · a∇u = f in Ω,
a∂u/∂n = g_{N} on Γ_{N}, u = g_{D} on Γ_{D},

where the coefficient a may depend on x ∈ Ω, Γ_{N} is the union of certain edges of ∂Ω and Γ_{D} is the union of the remaining edges. The corresponding weak formulation is:

Find u ∈ H^{1}(Ω) satisfying u = g_{D} on Γ_{D} and
∫_{Ω} a∇u · ∇v = ∫_{Ω} fv + ∫_{Γ_{N}} g_{N}v for all v ∈ H^{1}(Ω) vanishing on Γ_{D}.

We shall use MODULEF to solve this problem using continuous piecewise linear finite elements. Let S_{h} denote the space of all continuous piecewise linear functions subordinate to a given triangulation of Ω. The Galerkin approximation of our problem is

Find u_{h} ∈ S_{h} satisfying u_{h} = g_{D}^{I} on Γ_{D} and
∫_{Ω} a∇u_{h} · ∇v = ∫_{Ω} vf + ∫_{Γ_{N}} vg_{N} for all v ∈ S_{h} vanishing on Γ_{D},

where g_{D}^{I} denotes the piecewise linear interpolant of g_{D} on Γ_{D}. To complete the discretization we must specify quadrature rules for the three integrals entering the Galerkin approximation. For a triangle T with barycenter z_{0} and vertices z_{1}, z_{2}, and z_{3}, we define

∫_{T}^{[1]} w = area(T) · w(z_{0}), ∫_{T}^{[3]} w = area(T) · [w(z_{1}) + w(z_{2}) + w(z_{3})]/3.

Typeset by AMS-TEX

For a line segment e with endpoints z_{0} and z_{1} we set

∫_{e}^{[2]} w = length(e) · [w(z_{1}) + w(z_{2})]/2.

The fully discretized problem, the one we shall solve on the computer, is

Find u_{h} ∈ S_{h} satisfying u = g_{D}^{I} on Γ_{D} and
(2) ∑_{T} ∫_{T}^{[1]} a∇u_{h} · ∇v = ∑_{T} ∫_{T}^{[3]} vf + ∑_{e} ∫_{e}^{[2]} vg_{N}  for all v ∈ S_{h} vanishing on Γ_{D},

where the sum on T is over all triangles in the triangulation and the sum on e is over all triangle edges contained in Γ_{N}.

The solution of (2) with MODULEF can be broken into the six steps:

(1) Generation of the triangulation

(2) Description of the finite element space

(3) Generation of the elementary stiffness matrices

(4) Generation of the elementary load vector

(5) Assembly of the global stiffness matrix and load vector, and imposition of boundary conditions

(6) Solution of the resulting linear system There are MODULEF routines which perform steps 1, 5, and 6 in great generality. In fact, the weak formulation, finite element space, and quadrature rules we consider are treated by the element TRIA 2P1D in the thermal library of MODULEF, so we can use existing MODULEF modules to perform all six steps. How to do this was treated in the note “Example Implementation.” In this brochure we shall ignore the existence of the TRIA 2P1D element and explain how to write the routines to perform steps 3, 4, and 5, and how to interface them with existing MODULEF routines. MODULEF booklet #95, “Integration d’un Nouvel Element Fini,” treats this same subject, in greater generality and detail.

To incorporate an element into MODULEF the element should belong to an element library. There are presently two main libraries of finite elements in MODULEF: the thermal library THER and the elasticity library ELAS. (There are some smaller libraries as well.) Elements in the former library mostly pertain to scalar second order elliptic problems while those in the latter mostly pertain to elasticity problems. This is not a hard and fast rule, however, and is not in any way built into the software. In fact an element of any kind could be added to either library.^{1} We shall add our element to the thermal library. Each element in a library has a name which must consist of two words of four characters. We shall name our element TRIA MYP1.^{2}

1There is one major distinction between the thermal library and the elasticity library that is built into the software: the former usually uses single precision arithmetic and the latter double precision.

2This name doesn’t quite conform to MODULEF’s conventions since the conventional name for this element, TRIA 2P1D, is already taken.

2. Description of Interpolation. The mesh can be generated exactly as in “Example Implementation.” This results in a NOPO database describing a mesh of triangles with the boundary edges assigned reference numbers 1 and 2 according to whether they are contained in Γ_{D} or Γ_{N}. Once the mesh has been generated, we need to define our element interpolation. More precisely, we must describe the points, nodes, and degrees of freedom of our element on a reference triangle. This information will be incorporated into a MAIL database. The definition of the element interpolation is most commonly accomplished by calling the module COMAC2, which creates the mail database, passing to it the names of the element library (THER) and of the element (TRIA MYP1). COMAC2 then calls the routine DOBTHE (for elements in the thermal library; it calls DOBELA for elements in the elasticity library). DOBTHE in turn call a definition routine for the particular element TRIA MYP1. The definition routine is conventionally named TDMYP1 with “T” indicating the thermal library, “D” indicating a definition routine, and the remainder being the second word of the element name. Since our element is new, we shall have to supply the definition routine. The calling sequence for the element defintion routine is summarized in figure 1.

COMAC2

(MODULEF routine, creates MAIL database)

↓

DOBTHE

(MODULEF routine, calls definition routine, minor modifications needed for new element)

↓

TDMYP1

(Element definition routine, user-supplied)

Fig. 1. Calling sequence for the element definition routine.

The definition routine merely assigns values which encode the element interpolation to a pair of small arrays and calls a MODULEF utility routine (DOMAIL) to write them to appropriate locations in memory. It thus consists of only a few executable statements. There is a interactive MODULEF program, DEFIXX, which queries you for the information and writes the routine for you, thereby decreasing the chance of indexing errors.

The first array, MD, is integer. The main information encoded in it are the nodes and points of the element. Recall that a node is a point of the geometrical element associated with one or more degrees of freedom for an unknown of the problem. In our simple case there is only one unknown (the “temperature” u) and, since we are using continuous piecewise linear interpolation, its degrees of freedom are its values at the three vertices. Thus there are three nodes, and each is of the same type: it is associated to just the one unknown, temperature, and just one degree of freedom for that unknown, its value. The points of the element (in the MODULEF sense) determine its mapping onto a reference element; in fact they may be viewed as nodes for the mapping. In our case the mapping is just the usual affine mapping, which is determined by its values at the vertices. So the points coincide with the nodes. For a quadratic finite element, however, there would be six nodes but still just three points, so the points wouldn’t coincide with the nodes. For a curved finite element there might be more points, or the points might support degrees of freedom of the mapping to the reference element other that just its value (e.g., tangential derivatives). The entries of the array MD, in order, and their values for our element are as follows.

• length of this array minus 1 (determined by counting once the rest of the table has been completed)

28

• name of the library (character data to be converted to integer)

THER

• first word of element name (character data to be converted to integer)

TRIA

• second word of element name (character data to be converted to integer)

MYP1

• code number for the element (for a new element choose the next available number, thermal elements start at 100001)

100049

• code for element geometry (1=point, 2=segment, 3=triangle, 4=quadrilateral, …)

3

• number of nodes  3
• number of points  3
• number of associated arrays (usually none)  0
• dimension of the space  2
• number of unknowns  1
• number of types of nodes  1
• number of types of points  1

• 0 if the only degree of freedom of the reference map is its value at the points, 1 else

0

• 1 if the nodes and points coincide, 0 else  1

• for each unknown a description (character string of 4 words converted to integer)

TEMPERATURE

• for each type of node, the number of unknowns at that node

1

• for each type of node and for each unknown at that node, the number of the unknown and the number of degrees of freedom associated to it 1,1 • for each type of node, for each unknown at that

node, and for each degree of freedom associated with it, the code describing it (1=VN: value of the unknown, 2=DX: x derivative of the unknown, …; the list is long and can be extended)

1

• if points support additional degrees of freedom, they would come here

• for each node, its type  1,1,1
• for each point, its type  1,1,1

• if there were associated arrays, they would come here The necessary conversions of character data to integer can be accomplished by the

MODULEF utitility routines ICHAR4 and TRCHTA.

The second array, XD, simply gives the coordinates of the nodes on the reference element, followed by the coordinates of the points on the reference element. As for MD (and many arrays in MODULEF) the first element is an integer stating how many words follow. Thus the routine TDMYP1A is as follows.

                             subroutine tdmyp1(m)
                    c  definition routine for TRIA MYP1
                          dimension md(29),xd(13),xdp(12),m(*)
                          equivalence (xd(2),xdp(1)),(xd(1),ld)
                          data md/28,0,0,0,100049,3,3,3,0,2,1,1,1,0,1,0,
                         &         0,0,0,1,1,1,1,1,1,1,1,1,1/
                          data ld,xdp/12,0.,0.,1.,0.,0.,1.,0.,0.,1.,0.,0.,1./
                    c  put the character entries into md using Modulef utility routines
                          md(2) = ichar4('THER')
                          md(3) = ichar4('TRIA')
                          md(4) = ichar4('MYP1')
                          call trchta('TEMPERATURE',md(16),4)
                    c  now store the information
                          call domail(m,m,md,xd)
                          return
                          end

In order that DOBTHE can call TDMYP1, we must, of course, modify DOBTHE. There is a CHARACTER*8 array in DOBTHE called NB with the names of all the elements in the thermal library. Its dimension, a parameter named LNB, is the number of elements in the library. At the time of this writing this was 48. To add our element, we increase LNB to 49 and define NB(49) to be “TRIAMYP1.” The rest of DOBTHE mainly consists of a sequence of computed GOTO statements. We must add a branch to TDMYP1 as the destination point of the last of these. A listing of the relevant portion of DOBTHE follows. The changes made to incorporate MYP1 are underlined.

                          SUBROUTINE DOBTHE(M,I1,I2,NCOD,NRECH)

⋮

                          PARAMETER (LNB=\underline{49})

⋮

                          CHARACTER*8 TYPELE,NB(LNB)

                          DATA NB /
                         +   'TRIA2P1D','TRIA2P2D','TRIA2P2C','TRIA2P3D','TRIA2P3C',
                         +   'QUAD2Q1D','QUAD2Q2D','QUAD2Q2C','QUAD2Q3D','QUAD2Q3C',
                         +   'TETR3P1D','TETR3P2D','TETR3P2C','TETR3P3D','TETR3P3C',
                         +   'HEXA3Q1D','HEXA3Q2D','HEXA3Q2C','HEXA3Q3D','HEXA3Q3C',
                         +   'PENT3R1D','PENT3R2D','PENT3R2C','PENT3R3D','TETRM3T1',
                         +   'TRIAHT10','TRIAHT21','TRIAHT32','QUADHQ10','QUADHQ21',
                         +   '$$$$$$$$','TRI2THD2','TRI2TEQ2','TRIAMT10','QUADMQ10',
                         +   'TRIAMT21','QUADMQ21','QUADWACH','TRIAAP1D','QUADAQ1C',
                         +   'HEXAM3H1','TRIA2DCD','RESERVE ','RESERVE ','TRIAAP2C',
                         +   'QUADAQ2C','QUAD2I1D','QUADAI1D'\underline{,'TRIAMYP1'}/

⋮

                          GO TO (410,420,1111,1111,450,460,470,480\underline{,490}),II4
                    C
                      410 CALL TDM3H1(M)
                          RETURN

⋮

                      470 CALL TD2I1D(M)
                          RETURN
                      480 CALL TDAI1D(M)
                          RETURN
                      490 CALL TDMYP1(M)
                          \underline{RETURN}
                          END

3. Generation of the Elementary Stiffness Matrix. The elementary stiffness matrices are stored in a TAE database by the module THENEW (which replaces the older module THELAS). In the argument list for THENEW is an argument ONRAID, which gives an external routine. This routine should call the the elementary stiffness routine for the particular element. If the MODULEF program TNRAID is specified, the standard MODULEF thermal elements are used. In this case we supply our own routine ONRAID, which calls our own elementary stiffness routine TRMYP1 (note that “RAID” and the “R” in “TRMYP1” come from the French word raideur or stiffness). The main input parameters to the elementary stiffness routine give the coordinates of the vertices of the element, the number of the subdomain containing the element, and information about the coefficients of the differential equation on the element (exactly what information depends on the element; typically given are the values of the coefficients at the integration points). The main output parameter is an array, AE, containing the elementary stiffness matrix entries.

To compute the elementary stiffness matrix, the values of the coefficients of the differential equation at the integration points of the element must be known. In our case, this means that the value of a at the element barycenter is needed. There are three different ways we could have the user furnish this information. The first two possibilities are in a table (stored in the super-array), or via a subroutine, called MILIEV. In either case, THENEW will read this information for the element in question, and then pass it on to the elementary stiffness routine as an input parameter named SUR. The third possibility is that the elementary stiffness routine call a user supplied subroutine directly. The database MILI, created by MODULEF routine COMILI according to user supplied data cards, determines which means of evaluating the coefficients is used, and this information is passed to the elementary stiffness routine in the input parameter NOTELM. We shall use only the second possibility: the coefficients will be given in subroutine MILIEV.

Actually, the elementary stiffness routine is called twice for each element. On the first call, AE is not computed, but some information about it is returned: NMMAST which gives the length of the array AE which will be output, NSTOC which gives the option for its storage, NTYTA which gives the type of the array, NINDI which gives the number of indices of the array, and LSUR which gives the length of the input array SUR containing the coefficient data. In our case, the elementary stiffness matrix is a 3 × 3 symmetric matrix, which we shall store as a real vector of length 6. That is the elementary stiffness matrix is given as

AE(1) AE(2) AE(4)
AE(2) AE(3) AE(5)
AE(4) AE(5) AE(6)

Thus NMMAST=6, NSTOC=1 (symmetric matrix storage), NTYTA=2 (real values), NINDI=1, and LSUR=1 (since only the value of a single coefficient at the barycenter is to be passed in SUR). On the second call to the elementary stiffness routine, NMMAST and NSTOC must be returned again, as well as the array AE itself.

The routine TRMYP1 and the analogous routine TSMYP1 for computing the elementary load vectors share a common block labelled THELA0. This allows calculations common to both routines to be performed only once. Since, it could happen that in some circumstances one or the other of the elementary stiffness matrix or the elementary load vector is not computed, flags are used to indicate when the values in the common THELA0 have already been computed. By convention the flags are the first four variables in the common. These are named INDICV, INDICS, INDICL, and INDICP and are used to flag volumetric, superficial, linear, and point information respectively. In our case only INDICS is needed.

A listing of the elementary stiffness routine TRMYP1 follows.

                          subroutine trmyp1(nocal,notelm,nmmast,ntyta,nindi,nstoc,ndsde,
                         +                  x,y,sur,lsur,ae)
                    c**********************************************************************
                    c  This routine calculates the elementary stiffness matrix for the
                    c  TRIA MYP1 element.
                    c
                    c     Input parameters:
                    c     -----------------
                    c       nocal  - set to 1 on the first call, 2 on the second
                    c       notelm - options for coefficient evaluation
                    c       x      - array containing the x-coordinates of the nodes
                    c                of the triangle.
                    c       y      - array containing the y-coordinates of the nodes
                    c                of the triangle.
                    c       sur   - array containing coefficient information
                    c       ndsde - number of the subdomain containing the element
                    c

c     Output parameters:
c     ------------------
c       nmmast - length in words of the output array ae
c       nstoc -  type of storage used in the output array ae
c    nocal = 1:
c    ----------
c       ntyta  - type of the output array ae
c       nindi  - number of indices in the output array ae
c       lsur   - length in words of the input array sur
c    nocal = 2:
c    ----------
c       ae     - the elementary stiffness matrix
c
c**********************************************************************
      dimension ae(*),x(*),y(*),noref(*),sur(*),aret(*),notelm(*)
c****
c     Create a labeled common block so that information about
c     the current element can be shared among related routines,
c     (e.g., once the elementary stiffness matrix is computed for
c     an element, the quantities x21, y21, ..., delta can be reused
c     for the computation of the elementary load vector).
c****
      common/thela0/indicv,indics,indicl,indicp,x21,y21,x31,y31,x32,y32,
     +              delta,nww(1013)
c
      nstoc   = 1
      nmmast  = 6
      if ( nocal .le. 1 ) then
         ntyta   = 2
         nindi   = 1
         lsur    = 1
         return
      end if
c
c****
c     indics is a flag indicating whether the variables in common
c     /thela0/ have been initialized (indics=1) or not (indics=0).
c****
      if ( indics .le. 0 ) then
c****
c     Compute and store basic geometrical information about element.
c****
         x21 = x(2)-x(1)
         y21 = y(2)-y(1)
         x31 = x(3)-x(1)
         y31 = y(3)-y(1)
         x32 = x(3)-x(2)
         y32 = y(3)-y(2)
         delta = x21*y31 - x31*y21
c****
c     set indics to flag that the information has been computed
c****
         indics = 1
      end if
      alf1 = sur(1)/(2.*delta)
c****
c     calculation of ae
c****
      ae(1) = ( x32*x32 + y32*y32)*alf1

ae(2) = (-y31*y32 - x31*x32)*alf1
ae(3) = ( y31*y31 + x31*x31)*alf1
ae(4) = ( y21*y32 + x21*x32)*alf1
ae(5) = (-y21*y31 - x21*x31)*alf1
ae(6) = ( y21*y21 + x21*x21)*alf1
end

Next we must supply the routine ONRAID which will call TDMYP1. One way to do this is to copy and modify the MODULEF routine TNRAID. The changes are analogous to those we made to DOBTHE to have it call TDMYP1. In fact, this time its even easier: we just have to add one more branch to a computed GOTO statement.

      SUBROUTINE TNRAID(M,NDIM,NOCAL,NOTELM,NOTELF,NMMAST,NTYTA,NINDI,

⋮

GO TO (1001,1002,1003,1004,9999,1006,1007,1008,9999,9999,
+       1011,9999,1013,9999,9999,1016,1017,1018,9999,9999,
+       1021,9999,9999,9999,1025,1026,9999,9999,9999,9999,
+       9999,1032,1033,1034,1035,1036,1037,1038,1039,1040,
+       1041,1042,9999,9999,1045,1046,1047,1048\underline{,1049}),NCOD1

⋮

1049 CALL TRMYP1(NOCAL,NOTELM,NMMAST,NTYTA,NINDI,NSTOC,NDSDE,
+            X(1,1),X(1,2),SURFAC,LSURFA,AE)
\underline{GO TO 1000}
C
9999 WRITE (IMPRIM,99991) NCODE
STOP
C
1000 CONTINUE
IF (NOCAL .EQ. 2 .AND. IPAQUE .NE. 1) THEN
   CALL TYPMOT(NTYTA,NBREMO)
   INDICE(2,1) = NMMAST/NBREMO
   INDICE(2,4) = 1
END IF
END

This approach will give us access to all the existing MODULEF elements as well as our new one. (This would be useful, for example, if we mixed old and new elements on the same mesh.) However if we only want to use our own element, we can supply a trivial program in place of TNRAID. Here we have supplied the following simple program

subroutine onraid(m,ndim,nocal,notelm,notelf,nmmast,ntyta,nindi,
+          nstoc,x,ndi,npo,
+          volume,lvolum,surfac,lsurfa,arete,larete,point,lpoint,
+          volumf,lvoluf,surfaf,lsurff,aretf,laretf,pointf,lpoinf,
+          noref,iaxelo,naxelo,axeloc,lcpt,ae,ncode,nfpoba,ndsde,
+          nono,axenod,lxnod,lvect,trav,ltrav,indice)
c**************************************************************
c    This is a copy of MODULEF's program TNRAID, with calls to all
c    the existing elements removed, and replaced by a call to TRMYP1.
c    To incorporate the element into the MODULEF library, this
c    routine would not be used, and instead TNRAID would be modified
c    accordingly.
c**************************************************************
dimension m(*),ae(*),volume(*),surfac(*),arete(*), point(*),

     +           x(lvect*npo,*),volumf(*),surfaf(*),aretf(*),
     +           pointf(*),trav(*),indice(4,4)
      integer noref(*),notelm(*),notelf(*)
      common/paquet/ipaque
c
      call trmyp1(nocal,notelm,nmmast,ntyta,nindi,nstoc,ndsde,
     +            x(1,1),x(1,2),surfac,lsurfa,ae)
      if (nocal .eq. 2 .and. ipaque .ne. 1) then
         call typmot(ntyta,nbremo)
         indice(2,1) = nmmast/nbremo
         indice(2,4) = 1
      end if
      end

4. Generation of the Elementary Load Vectors. The treatment of the elementary load vectors is entirely analogous to that of the elementary stiffness matrices. Again these are stored in a TAE database by THENEW. This time the argument ONSECM (analogous to ONRAID discussed above) to THENEW specifies the routine which calls the load vector routine for the particular element. The load vector routine should be called TSMYP1, with the “S” standing for second membre or right hand side.

The following routine generates the elementary load vector for MYP1. We should say something more here before giving the listing.

       subroutine tsmyp1(nocal,notelf,nmmast,ntyta,nindi,nstoc,
     +                   ndsde,x,y,sur,lsur,aret,laret,noref,ndsm,be)
c*********************************************************************
c  This routine calculates the elementary right side for the
c  TRIA MYP1 element.
C
c   Input parameters:
c   -----------------
c
c       nocal - number of the call to this routine per
c               element: assumes the value of 1 or 2.
c       notelf - array of option parameters.
c       x - array containing the x-coordinates of the nodes
c           of the triangle.
c       y - array containing the y-coordinates of the nodes
c           of the triangle.
c       sur - array containing the forcing function information.
c       aret - array containing information on the traction on
c              neumann portion on the boundary.
c       ndsde - number of the subdomain.
c       noref - array containing the reference numbers of the faces.
c               only faces with nonzero reference numbers contribute
c               to calculation of boundary term.
c       ndsm - number of right hand sides.
c
c   Output parameters:
c   ------------------
c
c       nmmast - length of the array be (in words).
c       nstoc - type of storage of the el. r.h.s.;
c               -1 for non-symmetric.
c

c       nocal = 1:
c       ----------
c       ntyta - type of array;  2 for real single precision.
c       nindi - number of indices in be.
c       lsur - length of the sur array (in words).
c       laret - length of the aret array (in words).
c
c       nocal = 2:
c       ----------
c       be - the elementary right hand side.
c
C*********************************************************************
      dimension be(ndsm,*),x(*),y(*),sur(ndsm,*),aret(ndsm,*)
      integer noref(*),notelf(*)
c****
c     Create a labeled common block so that information about
c     the current element can be shared among related routines,
c     (e.g., once the elementary stiffness matrix is computed for
c     an element, the quantities x21, y21, ..., delta can be reused
c     for the computation of the elementary load vector).
c****
      common/thela0/indicv,indics,indicl,indicp,x21,y21,x31,y31,x32,y32,
     +              delta,nww(1013)
c
      nstoc   = -1
      nmmast  = 3*ndsm
      if ( nocal .le. 1 ) then
         ntyta   = 2
         nindi   = 2
         laret   =  ndsm
         lsur    = 3*ndsm
         return
      end if
c
c****
c     indics is a flag indicating whether the variables in common
c     /thela0/ have been initialized (indics=1) or not (indics=0).
c****
      if ( indics .le. 0 ) then
c****
c     Compute and store basic geometrical information about element.
c****
         x21 = x(2)-x(1)
         y21 = y(2)-y(1)
         x31 = x(3)-x(1)
         y31 = y(3)-y(1)
         x32 = x(3)-x(2)
         y32 = y(3)-y(2)
         delta = (x21*y31 - x31*y21)
         indics = 1
      end if
c****
c    Compute integral over triangle using three vertex values.
c****
      do 6 n=1,ndsm
         do 5 i=1,3
            be(n,i) = sur(n,i) * delta / 6.
    5    continue
    6 continue

c****
c    For each of the three sides of the triangle, check whether
c    that side lies on the boundary (has positive reference number),
c    and compute the contribution to the load vector on that edge
c    using the two endpoint values.
c****
      if( noref(1) .gt. 0 ) then
         do 8 n=1,ndsm
            temp   = sqrt( x21*x21 + y21*y21 ) / 2. * aret(n,1)
            be(n,1) = be(n,1) + temp
            be(n,2) = be(n,2) + temp
    8    continue
      end if
c
      if( noref(2) .gt. 0 ) then
         do 11 n=1,ndsm
            temp   = sqrt( x32*x32 + y32*y32 ) / 2. * aret(n,2)
            be(n,2) = be(n,2) + temp
            be(n,3) = be(n,3) + temp
   11    continue
      end if
c
      if( noref(3) .gt. 0 ) then
         do 14 n=1,ndsm
            temp   = sqrt( x31*x31 + y31*y31 ) / 2. * aret(n,3)
            be(n,1) = be(n,1) + temp
            be(n,3) = be(n,3) + temp
   14    continue
      end if
      end

The program ONSECM which we supply is completely analogous to ONRAID (and one could also make a more elaborate version based on MODULEF’s TNSECM, as discussed above).

      subroutine onsecm(m,ndim,nocal,notelf,notelm,nmmast,ntyta,nindi,
     +          nstoc,x,ndi,npo,
     +          volume,lvolum,surfac,lsurfa,arete,larete,point,lpoint,
     +          volumm,lvolmm,surfam,lsurfm,aretm,laretm,pointm,lpoinm,
     +          noref,iaxelo,naxelo,axeloc,lcpt,ndsm,
     +          be,ncode,nfpoba,ndsde,nthela,dudt,nono,axenod,lxnod,
     +          lvect,trav,ltrav,indice)
c**************************************************************
c    This is a copy of MODULEF's program TNSECM, with calls to all
c    the existing elements removed, and replaced by a call to TSMYP1.
c    To incorporate the element into the MODULEF library, this
c    routine would not be used, and instead TNSECM would be modified
c    accordingly.
c**************************************************************
      dimension m(*),be(*),volume(*),surfac(*),arete(*), point(*),
     +          x(npo*lvect,*), volumm(*),surfam(*),aretm(*),
     +          pointm(*),dudt(*),trav(ltrav),indice(4,4)
      integer noref(*),notelm(*),notelf(*)
      common/paquet/ipaque
c
 1049 call tsmyp1(nocal,notelf,nmmast,ntyta,nindi,nstoc,ndsde,
     +     x(1,1),x(1,2),surfac,lsurfa,arete,larete,noref,ndsm,be)

                          if( nocal .eq. 2 ) nmmast =nmmast / ndsm
                          if (nocal .eq. 2 .and. ipaque .ne. 1) then
                             call typmot(ntyta,nbremo)
                             indice(3,1) = ndsm
                             indice(3,2) = nmmast/nbremo
                             indice(3,4) = 2
                          end if
                          end

We have not furnished routines to compute the elementary mass matrices or flux vectors, although a complete element for the MODULEF library should have these. Routines TMMYP1 and TCMYP1 are left as an exercise to the reader.

Using the new element. In “Example Implementation” we carried out in full the solution of a particular boundary value problem of the form of (1) using MODULEF’s 2P1D thermal element. This was accomplished by three programs: MESH, which generated the mesh and interpolation, STFLD which generated the stiffness matrix and load vector, and SOLVE which solved the algebraic equations. Assuming that they are linked with the modified version of DOBTHE and the new routines TDMYP1, TRMYP1, TSMYP1, ONRAID, and ONSECM, the same three programs can be used with the new element with only very slight changes. The only change to MESH is in the one line which refers specifically to the element by NAME: the name 2P1D must of course be changed to MYP1. The only change needed in STFLD is the external routines TNRAID and TNSECM must be replace by ONRAID and ONSECM. No changes are needed in SOLVE. The data files for MESH, STFLD, and SOLVE can be taken to be identical with those used for the 2P1D computation.