From 79d34031339d135ad814f4cfee04bdf06d0fb9a2 Mon Sep 17 00:00:00 2001 From: amartin Date: Wed, 6 May 2020 16:42:07 +1000 Subject: [PATCH 1/5] Added a new, more general inner constructor to CartesianDiscreteModel --- src/Geometry/CartesianDiscreteModels.jl | 180 +++++++++++++++++- .../CartesianDiscreteModelsTests.jl | 3 + 2 files changed, 176 insertions(+), 7 deletions(-) diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl index c43565b2b..a9d23cfe5 100644 --- a/src/Geometry/CartesianDiscreteModels.jl +++ b/src/Geometry/CartesianDiscreteModels.jl @@ -22,6 +22,38 @@ struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D} _fill_cartesian_face_labeling!(labels,topo) new{D,T,F}(grid,topo,labels) end + + @doc """ + CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}, + cmin::CartesianIndex, + cmax::CartesianIndex) + + Builds a CartesianDiscreteModel object which represents a subgrid of + a (larger) grid represented by desc. This subgrid is described by its + D-dimensional minimum (cmin) and maximum (cmax) CartesianIndex + identifiers. + + Inner constructor + """ + function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}, + cmin::CartesianIndex, + cmax::CartesianIndex) where {D,T,F} + + suborigin = Tuple(desc.origin) .+ + (Tuple(cmin) .- 1) .* Tuple(desc.sizes) + subpartition = Tuple(cmax) .- Tuple(cmin) .+ 1 + subsizes = Tuple(desc.sizes) + subdesc = + CartesianDescriptor(suborigin, subsizes, subpartition, desc.map) + + grid = CartesianGrid(subdesc) + _grid = UnstructuredGrid(grid) + topo = UnstructuredGridTopology(_grid) + nfaces = [num_faces(topo, d) for d = 0:num_cell_dims(topo)] + labels = FaceLabeling(nfaces) + _fill_subgrid_cartesian_face_labeling!(labels,topo,subdesc,desc,cmin) + new{D,T,F}(grid, topo, labels) + end end """ @@ -34,6 +66,7 @@ function CartesianDiscreteModel(args...) CartesianDiscreteModel(desc) end + """ get_cartesian_descriptor(model::CartesianDiscreteModel) """ @@ -91,9 +124,14 @@ function _fill_cartesian_entities!(labels,topo) offset, interior_id,boundary_id,d,D) end - for d in 0:(D-2) - for j in (d+1):(D-1) - dface_to_jfaces = get_faces(topo,d,j) + _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id) + fill!(d_to_dface_to_entity[end],interior_id) +end + +function _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id) + for d = 0:(D-2) + for j = (d+1):(D-1) + dface_to_jfaces = get_faces(topo, d, j) dface_to_geolabel = d_to_dface_to_entity[d+1] jface_to_geolabel = d_to_dface_to_entity[j+1] _fix_dface_geolabels!( @@ -101,10 +139,11 @@ function _fill_cartesian_entities!(labels,topo) jface_to_geolabel, dface_to_jfaces.data, dface_to_jfaces.ptrs, - interior_id,boundary_id) + interior_id, + boundary_id, + ) end end - fill!(d_to_dface_to_entity[end],interior_id) end function _add_cartesian_tags!(labels,topo) @@ -151,7 +190,6 @@ function _generate_pre_geolabel_kernel!( cell_to_faces_data, cell_to_faces_ptrs, offset,boundary_id,max_ncells_around) - nfaces = length(face_to_geolabel) for face in 1:nfaces a = face_to_cells_ptrs[face]-1 @@ -172,7 +210,6 @@ function _generate_pre_geolabel_kernel!( face_to_geolabel[face] = boundary_id end end - end function _fix_dface_geolabels!( @@ -199,3 +236,132 @@ function _fix_dface_geolabels!( end end end + +function _fill_subgrid_cartesian_face_labeling!(labels,topo,subdesc,desc,cmin) + _fill_subgrid_cartesian_entities!(labels,topo,subdesc,desc,cmin) + _add_cartesian_tags!(labels,topo) +end + +function _fill_subgrid_cartesian_entities!(labels, topo, subdesc, desc, cmin) + D = num_cell_dims(topo) + d_to_dface_to_entity = labels.d_to_dface_to_entity + gcis = CartesianIndices(Tuple(desc.partition)) + subcis = CartesianIndices(Tuple(subdesc.partition)) + + polytope = first(get_polytopes(topo)) + face_labeling = labels + offsets = get_offsets(polytope) + interior_id = num_faces(polytope) + boundary_id = -1 + + minus_one_ci = CartesianIndex(tfill(-1, Val{D}())) + + polytope_d_face_to_jfaces = Matrix{Vector{Vector{Int}}}(undef, (D, D)) + for d = 0:(D-1) + for j = d+1:D-1 + polytope_d_face_to_jfaces[d+1, j+1] = get_faces(polytope, d, j) + end + end + + face_deltas = _find_ncube_face_neighbor_deltas(polytope) + for d = 0:(D-1) + face_to_cells = get_faces(topo, d, D) + cell_to_faces = get_faces(topo, D, d) + face_to_geolabel = face_labeling.d_to_dface_to_entity[d+1] + nfaces = length(face_to_geolabel) + for face_gid = 1:nfaces + cell_gid = face_to_cells.data[face_to_cells.ptrs[face_gid]] + a = cell_to_faces.ptrs[cell_gid] + b = cell_to_faces.ptrs[cell_gid+1] - 1 + gci = (subcis[cell_gid] + minus_one_ci) + cmin + + face_lid = -1 + for j = a:b + if (cell_to_faces.data[j] == face_gid) + face_lid = j - a + 1 + break + end + end + @assert face_lid != -1 + + face_lid += offsets[d+1] + # Check whether cell neighbour across face face_lid belongs to the + # global grid. If yes, the current face is actually at the interior + is_assigned_face_delta = isassigned(face_deltas, face_lid) + if (is_assigned_face_delta && (gci + face_deltas[face_lid]) in gcis) + face_to_geolabel[face_gid] = interior_id + else + cell_found = false + for j = d+1:D-1 + dface_to_jfaces = + polytope_d_face_to_jfaces[d+1, j+1][face_lid-offsets[d+1]] + cell_found = _is_there_interior_cell_across_higher_dim_faces( + dface_to_jfaces, + offsets[j+1], + gcis, + gci, + face_deltas, + ) + cell_found && break + end + if (cell_found) + # The current face is at least in two subgrid cells + face_to_geolabel[face_gid] = boundary_id + else + # The current face is only in one subgrid cell + face_to_geolabel[face_gid] = face_lid + end + end + end + end + _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id) + fill!(d_to_dface_to_entity[end], interior_id) +end + +function _is_there_interior_cell_across_higher_dim_faces( + dface_to_jfaces, + offset_j, + gcis, + gci, + face_deltas, +) + for k in dface_to_jfaces + jface_lid = k + offset_j + if (isassigned(face_deltas, jface_lid)) + if ((gci + face_deltas[jface_lid]) in gcis) + return true + end + end + end + return false +end + +""" + _find_ncube_face_neighbor_deltas(p::Polytope{D}) -> Vector{CartesianIndex} + + Given an n-cube type Polytope{D}, returns V=Vector{CartesianIndex} with as many + entries as the number of faces in the boundary of the Polytope. For an entry face_lid + in this vector, V[face_lid] returns what has to be added to the CartesianIndex of a + cell in order to obtain the CartesianIndex of the cell neighbour of K across the face F + with local ID face_lid. +""" +function _find_ncube_face_neighbor_deltas(p::Polytope{D}) where {D} + nfaces = num_faces(p) + delta_faces = Vector{CartesianIndex}(undef, nfaces - 1) + for face_lid = 1:nfaces-1 + delta_faces[face_lid] = _find_ncube_face_neighbor_delta(p, face_lid) + end + delta_faces +end + +function _find_ncube_face_neighbor_delta(p::Polytope{D}, face_lid) where {D} + @assert is_n_cube(p) + result = fill(0, D) + face = p.dface.nfaces[face_lid] + for d = 1:D + if (face.extrusion[d] == 0) + result[d] = (face.anchor[d] == 0) ? -1 : 1 + end + end + return CartesianIndex(Tuple(result)) +end diff --git a/test/GeometryTests/CartesianDiscreteModelsTests.jl b/test/GeometryTests/CartesianDiscreteModelsTests.jl index 66c048b4b..8728a48bf 100644 --- a/test/GeometryTests/CartesianDiscreteModelsTests.jl +++ b/test/GeometryTests/CartesianDiscreteModelsTests.jl @@ -43,6 +43,9 @@ test_discrete_model(model) model2 = from_dict(DiscreteModel,to_dict(model)) test_discrete_model(model2) +model3 = CartesianDiscreteModel(desc,CartesianIndex(2,2,2),CartesianIndex(3,3,3)) +test_discrete_model(model3) + #using Gridap.Visualization #writevtk(model2,"model2") From e4eda3cf9c768438b75cd7ebca348e5076a3ac4d Mon Sep 17 00:00:00 2001 From: amartin Date: Thu, 7 May 2020 10:33:37 +1000 Subject: [PATCH 2/5] Restricted _find_ncube_face_neighbor_deltas helper function arguments to type ExtrusionPolytope{D} --- src/Geometry/CartesianDiscreteModels.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl index a9d23cfe5..017c4f843 100644 --- a/src/Geometry/CartesianDiscreteModels.jl +++ b/src/Geometry/CartesianDiscreteModels.jl @@ -337,15 +337,15 @@ function _is_there_interior_cell_across_higher_dim_faces( end """ - _find_ncube_face_neighbor_deltas(p::Polytope{D}) -> Vector{CartesianIndex} + _find_ncube_face_neighbor_deltas(p::ExtrusionPolytope{D}) -> Vector{CartesianIndex} - Given an n-cube type Polytope{D}, returns V=Vector{CartesianIndex} with as many + Given an n-cube type ExtrusionPolytope{D}, returns V=Vector{CartesianIndex} with as many entries as the number of faces in the boundary of the Polytope. For an entry face_lid in this vector, V[face_lid] returns what has to be added to the CartesianIndex of a cell in order to obtain the CartesianIndex of the cell neighbour of K across the face F with local ID face_lid. """ -function _find_ncube_face_neighbor_deltas(p::Polytope{D}) where {D} +function _find_ncube_face_neighbor_deltas(p::ExtrusionPolytope{D}) where {D} nfaces = num_faces(p) delta_faces = Vector{CartesianIndex}(undef, nfaces - 1) for face_lid = 1:nfaces-1 @@ -354,7 +354,7 @@ function _find_ncube_face_neighbor_deltas(p::Polytope{D}) where {D} delta_faces end -function _find_ncube_face_neighbor_delta(p::Polytope{D}, face_lid) where {D} +function _find_ncube_face_neighbor_delta(p::ExtrusionPolytope{D}, face_lid) where {D} @assert is_n_cube(p) result = fill(0, D) face = p.dface.nfaces[face_lid] From a9e3d463ead74b417a8b2e384ee9fc3b3065bf1e Mon Sep 17 00:00:00 2001 From: amartin Date: Thu, 7 May 2020 20:07:29 +1000 Subject: [PATCH 3/5] Mitigated the effect of type instability in _fill_subgrid_cartesian_entities! introducing a function barrier (https://docs.julialang.org/en/v1/manual/performance-tips/index.html#kernel-functions-1) --- src/Geometry/CartesianDiscreteModels.jl | 120 +++++++++++++++--------- 1 file changed, 77 insertions(+), 43 deletions(-) diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl index 017c4f843..39e3e2bae 100644 --- a/src/Geometry/CartesianDiscreteModels.jl +++ b/src/Geometry/CartesianDiscreteModels.jl @@ -268,56 +268,90 @@ function _fill_subgrid_cartesian_entities!(labels, topo, subdesc, desc, cmin) face_to_cells = get_faces(topo, d, D) cell_to_faces = get_faces(topo, D, d) face_to_geolabel = face_labeling.d_to_dface_to_entity[d+1] - nfaces = length(face_to_geolabel) - for face_gid = 1:nfaces - cell_gid = face_to_cells.data[face_to_cells.ptrs[face_gid]] - a = cell_to_faces.ptrs[cell_gid] - b = cell_to_faces.ptrs[cell_gid+1] - 1 - gci = (subcis[cell_gid] + minus_one_ci) + cmin - - face_lid = -1 - for j = a:b - if (cell_to_faces.data[j] == face_gid) - face_lid = j - a + 1 - break - end + _generate_subgrid_pregeo_label!( + face_to_geolabel, + d, + D, + offsets, + polytope_d_face_to_jfaces, + minus_one_ci, + face_to_cells, + cell_to_faces, + subcis, + gcis, + face_deltas, + cmin, + interior_id, + boundary_id, + ) + end + _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id) + fill!(d_to_dface_to_entity[end], interior_id) +end + +function _generate_subgrid_pregeo_label!( + face_to_geolabel, + d, + D, + offsets, + polytope_d_face_to_jfaces, + minus_one_ci, + face_to_cells, + cell_to_faces, + subcis, + gcis, + face_deltas, + cmin, + interior_id, + boundary_id, +) + for face_gid = 1:length(face_to_geolabel) + cell_gid = face_to_cells.data[face_to_cells.ptrs[face_gid]] + a = cell_to_faces.ptrs[cell_gid] + b = cell_to_faces.ptrs[cell_gid+1] - 1 + gci = (subcis[cell_gid] + minus_one_ci) + cmin + + face_lid = -1 + for j = a:b + if (cell_to_faces.data[j] == face_gid) + face_lid = j - a + 1 + break + end + end + @assert face_lid != -1 + + face_lid += offsets[d+1] + # Check whether cell neighbour across face face_lid belongs to the + # global grid. If yes, the current face is actually at the interior + is_assigned_face_delta = isassigned(face_deltas, face_lid) + if (is_assigned_face_delta && (gci + face_deltas[face_lid]) in gcis) + face_to_geolabel[face_gid] = interior_id + else + cell_found = false + for j = d+1:D-1 + dface_to_jfaces = + polytope_d_face_to_jfaces[d+1, j+1][face_lid-offsets[d+1]] + cell_found = _is_there_interior_cell_across_higher_dim_faces( + dface_to_jfaces, + offsets[j+1], + gcis, + gci, + face_deltas, + ) + cell_found && break end - @assert face_lid != -1 - - face_lid += offsets[d+1] - # Check whether cell neighbour across face face_lid belongs to the - # global grid. If yes, the current face is actually at the interior - is_assigned_face_delta = isassigned(face_deltas, face_lid) - if (is_assigned_face_delta && (gci + face_deltas[face_lid]) in gcis) - face_to_geolabel[face_gid] = interior_id + if (cell_found) + # The current face is at least in two subgrid cells + face_to_geolabel[face_gid] = boundary_id else - cell_found = false - for j = d+1:D-1 - dface_to_jfaces = - polytope_d_face_to_jfaces[d+1, j+1][face_lid-offsets[d+1]] - cell_found = _is_there_interior_cell_across_higher_dim_faces( - dface_to_jfaces, - offsets[j+1], - gcis, - gci, - face_deltas, - ) - cell_found && break - end - if (cell_found) - # The current face is at least in two subgrid cells - face_to_geolabel[face_gid] = boundary_id - else - # The current face is only in one subgrid cell - face_to_geolabel[face_gid] = face_lid - end + # The current face is only in one subgrid cell + face_to_geolabel[face_gid] = face_lid end end end - _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id) - fill!(d_to_dface_to_entity[end], interior_id) end + function _is_there_interior_cell_across_higher_dim_faces( dface_to_jfaces, offset_j, From 7fb3ae88f8b4ba566d55603a7cb870cdb03cf48b Mon Sep 17 00:00:00 2001 From: amartin Date: Thu, 7 May 2020 20:49:42 +1000 Subject: [PATCH 4/5] Minor in _generate_subgrid_pregeo_label! --- src/Geometry/CartesianDiscreteModels.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl index 39e3e2bae..ee0ffe166 100644 --- a/src/Geometry/CartesianDiscreteModels.jl +++ b/src/Geometry/CartesianDiscreteModels.jl @@ -323,8 +323,7 @@ function _generate_subgrid_pregeo_label!( face_lid += offsets[d+1] # Check whether cell neighbour across face face_lid belongs to the # global grid. If yes, the current face is actually at the interior - is_assigned_face_delta = isassigned(face_deltas, face_lid) - if (is_assigned_face_delta && (gci + face_deltas[face_lid]) in gcis) + if ((gci + face_deltas[face_lid]) in gcis) face_to_geolabel[face_gid] = interior_id else cell_found = false From aa9e320e46dd57fe743196febd4cee571cf50eba Mon Sep 17 00:00:00 2001 From: amartin Date: Thu, 7 May 2020 22:06:54 +1000 Subject: [PATCH 5/5] Type instability fix in _find_ncube_face_neighbor_deltas --- src/Geometry/CartesianDiscreteModels.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl index ee0ffe166..b66301535 100644 --- a/src/Geometry/CartesianDiscreteModels.jl +++ b/src/Geometry/CartesianDiscreteModels.jl @@ -380,7 +380,7 @@ end """ function _find_ncube_face_neighbor_deltas(p::ExtrusionPolytope{D}) where {D} nfaces = num_faces(p) - delta_faces = Vector{CartesianIndex}(undef, nfaces - 1) + delta_faces = Vector{CartesianIndex{D}}(undef, nfaces - 1) for face_lid = 1:nfaces-1 delta_faces[face_lid] = _find_ncube_face_neighbor_delta(p, face_lid) end