Introduce unique entities in order to reduce redundant communication …

…and loops.

docs/Project.toml

@@ -1,10 +1,12 @@
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 FerriteDistributed = "570c3397-5fe4-4792-be0d-48dbf0d14605"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"

ext/FerriteDistributedHYPREAssembly.jl

@@ -5,7 +5,7 @@ module FerriteDistributedHYPREAssembly
 
 using FerriteDistributed
 # TODO remove me. These are merely hotfixes to split the extensions trasiently via an internal API.
-import FerriteDistributed: getglobalgrid, num_local_true_dofs, num_local_dofs, global_comm, interface_comm, global_rank, compute_owner, remote_entities, num_fields
+import FerriteDistributed: getglobalgrid, num_local_true_dofs, num_local_dofs, global_comm, interface_comm, global_rank, compute_owner, remote_entities, num_fields, local_entities
 using MPI
 using HYPRE
 using Base: @propagate_inbounds

ext/FerriteDistributedHYPREAssembly/conversion.jl

@@ -27,50 +27,53 @@ function FerriteDistributed.extract_local_part!(u::Vector{T}, uh::HYPREVector, d
     # TODO speed this up and better API
     dgrid = getglobalgrid(dh)
     for sv ∈ get_shared_vertices(dgrid)
-        lvi = sv.local_idx
         my_rank != compute_owner(dgrid, sv) && continue
-        for field_idx in 1:num_fields(dh)
-            if has_vertex_dofs(dh, field_idx, lvi)
-                local_dofs = vertex_dofs(dh, field_idx, lvi)
-                global_dofs = dh.ldof_to_gdof[local_dofs]
-                for receiver_rank ∈ keys(remote_entities(sv))
-                    for i ∈ 1:length(global_dofs)
-                        # Note that u already has the correct values for all locally owned dofs due to the loop above!
-                        gdof_value_send[receiver_rank][global_dofs[i]] = u[local_dofs[i]]
+        for lvi ∈ local_entities(sv)
+            for field_idx in 1:num_fields(dh)
+                if has_vertex_dofs(dh, field_idx, lvi)
+                    local_dofs = vertex_dofs(dh, field_idx, lvi)
+                    global_dofs = dh.ldof_to_gdof[local_dofs]
+                    for receiver_rank ∈ keys(remote_entities(sv))
+                        for i ∈ 1:length(global_dofs)
+                            # Note that u already has the correct values for all locally owned dofs due to the loop above!
+                            gdof_value_send[receiver_rank][global_dofs[i]] = u[local_dofs[i]]
+                        end
                     end
                 end
             end
         end
     end
 
     for se ∈ get_shared_edges(dgrid)
-        lei = se.local_idx
         my_rank != compute_owner(dgrid, se) && continue
-        for field_idx in 1:num_fields(dh)
-            if has_edge_dofs(dh, field_idx, lei)
-                local_dofs = edge_dofs(dh, field_idx, lei)
-                global_dofs = dh.ldof_to_gdof[local_dofs]
-                for receiver_rank ∈ keys(remote_entities(se))
-                    for i ∈ 1:length(global_dofs)
-                        # Note that u already has the correct values for all locally owned dofs due to the loop above!
-                        gdof_value_send[receiver_rank][global_dofs[i]] = u[local_dofs[i]]
+        for lei ∈ local_entities(se)
+            for field_idx in 1:num_fields(dh)
+                if has_edge_dofs(dh, field_idx, lei)
+                    local_dofs = edge_dofs(dh, field_idx, lei)
+                    global_dofs = dh.ldof_to_gdof[local_dofs]
+                    for receiver_rank ∈ keys(remote_entities(se))
+                        for i ∈ 1:length(global_dofs)
+                            # Note that u already has the correct values for all locally owned dofs due to the loop above!
+                            gdof_value_send[receiver_rank][global_dofs[i]] = u[local_dofs[i]]
+                        end
                     end
                 end
             end
         end
     end
 
     for sf ∈ get_shared_faces(dgrid)
-        lfi = sf.local_idx
         my_rank != compute_owner(dgrid, sf) && continue
-        for field_idx in 1:num_fields(dh)
-            if has_face_dofs(dh, field_idx, lfi)
-                local_dofs = face_dofs(dh, field_idx, lfi)
-                global_dofs = dh.ldof_to_gdof[local_dofs]
-                for receiver_rank ∈ keys(remote_entities(sf))
-                    for i ∈ 1:length(global_dofs)
-                        # Note that u already has the correct values for all locally owned dofs due to the loop above!
-                        gdof_value_send[receiver_rank][global_dofs[i]] = u[local_dofs[i]]
+        for lfi ∈ local_entities(sf)
+            for field_idx in 1:num_fields(dh)
+                if has_face_dofs(dh, field_idx, lfi)
+                    local_dofs = face_dofs(dh, field_idx, lfi)
+                    global_dofs = dh.ldof_to_gdof[local_dofs]
+                    for receiver_rank ∈ keys(remote_entities(sf))
+                        for i ∈ 1:length(global_dofs)
+                            # Note that u already has the correct values for all locally owned dofs due to the loop above!
+                            gdof_value_send[receiver_rank][global_dofs[i]] = u[local_dofs[i]]
+                        end
                     end
                 end
             end

ext/FerriteDistributedPartitionedArrays.jl

@@ -6,7 +6,7 @@ module FerriteDistributedPartitionedArrays
 using FerriteDistributed
 # TODO remove me. These are merely hotfixes to split the extensions trasiently via an internal API.
 import FerriteDistributed: getglobalgrid, num_global_dofs, num_local_true_dofs, num_local_dofs, global_comm, interface_comm, global_rank, compute_owner, remote_entities,
-    num_fields, getfieldnames, getdim
+    num_fields, getfieldnames, getdim, local_entities
 using MPI
 using PartitionedArrays
 using Base: @propagate_inbounds

ext/FerriteDistributedPartitionedArrays/assembler.jl

@@ -106,66 +106,31 @@ struct COOAssembler{T}
         ghost_dof_field_index_to_send = [Int[] for i ∈ 1:destination_len]
         ghost_dof_owner = [Int[] for i ∈ 1:destination_len] # corresponding owner
         ghost_dof_pivot_to_send = [Int[] for i ∈ 1:destination_len] # corresponding dof to interact with
-        for (pivot_vertex, pivot_shared_vertex) ∈ dgrid.shared_vertices
+        for pivot_shared_vertex ∈ get_shared_vertices(dgrid)
             # Start by searching shared entities which are not owned
             pivot_vertex_owner_rank = compute_owner(dgrid, pivot_shared_vertex)
-            pivot_cell_idx = pivot_vertex[1]
-
-            if my_rank != pivot_vertex_owner_rank
-                sender_slot = destination_index[pivot_vertex_owner_rank]
-
-                Ferrite.@debug println("$pivot_vertex may require synchronization (R$my_rank)")
-                # Note: We have to send ALL dofs on the element to the remote.
-                cell_dofs_upper_bound = (pivot_cell_idx == getncells(dh.grid)) ? length(dh.cell_dofs) : dh.cell_dofs_offset[pivot_cell_idx+1]
-                cell_dofs = dh.cell_dofs[dh.cell_dofs_offset[pivot_cell_idx]:cell_dofs_upper_bound]
-
-                for (field_idx, field_name) in zip(1:num_fields(dh), getfieldnames(dh))
-                    !has_vertex_dofs(dh, field_idx, pivot_vertex) && continue
-                    pivot_vertex_dofs = vertex_dofs(dh, field_idx, pivot_vertex)
-
-                    for d ∈ 1:dh.field_dims[field_idx]
-                        Ferrite.@debug println("  adding dof $(pivot_vertex_dofs[d]) to ghost sync synchronization on slot $sender_slot (R$my_rank)")
-
-                        # Extract dofs belonging to the current field
-                        #cell_field_dofs = cell_dofs[dof_range(dh, field_name)]
-                        #for cell_field_dof ∈ cell_field_dofs
-                        for cell_dof ∈ cell_dofs
-                            append!(ghost_dof_pivot_to_send[sender_slot], ldof_to_gdof[pivot_vertex_dofs[d]])
-                            append!(ghost_dof_to_send[sender_slot], ldof_to_gdof[cell_dof])
-                            append!(ghost_rank_to_send[sender_slot], ldof_to_rank[cell_dof])
-                            append!(ghost_dof_field_index_to_send[sender_slot], field_idx)
-                        end
-                    end
-                end
-            end
-        end
-
-        if dim > 1
-            for (pivot_face, pivot_shared_face) ∈ dgrid.shared_faces
-                # Start by searching shared entities which are not owned
-                pivot_face_owner_rank = compute_owner(dgrid, pivot_shared_face)
-                pivot_cell_idx = pivot_face[1]
+            for pivot_vi ∈ local_entities(pivot_shared_vertex)
+                pivot_cell_idx = pivot_vi[1]
+                if my_rank != pivot_vertex_owner_rank
+                    sender_slot = destination_index[pivot_vertex_owner_rank]
 
-                if my_rank != pivot_face_owner_rank
-                    sender_slot = destination_index[pivot_face_owner_rank]
-
-                    Ferrite.@debug println("$pivot_face may require synchronization (R$my_rank)")
+                    Ferrite.@debug println("$pivot_shared_vertex may require synchronization (R$my_rank)")
                     # Note: We have to send ALL dofs on the element to the remote.
                     cell_dofs_upper_bound = (pivot_cell_idx == getncells(dh.grid)) ? length(dh.cell_dofs) : dh.cell_dofs_offset[pivot_cell_idx+1]
                     cell_dofs = dh.cell_dofs[dh.cell_dofs_offset[pivot_cell_idx]:cell_dofs_upper_bound]
 
                     for (field_idx, field_name) in zip(1:num_fields(dh), getfieldnames(dh))
-                        !has_face_dofs(dh, field_idx, pivot_face) && continue
-                        pivot_face_dofs = face_dofs(dh, field_idx, pivot_face)
+                        !has_vertex_dofs(dh, field_idx, pivot_vi) && continue
+                        pivot_vertex_dofs = vertex_dofs(dh, field_idx, pivot_vi)
 
                         for d ∈ 1:dh.field_dims[field_idx]
-                            Ferrite.@debug println("  adding dof $(pivot_face_dofs[d]) to ghost sync synchronization on slot $sender_slot (R$my_rank)")
+                            Ferrite.@debug println("  adding dof $(pivot_vertex_dofs[d]) to ghost sync synchronization on slot $sender_slot (R$my_rank)")
 
                             # Extract dofs belonging to the current field
                             #cell_field_dofs = cell_dofs[dof_range(dh, field_name)]
                             #for cell_field_dof ∈ cell_field_dofs
                             for cell_dof ∈ cell_dofs
-                                append!(ghost_dof_pivot_to_send[sender_slot], ldof_to_gdof[pivot_face_dofs[d]])
+                                append!(ghost_dof_pivot_to_send[sender_slot], ldof_to_gdof[pivot_vertex_dofs[d]])
                                 append!(ghost_dof_to_send[sender_slot], ldof_to_gdof[cell_dof])
                                 append!(ghost_rank_to_send[sender_slot], ldof_to_rank[cell_dof])
                                 append!(ghost_dof_field_index_to_send[sender_slot], field_idx)
@@ -176,31 +141,33 @@ struct COOAssembler{T}
             end
         end
 
-        if dim > 2
-            for (pivot_edge, pivot_shared_edge) ∈ dgrid.shared_edges
+        if dim > 1
+            for pivot_shared_face ∈ get_shared_faces(dgrid)
                 # Start by searching shared entities which are not owned
-                pivot_edge_owner_rank = compute_owner(dgrid, pivot_shared_edge)
-                pivot_cell_idx = pivot_edge[1]
+                pivot_face_owner_rank = compute_owner(dgrid, pivot_shared_face)
+                pivot_fi = local_entities(pivot_shared_face)[1]
+                pivot_cell_idx = pivot_fi[1]
 
-                if my_rank != pivot_edge_owner_rank
-                    sender_slot = destination_index[pivot_edge_owner_rank]
+                if my_rank != pivot_face_owner_rank
+                    sender_slot = destination_index[pivot_face_owner_rank]
 
-                    Ferrite.@debug println("$pivot_edge may require synchronization (R$my_rank)")
+                    Ferrite.@debug println("$pivot_shared_face may require synchronization (R$my_rank)")
                     # Note: We have to send ALL dofs on the element to the remote.
                     cell_dofs_upper_bound = (pivot_cell_idx == getncells(dh.grid)) ? length(dh.cell_dofs) : dh.cell_dofs_offset[pivot_cell_idx+1]
                     cell_dofs = dh.cell_dofs[dh.cell_dofs_offset[pivot_cell_idx]:cell_dofs_upper_bound]
 
                     for (field_idx, field_name) in zip(1:num_fields(dh), getfieldnames(dh))
-                        !has_edge_dofs(dh, field_idx, pivot_edge) && continue
-                        pivot_edge_dofs = edge_dofs(dh, field_idx, pivot_edge)
+                        !has_face_dofs(dh, field_idx, pivot_fi) && continue
+                        pivot_face_dofs = face_dofs(dh, field_idx, pivot_fi)
 
                         for d ∈ 1:dh.field_dims[field_idx]
-                            Ferrite.@debug println("  adding dof $(pivot_edge_dofs[d]) to ghost sync synchronization on slot $sender_slot (R$my_rank)")
+                            Ferrite.@debug println("  adding dof $(pivot_face_dofs[d]) to ghost sync synchronization on slot $sender_slot (R$my_rank)")
+
                             # Extract dofs belonging to the current field
                             #cell_field_dofs = cell_dofs[dof_range(dh, field_name)]
                             #for cell_field_dof ∈ cell_field_dofs
                             for cell_dof ∈ cell_dofs
-                                append!(ghost_dof_pivot_to_send[sender_slot], ldof_to_gdof[pivot_edge_dofs[d]])
+                                append!(ghost_dof_pivot_to_send[sender_slot], ldof_to_gdof[pivot_face_dofs[d]])
                                 append!(ghost_dof_to_send[sender_slot], ldof_to_gdof[cell_dof])
                                 append!(ghost_rank_to_send[sender_slot], ldof_to_rank[cell_dof])
                                 append!(ghost_dof_field_index_to_send[sender_slot], field_idx)
@@ -211,6 +178,43 @@ struct COOAssembler{T}
             end
         end
 
+        if dim > 2
+            for pivot_shared_edge ∈ get_shared_edges(dgrid)
+                # Start by searching shared entities which are not owned
+                pivot_edge_owner_rank = compute_owner(dgrid, pivot_shared_edge)
+                for pivot_ei ∈ local_entities(pivot_shared_edge)
+                    pivot_cell_idx = pivot_ei[1]
+
+                    if my_rank != pivot_edge_owner_rank
+                        sender_slot = destination_index[pivot_edge_owner_rank]
+
+                        Ferrite.@debug println("$pivot_shared_edge may require synchronization (R$my_rank)")
+                        # Note: We have to send ALL dofs on the element to the remote.
+                        cell_dofs_upper_bound = (pivot_cell_idx == getncells(dh.grid)) ? length(dh.cell_dofs) : dh.cell_dofs_offset[pivot_cell_idx+1]
+                        cell_dofs = dh.cell_dofs[dh.cell_dofs_offset[pivot_cell_idx]:cell_dofs_upper_bound]
+
+                        for (field_idx, field_name) in zip(1:num_fields(dh), getfieldnames(dh))
+                            !has_edge_dofs(dh, field_idx, pivot_ei) && continue
+                            pivot_edge_dofs = edge_dofs(dh, field_idx, pivot_ei)
+
+                            for d ∈ 1:dh.field_dims[field_idx]
+                                Ferrite.@debug println("  adding dof $(pivot_edge_dofs[d]) to ghost sync synchronization on slot $sender_slot (R$my_rank)")
+                                # Extract dofs belonging to the current field
+                                #cell_field_dofs = cell_dofs[dof_range(dh, field_name)]
+                                #for cell_field_dof ∈ cell_field_dofs
+                                for cell_dof ∈ cell_dofs
+                                    append!(ghost_dof_pivot_to_send[sender_slot], ldof_to_gdof[pivot_edge_dofs[d]])
+                                    append!(ghost_dof_to_send[sender_slot], ldof_to_gdof[cell_dof])
+                                    append!(ghost_rank_to_send[sender_slot], ldof_to_rank[cell_dof])
+                                    append!(ghost_dof_field_index_to_send[sender_slot], field_idx)
+                                end
+                            end
+                        end
+                    end
+                end
+            end
+        end
+
         ghost_send_buffer_lengths = Int[length(i) for i ∈ ghost_dof_to_send]
         ghost_recv_buffer_lengths = zeros(Int, destination_len)
         MPI.Neighbor_alltoall!(UBuffer(ghost_send_buffer_lengths,1), UBuffer(ghost_recv_buffer_lengths,1), interface_comm(dgrid));

src/Entity.jl

@@ -0,0 +1,88 @@
+"""
+    VertexRepresentation
+
+We can identify an edge uniquely by the sorted node numbers associated with the end points.
+"""
+struct VertexRepresentation
+    node::Int
+end
+
+"""
+    EdgeRepresentation
+
+We can identify an edge uniquely by the sorted node numbers associated with the end points.
+"""
+struct EdgeRepresentation
+    a::Int
+    b::Int
+end
+
+function EdgeRepresentation(ab::Tuple{Int,Int}) 
+    se = Ferrite.sortedge_fast(ab)
+    return EdgeRepresentation(se[1], se[2])
+end
+
+"""
+    FaceRepresentation
+
+We can identify a face uniquely by 3 sorted node numbers associated with the vertices.
+3 points are sufficient, because 3 (non-aligned) points can uniquely describe a surface in 3D.
+"""
+struct FaceRepresentation
+    a::Int
+    b::Int
+    c::Int
+end
+
+
+function FaceRepresentation(ab::Tuple{Int,Int})
+    @warn "Fixme after https://github.com/Ferrite-FEM/Ferrite.jl/pull/789" maxlog=1
+    sf = Ferrite.sortface_fast(ab)
+    return FaceRepresentation(-1, sf[1], sf[2])
+end
+
+function FaceRepresentation(abc::Union{Tuple{Int,Int,Int}, Tuple{Int,Int,Int,Int}}) 
+    sf = Ferrite.sortface_fast(abc)
+    return FaceRepresentation(sf[1], sf[2], sf[3])
+end
+
+
+"""
+    Entity
+
+Supertype for geometric entities.
+"""
+abstract type Entity end
+
+# !!!THE STUFF BELOW IS CURRENTLY UNUSED!!!
+
+"""
+    Vertex <: Entity
+
+A shared vertex induced by a local vertex index and all remote vertex indices on all remote ranks.
+"""
+struct Vertex <: Entity
+    unique_local_representation::VertexRepresentation # Identify via node in grid
+    local_vertices::Vector{VertexIndex}
+end
+
+"""
+    Face <: Entity
+
+A face induced by a local face index and all remote face indices on all remote ranks.
+"""
+struct Face <: Entity
+    unique_representation::FaceRepresentation
+    local_faces::Pair{FaceIndex,FaceIndex}
+end
+
+"""
+    Edge <: Entity
+
+An edge induced by a local edge index and all remote edge indices on all remote ranks.
+"""
+struct Edge <: Entity
+    unique_representation::EdgeRepresentation # Identify via node in grid
+    local_edges::Vector{EdgeIndex}
+end
+

src/FerriteDistributed.jl

@@ -8,7 +8,7 @@ import Ferrite: get_coordinate_eltype, ScalarWrapper, @debug,
     nnodes_per_cell, n_components, get_grid, getdim,
     BoundaryIndex, FaceIndex, EdgeIndex, CellIndex, VertexIndex,
     AbstractTopology, EntityNeighborhood,
-    AbstractCell, boundaryfunction, faces, edges, vertices, nvertices, nfaces, nedges,
+    AbstractGrid, AbstractCell, boundaryfunction, faces, edges, vertices, nvertices, nfaces, nedges,
     cellnodes!, cellcoords!,
     getfieldnames, getfieldinterpolation, default_interpolation,
     reference_coordinates, value, getrefshape, dof_range, getfielddim,
@@ -18,6 +18,7 @@ include("utils.jl")
 
 include("CoverTopology.jl")
 
+include("Entity.jl")
 include("SharedEntity.jl")
 
 include("interface.jl")