update deps

Project.toml

@@ -11,6 +11,7 @@ GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
@@ -30,11 +31,10 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 FLoops = "0.2 - 0.3"
-FastGaussQuadrature = "0.5 - 0.6"
+FastGaussQuadrature = "0.5 - 1.6"
 FoldsThreads = "0.1 - 0.2"
 GSL = "1.0 - 1.2"
 IncompleteLU = "0.2 - 0.3"
-IterativeSolvers = "0.9"
 JLD2 = "0.4 - 0.5"
 LegendrePolynomials = "0.4 - 0.5"
 LinearMaps = "3.1 - 3.20"
@@ -43,7 +43,7 @@ OffsetArrays = "1.1 - 1.30"
 Primes = "0.5 - 0.6"
 ProgressMeter = "1.5 - 1.10"
 SpecialFunctions = "2.0 - 2.5"
-StaticArrays = "1.5 - 1.6"
+StaticArrays = "1.5 - 2.6"
 ThreadsX = "0.1 - 0.2"
 UnicodePlots = "3"
 julia = "1.8 - 1.10"

src/FastAlgorithm.jl

@@ -45,7 +45,7 @@ function getImpedanceOpt(geosInfo, bfsInfo)
     # 计算近场矩阵CSC
     Znear    =   calZnearCSC(leafLevel, geosInfo, bfsInfo)
     # 构建矩阵向量乘积算子
-    Zopt    =   MLMFAIterator(Znear, octree, geosInfo, bfsInfo)
+    Zopt    =   MLFMAIterator(Znear, octree, geosInfo, bfsInfo)
     # 返回
     Zopt
 end

src/MLFMA/AdjointIterateOnOctree.jl

@@ -358,7 +358,7 @@ end
 """
 计算远区矩阵的伴随矩阵向量乘积
 """
-function caladjZfarI!(adjZopt::Adjoint{ZT, MLMFAIterator{ZT, MT}}, IVec::AbstractArray{T}; setzero = true) where {ZT, T<:Number, MT<:Vector}
+function caladjZfarI!(adjZopt::Adjoint{ZT, MLFMAIterator{ZT, MT}}, IVec::AbstractArray{T}; setzero = true) where {ZT, T<:Number, MT<:Vector}
 
     Zopt = adjZopt.parent
     # 计算前置零

src/MLFMA/IterateOnOctree.jl

@@ -531,7 +531,7 @@ end
 """
 计算远区矩阵向量乘积
 """
-function calZfarI!(Zopt::MLMFAIterator{ZT, MT}, IVec::AbstractArray{T}; setzero = true) where {ZT, T<:Number, MT<:Vector}
+function calZfarI!(Zopt::MLFMAIterator{ZT, MT}, IVec::AbstractArray{T}; setzero = true) where {ZT, T<:Number, MT<:Vector}
 
     # 计算前置零
     setzero && fill!(Zopt.ZI, zero(T))

src/MLFMA/MLFMA.jl

@@ -4,7 +4,7 @@ include("OctreeInfo.jl")
 # 盒子内几何体 id 的计算函数
 include("SetGeoIDsInLevelCubes.jl")
 # 基于八叉树的迭代器
