[Nek5000-users] Mesh-morphing for bend pipe

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

Dear Neks,


I have been looking at the mesh-morphing for a bend pipe recently. I found the thread https://lists.mcs.anl.gov/mailman/htdig/nek5000-users/2011-June/001398.html on how to bend a straight pipe. The code worked perfectly for bending a straight pipe into a 90 degree bend shape. The problem is I need to bend my pipe at certain section, say I have a pipe with total length of 10 and I want to bend it in the second half, i.e. from 5-10. However, when I tried to do it in this way, the error of 'Vanishing Jacobian near Xth node of element XX' appeared and the simulation couldn't start. I can't see why it doesn't work.


Below is the code I'm using to bend the pipe (from 5-10) into a 90 degree bend. Hope someone could help me on this. Any suggestion would be appreciated. Thank you very much in advance.


Kind regards,


Tony


c-----------------------------------------------------------------------
      subroutine usrdat2()  ! This routine to modify mesh coordinates
      include 'SIZE'
      include 'TOTAL'

      n = nx1*ny1*nz1*nelv

c    First, rotate x into axial position
      do i=1,n
         x_original = xm1(i,1,1,1)
         y_original = ym1(i,1,1,1)
         z_original = zm1(i,1,1,1)

         xm1(i,1,1,1) =  z_original
         ym1(i,1,1,1) =  y_original
         zm1(i,1,1,1) = -x_original
      enddo

c    Second, bend pipe into 90 degree bend

      xmin = glmin(xm1,n)
      xmax = glmax(xm1,n)
      ymin = glmin(ym1,n)
      ymax = glmax(ym1,n)
      if (nid==0) print *,xmin,xmax,ymin,ymax            !which are 0, 10, -1, 1
      rad0 = 0.5*(ymax-ymin) ! Radius of initial pipe
      rad1 = 1.0            ! Radius of new pipe
      radm = 2.0            ! Major radius of torus

      do i=1,n
         x     = xm1(i,1,1,1)
         y     = ym1(i,1,1,1)

       if (x .gt. 5.0e0) then
         theta = 0.5*pi*(x-5.0e0)/(xmax-5.0e0)
         rad   = radm + 2.0e0*rad1*(y-ymin)/(ymax-ymin)

         ym1(i,1,1,1) = rad*cos(theta)
         xm1(i,1,1,1) = rad*sin(theta)+5.0e0
       end if

      end do

      return
      end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20150121/e9737b3d/attachment.html>