diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl index 935e3b9..0e1b602 100644 --- a/src/Distributed/DistributedDiscreteGeometries.jl +++ b/src/Distributed/DistributedDiscreteGeometries.jl @@ -4,41 +4,28 @@ end local_views(a::DistributedDiscreteGeometry) = a.geometries +# TODO: Is this really necessary? function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) where {Dc,Dp} @assert DomainStyle(φh) == ReferenceDomain() gids = get_cell_gids(model) values = map(local_views(φh),local_views(model),local_views(gids)) do φh, model, gids - # Maps from the no-ghost model to the original model own_model = remove_ghost_cells(model,gids) - own_to_local_node = get_face_to_parent_face(own_model,0) - local_to_own_node = find_inverse_index_map(own_to_local_node,num_nodes(model)) - own_to_local_cell = get_face_to_parent_face(own_model,Dc) - - # Cell-to-node map for the original model - # topo = get_grid_topology(model) - # c2n_map = get_faces(topo,Dc,0) - c2n_map = collect1d(get_cell_node_ids(model)) - - # Cell-wise node coordinates (in ReferenceDomain coordinates) - cell_reffe = get_cell_reffe(model) - cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) - - φh_data = CellData.get_data(φh) - T = return_type(testitem(CellData.get_data(φh)),testitem(testitem(cell_node_coords))) + own_cells = get_face_to_parent_face(own_model,Dc) + + trian = get_triangulation(φh) + cell_points = get_cell_points(trian) + cell_ids = get_cell_node_ids(own_model) + cell_values = φh(cell_points) + + T = eltype(testitem(cell_values)) values = Vector{T}(undef,num_nodes(own_model)) - touched = fill(false,num_nodes(model)) - cell_node_coords_cache = array_cache(cell_node_coords) - for cell in own_to_local_cell # For each owned cell - field = φh_data[cell] - node_coords = getindex!(cell_node_coords_cache,cell_node_coords,cell) - for (iN,node) in enumerate(c2n_map[cell]) # Go over local nodes - own_node = local_to_own_node[node] - if (own_node != 0) && !touched[node] # Compute value if suitable - values[own_node] = field(node_coords[iN]) - touched[node] = true - end - end + cell_ids_cache = array_cache(cell_ids) + cell_values_cache = array_cache(cell_values) + for (ocell,cell) in enumerate(own_cells) + ids = getindex!(cell_ids_cache,cell_ids,ocell) + vals = getindex!(cell_values_cache,cell_values,cell) + values[ids] .= vals end return values end @@ -56,7 +43,7 @@ function DiscreteGeometry(φh::CellField,model::DistributedDiscreteModel;name::S DistributedDiscreteGeometry(geometries) end -function get_geometry(a::AbstractArray{<:DiscreteGeometry}) +function distributed_geometry(a::AbstractArray{<:DiscreteGeometry}) DistributedDiscreteGeometry(a) end @@ -104,8 +91,8 @@ function distributed_embedded_triangulation( T, cutgeo::DistributedEmbeddedDiscretization, cutinorout, - geom::DistributedDiscreteGeometry) - + geom::DistributedDiscreteGeometry +) trians = map(local_views(cutgeo),local_views(geom)) do lcutgeo,lgeom T(lcutgeo,cutinorout,lgeom) end @@ -117,8 +104,8 @@ function distributed_aggregate( strategy::AggregateCutCellsByThreshold, cut::DistributedEmbeddedDiscretization, geo::DistributedDiscreteGeometry, - in_or_out=IN) - + in_or_out = IN +) bgmodel = get_background_model(cut) facet_to_inoutcut = compute_bgfacet_to_inoutcut(bgmodel,geo) _distributed_aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut) @@ -133,8 +120,8 @@ end function compute_bgfacet_to_inoutcut( cutter::Cutter, bgmodel::DistributedDiscreteModel, - geo::DistributedDiscreteGeometry) - + geo::DistributedDiscreteGeometry +) gids = get_cell_gids(bgmodel) bgf_to_ioc = map(local_views(bgmodel),local_views(gids),local_views(geo)) do model,gids,geo ownmodel = remove_ghost_cells(model,gids) diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl index 9ae968c..8bb2d6e 100644 --- a/src/Distributed/DistributedDiscretizations.jl +++ b/src/Distributed/DistributedDiscretizations.jl @@ -3,11 +3,12 @@ struct DistributedEmbeddedDiscretization{A,B} <: GridapType discretizations::A model::B function DistributedEmbeddedDiscretization( - d::AbstractArray{<:AbstractEmbeddedDiscretization}, - model::DistributedDiscreteModel) - A = typeof(d) + discretizations::AbstractArray{<:AbstractEmbeddedDiscretization}, + model::DistributedDiscreteModel + ) + A = typeof(discretizations) B = typeof(model) - new{A,B}(d,model) + new{A,B}(discretizations,model) end end @@ -16,12 +17,13 @@ local_views(a::DistributedEmbeddedDiscretization) = a.discretizations get_background_model(a::DistributedEmbeddedDiscretization) = a.model function get_geometry(a::DistributedEmbeddedDiscretization) - loc_geometries = map(get_geometry,local_views(a)) - get_geometry(loc_geometries) + geometries = map(get_geometry,local_views(a)) + distributed_geometry(geometries) end -function get_geometry(a::AbstractArray{<:CSG.Geometry}) - PartitionedArrays.getany(a) +# Neeed for dispatching between analytical geometries and distributed geometries +function distributed_geometry(geometries::AbstractArray{<:CSG.Geometry}) + PartitionedArrays.getany(geometries) end function cut(bgmodel::DistributedDiscreteModel,args...) @@ -127,8 +129,8 @@ end function compute_bgfacet_to_inoutcut( cutter::Cutter, bgmodel::DistributedDiscreteModel, - geo) - + geo +) gids = get_cell_gids(bgmodel) bgf_to_ioc = map(local_views(bgmodel),local_views(gids)) do model,gids ownmodel = remove_ghost_cells(model,gids) @@ -243,12 +245,11 @@ function change_bgmodel( DistributedEmbeddedDiscretization(cuts,model) end - function _change_bgmodels( cutgeo::DistributedEmbeddedDiscretization, model::DistributedDiscreteModel, - cell_to_newcell) - + cell_to_newcell +) map(local_views(cutgeo),local_views(model),cell_to_newcell) do c,m,c_to_nc change_bgmodel(c,m,c_to_nc) end @@ -256,8 +257,8 @@ end function _change_bgmodels( cutgeo::DistributedEmbeddedDiscretization, - model::DistributedDiscreteModel) - + model::DistributedDiscreteModel +) map(local_views(cutgeo),local_views(model)) do c,m change_bgmodel(c,m) end @@ -266,8 +267,8 @@ end function change_bgmodel( cut::EmbeddedDiscretization, newmodel::DiscreteModel, - cell_to_newcell=1:num_cells(get_background_model(cut))) - + cell_to_newcell=1:num_cells(get_background_model(cut)) +) ls_to_bgc_to_ioc = map(cut.ls_to_bgcell_to_inoutcut) do bgc_to_ioc new_bgc_to_ioc = Vector{Int8}(undef,num_cells(newmodel)) new_bgc_to_ioc[cell_to_newcell] = bgc_to_ioc @@ -289,10 +290,9 @@ end function change_bgmodel( cut::EmbeddedFacetDiscretization, newmodel::DiscreteModel, - facet_to_newfacet=1:num_facets(get_background_model(cut))) - + facet_to_newfacet=1:num_facets(get_background_model(cut)) +) nfacets = num_facets(newmodel) - ls_to_bgf_to_ioc = map(cut.ls_to_facet_to_inoutcut) do bgf_to_ioc new_bgf_to_ioc = Vector{Int8}(undef,nfacets) new_bgf_to_ioc[facet_to_newfacet] = bgf_to_ioc @@ -305,7 +305,8 @@ function change_bgmodel( subfacets, cut.ls_to_subfacet_to_inout, cut.oid_to_ls, - cut.geo) + cut.geo + ) end function change_bgmodel(cells::SubCellData,cell_to_newcell) diff --git a/src/LevelSetCutters/DiscreteGeometries.jl b/src/LevelSetCutters/DiscreteGeometries.jl index 23d8cbe..472fa27 100644 --- a/src/LevelSetCutters/DiscreteGeometries.jl +++ b/src/LevelSetCutters/DiscreteGeometries.jl @@ -34,7 +34,6 @@ function discretize(a::AnalyticalGeometry,point_to_coords::AbstractArray{<:Point end function discretize(a::AnalyticalGeometry,point_to_coords::Vector{<:Point}) - tree = get_tree(a) j_to_fun, oid_to_j = _find_unique_leaves(tree) j_to_ls = [ fun.(point_to_coords) for fun in j_to_fun ] @@ -48,13 +47,10 @@ function discretize(a::AnalyticalGeometry,point_to_coords::Vector{<:Point}) end newtree = replace_data(identity,conversion,tree) - DiscreteGeometry(newtree,point_to_coords) - end function _find_unique_leaves(tree) - i_to_fun = map(n->first(n.data),collect(Leaves(tree))) i_to_oid = map(objectid,i_to_fun) j_to_oid = unique(i_to_oid) @@ -65,31 +61,6 @@ function _find_unique_leaves(tree) j_to_fun, oid_to_j end -function _get_value_at_coords(φh::CellField,model::DiscreteModel{Dc,Dp}) where {Dc,Dp} - @assert DomainStyle(φh) == ReferenceDomain() - # Cell-to-node map for the original model - c2n_map = collect1d(get_cell_node_ids(model)) - - # Cell-wise node coordinates (in ReferenceDomain coordinates) - cell_reffe = get_cell_reffe(model) - cell_node_coords = lazy_map(get_node_coordinates,cell_reffe) - - # Get cell data - φh_data = CellData.get_data(φh) - T = return_type(testitem(CellData.get_data(φh)),testitem(testitem(cell_node_coords))) - values = Vector{T}(undef,num_nodes(model)) - cell_node_coords_cache = array_cache(cell_node_coords) - # Loop over cells - for cell in eachindex(c2n_map) - field = φh_data[cell] - node_coords = getindex!(cell_node_coords_cache,cell_node_coords,cell) - for (iN,node) in enumerate(c2n_map[cell]) - values[node] = field(node_coords[iN]) - end - end - return values -end - function DiscreteGeometry( point_to_value::AbstractVector,point_to_coords::AbstractVector;name::String="") data = (point_to_value,name,nothing) @@ -97,9 +68,14 @@ function DiscreteGeometry( DiscreteGeometry(tree,point_to_coords) end +# TODO: This assumes that the level set φh is 1st order, i.e that there is a 1-to-1 correspondence +# between nodes in the mesh and dofs in φh. +# Even if we allowed higher order, the cuts are always linear. Not only it would be a waste +# of time to use higher order, but cuts could actually be wrong. +# This might be developped in the future. function DiscreteGeometry( φh::CellField,model::DiscreteModel;name::String="") - point_to_value = _get_value_at_coords(φh,model) + point_to_value = get_free_dof_values(φh) point_to_coords = collect1d(get_node_coordinates(model)) DiscreteGeometry(point_to_value,point_to_coords;name) -end \ No newline at end of file +end