From c2801b8034e86eee469161f0fa4e03f97e78997f Mon Sep 17 00:00:00 2001 From: Carlos Brito Date: Tue, 21 Nov 2023 09:47:40 +0100 Subject: [PATCH] Come back to serial Connectivity --- examples/Geometry/Connectivity.cpp | 22 +- src/Rodin/Geometry/Connectivity.cpp | 348 ++++++++++------------------ src/Rodin/Geometry/Connectivity.h | 13 +- tests/benchmarks/CMakeLists.txt | 1 + tests/benchmarks/Connectivity.cpp | 82 +++++++ 5 files changed, 219 insertions(+), 247 deletions(-) create mode 100644 tests/benchmarks/Connectivity.cpp diff --git a/examples/Geometry/Connectivity.cpp b/examples/Geometry/Connectivity.cpp index 5748a761b..d4c74582e 100644 --- a/examples/Geometry/Connectivity.cpp +++ b/examples/Geometry/Connectivity.cpp @@ -4,32 +4,26 @@ * (See accompanying file LICENSE or copy at * https://www.boost.org/LICENSE_1_0.txt) */ +#include #include #include using namespace Rodin; using namespace Geometry; + int main(int, char**) { Mesh mesh; - mesh = mesh.UniformGrid(Polytope::Type::Triangle, 6, 6); + mesh = mesh.UniformGrid(Polytope::Type::Triangle, 1024, 1024); + + auto t0 = std::chrono::high_resolution_clock::now(); mesh.getConnectivity().compute(1, 2); mesh.getConnectivity().compute(2, 1); mesh.getConnectivity().compute(1, 1); - - mesh.save("miaow.mesh", IO::FileFormat::MEDIT); - - auto skin = mesh.skin(); - const auto& inc = skin.getConnectivity().getIncidence(1, 1); - for (size_t i = 0; i < inc.size(); i++) - { - std::cout << i << ": "; - for (const Index s : inc[i]) - std::cout << s << " "; - std::cout << std::endl; - } - skin.save("skin.mesh", IO::FileFormat::MEDIT); + auto t1 = std::chrono::high_resolution_clock::now(); + auto d = t1 - t0; + std::cout << std::chrono::duration_cast(d).count() << std::endl; return 0; } diff --git a/src/Rodin/Geometry/Connectivity.cpp b/src/Rodin/Geometry/Connectivity.cpp index cd68309fb..613e184a1 100644 --- a/src/Rodin/Geometry/Connectivity.cpp +++ b/src/Rodin/Geometry/Connectivity.cpp @@ -41,37 +41,20 @@ namespace Rodin::Geometry MeshConnectivity& MeshConnectivity::reserve(size_t d, size_t count) { assert(d < m_connectivity.size()); - m_index[d].write( - [&](auto& obj) - { - obj.left.rehash(count); - obj.right.rehash(count); - }); - m_geometry[d].write( - [&](auto& obj) - { - obj.reserve(count); - }); - m_connectivity[d][0].write( - [&](auto& obj) - { - obj.reserve(count); - }); + m_index[d].left.rehash(count); m_index[d].right.rehash(count); + m_geometry[d].reserve(count); + m_connectivity[d][0].reserve(count); return *this; } MeshConnectivity& MeshConnectivity::nodes(size_t count) { - m_count[0].write() = count; - m_gcount[Geometry::Polytope::Type::Point].write() = count; + m_count[0] = count; + m_gcount[Geometry::Polytope::Type::Point] = count; for (size_t i = 0; i < count; i++) { - m_index[0].write( - [&](auto& obj) - { - auto p = obj.left.insert({ IndexArray{{ i }}, i }); - assert(p.second); - }); + auto p = m_index[0].left.insert({ IndexArray{{ i }}, i }); + assert(p.second); } return *this; } @@ -83,31 +66,14 @@ namespace Rodin::Geometry const size_t d = Polytope::getGeometryDimension(t); assert(d > 0); assert(d <= m_maximalDimension); - std::pair res; - m_index[d].write( - [&](PolytopeIndex& obj) - { - size_t countd; - m_count[d].read([&](const auto& obj) { countd = obj; }); - res = obj.left.insert({ in, countd }); - }); - const auto& it = res.first; - const bool& inserted = res.second; + auto [it, inserted] = m_index[d].left.insert({ in, m_count[d] }); if (inserted) { - m_connectivity[d][0].write( - [&](auto& obj) - { - obj.emplace_back().insert_unique(it->first.begin(), it->first.end()); - }); - m_geometry[d].write( - [&](auto& obj) - { - obj.push_back(t); - }); - m_count[d].write([](auto& obj){ obj += 1; }); - m_gcount[t].write([](auto& obj){ obj += 1; }); - m_dirty[d][0].write() = false; + m_connectivity[d][0].emplace_back().insert_unique(it->first.begin(), it->first.end()); + m_geometry[d].push_back(t); + m_count[d] += 1; + m_gcount[t] += 1; + m_dirty[d][0] = false; } return *this; } @@ -119,63 +85,37 @@ namespace Rodin::Geometry const size_t d = Polytope::getGeometryDimension(t); assert(d > 0); assert(d <= m_maximalDimension); - - size_t countd; - m_count[d].read([&](const auto& obj) { countd = obj; }); - - std::pair res; - m_index[d].write( - [&](auto& obj) - { - res = obj.left.insert({ std::move(in), countd }); - }); - const auto& it = res.first; - const bool& inserted = res.second; + auto [it, inserted] = m_index[d].left.insert({ std::move(in), m_count[d] }); if (inserted) { - m_connectivity[d][0].write( - [&](auto& obj) - { - obj.emplace_back().insert_unique(it->first.begin(), it->first.end()); - }); - m_geometry[d].write( - [&](auto& obj) - { - obj.push_back(t); - }); - m_count[d].write([](auto& obj){ obj += 1; }); - m_gcount[t].write([](auto& obj){ obj += 1; }); - m_dirty[d][0].write() = false; + m_connectivity[d][0].emplace_back().insert_unique(it->first.begin(), it->first.end()); + m_geometry[d].push_back(t); + m_count[d] += 1; + m_gcount[t] += 1; + m_dirty[d][0] = false; } return *this; } const MeshConnectivity::PolytopeIndex& MeshConnectivity::getIndexMap(size_t dim) const { - return m_index[dim].read(); + return m_index[dim]; } const std::optional MeshConnectivity::getIndex(size_t dim, const IndexArray& key) const { - bool found; - PolytopeIndex::left_const_iterator it; - m_index[dim].read( - [&](const auto& obj) - { - it = obj.left.find(key); - found = (it == obj.left.end()); - }); - if (found) - return it->second; - else + auto it = m_index[dim].left.find(key); + if (it == m_index[dim].left.end()) return {}; + else + return it->second; } const Incidence& MeshConnectivity::getIncidence(size_t d, size_t dp) const { assert(d < m_connectivity.size()); assert(dp < m_connectivity[d].size()); - return m_connectivity[d][dp].read(); + return m_connectivity[d][dp]; } const IndexSet& MeshConnectivity::getIncidence( @@ -184,27 +124,25 @@ namespace Rodin::Geometry const auto& [d, dp] = p; assert(d < m_connectivity.size()); assert(dp < m_connectivity[d].size()); - assert(idx < m_connectivity[d][dp].read().size()); - return m_connectivity[d][dp].read()[idx]; + assert(idx < m_connectivity[d][dp].size()); + return m_connectivity[d][dp][idx]; } size_t MeshConnectivity::getCount(size_t dim) const { - return m_count[dim].read(); + return m_count[dim]; } size_t MeshConnectivity::getCount(Polytope::Type g) const { - return m_gcount[g].read(); + return m_gcount[g]; } size_t MeshConnectivity::getMeshDimension() const { for (int i = m_count.size() - 1; i >= 0; i--) { - bool b; - m_count[i].read([&](const auto& obj) { b = (obj > 0); }); - if (b) + if (m_count[i] > 0) return i; } return 0; @@ -215,16 +153,12 @@ namespace Rodin::Geometry if (d == 0) return Polytope::Type::Point; else - { - Polytope::Type t; - m_geometry[d].read([&](const auto& obj) { t = obj[idx]; }); - return t; - } + return m_geometry[d][idx]; } const Array& MeshConnectivity::getPolytope(size_t d, Index idx) const { - return m_index[d].read().right.at(idx); + return m_index[d].right.at(idx); } MeshConnectivity& MeshConnectivity::setIncidence(const std::pair& p, Incidence&& inc) @@ -232,7 +166,7 @@ namespace Rodin::Geometry const auto& [d, dp] = p; assert(d < m_connectivity.size()); assert(dp < m_connectivity[d].size()); - m_connectivity[d][dp].write() = std::move(inc); + m_connectivity[d][dp] = std::move(inc); return *this; } @@ -241,18 +175,18 @@ namespace Rodin::Geometry const size_t D = getMeshDimension(); if (d == D && dp == 0) return *this; - if (m_dirty[D][D].read()) + if (m_dirty[D][D]) transpose(0, D).intersection(D, D, 0); - assert(!m_dirty[D][D].read()); - if (d != D && d != 0 && (m_dirty[D][d].read() || m_dirty[d][0].read())) + assert(!m_dirty[D][D]); + if (d != D && d != 0 && (m_dirty[D][d] || m_dirty[d][0])) build(d); - assert(!m_dirty[D][d].read()); - assert(!m_dirty[d][0].read() || d == D || d == 0); - if (dp != D && dp != 0 && (m_dirty[D][dp].read() || m_dirty[dp][0].read())) + assert(!m_dirty[D][d]); + assert(!m_dirty[d][0] || d == D || d == 0); + if (dp != D && dp != 0 && (m_dirty[D][dp] || m_dirty[dp][0])) build(dp); - assert(!m_dirty[D][dp].read()); - assert(!m_dirty[dp][0].read() || dp == D || dp == 0); - if (m_dirty[d][dp].read()) + assert(!m_dirty[D][dp]); + assert(!m_dirty[dp][0] || dp == D || dp == 0); + if (m_dirty[d][dp]) { if (d < dp) { @@ -268,7 +202,7 @@ namespace Rodin::Geometry compute(d, dpp).compute(dpp, dp).intersection(d, dp, dpp); } } - m_dirty[d][dp].write() = false; + m_dirty[d][dp] = false; return *this; } @@ -277,50 +211,32 @@ namespace Rodin::Geometry const size_t D = getMeshDimension(); assert(d > 0); assert(d < D); - assert(!m_dirty[D][0].read()); - assert(!m_dirty[D][D].read()); - for (Index i = 0; i < m_count[D].read(); i++) + assert(!m_dirty[D][0]); + assert(!m_dirty[D][D]); + for (Index i = 0; i < m_count[D]; i++) { IndexSet s; std::vector subpolytopes; local(subpolytopes, d, i); - for (auto& gv : subpolytopes) + for (auto& [geometry, vertices] : subpolytopes) { - const auto& geometry = gv.geometry; - auto& vertices = gv.vertices; - size_t countd; - m_count[d].read([&](const auto& obj) { countd = obj; }); - std::pair insert; - m_index[d].write([&](auto& obj) { insert = obj.left.insert({ std::move(vertices), countd }); }); - const PolytopeIndex::left_const_iterator it = insert.first; + auto insert = m_index[d].left.insert({ std::move(vertices), m_count[d] }); + const PolytopeIndex::left_iterator it = insert.first; const bool inserted = insert.second; - const auto& v = it->first; - const auto& idx = it->second; + const auto& [v, idx] = *it; if (inserted) { - m_geometry[d].write( - [&](auto& obj) - { - obj.push_back(geometry); - }); - m_connectivity[d][0].write( - [&](auto& obj) - { - obj.emplace_back().insert_unique(v.begin(), v.end()); - }); + m_geometry[d].push_back(geometry); + m_connectivity[d][0].emplace_back().insert_unique(v.begin(), v.end()); } - m_count[d].write([&](auto& obj) { obj += inserted && !(d == D || d == 0); } ); - m_gcount[geometry].write([&](auto& obj) { obj += inserted && !(d == D || d == 0); }); + m_count[d] += inserted && !(d == D || d == 0); + m_gcount[geometry] += inserted && !(d == D || d == 0); s.insert(idx); } - m_connectivity[D][d].write( - [&](auto& obj) - { - obj.push_back(std::move(s)); - }); + m_connectivity[D][d].push_back(std::move(s)); } - m_dirty[D][d].write() = false; - m_dirty[d][0].write() = false; + m_dirty[D][d] = false; + m_dirty[d][0] = false; return *this; } @@ -329,75 +245,55 @@ namespace Rodin::Geometry assert(d < dp); assert(d < m_connectivity.size()); assert(dp < m_connectivity[d].size()); - size_t countd; - m_count[d].read([&](const auto& obj) { countd = obj; }); - m_connectivity[d][dp].write([&](auto& obj) { obj.resize(countd); }); - size_t countdp; - m_count[dp].read([&](const auto& obj) { countdp = obj; }); - for (Index j = 0; j < countdp; j++) + m_connectivity[d][dp].resize(m_count[d]); + for (Index j = 0; j < m_count[dp]; j++) { - const IndexSet* isj = nullptr; - m_connectivity[dp][d].read([&](const auto& obj) { isj = &obj[j]; }); - for (Index i : *isj) - m_connectivity[d][dp].write([&](auto& obj) { obj[i].insert(j); }); + for (Index i : m_connectivity[dp][d][j]) + m_connectivity[d][dp][i].insert(j); } - m_dirty[d][dp].write() = false; + m_dirty[d][dp] = false; return *this; } MeshConnectivity& MeshConnectivity::intersection(size_t d, size_t dp, size_t dpp) { assert(d >= dp); - size_t countd; - m_count[d].read([&](const auto& obj) { countd = obj; }); - m_connectivity[d][dp].write([&](auto& obj) { obj.resize(countd); }); - for (Index i = 0; i < countd; i++) + m_connectivity[d][dp].resize(m_count[d]); + for (Index i = 0; i < m_count[d]; i++) { - m_connectivity[d][dpp].read( - [&](const auto& obj) { assert(i < obj.size()); }); - const IndexSet* isi = nullptr; - m_connectivity[d][dpp].read([&](const auto& obj) { isi = &obj[i]; }); - for (Index k : *isi) + assert(i < m_connectivity[d][dpp].size()); + for (Index k : m_connectivity[d][dpp][i]) { - m_connectivity[dpp][dp].read( - [&](const auto& obj) { assert(k < obj.size()); }); - const IndexSet* isk = nullptr; - m_connectivity[dpp][dp].read([&](const auto& obj) { isk = &obj[k]; }); - for (Index j : *isk) + assert(k < m_connectivity[dpp][dp].size()); + for (Index j : m_connectivity[dpp][dp][k]) { - m_connectivity[d][0].read( - [&](const auto& obj) { assert(i < obj.size()); }); - m_connectivity[dp][0].read( - [&](const auto& obj) { assert(j < obj.size()); }); - const IndexSet* isd0i = nullptr; - m_connectivity[d][0].read([&](const auto& obj) { isd0i = &obj[i]; }); - const IndexSet* isdp0j = nullptr; - m_connectivity[dp][0].read([&](const auto& obj) { isdp0j = &obj[j]; }); + assert(i < m_connectivity[d][0].size()); + assert(j < m_connectivity[dp][0].size()); + const auto& d0i = m_connectivity[d][0][i]; + const auto& d0j = m_connectivity[dp][0][j]; if ((d == dp && i != j) || - (d > dp && std::includes(isd0i->begin(), isd0i->end(), isdp0j->begin(), isdp0j->end()))) + (d > dp && std::includes(d0i.begin(), d0i.end(), d0j.begin(), d0j.end()))) { - m_connectivity[d][dp].write([&](auto& obj) { obj[i].insert(j); }); + m_connectivity[d][dp][i].insert(j); } } } } - m_dirty[d][dp].write() = false; + m_dirty[d][dp] = false; return *this; } void MeshConnectivity::local(std::vector& out, size_t dim, Index i) { const size_t D = getMeshDimension(); - const IndexArray* p = nullptr; - m_index[D].read([&](const auto& obj) { p = &obj.right.at(i); }); - Polytope::Type geometryD; - m_geometry[D].read([&](const auto& obj) { geometryD = obj[i]; }); - switch (geometryD) + + const auto& p = m_index[D].right.at(i); + switch (m_geometry[D][i]) { case Polytope::Type::Point: { assert(dim == 0); - assert(p->size() == 0); + assert(p.size() == 0); out.resize(1); out[0].geometry = Polytope::Type::Point; out[0].vertices.resize(1); @@ -407,22 +303,22 @@ namespace Rodin::Geometry case Polytope::Type::Segment: { assert(dim <= 1); - assert(p->size() == 2); + assert(p.size() == 2); if (dim == 0) { out.resize(2); out[0].geometry = Polytope::Type::Point; out[0].vertices.resize(1); - out[0].vertices.coeffRef(0) = p->coeff(0); + out[0].vertices.coeffRef(0) = p.coeff(0); out[1].geometry = Polytope::Type::Point; out[1].vertices.resize(1); - out[1].vertices.coeffRef(0) = p->coeff(1); + out[1].vertices.coeffRef(0) = p.coeff(1); } else if (dim == 1) { out.resize(1); out[0].geometry = Polytope::Type::Segment; - out[0].vertices = *p; + out[0].vertices = p; } else { @@ -434,21 +330,21 @@ namespace Rodin::Geometry case Polytope::Type::Triangle: { assert(dim <= 2); - assert(p->size() == 3); + assert(p.size() == 3); if (dim == 0) { out.resize(3); out[0].geometry = Polytope::Type::Point; out[0].vertices.resize(1); - out[0].vertices.coeffRef(0) = p->coeff(0); + out[0].vertices.coeffRef(0) = p.coeff(0); out[1].geometry = Polytope::Type::Point; out[1].vertices.resize(1); - out[1].vertices.coeffRef(0) = p->coeff(1); + out[1].vertices.coeffRef(0) = p.coeff(1); out[2].geometry = Polytope::Type::Point; out[2].vertices.resize(1); - out[2].vertices.coeffRef(0) = p->coeff(2); + out[2].vertices.coeffRef(0) = p.coeff(2); } else if (dim == 1) { @@ -456,24 +352,24 @@ namespace Rodin::Geometry out[0].geometry = Polytope::Type::Segment; out[0].vertices.resize(2); - out[0].vertices.coeffRef(0) = p->coeff(0); - out[0].vertices.coeffRef(1) = p->coeff(1); + out[0].vertices.coeffRef(0) = p.coeff(0); + out[0].vertices.coeffRef(1) = p.coeff(1); out[1].geometry = Polytope::Type::Segment; out[1].vertices.resize(2); - out[1].vertices.coeffRef(0) = p->coeff(1); - out[1].vertices.coeffRef(1) = p->coeff(2); + out[1].vertices.coeffRef(0) = p.coeff(1); + out[1].vertices.coeffRef(1) = p.coeff(2); out[2].geometry = Polytope::Type::Segment; out[2].vertices.resize(2); - out[2].vertices.coeffRef(0) = p->coeff(0); - out[2].vertices.coeffRef(1) = p->coeff(2); + out[2].vertices.coeffRef(0) = p.coeff(0); + out[2].vertices.coeffRef(1) = p.coeff(2); } else if (dim == 2) { out.resize(1); out[0].geometry = Polytope::Type::Triangle; - out[0].vertices = *p; + out[0].vertices = p; } else { @@ -485,27 +381,27 @@ namespace Rodin::Geometry case Polytope::Type::Quadrilateral: { assert(dim <= 2); - assert(p->size() == 4); + assert(p.size() == 4); if (dim == 0) { out.resize(4); - out[0] = { Polytope::Type::Point, { { (*p)(0) } } }; - out[1] = { Polytope::Type::Point, { { (*p)(1) } } }; - out[2] = { Polytope::Type::Point, { { (*p)(2) } } }; - out[3] = { Polytope::Type::Point, { { (*p)(3) } } }; + out[0] = { Polytope::Type::Point, { { p(0) } } }; + out[1] = { Polytope::Type::Point, { { p(1) } } }; + out[2] = { Polytope::Type::Point, { { p(2) } } }; + out[3] = { Polytope::Type::Point, { { p(3) } } }; } else if (dim == 1) { out.resize(4); - out[0] = { Polytope::Type::Segment, { { (*p)(0), (*p)(1) } } }; - out[1] = { Polytope::Type::Segment, { { (*p)(1), (*p)(3) } } }; - out[2] = { Polytope::Type::Segment, { { (*p)(3), (*p)(2) } } }; - out[3] = { Polytope::Type::Segment, { { (*p)(2), (*p)(0) } } }; + out[0] = { Polytope::Type::Segment, { { p(0), p(1) } } }; + out[1] = { Polytope::Type::Segment, { { p(1), p(3) } } }; + out[2] = { Polytope::Type::Segment, { { p(3), p(2) } } }; + out[3] = { Polytope::Type::Segment, { { p(2), p(0) } } }; } else if (dim == 2) { out.resize(1); - out[0] = { Polytope::Type::Quadrilateral, *p }; + out[0] = { Polytope::Type::Quadrilateral, p }; } else { @@ -517,37 +413,37 @@ namespace Rodin::Geometry case Polytope::Type::Tetrahedron: { assert(dim <= 3); - assert(p->size() == 4); + assert(p.size() == 4); if (dim == 0) { out.resize(4); - out[0] = { Polytope::Type::Point, { { (*p)(0) } } }; - out[1] = { Polytope::Type::Point, { { (*p)(1) } } }; - out[2] = { Polytope::Type::Point, { { (*p)(2) } } }; - out[3] = { Polytope::Type::Point, { { (*p)(3) } } }; + out[0] = { Polytope::Type::Point, { { p(0) } } }; + out[1] = { Polytope::Type::Point, { { p(1) } } }; + out[2] = { Polytope::Type::Point, { { p(2) } } }; + out[3] = { Polytope::Type::Point, { { p(3) } } }; } else if (dim == 1) { out.resize(6); - out[0] = { Polytope::Type::Segment, { { (*p)(0), (*p)(1) } } }; - out[1] = { Polytope::Type::Segment, { { (*p)(1), (*p)(2) } } }; - out[2] = { Polytope::Type::Segment, { { (*p)(0), (*p)(2) } } }; - out[3] = { Polytope::Type::Segment, { { (*p)(0), (*p)(3) } } }; - out[4] = { Polytope::Type::Segment, { { (*p)(1), (*p)(3) } } }; - out[5] = { Polytope::Type::Segment, { { (*p)(2), (*p)(3) } } }; + out[0] = { Polytope::Type::Segment, { { p(0), p(1) } } }; + out[1] = { Polytope::Type::Segment, { { p(1), p(2) } } }; + out[2] = { Polytope::Type::Segment, { { p(0), p(2) } } }; + out[3] = { Polytope::Type::Segment, { { p(0), p(3) } } }; + out[4] = { Polytope::Type::Segment, { { p(1), p(3) } } }; + out[5] = { Polytope::Type::Segment, { { p(2), p(3) } } }; } else if (dim == 2) { out.resize(4); - out[0] = { Polytope::Type::Triangle, { { (*p)(0), (*p)(1), (*p)(2) } } }; - out[1] = { Polytope::Type::Triangle, { { (*p)(0), (*p)(1), (*p)(3) } } }; - out[2] = { Polytope::Type::Triangle, { { (*p)(0), (*p)(2), (*p)(3) } } }; - out[3] = { Polytope::Type::Triangle, { { (*p)(1), (*p)(2), (*p)(3) } } }; + out[0] = { Polytope::Type::Triangle, { { p(0), p(1), p(2) } } }; + out[1] = { Polytope::Type::Triangle, { { p(0), p(1), p(3) } } }; + out[2] = { Polytope::Type::Triangle, { { p(0), p(2), p(3) } } }; + out[3] = { Polytope::Type::Triangle, { { p(1), p(2), p(3) } } }; } else if (dim == 3) { out.resize(1); - out[0] = { Polytope::Type::Tetrahedron, *p }; + out[0] = { Polytope::Type::Tetrahedron, p }; } else { @@ -565,8 +461,8 @@ namespace Rodin::Geometry { assert(d < m_connectivity.size()); assert(dp < m_connectivity[d].size()); - m_dirty[d][dp].write() = true; - m_connectivity[d][dp].write().clear(); + m_dirty[d][dp] = true; + m_connectivity[d][dp].clear(); return *this; } } diff --git a/src/Rodin/Geometry/Connectivity.h b/src/Rodin/Geometry/Connectivity.h index 0c67a9a5a..083c383b8 100644 --- a/src/Rodin/Geometry/Connectivity.h +++ b/src/Rodin/Geometry/Connectivity.h @@ -15,7 +15,6 @@ #include #include "Rodin/Array.h" -#include "Rodin/Threads/Shared.h" #include "ForwardDecls.h" @@ -132,12 +131,12 @@ namespace Rodin::Geometry private: size_t m_maximalDimension; - std::vector> m_count; - GeometryIndexed> m_gcount; - std::vector> m_index; - std::vector>> m_dirty; - std::vector>> m_geometry; - std::vector>> m_connectivity; + std::vector m_count; + GeometryIndexed m_gcount; + std::vector m_index; + std::vector> m_dirty; + std::vector> m_geometry; + std::vector> m_connectivity; }; } diff --git a/tests/benchmarks/CMakeLists.txt b/tests/benchmarks/CMakeLists.txt index 78ef89e72..582da39fc 100644 --- a/tests/benchmarks/CMakeLists.txt +++ b/tests/benchmarks/CMakeLists.txt @@ -3,6 +3,7 @@ set(RodinBenchmarks_SRCS Poisson.cpp MeshIO.cpp UniformGrid.cpp + Connectivity.cpp ) add_executable(RodinBenchmarks ${RodinBenchmarks_SRCS}) diff --git a/tests/benchmarks/Connectivity.cpp b/tests/benchmarks/Connectivity.cpp new file mode 100644 index 000000000..9b627ffda --- /dev/null +++ b/tests/benchmarks/Connectivity.cpp @@ -0,0 +1,82 @@ +/* + * Copyright Carlos BRITO PACHECO 2021 - 2023. + * Distributed under the Boost Software License, Version 1.0. + * (See accompanying file LICENSE or copy at + * https://www.boost.org/LICENSE_1_0.txt) + */ +#include + +#include +#include + +using namespace Rodin; +using namespace Rodin::Geometry; +using namespace Rodin::Variational; + +namespace RodinBenchmark +{ + struct Connectivity : public benchmark::Fixture + { + public: + template + void test(Mesh& mesh, benchmark::State& st) + { + const size_t D = mesh.getDimension(); + for (auto _ : st) + { + st.ResumeTiming(); + for (size_t d = 0; d <= D; d++) + for (size_t dp = 0; dp <= D; dp++) + mesh.getConnectivity().compute(d, dp); + st.PauseTiming(); + for (size_t d = 0; d <= D; d++) + { + for (size_t dp = 0; dp <= D; dp++) + { + if (!(d == D && dp == 0)) + mesh.getConnectivity().clear(d, dp); + } + } + } + } + + void SetUp(const benchmark::State&) + {} + + void TearDown(const benchmark::State&) + {} + }; + + + BENCHMARK_F(Connectivity, Triangular_16x16)(benchmark::State& st) + { + auto mesh = SerialMesh::UniformGrid(Polytope::Type::Triangle, 16, 16); + test(mesh, st); + } + + BENCHMARK_F(Connectivity, Triangular_64x64)(benchmark::State& st) + { + auto mesh = SerialMesh::UniformGrid(Polytope::Type::Triangle, 32, 32); + test(mesh, st); + } + + BENCHMARK_F(Connectivity, Triangular_128x128)(benchmark::State& st) + { + auto mesh = SerialMesh::UniformGrid(Polytope::Type::Triangle, 128, 128); + test(mesh, st); + } + + + BENCHMARK_F(Connectivity, Triangular_256x256)(benchmark::State& st) + { + auto mesh = SerialMesh::UniformGrid(Polytope::Type::Triangle, 256, 256); + test(mesh, st); + } + + BENCHMARK_F(Connectivity, Triangular_512x512)(benchmark::State& st) + { + auto mesh = SerialMesh::UniformGrid(Polytope::Type::Triangle, 512, 512); + test(mesh, st); + } +} +