Commits

Nathan Collier  committed d6074fe

removed all data copies

  • Participants
  • Parent commits 706eaa6

Comments (0)

Files changed (2)

File src/pyfly/uvlm.f90

   use base
   implicit none
   integer(kind=4), intent(in)  :: nr,nc,symaxis
-  real   (kind=8), intent(in)  :: grid(nr+1,nc+1,3),normals(nr,nc,3),cutoff,symcord
+  real   (kind=8), intent(in)  :: grid(3,nc+1,nr+1),normals(3,nc,nr),cutoff,symcord
   real   (kind=8), intent(out) :: A(nr*nc,nr*nc)
-  real   (kind=8)              :: mid(3),vel(3),mirror(nr+1,nc+1,3)
+  real   (kind=8)              :: mid(3),vel(3),mirror(3,nc+1,nr+1)
   integer(kind=4)              :: i,j,k,l
   if(symaxis .ge. 0) then
      mirror = grid
-     mirror(:,:,symaxis+1) = -grid(:,:,symaxis+1) + 2.0d0*symcord
+     mirror(symaxis+1,:,:) = -grid(symaxis+1,:,:) + 2.0d0*symcord
   end if
   do i=1,nr
      do j=1,nc
-        mid = 0.25d0*(grid(i,j,:)+grid(i+1,j,:)+grid(i,j+1,:)+grid(i+1,j+1,:))
+        mid = 0.25d0*(grid(:,j,i)+grid(:,j,i+1)+grid(:,j+1,i)+grid(:,j+1,i+1))
         do k=1,nr
            do l=1,nc
-              vel =       BiotSavartLaw(grid(k  ,l  ,:),grid(k  ,l+1,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid(k  ,l+1,:),grid(k+1,l+1,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid(k+1,l+1,:),grid(k+1,l  ,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid(k+1,l  ,:),grid(k  ,l  ,:),mid,cutoff)
-              A((i-1)*nc+j,(k-1)*nc+l) = DotProduct(vel,normals(i,j,:))
+              vel =     + BiotSavartLaw(grid(:,l  ,k  ),grid(:,l+1,k  ),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l+1,k  ),grid(:,l+1,k+1),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l+1,k+1),grid(:,l  ,k+1),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l  ,k+1),grid(:,l  ,k  ),mid,cutoff)
+              A((k-1)*nc+l,(i-1)*nc+j) = DotProduct(vel,normals(:,j,i))
            end do
         end do
         if(symaxis .ge. 0) then
            do k=1,nr
               do l=1,nc
-                 vel =       BiotSavartLaw(mirror(k  ,l  ,:),mirror(k  ,l+1,:),mid,cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k  ,l+1,:),mirror(k+1,l+1,:),mid,cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k+1,l+1,:),mirror(k+1,l  ,:),mid,cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k+1,l  ,:),mirror(k  ,l  ,:),mid,cutoff)
+                 vel =     + BiotSavartLaw(mirror(:,l  ,k  ),mirror(:,l+1,k  ),mid,cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l+1,k  ),mirror(:,l+1,k+1),mid,cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l+1,k+1),mirror(:,l  ,k+1),mid,cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l  ,k+1),mirror(:,l  ,k  ),mid,cutoff)
                  vel = -vel
-                 A((i-1)*nc+j,(k-1)*nc+l) = A((i-1)*nc+j,(k-1)*nc+l) + DotProduct(vel,normals(i,j,:))
+                 A((k-1)*nc+l,(i-1)*nc+j) = A((k-1)*nc+l,(i-1)*nc+j) + DotProduct(vel,normals(:,j,i))
               end do
            end do
         end if
   use base
   implicit none
   integer(kind=4), intent(in)  :: nr1,nc1,nr2,nc2,symaxis2
