c This is an example 'main' program for calling MGGHAT
c
      program main
c
c In 'commons', parameter statements are used to set the dimension for
c the arrays.  Also, all of the program parameters are passed through it.
c
      include 'commons'
c
c solve the problem
c
      call mgghat
      stop
      end
c
c -------- pde
c
      subroutine pde(x,y,p,q,r,f)
      real x,y,p,q,r,f
c
c return the values of the pde coefficents at (x,y)
c pde is
c
c    -( p(x,y) * u  )  -( q(x,y) * u  ) + r(x,y) * u = f(x,y)
c                 x  x              y  y
c
c NOTE: BE CAREFUL TO GET THE SIGNS RIGHT
c e.g. p=q=1. means f=-(uxx+uyy)
c
      ipower=5
      p=1.
      q=1.
      r = 0.
      f=-ipower*(ipower-1.)*(x**(ipower-2)+y**(ipower-2))
c
      return
      end
c
c -------- bcond
c
      subroutine bcond(x,y,ipiece,c,g,itype)
      real x,y,c,g
      integer ipiece,itype
c
c returns boundary condition coefficients at (x,y)
c
c boundary condition is either
c
c    u  + c(x,y)*u = g(x,y)   or   u = g(x,y)
c     n
c
c (the natural b.c. is the Neumann b.c. because p=q=1)
c
c ipiece indicates the boundary piece from which the boundary condition
c is determined. itype must be set to the type of boundary condition for
c that piece (the type cannot change within a piece).
c itype = 1 Dirichlet (second condition above)
c         2 Neuman (first condition with c = 0)
c         3 Mixed (first condition with c != 0)
c
c (these pieces assume the assignment in the example inittr for
c  a rectangular domain)
c
      if (ipiece.eq.1) then
c left side; Dirichlet b.c.; U=g
         itype = 1
         c = 0.
         g = true(x,y)
      elseif (ipiece.eq.2) then
c bottom; Mixed b.c.; Un + U = g
         itype = 3
         c = 1.
         g = true(x,y) - truey(x,y)
      elseif (ipiece.eq.3) then
c right side; Neuman; Un = g
         itype = 2
         c = 0.
         g = truex(x,y)
      else
c top; dirichlet b.c.; U=g
         itype = 1
         c = 0.
         g = true(x,y)
      endif
c
      return
      end
c
c -------- true
c
      real function true(x,y)
      real x,y
      ipower=5
      true = x**ipower + y**ipower
      return
      end
c
c -------- truex
c
      real function truex(x,y)
      real x,y
      ipower = 5
      truex = ipower*x**(ipower-1)
      return
      end 
c
c -------- truey
c
      real function truey(x,y)
      real x,y
      ipower = 5
      truey = ipower*y**(ipower-1)
      return
      end
c
c -------- inittr
c
      subroutine inittr
      include 'commons'
c
c rectangular domain (ax,bx) X (ay,by) with an ngridx X ngridy grid
c
      ax = 0.
      bx = 1.
      ay = 0.
      by = 1.
      ngridx = 2
      ngridy = 2
c
      if (outlev.ge.3) write(ioutpt,101) ax,bx,ay,by
  101 format(' begin initializing triangulation'/
     .       ' triangulation for rectangle (',f8.2,
     .       ',',f8.2,') X (',f8.2,',',f8.2,')')
c
c set initial triangulation
c
c this is a user provided routine to define the domain and
c the initial triangulation
c
c the user provides:
c
c   nvert   - number of vertices
c   ntri    - number of triangles
c   xvert(1..nvert),yvert(1..nvert) - x and y coordinates
c                                     of the vertices
c   vertex(1..3,1..ntri) - the 3 vertices of each triangle
c   neigh(i,1..ntri) - which piece of the boundary contains the triangle
c                      side opposite vertex i. Need not be set if
c                      that side is not on the domain boundary.
c
c the peak of each triangle is the third vertex, i.e.,
c  peak(triangle)=vertex(3,triangle)
c
c assumptions on the triangulation are:
c
c 1) each vertex is in the same position (first, second, or third)
c    in every triangle that contains it, e.g., if vertex 2 is
c    the first vertex of triangle 3 (vertex(1,3)=2). vertex 2 must
c    also be the first vertex of any other triangle containing it.
c    thus we can have vertex(1,4)=2, and cannot have vertex(2,4)=2
c
c 2) vertices which are first or second vertices can be in at most
c    8 triangles.  vertices which are third vertices (peaks) can
c    be in at most 4 triangles
c
c This version of the routine triangulates the rectangle (ax,bx)X(ay,by).
c ngridx and ngridy specify the number of intervals in each dimension.
c The boundary pieces (ipiece for subroutine bcond) are:
c   1 - left
c   2 - bottom
c   3 - right
c   4 - top
c
c for the case of a 4X4 grid, the vertex and triangle numbers are:
c
c  5---10---15---20---25
c  |\   |   /|\   |   /|
c  | \ 8|15/ | \24|31/ |
c  | 7\ | /16|23\ | /32|
c  |   \|/   |   \|/   |
c  4----9---14---19---24
c  |   /|\   |   /|\   |
c  | 5/ | \14|21/ | \30|
c  | /6 |13\ | /22|29\ |
c  |/   |   \|/   |   \|
c  3----8---13---18---23
c  |\   |   /|\   |   /|
c  | \4 |11/ | \20|27/ |
c  | 3\ | /12|19\ | /28|
c  |   \|/   |   \|/   |
c  2----7---12---17---22
c  |   /|\   |   /|\   |
c  | 1/ | \10|17/ | \26|
c  | /2 | 9\ | /18|25\ |
c  |/   |   \|/   |   \|
c  1----6---11---16---21
c
      dy = by-ay
      dx = bx-ax
      dy = dy/ngridy
      dx = dx/ngridx
