C******************************************************************* C** C** v e m e 0 2 e x m 0 4 : C** C** velocity driven diffusion on the 3-dimensional domain. For the C** solution the nonlinear solver is used. The mesh is read from C** an I-DEAS universal file. C** C** by L. Grosz Karlsruhe, Jan. 1995 C** C******************************************************************* C** C** The problem is the velocity driven diffusion problem, which C** is solved by the nonlinear solver veme02. The domain is an C** 3-dimensional arbitrary domain. An all boundaries Neuman C** boundary conditions are prescribed and on one point a C** Dirichlet condition is set. Using the notations in C** equation the problem is given by the functional equation: C** C** Dirichlet conditions: C** u1=b C** C** linear functional equation: F{u}(v)=0 C** C** with F{u}(v)= C** volume{v1x1 * u1x1 + v1x2 * u1x2 + v1x3 * u1x3 + C** v1 * ( w1 * u1x1 + w2 * u1x2 + w3 * u1x3 - f) } C** + area{v1 * g} C** C** The functions b, f and g are selected so that u1=x3 C** is the exact solution of this problem. We set w1=w2=0 C** and w3=16*x1*x2*(1-x1)*(1-x2). C** C** An example of an I-DEAS universal file you get in the data set C** cube.unv. The domain is the unit cube with center C** (0,0,0) and edge length 1. The mesh uses tetrahedron elements C** of order 2. cube.unv is the universal file. C** C**----------------------------------------------------------------- C** PROGRAM VEMEXM C** C**----------------------------------------------------------------- C** IMPLICIT NONE include 'bytes.h' C** C**----------------------------------------------------------------- C** C** some parameters which may be chanced: C** C** MAXNN = maximal number of nodes per process C** MAXNE = maximal number of elements per process C** MAXNG = maximal number of groups C** INPUT = name of the I-DEAS universal file C** STORE = total storage of process in Mbytes. C** INTEGER MAXNE,MAXNN,MAXNG,STORE CHARACTER*80 INPUT PARAMETER (MAXNN=4000, & MAXNE=2000, & MAXNG=5, & INPUT='cube.unv', & STORE=25) C** C**----------------------------------------------------------------- C** C** the parameters are explained in mesh. C** INTEGER NK,DIM,MESH,LOUT PARAMETER (NK=1, & DIM=3, & MESH=500, & LOUT=6) C** C**----------------------------------------------------------------- C** C** the length of the array for the mesh are set: C** INTEGER LU,LNODN,LNOD,LNOPRM,LNEK,LRPARM,LIPARM, & LDNOD,LIDPRM,LRDPRM,LIVEM,LRVEM,LLVEM,LBIG PARAMETER (LU =MAXNN*NK, & LNODN =MAXNN, & LNOD =MAXNN*DIM, & LNOPRM=1, & LNEK=20*MAXNE, & LIPARM=2*MAXNE, & LRPARM=1, & LDNOD =2*NK, & LIDPRM=NK, & LRDPRM=1, & LIVEM =MESH+800+LU+LDNOD/2, & LLVEM =500, & LRVEM =60+2*LU) C** C**----------------------------------------------------------------- C** C** RBIG should be as large as possible: the available C** storage STORE is reduced by all allocated array. C** the remaining storage is reserved for RBIG. C** PARAMETER ( LBIG=(STORE * 1 000 000)/IREAL & - (3*LU+LNOD+LNOPRM+LRPARM+LRDPRM) & - (LIVEM+LNODN+LNEK+LIPARM+LDNOD+LIDPRM)/RPI ) C** C**----------------------------------------------------------------- C** C** variables and arrays : C** -------------------- C** DOUBLE PRECISION T,NOD(LNOD),NOPARM(LNOPRM),RPARM(LRPARM), & RDPARM(LRDPRM),RBIG(LBIG),U(LU),RVEM(LRVEM), & EEST(LU),ERRG(LU),NRMERR(NK) INTEGER IVEM(LIVEM),NODNUM(LNODN),NEK(LNEK), & IPARM(LIPARM),DNOD(LDNOD),IDPARM(LIDPRM), & IBIG(RPI*LBIG) LOGICAL MASKL(NK,NK,MAXNG),MASKF(NK,MAXNG),LVEM(LLVEM) C*** INTEGER MYPROC,INFO,OUTFLG,GINFO,GINFO1,CLASS,NGROUP,G CHARACTER*80 NAME C*** EXTERNAL VEM630,VEM500 EXTERNAL DUMMY,USRFU,USERF,USERC,USERB C** C**----------------------------------------------------------------- C** C** The equivalence of RBIG and IBIG is very important : C** EQUIVALENCE (RBIG,IBIG) C** C**----------------------------------------------------------------- C** C** get task ids : C** NAME='a.out' CALL COMBGN(IVEM(200),MYPROC,LIVEM-203,IVEM(204),NAME,INFO) IF (INFO.NE.0) GOTO 9999 IVEM(201)=MYPROC IVEM(202)=0 IVEM(203)=IVEM(204) C** C**----------------------------------------------------------------- C** C** a protocol is printed only on process 1 : C** IF (MYPROC.EQ.1) THEN OUTFLG=1 ELSE OUTFLG=0 ENDIF C** C**----------------------------------------------------------------- C** C**** the parameters are copied into IVEM : C** ----------------------------------- C** IVEM(1)=MESH IVEM(MESH+ 2)=NK IVEM(MESH+ 3)=DIM C** C**----------------------------------------------------------------- C** C**** read the mesh from a universal file : C** ----------------------------------- C** IVEM(120)=LOUT IVEM(121)=OUTFLG IVEM(122)=9 IVEM(124)=1 IF (MYPROC.EQ.1) OPEN(9,FILE=INPUT,STATUS= 'UNKNOWN', & FORM='FORMATTED') CALL IDEVEM (LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG) CLOSE(9) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** print mesh on processor 1 C** ------------------------- C** IVEM(20)=LOUT IVEM(21)=0000*OUTFLG IVEM(22)=2 CALL VEMU01(LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM, & LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** distribute mesh : C** ---------------- C** IVEM(80)=LOUT IVEM(81)=OUTFLG IVEM(51)=2 CALL VEMDIS (LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM, & LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM, & LBIG,RBIG,IBIG) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** set masks : C** --------- C** NGROUP=IVEM(IVEM(1)+4) GINFO=IVEM(IVEM(1)+21)+IVEM(1) GINFO1=IVEM(IVEM(1)+22) DO 400 G=1,NGROUP MASKF(1,G)=.FALSE. MASKL(1,1,G)=.FALSE. CLASS=IVEM(GINFO+GINFO1*(G-1)+3) IF (CLASS.EQ.DIM) THEN MASKF(1,G)=.TRUE. MASKL(1,1,G)=.TRUE. ELSEIF (CLASS.EQ.DIM-1) THEN MASKF(1,G)=.TRUE. ENDIF 400 CONTINUE C** C**----------------------------------------------------------------- C** C**** call of VECFEM : C** -------------- C** OPEN(10,FORM='UNFORMATTED',STATUS='SCRATCH') OPEN(11,FORM='UNFORMATTED',STATUS='SCRATCH') OPEN(12,FORM='UNFORMATTED',STATUS='SCRATCH') LVEM(1)=.FALSE. LVEM(4)=.FALSE. LVEM(5)=.FALSE. LVEM(6)=.TRUE. LVEM(7)=.TRUE. LVEM(8)=.TRUE. LVEM(9)=.FALSE. LVEM(10)=.TRUE. LVEM(11)=.FALSE. RVEM(1)=0 RVEM(3)=1.D-2 RVEM(10)=1.D-8 IVEM(3)=0 IVEM(10)=10 IVEM(11)=11 IVEM(12)=12 IVEM(40)=LOUT IVEM(41)=50*OUTFLG IVEM(45)=500 IVEM(46)=0 IVEM(60)=0 IVEM(70)=10 IVEM(71)=11 IVEM(72)=10 000 CALL VEME02 (T,LU,U,EEST,LIVEM,IVEM,LLVEM,LVEM,LRVEM,RVEM, & LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG, & MASKL,MASKF,USERB,USRFU,USERF,VEM500,VEM630) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** compute the error on the geometrical nodes : C** ------------------------------------------ C** IVEM(4)=30 IVEM(30)=LOUT IVEM(31)=OUTFLG*0 IVEM(32)=MAXNN IVEM(33)=NK CALL VEMU05 (T,LU,ERRG,LU,U,LIVEM,IVEM, & LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG, & USERC) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** print the error and its norm : C** ---------------------------- C** IVEM(23)=LOUT IVEM(24)=OUTFLG IVEM(25)=IVEM(32) IVEM(26)=IVEM(33) CALL VEMU13 (LU,ERRG,NRMERR,LIVEM,IVEM, & LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** 9999 CALL COMEND(IVEM(200),INFO) E N D SUBROUTINE USERB(T,COMPU,RHS, & NRSDP,RSDPRM,NRVDP,RVDP1,RVDPRM, & NISDP,ISDPRM,NIVDP,IVDP1,IVDPRM, & NDC,DIM,X,NOP,NOPARM,B) C** C******************************************************************* C** C** the routine USERB sets the Dirichlet boundary conditions, C** see userb. here the exact solution x3 is prescribed. C** C******************************************************************* C** INTEGER COMPU,RHS,NRSDP,NRVDP,RVDP1,NISDP,NIVDP,IVDP1, & NDC,DIM,NOP DOUBLE PRECISION T,RSDPRM(NRSDP),RVDPRM(RVDP1,NRVDP), & X(NDC,DIM),NOPARM(NDC,NOP),B(NDC) INTEGER ISDPRM(NISDP),IVDPRM(IVDP1,NIVDP) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** IF (COMPU.EQ.1) THEN DO 10 Z=1,NDC B(Z) = X(Z,3) 10 CONTINUE ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERB-------------------------------------------------- E N D SUBROUTINE USRFU(T,GROUP,CLASS,COMPV,COMPU,LAST, & NELIS,L,DIM,X,TAU,NK,U,DUDX, & LT,UT,DUTDX,NOP,NOPARM,DNOPDX, & NRSP,RSPARM,NRVP,RVP1,RVPARM, & NISP,ISPARM,NIVP,IVP1,IVPARM, & F1UX,F1U,F0UX,F0U) C** C******************************************************************* C** C** the routine USRFU sets the Frechet derivative of the linear C** form F, see usrfu: C** C******************************************************************* C** INTEGER GROUP,CLASS,COMPV,COMPU,LAST,NELIS,L,LT, & DIM,NK,NOP,NRSP,RVP1,NRVP,NISP,IVP1,NIVP DOUBLE PRECISION T,X(L,DIM),TAU(L,DIM,CLASS),U(L,NK),UT(LT,NK), & DUDX(L,NK,CLASS),DUTDX(LT,NK,CLASS), & NOPARM(L,NOP),DNOPDX(L,NOP,CLASS), & RSPARM(NRSP),RVPARM(RVP1,NRVP), & F1UX(L,CLASS,CLASS),F1U(L,CLASS),F0UX(L,CLASS), & F0U(L) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** --------------------- C** C** the coefficients for the area integration : C** IF (CLASS.EQ.3) THEN IF ((COMPV.EQ.1).AND.(COMPU.EQ.1)) THEN DO 112 Z=1,NELIS F1UX(Z,1,1)=1. F1UX(Z,2,2)=1. F1UX(Z,3,3)=1. F0UX(Z,3)=16*X(Z,1)*X(Z,2)*(1.-X(Z,1))*(1.-X(Z,2)) 112 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USRFU-------------------------------------------------- E N D SUBROUTINE USERF (T,GROUP,CLASS,COMPV,RHS,LAST, & NELIS,L,DIM,X,TAU,NK,U,DUDX, & LT,UT,DUTDX,NOP,NOPARM,DNOPDX, & NRSP,RSPARM,NRVP,RVP1,RVPARM, & NISP,ISPARM,NIVP,IVP1,IVPARM, & F1,F0) C** C******************************************************************* C** C** the routine USERF sets the coefficients of the linear form F, C** see userf: C** C** It is f=-16*x1*x2*(1-x1)*(1-x2) * u1x3 C** g=-(n1*u1x1+n2*u1x1+n3*u1x3)=-n3 C** C** computed by the exact solution u1=x3. (n1,n2,n3) is the outer C** normal field, which is computed by the tangential field: C** n3=(tau11*tau23-tau21*tau12)/norm , norm is the normalization. C** C******************************************************************* C** INTEGER GROUP,CLASS,COMPV,RHS,LAST,NELIS,L,LT,DIM,NK,NOP, & NRSP,RVP1,NRVP,NISP,IVP1,NIVP DOUBLE PRECISION T,X(L,DIM),TAU(L,DIM,CLASS),U(L,NK),UT(LT,NK), & DUDX(L,NK,CLASS),DUTDX(LT,NK,CLASS), & NOPARM(L,NOP),DNOPDX(L,NOP,CLASS), & RSPARM(NRSP),RVPARM(RVP1,NRVP), & F1(L,CLASS),F0(L) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z DOUBLE PRECISION NORM C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** C** the coefficients for the volume integration : C** IF (CLASS.EQ.3) THEN IF (COMPV.EQ.1) THEN DO 12 Z=1,NELIS F1(Z,1)=DUDX(Z,1,1) F1(Z,2)=DUDX(Z,1,2) F1(Z,3)=DUDX(Z,1,3) F0(Z)=16*X(Z,1)*X(Z,2)*(1.-X(Z,1))*(1.-X(Z,2))* & (DUDX(Z,1,3)-1.) 12 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C** the coefficients for the area integration : C** the local numbering of the surface elements in I-DEAS C** makes it necessary to change the sign of one tangential C** direction to get the outer normal field. C** IF (CLASS.EQ.2) THEN IF (COMPV.EQ.1) THEN DO 11 Z=1,NELIS NORM = (TAU(Z,2,1)*TAU(Z,3,2)-TAU(Z,3,1)*TAU(Z,2,2))**2 & + (TAU(Z,1,1)*TAU(Z,3,2)-TAU(Z,3,1)*TAU(Z,1,2))**2 & + (TAU(Z,1,1)*TAU(Z,2,2)-TAU(Z,2,1)*TAU(Z,1,2))**2 F0(Z)= (TAU(Z,1,1)*TAU(Z,2,2)-TAU(Z,2,1)*TAU(Z,1,2)) & /SQRT(NORM) 11 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERF------------------------------------------------- E N D SUBROUTINE USERC(T,GROUP,LAST,NELIS, & NRSP,RSPARM,NRVP,RVP1,RVPARM, & NISP,ISPARM,NIVP,IVP1,IVPARM, & L,DIM,X,NK,U,DUDX,NOP,NOPARM,DNOPDX,N,CU) C** C******************************************************************* C** C** the routine USERC computes in this case the error of the C** computed solution, see userc. C** C******************************************************************* C** INTEGER GROUP,LAST,NELIS,L,DIM,NK,N, & NRSP,RVP1,NRVP,NISP,IVP1,NIVP,NOP DOUBLE PRECISION T,X(L,DIM),U(L,NK),DUDX(L,NK,DIM), & RSPARM(NRSP),RVPARM(RVP1,NRVP), & NOPARM(L,NOP),DNOPDX(L,NOP,DIM),CU(L,N) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** DO 10 Z=1,NELIS CU(Z,1) = ABS( U(Z,1) - X(Z,3) ) 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERC-------------------------------------------------- E N D