Merge pull request #20 from gridap/background_tet_mesh

Background tetrahedral grid

Project.toml

@@ -23,7 +23,7 @@ Downloads = "1.4"
 FileIO = "1.6"
 FillArrays = "0.11, 0.13"
 GraphRecipes = "0.5"
-Gridap = "0.17"
+Gridap = "0.17.16"
 GridapEmbedded = "0.8"
 IterativeSolvers = "0.9"
 MeshIO = "0.4"

src/Embedded.jl

@@ -73,37 +73,38 @@ function _cut_stl(model::DiscreteModel,geom::STLGeometry;kwargs...)
   (bgcell_to_ioc,subcells,cell_to_io,subfacets,face_to_io,oid_to_ls),bgface_to_ioc
 end
 
+function send_to_ref_space(
+  model::DiscreteModel,
+  cell_to_bgcell::Vector,
+  subgrid::Grid)
+
+  send_to_ref_space(get_grid(model),cell_to_bgcell,subgrid)
+end
+
 function send_to_ref_space(grid::Grid,cell_to_bgcell::Vector,subgrid::Grid)
-  send_to_ref_space(
-    grid,
-    cell_to_bgcell,
-    get_cell_node_ids(subgrid),
-    get_node_coordinates(subgrid))
+  bgcell_map = get_cell_map(grid)
+  bgcell_invmap = lazy_map(inverse_map,bgcell_map)
+  cell_invmap = lazy_map(Reindex(bgcell_invmap),cell_to_bgcell)
+  cell_nodes = get_cell_node_ids(subgrid)
+  node_coords = get_node_coordinates(subgrid)
+  node_cell = _first_inverse_index_map(cell_nodes,num_nodes(subgrid))
+  node_invmap = lazy_map(Reindex(cell_invmap),node_cell)
+  node_rcoords = lazy_map(evaluate,node_invmap,node_coords)
+  _collect(node_rcoords)
 end
 
-function send_to_ref_space(
-  grid::Grid,
-  cell_to_bgcell::Vector,
-  cell_to_points::AbstractVector,
-  points::Vector)
-
-  rpoints = zero(points)
-  cache = array_cache(cell_to_points)
-  bgcache = array_cache(get_cell_node_ids(grid))
-  for i in 1:length(cell_to_points)
-    for p in getindex!(cache,cell_to_points,i)
-      rpoints[p] = send_to_ref_space!(bgcache,grid,cell_to_bgcell[i],points[p])
+function _first_inverse_index_map(a_to_b,nb)
+  na = length(a_to_b)
+  T = eltype(eltype(a_to_b))
+  b_to_a = ones(T,nb)
+  c = array_cache(a_to_b)
+  for a in 1:na
+    bs = getindex!(c,a_to_b,a)
+    for b in bs
+      b_to_a[b] = a
     end
   end
-  rpoints
-end
-
-function send_to_ref_space!(cache,grid::Grid,cell::Integer,point::Point)
-  cell_nodes = getindex!(cache,get_cell_node_ids(grid),cell)
-  node_coordinates = get_node_coordinates(grid)
-  pmin = node_coordinates[cell_nodes[1]]
-  pmax = node_coordinates[cell_nodes[end]]
-  Point(Tuple(point-pmin)./Tuple(pmax-pmin))
+  b_to_a
 end
 
 function _normals(geom::STLGeometry,face_to_stlface)
@@ -130,3 +131,13 @@ function check_requisites(geo::STLGeometry,bgmodel::DiscreteModel;kwargs...)
   stl = get_stl(geo)
   check_requisites(stl,bgmodel)
 end
+
+function _collect(a::LazyArray)
+  c = array_cache(a)
+  T = eltype(a)
+  b = zeros(T,size(a))
+  for i in eachindex(a)
+    b[i] = getindex!(c,a,i)
+  end
+  b
+end

src/STLs.jl

@@ -223,6 +223,7 @@ function get_cell!(c,grid::Grid,i)
   p = get_polytope( only(get_reffes(grid)) )
   get_cell!(c,T,X,p,i)
 end
+
 function get_facet_cache(model::DiscreteModel)
   Dc = num_dims(model)
   get_dface_cache(model,Dc-1)
@@ -240,9 +241,11 @@ end
 
 function get_dface!(c,model::DiscreteModel,i,::Val{d}) where d
   topo = get_grid_topology(model)
+  Dc = num_dims(model)
+  PT = ExtrusionPolytope{Dc}
   T = get_face_vertices(topo,d)
   X = get_vertex_coordinates(topo)
-  p = get_polytope( only(get_reffes(model)) )
+  p = get_polytope( only(get_reffes(model)) )::PT
   get_dface!(c,T,X,p,i,Val{d}())
 end
 
@@ -867,6 +870,18 @@ end
 
 min_height(model::DiscreteModel) = min_height(get_grid(model))
 
+function max_length(grid::Grid)
+  max_len = 0.0
+  c = get_cell_cache(grid)
+  for i in 1:num_cells(grid)
+    cell = get_cell!(c,grid,i)
+    max_len = max(max_len,max_length(cell))
+  end
+  max_len
+end
+
+max_length(model::DiscreteModel) = max_length(get_grid(model))
+
 # Cartesian Grid
 
 function get_bounding_box(grid::CartesianGrid)
@@ -886,8 +901,13 @@ function measure(a::CartesianGrid{D,T,typeof(identity)},mask) where {D,T}
   cell_vol * count(mask)
 end
 
+function compute_cell_to_facets(model_a::DiscreteModel,grid_b)
+  grid_a = get_grid(model_a)
+  compute_cell_to_facets(grid_a,grid_b)
+end
 
