Developing.ru

Код: Выделить всё

      subroutine velocity(num_step, third_name, v11, v22, v33)
      character(16) third_name
      real*8 vx, vy, vz, p
      real*8 v1(512), v2(512), v3(512)
      real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8)
      integer size_file, num_step(512), m, l, r, s, q
      open(unit = 2, file = third_name, form='unformatted',
     &  status='old')
      read(2) size_file
c Я считываю данные из файла и записываю их так, чтобы применить метод Лагранжа. 
      do l = 1, size_file
          read(2) vx, vy, vz, p
          do m = 1, 512
            if(l.eq.num_step(m)) then
               v1(m) = vx
               v2(m) = vy
               v3(m) = vz
            endif
          enddo
      enddo
 close(2)
c Отсортированный массив трёх компонет скорости
      do m = 1, 1
        do r = 1, 8
          do s = 1, 8
            do q = 1, 8
              v11(m,r,s,q) = v1(q+8*(s-1)+64*(r-1)+512*(m-1))
              v22(m,r,s,q) = v2(q+8*(s-1)+64*(r-1)+512*(m-1))
              v33(m,r,s,q) = v3(q+8*(s-1)+64*(r-1)+512*(m-1))
c              write(3,*) v11(m,r,s,q), v22(m,r,s,q), v33(m,r,s,q)
            enddo
          enddo
        enddo
      enddo
      end
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
 Вычисление среднего значения скорости для каждой точки 
      subroutine av_velocity(xx, yy, zz, pox, poy, poz, v11, v22, v33)
      real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8), lag(1,8,8,8)
      integer i, j, k, n, xstep(1), ystep(1), zstep(1)
      real*8 vx(1), vy(1), vz(1), xx(1), yy(1), zz(1)
      real*8 pox(1,8,8,8), poy(1,8,8,8), poz(1,8,8,8)
      call place(xx, yy, zz, xstep, ystep, zstep)
      vx = 0
      vy = 0
      vz = 0
      lag = 1
      do n = 1, 1
        do i = 1, 8
          do j = 1, 8
            do k = 1, 8
c              write(*,*) v11(n,i,j,k), v22(n,i,j,k), v33(n,i,j,k)
c Применние интерполяционных полиномов Лагранжа
              call lagrange(xx(n),yy(n),zz(n),pox(n,:,:, :) ,poy(n,:,:, :) ,
     & poz(n,:,:, :) , lag(n,:,:, :) , i, j, k)
              vx(n) = vx(n) + v11(n,i,j,k) * lag(n,i,j,k)
              vy(n) = vy(n) + v22(n,i,j,k) * lag(n,i,j,k)
              vz(n) = vz(n) + v33(n,i,j,k) * lag(n,i,j,k)
            enddo
          enddo
        enddo
      enddo
      call place(xx, yy, zz, xstep, ystep, zstep)
      do n = 1, 1
c          write(3,*) xx(n), yy(n), zz(n), vx(n), vy(n), vz(n)
      enddo
      end
c----------------------------------------------------------------------------
c Нахождение зависимости эволюции средних скоростей
c частиц в зависимости отвремени 
      subroutine v_time
 
      integer i, num_step(512)
      character(6) first_name
      character(16) third_name
      real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8), poz(1,8,8,8)
      real*8 xx(1), yy(1), zz(1),  pox(1,8,8,8), poy(1,8,8,8)
c      call polpoint(lag, num_step)
 
      third_name = 'Velo_total_00001'
      first_name = '5_file'
      open(unit = 3, file = first_name)
      call polpoint(pox, poy, poz, num_step)
      do i = 1, 1
         write(third_name(12:16),'(I5.5)') i
         call velocity(num_step, third_name, v11, v22, v33)
         call av_velocity(xx, yy, zz, pox, poy, poz, v11, v22, v33)
         write(3,*) " "
      enddo
      close(3)
  end

Код: Выделить всё

program time
 
      call v_time
 
      end
c------------------------------------------------------------------------
c Определение  номера элемента, в который попадёт точка   
      subroutine place(xx,yy,zz,xstep,ystep,zstep)
      real*8 zz(1), xx(1), yy(1), step, h, w, s
      integer m, n, k, l, xstep(1), ystep(1), zstep(1)
c Задание точек, в которых нужно найти скорость жидкости
      step = 0.1
c Размеры исходной области
      h = 4 * acos(-1.0)/3
      w = 1
      s = 4 * acos(-1.0)
      do m = 1, 1
        xx(1) = 0.0
        yy(1) = 0.0
        zz(1) = 0.0