-  real   (kind=8), intent(in)  :: grid1(nr1+1,nc1+1,3),normals1(nr1,nc1,3),cutoff,symcord2
-  real   (kind=8), intent(in)  :: grid2(nr2+1,nc2+1,3)
-  real   (kind=8), intent(out) :: A(nr1*nc1,nr2*nc2)
-  real   (kind=8)              :: mid(3),vel(3),mirror(nr2+1,nc2+1,3)
+  real   (kind=8), intent(in)  :: grid1(3,nc1+1,nr1+1),normals1(3,nc1,nr1),cutoff,symcord2
+  real   (kind=8), intent(in)  :: grid2(3,nc2+1,nr2+1)
+  real   (kind=8), intent(out) :: A(nr2*nc2,nr1*nc1)
+  real   (kind=8)              :: mid(3),vel(3),mirror(3,nc2+1,nr2+1)
   integer(kind=4)              :: i,j,k,l
   if(symaxis2 .ge. 0) then
      mirror = grid2
-     mirror(:,:,symaxis2+1) = -grid2(:,:,symaxis2+1) + 2.0d0*symcord2
+     mirror(symaxis2+1,:,:) = -grid2(symaxis2+1,:,:) + 2.0d0*symcord2
   end if
   do i=1,nr1
      do j=1,nc1
-        mid = 0.25d0*(grid1(i,j,:)+grid1(i+1,j,:)+grid1(i,j+1,:)+grid1(i+1,j+1,:))
+        mid = 0.25d0*(grid1(:,j,i)+grid1(:,j,i+1)+grid1(:,j+1,i)+grid1(:,j+1,i+1))
         do k=1,nr2
            do l=1,nc2
-              vel =       BiotSavartLaw(grid2(k  ,l  ,:),grid2(k  ,l+1,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid2(k  ,l+1,:),grid2(k+1,l+1,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid2(k+1,l+1,:),grid2(k+1,l  ,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid2(k+1,l  ,:),grid2(k  ,l  ,:),mid,cutoff)
-              A((i-1)*nc1+j,(k-1)*nc2+l) = DotProduct(vel,normals1(i,j,:))
+              vel =     + BiotSavartLaw(grid2(:,l  ,k  ),grid2(:,l+1,k  ),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid2(:,l+1,k  ),grid2(:,l+1,k+1),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid2(:,l+1,k+1),grid2(:,l  ,k+1),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid2(:,l  ,k+1),grid2(:,l  ,k  ),mid,cutoff)
+              A((k-1)*nc2+l,(i-1)*nc1+j) = DotProduct(vel,normals1(:,j,i))
            end do
         end do
         if(symaxis2 .ge. 0) then
            do k=1,nr2
               do l=1,nc2
-                 vel =       BiotSavartLaw(mirror(k  ,l  ,:),mirror(k  ,l+1,:),mid,cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k  ,l+1,:),mirror(k+1,l+1,:),mid,cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k+1,l+1,:),mirror(k+1,l  ,:),mid,cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k+1,l  ,:),mirror(k  ,l  ,:),mid,cutoff)
+                 vel =     + BiotSavartLaw(mirror(:,l  ,k  ),mirror(:,l+1,k  ),mid,cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l+1,k  ),mirror(:,l+1,k+1),mid,cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l+1,k+1),mirror(:,l  ,k+1),mid,cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l  ,k+1),mirror(:,l  ,k  ),mid,cutoff)
                  vel = -vel
-                 A((i-1)*nc1+j,(k-1)*nc2+l) = A((i-1)*nc1+j,(k-1)*nc2+l) + DotProduct(vel,normals1(i,j,:))
+                 A((k-1)*nc2+l,(i-1)*nc1+j) = A((k-1)*nc2+l,(i-1)*nc1+j) + DotProduct(vel,normals1(:,j,i))
               end do
            end do
         end if
   end do
 end subroutine Grid2OnGrid1
 
-function VelocityAtPoints(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,cutoff) result (velocity)
+subroutine VelocityAtPoints(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,cutoff,velocity)
   use base
   implicit none
   integer(kind=4), intent(in)  :: nr,nc,symaxis,npx,npy