-include("MLMFAIterators.jl")
+include("MLFMAIterators.jl")
 
 """
 计算基函数中心的数组，用于方便混合基函数使用时的情况

src/MLFMA/MLFMAIterators.jl

@@ -0,0 +1,159 @@
+## 本文件的函数用于构建迭代算子
+include("AggregateOnLevel.jl")
+
+"""
+保存 MLFMA 相关信息的结构体
+"""
+mutable struct MLFMAIterator{T<:Number, VT<:AbstractVector}
+    octree          ::OctreeInfo
+    Znear           ::AbstractMatrix{T}
+    vsCellsInfo     ::CellT where {CellT<:AbstractVector}
+    bfsInfo         ::BFT where {BFT<:AbstractVector}
+    aggSBF          ::AbstractArray{T}
+    disaggSBF       ::AbstractArray{T}
+    ZI              ::VT
+end
+
+# 各种可作为线性算子的集合
+LinearMapType{T}   =   Union{AbstractMatrix{T}, MLFMAIterator{T, VT}, LinearMap{T}} where {T<:Number, VT}
+
+Base.eltype(::MLFMAIterator{T, VT}) where {T, VT} = T
+Base.size(opt::MLFMAIterator) = size(opt.Znear)
+Base.size(opt::MLFMAIterator, ind...) = size(opt.Znear, ind...)
+
+function Base.getproperty(obj::MLFMAIterator, sym::Symbol)
+    if sym === :leafLevel
+        return obj.octree.levels[obj.nLevels]
+    elseif sym === :nLevels
+        return obj.octree.nLevels
+    elseif sym === :levels
+        return obj.octree.levels
+    else # fallback to getfield
+        return getfield(obj, sym)
+    end
+end
+
+# 多线程迭代函数
+include("IterateOnOctree.jl")
+
+
+"""
+实现矩阵向量乘积，并封装为线性算子
+"""
+function MLFMAIterator(Znear, octree::OctreeInfo{FT, LT}, 
+    vsCellsInfo::Vector, bfsInfo::Vector) where {FT<:Real, LT<:LevelInfo}
+
+    # 叶层ID
+    nLevels     =   octree.nLevels
+    # 层信息
+    levels      =   octree.levels
+    # 叶层
+    leafLevel   =   octree.levels[nLevels]
+
+    # 基函数数量
+    nbf         =   getNUnknown(bfsInfo)
+
+    # 预先计算叶层聚合项，内存占用为该层采样点数 nPoles × Nbf， 因此仅在内存充足时使用
+    @clock "叶层聚合项计算" begin
+        aggSBF, disaggSBF   =   getAggSBFOnLevel(leafLevel, vsCellsInfo, bfsInfo)
+    end
+
+    # 给各层的聚合项、解聚项预分配内存
+    memoryAllocationOnLevels!(nLevels, levels)
+
+    # 给矩阵向量乘积预分配内存
+    ZI  =   zeros(Complex{FT}, nbf)
+
+    Zopt =  MLFMAIterator{Complex{FT}, Vector}(octree, Znear, vsCellsInfo, bfsInfo, aggSBF, disaggSBF, ZI)
+
+    record_memorys(Zopt)
+
+    return Zopt
+
+end
+
+"""
+    LinearAlgebra.mul!(y, Zopt::MLFMAIterator, x)
+
+    重载以实现矩阵向量乘积计算
+TBW
+"""
+function LinearAlgebra.mul!(y, Zopt::MLFMAIterator, x)
+    # near
+    fill!(y, 0)
+    mul!(y, Zopt.Znear, x)
+    # far
+    calZfarI!(Zopt, x)
+    # cumsum results
+    axpy!(1, Zopt.ZI, y)
+
+    return y
+
+end
+
+"""
+    LinearAlgebra.mul!(y, Zopt::MLFMAIterator, x)
+
+    重载以实现矩阵向量乘积计算
+TBW
+"""
+function LinearAlgebra.mul!(y::AbstractVector, Zopt::MLFMAIterator, x::AbstractVector, α::Number, β::Number)
+    # near
+    mul!(y, Zopt.Znear, x, α, β)
+    # far
+    calZfarI!(Zopt, x)
+    # cumsum results
+    axpy!(α, Zopt.ZI, y)
+
+    return y
+
+end
+
+Base.:*(opt::MLFMAIterator, x) = mul!(deepcopy(opt.ZI), opt, x)
+
+# 实现矩阵的伴随算子矩阵向量相乘 
+Base.adjoint(opt::MLFMAIterator) = Adjoint(opt)
+Base.size(adjopt::Adjoint{T, MLFMAIterator{T, V}}) where {T, V} = size(adjoint(adjopt.parent.Znear))
+Base.size(adjopt::Adjoint{T, MLFMAIterator{T, V}}, ind...) where {T, V}  = size(adjoint(adjopt.parent.Znear), ind...)
+Base.:*(adjopt::Adjoint{T, MLFMAIterator{T, V}}, x::AbstractVector{S}) where {T, V, S} = mul!(deepcopy(adjopt.parent.ZI), adjopt, x)
+
+include("AdjointIterateOnOctree.jl")
+
+"""
+    LinearAlgebra.mul!(y, Zopt::MLFMAIterator, x)
+
+    重载以实现矩阵向量乘积计算
+TBW
+"""
+function LinearAlgebra.mul!(y, adjZopt::Adjoint{T, MLFMAIterator{T, V}}, x) where {T, V}
+    Zopt = adjZopt.parent
+    # near
+    fill!(y, 0)
+    mul!(y, adjoint(Zopt.Znear), x)
+    # far
+    caladjZfarI!(adjZopt, x)
+    # cumsum results
+    axpy!(1, Zopt.ZI, y)
+
+    return y
+
+end
+
+"""
+    LinearAlgebra.mul!(y, Zopt::MLFMAIterator, x)
+
+    重载以实现矩阵向量乘积计算
+TBW
+"""
+function LinearAlgebra.mul!(y::AbstractVector, adjZopt::Adjoint{T, MLFMAIterator{T, V}}, x::AbstractVector, α::Number, β::Number) where {T,V}
+    Zopt = adjZopt.parent
+    # near adjoint
+    mul!(y, adjoint(Zopt.Znear), x, α, β)
+    # far
+    caladjZfarI!(adjZopt, x)
+    # cumsum results
+    axpy!(α, Zopt.ZI, y)
+
+    return y
+
+end

src/MLFMA/OctreeInfo.jl

@@ -211,6 +211,7 @@ function reOrderBasisFunctionAndGeoInfo!(bfReOrderID::Vector{IT}, geosInfo::Vect
     # 再更新几何单元内的基函数 id
     @threads for geo in geosInfo
         for ii in eachindex(geo.inBfsID)
+            0 == geo.inBfsID[ii] && continue
             geo.inBfsID[ii]  = oldNewIDpair[geo.inBfsID[ii], 2]
         end
     end # cube

src/MLFMA/Znear/ZnearCFIE.jl

@@ -16,7 +16,7 @@ function calZnearCSCCFIE!(level, trianglesInfo::Vector{TriangleInfo{IT, FT}},
     # 叶层盒子数量
     nCubes      =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (RWG)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (RWG)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)

src/MLFMA/Znear/ZnearChunks.jl

@@ -3,14 +3,14 @@ abstract type ZNEARCHUNK{T} <: AbstractMatrix{T} end
 """
 创建近场矩阵结构体，所包含的数据为所有盒子内的近场矩阵元，多线程版本
 ```
-m::Int，行数
-n::Int，列数
-nChunks::Int，矩阵块儿数
-chunks::Vector{ZnearChunksStruct{T}}，矩阵
-lmul::Vector{T}，用于左乘其它矩阵、向量的临时数组，大小与列数相同
-lmuld::Vector{T}，用于左乘其它矩阵、向量的临时分布式数组，大小与列数相同，默认不分配
-rmul::Vector{T}，用于右乘其它矩阵、向量的临时数组，大小与行数相同
-lmuld::Vector{T}，用于左乘其它矩阵、向量的临时分布式数组，大小与列数相同，默认不分配
+m::Int                  行数
+n::Int                  列数
+nChunks::Int            矩阵块儿数
+chunks::Vector{ZnearChunksStruct{T}}    矩阵
+lmul::Vector{T}         用于左乘其它矩阵、向量的临时数组，大小与列数相同
+lmuld::Vector{T}        用于左乘其它矩阵、向量的临时分布式数组，大小与列数相同，默认不分配
+rmul::Vector{T}         用于右乘其它矩阵、向量的临时数组，大小与行数相同
+lmuld::Vector{T}        用于左乘其它矩阵、向量的临时分布式数组，大小与列数相同，默认不分配
 ```
 """
 mutable struct ZnearChunksStruct{T<:Number} <: ZNEARCHUNK{T}
@@ -25,7 +25,9 @@ mutable struct ZnearChunksStruct{T<:Number} <: ZNEARCHUNK{T}
 end
 
 """