c        xx(m + 1) = xx(m) + step
c        yy(m + 1) = yy(m) + step
c        zz(m + 1) = zz(m) + step
      enddo
c Определение номера элемента по оси z, в который попадёт точка
      do m = 1, 1
        do n = 0, 24
           if(zz(m).ge.(n*h/25) .and. (zz(m).le.((n + 1)*h/25)))
     & then
             zstep(m) = n + 1
           endif
        enddo
      enddo
c Определение номера элемента по оси x, в который попадёт точка
      do m = 1, 1
        do l = 0, 29
           if(xx(m).ge.(l*s/30) .and. (xx(m).le.((l + 1)*s/30)))
     & then
             xstep(m) = l + 1
           endif
        enddo
      enddo
c Определение номера элемента по оси y, в который попадёт точка
       do m = 1, 1
        do k = 0, 9
           if(yy(m).ge.(k*w/10) .and. (yy(m).le.((k + 1)*w/10)))
     & then
             ystep(m) = k + 1
           endif
       enddo
       enddo
      end
c--------------------------------------------------------------------------
c Вычисление базисных интерполяционных полиномов Лагранжа
      subroutine lagrange(xx, yy, zz, x, y, z, l, i, j, k)
      integer i, j, k, s, q, r
      real*8  zz, xx, yy, x(8,8,8), y(8,8,8), z(8,8,8), l(8,8,8)
      l = 1
      do r = 1, 8
        do s = 1, 8
          do q = 1, 8
            if(i.ne.r .and. j.ne.s .and. k.ne.q) then
              l(i, j, k) = l(i, j, k) *
     & * (xx - x(r,s,q))/(x(i,j,k) - x(r,s,q))
     & * (yy - y(r,s,q))/(y(i,j,k) - y(r,s,q))
     & * (zz - z(r,s,q))/(z(i,j,k) - z(r,s,q))
            else
               l(i, j, k) = l(i, j, k) * 1
            endif
 enddo
        enddo
      enddo
c      write(*,*) l(i,j,k)
      end
c--------------------------------------------------------------------------
c Сортировка точек в элементе с выделенной точкой
c для применения базисных Полиномов
      subroutine polpoint(pox, poy, poz, num_step)
      integer xstep(1), ystep(1), zstep(1)
      integer l, n, num_step(512)
      integer r, q, s
      integer size_file, n1, el, gel, ex, ey, ez, i, j, k
      real*8 x, y, z, zz(1), xx(1), yy(1), t(512)
      real*8 u(512), v(512), pox(1,8,8,8)
      real*8 poy(1,8,8,8), poz(1,8,8,8)
      character(10) second_name
c Файл с расположением точек 
      second_name = 'Geom_total'
 
      call place(xx,yy,zz,xstep,ystep,zstep)
 
      open(unit = 1, file = second_name, form = 'unformatted',
     & status='old')
 
      read(1) size_file
c Я считываю данные из файла и записываю их так, чтобы применить метод Лагранжа. 
      do l = 1, size_file
        read(1) n1, el, gel, ex, ey, ez, i, j, k, x,
     &  y,  z
        do n = 1, 1
          if(ex.eq.xstep(n) .and. (ey-10).eq.ystep(n)) then
            if(ez.eq.zstep(n)) then
              do r = 1, 8
                if(i.eq.r) then
                  do s = 1, 8
                    if(j.eq.s) then
                      do q = 1, 8
                        if(k.eq.q) then
                          t(q+8*(s-1)+64*(r-1)+512*(n-1)) = x
                          u(q+8*(s-1)+64*(r-1)+512*(n-1)) = y
                          v(q+8*(s-1)+64*(r-1)+512*(n-1)) = z
                          num_step(q+8*(s-1)+64*(r-1)+512*(n-1)) = l
                        endif
                      enddo
                    endif
                  enddo
                endif
              enddo
            endif
          endif
        enddo
      enddo
      close(1)
c Отсортированные точки в спекиральном элементе
      do n = 1, 1
        do r = 1, 8
          do s = 1, 8
            do q = 1, 8
              pox(n,r,s,q) = t(q+8*(s-1)+64*(r-1)+512*(n-1))
              poy(n,r,s,q) = u(q+8*(s-1)+64*(r-1)+512*(n-1))
              poz(n,r,s,q) = v(q+8*(s-1)+64*(r-1)+512*(n-1))
c              write(3,*) pox(m,r,s,q), poy(m,r,s,q), poz(m,r,s,q)
            enddo
          enddo
        enddo
      enddo
      end