subroutine generate(pri,isize,ierr)
*
*********************************************************
*     mesh generation				
*********************************************************
*
      implicit double precision (a-h,o-z)
      dimension pri(*)
      common /files/ infile,iofile,iscree
*
      write(6,*) 
      write(6,*) ' Start mesh generation:'
*
*********************************************************
*     generate triangular mesh						
*********************************************************
*
      call gener1(infile,iscree,pri,isize,ierr)
      if (ierr.ne.0) go to 9000
*
 9000 continue
      write(6,*) ' End mesh generation'
      return
*
      end


      subroutine gener1(infile,iw,pri,isize,ierr)
*
*********************************************************
*     pave the domain with triangles
*********************************************************
*
      implicit double precision (a-h,o-z)
      dimension pri(*), rquad(3,10000)
      common /misc/ maxpo,numpo,numse,numint,
     & maxsb,numbo,nsab,numms,numcseg,nexmax
      common /point/ lacoep,laiseg,laibnd,lalbnd,
     & larpbg,laics,lacseg
*
*********************************************************
*     maxb - max number of boundaries				
*     maxp - max number of points 					
*********************************************************
*
      call gener11(iw,3,numse,pri(laiseg),
     i 4,numbo,pri(laibnd),pri(lalbnd),
     o maxb,maxp,ierr)
      if (ierr.ne.0) go to 255
*
*********************************************************
*     create work array that can store (u,v) for all points on all
*     boundaries on a surface. create work array that can store
*     pointers for each boundary to the aforementioned array and pointers
*     to segment list for each boundary on a surface.
*     create work array that can store point numbers for all points on
*     all boundaries on a surface and segment numbers for all segments
*     on all boundaries on a surface.
*********************************************************
*
      call dimpri(pri,maxb,1,4,lalpbn,1)
      call dimpri(pri,maxb,1,4,lalpbs,1)
      call dimpri(pri,maxp,1,4,laidpt,1)
      call dimpri(pri,maxp,1,4,laidsg,1)
*
*********************************************************
*     set up globally defined arrays		 
*********************************************************
*
      mheap=isize/100
      mface=isize/40
      mele=isize/100
      mnode=isize/100
*
      call dimpri(pri,2,mnode,8,lacoeg,1)
      call dimpri(pri,7,mnode,4,lalqua,1)
      call dimpri(pri,3,mele,4,laitno,1)
      call dimpri(pri,3,mele,4,laited,1)
      call dimpri(pri,3,mface,4,lalfac,1)
      call dimpri(pri,4,mface,8,lavfac,1)
      call dimpri(pri,1,mface,4,lalpoi,1)
      call dimpri(pri,3,mface,4,lalfap,1)
      call dimpri(pri,1,mface,8,larfac,1)
      call dimpri(pri,1,mface,4,lalhea,1)
      call dimpri(pri,4,mnode,4,lanode,1)
      call dimpri(pri,mface,1,4,lagseg,1)
      call dimpri(pri,mface,1,4,laiare,1)
      call dimpri(pri,1,mnode,4,laibou,1)
*
      call backgr(numpo,pri(lacoep),pri(lalqua),rquad)
      hmax=rquad(3,1)/25.d0
*
      npoin=numpo
*
      nquad=0
      nele=0
      nface=0
*
*********************************************************
*     fill work array with point numbers and pointers 
*********************************************************
*
      call gener12(iw,3,numse,pri(laiseg),pri(lacoep),
     & 4,numbo,pri(laibnd),pri(lalbnd),
     & pri(lalpbn),pri(lalpbs),pri(laidsg),pri(laidpt),
     & numseg,numpo,ierr)
*
      if (ierr.ne.0) go to 255
*
      call phcoeg(pri(lacoeg),pri(lacoep),numpo,
     & pri(larpbg),numms,hmax,hmin,xmax,xmin,ymax,ymin)
*
      dl=dmax1(xmax-xmin,ymax-ymin)
      hloc=dl/10.d0
*
      if (ierr.ne.0) go to 255
*
*
*********************************************************
*     triangulate the boundary			
*********************************************************
*
*********************************************************
*     file to store the boundary normal vectors
*********************************************************
*
      ns=17
      open(ns,file='normal.scratch',access='sequential',
     & form='unformatted',status='unknown')