-  real   (kind=8), intent(in)  :: grid(nr+1,nc+1,3),circ(nr,nc),symcord,cutoff
-  real   (kind=8), intent(in)  :: points(npx,npy,3)
-  real   (kind=8)              :: velocity(npx,npy,3),vel(3),mirror(nr+1,nc+1,3)
+  real   (kind=8), intent(in)  :: grid(3,nc+1,nr+1),circ(nc,nr),symcord,cutoff
+  real   (kind=8), intent(in)  :: points(3,npy,npx)
+  real   (kind=8), intent(out) :: velocity(3,npy,npx)
+  real   (kind=8)              :: vel(3),mirror(3,nc+1,nr+1)
+  integer(kind=4)              :: i,j,k,l
+  if(symaxis .ge. 0) then
+     mirror = grid
+     mirror(symaxis+1,:,:) = -grid(symaxis+1,:,:) + 2.0d0*symcord
+  end if
+  do i=1,npx
+     do j=1,npy
+        do k=1,nr
+           do l=1,nc
+              vel =     + BiotSavartLaw(grid(:,l  ,k  ),grid(:,l+1,k  ),points(:,j,i),cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l+1,k  ),grid(:,l+1,k+1),points(:,j,i),cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l+1,k+1),grid(:,l  ,k+1),points(:,j,i),cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l  ,k+1),grid(:,l  ,k  ),points(:,j,i),cutoff)
+              velocity(:,j,i) = vel * circ(l,k)
+           end do
+        end do
+        if(symaxis .ge. 0) then
+           do k=1,nr
+              do l=1,nc
+                 vel =     + BiotSavartLaw(mirror(:,l  ,k  ),mirror(:,l+1,k  ),points(:,j,i),cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l+1,k  ),mirror(:,l+1,k+1),points(:,j,i),cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l+1,k+1),mirror(:,l  ,k+1),points(:,j,i),cutoff)
+                 vel = vel + BiotSavartLaw(mirror(:,l  ,k+1),mirror(:,l  ,k  ),points(:,j,i),cutoff)
+                 velocity(:,j,i) = -vel * circ(l,k)
+              end do
+           end do
+        end if
+     end do
+  end do
+end subroutine VelocityAtPoints
+
+subroutine NormalVelocityAtCenter(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,normal,cutoff,velocity)
+  use base
+  implicit none
+  integer(kind=4), intent(in)  :: nr,nc,symaxis,npx,npy
+  real   (kind=8), intent(in)  :: grid(3,nc+1,nr+1),circ(nc,nr),symcord,cutoff
+  real   (kind=8), intent(in)  :: points(3,npy+1,npx+1),normal(3,npy,npx)
+  real   (kind=8), intent(out) :: velocity(npy,npx)
+  real   (kind=8)              :: vel(3),mirror(3,nc+1,nr+1),mid(3)
   integer(kind=4)              :: i,j,k,l
   if(symaxis .ge. 0) then
      mirror = grid
      mirror(:,:,symaxis+1) = -grid(:,:,symaxis+1) + 2.0d0*symcord
   end if
+  velocity = 0.
   do i=1,npx
      do j=1,npy
+        mid = 0.25d0*(points(:,j,i)+points(:,j,i+1)+points(:,j+1,i)+points(:,j+1,i+1))
         do k=1,nr
            do l=1,nc
