C******************************************************************* C** C** v e m e 0 2 e x m 0 5 C** C** thermal diffusion with temperature-dependent material C** coefficients on a 2-dimensional ring segment. The C** isoparametrical mesh is generated distributed on the C** processors. C** C** by L. Grosz Karlsruhe, Jan. 1995 C** C******************************************************************* C** C** The problem is equation of the Poisson type with a diffusion, C** coefficient, which depends on the solution. The equation C** is solved by the nonlinear solver veme02. The domain C** is a quarter segment of a ring with inner radius 1 and outer C** radius 2. On the arc of the circular with radius 1 and the C** x1-axis (boundaries 1 and 2) Neuman boundary conditions are C** prescribed and on the arc of the circular with radius 2 and C** the x2-axis (boundaries 3 and 4) Dirichlet conditions C** are set. Using the notations in equation the problem is C** given by the functional equation: C** C** Dirichlet conditions: C** u1=b C** C** nonlinear functional Equation: F{u}(v)=0 C** C** with F{u}(v)= C** area{ (1+theta*u1)*(v1x1 * u1x1 + v1x2 * u1x2)+ f*v1 } C** + line{v1 * g} C** C** where theta is a given real constant. The functions b, f and C** g are selected so u1=cos(alpha*(x1^2+x2)^2)*x1 is the C** exact solution of the problem (alpha is a real constants). C** For fixed finite element mesh the quality of the FEM C** approximation will become worse for increasing alpha. C** C** The ring segment is subdivided into quadrilateral elements of C** order 2. Therefore the boundary is subdivided into line C** elements of order 2. The mesh is generated distributed onto C** the processes. The error of the computed solution C** approximation is calculated. C** PROGRAM VEMEXM C** C**----------------------------------------------------------------- C** IMPLICIT NONE include 'bytes.h' C** C**----------------------------------------------------------------- C** C** some parameters which may be chanced: C** C** ELEM1 = number of elements in angle direction, C** in radial direction .64*ELEM1 elements will be C** generated to get nearly quadratical elements, but C** only about ELEM1/NPROC on this process. C** STORE = total storage of process in Mbytes. C** INTEGER NPROC,ELEM1,STORE PARAMETER (NPROC=1, & ELEM1=40, & STORE=50) C** C**----------------------------------------------------------------- C** C** ELEM2 = number of elements in radial direction on process C** N1,N2 = number of nodes in angle/radial-direction on the C** process C** C** other parameters are explained in mesh. C** INTEGER N1,N2,NK,NGROUP,DIM,MESH,GINFO,GINFO1,DINFO,DINFO1, & ELEM2,LOUT PARAMETER (ELEM2=(0.64*ELEM1+NPROC-1)/NPROC, & N1=2*ELEM1+1, & N2=2*ELEM2+1, & NK=1, & NGROUP=2, & DIM=2, & MESH =210+NPROC, & GINFO =30, & GINFO1=23+2*NK, & DINFO =GINFO+GINFO1*NGROUP, & DINFO1=17, & LOUT=6) C** C**----------------------------------------------------------------- C** C** the length of the array for the mesh are set: C** they are a little bit greater than actual used in the C** mesh generation. this is necessary for the mesh distribution. C** INTEGER NN,LU,LNODN,LNOD,LNOPRM,LNEK,LRPARM,LIPARM, & LDNOD,LIDPRM,LRDPRM,LIVEM,LRVEM,LLVEM,LBIG PARAMETER (NN=N1*N2*1.5, & LU =NN*NK, & LNODN =NN, & LNOD =NN*DIM, & LNOPRM=1, & LNEK=(8*(ELEM1*ELEM2+1)+6*(ELEM1+ELEM2+1))*3, & LIPARM=(ELEM1*ELEM2+2*(ELEM1+ELEM2)+2)*1.5, & LRPARM=1, & LDNOD =4*(N1+N2)*NK*1.5, & LIDPRM=2*(N1+N2)*NK*1.5, & LRDPRM=1, & LIVEM =MESH+DINFO+DINFO1*NK+600+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,NGROUP),MASKF(NK,NGROUP),LVEM(LLVEM) C*** INTEGER MYPROC,INFO,OUTFLG,NDNUM0,HERE,S,NE1,ADGEO1, & NE2,ADGEO2,ADIVP2,ADIVP1,NE0,NDC1,ADDCG1,ELNUM0 INTEGER Z1,Z2 DOUBLE PRECISION R0,PIH,THETA,PHI,R,ALPHA CHARACTER*80 NAME C*** EXTERNAL VEM630,VEM500 EXTERNAL DUMMY,USERB,USRFU,USERF,USERC C** C**----------------------------------------------------------------- C** C** The equivalence of RBIG and IBIG is very important : C** EQUIVALENCE (RBIG,IBIG) C** C**----------------------------------------------------------------- C** C** The COMMON block is used to transport the constants into the C** subroutines USERB, USERF, USRFU, USERC: C** COMMON /PROP/THETA,ALPHA C** THETA=.3 ALPHA=1 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) IF (NPROC.NE.IVEM(200)) THEN PRINT*,'Set NPROC=',IVEM(200) GOTO 9999 ENDIF 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+ 1)=N1*N2 IVEM(MESH+ 2)=NK IVEM(MESH+ 3)=DIM IVEM(MESH+ 4)=NGROUP IVEM(MESH+ 5)=NN IVEM(MESH+13)=NN IVEM(MESH+14)=0 IVEM(MESH+15)=0 IVEM(MESH+18)=0 IVEM(MESH+21)=GINFO IVEM(MESH+22)=GINFO1 IVEM(MESH+23)=DINFO IVEM(MESH+24)=DINFO1 C** C**----------------------------------------------------------------- C** C** This process generates the nodes in the subdomain with a C** distance between R0 and R0+1/NPROC from the origin starting C** with the node id number NDNUM0. The nodes with distance R0 C** are also generated on process MYPROC-1 and the nodes with C** distance R0+1/NPROC are also generated on process MYPROC+1. C** The first element generated on the process gets the element C** id number ELNUM0. C** R0=1.+DBLE(MYPROC-1)/DBLE(IVEM(200)) NDNUM0=(MYPROC-1)*N1*(N2-1)+1 ELNUM0=(MYPROC-1)*(ELEM1*ELEM2+ELEM1+2*ELEM2)+1 C** C**----------------------------------------------------------------- C** C**** the generation of the geometrical nodes : C** --------------------------------------- C** PIH=ATAN(1.)*2 DO 10 Z2=1,N2 DO 10 Z1=1,N1 PHI=1-DBLE(Z1-1)/DBLE(N1-1) R=R0+DBLE(Z2-1)/DBLE(IVEM(200)*(N2-1)) NOD(Z1+N1*(Z2-1) )=R*COS(PIH*PHI) NOD(Z1+N1*(Z2-1)+NN)=R*SIN(PIH*PHI) NODNUM(Z1+N1*(Z2-1) )=Z1+N1*(Z2-1)+NDNUM0-1 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** the generation of the elements : C** ------------------------------- C** C** The domain is covered by quadrilateral elements of order 2 C** and consequently the boundary is described by line elments C** of order 2. The following picture illustrates the C** construction of the quadrilateral element with lower S C** and its boundary elements which are only generated C** if they are a subset of the boundaries 1 or 2: C** C** S+2*N1 S+2*N1+1 S+2*N1+2 C** 4-----7-----3 2 C** | | | C** S+N1 8 6 3 S+N1+2 C** | | | C** 1-----5-----4 1 C** 1-----3-----2 C** S S+1 S+2 C** C** ADGEO1 defines the start address of the quadrilateral C** elements in NEK and ADIVP1 defines the start address of C** the element id number assigned to the elements. The element C** id number is unique over all processes. NE1 is the number of C** quadrilateral elements generated on the process. HERE gives C** the address of the element in NEK, which lowest vertex has C** the node id S over all processes. C** ADGEO1=1 ADIVP1=1 NE1=ELEM1*ELEM2 DO 20 Z2=1,ELEM2 DO 20 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z2-1)+ADGEO1-1 S =2*(Z1-1) + 2*(Z2-1) * N1 + NDNUM0 IPARM(ADIVP1-1+Z1+ELEM1*(Z2-1))=Z1+ELEM1*(Z2-1)+ ELNUM0-1 NEK(HERE )=S NEK(HERE+ NE1)=S + 2 NEK(HERE+2*NE1)=S + 2 *N1 + 2 NEK(HERE+3*NE1)=S + 2 *N1 NEK(HERE+4*NE1)=S + 1 NEK(HERE+5*NE1)=S + N1 + 2 NEK(HERE+6*NE1)=S + 2 *N1 + 1 NEK(HERE+7*NE1)=S + N1 20 CONTINUE C** C** ADGEO2 defines the start address of the line elements C** in NEK and ADIVP2 defines the start address of the C** element id number assigned to the elements. The entries 1 to C** 8*NE1 in NEK and 1 to NE1 in IPARM are already used by C** the elements in group 1. NE2 is the number of C** line elements generated on the process, where the C** elements on boundary 1 are only generated on process 1. C** HERE gives the address of the element in NEK, which is C** a boundary element of the quadrilateral element with lowest C** node id S. C** ADGEO2=ADGEO1+8*NE1 ADIVP2=ADIVP1+NE1 NE2=ELEM2 IF (MYPROC.EQ.1) NE2=NE2+ELEM1 C** C** these are the line elements on boundary 1: C** only on process 1. NE0 counts the already generated line C** elements C** NE0=0 IF (MYPROC.EQ.1) THEN DO 31 Z1=1,ELEM1 HERE=Z1+NE0+ADGEO2-1 S =2*(Z1-1) + NDNUM0 IPARM(ADIVP2-1+Z1+NE0)=Z1+NE0+ELEM1*ELEM2+ELNUM0-1 NEK(HERE )=S NEK(HERE+ NE2)=S + 2 NEK(HERE+2*NE2)=S + 1 31 CONTINUE NE0=NE0+ELEM1 ENDIF C** C** these are the line elements on boundary 2: C** DO 32 Z2=1,ELEM2 HERE=Z2+NE0+ADGEO2-1 S = 2*(ELEM1-1) + 2*(Z2-1) * N1 + NDNUM0 IPARM(ADIVP2-1+Z2+NE0)=Z2+NE0+ELEM1*ELEM2+ELNUM0-1 NEK(HERE )=S + 2 NEK(HERE+ NE2)=S + 2 + N1*2 NEK(HERE+2*NE2)=S + 2 + N1 32 CONTINUE NE0=NE0+ELEM2 C** C**----------------------------------------------------------------- C** C** the start addresses, etc are written to IVEM: C** C** group 1: quadrilateral elements C** IVEM(MESH+GINFO ) = NE1 IVEM(MESH+GINFO+ 2) = 4 IVEM(MESH+GINFO+ 3) = 2 IVEM(MESH+GINFO+ 8) = 0 IVEM(MESH+GINFO+11) = 0 IVEM(MESH+GINFO+13) = 0 IVEM(MESH+GINFO+14) = ADIVP1 IVEM(MESH+GINFO+15) = NE1 IVEM(MESH+GINFO+16) = 1 IVEM(MESH+GINFO+20) = ADGEO1 IVEM(MESH+GINFO+21) = NE1 IVEM(MESH+GINFO+23) = 8 C** C** group 2: line elements C** IVEM(MESH+GINFO+GINFO1 ) = NE2 IVEM(MESH+GINFO+GINFO1+ 2) = 2 IVEM(MESH+GINFO+GINFO1+ 3) = 1 IVEM(MESH+GINFO+GINFO1+ 8) = 0 IVEM(MESH+GINFO+GINFO1+11) = 0 IVEM(MESH+GINFO+GINFO1+13) = 0 IVEM(MESH+GINFO+GINFO1+14) = ADIVP2 IVEM(MESH+GINFO+GINFO1+15) = NE2 IVEM(MESH+GINFO+GINFO1+16) = 1 IVEM(MESH+GINFO+GINFO1+20) = ADGEO2 IVEM(MESH+GINFO+GINFO1+21) = NE2 IVEM(MESH+GINFO+GINFO1+23) = 3 C** C**----------------------------------------------------------------- C** C** generation of the nodes with Dirichlet conditions : C** ------------------------------------------------- C** C** NDC1 counts the Dirichlet condition generated for C** component 1. ADDCGC gives the start address of C** node id in DNOD and the vector integer parameter in IPARM, C** which is the id number of the boundary of the node. C** NDC1=0 ADDCG1=1 C** C** for component 1 the boundary 4 gets Dirichlet conditions: C** DO 43 Z2=1,N2 DNOD (ADDCG1+NDC1-1+Z2)= (Z2-1)*N1+NDNUM0 IDPARM(ADDCG1+NDC1-1+Z2)= 4 43 CONTINUE NDC1=NDC1+N2 C** C** for component 1 the boundary 3 gets Dirichlet conditions: C** (only on processor 1) C** IF (MYPROC.EQ.IVEM(200)) THEN DO 44 Z1=1,N1 DNOD (ADDCG1+NDC1-1+Z1)=N1*(N2-1)+Z1+NDNUM0-1 IDPARM(ADDCG1+NDC1-1+Z1)=3 44 CONTINUE NDC1=NDC1+N1 ENDIF C** C**----------------------------------------------------------------- C** C** the start addresses, etc are written to IVEM: C** C** component 1: C** IVEM(MESH+DINFO ) = NDC1 IVEM(MESH+DINFO+ 2) = ADDCG1 IVEM(MESH+DINFO+ 4) = 0 IVEM(MESH+DINFO+ 7) = 0 IVEM(MESH+DINFO+ 9) = 0 IVEM(MESH+DINFO+10) = ADDCG1 IVEM(MESH+DINFO+11) = NDC1 IVEM(MESH+DINFO+12) = 1 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**** 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(3)=1.D-3 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)=1 IVEM(72)=10 000 MASKL(1,1,1)=.TRUE. MASKL(1,1,2)=.FALSE. MASKF(1,1)=.TRUE. MASKF(1,2)=.TRUE. 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 IVEM(32)=NN 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. C** C******************************************************************* C** IMPLICIT NONE 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 DOUBLE PRECISION THETA,ALPHA COMMON /PROP/THETA,ALPHA C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** IF (COMPU.EQ.1) THEN DO 10 Z=1,NDC B(Z) = COS(ALPHA*(X(Z,1)**2+X(Z,2)**2))*X(Z,1) 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 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, C** see userf: C** C******************************************************************* C** IMPLICIT NONE 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 THETA,ALPHA,U0,H1,U0X1,U0X2,K COMMON /PROP/THETA,ALPHA C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** C** the coefficients for the area integration : C** IF (CLASS.EQ.2) THEN IF (COMPV.EQ.1) THEN DO 12 Z=1,NELIS K=1+THETA*U(Z,1) F1(Z,1)=K*DUDX(Z,1,1) F1(Z,2)=K*DUDX(Z,1,2) H1=ALPHA*(X(Z,1)**2+X(Z,2)**2) U0=COS(H1)*X(Z,1) K=1+THETA*U0 F0(Z)=(-2*THETA*SIN(H1)*ALPHA*X(Z,1)**2+THETA*COS(H1))* & (-2*SIN(H1)*ALPHA*X(Z,1)**2+COS(H1))+ & K*(-4*COS(H1)*ALPHA**2*X(Z,1)**3-6*SIN(H1)*ALPHA*X(Z,1))+ & 4*THETA*SIN(H1)**2*ALPHA**2*X(Z,2)**2*X(Z,1)**2- & 4*K*COS(H1)*ALPHA**2*X(Z,2)**2*X(Z,1)- & 2*K*SIN(H1)*ALPHA*X(Z,1) 12 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C** the coefficients for the line integration : C** IF (CLASS.EQ.1) THEN IF (COMPV.EQ.1) THEN DO 11 Z=1,NELIS H1=ALPHA*(X(Z,1)**2+X(Z,2)**2) U0X1=-2*ALPHA*SIN(H1)*X(Z,1)**2+COS(H1) U0X2=-2*ALPHA*SIN(H1)*X(Z,1)*X(Z,2) K=1+THETA*COS(H1)*X(Z,1) F0(Z)=-K*(U0X1*TAU(Z,2,1)-U0X2*TAU(Z,1,1)) 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 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 F with C** respect of u, see usrfu: C** C******************************************************************* C** IMPLICIT NONE 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 DOUBLE PRECISION THETA,ALPHA,K COMMON /PROP/THETA,ALPHA C** C**----------------------------------------------------------------- C** C**** start of calculation : C** --------------------- C** C** the coefficients for the area integration : C** IF (CLASS.EQ.2) THEN IF ((COMPV.EQ.1).AND.(COMPU.EQ.1)) THEN DO 112 Z=1,NELIS K=1+THETA*U(Z,1) F1UX(Z,1,1)=K F1UX(Z,2,2)=K F1U(Z,1)=THETA*DUDX(Z,1,1) F1U(Z,2)=THETA*DUDX(Z,1,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 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** IMPLICIT NONE 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 DOUBLE PRECISION THETA,ALPHA COMMON /PROP/THETA,ALPHA C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** DO 10 Z=1,NELIS CU(Z,1) = ABS( U(Z,1) - & COS(ALPHA*(X(Z,1)**2+X(Z,2)**2))*X(Z,1)) 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERC-------------------------------------------------- E N D