new grid glueing algorithm using BinnedPointList

src/binnedpointlist.jl

@@ -213,6 +213,16 @@ function Base.insert!(bpl::BinnedPointList{T}, p) where {T}
     end
 end
 
+function findpoint(bpl::BinnedPointList{T}, p) where {T}
+    _rebin_all_points!(bpl)
+    _bin_of_point!(bpl, p)
+    if reduce(*, bpl.current_bin) > 0
+        return _findpoint(bpl, bpl.bins[bpl.current_bin...], p)
+    else
+        return _findpoint(bpl, bpl.unbinned, p)
+    end
+end
+
 Base.insert!(bpl::BinnedPointList{T}, x::Number) where {T} = insert!(bpl, (x))
 Base.insert!(bpl::BinnedPointList{T}, x::Number, y::Number) where {T} = insert!(bpl, (x, y))
 Base.insert!(bpl::BinnedPointList{T}, x::Number, y::Number, z::Number) where {T} = insert!(bpl, (x, y, z))

src/simplexgrid.jl

@@ -10,7 +10,7 @@ function simplexgrid(coord::Array{Tc,2},
 
 Create  d-dimensional simplex grid from five arrays.
 
--    coord: d ``\\times`` n_points matrix of coordinates
+-        coord: d ``\\times`` n_points matrix of coordinates
 -    cellnodes: d+1 ``\\times`` n_tri matrix of triangle - point incidence
 -  cellregions: n_tri vector of cell region markers
 -   bfacenodes: d ``\\times`` n_bf matrix of boundary facet - point incidences
@@ -749,23 +749,79 @@ function glue(g1::ExtendableGrid,g2::ExtendableGrid;
       end
     end
 
+    function barycenter!(bc,bfnodes,coord,iface)
+        for idim=1:dim
+            bc[idim]=0.0
+            for jdim=1:dim
+                bc[idim]+=coord[idim,bfnodes[jdim,iface]]/dim
+            end
+        end
+        bc
+    end
+
+    nbfreg1=num_bfaceregions(g1)
+    nbfreg2=num_bfaceregions(g2)
+    breg1used=zeros(Bool,nbfreg1)
+    breg2used=zeros(Bool,nbfreg2)
+    for reg in g1regions
+        breg1used[reg]=true
+    end
+    for reg in g2regions
+        breg2used[reg]=true
+    end
 
-    use_region(ireg,regionlist) = findfirst(i->i==ireg,regionlist)!=nothing
-    # Run over all pairs of boundary faces and try to match them
-    for ibf1=1:nbf1
-        if use_region(bfreg1[ibf1],g1regions)
-            for ibf2=1:nbf2
-                if use_region(bfreg2[ibf2],g2regions)
+    if false
+        # Run over all pairs of boundary faces and try to match them
+        # This can scale catastrophically...
+        for ibf1=1:nbf1
+            if breg1used[bfreg1[ibf1]]
+                for ibf2=1:nbf2
+                    if breg2used[bfreg2[ibf2]]
+                        match_faces(ibf1,ibf2)
+                    end
+                end
+            end
+        end
+    else
+        # Use a binned point list for speed up.
+        # Facets are found through their barycenters.
+        bcenter=zeros(dim)
+        bpl=BinnedPointList(eltype(coord1),dim)
+        nbf1used=0
+
+        for ibf1=1:nbf1
+            if breg1used[bfreg1[ibf1]]
+                nbf1used+=1
+            end
+        end
+
+        # marker for facet numbers
+        bf1used=zeros(Int,nbf1used)
+
+        # insert all barycenters of used bfaces from grid 1 into the binned pointlist
+        for ibf1=1:nbf1
+            if breg1used[bfreg1[ibf1]]
+                i=insert!(bpl,barycenter!(bcenter,bfn1,coord1,ibf1))
+                bf1used[i]=ibf1
+            end
+        end
+
+        # loop over used grid2 bfaces and find if barycenters match with a point in
+        # the binned pointlist
+        for ibf2=1:nbf2
+            if breg2used[bfreg2[ibf2]]
+                i=findpoint(bpl,barycenter!(bcenter,bfn2,coord2,ibf2))
+                # if a match has been found, identify facets
+                if i>0
+                    ibf1=bf1used[i]
+                    # use "old" matching, double check at once
                     match_faces(ibf1,ibf2)
                 end
             end
         end
     end
 
     @info "glue: matches found: nodes $(n_matching_nodes) bfaces: $(n_matching_faces)"
-
-
-
 
     creg1=g1[CellRegions]
     creg2=g2[CellRegions]