-ZnearChunksStruct 类的初始化函数，将 lumld 和 rmuld 初始化为
+    ZnearChunksStruct{T}(chunks; m, n) where {T<:Number}
+
+ZnearChunksStruct 类的初始化函数。
 """
 function ZnearChunksStruct{T}(chunks; m, n) where {T<:Number}
 
@@ -71,6 +73,8 @@ function getindex(Z::T, i1::Int, i2::Int) where{T<:ZNEARCHUNK}
 end
 
 """
+    setindex!(Z::T, x, i1::Int, i2::Int) where{T<:ZNEARCHUNK}
+
 重载 setindex! 函数
 """
 function setindex!(Z::T, x, i1::Int, i2::Int) where{T<:ZNEARCHUNK}
@@ -87,7 +91,9 @@ end
 
 
 """
-初始化 阻抗矩阵 左乘 向量 乘积的 分布式数组
+    initialZchunksMulV!(Z::T) where{T<:ZnearChunksStruct}
+
+初始化 阻抗矩阵 左乘 向量 乘积的 分布式数组。
 """
 function initialZchunksMulV!(Z::T) where{T<:ZnearChunksStruct}
 
@@ -97,6 +103,8 @@ function initialZchunksMulV!(Z::T) where{T<:ZnearChunksStruct}
 end
 
 """