*
      call boundary(
     i maxpo,pri(laiseg),numbo,
     i pri(lalpbn),pri(laidpt),pri(lalpbs),pri(laidsg),
     i hmax,hmin,numpo,numseg,numms,pri(larpbg),
     i pri(laics),pri(lacseg),numcseg,
     m nface,pri(lalfac),pri(lavfac),pri(lacoeg),
     m pri(lagseg),pri(laibou),pri(laiare),pri(lalqua),
     m nquad,rquad,npoin,maxseg)
*
      call lfapc(pri(lalfac),nface,npoin,
     x pri(lalpoi),pri(lalfap),nfapo)
*
      nquad=nquad
      npbou=npoin
      nfbou=nface
      nfacv=nface
*
      call arheap(pri(laiare),pri(lalfac),pri(lacoeg),
     x pri(lalhea),pri(larfac),nheap)
*
*********************************************************
*     triangulate the interior			
*********************************************************
*
*      mele=300
      call gridc(nheap,nface,pri(lalfac),pri(lacoeg),pri(lalhea),
     x pri(larfac),pri(lalpoi),pri(lalfap),nfapo,pri(laitno),
     x pri(lalqua),npoin,rquad,nquad,nele,mele,ierr,
     x pri(laibou),hmax,hmin,iums,numms,pri(larpbg))
*
      if(ierr.ne.0) goto 255
*
      call phsmo(pri(lalpoi),pri(lalfap),
     x     pri(lalfac),pri(lacoeg),pri(laibou),npoin)
*
*********************************************************
*     write to file				
*********************************************************
*
  255 continue
      if(ierr.eq.12) ierr=0
*
      call wrtr(pri,ierr,pri(lacoeg),nele,npoin,
     $     pri(laitno),pri(laited),pri(lavfac),nfacv,
     $     pri(lagseg),pri(lalfac),maxseg,infile)
*
      return
      end
*
*
*
      integer function inpi(in)
       inpi=in
      return
      end
*
*
*
      subroutine gener11(
     i iw,nrowse,numse,nseg,
     i nrowbn,numbo,nbnd,lbnd,
     o maxb,maxp,ierr)
*
*********************************************************
*     get number of boundaries 						
*     get max number of points on all boundaries	 
*********************************************************
*
      implicit double precision (a-h,o-z)
      dimension nseg(nrowse,*), nbnd(nrowbn,*), lbnd(*)
*
      maxp = 0
      maxb = numbo
*
*********************************************************
*     count number of points on all segments 		
*     on all boundaries			
*********************************************************
*
      mp = 0
      do 40 ib = 1, numbo
*
      mseg = nbnd(2,ib)
      lp3  = nbnd(3,ib) - 1
*
*********************************************************
*     add (number of points)-1 for each segment 	
*     on boundary ib	
*********************************************************
*
      do 20 i = 1, mseg
      name = lbnd(lp3+i)
      do 15 j = 1, numse
      lp4 = j
      if (nseg(1,j).eq.name) go to 16
   15 continue
      go to 8020
   16 continue
      mp = mp + 2
   20 continue
*
   40 continue
      if (mp.gt.maxp) maxp = mp
*
      go to 9000
*
 8020 continue
      ierr = 112
      write(6,*) ' Error: a segment has not been defined, no.',name
c
 9000 continue
      return
      end
*
*
*
      subroutine gener12(
     i     iw,nrowse,numse,nseg,coep,
     i     nrowbn,numbo,nbnd,lbnd,
     o     lpbnd,lpbnds,idseg,idpnt,numseg,numpo,ierr)
*     
*********************************************************
*     load boundary information 	
*********************************************************
*     
      implicit double precision (a-h,o-z)
      logical icha,iclose
      dimension nseg(nrowse,*), coep(2,*)
      dimension nbnd(nrowbn,*), lbnd(*)
      dimension lpbnd(*), lpbnds(*), idseg(*), idpnt(*)
*     
      lps = 1
      lpp = 1
*     
*********************************************************
*     for each boundary		
*********************************************************
*     
      do 400 ib = 1, numbo
         lpps = lpp
         lpss = lps
         iib = ib
         lpbnds(ib) = lps
         lpbnd(ib)  = lpp
         mseg = nbnd(2,ib)
         lp3  = nbnd(3,ib) - 1
         ityp  = nbnd(4,ib)
