Added is_partition

Author: JordiManyer

Project.toml

@@ -14,6 +14,7 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
@@ -49,6 +50,7 @@ FastGaussQuadrature = "0.4.5, 1"
 FileIO = "1.5"
 FillArrays = "1.11"
 ForwardDiff = "0.10.14, 1"
+Graphs = "1"
 JLD2 = "0.5"
 JSON = "0.21"
 LineSearches = "7.0.1"

src/Adaptivity/PolytopalRefinement.jl

@@ -8,6 +8,7 @@
 # are correctly oriented.
 
 function coarsen(model::Geometry.PolytopalDiscreteModel,ptopo::Geometry.PatchTopology; return_glue=false)
+  @check Geometry.is_partition(ptopo) "The patch topology is not a valid partition of the model"
   new_polys, new_connectivity = generate_patch_polytopes(model,ptopo)
 
   vertex_coords = get_vertex_coordinates(get_grid_topology(model))
@@ -84,7 +85,7 @@ function generate_patch_polytopes(
 
   new_connectivity = Vector{Vector{Int32}}(undef, npatches)
   new_polys = Vector{GeneralPolytope{D}}(undef, npatches)
-  Threads.@threads for patch in 1:npatches
+  for patch in 1:npatches
     cells = view(patch_cells,patch)
     if isone(length(cells))
       new_polys[patch] = polys[first(cells)]

src/Geometry/Geometry.jl

@@ -12,6 +12,8 @@ using FillArrays
 using LinearAlgebra: ⋅
 using Statistics: mean
 using DataStructures: SortedSet
+using SparseArrays: sparse
+using Graphs: Graph, is_connected
 
 using Gridap.Helpers
 using Gridap.Arrays

src/Geometry/GridTopologies.jl

@@ -556,6 +556,75 @@ function restrict(topo::GridTopology,parent_cell_to_mask::AbstractVector{Bool})
   restrict(topo, cell_to_parent_cell)
 end
 
+"""
+    compute_graph(topo::GridTopology, d_node::Integer, d_edge::Integer; self_loops=true, Tv=Int8)
+
+Computes the adjacency graph between the `d_node` entities w.r.t the `d_edge` entities. The 
+graph is returned as a sparse matrix, where nonzero entries indicate an adjacency between the nodes.
+
+For instance for a D-dimensional topology and pairs `(d_node,d_edge)` we have: 
+
+- `(1,0)`: graph of edge-to-edge adjacencies where two edges (dim 1) are connected if they share a node (dim 0).
+- `(D,D-1)`: graph of cell-to-cell adjacencies where two cells (dim D) are connected if they share a facet (dim D-1).
+
+Weights are defaulted to `one(Tv)`. If `self_loops` is `true`, self-loops (diagonal entries) are included in the graph.
+
+"""
+function compute_graph(topo::GridTopology, d_node::Integer, d_edge::Integer; self_loops=true, Tv=Int8)
+  edge_to_nodes = Geometry.get_faces(topo, d_edge, d_node)
+  node_to_edges = Geometry.get_faces(topo, d_node, d_edge)
+  compute_graph(node_to_edges, edge_to_nodes; self_loops, Tv)
+end
+
+function compute_graph(
+  node_to_edges::Table,edge_to_nodes::Table;Tv=Int8,self_loops=false
+)
+  nodes = eachindex(node_to_edges)
+  node_to_lnode = eachindex(node_to_edges)
+  return compute_graph(
+    nodes,node_to_edges,edge_to_nodes;node_to_lnode,Tv,self_loops
+  )
+end
+
+function compute_graph(
+  nodes,node_to_edges::Table,edge_to_nodes::Table;
+  node_to_lnode = Dict{Int32,Int32}(n => i for (i,n) in enumerate(nodes)),
+  self_loops=true,
+  Tv = Int8,
+)
+  ndata = 0
+  for node in nodes
+    ndata += self_loops
+    for edge in view(node_to_edges,node)
+      for nbor in view(edge_to_nodes,edge)
+        ndata += (nbor != node) && (nbor in nodes)
+      end
+    end
+  end
+  I = zeros(eltype(nodes),ndata)
+  J = zeros(eltype(nodes),ndata)
+  p = 0
+  for (lnode,node) in enumerate(nodes)
+    if self_loops
+      p += 1
+      I[p] = lnode
+      J[p] = lnode
+    end
+    for edge in view(node_to_edges,node)
+      for nbor in view(edge_to_nodes,edge)
+        if (nbor != node) && (nbor in nodes)
+          p += 1
+          I[p] = lnode
+          J[p] = node_to_lnode[nbor]
+        end
+      end
+    end
+  end
+  V = ones(Tv,ndata)
+  n = length(nodes)
+  return sparse(I,J,V,n,n)
+end
+
 # Helpers
 
 function  _compute_cell_perm_indices!(

src/Geometry/PatchTriangulations.jl

@@ -101,6 +101,79 @@ function patch_extend(ptopo::PatchTopology{Dc},patch_to_data,Df=Dc) where Dc
   return pface_to_data
 end
 
+"""
+    is_cover(ptopo::PatchTopology)
+    is_cover(topo::GridTopology, patch_cells::Table)
+
+Returns `true` if the given patch topology is a cover of the underlying topology, i.e 
+if every cell in the topology is contained in at least one patch.
+"""
+function is_cover(ptopo::PatchTopology)
+  patch_cells = Geometry.get_patch_cells(ptopo)
+  is_cover(ptopo.topo, patch_cells)
+end
+
+function is_cover(topo::GridTopology, patch_cells)
+  cache = array_cache(patch_cells)
+  cell_is_covered = fill(false, num_cells(topo))
+  for patch in eachindex(patch_cells)
+    cells = getindex!(cache, patch_cells, patch)
+    cell_is_covered[cells] .= true
+  end
+  return all(cell_is_covered)
+end
+
+"""
+  is_partition(ptopo::PatchTopology; kwargs...)
+  is_partition(topo::GridTopology, patch_cells::Table; fail_fast = true)
+
+Check if the given patch topology is a valid partition of the underlying topology.
+
+To be a valid partition, the patches
+
+  - must cover the whole topology
+  - must be disjoint (i.e non-overlapping)
+  - must be connected (a patch cannot be split)
+
+If `fail_fast` is `true`, the function will exit as soon as it finds a patch that is not connected.
+Otherwise, it will check all patches and print a warning with the indices of the bad patches.
+
+"""
+function is_partition(ptopo::PatchTopology; kwargs...)
+  patch_cells = Geometry.get_patch_cells(ptopo)
+  is_partition(ptopo.topo, patch_cells; kwargs...)
+end
+
+function is_partition(topo::GridTopology, patch_cells; fail_fast = true)
+  if !is_cover(topo, patch_cells) || !isequal(sum(length, patch_cells), num_cells(topo))
+    return false # Not a cover or not disjoint
+  end
+
+  # Check patch connectivity
+  # We could return false at the first bad patch, but I'd rather get 
+  # all bad patch indices for debugging purposes.
+  D = num_cell_dims(topo)
+  cell_to_facets = get_faces(topo,D,D-1)
+  facet_to_cells = get_faces(topo,D-1,D)
+
+  cache = array_cache(patch_cells)
+  bad_patches = Int[]
+  for patch in eachindex(patch_cells)
+    cells = getindex!(cache, patch_cells, patch)
+    G = Graph(compute_graph(Set(cells), cell_to_facets, facet_to_cells))
+    if !is_connected(G)
+      push!(bad_patches, patch)
+      fail_fast && return false # Stop at the first bad patch
+    end
+  end
+
+  if !isempty(bad_patches)
+    @warn "The following patches are not connected: $(bad_patches)"
+    return false
+  end
+  return true
+end
+
 function generate_patch_faces(ptopo::PatchTopology{Dc},Df) where Dc
   cell_to_faces = get_faces(ptopo.topo,Dc,Df)
   patch_cells = get_patch_cells(ptopo)