+    initialVMulZchunks!(Z::T) where{T<:ZnearChunksStruct}
+
 初始化 阻抗矩阵 右乘 向量 乘积的 分布式数组
 """
 function initialVMulZchunks!(Z::T) where{T<:ZnearChunksStruct}
@@ -107,6 +115,8 @@ end
 
 
 """
+    Base.:*(Z::T, x::AbstractVector) where{T<:ZnearChunksStruct}
+
 实现左乘其它向量
 """
 function Base.:*(Z::T, x::AbstractVector) where{T<:ZnearChunksStruct}
@@ -139,7 +149,7 @@ end
 """
     get_chunks_minmax_col(matchunks)
 
-TBW
+
 """
 function get_chunks_minmax_col(matchunks)::UnitRange{Int64}
 
@@ -151,6 +161,8 @@ function get_chunks_minmax_col(matchunks)::UnitRange{Int64}
 end
 
 """
+    Base.:*(Z::ZNEARCHUNK{T}, mat::AbstractMatrix) where{T<:Number}
+
 实现左乘其它矩阵
 """
 function Base.:*(Z::ZNEARCHUNK{T}, mat::AbstractMatrix) where{T<:Number}

src/MLFMA/Znear/ZnearEFIE.jl

@@ -14,7 +14,7 @@ function calZnearCSCEFIE!(level, trianglesInfo::Vector{TriangleInfo{IT, FT}},
     # 叶层盒子数量
     nCubes      =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (RWG)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (RWG)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)
@@ -188,7 +188,7 @@ function calZnearCSCEFIE!(level, tetrasInfo::AbstractVector{TetrahedraInfo{IT, F
     # 叶层盒子数量
     nCubes  =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (SWG)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (SWG)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)
@@ -403,7 +403,7 @@ function calZnearCSCEFIE!(level, tetrasInfo::AbstractVector{TetrahedraInfo{IT, F
     # 叶层盒子数量
     nCubes      =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (PWC)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (PWC)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)
@@ -548,7 +548,7 @@ function calZnearCSCEFIEnew!(level, hexasInfo::AbstractVector{HexahedraInfo{IT,
     # 叶层盒子数量
     nCubes      =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (RBF)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (RBF)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)
@@ -710,7 +710,7 @@ function calZnearCSCEFIE!(level, hexasInfo::AbstractVector{HexahedraInfo{IT, FT,
     # 线程锁防止对同一数据写入出错
     lockZ   =   SpinLock()
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (RBF)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (RBF)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)
@@ -872,7 +872,7 @@ function calZnearCSCEFIE!(level, hexasInfo::AbstractVector{HexahedraInfo{IT, FT,
     # 叶层盒子数量
     nCubes      =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (PWC)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (PWC)...", dt = 1)
     # 对盒子循环计算
     @threads for iCube in 1:nCubes
         next!(pmeter)
@@ -1025,7 +1025,7 @@ function calZnearCSCEFIE!(  level, geos1Info::AbstractVector{VT1}, geos2Info::Ab
     # 叶层盒子数量
     nCubes  =   cubesIndices.stop
     # Progress Meter
-    pmeter  =   Progress(nCubes; desc = "Calculating Znear (PWC + PWC)...", dt = 1)
+    pmeter  =   Progress(nCubes; desc = "Znear (PWC + PWC)...", dt = 1)
     # 对盒子循环计算@threads
     @floop WorkStealingEx(; basesize = 1 ) for iCube in 1:nCubes
         next!(pmeter)