C******************************************************************* C** C** v e m g e n 3 d : C** C** generation of a subdivision of a 3-dimensional domain with C** well-known boundaries into hexahedron elements. The program C** can easily changed for other problems. The mesh data is read C** by vemu02. C** C** by L. Grosz Karlsruhe, Oct. 1994 C** C** C******************************************************************* C** B4 C** (0,1,1) / (1,1,1) C** *----------* C** / : B6 / | C** / : / | C** B5 *----------* B3 | C** \ | *----|-----* (1,1,0) C** | / | / C** | / B2 | / C** *----------* C** xi3 ^ xi2 \ (1,0,0) C** | / B1 C** --xi1 C** C** This FORTRAN program generates an order 2 subdivison of a C** three dimensional domain into hexahedron elements. The domain C** has a parametrical representation whose domain is the unit C** cube [0,1]^3 in the variables (xi1,xi2,xi3). There are six C** boundary portions numbered from 1 to 6. Boundary 1 is C** represented by xi3=0, boundary 2 is represented by xi2=0, C** boundary 3 is represented by xi1=1, boundary 4 is represented C** by xi2=1, boundary 5 is represented by xi1=0 and boundary 6 C** is reprsented by xi3=1. The parametrical representation C** of the domain is computed from the parametrical representions C** F1 and F6 of boundary 1 and 6 by a linear interpolation. C** The code can be adapted to other domains by entering suitable C** statements for F1 and F6 into the node generation. C** C** The boundaries are subdivied into quadrilateral elements of C** order 2 and additionally all boundary nodes are specified as C** nodes with Dirichlet conditions. Since the Dirichlet C** conditions has the higher priority, the line elements will be C** useless 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** ELEM3 = number of elements in x3 direction C** NK = number of solution components C** MESH = name of the mesh file C** INTEGER ELEM1,ELEM2,ELEM3,NK CHARACTER*80 MESH PARAMETER (NK=4) C*** INTEGER Z1,Z2,Z3,N1,N2,N3,S,ELMID,D,IND DOUBLE PRECISION X1,X2,X3,XI1,XI2,XI3,F11,F12,F13,F61,F62,F63,PI, & ZERO,VALUE ZERO=0 IND=1 C** C**----------------------------------------------------------------- C** PRINT*,'Enter number of elements in x1- , x2- and x3- direction:' READ(5,*) ELEM1,ELEM2,ELEM3 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]^3 of the C** parametrical representation of the actual domain, where Ni C** is the number points in XIi direction (i=1,2,3): C** N1=2*ELEM1+1 N2=2*ELEM2+1 N3=2*ELEM3+1 WRITE(99,*) N1*N2*N3 C** DO 10 Z3=1,N3 DO 10 Z2=1,N2 DO 10 Z1=1,N1 XI1=DBLE(Z1-1)/DBLE(N1-1) XI2=DBLE(Z2-1)/DBLE(N2-1) XI3=DBLE(Z3-1)/DBLE(N3-1) C** C**----------------------------------------------------------------- C** C** here you can change the shape of the domain : C** C** C** F6(0,1) *----------* F6(1,1) C** / : B6 / | C** / : / | if you set : C** F6(0,0) *----------* F6(1,0) C** | *----|-----* F1(1,1) F11=XI1 F12=XI2 F13=0 C** | / | / F61=XI1 F62=XI2 F63=1 C** | / | / C** F1(0,0) *----------* F1(1,0) you will get the unit C** / cube [0,1]^3. C** B1 C** C** Here we generate a curved channel with the unit square C** as sectional area. There are straight inlet and outlet C** of length 1. The inner radius of the curved section is C** 1 and the outer radius is 2. C** C** ------, C** / , C** 3 -|----- , , C** | ' , C** 2 -|----- . . , , C** . , | C** 1 -- | | | C** / | |/ C** 0 -| | ------ C** 0 1 2 3 C** PI=4.*ATAN(1.D0) IF (XI1.LE.1.D0/3) THEN F11=XI1*3 F12=XI2 F13=2. F61=XI1*3 F62=XI2 F63=3. ELSEIF (XI1.LE.2.D0/3) THEN F11=1.D0+SIN((XI1*3-1)*PI/2) F12=XI2 F13=1.D0+COS((XI1*3-1)*PI/2) F61=1.+2.D0*SIN((XI1*3-1)*PI/2) F62=XI2 F63=1.+2.D0*COS((XI1*3-1)*PI/2) ELSE F11=2. F12=XI2 F13=3*(1.-XI1) F61=3. F62=XI2 F63=3*(1.-XI1) ENDIF C** C**----------------------------------------------------------------- C** X1=(1.-XI3)*F11+F61*XI3 X2=(1.-XI3)*F12+F62*XI3 X3=(1.-XI3)*F13+F63*XI3 WRITE(99,*) Z1+N1*(Z2-1)+N1*N2*(Z3-1),X1,X2,X3 10 CONTINUE C** PRINT*,'written nodes : ',N1*N2*N3 C** C**----------------------------------------------------------------- C** C**** the generation of the elements : C** ------------------------------- C** C** The domain is covered by hexahedron elements of order 2 C** and consequently the boundaries are described by C** quadrilateral elements of order 2. The succesion of the C** nodes in the element is defined in vemu02 and vembuild(3). C** The lowest node id in an element is S. C** C** There are two different element types: C** WRITE(99,*) 2 C** C** These are the hexahedron 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,*) ELEM1*ELEM2*ELEM3,3,8,20 DO 20 Z3=1,ELEM3 DO 20 Z2=1,ELEM2 DO 20 Z1=1,ELEM1 S=2*(Z1-1)+2*(Z2-1)*N1+2*N1*N2*(Z3-1)+1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*(Z3-1) WRITE(99,'(22I5)') ELMID,IND, & S,S+2,S+2*N1+2,S+2*N1, & S+2*N1*N2,S+2*N1*N2+2,S+2*N1*N2+2*N1+2,S+2*N1*N2+2*N1, & S+1,S+N1+2,S+2*N1+1,S+N1, & S+N1*N2,S+N1*N2+2,S+N1*N2+2*N1+2,S+N1*N2+2*N1, & S+2*N1*N2+1,S+2*N1*N2+N1+2,S+2*N1*N2+2*N1+1,S+2*N1*N2+N1 20 CONTINUE C** PRINT*,'written hexahedron elements : ',ELEM1*ELEM2*ELEM3 C** C**-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- C** C** These are the quadrilateral elements: C** WRITE(99,*) 2*(ELEM1*ELEM2+ELEM1*ELEM3+ELEM2*ELEM3),2,4,8 C** DO 31 Z2=1,ELEM2 DO 31 Z1=1,ELEM1 C** C**** elements on boundary 1: C** S=2*(Z1-1)+2*(Z2-1)*N1+1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*ELEM3 WRITE(99,*) ELMID,IND, & S,S+2*N1,S+2*N1+2,S+2,S+N1,S+2*N1+1,S+N1+2,S+1 C** C**** elements on boundary 6: C** S=2*(Z1-1)+2*(Z2-1)*N1+2*N1*N2*ELEM3+1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2+ELEM1*ELEM2*ELEM3 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 31 CONTINUE DO 32 Z3=1,ELEM3 DO 32 Z2=1,ELEM2 C** C**** elements on boundary 5: C** S=2*(Z2-1)*N1+2*N1*N2*(Z3-1)+1 ELMID=Z2+ELEM2*(Z3-1)+2*ELEM1*ELEM2+ELEM1*ELEM2*ELEM3 WRITE(99,*) ELMID,IND, & S,S+2*N1*N2,S+2*N1*N2+2*N1,S+2*N1, & S+N1*N2,S+2*N1*N2+N1,S+N1*N2+2*N1,S+N1 C** C**** elements on boundary 3: C** S=2*ELEM1+2*(Z2-1)*N1+2*N1*N2*(Z3-1)+1 ELMID=Z2+ELEM2*(Z3-1)+ & ELEM2*ELEM3+2*ELEM1*ELEM2+ELEM1*ELEM2*ELEM3 WRITE(99,*) ELMID,IND, & S,S+2*N1,S+2*N1*N2+2*N1,S+2*N1*N2, & S+N1,S+N1*N2+2*N1,S+2*N1*N2+N1,S+N1*N2 32 CONTINUE DO 33 Z3=1,ELEM3 DO 33 Z1=1,ELEM1 C** C**** elements on boundary 2: C** S=2*(Z1-1)+2*N1*N2*(Z3-1)+1 ELMID=Z3+ELEM3*(Z1-1)+ & 2*ELEM2*ELEM3+2*ELEM1*ELEM2+ELEM1*ELEM2*ELEM3 WRITE(99,*) ELMID,IND, & S,S+2,S+2*N2*N1+2,S+2*N1*N2, & S+1,S+N1*N2+2,S+2*N2*N1+1,S+N1*N2 C** C**** elements on boundary 4: C** S=2*(Z1-1)+2*ELEM2*N1+2*N1*N2*(Z3-1)+1 ELMID=Z3+ELEM3*(Z1-1)+ & ELEM1*ELEM3+2*ELEM2*ELEM3+2*ELEM1*ELEM2+ELEM1*ELEM2*ELEM3 WRITE(99,*) ELMID,IND, & S,S+2*N1*N2,S+2*N1*N2+2,S+2, & S+N1*N2,S+2*N2*N1+1,S+N1*N2+2,S+1 33 CONTINUE C** PRINT*,'written quadrilateral elements : ', & 2*(ELEM1*ELEM2+ELEM1*ELEM3+ELEM2*ELEM3) C** C**----------------------------------------------------------------- C** C** generation of the nodes with Dirichlet conditions : C** ------------------------------------------------- C** C** Here we generate the Dirichlet conditions for 3D Navier- C** Stokes equation. The components 1,2,3 are the velocity C** compontents in x1,x2,x3 direction. The velocity is prescribed C** on the total boundary. The fourth component is the pressure, C** which is set on a single point. The real parameter is used C** to specify the velocity at the inlet (v1,v2,v3)=(val,0,0,0) C** and the outlet (v1,v2,v3)=(0,val,0), where val is the C** parabolic bubble over the sectional area. The program can be C** easily adapted to other nodes selections and other input or C** output profiles. C** WRITE(99,*) NK C** C** First the velocity components 1,2,3 : C** DO 40 D=1,NK-1 C** WRITE(99,*) 2*((N1-1)*(N2-1)+(N1-1)*(N3-1)+(N2-1)*(N3-1)+1) C** DO 41 Z2=1,N2-1 DO 41 Z1=1,N1-1 WRITE(99,*) Z1+N1*(Z2-1),ZERO WRITE(99,*) N1*N2*N3+1-Z1-N1*(Z2-1),ZERO 41 CONTINUE DO 42 Z3=1,N3-1 DO 42 Z1=1,N1-1 WRITE(99,*) N1*N2*(N3-1)-N1*N2*(Z3-1)+Z1,ZERO WRITE(99,*) N1*N2*Z3-Z1+1,ZERO 42 CONTINUE DO 43 Z3=1,N3-1 DO 43 Z2=1,N2-1 XI2=DBLE(Z2-1)/DBLE(N2-1) XI3=DBLE(Z3-1)/DBLE(N3-1) VALUE=XI2*(XI2-1)*XI3*(XI3-1)*16 IF (D.EQ.1) THEN WRITE(99,*) N1*N2*N3-N1*N2*(Z3-1)-N1*Z2+1,VALUE WRITE(99,*) N1*N2*(Z3-1)+N1*Z2,ZERO ELSEIF (D.EQ.3) THEN WRITE(99,*) N1*N2*N3-N1*N2*(Z3-1)-N1*Z2+1,ZERO WRITE(99,*) N1*N2*(Z3-1)+N1*Z2,-VALUE ELSE WRITE(99,*) N1*N2*N3-N1*N2*(Z3-1)-N1*Z2+1,ZERO WRITE(99,*) N1*N2*(Z3-1)+N1*Z2,ZERO ENDIF 43 CONTINUE WRITE(99,*) N1*(N2-1)+1,ZERO WRITE(99,*) N1*N2*(N3-1)+N1,ZERO C** PRINT*,2*((N1-1)*(N2-1)+(N1-1)*(N3-1)+(N2-1)*(N3-1)+1), & ' Dirichlet conditions for component ',d 40 CONTINUE C** C** Only one Dirichlet condition for the fourth component: C** WRITE(99,*) 1 WRITE(99,*) 2*((ELEM1+1)/2-1)+2*((ELEM2+1)/2-1)*N1+ & 2*N1*N2*((ELEM2+1)/2-1)+1,ZERO PRINT*,1,' Dirichlet condition for component ',NK C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C**