Merge branch 'compute-area-nek-expansion' into 'develop'

Compute area nek expansion

See merge request !18

examples/integration_example/integration_example.usr

@@ -20,14 +20,16 @@ c-----------------------------------------------------------------------
       return
       end
 c-----------------------------------------------------------------------
+! Sets the force term in the momentum equation. Use this to set body
+! forces such as gravity. Note that ffx, ffy, and ffz will later be
+! multiplied by the density.
       subroutine userf  (ix,iy,iz,eg)
       include 'SIZE'
       include 'TOTAL'
       include 'NEKUSE'
-      integer e,f,eg
       ffx=0.0
       ffy=0.0
-      ffz=0.0     
+      ffz=0.0
       return
       end
 c-----------------------------------------------------------------------
@@ -305,42 +307,45 @@ C=======================================================================
       parameter (nf_max=100)
 
       include 'SIZE'
-      include 'TOTAL'       
+      include 'TOTAL'
       common/expansion_tdata/n_legendre, m_fourier
-      common/expansion_tcoef/coeff_tij(nl_max,nf_max)  
+      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
       integer e,f
 
-      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt) 
-      real*8 sint, sint1
+      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
+      real*8 sint, sint1, sarea, sarea1
+      real*8 pi
 
+      pi=4.0*atan(1.0)
       ntot=nx1*ny1*nz1*nelt
 
-      zmin=glmin(zm1,ntot)
-      zmax=glmax(zm1,ntot)
-
       call rzero(fmode,ntot)
-         do i0=1,n_legendre
-           do j0=1,m_fourier
-             call nek_mode(fmode,i0,j0)
-             sint1=0.0 
-             do e=1,nelt
-               do f=1,6
-                 sint=0.0
-                 if (cbc(f,e,1).eq.'W  ') then
-                   call col3(cache,fmode,t,ntot)
-                   call surface_int(sint,sarea,cache,e,f)
-                   sint1=sint1+sint
-                 endif
-               enddo
-             enddo
+      do i0=1,n_legendre
+        do j0=1,m_fourier
+          call nek_mode(fmode,i0,j0)
+          sarea1=0.0
+          sint1= 0.0
+          do e=1,nelt
+            do f=1,6
+              sint=0.0
+              sarea=0.0
+              if (cbc(f,e,1).eq.'W  ') then
+                call col3(cache,fmode,t,ntot)
+                call surface_int(sint,sarea,cache,e,f)
+                sint1=sint1+sint
+                sarea1=sarea1+sarea
+              endif
+            enddo
+          enddo
           call  gop(sint1,wtmp,'+  ',1)
-!                                   
-          coeff_tij(i0,j0)=sint1*2.0/(0.5*(zmax-zmin))
-!         Note that R=0.5 here!!!!
-        enddo 
+          call  gop(sarea1,wtmp,'+  ',1)
+!
+          coeff_tij(i0,j0)=sint1*4.0*pi/sarea1
+!
+        enddo
       enddo
 
-c     For Testing 
+c     For Testing
       call nek_testp()
       return
       end
@@ -535,41 +540,47 @@ c-----------------------
       return
       end
 c-----------------------------------------------------------------------
+! calculates Legendre polynomials using a recurrence relationship. If
+! n > the maximum Legendre order, the function returns 0.0.
       function pl_leg(x,n)
-!======================================
-!    calculates Legendre polynomials Pn(x)
-!    using the recurrence relation
-!    if n > 100 the function retuns 0.0
-!======================================
+      parameter (nl_max=100)
       real*8 pl,pl_leg
       real*8 x
       real*8 pln(0:n)
       integer n, k
+
       pln(0) = 1.0
       pln(1) = x
+
       if (n.le.1) then
         pl = pln(n)
-          else
-            do k=1,n-1
-              pln(k+1)=
-     &  ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
-            end do
-            pl = pln(n)
+      else if (n.le.nl_max) then
+        do k=1,n-1
+          pln(k+1) = ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
+        end do
+        pl = pln(n)
+      else
+        pl = 0.0
       end if
-      pl_leg=pl*sqrt(dble(2*n+1)/2.0)
-      return 
+
+      pl_leg = pl*sqrt(dble(2*n+1)/2.0)
+
+      return
       end
-!======================================
+c-----------------------------------------------------------------------
+! calculates Fourier polynomials Fn(x)
       function fl_four(x,n)
       real*8 fl_four,pi
       real*8 x, A
       integer n, k
+
       pi=4.0*atan(1.0)
-      A=1.0/sqrt(2*pi)
-      if (n.gt.0) A=1.0/sqrt(pi)  
-      fl_four=A*cos(n*x)    
-c       fl_four=1.0/sqrt(pi)
-      return 
+      A=1.0/sqrt(pi)
+
+      if (n.eq.0) A=1.0/sqrt(2*pi)
+      fl_four=A*cos(n*x)
+
+      return
       end
 c-----------------------------------------------------------------------
       subroutine heat_balance(fflux)