-              vel =       BiotSavartLaw(grid(k  ,l  ,:),grid(k  ,l+1,:),points(i,j,:),cutoff)
-              vel = vel + BiotSavartLaw(grid(k  ,l+1,:),grid(k+1,l+1,:),points(i,j,:),cutoff)
-              vel = vel + BiotSavartLaw(grid(k+1,l+1,:),grid(k+1,l  ,:),points(i,j,:),cutoff)
-              vel = vel + BiotSavartLaw(grid(k+1,l  ,:),grid(k  ,l  ,:),points(i,j,:),cutoff)
-              velocity(i,j,:) = vel * circ(k,l)
-           end do
-        end do
-        if(symaxis .ge. 0) then
-           do k=1,nr
-              do l=1,nc
-                 vel =       BiotSavartLaw(mirror(k  ,l  ,:),mirror(k  ,l+1,:),points(i,j,:),cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k  ,l+1,:),mirror(k+1,l+1,:),points(i,j,:),cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k+1,l+1,:),mirror(k+1,l  ,:),points(i,j,:),cutoff)
-                 vel = vel + BiotSavartLaw(mirror(k+1,l  ,:),mirror(k  ,l  ,:),points(i,j,:),cutoff)
-                 velocity(i,j,:) = -vel * circ(k,l)
-              end do
-           end do
-        end if
-     end do
-  end do
-end function VelocityAtPoints
-
-function NormalVelocityAtCenter(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,normal,cutoff) result (velocity)
-  use base
-  implicit none
-  integer(kind=4), intent(in)  :: nr,nc,symaxis,npx,npy
-  real   (kind=8), intent(in)  :: grid(nr+1,nc+1,3),circ(nr,nc),symcord,cutoff
-  real   (kind=8), intent(in)  :: points(npx+1,npy+1,3),normal(npx,npy,3)
-  real   (kind=8)              :: velocity(npx,npy),vel(3),mirror(nr+1,nc+1,3),mid(3)
-  integer(kind=4)              :: i,j,k,l
-  if(symaxis .ge. 0) then
-     mirror = grid
-     mirror(:,:,symaxis+1) = -grid(:,:,symaxis+1) + 2.0d0*symcord
-  end if
-  do i=1,npx
-     do j=1,npy
-        mid = 0.25d0*(points(i,j,:)+points(i+1,j,:)+points(i,j+1,:)+points(i+1,j+1,:))
-        do k=1,nr
-           do l=1,nc
-              vel =       BiotSavartLaw(grid(k  ,l  ,:),grid(k  ,l+1,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid(k  ,l+1,:),grid(k+1,l+1,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid(k+1,l+1,:),grid(k+1,l  ,:),mid,cutoff)
-              vel = vel + BiotSavartLaw(grid(k+1,l  ,:),grid(k  ,l  ,:),mid,cutoff)
-              velocity(i,j) = circ(k,l)*DotProduct(vel,normal(i,j,:))
+              vel =     + BiotSavartLaw(grid(:,l  ,k  ),grid(:,l+1,k  ),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l+1,k  ),grid(:,l+1,k+1),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l+1,k+1),grid(:,l  ,k+1),mid,cutoff)
+              vel = vel + BiotSavartLaw(grid(:,l  ,k+1),grid(:,l  ,k  ),mid,cutoff)
+              velocity(j,i) = velocity(j,i) + circ(k,l)*DotProduct(vel,normal(:,j,i))
            end do
         end do
         if(symaxis .ge. 0) then
                  vel = vel + BiotSavartLaw(mirror(k  ,l+1,:),mirror(k+1,l+1,:),mid,cutoff)
                  vel = vel + BiotSavartLaw(mirror(k+1,l+1,:),mirror(k+1,l  ,:),mid,cutoff)
                  vel = vel + BiotSavartLaw(mirror(k+1,l  ,:),mirror(k  ,l  ,:),mid,cutoff)
-                 velocity(i,j) = -circ(k,l)*DotProduct(vel,normal(i,j,:))
+                 velocity(j,i) = velocity(j,i) + circ(k,l)*DotProduct(vel,normal(:,j,i))
               end do
            end do
         end if
      end do
   end do
-end function NormalVelocityAtCenter
+end subroutine NormalVelocityAtCenter
 
 end module influence

File src/pyfly/uvlm.pyf

 
      subroutine GridOnSelf(nr,nc,grid,normals,symaxis,symcord,cutoff,A)
        implicit none
-       integer(kind=4), intent(hide) :: nr,nc
-       integer(kind=4), optional     :: symaxis = -1
-       real   (kind=8), optional     :: symcord = 0.
-       real   (kind=8), intent(in)   :: grid(nr+1,nc+1,3)
-       real   (kind=8), intent(in)   :: normals(nr,nc,3)
-       real   (kind=8), optional     :: cutoff = 1.0e-9
-       real   (kind=8), intent(out)  :: A(nr*nc,nr*nc)
+       integer(kind=4), intent(hide)  :: nr,nc
+       integer(kind=4), optional      :: symaxis = -1
+       real   (kind=8), optional      :: symcord = 0.
+       real   (kind=8), intent(c,in)  :: grid(nr+1,nc+1,3)
+       real   (kind=8), intent(c,in)  :: normals(nr,nc,3)
+       real   (kind=8), optional      :: cutoff = 1.0e-9
+       real   (kind=8), intent(c,out) :: A(nr*nc,nr*nc)
      end subroutine GridOnSelf
 
      subroutine Grid2OnGrid1(nr1,nc1,grid1,normals1,nr2,nc2,grid2,symaxis2,symcord2,cutoff,A)
        implicit none
