[Nek5000-users] How to identify boundary faces ?

[nek5000-users at lists.mcs.anl.gov]

Hi Praveen,

Below is one approach... not tested, but it should work.

In general, I would use facint_v() from navier5.f to get the
average of (say) xm1 on a particular face.

I didn't do that here because we also needed to verify that the
face was internal.   Note that 'E  ' is not a guaranteed way to
check because there are some cases where CBC is not defined as 'E  ',
even when it is in fact an E-E boundary condition.   The reason for
this is that the E-E information is no longer used in the nek5000
simulation code.  This information comes instead from the .map file.

Paul



c-----------------------------------------------------------------------
       subroutine usrdat2

       include 'SIZE'
       include 'TOTAL'

       character*3 cbv
       integer e,f,fmid(6)
       real xmin, xmax, ymax

       common /scrns/ one(lx1,ly1,lz1,lelt)


       n = nx1*ny1*nz1*nelt
       call rone (one,n)

       ifield = 1
       call dssum(one,nx1,ny1,nz1)  ! Identify shared faces

       im = (1+nx1)/2               ! Face midpoints
       jm = (1+ny1)/2
       km = (1+nz1)/2

       fmid(4) =   1 + nx1*( jm-1) + nx1*ny1*( km-1)   !  x = -1
       fmid(2) = nx1 + nx1*( jm-1) + nx1*ny1*( km-1)   !  x =  1
       fmid(1) =  im + nx1*(  1-1) + nx1*ny1*( km-1)   !  y = -1
       fmid(3) =  im + nx1*(ny1-1) + nx1*ny1*( km-1)   !  y =  1
       fmid(5) =  im + nx1*( jm-1) + nx1*ny1*(  1-1)   !  z = -1
       fmid(6) =  im + nx1*( jm-1) + nx1*ny1*(nz1-1)   !  z =  1

       xmin = -5.0
       xmax = +10.0
       ymax =  3.0
       tol  =  1.e-5

       do e = 1,nelt      ! set boundary conditions
       do f = 1,2*ndim
          count = one(fmid(f),1,1,e)
          xf    = xm1(fmid(f),1,1,e)
          yf    = ym1(fmid(f),1,1,e)
          zf    = zm1(fmid(f),1,1,e)

          if (count.lt.1.1) then       ! This is boundary edge
             if(abs(xf-xmin).lt.tol)then
                 cbc(f,e,1) = 'v  '   ! inlet
             else if(abs(xf-xmax).lt.tol)then
                 cbc(f,e,1) = 'O  '   ! outlet
             else if(abs(yf-ymax).lt.tol)then
                 cbc(f,e,1) = 'SYM'   ! top wall
             else
                 cbc(f,e,1) = 'W  '   ! bottom wall
             endif
          endif
       enddo
       enddo

       return
       end
c-----------------------------------------------------------------------





On Wed, 30 Jan 2013, nek5000-users at lists.mcs.anl.gov wrote:

> Hello
> I want to identify boundary faces to set their boundary condition type.
> Following some of the examples, I created the following function. To find
> the location of the face, I need to know the coordinates (xf,yf) of the
> face center using which I can figure out where the face lies. Is this
> possible to find ?
>
> Thanks
> praveen
>
> C=======================================================================
>      subroutine usrdat2
>
>      include 'SIZE'
>      include 'TOTAL'
>
>      parameter(XTOL=1e-10)
>      character*3 cbv
>      integer e, f
>      real xmin, xmax, ymax
>
>      xmin = -5.0
>      xmax = +10.0
>      ymax =  3.0
>
> c set boundary conditions
>      do e = 1,nelt
>         do f = 1,2*ndim
>            cbv = cbc(f,e,1)
>            if (cbv.ne.'E  ') then       ! This is boundary edge
> c TODO: Get mid-point coordinate of face f = (xf,yf)
>               if(abs(xf-xmin).lt.XTOL)then
>                   cbc(f,e,1) = 'v  '   ! inlet
>               else if(abs(xf-xmax).lt.XTOL)then
>                   cbc(f,e,1) = 'O  '   ! outlet
>               else if(abs(yf-ymax).lt.XTOL)then
>                   cbc(f,e,1) = 'SYM'   ! top wall
>               else
>                   cbc(f,e,1) = 'W  '   ! bottom wall
>               endif
>            endif
>         enddo
>      enddo
>
>      return
>      end
>