|
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** |