C******************************************************************* C** C** v e m g e n 2 d : C** C** generation of a triangulation of 2-dimensional domain C** with well-known boundaries. The program can easily changed for C** other problems. The mesh data are read by vemu02. C** C** by L. Grosz Karlsruhe, Oct. 1994 C** C******************************************************************* C** C** B3 C** (0,1) (1,1) C** *------* C** | | C** B4 | | B2 C** | | C** *------* C** xi2 ^ (0,0) (1,0) C** | B1 C** --xi1 C** C** C** This FORTRAN program generates an order 2 subdivision of a C** two dimensional domain into triangle elements. The domain C** has a parametrical representation whose domain is the unit C** cube [0,1]^2 in the variables (xi1,xi2). There are four C** boundary portions numbered from 1 to 4. Boundary 1 is C** represented by xi2=0, boundary 2 is represented by xi1=1, C** boundary 3 is represented by xi2=1 and boundary 4 is C** represented by xi1=0. The parametrical representation of C** the domain is computed from the parametrical representation C** F1 and F4 of boundary 1 and 4 by a linear interpolation. C** The code can be adapted to other domains by entering suitable C** statements for F1 and F4 into the node generation. C** C** The boundaries are subdivided into line elements of order 2 C** and additionally all boundary nodes are specified as nodes C** with Dirichlet conditions. Since the Dirichlet conditions C** has the higher priority, the line elements will be useless C** in the vecfem application. But if you want to activate C** special boundary elements, you can remove the C** corresponding boundary nodes from the generation of the C** 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,IND DOUBLE PRECISION X1,X2,X3,XI1,XI2,XI,F11,F12,F31,F32 IND=1 C** C**----------------------------------------------------------------- C** PRINT*,'Enter number of elements in x1- and x2- direction:' READ(5,*) ELEM1,ELEM2 PRINT*,'Name of the mesh file :' READ(5,'(80A)') MESH C** C**----------------------------------------------------------------- C** C** open output file : C** 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 of the C** parametrical representation of the actual domain, where Ni C** is the number points in XIi direction (i=1,2): C** N1=2*ELEM1+1 N2=2*ELEM2+1 WRITE(99,*) N1*N2 C** DO 10 Z2=1,N2 DO 10 Z1=1,N1 XI1=DBLE(Z1-1)/DBLE(N1-1) XI2=DBLE(Z2-1)/DBLE(N2-1) C** XI=XI1 C** C**----------------------------------------------------------------- C** C** here you can change the shape of the domain : C** C** B3 C** F3(0) F3(1) if you set : C** *-------* C** | | F11=XI F12=0 C** | | F31=XI F32=1 C** | | C** *-------* you will get the unit cube [0,1]^2. C** F1(0) F1(1) C** B1 C** F11=XI F12=0 F31=XI F32=1 C** C**----------------------------------------------------------------- C** X1=(1.-XI2)*F11+F31*XI2 X2=(1.-XI2)*F12+F32*XI2 X3=0 WRITE(99,*) Z1+N1*(Z2-1),X1,X2,X3 10 CONTINUE C** PRINT*,'written nodes : ',N1*N2 C** C**----------------------------------------------------------------- C** C**** the generation of the elements : C** ------------------------------- C** C** The domain is covered by triangle elements of order 2 C** and consequently the boundaries are described by line elments C** of order 2. The following picture illustrates the C** construction of the triangle elements with lower node S C** and their boundary elements which are only generated C** if they are a subset of the boundary of the domain: C** C** S+2*N1 S+2*N1+1 S+2*N1+2 C** 3-----2-----1 C** 1 3-----6-----1 2 C** | | \ | | C** | | \ | | C** S+N1 3 6 5 4 3 S+N1+2 C** | | \ | | C** | | \ | | C** 2 1-----4-----2 1 C** 1-----3-----2 C** S S+1 S+2 C** C** There are two different element types: C** WRITE(99,*) 2 C** C** These are the triangle elements: C** C** ELMID is the element id number, which allows to recognize C** an element after the distribution to the processors. C** WRITE(99,*) 2*ELEM1*ELEM2,2,3,6 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,S+1,S+N1+1,S+N1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2 WRITE(99,*) ELMID,IND, & S+2*N1+2,S+2,S+2*N1,S+N1+2,S+N1+1,S+2*N1+1 20 CONTINUE C** PRINT*,'written triangle elements : ',2*ELEM1*ELEM2 C** C**-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- C** C** These are the line elements: C** C** The orientation of the line elements is counterclockwise. C** WRITE(99,*) 2*(ELEM1+ELEM2),1,2,3 C** DO 31 Z1=1,ELEM1 C** C**** elements on boundary 1: C** S=2*(Z1-1)+1 ELMID=Z1+2*ELEM1*ELEM2 WRITE(99,*) ELMID,IND,S,S+2,S+1 C** C**** elements on boundary 3: C** S=2*(Z1-1)+2*(ELEM2-1)*N1+2*N1+1 ELMID=Z1+ELEM1+2*ELEM1*ELEM2 WRITE(99,*) ELMID,IND,S+2,S,S+1 C** 31 CONTINUE DO 32 Z2=1,ELEM2 C** C**** elements on boundary 2: C** S=2*(ELEM1-1)+2*(Z2-1)*N1+3 ELMID=Z2+ELEM2+2*ELEM1+2*ELEM1*ELEM2 WRITE(99,*) ELMID,IND,S,S+2*N1,S+N1 C** C**** elements on boundary 4: C** S=2*(Z2-1)*N1+1 ELMID=Z2+2*ELEM1+2*ELEM1*ELEM2 WRITE(99,*) ELMID,IND,S+2*N1,S,S+N1 32 CONTINUE C** PRINT*,'written line elements : ',2*(ELEM1+ELEM2) C** C**----------------------------------------------------------------- C** C** generation of the nodes with Dirichlet conditions : C** ------------------------------------------------- C** C** The real parameter is used to notify the boundary id number C** where the node with Dirichlet condition belongs to. Here C** all nodes on the boundary gets a Dirichlet condition for C** all components, so that the boundary elements are useless. C** It is easy to change the program, so that special boundary C** portions get no Dirichlet conditions for a component. Then C** boundary condition of the Neuman type can be considered. C** WRITE(99,*) NK C** DO 40 D=1,NK C** WRITE(99,*) 2*(N1+N2-2) C** DO 41 Z1=1,N1-1 WRITE(99,*) Z1,DBLE(1) WRITE(99,*) N1*N2+1-Z1,DBLE(3) 41 CONTINUE DO 42 Z2=1,N2-1 WRITE(99,*) Z2*N1,DBLE(2) WRITE(99,*) (N2-Z2)*N1+1,DBLE(4) 42 CONTINUE C** PRINT*,2*(N1+N2-2),' Dirichlet conditions for component ',d 40 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C**