! Program to test integration of Tetra1 ! ------------------------------------- ! We compute numerically the integrals of some functions for which we can ! get the exact solution. ! From the values of those functions and their first derivatives (v) at the ! nodes of a random tetrahedron (coor), we compute the integrals by using ! the shape functions (phi) defined at a set of integration points (point) ! and the related weights (weight). ! This does only check the values of the shape functions (stored into phi(1,...)), ! not their derivatives stored into phi(2-4,...). IMPLICIT INTEGER*4 (i-n), REAL*8 (a-h,o-z) DIMENSION coor(3,4),point(3,14),phi(4,16,14),weight(14),v(4),cu(3,4),cu1(3,4) REAL*4 xx DATA cu /0.d0,0.d0,0.d0, 1.d0,0.d0,0.d0, 0.d0,1.d0,0.d0, 0.d0,0.d0,1.d0/ DATA cu1/1.d0,0.d0,0.d0, 1.d0,1.d0,0.d0, 0.d0,1.d0,0.d0, 0.d0,0.d0,1.d0/ OPEN (UNIT=10,file='resu1.txt',STATUS='UNKNOWN') WRITE (10,1) 1 FORMAT (' Hermite (Délèze) tetrahedron interpolation with 16 parameters'/ & ' Check integration of quadratic terms.'/) kcas = 1 2 IF (kcas.EQ.1) THEN coor = cu WRITE (10,3) kcas 3 FORMAT (/' Case ',i2,': Unit tetrahedron'/ & ' =========================') ELSEIF (kcas.EQ.2) THEN coor = cu1 WRITE (10,4) kcas 4 FORMAT (/' Case ',i2,': Unit rotated tetrahedron'/ & ' ================================') ELSEIF (kcas.EQ.3) THEN iran = 1. DO i = 1,3 DO j = 1,4 xx = RAND (iran); iran = 0 coor(i,j) = xx*10.d0 ENDDO ENDDO WRITE (10,5) kcas 5 FORMAT (/' Case ',i2,': Random tetrahedron'/ & ' ===========================') ENDIF x1 = coor(1,1); y1 = coor(2,1); z1 = coor(3,1) x2 = coor(1,2); y2 = coor(2,2); z2 = coor(3,2) x3 = coor(1,3); y3 = coor(2,3); z3 = coor(3,3) x4 = coor(1,4); y4 = coor(2,4); z4 = coor(3,4) WRITE (10,'('' x1,y1,z1 = '',1P3D15.6)') x1,y1,z1 WRITE (10,'('' x2,y2,z2 = '',1P3D15.6)') x2,y2,z2 WRITE (10,'('' x3,y3,z3 = '',1P3D15.6)') x3,y3,z3 WRITE (10,'('' x4,y4,z4 = '',1P3D15.6)') x4,y4,z4 ! volume vol = ABS ( (x2-x1)*((y3-y1)*(z4-z1)-(y4-y1)*(z3-z1)) + & (x3-x1)*((y4-y1)*(z2-z1)-(y2-y1)*(z4-z1)) + & (x4-x1)*((y2-y1)*(z3-z1)-(y3-y1)*(z2-z1))) vol = vol/6.d0 WRITE (10,'('' Volume = '',1PD15.6)') vol ! centroid xm = (coor(1,1)+coor(1,2)+coor(1,3)+coor(1,4))/4.d0 ym = (coor(2,1)+coor(2,2)+coor(2,3)+coor(2,4))/4.d0 zm = (coor(3,1)+coor(3,2)+coor(3,3)+coor(3,4))/4.d0 ! shape functions npi = 5 ivers = 1 10 WRITE (10,12) npi 12 FORMAT (/' Nb. of integration points = ',i3) CALL FElib_Tetra_Shape (coor,npi,point,ivers,phi,weight) ! check integration of u(x,y,z) = 1 sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = 1.d0; v(2) = 0.d0; v(3) = 0.d0; v(4) = 0.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = vol WRITE (10,*) 'Int (1) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = x-xm sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = coor(1,j)-xm; v(2) = 1.d0; v(3) = 0.d0; v(4) = 0.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 WRITE (10,*) 'Int (x-xm) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = y-ym sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = coor(2,j)-ym; v(2) = 0.d0; v(3) = 1.d0; v(4) = 0.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 WRITE (10,*) 'Int (y-ym) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = z-zm sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = coor(3,j)-zm; v(2) = 0.d0; v(3) = 0.d0; v(4) = 1.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 WRITE (10,*) 'Int (z-zm) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = (x-xm)**2 sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = (coor(1,j)-xm)**2; v(2) = 2*(coor(1,j)-xm); v(3) = 0.d0 v(4) = 0.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 DO i = 1,4 sume = sume+(coor(1,i)-xm)**2 ENDDO sume = sume*vol/20.d0 WRITE (10,*) 'Int ((x-xm)**2) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = (y-ym)**2 sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = (coor(2,j)-ym)**2; v(2) = 0.d0; v(3) = 2*(coor(2,j)-ym); v(4) = 0.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 DO i = 1,4 sume = sume+(coor(2,i)-ym)**2 ENDDO sume = sume*vol/20.d0 WRITE (10,*) 'Int ((y-ym)**2) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = (z-zm)**2 sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = (coor(3,j)-zm)**2; v(2) = 0.d0; v(3) = 0.d0 v(4) = 2*(coor(3,j)-zm); j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 DO i = 1,4 sume = sume+(coor(3,i)-zm)**2 ENDDO sume = sume*vol/20.d0 WRITE (10,*) 'Int ((z-zm)**2) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = (x-xm)*(y-ym) sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = (coor(1,j)-xm)*(coor(2,j)-ym); v(2) = (coor(2,j)-ym) v(3) = (coor(1,j)-xm); v(4) = 0.d0 j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 DO i = 1,4 sume = sume+(coor(1,i)-xm)*(coor(2,i)-ym) ENDDO sume = sume*vol/20.d0 WRITE (10,*) 'Int ((x-xm)*(y-ym)) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = (x-xm)*(z-zm) sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = (coor(1,j)-xm)*(coor(3,j)-zm); v(2) = (coor(3,j)-zm) v(3) = 0.d0; v(4) = (coor(1,j)-xm) j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 DO i = 1,4 sume = sume+(coor(1,i)-xm)*(coor(3,i)-zm) ENDDO sume = sume*vol/20.d0 WRITE (10,*) 'Int ((x-xm)*(z-zm)) = ',sum WRITE (10,*) 'exact = ',sume ! check integration of u(x,y,z) = (y-ym)*(z-zm) sum = 0.d0 DO i = 1,npi DO j = 1,4 v(1) = (coor(2,j)-ym)*(coor(3,j)-zm); v(2) = 0.d0 v(3) = (coor(3,j)-zm); v(4) = (coor(2,j)-ym) j1 = 4*(j-1)+1 sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ & phi(1,j1+3,i)*v(4))*weight(i) ENDDO ENDDO sume = 0.d0 DO i = 1,4 sume = sume+(coor(2,i)-ym)*(coor(3,i)-zm) ENDDO sume = sume*vol/20.d0 WRITE (10,*) 'Int ((y-xm)*(z-ym)) = ',sum WRITE (10,*) 'exact = ',sume ! next integration scheme IF (npi.EQ.5) THEN npi = 10; GOTO 10 ELSEIF (npi.EQ.10) THEN npi = 11; GOTO 10 ELSEIF (npi.EQ.11) THEN npi = 14; GOTO 10 ENDIF ! next tetra IF (kcas.EQ.1) THEN kcas = 2; GOTO 2 ELSEIF (kcas.EQ.2) THEN kcas = 3; GOTO 2 ENDIF CLOSE (UNIT=10) END