From 5212b935eeaccec3c1260f6c4c35c04d6847a4e5 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 15 Aug 2024 15:34:02 -0400 Subject: [PATCH 01/34] WIP: form patches --- src/CMakeLists.txt | 1 + src/formPatches.cpp | 45 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) create mode 100644 src/formPatches.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b9e15ab0..0582fc5f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -391,6 +391,7 @@ if(BUILD_TESTING) osh_add_exe(arrayops_test) osh_add_exe(sort_test) osh_add_exe(unit_math) + osh_add_exe(formPatches) if (Omega_h_USE_Kokkos) osh_add_exe(bbox_reduce_test) osh_add_exe(initKokkosAndLib) diff --git a/src/formPatches.cpp b/src/formPatches.cpp new file mode 100644 index 00000000..6a0a26c0 --- /dev/null +++ b/src/formPatches.cpp @@ -0,0 +1,45 @@ +#include +#include +#include + +#include + +using namespace Omega_h; + +/** + * \brief expand the patches + * \param m (in) mesh of simplices + * \param patches (in) graph of key entities to elements + * \param bridgeDim (in) the entity dimension used for expansion via second + * order element-to-element adjacencies + * \return an expanded graph from key entities to elements +*/ +Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { + OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); + return Graph(); +} + +/** + * \brief form a patch of at least minPatchSize elements surrounding each entity + * of dimension keyDim + * \param m (in) mesh of simplices + * \param keyDim (in) the dimension of mesh entities that the patches are + * created around + * \param minPatchSize (in) the minimum number of elements in each patch + * \return a graph whose source nodes are the entities of keyDim dimension, and + * edges are connecting to elements in the patch + */ +Graph formPatches(Mesh& m, LO keyDim, Int minPatchSize) { + OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); + OMEGA_H_CHECK(minPatchSize > 0); + return Graph(); +} + +int main(int argc, char** argv) { + auto lib = Library(&argc, &argv); + OMEGA_H_CHECK(argc == 3); + Mesh mesh(&lib); + binary::read(argv[1], lib.world(), &mesh); + auto patches = formPatches(mesh, VERT, 3); + vtk::write_parallel(argv[2], &mesh, dim); +} From 46508677ca8d5f07bd95fffab71016900ea4cfc6 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 15 Aug 2024 16:22:31 -0400 Subject: [PATCH 02/34] compiles --- src/formPatches.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 6a0a26c0..8dc10fad 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -15,7 +15,7 @@ using namespace Omega_h; * \return an expanded graph from key entities to elements */ Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { - OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); + OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); return Graph(); } @@ -41,5 +41,5 @@ int main(int argc, char** argv) { Mesh mesh(&lib); binary::read(argv[1], lib.world(), &mesh); auto patches = formPatches(mesh, VERT, 3); - vtk::write_parallel(argv[2], &mesh, dim); + vtk::write_parallel(argv[2], &mesh, mesh.dim()); } From 5852895d4a6aaed68a03285953b361a074b724b4 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 16 Aug 2024 11:22:39 -0400 Subject: [PATCH 03/34] add expansion loop doesn't support more than one expansion round for now --- src/formPatches.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 8dc10fad..89cb9e55 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -6,6 +6,14 @@ using namespace Omega_h; +bool patchSufficient(Graph patches, Int minPatchSize) { + //find min degree for each patch + //if( minDegree < minPatchSize) + // return false; + //else + return true; +} + /** * \brief expand the patches * \param m (in) mesh of simplices @@ -32,6 +40,14 @@ Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { Graph formPatches(Mesh& m, LO keyDim, Int minPatchSize) { OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); OMEGA_H_CHECK(minPatchSize > 0); + auto patches = Graph(); + for(Int bridgeDim = m.dim()-1; bridgeDim >= 0; bridgeDim--) { + auto bridgePatches = expandPatches(m, patches, bridgeDim); + if( patchSufficient(bridgePatches, minPatchSize) ) { + return bridgePatches; + } + } + assert(false); return Graph(); } From ea133cd0a60b45941f95e29a37d54388c14fca37 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 16 Aug 2024 22:07:16 -0400 Subject: [PATCH 04/34] wip --- src/formPatches.cpp | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 89cb9e55..b0c98f2a 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -1,11 +1,17 @@ #include #include #include +#include #include using namespace Omega_h; +Graph getElmToElm2ndOrderAdj(Mesh& m, Int bridgeDim) { + OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); + return Graph(); +} + bool patchSufficient(Graph patches, Int minPatchSize) { //find min degree for each patch //if( minDegree < minPatchSize) @@ -14,6 +20,28 @@ bool patchSufficient(Graph patches, Int minPatchSize) { return true; } +Graph adj_segment_sort(Graph& g) { + using ExecSpace = Kokkos::DefaultExecutionSpace; + using TeamPol = Kokkos::TeamPolicy; + using TeamMem = typename TeamPol::member_type; + auto offsets = g.a2ab; + auto elms = g.ab2b.view(); + Kokkos::View elms_v("elms", elms.size()); + Kokkos::deep_copy(elms_v, elms); + auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { + //Sort a row of A using the whole team. + auto i = t.league_rank(); + auto patch = Kokkos::subview(elms_v, Kokkos::make_pair(offsets[i], offsets[i+1])); + Kokkos::Experimental::sort_team(t, patch); + }; + Kokkos::parallel_for(TeamPol(offsets.size()-1, Kokkos::AUTO()), segment_sort); + return Graph(offsets,elms_w);//need a LOs instead of elms_w... +} + +Graph remove_duplicate_edges(Graph g) { + return Graph(); +} + /** * \brief expand the patches * \param m (in) mesh of simplices @@ -24,7 +52,12 @@ bool patchSufficient(Graph patches, Int minPatchSize) { */ Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); - return Graph(); + auto patch_elms = patches.ab2b; + auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); + auto expanded = unmap_graph(patch_elms, adjElms); + adj_segment_sort(expanded); + auto dedup = remove_duplicate_edges(expanded); + return dedup; } /** From 4133ccd6f1bb767e971bb9845f321d63f32c76b8 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 06:27:01 -0400 Subject: [PATCH 05/34] segment_sort compiles --- src/formPatches.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index b0c98f2a..d5277ca1 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -25,17 +25,17 @@ Graph adj_segment_sort(Graph& g) { using TeamPol = Kokkos::TeamPolicy; using TeamMem = typename TeamPol::member_type; auto offsets = g.a2ab; - auto elms = g.ab2b.view(); - Kokkos::View elms_v("elms", elms.size()); - Kokkos::deep_copy(elms_v, elms); + auto elms_r = g.ab2b.view(); //read only + Kokkos::View elms("elms", elms.size()); + Kokkos::deep_copy(elms, elms_r); auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { //Sort a row of A using the whole team. auto i = t.league_rank(); - auto patch = Kokkos::subview(elms_v, Kokkos::make_pair(offsets[i], offsets[i+1])); + auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1])); Kokkos::Experimental::sort_team(t, patch); }; Kokkos::parallel_for(TeamPol(offsets.size()-1, Kokkos::AUTO()), segment_sort); - return Graph(offsets,elms_w);//need a LOs instead of elms_w... + return Graph(offsets,Write(elms)); } Graph remove_duplicate_edges(Graph g) { From c8fe9d9401e91e951707358a1f4370d521058c24 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 06:30:27 -0400 Subject: [PATCH 06/34] nodiscard attributes --- src/formPatches.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index d5277ca1..361ce08f 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -7,12 +7,12 @@ using namespace Omega_h; -Graph getElmToElm2ndOrderAdj(Mesh& m, Int bridgeDim) { +[[nodiscard]] Graph getElmToElm2ndOrderAdj(Mesh& m, Int bridgeDim) { OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); return Graph(); } -bool patchSufficient(Graph patches, Int minPatchSize) { +[[nodiscard]] bool patchSufficient(Graph patches, Int minPatchSize) { //find min degree for each patch //if( minDegree < minPatchSize) // return false; @@ -20,7 +20,7 @@ bool patchSufficient(Graph patches, Int minPatchSize) { return true; } -Graph adj_segment_sort(Graph& g) { +[[nodiscard]] Graph adj_segment_sort(Graph& g) { using ExecSpace = Kokkos::DefaultExecutionSpace; using TeamPol = Kokkos::TeamPolicy; using TeamMem = typename TeamPol::member_type; @@ -38,7 +38,7 @@ Graph adj_segment_sort(Graph& g) { return Graph(offsets,Write(elms)); } -Graph remove_duplicate_edges(Graph g) { +[[nodiscard]] Graph remove_duplicate_edges(Graph g) { return Graph(); } @@ -50,13 +50,13 @@ Graph remove_duplicate_edges(Graph g) { * order element-to-element adjacencies * \return an expanded graph from key entities to elements */ -Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { +[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); auto patch_elms = patches.ab2b; auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); auto expanded = unmap_graph(patch_elms, adjElms); - adj_segment_sort(expanded); - auto dedup = remove_duplicate_edges(expanded); + auto sorted = adj_segment_sort(expanded); + auto dedup = remove_duplicate_edges(sorted); return dedup; } @@ -70,7 +70,7 @@ Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { * \return a graph whose source nodes are the entities of keyDim dimension, and * edges are connecting to elements in the patch */ -Graph formPatches(Mesh& m, LO keyDim, Int minPatchSize) { +[[nodiscard]] Graph formPatches(Mesh& m, LO keyDim, Int minPatchSize) { OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); OMEGA_H_CHECK(minPatchSize > 0); auto patches = Graph(); From be07aeabcda105c25dead9015a64781e4340539c Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 06:38:32 -0400 Subject: [PATCH 07/34] patchSufficient compiles --- src/formPatches.cpp | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 361ce08f..ade28c7d 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -1,7 +1,9 @@ #include #include #include -#include +#include //parallel_for +#include //get_min +#include //sort_team #include @@ -13,10 +15,16 @@ using namespace Omega_h; } [[nodiscard]] bool patchSufficient(Graph patches, Int minPatchSize) { - //find min degree for each patch - //if( minDegree < minPatchSize) - // return false; - //else + const auto num_patches = patches.a2ab.size()-1; + auto offsets = patches.a2ab; + Write degree(num_patches); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { + degree[i] = offsets[i+1]-offsets[i]; + }); + auto minDegree = get_min(read(degree)); + if( minDegree < minPatchSize) + return false; + else return true; } From f5e129d9b7f1c34233989d415af8ef8a09cfdaeb Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 07:09:07 -0400 Subject: [PATCH 08/34] 2nd order adj wip --- src/formPatches.cpp | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index ade28c7d..405be7c1 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -10,7 +10,24 @@ using namespace Omega_h; [[nodiscard]] Graph getElmToElm2ndOrderAdj(Mesh& m, Int bridgeDim) { - OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); + const auto elmDim = m.dim(); + OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < elmDim); + auto elm2Bridge = m.ask_down(elmDim, bridgeDim); + auto e2bOffsets = elm2Bridge.a2ab; + auto e2bValues = elm2Bridge.ab2b; + //get bridge-to-elm + auto bridge2Elm = m.ask_up(bridgeDim, elmDim); + //traverse elm-to-bridge and bridge-to-elm to form elm-to-elm + //are there map and graph functions to compute this? + //first count how many elms there are per elm + Write e2eDegree(m.nelems()); + auto b2eOffsets = bridge2Elm.a2ab; + parallel_for(elm, OMEGA_H_LAMBDA(LO i) { + //HERE loop over e2bOffsets + //get bridge + e2eDegree[i] = b2eOffsets[bridge+1]-b2eOffsets[bridge]; //counts duplicates! + }); + return Graph(); } From 29da527eb1bc906bbfb322e4dfc224cbe7b959b3 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 15:55:04 -0400 Subject: [PATCH 09/34] remove duplicates compiles --- src/formPatches.cpp | 67 ++++++++++++++++++++++++++------------------- 1 file changed, 39 insertions(+), 28 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 405be7c1..0a017576 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -9,9 +9,42 @@ using namespace Omega_h; +[[nodiscard]] Graph adj_segment_sort(Graph& g) { + using ExecSpace = Kokkos::DefaultExecutionSpace; + using TeamPol = Kokkos::TeamPolicy; + using TeamMem = typename TeamPol::member_type; + auto offsets = g.a2ab; + auto elms_r = g.ab2b.view(); //read only + Kokkos::View elms("elms", elms.size()); + Kokkos::deep_copy(elms, elms_r); + auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { + auto i = t.league_rank(); + auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1])); + Kokkos::Experimental::sort_team(t, patch); + }; + Kokkos::parallel_for(TeamPol(offsets.size()-1, Kokkos::AUTO()), segment_sort); + return Graph(offsets,Write(elms)); +} + +[[nodiscard]] Graph remove_duplicate_edges(Graph g) { + auto offsets = g.a2ab; + auto values = g.ab2b; + Write keep(values.size()); + auto markDups = OMEGA_H_LAMBDA(LO i) { + keep[offsets[i]] = 1; + for(int j=offsets[i]+1; j= 0 && bridgeDim < elmDim); + if(bridgeDim == elmDim-1) + return m.ask_dual(); auto elm2Bridge = m.ask_down(elmDim, bridgeDim); auto e2bOffsets = elm2Bridge.a2ab; auto e2bValues = elm2Bridge.ab2b; @@ -22,12 +55,11 @@ using namespace Omega_h; //first count how many elms there are per elm Write e2eDegree(m.nelems()); auto b2eOffsets = bridge2Elm.a2ab; - parallel_for(elm, OMEGA_H_LAMBDA(LO i) { - //HERE loop over e2bOffsets - //get bridge - e2eDegree[i] = b2eOffsets[bridge+1]-b2eOffsets[bridge]; //counts duplicates! - }); - +// parallel_for(elm, OMEGA_H_LAMBDA(LO i) { +// //HERE loop over e2bOffsets +// //get bridge +// e2eDegree[i] = b2eOffsets[bridge+1]-b2eOffsets[bridge]; //counts duplicates! +// }); return Graph(); } @@ -45,27 +77,6 @@ using namespace Omega_h; return true; } -[[nodiscard]] Graph adj_segment_sort(Graph& g) { - using ExecSpace = Kokkos::DefaultExecutionSpace; - using TeamPol = Kokkos::TeamPolicy; - using TeamMem = typename TeamPol::member_type; - auto offsets = g.a2ab; - auto elms_r = g.ab2b.view(); //read only - Kokkos::View elms("elms", elms.size()); - Kokkos::deep_copy(elms, elms_r); - auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { - //Sort a row of A using the whole team. - auto i = t.league_rank(); - auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1])); - Kokkos::Experimental::sort_team(t, patch); - }; - Kokkos::parallel_for(TeamPol(offsets.size()-1, Kokkos::AUTO()), segment_sort); - return Graph(offsets,Write(elms)); -} - -[[nodiscard]] Graph remove_duplicate_edges(Graph g) { - return Graph(); -} /** * \brief expand the patches @@ -98,7 +109,7 @@ using namespace Omega_h; [[nodiscard]] Graph formPatches(Mesh& m, LO keyDim, Int minPatchSize) { OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); OMEGA_H_CHECK(minPatchSize > 0); - auto patches = Graph(); + auto patches = m.ask_up(VERT,m.dim()); for(Int bridgeDim = m.dim()-1; bridgeDim >= 0; bridgeDim--) { auto bridgePatches = expandPatches(m, patches, bridgeDim); if( patchSufficient(bridgePatches, minPatchSize) ) { From 25baf30f431ad218a1d697943bbaf050285795e7 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 16:29:04 -0400 Subject: [PATCH 10/34] 2nd order adj compiles --- src/formPatches.cpp | 44 ++++++++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 0a017576..d6865074 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -3,13 +3,14 @@ #include #include //parallel_for #include //get_min +#include //offset_scan #include //sort_team #include using namespace Omega_h; -[[nodiscard]] Graph adj_segment_sort(Graph& g) { +[[nodiscard]] Graph adj_segment_sort(Graph& g) { //not tested using ExecSpace = Kokkos::DefaultExecutionSpace; using TeamPol = Kokkos::TeamPolicy; using TeamMem = typename TeamPol::member_type; @@ -26,7 +27,7 @@ using namespace Omega_h; return Graph(offsets,Write(elms)); } -[[nodiscard]] Graph remove_duplicate_edges(Graph g) { +[[nodiscard]] Graph remove_duplicate_edges(Graph g) { //not tested auto offsets = g.a2ab; auto values = g.ab2b; Write keep(values.size()); @@ -50,17 +51,32 @@ using namespace Omega_h; auto e2bValues = elm2Bridge.ab2b; //get bridge-to-elm auto bridge2Elm = m.ask_up(bridgeDim, elmDim); - //traverse elm-to-bridge and bridge-to-elm to form elm-to-elm - //are there map and graph functions to compute this? - //first count how many elms there are per elm - Write e2eDegree(m.nelems()); auto b2eOffsets = bridge2Elm.a2ab; -// parallel_for(elm, OMEGA_H_LAMBDA(LO i) { -// //HERE loop over e2bOffsets -// //get bridge -// e2eDegree[i] = b2eOffsets[bridge+1]-b2eOffsets[bridge]; //counts duplicates! -// }); - return Graph(); + auto b2eValues = bridge2Elm.ab2b; + //traverse elm-to-bridge and bridge-to-elm to form elm-to-elm + //first count how many adj elms there are per elm + Write e2eDegree(m.nelems(),0); + parallel_for(m.nelems(), OMEGA_H_LAMBDA(LO i) { + for(int j=e2bOffsets[i]; j e2eValues(e2eOffsets.last()); + parallel_for(m.nelems(), OMEGA_H_LAMBDA(LO i) { + auto pos = e2eOffsets[i]; + for(int j=e2bOffsets[i]; j= 0 && bridgeDim < m.dim()); auto patch_elms = patches.ab2b; From 55fa0fa263a7a52d153ac03aab1e9bd182eb1139 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 17 Aug 2024 16:46:02 -0400 Subject: [PATCH 11/34] down adj doesn't have an offset array --- src/formPatches.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index d6865074..4458f81a 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -4,6 +4,7 @@ #include //parallel_for #include //get_min #include //offset_scan +#include //simplex_degree #include //sort_team #include @@ -16,7 +17,7 @@ using namespace Omega_h; using TeamMem = typename TeamPol::member_type; auto offsets = g.a2ab; auto elms_r = g.ab2b.view(); //read only - Kokkos::View elms("elms", elms.size()); + Kokkos::View elms("elms", elms_r.size()); Kokkos::deep_copy(elms, elms_r); auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { auto i = t.league_rank(); @@ -47,7 +48,7 @@ using namespace Omega_h; if(bridgeDim == elmDim-1) return m.ask_dual(); auto elm2Bridge = m.ask_down(elmDim, bridgeDim); - auto e2bOffsets = elm2Bridge.a2ab; + auto e2bDegree = simplex_degree(elmDim, bridgeDim); auto e2bValues = elm2Bridge.ab2b; //get bridge-to-elm auto bridge2Elm = m.ask_up(bridgeDim, elmDim); @@ -57,7 +58,7 @@ using namespace Omega_h; //first count how many adj elms there are per elm Write e2eDegree(m.nelems(),0); parallel_for(m.nelems(), OMEGA_H_LAMBDA(LO i) { - for(int j=e2bOffsets[i]; j e2eValues(e2eOffsets.last()); parallel_for(m.nelems(), OMEGA_H_LAMBDA(LO i) { auto pos = e2eOffsets[i]; - for(int j=e2bOffsets[i]; j Date: Sun, 18 Aug 2024 13:02:24 -0400 Subject: [PATCH 12/34] rework expand patches The elmsInPatch-to-2ndAdjElms graph was being created and returned. It was not 'merging' those subgraphs into the patch-to-elmsInPatch graph that needs to be returned. --- src/formPatches.cpp | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 4458f81a..dadf32a0 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -108,8 +108,23 @@ using namespace Omega_h; auto patch_elms = patches.ab2b; auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); auto expanded = unmap_graph(patch_elms, adjElms); - auto sorted = adj_segment_sort(expanded); - auto dedup = remove_duplicate_edges(sorted); + //patchDegree = 0 + //for each patch (i.e., formed by a vertex) + // loop over patches.offsets + // patchElm = patches.values[j] + // patchDegree += expanded.offsets[j+1]-expanded.offsets[j] + //offsets = offset_scan(patchDegree) + //Write mergedPatches(offsets.last()) + //for each patch (i.e., formed by a vertex) + // idx = 0 + // loop over patches.offsets + // patchElm = patches.values[j] + // loop over expanded.offsets[patchElm] + // mergedPatches[idx] = expanded.values[k] + // idx++; + // Graph merged(offsets,mergedPatches) + // auto sorted = adj_segment_sort(merged); + // auto dedup = remove_duplicate_edges(sorted); return dedup; } @@ -128,7 +143,7 @@ using namespace Omega_h; OMEGA_H_CHECK(minPatchSize > 0); auto patches = m.ask_up(VERT,m.dim()); for(Int bridgeDim = m.dim()-1; bridgeDim >= 0; bridgeDim--) { - auto bridgePatches = expandPatches(m, patches, bridgeDim); + auto bridgePatches = expandPatches(m, patches, bridgeDim); //FIXME a2ab.size is 25... something wrong here, should be 10 if( patchSufficient(bridgePatches, minPatchSize) ) { return bridgePatches; } From aac6cf0112917849bbef56c23fdf70ef20dcdc47 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 19 Aug 2024 12:40:53 -0400 Subject: [PATCH 13/34] less broken expand patches use graph functions nnodes and nedges --- src/formPatches.cpp | 70 ++++++++++++++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 23 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index dadf32a0..ac0bdff6 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -24,20 +24,21 @@ using namespace Omega_h; auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1])); Kokkos::Experimental::sort_team(t, patch); }; - Kokkos::parallel_for(TeamPol(offsets.size()-1, Kokkos::AUTO()), segment_sort); + Kokkos::parallel_for(TeamPol(g.nnodes(), Kokkos::AUTO()), segment_sort); return Graph(offsets,Write(elms)); } [[nodiscard]] Graph remove_duplicate_edges(Graph g) { //not tested auto offsets = g.a2ab; auto values = g.ab2b; - Write keep(values.size()); + Write keep(g.nedges()); auto markDups = OMEGA_H_LAMBDA(LO i) { keep[offsets[i]] = 1; for(int j=offsets[i]+1; j degree(num_patches); parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { @@ -103,28 +104,37 @@ using namespace Omega_h; * order element-to-element adjacencies * \return an expanded graph from key entities to elements */ +//FIXME ideally, this used the Omega_h_map and Omega_h_graph functions [[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); - auto patch_elms = patches.ab2b; auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); - auto expanded = unmap_graph(patch_elms, adjElms); - //patchDegree = 0 - //for each patch (i.e., formed by a vertex) - // loop over patches.offsets - // patchElm = patches.values[j] - // patchDegree += expanded.offsets[j+1]-expanded.offsets[j] - //offsets = offset_scan(patchDegree) - //Write mergedPatches(offsets.last()) - //for each patch (i.e., formed by a vertex) - // idx = 0 - // loop over patches.offsets - // patchElm = patches.values[j] - // loop over expanded.offsets[patchElm] - // mergedPatches[idx] = expanded.values[k] - // idx++; - // Graph merged(offsets,mergedPatches) - // auto sorted = adj_segment_sort(merged); - // auto dedup = remove_duplicate_edges(sorted); + auto adjElms_offsets = adjElms.a2ab; + auto adjElms_elms = adjElms.ab2b; + const auto num_patches = patches.nnodes(); + auto patch_offsets = patches.a2ab; + auto patch_elms = patches.ab2b; + Write degree(num_patches); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { + degree[patch] = 0; + for(int j=patch_offsets[patch]; j patchExpDup_elms(patchExpDup_offsets.last()); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { + auto idx = 0; + for(int j=patch_offsets[patch]; j= 0 && keyDim < m.dim()); OMEGA_H_CHECK(minPatchSize > 0); auto patches = m.ask_up(VERT,m.dim()); + if( patchSufficient(patches, minPatchSize) ) { + return patches; + } for(Int bridgeDim = m.dim()-1; bridgeDim >= 0; bridgeDim--) { - auto bridgePatches = expandPatches(m, patches, bridgeDim); //FIXME a2ab.size is 25... something wrong here, should be 10 + auto bridgePatches = expandPatches(m, patches, bridgeDim); if( patchSufficient(bridgePatches, minPatchSize) ) { return bridgePatches; } @@ -157,6 +170,17 @@ int main(int argc, char** argv) { OMEGA_H_CHECK(argc == 3); Mesh mesh(&lib); binary::read(argv[1], lib.world(), &mesh); + { + Graph g({0,2,6},{1,1,2,3,7,7}); + auto res = remove_duplicate_edges(g); + Graph expected({0,1,4},{1,2,3,7}); + OMEGA_H_CHECK(res == expected); + } + Graph g({0,6},{1,1,1,1,1,1}); + auto res = remove_duplicate_edges(g); + Graph expected({0,1},{1}); + OMEGA_H_CHECK(res == expected); + } auto patches = formPatches(mesh, VERT, 3); vtk::write_parallel(argv[2], &mesh, mesh.dim()); } From fe59d80e7b4d416bb3e13f37d02a3bf966d58e35 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 19 Aug 2024 12:49:48 -0400 Subject: [PATCH 14/34] missing bracket --- src/formPatches.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index ac0bdff6..7f33b7d6 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -176,6 +176,7 @@ int main(int argc, char** argv) { Graph expected({0,1,4},{1,2,3,7}); OMEGA_H_CHECK(res == expected); } + { Graph g({0,6},{1,1,1,1,1,1}); auto res = remove_duplicate_edges(g); Graph expected({0,1},{1}); From ad7e0567141507d057bea87c214769732ad7b8dd Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 19 Aug 2024 12:58:00 -0400 Subject: [PATCH 15/34] test graph sort --- src/formPatches.cpp | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 7f33b7d6..1492b958 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -11,7 +11,7 @@ using namespace Omega_h; -[[nodiscard]] Graph adj_segment_sort(Graph& g) { //not tested +[[nodiscard]] Graph adj_segment_sort(Graph& g) { using ExecSpace = Kokkos::DefaultExecutionSpace; using TeamPol = Kokkos::TeamPolicy; using TeamMem = typename TeamPol::member_type; @@ -28,7 +28,7 @@ using namespace Omega_h; return Graph(offsets,Write(elms)); } -[[nodiscard]] Graph remove_duplicate_edges(Graph g) { //not tested +[[nodiscard]] Graph remove_duplicate_edges(Graph g) { auto offsets = g.a2ab; auto values = g.ab2b; Write keep(g.nedges()); @@ -165,11 +165,22 @@ using namespace Omega_h; return Graph(); } -int main(int argc, char** argv) { - auto lib = Library(&argc, &argv); - OMEGA_H_CHECK(argc == 3); - Mesh mesh(&lib); - binary::read(argv[1], lib.world(), &mesh); +void testGraphSort() { + { + Graph g({0,6},{1,3,2,4,4,3}); + auto res = adj_segment_sort(g); + Graph expected({0,6},{1,2,3,3,4,4}); + OMEGA_H_CHECK(res == expected); + } + { + Graph g({0,4,6},{4,3,2,2,101,100}); + auto res = adj_segment_sort(g); + Graph expected({0,4,6},{2,2,3,4,100,101}); + OMEGA_H_CHECK(res == expected); + } +} + +void testGraphDuplicateRemoval() { { Graph g({0,2,6},{1,1,2,3,7,7}); auto res = remove_duplicate_edges(g); @@ -182,6 +193,15 @@ int main(int argc, char** argv) { Graph expected({0,1},{1}); OMEGA_H_CHECK(res == expected); } +} + +int main(int argc, char** argv) { + auto lib = Library(&argc, &argv); + testGraphSort(); + testGraphDuplicateRemoval(); + OMEGA_H_CHECK(argc == 3); + Mesh mesh(&lib); + binary::read(argv[1], lib.world(), &mesh); auto patches = formPatches(mesh, VERT, 3); vtk::write_parallel(argv[2], &mesh, mesh.dim()); } From f56c0a857620f2024ae9ce4345c8bd1a592037de Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 19 Aug 2024 12:59:17 -0400 Subject: [PATCH 16/34] sort test --- src/formPatches.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 1492b958..697772d1 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -173,9 +173,9 @@ void testGraphSort() { OMEGA_H_CHECK(res == expected); } { - Graph g({0,4,6},{4,3,2,2,101,100}); + Graph g({0,4,6},{4,3,2,2,1,0}); auto res = adj_segment_sort(g); - Graph expected({0,4,6},{2,2,3,4,100,101}); + Graph expected({0,4,6},{2,2,3,4,0,1}); OMEGA_H_CHECK(res == expected); } } From 67a064b207ace8f5b40243d7c495f3dc13e263c5 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 19 Aug 2024 16:43:53 -0400 Subject: [PATCH 17/34] debugging helper functions --- src/formPatches.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 697772d1..c0be7451 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -11,6 +11,30 @@ using namespace Omega_h; +void render(Mesh& m, Graph patches, std::string suffix) { + const auto num_patches = patches.nnodes(); + auto offsets = patches.a2ab; + Write degree(num_patches); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { + degree[i] = offsets[i+1]-offsets[i]; + }); + m.add_tag(VERT, "patchDegree", 1, read(degree)); + std::string name = "patch_" + suffix + ".vtk"; + vtk::write_parallel(name, &m, m.dim()); +} + +void writeGraph(Graph g) { + HostRead offsets(g.a2ab); + HostRead values(g.ab2b); + for(int node = 0; node < g.nnodes(); node++) { + std::cout << node << ": "; + for(int edge = offsets[node]; edge < offsets[node+1]; edge++) { + std::cout << values[edge] << " "; + } + std::cout << "\n"; + } +} + [[nodiscard]] Graph adj_segment_sort(Graph& g) { using ExecSpace = Kokkos::DefaultExecutionSpace; using TeamPol = Kokkos::TeamPolicy; @@ -133,6 +157,7 @@ using namespace Omega_h; } }); Graph patchExpDup(patchExpDup_offsets,patchExpDup_elms); + writeGraph(patchExpDup); auto sorted = adj_segment_sort(patchExpDup); auto dedup = remove_duplicate_edges(sorted); return dedup; @@ -152,11 +177,15 @@ using namespace Omega_h; OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); OMEGA_H_CHECK(minPatchSize > 0); auto patches = m.ask_up(VERT,m.dim()); + writeGraph(patches); + render(m,patches,"init"); + //TODO modify patchSufficient to return a mask of patches that need to be expanded if( patchSufficient(patches, minPatchSize) ) { return patches; } for(Int bridgeDim = m.dim()-1; bridgeDim >= 0; bridgeDim--) { auto bridgePatches = expandPatches(m, patches, bridgeDim); + render(m,bridgePatches,std::to_string(bridgeDim)); if( patchSufficient(bridgePatches, minPatchSize) ) { return bridgePatches; } From 1069db3aca408f03e005060821e2a5b2fd7c72f7 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 20 Aug 2024 10:20:44 -0400 Subject: [PATCH 18/34] fill the patch correctly --- src/formPatches.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index c0be7451..81ed57ee 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -148,7 +148,7 @@ void writeGraph(Graph g) { auto patchExpDup_offsets = offset_scan(read(degree)); Write patchExpDup_elms(patchExpDup_offsets.last()); parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { - auto idx = 0; + auto idx = patchExpDup_offsets[patch]; for(int j=patch_offsets[patch]; j Date: Tue, 20 Aug 2024 12:19:55 -0400 Subject: [PATCH 19/34] wip - don't expand patches that are of sufficient size --- src/formPatches.cpp | 57 +++++++++++++++++++++++++-------------------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 81ed57ee..c5ea8d85 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -105,18 +105,14 @@ void writeGraph(Graph g) { return remove_duplicate_edges(e2e_dups); } -[[nodiscard]] bool patchSufficient(Graph patches, Int minPatchSize) { +[[nodiscard]] Read patchSufficient(Graph patches, Int minPatchSize) { const auto num_patches = patches.nnodes(); auto offsets = patches.a2ab; - Write degree(num_patches); + Write done(num_patches); parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { - degree[i] = offsets[i+1]-offsets[i]; + done[i] = ((offsets[i+1]-offsets[i]) >= minPatchSize); }); - auto minDegree = get_min(read(degree)); - if( minDegree < minPatchSize) - return false; - else - return true; + return read(done); } @@ -129,7 +125,7 @@ void writeGraph(Graph g) { * \return an expanded graph from key entities to elements */ //FIXME ideally, this used the Omega_h_map and Omega_h_graph functions -[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim) { +[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim, Read patchDone) { OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); auto adjElms_offsets = adjElms.a2ab; @@ -139,27 +135,38 @@ void writeGraph(Graph g) { auto patch_elms = patches.ab2b; Write degree(num_patches); parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { - degree[patch] = 0; - for(int j=patch_offsets[patch]; j patchExpDup_elms(patchExpDup_offsets.last()); parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { auto idx = patchExpDup_offsets[patch]; - for(int j=patch_offsets[patch]; j= 0; bridgeDim--) { - auto bridgePatches = expandPatches(m, patches, bridgeDim); + std::cout << "expanding via bridge " << bridgeDim << "\n"; + auto bridgePatches = expandPatches(m, patches, bridgeDim, patchDone); render(m,bridgePatches,std::to_string(bridgeDim)); - if( patchSufficient(bridgePatches, minPatchSize) ) { - return bridgePatches; - } + patchDone = patchSufficient(bridgePatches, minPatchSize); + if( get_min(patchDone) == 1 ) + return patches; } assert(false); //fails here return Graph(); From ed1f6ae3629a8fd72c08aa3451e1c15f24f40f91 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 20 Aug 2024 13:59:02 -0400 Subject: [PATCH 20/34] works on 2x2 test mesh of tris --- src/formPatches.cpp | 43 +++++++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index c5ea8d85..045d7d86 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -23,9 +23,10 @@ void render(Mesh& m, Graph patches, std::string suffix) { vtk::write_parallel(name, &m, m.dim()); } -void writeGraph(Graph g) { +void writeGraph(Graph g, std::string name="") { HostRead offsets(g.a2ab); HostRead values(g.ab2b); + std::cout << "== " << name << " ==\n"; for(int node = 0; node < g.nnodes(); node++) { std::cout << node << ": "; for(int edge = offsets[node]; edge < offsets[node+1]; edge++) { @@ -115,19 +116,24 @@ void writeGraph(Graph g) { return read(done); } +template +void writeArray(T arr, std::string prefix) { + std::cout << prefix << "\n"; + for(int i=0; i< arr.size(); i++) { + std::cout << i << ": " << arr[i] << "\n"; + } +} /** * \brief expand the patches * \param m (in) mesh of simplices * \param patches (in) graph of key entities to elements - * \param bridgeDim (in) the entity dimension used for expansion via second - * order element-to-element adjacencies + * \param adjElms (in) second order element-to-element adjacencies + * used for expansion * \return an expanded graph from key entities to elements */ //FIXME ideally, this used the Omega_h_map and Omega_h_graph functions -[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Int bridgeDim, Read patchDone) { - OMEGA_H_CHECK(bridgeDim >= 0 && bridgeDim < m.dim()); - auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); +[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Graph adjElms, Read patchDone) { auto adjElms_offsets = adjElms.a2ab; auto adjElms_elms = adjElms.ab2b; const auto num_patches = patches.nnodes(); @@ -136,8 +142,6 @@ void writeGraph(Graph g) { Write degree(num_patches); parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { degree[patch] = patch_offsets[patch+1] - patch_offsets[patch]; - }); - parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { if(!patchDone[patch]) { for(int j=patch_offsets[patch]; j patchExpDup_elms(patchExpDup_offsets.last()); parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { auto idx = patchExpDup_offsets[patch]; - for(int j=patch_offsets[patch]; j= 0; bridgeDim--) { - std::cout << "expanding via bridge " << bridgeDim << "\n"; - auto bridgePatches = expandPatches(m, patches, bridgeDim, patchDone); - render(m,bridgePatches,std::to_string(bridgeDim)); - patchDone = patchSufficient(bridgePatches, minPatchSize); + const auto bridgeDim = m.dim()-1; + auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); + writeGraph(adjElms, "adjElms"); + for(Int iter = 0; iter < 10; iter++) { + std::cout << iter << " expanding via bridge " << bridgeDim << "\n"; + patches = expandPatches(m, patches, adjElms, patchDone); + render(m,patches,std::to_string(bridgeDim)); + patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) return patches; } From 5dd36fd0a6a7a81183f652b92ed27febc2987f62 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 20 Aug 2024 22:16:26 -0400 Subject: [PATCH 21/34] verbosity control --- src/formPatches.cpp | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 045d7d86..d89cd84d 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -11,6 +11,8 @@ using namespace Omega_h; +bool verbose = false; + void render(Mesh& m, Graph patches, std::string suffix) { const auto num_patches = patches.nnodes(); auto offsets = patches.a2ab; @@ -24,6 +26,7 @@ void render(Mesh& m, Graph patches, std::string suffix) { } void writeGraph(Graph g, std::string name="") { + if(verbose!=2) return; HostRead offsets(g.a2ab); HostRead values(g.ab2b); std::cout << "== " << name << " ==\n"; @@ -116,14 +119,6 @@ void writeGraph(Graph g, std::string name="") { return read(done); } -template -void writeArray(T arr, std::string prefix) { - std::cout << prefix << "\n"; - for(int i=0; i< arr.size(); i++) { - std::cout << i << ": " << arr[i] << "\n"; - } -} - /** * \brief expand the patches * \param m (in) mesh of simplices @@ -193,14 +188,16 @@ void writeArray(T arr, std::string prefix) { auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); writeGraph(adjElms, "adjElms"); for(Int iter = 0; iter < 10; iter++) { - std::cout << iter << " expanding via bridge " << bridgeDim << "\n"; + if(verbose==2) std::cout << iter << " expanding via bridge " << bridgeDim << "\n"; patches = expandPatches(m, patches, adjElms, patchDone); render(m,patches,std::to_string(bridgeDim)); patchDone = patchSufficient(patches, minPatchSize); - if( get_min(patchDone) == 1 ) + if( get_min(patchDone) == 1 ) { + if(verbose) std::cout << "iterations: " << iter << "\n"; return patches; + } } - assert(false); //fails here + assert(false); return Graph(); } @@ -238,9 +235,10 @@ int main(int argc, char** argv) { auto lib = Library(&argc, &argv); testGraphSort(); testGraphDuplicateRemoval(); - OMEGA_H_CHECK(argc == 3); + OMEGA_H_CHECK(argc == 4); Mesh mesh(&lib); binary::read(argv[1], lib.world(), &mesh); + verbose = std::stoi(argv[3]); auto patches = formPatches(mesh, VERT, 3); vtk::write_parallel(argv[2], &mesh, mesh.dim()); } From fb172743b7175ebcc9fda56646c018598d19d219 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 11:22:37 -0400 Subject: [PATCH 22/34] add box tests with expected results --- src/formPatches.cpp | 73 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 68 insertions(+), 5 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index d89cd84d..2adcc627 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -6,12 +6,13 @@ #include //offset_scan #include //simplex_degree #include //sort_team +#include //build_box #include using namespace Omega_h; -bool verbose = false; +int verbose = 0; void render(Mesh& m, Graph patches, std::string suffix) { const auto num_patches = patches.nnodes(); @@ -26,7 +27,7 @@ void render(Mesh& m, Graph patches, std::string suffix) { } void writeGraph(Graph g, std::string name="") { - if(verbose!=2) return; + if(verbose<2) return; HostRead offsets(g.a2ab); HostRead values(g.ab2b); std::cout << "== " << name << " ==\n"; @@ -37,6 +38,19 @@ void writeGraph(Graph g, std::string name="") { } std::cout << "\n"; } + if(verbose<3) return; + std::cout << "offsets = {"; + for(int i = 0; i < offsets.size(); i++) { + std::cout << offsets[i]; + if(i!=offsets.size()-1) std::cout << ","; + } + std::cout << "};\n"; + std::cout << "values = {"; + for(int i = 0; i < values.size(); i++) { + std::cout << values[i]; + if(i!=values.size()-1) std::cout << ","; + } + std::cout << "};\n"; } [[nodiscard]] Graph adj_segment_sort(Graph& g) { @@ -188,12 +202,12 @@ void writeGraph(Graph g, std::string name="") { auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); writeGraph(adjElms, "adjElms"); for(Int iter = 0; iter < 10; iter++) { - if(verbose==2) std::cout << iter << " expanding via bridge " << bridgeDim << "\n"; + if(verbose>=2) std::cout << iter << " expanding via bridge " << bridgeDim << "\n"; patches = expandPatches(m, patches, adjElms, patchDone); render(m,patches,std::to_string(bridgeDim)); patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) { - if(verbose) std::cout << "iterations: " << iter << "\n"; + if(verbose>=1) std::cout << "iterations: " << iter << "\n"; return patches; } } @@ -231,14 +245,63 @@ void testGraphDuplicateRemoval() { } } +void test2x2(Omega_h::Library& lib) { + auto world = lib.world(); + const auto x = 2.0; + const auto y = 2.0; + const auto z = 0.0; + const auto nx = 2; + const auto ny = 2; + const auto nz = 0; + const auto symmetric = false; + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + const auto minPatchSize = 3; + auto patches = formPatches(mesh, VERT, minPatchSize); + Graph expected( + {0,4,7,11,17,20,24,27,30,34}, + {0,1,2,6,1,2,4,1,2,4,5,0,1,2,3,5,6,1,4,5,1,3,5,6,3,6,7,0,6,7,0,3,6,7}); + OMEGA_H_CHECK(patches == expected); +} + +void test1x5(Omega_h::Library& lib) { + auto world = lib.world(); + const auto x = 1.0; + const auto y = 5.0; + const auto z = 0.0; + const auto nx = 1; + const auto ny = 5; + const auto nz = 0; + const auto symmetric = false; + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + { + const auto minPatchSize = 3; + auto patches = formPatches(mesh, VERT, minPatchSize); + Graph expected( + {0,3,6,9,12,15,18,21,24,27,30,33,36}, + {0,1,2,1,2,3,0,1,2,0,1,2,1,3,4,3,4,6,5,6,9,4,5,6,7,8,9,7,8,9,7,8,9,5,7,9}); + OMEGA_H_CHECK(patches == expected); + } + { + const auto minPatchSize = 4; + auto patches = formPatches(mesh, VERT, minPatchSize); + Graph expected( + {0,4,9,13,17,22,27,32,37,41,45,49,54}, + {0,1,2,3,0,1,2,3,4,0,1,2,3,0,1,2,3,1,2,3,4,6,1,3, + 4,5,6,4,5,6,7,9,3,4,5,6,9,5,7,8,9,5,7,8,9,5,7,8,9,5,6,7,8,9}); + OMEGA_H_CHECK(patches == expected); + } +} + int main(int argc, char** argv) { auto lib = Library(&argc, &argv); testGraphSort(); testGraphDuplicateRemoval(); OMEGA_H_CHECK(argc == 4); + verbose = std::stoi(argv[3]); + test2x2(lib); + test1x5(lib); Mesh mesh(&lib); binary::read(argv[1], lib.world(), &mesh); - verbose = std::stoi(argv[3]); auto patches = formPatches(mesh, VERT, 3); vtk::write_parallel(argv[2], &mesh, mesh.dim()); } From 5f6baceac9dac3f1c48c30a7523d8e6806f0e9bf Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 13:22:13 -0400 Subject: [PATCH 23/34] only using elm-to-elm via edges(2d) --- src/formPatches.cpp | 49 ++++++--------------------------------------- 1 file changed, 6 insertions(+), 43 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 2adcc627..89253bea 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -75,7 +75,7 @@ void writeGraph(Graph g, std::string name="") { auto values = g.ab2b; Write keep(g.nedges()); auto markDups = OMEGA_H_LAMBDA(LO i) { - keep[offsets[i]] = 1; + keep[offsets[i]] = 1; //always keep the first edge in the segment for(int j=offsets[i]+1; j= 0 && bridgeDim < elmDim); - if(bridgeDim == elmDim-1) - return m.ask_dual(); - auto elm2Bridge = m.ask_down(elmDim, bridgeDim); - auto e2bDegree = simplex_degree(elmDim, bridgeDim); - auto e2bValues = elm2Bridge.ab2b; - //get bridge-to-elm - auto bridge2Elm = m.ask_up(bridgeDim, elmDim); - auto b2eOffsets = bridge2Elm.a2ab; - auto b2eValues = bridge2Elm.ab2b; - //traverse elm-to-bridge and bridge-to-elm to form elm-to-elm - //first count how many adj elms there are per elm - Write e2eDegree(m.nelems(),0); - parallel_for(m.nelems(), OMEGA_H_LAMBDA(LO i) { - for(int j=i*e2bDegree; j<(i+1)*e2bDegree; j++) { - auto bridge = e2bValues[j]; - e2eDegree[i] += b2eOffsets[bridge+1]-b2eOffsets[bridge]; - } - }); - auto e2eOffsets = offset_scan(read(e2eDegree)); - Write e2eValues(e2eOffsets.last()); - parallel_for(m.nelems(), OMEGA_H_LAMBDA(LO i) { - auto pos = e2eOffsets[i]; - for(int j=i*e2bDegree; j<(i+1)*e2bDegree; j++) { - auto bridge = e2bValues[j]; - for(int k=b2eOffsets[bridge]; k patchSufficient(Graph patches, Int minPatchSize) { const auto num_patches = patches.nnodes(); auto offsets = patches.a2ab; @@ -182,6 +144,8 @@ void writeGraph(Graph g, std::string name="") { /** * \brief form a patch of at least minPatchSize elements surrounding each entity * of dimension keyDim + * \remark the patch is expanded via 2nd order adjacencies using meshDim-1 as + * the bridge entity (e.g., faces for 3d, edges for 2d) * \param m (in) mesh of simplices * \param keyDim (in) the dimension of mesh entities that the patches are * created around @@ -198,13 +162,12 @@ void writeGraph(Graph g, std::string name="") { auto patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) return patches; - const auto bridgeDim = m.dim()-1; - auto adjElms = getElmToElm2ndOrderAdj(m, bridgeDim); + auto adjElms = m.ask_dual(); writeGraph(adjElms, "adjElms"); for(Int iter = 0; iter < 10; iter++) { - if(verbose>=2) std::cout << iter << " expanding via bridge " << bridgeDim << "\n"; + if(verbose>=2) std::cout << iter << " expanding patch\n"; patches = expandPatches(m, patches, adjElms, patchDone); - render(m,patches,std::to_string(bridgeDim)); + render(m,patches,std::to_string(iter)); patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) { if(verbose>=1) std::cout << "iterations: " << iter << "\n"; From 12abc83b2d46d89dd9fe9e69f92aac1c1faa5ff0 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 13:49:54 -0400 Subject: [PATCH 24/34] remove debug helper functions --- src/formPatches.cpp | 43 ------------------------------------------- 1 file changed, 43 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 89253bea..c37746af 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -14,45 +14,6 @@ using namespace Omega_h; int verbose = 0; -void render(Mesh& m, Graph patches, std::string suffix) { - const auto num_patches = patches.nnodes(); - auto offsets = patches.a2ab; - Write degree(num_patches); - parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { - degree[i] = offsets[i+1]-offsets[i]; - }); - m.add_tag(VERT, "patchDegree", 1, read(degree)); - std::string name = "patch_" + suffix + ".vtk"; - vtk::write_parallel(name, &m, m.dim()); -} - -void writeGraph(Graph g, std::string name="") { - if(verbose<2) return; - HostRead offsets(g.a2ab); - HostRead values(g.ab2b); - std::cout << "== " << name << " ==\n"; - for(int node = 0; node < g.nnodes(); node++) { - std::cout << node << ": "; - for(int edge = offsets[node]; edge < offsets[node+1]; edge++) { - std::cout << values[edge] << " "; - } - std::cout << "\n"; - } - if(verbose<3) return; - std::cout << "offsets = {"; - for(int i = 0; i < offsets.size(); i++) { - std::cout << offsets[i]; - if(i!=offsets.size()-1) std::cout << ","; - } - std::cout << "};\n"; - std::cout << "values = {"; - for(int i = 0; i < values.size(); i++) { - std::cout << values[i]; - if(i!=values.size()-1) std::cout << ","; - } - std::cout << "};\n"; -} - [[nodiscard]] Graph adj_segment_sort(Graph& g) { using ExecSpace = Kokkos::DefaultExecutionSpace; using TeamPol = Kokkos::TeamPolicy; @@ -137,7 +98,6 @@ void writeGraph(Graph g, std::string name="") { Graph patchExpDup(patchExpDup_offsets,patchExpDup_elms); auto sorted = adj_segment_sort(patchExpDup); auto dedup = remove_duplicate_edges(sorted); - writeGraph(dedup, "dedup"); return dedup; } @@ -157,8 +117,6 @@ void writeGraph(Graph g, std::string name="") { OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); OMEGA_H_CHECK(minPatchSize > 0); auto patches = m.ask_up(VERT,m.dim()); - writeGraph(patches); - render(m,patches,"init"); auto patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) return patches; @@ -167,7 +125,6 @@ void writeGraph(Graph g, std::string name="") { for(Int iter = 0; iter < 10; iter++) { if(verbose>=2) std::cout << iter << " expanding patch\n"; patches = expandPatches(m, patches, adjElms, patchDone); - render(m,patches,std::to_string(iter)); patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) { if(verbose>=1) std::cout << "iterations: " << iter << "\n"; From 98f76bc736496fb98148af4848134dbb86d0c155 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 13:50:11 -0400 Subject: [PATCH 25/34] better upper bound on iterations --- src/formPatches.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index c37746af..1a17b70c 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -121,8 +121,10 @@ int verbose = 0; if( get_min(patchDone) == 1 ) return patches; auto adjElms = m.ask_dual(); - writeGraph(adjElms, "adjElms"); - for(Int iter = 0; iter < 10; iter++) { + //assuming each iteration adds at least one element then + //the minPatchSize is a conservative upper bound on the + //iteration count + for(Int iter = 0; iter < minPatchSize; iter++) { if(verbose>=2) std::cout << iter << " expanding patch\n"; patches = expandPatches(m, patches, adjElms, patchDone); patchDone = patchSufficient(patches, minPatchSize); From db8ce1636c5442b24d8071c7531cd43e17028036 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 13:58:27 -0400 Subject: [PATCH 26/34] rename, return empty graph if failure --- src/formPatches.cpp | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 1a17b70c..92b1c108 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -102,19 +102,18 @@ int verbose = 0; } /** - * \brief form a patch of at least minPatchSize elements surrounding each entity - * of dimension keyDim + * \brief form a patch of at least minPatchSize elements surrounding each mesh + * vertex * \remark the patch is expanded via 2nd order adjacencies using meshDim-1 as * the bridge entity (e.g., faces for 3d, edges for 2d) * \param m (in) mesh of simplices - * \param keyDim (in) the dimension of mesh entities that the patches are - * created around * \param minPatchSize (in) the minimum number of elements in each patch - * \return a graph whose source nodes are the entities of keyDim dimension, and - * edges are connecting to elements in the patch + * \return a graph whose source nodes are mesh vertices, and + * edges are connecting to elements in the patch + * OR + * an empty graph upon failure */ -[[nodiscard]] Graph formPatches(Mesh& m, LO keyDim, Int minPatchSize) { - OMEGA_H_CHECK(keyDim >= 0 && keyDim < m.dim()); +[[nodiscard]] Graph formVertexPatches(Mesh& m, Int minPatchSize) { OMEGA_H_CHECK(minPatchSize > 0); auto patches = m.ask_up(VERT,m.dim()); auto patchDone = patchSufficient(patches, minPatchSize); @@ -133,7 +132,6 @@ int verbose = 0; return patches; } } - assert(false); return Graph(); } @@ -178,7 +176,7 @@ void test2x2(Omega_h::Library& lib) { const auto symmetric = false; auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); const auto minPatchSize = 3; - auto patches = formPatches(mesh, VERT, minPatchSize); + auto patches = formVertexPatches(mesh, minPatchSize); Graph expected( {0,4,7,11,17,20,24,27,30,34}, {0,1,2,6,1,2,4,1,2,4,5,0,1,2,3,5,6,1,4,5,1,3,5,6,3,6,7,0,6,7,0,3,6,7}); @@ -197,7 +195,7 @@ void test1x5(Omega_h::Library& lib) { auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); { const auto minPatchSize = 3; - auto patches = formPatches(mesh, VERT, minPatchSize); + auto patches = formVertexPatches(mesh, minPatchSize); Graph expected( {0,3,6,9,12,15,18,21,24,27,30,33,36}, {0,1,2,1,2,3,0,1,2,0,1,2,1,3,4,3,4,6,5,6,9,4,5,6,7,8,9,7,8,9,7,8,9,5,7,9}); @@ -205,7 +203,7 @@ void test1x5(Omega_h::Library& lib) { } { const auto minPatchSize = 4; - auto patches = formPatches(mesh, VERT, minPatchSize); + auto patches = formVertexPatches(mesh, minPatchSize); Graph expected( {0,4,9,13,17,22,27,32,37,41,45,49,54}, {0,1,2,3,0,1,2,3,4,0,1,2,3,0,1,2,3,1,2,3,4,6,1,3, @@ -224,6 +222,7 @@ int main(int argc, char** argv) { test1x5(lib); Mesh mesh(&lib); binary::read(argv[1], lib.world(), &mesh); - auto patches = formPatches(mesh, VERT, 3); + auto patches = formVertexPatches(mesh, 3); + OMEGA_H_CHECK(patches.nnodes() != 0); vtk::write_parallel(argv[2], &mesh, mesh.dim()); } From 80675f5eb2125fd515d67dbe759bb44dad2f3dd4 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 14:55:13 -0400 Subject: [PATCH 27/34] remove verbosity flag --- src/formPatches.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/formPatches.cpp b/src/formPatches.cpp index 92b1c108..fefec1df 100644 --- a/src/formPatches.cpp +++ b/src/formPatches.cpp @@ -12,8 +12,6 @@ using namespace Omega_h; -int verbose = 0; - [[nodiscard]] Graph adj_segment_sort(Graph& g) { using ExecSpace = Kokkos::DefaultExecutionSpace; using TeamPol = Kokkos::TeamPolicy; @@ -124,11 +122,9 @@ int verbose = 0; //the minPatchSize is a conservative upper bound on the //iteration count for(Int iter = 0; iter < minPatchSize; iter++) { - if(verbose>=2) std::cout << iter << " expanding patch\n"; patches = expandPatches(m, patches, adjElms, patchDone); patchDone = patchSufficient(patches, minPatchSize); if( get_min(patchDone) == 1 ) { - if(verbose>=1) std::cout << "iterations: " << iter << "\n"; return patches; } } @@ -216,8 +212,7 @@ int main(int argc, char** argv) { auto lib = Library(&argc, &argv); testGraphSort(); testGraphDuplicateRemoval(); - OMEGA_H_CHECK(argc == 4); - verbose = std::stoi(argv[3]); + OMEGA_H_CHECK(argc == 3); test2x2(lib); test1x5(lib); Mesh mesh(&lib); From 6bbd0abbd6c9dd312173210c4d240b7b3ace1241 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 16:26:26 -0400 Subject: [PATCH 28/34] move and rename patch function to mesh class given the use of kokkos team sorting, the patch api is only supported when kokkos is enabled --- src/CMakeLists.txt | 3 +- src/Omega_h_mesh.hpp | 15 +++++ src/Omega_h_patches.cpp | 134 ++++++++++++++++++++++++++++++++++++++++ src/test_patches.cpp | 63 +++++++++++++++++++ 4 files changed, 214 insertions(+), 1 deletion(-) create mode 100644 src/Omega_h_patches.cpp create mode 100644 src/test_patches.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0582fc5f..f36473d8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -71,6 +71,7 @@ set(Omega_h_SOURCES Omega_h_owners.cpp Omega_h_parser.cpp Omega_h_parser_graph.cpp + Omega_h_patches.cpp Omega_h_pool.cpp Omega_h_print.cpp Omega_h_profile.cpp @@ -391,7 +392,7 @@ if(BUILD_TESTING) osh_add_exe(arrayops_test) osh_add_exe(sort_test) osh_add_exe(unit_math) - osh_add_exe(formPatches) + osh_add_exe(test_patches) if (Omega_h_USE_Kokkos) osh_add_exe(bbox_reduce_test) osh_add_exe(initKokkosAndLib) diff --git a/src/Omega_h_mesh.hpp b/src/Omega_h_mesh.hpp index 30dac633..18f75222 100644 --- a/src/Omega_h_mesh.hpp +++ b/src/Omega_h_mesh.hpp @@ -167,6 +167,21 @@ class Mesh { std::string const& name, Read array); friend class ScopedChangeRCFieldsToMesh; + /** + * \brief form a patch of at least minPatchSize elements surrounding each mesh + * vertex + * \remark the patch is expanded via 2nd order adjacencies using meshDim-1 as + * the bridge entity (e.g., faces for 3d, edges for 2d) + * \param m (in) mesh of simplices + * \param minPatchSize (in) the minimum number of elements in each patch + * \return a graph whose source nodes are mesh vertices, and + * edges are connecting to elements in the patch + * OR + * an empty graph upon failure + */ + [[nodiscard]] Graph get_vtx_patches(Int minPatchSize); + + private: bool change_all_rcFieldsToMesh(); bool change_all_rcFieldsTorc(); diff --git a/src/Omega_h_patches.cpp b/src/Omega_h_patches.cpp new file mode 100644 index 00000000..bffe978d --- /dev/null +++ b/src/Omega_h_patches.cpp @@ -0,0 +1,134 @@ +#include +#include +#include +#include //parallel_for +#include //get_min +#include //offset_scan +#include //simplex_degree +#include //sort_team +#include //build_box + +#include + +using namespace Omega_h; + +#if defined(OMEGA_H_USE_KOKKOS) + +namespace { +[[nodiscard]] Graph adj_segment_sort(Graph& g) { + using ExecSpace = Kokkos::DefaultExecutionSpace; + using TeamPol = Kokkos::TeamPolicy; + using TeamMem = typename TeamPol::member_type; + auto offsets = g.a2ab; + auto elms_r = g.ab2b.view(); //read only + Kokkos::View elms("elms", elms_r.size()); + Kokkos::deep_copy(elms, elms_r); + auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { + auto i = t.league_rank(); + auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1])); + Kokkos::Experimental::sort_team(t, patch); + }; + Kokkos::parallel_for(TeamPol(g.nnodes(), Kokkos::AUTO()), segment_sort); + return Graph(offsets,Write(elms)); +} + +[[nodiscard]] Graph remove_duplicate_edges(Graph g) { + auto offsets = g.a2ab; + auto values = g.ab2b; + Write keep(g.nedges()); + auto markDups = OMEGA_H_LAMBDA(LO i) { + keep[offsets[i]] = 1; //always keep the first edge in the segment + for(int j=offsets[i]+1; j patchSufficient(Graph patches, Int minPatchSize) { + const auto num_patches = patches.nnodes(); + auto offsets = patches.a2ab; + Write done(num_patches); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { + done[i] = ((offsets[i+1]-offsets[i]) >= minPatchSize); + }); + return read(done); +} + +/** + * \brief expand the patches + * \param m (in) mesh of simplices + * \param patches (in) graph of key entities to elements + * \param adjElms (in) second order element-to-element adjacencies + * used for expansion + * \return an expanded graph from key entities to elements +*/ +//FIXME ideally, this used the Omega_h_map and Omega_h_graph functions +[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Graph adjElms, Read patchDone) { + auto adjElms_offsets = adjElms.a2ab; + auto adjElms_elms = adjElms.ab2b; + const auto num_patches = patches.nnodes(); + auto patch_offsets = patches.a2ab; + auto patch_elms = patches.ab2b; + Write degree(num_patches); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { + degree[patch] = patch_offsets[patch+1] - patch_offsets[patch]; + if(!patchDone[patch]) { + for(int j=patch_offsets[patch]; j patchExpDup_elms(patchExpDup_offsets.last()); + parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { + auto idx = patchExpDup_offsets[patch]; + for(int j=patch_offsets[patch]; j 0); + auto patches = ask_up(VERT,dim()); + auto patchDone = patchSufficient(patches, minPatchSize); + if( get_min(patchDone) == 1 ) + return patches; + auto adjElms = ask_dual(); + //assuming each iteration adds at least one element then + //the minPatchSize is a conservative upper bound on the + //iteration count + for(Int iter = 0; iter < minPatchSize; iter++) { + patches = expandPatches(*this, patches, adjElms, patchDone); + patchDone = patchSufficient(patches, minPatchSize); + if( get_min(patchDone) == 1 ) { + return patches; + } + } + return Graph(); +} + +#else +[[nodiscard]] Graph Mesh::get_vtx_patches(Int minPatchSize) { + auto message = "Patch creation requires Kokkos support. Please rebuild with Kokkos enabled.\n"; + if( library_->self() == 0 ) { + std::cerr << message; + } + return Graph(); +} +#endif diff --git a/src/test_patches.cpp b/src/test_patches.cpp new file mode 100644 index 00000000..67daeaf9 --- /dev/null +++ b/src/test_patches.cpp @@ -0,0 +1,63 @@ +#include +#include +#include +#include //build_box + +#include + +using namespace Omega_h; + +void test2x2(Omega_h::Library& lib) { + auto world = lib.world(); + const auto x = 2.0; + const auto y = 2.0; + const auto z = 0.0; + const auto nx = 2; + const auto ny = 2; + const auto nz = 0; + const auto symmetric = false; + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + const auto minPatchSize = 3; + auto patches = mesh.get_vtx_patches(minPatchSize); + Graph expected( + {0,4,7,11,17,20,24,27,30,34}, + {0,1,2,6,1,2,4,1,2,4,5,0,1,2,3,5,6,1,4,5,1,3,5,6,3,6,7,0,6,7,0,3,6,7}); + OMEGA_H_CHECK(patches == expected); +} + +void test1x5(Omega_h::Library& lib) { + auto world = lib.world(); + const auto x = 1.0; + const auto y = 5.0; + const auto z = 0.0; + const auto nx = 1; + const auto ny = 5; + const auto nz = 0; + const auto symmetric = false; + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + { + const auto minPatchSize = 3; + auto patches = mesh.get_vtx_patches(minPatchSize); + Graph expected( + {0,3,6,9,12,15,18,21,24,27,30,33,36}, + {0,1,2,1,2,3,0,1,2,0,1,2,1,3,4,3,4,6,5,6,9,4,5,6,7,8,9,7,8,9,7,8,9,5,7,9}); + OMEGA_H_CHECK(patches == expected); + } + { + const auto minPatchSize = 4; + auto patches = mesh.get_vtx_patches(minPatchSize); + Graph expected( + {0,4,9,13,17,22,27,32,37,41,45,49,54}, + {0,1,2,3,0,1,2,3,4,0,1,2,3,0,1,2,3,1,2,3,4,6,1,3, + 4,5,6,4,5,6,7,9,3,4,5,6,9,5,7,8,9,5,7,8,9,5,7,8,9,5,6,7,8,9}); + OMEGA_H_CHECK(patches == expected); + } +} + +int main(int argc, char** argv) { + auto lib = Omega_h::Library(&argc, &argv); + OMEGA_H_CHECK(argc == 1); + test2x2(lib); + test1x5(lib); + return 0; +} From 62607769b0d89063672d8b2cccaf53939d99ce1a Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 21 Aug 2024 16:30:43 -0400 Subject: [PATCH 29/34] add ctest --- src/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f36473d8..ac084c61 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -392,10 +392,11 @@ if(BUILD_TESTING) osh_add_exe(arrayops_test) osh_add_exe(sort_test) osh_add_exe(unit_math) - osh_add_exe(test_patches) if (Omega_h_USE_Kokkos) osh_add_exe(bbox_reduce_test) osh_add_exe(initKokkosAndLib) + osh_add_exe(test_patches) + test_func(run_test_patches 1 ./test_patches) endif() set(TEST_EXES ${TEST_EXES} unit_math) test_basefunc(run_unit_math 1 ./unit_math) From 64c5bbabf4c831355157242e15c9f72c728612e2 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 30 Aug 2024 12:06:04 -0400 Subject: [PATCH 30/34] remove old test implementation moved to Omega_h_patches.[h|c]pp --- src/formPatches.cpp | 223 -------------------------------------------- 1 file changed, 223 deletions(-) delete mode 100644 src/formPatches.cpp diff --git a/src/formPatches.cpp b/src/formPatches.cpp deleted file mode 100644 index fefec1df..00000000 --- a/src/formPatches.cpp +++ /dev/null @@ -1,223 +0,0 @@ -#include -#include -#include -#include //parallel_for -#include //get_min -#include //offset_scan -#include //simplex_degree -#include //sort_team -#include //build_box - -#include - -using namespace Omega_h; - -[[nodiscard]] Graph adj_segment_sort(Graph& g) { - using ExecSpace = Kokkos::DefaultExecutionSpace; - using TeamPol = Kokkos::TeamPolicy; - using TeamMem = typename TeamPol::member_type; - auto offsets = g.a2ab; - auto elms_r = g.ab2b.view(); //read only - Kokkos::View elms("elms", elms_r.size()); - Kokkos::deep_copy(elms, elms_r); - auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) { - auto i = t.league_rank(); - auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1])); - Kokkos::Experimental::sort_team(t, patch); - }; - Kokkos::parallel_for(TeamPol(g.nnodes(), Kokkos::AUTO()), segment_sort); - return Graph(offsets,Write(elms)); -} - -[[nodiscard]] Graph remove_duplicate_edges(Graph g) { - auto offsets = g.a2ab; - auto values = g.ab2b; - Write keep(g.nedges()); - auto markDups = OMEGA_H_LAMBDA(LO i) { - keep[offsets[i]] = 1; //always keep the first edge in the segment - for(int j=offsets[i]+1; j patchSufficient(Graph patches, Int minPatchSize) { - const auto num_patches = patches.nnodes(); - auto offsets = patches.a2ab; - Write done(num_patches); - parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) { - done[i] = ((offsets[i+1]-offsets[i]) >= minPatchSize); - }); - return read(done); -} - -/** - * \brief expand the patches - * \param m (in) mesh of simplices - * \param patches (in) graph of key entities to elements - * \param adjElms (in) second order element-to-element adjacencies - * used for expansion - * \return an expanded graph from key entities to elements -*/ -//FIXME ideally, this used the Omega_h_map and Omega_h_graph functions -[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Graph adjElms, Read patchDone) { - auto adjElms_offsets = adjElms.a2ab; - auto adjElms_elms = adjElms.ab2b; - const auto num_patches = patches.nnodes(); - auto patch_offsets = patches.a2ab; - auto patch_elms = patches.ab2b; - Write degree(num_patches); - parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { - degree[patch] = patch_offsets[patch+1] - patch_offsets[patch]; - if(!patchDone[patch]) { - for(int j=patch_offsets[patch]; j patchExpDup_elms(patchExpDup_offsets.last()); - parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) { - auto idx = patchExpDup_offsets[patch]; - for(int j=patch_offsets[patch]; j 0); - auto patches = m.ask_up(VERT,m.dim()); - auto patchDone = patchSufficient(patches, minPatchSize); - if( get_min(patchDone) == 1 ) - return patches; - auto adjElms = m.ask_dual(); - //assuming each iteration adds at least one element then - //the minPatchSize is a conservative upper bound on the - //iteration count - for(Int iter = 0; iter < minPatchSize; iter++) { - patches = expandPatches(m, patches, adjElms, patchDone); - patchDone = patchSufficient(patches, minPatchSize); - if( get_min(patchDone) == 1 ) { - return patches; - } - } - return Graph(); -} - -void testGraphSort() { - { - Graph g({0,6},{1,3,2,4,4,3}); - auto res = adj_segment_sort(g); - Graph expected({0,6},{1,2,3,3,4,4}); - OMEGA_H_CHECK(res == expected); - } - { - Graph g({0,4,6},{4,3,2,2,1,0}); - auto res = adj_segment_sort(g); - Graph expected({0,4,6},{2,2,3,4,0,1}); - OMEGA_H_CHECK(res == expected); - } -} - -void testGraphDuplicateRemoval() { - { - Graph g({0,2,6},{1,1,2,3,7,7}); - auto res = remove_duplicate_edges(g); - Graph expected({0,1,4},{1,2,3,7}); - OMEGA_H_CHECK(res == expected); - } - { - Graph g({0,6},{1,1,1,1,1,1}); - auto res = remove_duplicate_edges(g); - Graph expected({0,1},{1}); - OMEGA_H_CHECK(res == expected); - } -} - -void test2x2(Omega_h::Library& lib) { - auto world = lib.world(); - const auto x = 2.0; - const auto y = 2.0; - const auto z = 0.0; - const auto nx = 2; - const auto ny = 2; - const auto nz = 0; - const auto symmetric = false; - auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); - const auto minPatchSize = 3; - auto patches = formVertexPatches(mesh, minPatchSize); - Graph expected( - {0,4,7,11,17,20,24,27,30,34}, - {0,1,2,6,1,2,4,1,2,4,5,0,1,2,3,5,6,1,4,5,1,3,5,6,3,6,7,0,6,7,0,3,6,7}); - OMEGA_H_CHECK(patches == expected); -} - -void test1x5(Omega_h::Library& lib) { - auto world = lib.world(); - const auto x = 1.0; - const auto y = 5.0; - const auto z = 0.0; - const auto nx = 1; - const auto ny = 5; - const auto nz = 0; - const auto symmetric = false; - auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); - { - const auto minPatchSize = 3; - auto patches = formVertexPatches(mesh, minPatchSize); - Graph expected( - {0,3,6,9,12,15,18,21,24,27,30,33,36}, - {0,1,2,1,2,3,0,1,2,0,1,2,1,3,4,3,4,6,5,6,9,4,5,6,7,8,9,7,8,9,7,8,9,5,7,9}); - OMEGA_H_CHECK(patches == expected); - } - { - const auto minPatchSize = 4; - auto patches = formVertexPatches(mesh, minPatchSize); - Graph expected( - {0,4,9,13,17,22,27,32,37,41,45,49,54}, - {0,1,2,3,0,1,2,3,4,0,1,2,3,0,1,2,3,1,2,3,4,6,1,3, - 4,5,6,4,5,6,7,9,3,4,5,6,9,5,7,8,9,5,7,8,9,5,7,8,9,5,6,7,8,9}); - OMEGA_H_CHECK(patches == expected); - } -} - -int main(int argc, char** argv) { - auto lib = Library(&argc, &argv); - testGraphSort(); - testGraphDuplicateRemoval(); - OMEGA_H_CHECK(argc == 3); - test2x2(lib); - test1x5(lib); - Mesh mesh(&lib); - binary::read(argv[1], lib.world(), &mesh); - auto patches = formVertexPatches(mesh, 3); - OMEGA_H_CHECK(patches.nnodes() != 0); - vtk::write_parallel(argv[2], &mesh, mesh.dim()); -} From c74c8f4156764cc99677d7b6bf61d729690c7dc0 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 30 Aug 2024 12:47:08 -0400 Subject: [PATCH 31/34] cleanup --- src/Omega_h_patches.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/Omega_h_patches.cpp b/src/Omega_h_patches.cpp index bffe978d..12215f42 100644 --- a/src/Omega_h_patches.cpp +++ b/src/Omega_h_patches.cpp @@ -1,18 +1,12 @@ -#include -#include #include #include //parallel_for #include //get_min #include //offset_scan -#include //simplex_degree -#include //sort_team -#include //build_box - -#include using namespace Omega_h; #if defined(OMEGA_H_USE_KOKKOS) +#include //sort_team namespace { [[nodiscard]] Graph adj_segment_sort(Graph& g) { @@ -65,7 +59,7 @@ namespace { * used for expansion * \return an expanded graph from key entities to elements */ -//FIXME ideally, this used the Omega_h_map and Omega_h_graph functions +//TODO use Omega_h_map and Omega_h_graph functions [[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Graph adjElms, Read patchDone) { auto adjElms_offsets = adjElms.a2ab; auto adjElms_elms = adjElms.ab2b; @@ -125,9 +119,11 @@ namespace { #else [[nodiscard]] Graph Mesh::get_vtx_patches(Int minPatchSize) { - auto message = "Patch creation requires Kokkos support. Please rebuild with Kokkos enabled.\n"; + const auto message = "get_vtx_patches requires Kokkos. Please rebuild with Kokkos enabled.\n"; + static_assert(false,message); if( library_->self() == 0 ) { std::cerr << message; + exit(EXIT_FAILURE); } return Graph(); } From 5577fc7e2139ccd1574ad669d19a70cb06cad34c Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 30 Aug 2024 13:06:16 -0400 Subject: [PATCH 32/34] don't define patch api when kokkos disabled --- src/Omega_h_mesh.hpp | 2 ++ src/Omega_h_patches.cpp | 11 ----------- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/src/Omega_h_mesh.hpp b/src/Omega_h_mesh.hpp index 18f75222..daba9e44 100644 --- a/src/Omega_h_mesh.hpp +++ b/src/Omega_h_mesh.hpp @@ -167,6 +167,7 @@ class Mesh { std::string const& name, Read array); friend class ScopedChangeRCFieldsToMesh; + #if defined(OMEGA_H_USE_KOKKOS) /** * \brief form a patch of at least minPatchSize elements surrounding each mesh * vertex @@ -180,6 +181,7 @@ class Mesh { * an empty graph upon failure */ [[nodiscard]] Graph get_vtx_patches(Int minPatchSize); + #endif private: diff --git a/src/Omega_h_patches.cpp b/src/Omega_h_patches.cpp index 12215f42..c2ce3316 100644 --- a/src/Omega_h_patches.cpp +++ b/src/Omega_h_patches.cpp @@ -116,15 +116,4 @@ namespace { } return Graph(); } - -#else -[[nodiscard]] Graph Mesh::get_vtx_patches(Int minPatchSize) { - const auto message = "get_vtx_patches requires Kokkos. Please rebuild with Kokkos enabled.\n"; - static_assert(false,message); - if( library_->self() == 0 ) { - std::cerr << message; - exit(EXIT_FAILURE); - } - return Graph(); -} #endif From f03a1b85a256f40a6c40969806fefc7bfbc0b853 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 30 Aug 2024 14:29:57 -0400 Subject: [PATCH 33/34] add 3d and mpi tests --- src/CMakeLists.txt | 6 ++++- src/test_patches.cpp | 52 +++++++++++++++++++++++++++++++++++++------- 2 files changed, 49 insertions(+), 9 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ac084c61..d08c8912 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -396,7 +396,11 @@ if(BUILD_TESTING) osh_add_exe(bbox_reduce_test) osh_add_exe(initKokkosAndLib) osh_add_exe(test_patches) - test_func(run_test_patches 1 ./test_patches) + if(Omega_h_USE_MPI) + test_func(run_test_patches_par 4 ./test_patches) + else() + test_func(run_test_patches 1 ./test_patches) + endif() endif() set(TEST_EXES ${TEST_EXES} unit_math) test_basefunc(run_unit_math 1 ./unit_math) diff --git a/src/test_patches.cpp b/src/test_patches.cpp index 67daeaf9..1f065159 100644 --- a/src/test_patches.cpp +++ b/src/test_patches.cpp @@ -7,8 +7,8 @@ using namespace Omega_h; -void test2x2(Omega_h::Library& lib) { - auto world = lib.world(); +void test2x2(Omega_h::CommPtr comm) { + OMEGA_H_CHECK(comm->size() == 1); const auto x = 2.0; const auto y = 2.0; const auto z = 0.0; @@ -16,7 +16,7 @@ void test2x2(Omega_h::Library& lib) { const auto ny = 2; const auto nz = 0; const auto symmetric = false; - auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); const auto minPatchSize = 3; auto patches = mesh.get_vtx_patches(minPatchSize); Graph expected( @@ -25,8 +25,8 @@ void test2x2(Omega_h::Library& lib) { OMEGA_H_CHECK(patches == expected); } -void test1x5(Omega_h::Library& lib) { - auto world = lib.world(); +void test1x5(Omega_h::CommPtr comm) { + OMEGA_H_CHECK(comm->size() == 1); const auto x = 1.0; const auto y = 5.0; const auto z = 0.0; @@ -34,7 +34,7 @@ void test1x5(Omega_h::Library& lib) { const auto ny = 5; const auto nz = 0; const auto symmetric = false; - auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); { const auto minPatchSize = 3; auto patches = mesh.get_vtx_patches(minPatchSize); @@ -54,10 +54,46 @@ void test1x5(Omega_h::Library& lib) { } } +void test3D(Omega_h::CommPtr comm) { + OMEGA_H_CHECK(comm->size() == 1); + const auto x = 1.0; + const auto y = 1.0; + const auto z = 1.0; + const auto nx = 2; + const auto ny = 2; + const auto nz = 2; + const auto symmetric = false; + auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + const auto minPatchSize = 3; + auto patches = mesh.get_vtx_patches(minPatchSize); +} + +void testPar(Omega_h::CommPtr comm) { + OMEGA_H_CHECK(comm->size() == 4); + const auto x = 1.0; + const auto y = 1.0; + const auto z = 1.0; + const auto nx = 4; + const auto ny = 4; + const auto nz = 4; + const auto symmetric = false; + auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); + mesh.set_parting(OMEGA_H_GHOSTED); + const auto minPatchSize = 4; + auto patches = mesh.get_vtx_patches(minPatchSize); +} + int main(int argc, char** argv) { auto lib = Omega_h::Library(&argc, &argv); + auto world = lib.world(); OMEGA_H_CHECK(argc == 1); - test2x2(lib); - test1x5(lib); + if (world->rank() == 0) { + test2x2(lib.self()); + test1x5(lib.self()); + test3D(lib.self()); + } + if (world->size() == 4) { + testPar(world); + } return 0; } From 1beea4c20407fcd87728ffb1b4741f3ac33e0161 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 30 Aug 2024 15:43:33 -0400 Subject: [PATCH 34/34] fix parallel test --- src/test_patches.cpp | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/test_patches.cpp b/src/test_patches.cpp index 1f065159..a282a3f8 100644 --- a/src/test_patches.cpp +++ b/src/test_patches.cpp @@ -3,8 +3,6 @@ #include #include //build_box -#include - using namespace Omega_h; void test2x2(Omega_h::CommPtr comm) { @@ -66,21 +64,35 @@ void test3D(Omega_h::CommPtr comm) { auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); const auto minPatchSize = 3; auto patches = mesh.get_vtx_patches(minPatchSize); + OMEGA_H_CHECK(patches.nnodes() != -1); } void testPar(Omega_h::CommPtr comm) { OMEGA_H_CHECK(comm->size() == 4); const auto x = 1.0; const auto y = 1.0; - const auto z = 1.0; + const auto z = 0.0; const auto nx = 4; const auto ny = 4; - const auto nz = 4; + const auto nz = 0; const auto symmetric = false; auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric); mesh.set_parting(OMEGA_H_GHOSTED); - const auto minPatchSize = 4; - auto patches = mesh.get_vtx_patches(minPatchSize); + { + const auto min_elems = comm->allreduce(mesh.nelems(), OMEGA_H_MIN); + const auto minPatchSize = min_elems+1; + auto patches = mesh.get_vtx_patches(minPatchSize); + if( mesh.nelems() < minPatchSize ) { + OMEGA_H_CHECK(patches.nnodes() == -1); //expected to fail + } else { + OMEGA_H_CHECK(patches.nnodes() != -1); //expected to pass + } + } + { + const auto minPatchSize = 6; + auto patches = mesh.get_vtx_patches(minPatchSize); + OMEGA_H_CHECK(patches.nnodes() != -1); + } } int main(int argc, char** argv) {