У меня есть 3D-матрица (размерность nx, nz, ny), которая соответствует физической области.Эта матрица содержит непрерывное поле от -1 (фаза 1) до +1 (фаза 2);интерфейс между этими двумя фазами - уровень этого поля 0.
Теперь я хочу эффективно рассчитать функцию расстояния со знаком для интерфейса для каждой точки в домене.
Я попробовал двавозможности (sgn - это знак моего поля, со значениями + 1,0, -1, xyz содержит сетку в виде тройки x, y, z в каждой точке, а dist - это функция расстояния со знаком, которую я хочу вычислить).
double precision, dimension(nx,nz,ny) :: dist,sgn,eudist
integer :: i,j,k
double precision :: seed,posit,tmp(nx)
do j=1,ny
do k=1,nz
do i=1,nx
seed=sgn(i,k,j)
! look for interface
eudist=(xyz(:,:,:,1)-x(i))**2+(xyz(:,:,:,2)-z(k))**2+(xyz(:,:,:,3)-y(j))**2
! find min within mask
posit=minval(eudist,seed*sgn.le.0)
! tmp fits in cache, small speed-up
tmp(i)=-seed*dsqrt(posit)
enddo
dist(:,k,j)=tmp
enddo
enddo
Я также попробовал вторую версию, которая очень похожа на предыдущую, но она рассчитывает евклидово расстояние только в подмножестве всей матрицы.В этой второй версии есть некоторое ускорение, но оно все еще слишком медленное.Я хотел бы знать, есть ли более эффективный способ вычисления функции расстояния.
Вторая версия:
double precision, dimension(nx,nz,ny) :: dist,sgn
double precision, allocatable, dimension(:,:,:) :: eudist
integer :: i,j,k , ii,jj,kk
integer :: il,iu,jl,ju,kl,ku
double precision :: seed, deltax,deltay,deltaz,tmp(nx)
deltax=max(int(nx/4),1)
deltay=max(int(ny/4),1)
deltaz=max(int(nz/2),1)
allocate(eudist(2*deltax+1,2*deltaz+1,2*deltay+1))
do j=1,ny
do k=1,nz
do i=1,nx
! look for closest point in box 2*deltax+1,2*deltaz+1,2*deltay+1
il=max(1,i-deltax)
iu=min(nx,i+deltax)
jl=max(1,j-deltay)
ju=min(ny,j+deltay)
kl=max(1,k-deltaz)
ku=min(nz,k+deltaz)
eudist(:,1:ku-kl+1,:)=(xyz(il:iu,kl:ku,jl:ju,1)-x(i))**2 &
& +(xyz(il:iu,kl:ku,jl:ju,2)-z(k))**2 &
& +(xyz(il:iu,kl:ku,jl:ju,3)-y(j))**2
seed=sgn(i,k,j)
tmp(i)=minval(eudist(:,1:ku-kl+1,:),seed*sgn(il:iu,kl:ku,jl:ju).le.0)
tmp(i)=-seed*dsqrt(tmp(i))
enddo
dist(:,k,j)=tmp
enddo
enddo
eudist: евклидово расстояние между точкой i, k, j и любымдругая точка в блоке 2 * deltax + 1,2 * deltaz + 1,2 * deltay + 1 с центром в i, k, j.Это снижает вычислительные затраты, так как расстояние вычисляется только в подмножестве всей сетки (здесь я предполагаю, что подмножество достаточно велико, чтобы содержать точку раздела).
После предположения Владимира (x, y,z - оси, определяющие положение сетки, xyz (i, k, j) = (x (i), z (k), y (j))):
double precision, dimension(nx,nz,ny) :: dist,sgn
double precision :: x(nx), y(ny), z(nz)
double precision, allocatable, dimension(:,:,:) :: eudist
double precision, allocatable, dimension(:) :: xd,yd,zd
integer :: i,j,k , ii,jj,kk
integer :: il,iu,jl,ju,kl,ku
double precision :: seed, deltax,deltay,deltaz,tmp(nx)
deltax=max(int(nx/4),1)
deltay=max(int(ny/4),1)
deltaz=max(int(nz/2),1)
allocate(eudist(2*deltax+1,2*deltaz+1,2*deltay+1))
allocate(xd(2*deltax+1))
allocate(yd(2*deltay+1))
allocate(zd(2*deltaz+1))
do j=1,ny
do k=1,nz
do i=1,nx
! look for closest point in box 2*deltax+1,2*deltaz+1,2*deltay+1
il=max(1,i-deltax)
iu=min(nx,i+deltax)
jl=max(1,j-deltay)
ju=min(ny,j+deltay)
kl=max(1,k-deltaz)
ku=min(nz,k+deltaz)
do ii=1,iu-il+1
xd(ii)=(xyz(il+ii-1)-x(i))**2
end do
do jj=1,ju-jl+1
yd(jj)=(y(jj+jl-1)-y(j))**2
end do
do kk=1,ku-kl+1
zd(kk)=(z(kk+kl-1)-z(k))**2
end do
do jj=1,ju-jl+1
do kk=1,ku-kl+1
do ii=1,iu-il+1
eudist(ii,kk,jj)=xd(ii)+yd(jj)+zd(kk)
enddo
enddo
enddo
seed=sgn(i,k,j)
tmp(i)=minval(eudist(:,1:ku-kl+1,:),seed*sgn(il:iu,kl:ku,jl:ju).le.0)
tmp(i)=-seed*dsqrt(tmp(i))
enddo
dist(:,k,j)=tmp
enddo
enddo
РЕДАКТИРОВАТЬ: дополнительная информация опроблема под рукой.Сетка - это ортогональная сетка, сопоставленная с матрицей.Количество точек этой сетки составляет порядка 1000 в каждом направлении (всего около 1 миллиарда точек).
Моя цель - перейти от функции знака (+ 1,0, -1) кФункция расстояния со знаком во всей сетке эффективным способом.