-function compute_cell_to_facets(grid::CartesianGrid,stl::STL)
+
+function compute_cell_to_facets(grid::CartesianGrid,stl)
   CELL_EXPANSION_FACTOR = 1e-3
   desc = get_cartesian_descriptor(grid)
   @assert length(get_reffes(grid)) == 1
@@ -920,21 +940,120 @@ function compute_cell_to_facets(grid::CartesianGrid,stl::STL)
       end
     end
   end
-  cell_to_stl_facets = [ Int32[] for _ in 1:num_cells(grid) ]
-  for (cells,stl_facets) in zip(thread_to_cells,thread_to_stl_facets)
-    for (cell,stl_facet) in zip(cells,stl_facets)
-      push!(cell_to_stl_facets[cell],stl_facet)
+  ncells = num_cells(grid)
+  assemble_threaded_sparse_map(thread_to_cells,thread_to_stl_facets,ncells)
+end
+
+function compute_cell_to_facets(a::UnstructuredGrid,b)
+  tmp = cartesian_bounding_model(a)
+  tmp_to_a = compute_cell_to_facets(tmp,a)
+  tmp_to_b = compute_cell_to_facets(tmp,b)
+  a_to_tmp = inverse_index_map(tmp_to_a,num_cells(a))
+  a_to_b = compose_index_map(a_to_tmp,tmp_to_b)
+  a_to_b = filter_intersections(a,b,a_to_b)
+  a_to_b
+end
+
+function filter_intersections(a::Grid,b,a_to_b)
+  CELL_EXPANSION_FACTOR = 1e-3
+  n = Threads.nthreads()
+  t_as = map(_->Int32[],1:n)
+  t_bs = map(_->Int32[],1:n)
+  t_c = map(_->(get_cell_cache(a),get_cell_cache(b),array_cache(a_to_b)),1:n)
+  Threads.@threads for ai in 1:num_cells(a)
+    i = Threads.threadid()
+    ca,cb,c = t_c[i]
+    as = t_as[i]
+    bs = t_bs[i]
+    cell_a = get_cell!(ca,a,ai)
+    pmin,pmax = get_bounding_box(cell_a)
+    Δ = norm(pmin-pmax) * CELL_EXPANSION_FACTOR
+    cell_a = expand_face(cell_a,Δ)
+    for bi in getindex!(c,a_to_b,ai)
+      cell_b = get_cell!(cb,b,bi)
+      if has_intersection(cell_a,cell_b)
+        push!(as,ai)
+        push!(bs,bi)
+      end
+    end
+  end
+  assemble_threaded_sparse_map(t_as,t_bs,num_cells(a))
+end
+
+function cartesian_bounding_model(grid::Grid)
+  h = max_length(grid)
+  pmin,pmax = get_bounding_box(grid)
+  s = Tuple(pmax-pmin)
+  partition = Int.( ceil.( s ./ h  ))
+  CartesianDiscreteModel(pmin,pmax,partition)
+end
+
+function inverse_index_map(
+  a_to_b::AbstractVector{<:AbstractVector},
+  nb=maximum(maximum(a_to_b)))
+
+  na = length(a_to_b)
+  as = Int32[]
+  bs = Int32[]
+  cache = array_cache(a_to_b)
+  for a in 1:na
+    for b in getindex!(cache,a_to_b,a)
+      push!(as,a)
+      push!(bs,b)
+    end
+  end
+  assemble_sparse_map(bs,as,nb)
+end
+
+function compose_index_map(
+  a_to_b::AbstractVector{<:AbstractVector},
+  b_to_c::AbstractVector{<:AbstractVector})
+
+  na = length(a_to_b)
+  as = Int32[]
+  cs = Int32[]
+  cache_ab = array_cache(a_to_b)
+  cache_bc = array_cache(b_to_c)
+  for a in 1:na
+    for b in getindex!(cache_ab,a_to_b,a)
+      for c in getindex!(cache_bc,b_to_c,b)
+        push!(as,a)
+        push!(cs,c)
+      end
+    end
+  end
+  assemble_sparse_map(as,cs,na)
+end
+
+function assemble_threaded_sparse_map(tas,tbs,na=maximum(maximum,tas))
+  a_to_b = map(_->Int32[],1:na)
+  for (as,bs) in zip(tas,tbs)
+    assemble_sparse_map!(a_to_b,as,bs)
+  end
+  a_to_b
+end
+
+function assemble_sparse_map(as,bs,na=maximum(as))
+  a_to_b = map(_->Int32[],1:na)
+  assemble_sparse_map!(a_to_b,as,bs)
+  a_to_b
+end
+
+function assemble_sparse_map!(a_to_b,as,bs)
+  for (a,b) in zip(as,bs)
+    if b ∉ a_to_b[a]
+      push!(a_to_b[a],b)
     end
   end
-  cell_to_stl_facets
+  a_to_b
 end
 
 function get_cells_around(desc::CartesianDescriptor{D},pmin::Point,pmax::Point) where D
   cmin,_ = get_cell_bounds(desc,pmin)
   _,cmax = get_cell_bounds(desc,pmax)
   cmin = CartesianIndices(desc.partition)[cmin]
   cmax = CartesianIndices(desc.partition)[cmax]
-  ranges = ntuple( i -> cmin.I[i] : cmax.I[i], Val{D}() )
+  ranges = ntuple( i -> cmin.I[i]::Int : cmax.I[i]::Int, Val{D}() )
   CartesianIndices( ranges )
 end