我正在努力在有限元問題中實作周期性邊界條件，并且我想將邊界 A 上的節點與邊界 B 上的節點配對，給定一個將邊界 A 與邊界 B 對齊的向量trans。邊界 A 中的節點在串列中給出g1; 在 B 他們是g2。節點坐標在 中查找mesh%nodes(:,nodenum)。
我決定通過為每個節點創建一個距離矩陣來做到這一點，我意識到這不是最有效的方法，老實說，我不希望通過優化這個演算法來節省大量時間。這個問題更具學術性。
我知道 Fortran 以列優先順序存盤，另一方面，陣列將是對稱的，當陣列完成時，我想獲取它的列切片以找到最近的節點。所以問題是應該如何填充這個？
這是我天真的嘗試。
subroutine autopair_nodes_in_groups(mesh, g1, g2, pairs, trans)
type(meshdata) :: mesh
integer(kind=sp) :: i,j
integer(kind=sp),dimension(:) :: g1,g2
integer(kind=sp),dimension(:,:) :: pairs
real(kind=dp) :: trans(3) !xyz translate
real(kind=dp) :: dist_mat(size(g1),size(g2))
real(kind=dp) :: p1(3), p2(3)
dist_mat = -1.0_wp
! make a distance matrix
do j=1,size(g2)
p2 = mesh%nodes(1:3,g2(j))-trans
do i=1,j
p1 = mesh%nodes(1:3,g1(i))
dist_mat(i,j) = norm2(p1-p2) !equivalent to norm2(n1pos trans-n2pos)
if (i.ne.j) dist_mat(j,i) = dist_mat(i,j) !fill symmetry
end do
end do
! Remainder of routine to find nearest nodes
end subroutine autopair_nodes_in_groups
據我所知，這個問題在記憶體訪問方面是有效的，直到一個對稱陣列。
uj5u.com熱心網友回復：
要進行快速最近鄰搜索，您應該實作具有搜索復雜度 O(log(N)) 的樹結構，而不是查看所有 O(N^2) 的點到點距離。
無論如何，關于對稱矩陣處理，您將擁有：
! Storage size of a symmetric matrix
elemental integer function sym_size(N)
integer, intent(in) :: N
sym_size = (N*(N 1))/2
end function sym_size
! Compute pointer to an element in the symmetric matrix array
elemental integer function sym_ptr(N,i,j)
integer, intent(in) :: N,i,j
integer :: irow,jcol
! Column-major storage
irow = merge(i,j,i>=j)
jcol = merge(j,i,i>=j)
! Skip previous columns
sym_ptr = N*(jcol-1)-((jcol-1)*(jcol-2))/2
! Locate in current column
sym_ptr = sym_ptr (irow-jcol 1)
end function sym_ptr
然后做你的作業：
N = size(g2)
allocate(sym_dist_mat(sym_size(N)))
do j=1,size(g2)
p2 = mesh%nodes(1:3,g2(j))-trans
do i=j,size(g2)
p1 = mesh%nodes(1:3,g1(i))
sym_dist_mat(sym_ptr(N,i,j)) = norm2(p1-p2)
end do
end do
minloc 函式應該看起來像這樣（未經測驗）：
! minloc-style
pure function symmetric_minloc(N,matrix,dim) result(loc_min)
integer, intent(in) :: N
real(8), intent(in) :: matrix(N*(N 1)/2)
integer, intent(in) :: dim
real(8) :: dim_min(N),min_column
integer :: loc_min(N)
select case (dim)
! Row/column does not matter: it's symmetric!
case (1,2)
dim_min = huge(dim_row)
loc_min = -1
ptr = 0
do j=1,N
! Diagonal element m(j,j)
ptr=ptr 1
min_column = matrix(ptr)
if (min_column<=dim_min(j)) then
loc_min(j) = j
dim_min(j) = min_column
end if
! Lower-diagonal part at column j
do i=j 1,N
ptr=ptr 1
! Applies to both this column,
if (matrix(ptr)<=dim_min(j)) then
loc_min(j) = i
dim_min(j) = matrix(ptr)
end if
! And the i-th column
if (matrix(ptr)<=dim_min(i)) then
loc_min(i) = j
dim_min(i) = matrix(ptr)
end if
end do
end do
case default
! Invalid dimension
loc_min = -1
end select
end function symmetric_minloc
轉載請註明出處，本文鏈接：https://www.uj5u.com/caozuo/467612.html
下一篇：根據組中其他值的平均值創建新列