*     
*********************************************************
*     for each segment on this boundary				
*********************************************************
*     
         do 300 i = 1, mseg
            iclose=.true.
            name = lbnd(lp3+i)
            do j = 1, numse
               lp4 = j
               if (nseg(1,j).eq.name) go to 16
            enddo
            go to 8020
 16         continue
*     
*     store segment number 
*     
            idseg(lps) = lp4
            lps = lps + 1
*     
*********************************************************
*     the points on segment	
*********************************************************
*     
            npiol1=npi1
            npiol2=npi2
            npi1 = nseg(2,lp4)
            npi2 = nseg(3,lp4)
*     
*********************************************************
*     find first and last point on the following segment*
*********************************************************
*     
            k = i + 1
            if (k.gt.mseg) k = 1
            name = lbnd(lp3+k)
            do 25 j = 1, numse
               lp5 = j
               if (nseg(1,j).eq.name) go to 26
 25         continue
            go to 8020
 26         continue
*     
            npk1 = nseg(2,lp5)
            npk2 = nseg(3,lp5)
*     
*********************************************************
*     determine start point npis = npi1 or npi2 	
*     (whichever is not equal to npk1 or npk2)		
*********************************************************
*     
            nbeg = 2
            if (npi1.ne.npk1.and.npi1.ne.npk2) go to 50
            nbeg = 3
 50         continue
*     
*********************************************************
*     store the first point of each segment			
*********************************************************
*     
            idpnt(lpp) = nseg(nbeg,lp4)
            itest=0
            if(mseg.le.2.and.i.eq.mseg) itest=1
            if ((npi1.ne.npk1.and.npi1.ne.npk2).and.
     &           (npi2.ne.npk1.and.npi2.ne.npk2)) itest=1
*
            if(itest.eq.1) then
               npp=npiol1
               if(npiol1.ne.npi1.and.npiol1.ne.npi2) npp=npiol2
               idpnt(lpp)=npp
               lpp=lpp+1
               idpnt(lpp)=-npi1
               if(npi1.eq.npp) idpnt(lpp)=-npi2
               iclose=.false.
            endif
            lpp = lpp + 1
*     
 300     continue
*     
*********************************************************
*     check boundary type and direction				
*********************************************************
*     
         npoin = lpp - lpps
         call checkdir(
     i        npoin,coep,idpnt(lpps),
     o        area)
*     
         icha = .false.
         if(area.lt.0.d0.and.ityp.eq.1)  icha = .true.
         if(area.gt.0.d0.and.ityp.eq.2)  icha = .true.
         if(ityp.eq.4)  icha = .true.
*         call dwmatg(idpnt(lpps),1,lpp-1,0,'i    ','idpnt before')
*         call dwmatg(idseg(lpss),1,lps-1,0,'i    ','idseg before')
         if(icha) then
*     
*********************************************************
*     change direction of boundary					
*********************************************************
*     
            nend=lpp-1
            nstart=lpps+1
            if(iclose) then
               nstart=lpps+1
            else
               nstart=lpps
            endif
            nmid=(nend+nstart)/2
            i1=0
            do 305 i=nstart,nmid
               i1=i1+1
               ik=idpnt(i)
               idpnt(i)=idpnt(nend-i1+1)
               idpnt(nend-i1+1)=ik
 305        continue
            if(idpnt(nstart).lt.0) then
               idpnt(nstart)=-idpnt(nstart)
               idpnt(nend)=-idpnt(nend)
            endif
            nend=lps-1
            nstart=lpss
            nmid=(nend+nstart)/2
            i1=0
            do 306 i=nstart,nmid
               i1=i1+1
               ik=idseg(i)
               idseg(i)=idseg(nend-i1+1)
               idseg(nend-i1+1)=ik
 306        continue
         endif
*     
 400  continue
      go to 9000
*     
 8020 continue
      ierr = 112
      write(6,*) ' Error: a segment has not been defined, no.',name
      go to 9000
 8030 continue
      ierr = 113
      write(6,*) 
     $     ' Error: segments do not form ',
     $     'an unbroken chain, boundary no.',ib

*     
 9000 continue
      numseg= lps - 1
      numpo = lpp - 1
