[Nek5000-users] How to identify boundary faces ?

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

On Wed, Jan 30, 2013 at 10:13 AM, <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
>

Hello

By looking at this

http://nek5000.mcs.anl.gov/index.php/Element_Definition

I made the following code. This is a 2d problem and in blue part below I am
getting face midpoint coordinates. I am not sure I have got the numbering
of the faces correct.

Thanks
praveen

C=======================================================================
      subroutine usrdat2

      include 'SIZE'
      include 'TOTAL'

      parameter(XTOL=1e-10)
      character*3 cbv
      integer e, f
      real xmin, xmax, ymax, xf(4), yf(4)

      xmin = -5.0
      xmax = +10.0
      ymax =  3.0

c set boundary conditions
      do e = 1,nelv
         xf(3) = 0.5*(xm1(1,  1,  1,e) + xm1(nx1,1,  1,e))
         xf(2) = 0.5*(xm1(nx1,1,  1,e) + xm1(nx1,ny1,1,e))
         xf(4) = 0.5*(xm1(1,  ny1,1,e) + xm1(nx1,ny1,1,e))
         xf(1) = 0.5*(xm1(1,  ny1,1,e) + xm1(1,  1,  1,e))
         yf(3) = 0.5*(ym1(1,  1,  1,e) + ym1(nx1,1,  1,e))
         yf(2) = 0.5*(ym1(nx1,1,  1,e) + ym1(nx1,ny1,1,e))
         yf(4) = 0.5*(ym1(1,  ny1,1,e) + ym1(nx1,ny1,1,e))
         yf(1) = 0.5*(ym1(1,  ny1,1,e) + ym1(1,  1,  1,e))
         do f = 1,2*ndim
            cbv = cbc(f,e,1)
            if (cbv.ne.'E  ') then       ! This is boundary edge
               if(abs(xf(f)-xmin).lt.XTOL)then
                   cbc(f,e,1) = 'v  '   ! inlet
               else if(abs(xf(f)-xmax).lt.XTOL)then
                   cbc(f,e,1) = 'O  '   ! outlet
               else if(abs(yf(f)-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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20130130/b5973751/attachment.html>