c
c set number of triangles and vertices
c
      nvert = (ngridx+1)*(ngridy+1)
      ntri  = 2*ngridx*ngridy
c
c set coordinates of vertices
c
      k=0
      do 11 i=1,ngridx+1
         xtemp=ax+(i-1)*dx
         if (i.eq.ngridx+1) xtemp=bx
         do 10 j=1,ngridy+1
            k=k+1
            xvert(k)=xtemp
            yvert(k)=ay+(j-1)*dy
            if (j.eq.ngridy+1) yvert(k)=by
   10    continue
   11 continue
c
c set vertices and boundary pieces of triangles
c
      do 21 i=1,ngridx,2
         do 20 j=1,ngridy,2
            ivbase=(i-1)*(ngridy+1)+j
            itbase=2*((i-1)*ngridy+j)-1
            vertex(1,itbase) = ivbase+ngridy+2
            vertex(2,itbase) = ivbase
            vertex(3,itbase) = ivbase+1
            if (i.eq.1) neigh(1,itbase) = -1
            if (j.eq.ngridy) neigh(2,itbase) = -4
            vertex(1,itbase+1) = ivbase+ngridy+2
            vertex(2,itbase+1) = ivbase
            vertex(3,itbase+1) = ivbase+ngridy+1
            if (j.eq.1) neigh(1,itbase+1) = -2
            if (i.eq.ngridx) neigh(2,itbase+1) = -3
            if (j.ne.ngridy) then
               vertex(1,itbase+2) = ivbase+ngridy+2
               vertex(2,itbase+2) = ivbase+2
               vertex(3,itbase+2) = ivbase+1
               if (i.eq.1) neigh(1,itbase+2) = -1
               vertex(1,itbase+3) = ivbase+ngridy+2
               vertex(2,itbase+3) = ivbase+2
               vertex(3,itbase+3) = ivbase+ngridy+3
               if (j.eq.ngridy-1) neigh(1,itbase+3) = -4
               if (i.eq.ngridx) neigh(2,itbase+3) = -3
            endif
            if (i.ne.ngridx) then
               vertex(1,itbase+2*ngridy) = ivbase+ngridy+2
               vertex(2,itbase+2*ngridy) = ivbase+2*ngridy+2
               vertex(3,itbase+2*ngridy) = ivbase+ngridy+1
               if (j.eq.1) neigh(1,itbase+2*ngridy) = -2
               vertex(1,itbase+2*ngridy+1) = ivbase+ngridy+2
               vertex(2,itbase+2*ngridy+1) = ivbase+2*ngridy+2
               vertex(3,itbase+2*ngridy+1) = ivbase+2*ngridy+3
               if (i.eq.ngridx-1) neigh(1,itbase+2*ngridy+1) = -3
               if (j.eq.ngridy) neigh(2,itbase+2*ngridy+1) = -4
               if (j.ne.ngridy) then
                  vertex(1,itbase+2*ngridy+2) = ivbase+ngridy+2
                  vertex(2,itbase+2*ngridy+2) = ivbase+2*ngridy+4
                  vertex(3,itbase+2*ngridy+2) = ivbase+ngridy+3
                  if (j.eq.ngridy-1) neigh(1,itbase+2*ngridy+2) = -4
                  vertex(1,itbase+2*ngridy+3) = ivbase+ngridy+2
                  vertex(2,itbase+2*ngridy+3) = ivbase+2*ngridy+4
                  vertex(3,itbase+2*ngridy+3) = ivbase+2*ngridy+3
                  if (i.eq.ngridx-1) neigh(1,itbase+2*ngridy+3) = -3
               endif
            endif
   20    continue
   21 continue
c
      if (outlev.ge.3) then
         write(ioutpt,102) ngridx,ngridy,ntri,nvert
  102    format(' grid lines   ',i3,' X',i3/
     1          ' triangles    ',i5/
     2          ' vertices     ',i5/
     3          ' initial triangulation complete')
      endif
      return
      end