*      call dwmatg(idpnt,1,numpo,0,'i    ','idpnt')
*      call dwmatg(idseg,1,numseg,0,'i    ','idseg')
*      stop
      return
      end


      subroutine wrtr(pri,ierr,coeg,nele,nno,itri,itedge,
     $     vface,nface,lseg,lface,maxs,infile)
*     
*********************************************************
*     write mesh on file "solu.scratch"					
*********************************************************
*     
      implicit double precision (a-h,o-z)
*     
      dimension pri(*),coeg(2,*),itri(3,*),lface(3,*),lseg(2,*)
      dimension itedge(3,*),vface(4,*),iw(10)
      real x(2)
*     
*
*********************************************************
*     create edge connections   
*********************************************************
*
      call dimpri(pri,4,nno*4,4,laa,1)
      call nodlis(nno,nele,pri(laa),itri,nlist)
      call edgcr(nele,pri(laa),itri,itedge)
      call edgfac(itri,itedge,lface,nele,nface)
      call dimpri(pri,4,nno*4,4,laa,-1)
*
      open(unit=61,file='solu.scratch',status='unknown')
      open(unit=62,file='solu2.scratch',status='unknown')
      write(61,*) nele,nno
      write(61,*) ((itri(i,j),i=1,3),j=1,nele)
      write(61,*) ((coeg(i,j),i=1,2),j=1,nno)
      close(unit=61)
*******************************************
      write(62,*) nele,nno,nface
      write(62,*) ((itedge(i,j),i=1,3),j=1,nele)
      write(62,*) ((vface(i,j),i=1,4),j=1,nface)
      close(unit=62)
*******************************************
      backspace(infile) 
      write(infile,90) 
      write(infile,*) ierr
      if(ierr.ne.0) return
      write(infile,100) 
      write(infile,*) nno,nele
      do 10 i=1,nno
         x(1)=coeg(1,i)
         x(2)=coeg(2,i)
         write(infile,*) (x(j),j=1,2)
 10   continue
      do 20 i=1,nele
         write(infile,*) (itri(j,i),j=1,3)
 20   continue
*     
*     write points on segment
*     
      ns=78
      open(ns,file='temp.scratch',access='sequential',
     &     form='unformatted',status='unknown')
*     
      i1=0
      do 30 iseg=1,maxs
         if(lseg(1,iseg).gt.0) then
            i1=i1+1
            npoin=lseg(1,iseg)+1
            write(ns) iseg,npoin
            istart=lseg(2,iseg)
*     
            ich=1
            nod1=lface(1,istart)
            n1=lface(1,istart+1)
            n2=lface(2,istart+1)
            if(nod1.eq.n1.or.nod1.eq.n2) ich=2
            do 25 i=1,npoin-1
               nod1=lface(ich,istart+i-1)
               write(ns) nod1
 25         continue
*     
            ich=ich-1
            if(ich.eq.0) ich=2
            nod1=lface(ich,istart+npoin-2)
            write(ns) nod1
         endif
 30   continue
      rewind(ns)
*     
      write(infile,150) 
      write(infile,*) i1
      do 40 i=1,i1
         read(ns) iseg,npoin
         write(infile,*) iseg,npoin
         nlin = npoin/10
         nextra = npoin - nlin*10
         if(nlin.eq.0) goto 33
         do 31 j=1,nlin
            do 32 k=1,10
               read(ns) ipoin
               iw(k)=ipoin
 32         continue
            write(infile,*) (iw(k),k=1,10)
 31      continue
 33      continue
         if(nextra.eq.0) goto 35
         do 34 j=1,nextra
            read(ns) ipoin
            iw(j)=ipoin
 34      continue
         write(infile,*) (iw(k),k=1,nextra)
 35      continue
 40   continue
      close(unit=ns,status='delete')
*     
 160  continue
      write(infile,200) 
      close (infile,status='keep',err=8020)
      goto 9000
 8020 continue
      ierr = 3
      write(6,*) ' Error: failed to close input file'
      go to 9000
 9000 continue
      return
 90   format('ERROR')
 100  format('MESH DATA')
 150  format('NODES ON SEGMENTS')
 200  format('END OF DATA')
      end


      subroutine phcoeg(coeg,coor,nno,rpbg,numms,hmax,hmin,
     & xmax,xmin,ymax,ymin)
      implicit double precision (a-h,o-z)
