Feature/blockassembler

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.1"
 
 [deps]
 BEAST = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
+BlockSparseMatrices = "1ef78180-cc57-4e70-b6ec-c3e471307147"
 ClusterTrees = "5100927d-02aa-593a-b4f9-7235df19f0db"
 CompScienceMeshes = "3e66a162-7b8c-5da0-b8f8-124ecd2c3ae1"
 ExaFMMt = "add1291b-d076-47cc-8667-64576de04ee0"
@@ -18,6 +19,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 
 [compat]
 BEAST = "2.4.0"

src/blockassembler/galerkinblkassembler.jl

@@ -0,0 +1,74 @@
+using ClusterTrees
+using ThreadsX
+
+function assemble(
+    op::BEAST.AbstractOperator,
+    space::BEAST.Space,
+    tree::ClusterTrees.BlockTrees.BlockTree{T},
+    interactions::Vector{Tuple{I, I}},
+    ::Type{K};
+    quadstrat=BEAST.defaultquadstrat(op, space, space),
+    multithreading=true
+) where {I,K,T}
+
+    @views assembler = BEAST.blockassembler(
+        op, space, space, quadstrat=quadstrat
+    )
+    @views function nearassembler(Z, tdata, sdata)
+        @views store(v,m,n) = (Z[m,n] += v)
+        assembler(tdata,sdata,store)
+    end
+
+    nears = sort(interactions)
+    snears = Tuple{Int, Vector{Int}}[]
+    selfs = Tuple{Int, Int}[]
+    for near in nears
+        if near[2] > near[1]
+            if length(snears)!= 0 && snears[end][1] == near[1]
+                push!(snears[end][2], near[2])
+            else
+                push!(snears, (near[1], [near[2]]))
+            end
+        elseif near[1] == near[2]
+            push!(selfs, near)
+        end
+    end
+
+    _foreach = multithreading ? ThreadsX.foreach : Base.foreach
+    blocks = Vector{BlockSparseMatrices.DenseMatrixBlock{K, Matrix{K}, Vector{I}}}(
+        undef, length(snears)
+    )
+    diagonals = Vector{BlockSparseMatrices.DenseMatrixBlock{K, Matrix{K}, Vector{I}}}(
+        undef, length(selfs)
+    )
+
+    _foreach(enumerate(snears)) do (idx, snear) 
+        testidcs = value(tree.test_cluster, snear[1])
+        trialidcs = value(tree.trial_cluster, snear[2])
+        matrix = zeros(K, length(testidcs), length(trialidcs))
+        assembler(matrix, testidcs, trialidcs)
+
+        blocks[idx] =  BlockSparseMatrices.DenseMatrixBlock{K, Matrix{K}, Vector{Int}}(
+            matrix,
+            testidcs,
+            trialidcs
+        )
+    end
+
+    _foreach(enumerate(selfs)) do (idx, self) 
+        testidcs = value(tree.test_cluster, self[1])
+        trialidcs = value(tree.trial_cluster, self[2])
+        matrix = zeros(K, length(testidcs), length(trialidcs))
+        assembler(matrix, testidcs, trialidcs)
+
+        diagonals[idx] =  BlockSparseMatrices.DenseMatrixBlock{K, Matrix{K}, Vector{Int}}(
+            matrix,
+            testidcs,
+            trialidcs
+        )
+    end
+
+    return SymmetricBlockMatrix(diagonals, blocks, (
+        tree.test_cluster.num_elements, tree.trial_cluster.num_elements
+    ))
+end

src/blockassembler/petrovgalerkinblkassembler.jl

@@ -0,0 +1,53 @@
+using ClusterTrees
+using ThreadsX
+
+function assemble(
+    op::BEAST.AbstractOperator,
+    testspace::BEAST.Space,
+    trialspace::BEAST.Space,
+    tree::ClusterTrees.BlockTrees.BlockTree{T},
+    interactions::Vector{Tuple{I, I}},
+    ::Type{K};
+    quadstrat=BEAST.defaultquadstrat(op, testspace, trialspace),
+    multithreading=true
+) where {I,K,T}
+
+    @views nearblkassembler = BEAST.blockassembler(
+        op, testspace, trialspace, quadstrat=quadstrat
+    )
+    @views function nearassembler(Z, tdata, sdata)
+        @views store(v,m,n) = (Z[m,n] += v)
+        nearblkassembler(tdata,sdata,store)
+    end
+
+    snears = Tuple{Int, Vector{Int}}[]
+    for near in interactions
+        if length(snears)!= 0 && snears[end][1] == near[1] 
+            push!(snears[end][2], near[2])
+        else
+            push!(snears, (near[1], [near[2]]))
+        end
+    end
+
+    _foreach = multithreading ? ThreadsX.foreach : Base.foreach
+    blocks = Vector{BlockSparseMatrices.DenseMatrixBlock{K, Matrix{K}, Vector{I}}}(
+        undef, length(snears)
+    )
+
+    _foreach(enumerate(snears)) do (idx, snear) 
+        testidcs = value(tree.test_cluster, snear[1])
+        trialidcs = value(tree.trial_cluster, snear[2])
+        matrix = zeros(K, length(testidcs), length(trialidcs))
+        assembler(matrix, testidcs, trialidcs)
+
+        blocks[idx] =  BlockSparseMatrices.DenseMatrixBlock{K, Matrix{K}, Vector{Int}}(
+            matrix,
+            testidcs,
+            trialidcs
+        )
+    end
+
+    return BlockSparseMatrix(
+        blocks, (tree.test_cluster.num_elements, tree.trial_cluster.num_elements)
+    )
+end