C*******************************************************************
C**
C**    v e m g e n 2 d q:
C**
C**  subdivision of the domain [0,1]^2 into quadrilateral elements.
C**  The generated mesh is used by example vembldexm10. The mesh
C**  data are read by vemu02.
C**
C**   by L. Grosz                           Karlsruhe, June 1995
C**
C*******************************************************************
C**
C**              (0,1)    (1,1)
C**                 *------*<--boundary X2=1
C**                 |      |
C**                 |      |
C**                 |      |
C**                 *------*<--boundary X2=0
C**  x2  ^       (0,0)    (1,0)
C**      |
C**      -->x1
C**
C**
C**  This FORTRAN program generates an order 2 subdivision of the
C**  two dimensional unit cube into quadrilateral elements.
C**  The boundary portion X2=1 is subdivided into line elements of
C**  order 2 and all nodes on the boundary X2=0 are specified as
C**  nodes with Dirichlet conditions.
C**
      PROGRAM VEMEXM
C**
C**-----------------------------------------------------------------
C**
      IMPLICIT NONE
C**
C**-----------------------------------------------------------------
C**
C**    some parameters which may be chanced:
C**
C**    ELEM1 = number of elements in x1 direction
C**    ELEM2 = number of elements in x2 direction
C**    NK    = number of solution components
C**    MESH  = name of the mesh file
C**
      INTEGER        ELEM1,ELEM2,NK
      CHARACTER*80   MESH
      PARAMETER (NK=1)
C***
      INTEGER          Z1,Z2,N1,N2,S,ELMID,D,NDEG,NODNUM,NDC,IND
      DOUBLE PRECISION NOD1,NOD2,NOD3
      IND=1
C**
C**-----------------------------------------------------------------
C**
      PRINT*,'Enter number of elements in x1- and x2- direction:'
      READ(5,*) ELEM1,ELEM2
C**
C**-----------------------------------------------------------------
C**
C**   open output file :
C**
      MESH="mesh.vem"
      OPEN (99,FILE=MESH,STATUS= 'UNKNOWN',FORM='FORMATTED')
      PRINT*,'opened file : ',MESH
C**
C**-----------------------------------------------------------------
C**
C**** the generation of the geometrical nodes :
C**   ---------------------------------------
C**
C**  The grid is rectangular in the domain [0,1]^2, where Ni
C**  is the number points in Xi direction (i=1,2 ):
C**
      N1=2*ELEM1+1
      N2=2*ELEM2+1
      NDEG=N1*N2
      WRITE(99,*) NDEG
C**
      DO 10 Z2=1,N2
        DO 10 Z1=1,N1
          NODNUM=Z1+N1*(Z2-1)
          NOD1=DBLE(Z1-1)/DBLE(N1-1)
          NOD2=DBLE(Z2-1)/DBLE(N2-1)
	  NOD3=0
          WRITE(99,*) NODNUM,NOD1,NOD2,NOD3
 10   CONTINUE
C**
      PRINT*,'written nodes : ',NDEG
C**
C**-----------------------------------------------------------------
C**
C**** the generation of the elements :
C**   -------------------------------
C**
C**
C**   The domain is covered by quadrilateral elements of order 2
C**   and consequently the upper boundary of the domain is
C**   described by line elments of order 2.
C**
C**   There are two different element types:
C**
      WRITE(99,*) 2
C**
C**   These are the quadrilateral elements:
C**   The following picture illustrates the construction of the
C**   quadrilateral element with lower node S :
C**
C**       S+2*N1   S+2*N1+1    S+2*N1+2
C**             4-----7-----3
C**             |           |
C**             |           |
C**      S+N1   8           6   S+N1+2
C**             |           |
C**             |           |
C**             1-----5-----2
C**            S     S+1        S+2
C**
C**   ELMID is the element id number, which allows to recognize
C**   an element after the distribution to the processors.
C**
      WRITE(99,*) ELEM1*ELEM2,2,4,8
      DO 20 Z2=1,ELEM2
        DO 20 Z1=1,ELEM1
          S=2*(Z1-1)+2*(Z2-1)*N1+1
	  ELMID=Z1+ELEM1*(Z2-1)
	  WRITE(99,*) ELMID,IND,S,S+2,S+2*N1+2,S+2*N1,
     &                      S+1,S+N1+2,S+2*N1+1,S+N1
 20   CONTINUE
C**
      PRINT*,'written quadrilateral elements : ',ELEM1*ELEM2
C**
C**-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
C**
C**   These are the line elements:
C**   The following picture illustrates the construction of the
C**   line element with right node S :
C**
C**             2-----3-----1
C**            S     S+1     S+2
C**
C**   The orientation of the line elements is counterclockwise.
C**
      WRITE(99,*) ELEM1,1,2,3
      DO 31 Z1=1,ELEM1
        S=2*(Z1-1)+(N2-1)*N1+1
	ELMID=Z1+ELEM1*ELEM2
	WRITE(99,*) ELMID,IND,S+2,S,S+1
 31   CONTINUE
C**
      PRINT*,'written line elements : ',ELEM1
C**
C**-----------------------------------------------------------------
C**
C**   generation of the nodes with Dirichlet conditions :
C**   -------------------------------------------------
C**
C**   All nodes in the lower boundary X2=0 get a Dirichlet
C**   condition. The prescribed value is 10. for all nodes.
C**
      WRITE(99,*) NK
C**
      DO 40 D=1,NK
	NDC=N1
        WRITE(99,*) NDC
        DO 41 Z1=1,N1
          WRITE(99,*) Z1,10.
   41   CONTINUE
C**
        PRINT*,N1,' Dirichlet conditions for component ',d
40    CONTINUE
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C**   ------------------
C**