-       integer(kind=4), intent(hide) :: nr1,nc1,nr2,nc2
-       integer(kind=4), optional     :: symaxis2 = -1
-       real   (kind=8), optional     :: symcord2 = 0.
-       real   (kind=8), intent(in)   :: grid1(nr1+1,nc1+1,3)
-       real   (kind=8), intent(in)   :: grid2(nr2+1,nc2+1,3)
-       real   (kind=8), intent(in)   :: normals1(nr1,nc1,3)
-       real   (kind=8), optional     :: cutoff = 1.0e-9
-       real   (kind=8), intent(out)  :: A(nr1*nc1,nr2*nc2)
+       integer(kind=4), intent(hide)  :: nr1,nc1,nr2,nc2
+       integer(kind=4), optional      :: symaxis2 = -1
+       real   (kind=8), optional      :: symcord2 = 0.
+       real   (kind=8), intent(c,in)  :: grid1(nr1+1,nc1+1,3)
+       real   (kind=8), intent(c,in)  :: grid2(nr2+1,nc2+1,3)
+       real   (kind=8), intent(c,in)  :: normals1(nr1,nc1,3)
+       real   (kind=8), optional      :: cutoff = 1.0e-9
+       real   (kind=8), intent(c,out) :: A(nr1*nc1,nr2*nc2)
      end subroutine Grid2OnGrid1
 
-     function VelocityAtPoints(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,cutoff) result (velocity)
+     subroutine VelocityAtPoints(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,cutoff,velocity)
        implicit none
        integer(kind=4), intent(hide)  :: nr,nc,npx,npy
        integer(kind=4), optional      :: symaxis = -1
        real   (kind=8), optional      :: symcord = 0.,cutoff = 1.0e-9
-       real   (kind=8), intent(in)    :: grid(nr+1,nc+1,3),circ(nr,nc)
-       real   (kind=8), intent(in)    :: points(npx,npy,3)
-       real   (kind=8), intent(out)   :: velocity(npx,npy,3)
-     end function VelocityAtPoints
+       real   (kind=8), intent(c,in)  :: grid(nr+1,nc+1,3),circ(nr,nc)
+       real   (kind=8), intent(c,in)  :: points(npx,npy,3)
+       real   (kind=8), intent(c,out) :: velocity(npx,npy,3)
+     end subroutine VelocityAtPoints
 
-     function NormalVelocityAtCenter(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,normal,cutoff) result (velocity)
+     subroutine NormalVelocityAtCenter(nr,nc,grid,circ,symaxis,symcord,npx,npy,points,normal,cutoff,velocity)
        implicit none
-       integer(kind=4), intent(hide)  :: nr,nc,npx,npy,npz
+       integer(kind=4), intent(hide)  :: nr,nc,npx,npy
        integer(kind=4), optional      :: symaxis = -1
        real   (kind=8), optional      :: symcord = 0.,cutoff = 1.0e-9
-       real   (kind=8), intent(in)    :: grid(nr+1,nc+1,3),circ(nr,nc)
-       real   (kind=8), intent(in)    :: points(npx+1,npy+1,3),normal(npx,npy,3)
-       real   (kind=8), intent(out)   :: velocity(npx,npy)
-     end function NormalVelocityAtCenter
+       real   (kind=8), intent(c,in)  :: grid(nr+1,nc+1,3),circ(nr,nc)
+       real   (kind=8), intent(c,in)  :: points(npx+1,npy+1,3),normal(npx,npy,3)
+       real   (kind=8), intent(c,out) :: velocity(npx,npy)
+     end subroutine NormalVelocityAtCenter
 
    end module influence