*
      dimension coeg(2,*),coor(2,*),rpbg(numms,*)
*
      xmax=-1.d9
      ymax=-1.d9
      xmin=1.d9
      ymin=1.d9
*
      do 10 i=1,nno
         u = coor(1,i)
         v = coor(2,i)
         coeg(1,i) = u
         coeg(2,i) = v
         xmax=dmax1(xmax,u)
         xmin=dmin1(xmin,u)
         ymax=dmax1(ymax,v)
         ymin=dmin1(ymin,v)
 10   continue
*
      hmax=-9.d10
      hmin=9.d10
*
      do 20 i = 1, numms
         hmax=dmax1(hmax,rpbg(i,3))
         hmin=dmin1(hmin,rpbg(i,3))
   20 continue
      return
      end
      subroutine nodlis(nno,nele,nodel,itri,nlist)
*
      implicit double precision (a-h,o-z)
      dimension nodel(4,*),itri(3,*)
*
      imax=nno+1
      nlist=nno
      do 50 iel=1,nele
      do 40 i=1,3
      inod=itri(i,iel)
10    continue
      ik=nodel(4,inod)
      if(ik.lt.0) inod=-nodel(4,inod)
      if(ik.lt.0) goto 10
      if(ik.eq.3) goto 20
      nodel(ik+1,inod)=iel
      nodel(4,inod)=ik+1
      goto 40
20    continue
      nodel(4,inod)=-imax
      nodel(1,imax)=iel
      nodel(4,imax)=1
      nlist=imax
      imax=imax+1
40    continue
50    continue
      return
      end
      subroutine edgcr(nele,nodel,nodes,nelcon)
*     
      implicit double precision (a-h,o-z)
      dimension nodel(4,*),nodes(3,*),nelcon(3,*)
*     
*********************************************************
*     create edge connections   
*********************************************************
*     
      do 50 iel=1,nele
         iv1=nodes(1,iel)
         iv2=nodes(2,iel)
         iv3=nodes(3,iel)
         do 40 i=1,3
*     
*********************************************************
*     first loop for iv2   
*********************************************************
*     
            inod1=iv2
 10         continue
            ik1=nodel(4,inod1)
            inr1=ik1
            if(ik1.lt.0) inr1=3
            do 30 ites1=1,inr1
               iel1=nodel(ites1,inod1)
               if(iel1.eq.iel) goto 30
*     
*********************************************************
*     second loop for iv3  
*********************************************************
*     
               inod2=iv3
 20            continue
               ik2=nodel(4,inod2)
               inr2=ik2
               if(ik2.lt.0) inr2=3
               do ites2=1,inr2
                  iel2=nodel(ites2,inod2)
                  if(iel1.eq.iel2) goto 35
               enddo
               if(ik2.lt.0) inod2=-ik2
               if(ik2.lt.0) goto 20
 30         continue
            if(ik1.lt.0) inod1=-ik1
            if(ik1.lt.0) goto 10
            nelcon(i,iel)=-1
            goto 36
 35         continue
            nelcon(i,iel)=iel1
 36         ii=iv1
            iv1=iv2
            iv2=iv3
            iv3=ii
 40      continue
 50   continue
      return
      end
      subroutine edgfac(itri,itedge,lface,nele,nface)
      dimension itri(3,*),itedge(3,*),lface(3,*)
*
*     Give the right boundary segment number to itedge(3,*)
*
      do i=1,nele
         do j=1,3
            if(itedge(j,i).eq.-1) then
*
               do iface=1,nface
                  nodco=0
                  do  k=1,3
                     inod=itri(k,i)
                     do  l=1,2
                        inod2=lface(l,iface)
                        if(inod.eq.inod2) nodco=nodco+1
                     enddo
                  enddo
                  if (nodco.ne.2) goto 9000
                  do 10 k=1,3
                     nodnr=itri(k,i)
                     if(nodnr.eq.lface(1,iface).or.
     $                    nodnr.eq.lface(2,iface)) goto 10
                     goto 9000
 10               continue
 9000             continue
                  if(itri(j,i).eq.nodnr.and.nodco.eq.2) goto 20
               enddo
 20            continue
               itedge(j,i)=-iface
            endif
*
         enddo
      enddo
      return
      end