Skip to content

Commit

Permalink
Begin work on Threaded assembly
Browse files Browse the repository at this point in the history
  • Loading branch information
@cbritopacheco
cbritopacheco committed Nov 10, 2023
1 parent cce88cf commit fa5c939
Show file tree
Hide file tree
Showing 17 changed files with 543 additions and 177 deletions.
2 changes: 1 addition & 1 deletion examples/BoundaryOptimization/BoundaryHeatConductivity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ static constexpr Geometry::Attribute SigmaN = 2;
static constexpr size_t maxIt = 10000;

static constexpr Scalar epsilon = 0.001;
static constexpr Scalar ell = 5;
static constexpr Scalar ell = 1;
static constexpr Scalar radius = 0.02;
static constexpr Scalar tgv = std::numeric_limits<float>::max();

Expand Down
9 changes: 9 additions & 0 deletions examples/Models/Distance/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,12 @@
add_executable(FMM FMM.cpp)
target_link_libraries(FMM
PUBLIC
Rodin::Geometry
Rodin::Solver
Rodin::Variational
Rodin::External::MMG
)

add_executable(SpaldingTucker SpaldingTucker.cpp)
target_link_libraries(SpaldingTucker
PUBLIC
Expand Down
58 changes: 58 additions & 0 deletions examples/Models/Distance/FMM.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
/*
* Copyright Carlos BRITO PACHECO 2021 - 2022.
* 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 <Rodin/Math.h>
#include <Rodin/Solver.h>
#include <Rodin/Geometry.h>
#include <Rodin/Variational.h>
#include <Rodin/Models/Distance/FMM.h>

#include <RodinExternal/MMG.h>

using namespace Rodin;
using namespace Rodin::External;
using namespace Rodin::Geometry;
using namespace Rodin::Variational;

static constexpr Attribute boundary = 2;

int main(int, char**)
{
// Build a mesh
Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, 32, 32);
mesh.scale(1. / (31));
mesh.getConnectivity().compute(2, 0);

// Distance
{
P1 vh(mesh);
GridFunction dist(vh);
dist = [&](const Point& p)
{
Scalar d = (p - Math::SpatialVector{{0.75, 0.25}}).norm() - 0.05;
d = std::min(d, (p - Math::SpatialVector{{0.25, 0.25}}).norm() - 0.25);
d = std::min(d, (p - Math::SpatialVector{{0.75, 0.75}}).norm() - 0.1);
return d;
};

mesh = MMG::ImplicitDomainMesher().setBoundaryReference(5)
.setAngleDetection().setHMax(0.02).discretize(dist);
mesh.save("implicit.mesh", IO::FileFormat::MEDIT);
}

std::cout << "implicit done\n";
mesh.getConnectivity().compute(1, 2);

P1 vh(mesh);
auto dist = Models::Distance::FMM()(5, 3, vh);

dist.save("miaow.sol", IO::FileFormat::MEDIT);
mesh.save("miaow.mesh", IO::FileFormat::MEDIT);

return 0;
}

2 changes: 2 additions & 0 deletions src/Rodin/Configure.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,8 @@
*/
#cmakedefine RODIN_THIRD_PARTY_DIR "@RODIN_THIRD_PARTY_DIR@"

#cmakedefine RODIN_THREAD_SAFE

/**
* @ingroup RodinDirectives
* @brief Indicates the maximal space dimension.
Expand Down
139 changes: 90 additions & 49 deletions src/Rodin/Geometry/Mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,36 @@ namespace Rodin::Geometry
}

// ---- Mesh<Context::Serial> ----------------------------------------------
Mesh<Context::Serial>::Mesh(const Mesh& other)
: m_sdim(other.m_sdim),
m_vertices(other.m_vertices),
m_connectivity(other.m_connectivity),
m_attributeIndex(other.m_attributeIndex),
m_attributes(other.m_attributes)
{}

Mesh<Context::Serial>::Mesh(Mesh&& other)
: m_sdim(std::move(other.m_sdim)),
m_vertices(std::move(other.m_vertices)),
m_connectivity(std::move(other.m_connectivity)),
m_attributeIndex(std::move(other.m_attributeIndex)),
m_transformationIndex(std::move(other.m_transformationIndex)),
m_attributes(std::move(other.m_attributes))
{}

Mesh<Context::Serial>& Mesh<Context::Serial>::operator=(Mesh&& other)
{
Parent::operator=(std::move(other));
m_sdim = std::move(other.m_sdim);
m_vertices = std::move(other.m_vertices);
m_connectivity = std::move(other.m_connectivity);
m_attributeIndex = std::move(other.m_attributeIndex);
m_transformationIndex = std::move(other.m_transformationIndex);
m_attributes = std::move(other.m_attributes);
return *this;
}


Mesh<Context::Serial>&
Mesh<Context::Serial>::load(const boost::filesystem::path& filename, IO::FileFormat fmt)
{
Expand Down Expand Up @@ -174,28 +204,31 @@ namespace Rodin::Geometry
}

Mesh<Context::Serial>& Mesh<Context::Serial>::trace(
const Map<std::pair<Attribute, Attribute>, Attribute>& interface)
const Map<std::pair<Attribute, Attribute>, Attribute>& interface, const FlatSet<Attribute>& attrs)
{
const size_t D = getDimension();
RODIN_GEOMETRY_MESH_REQUIRE_INCIDENCE(D - 1, D);
const auto& conn = getConnectivity();
for (auto it = getFace(); it; ++it)
{
assert(it->getDimension() == D - 1);
const auto& inc = conn.getIncidence({ D - 1, D }, it->getIndex());
assert(inc.size() == 2);
auto el1 = getCell(*inc.begin());
auto el2 = getCell(*std::next(inc.begin()));
auto find = interface.find({ el1->getAttribute(), el2->getAttribute() });
if (find != interface.end())
if (attrs.size() == 0 || attrs.count(it->getAttribute()))
{
setAttribute({ D - 1, it->getIndex() }, find->second);
}
else
{
find = interface.find({ el2->getAttribute(), el1->getAttribute() });
assert(it->getDimension() == D - 1);
const auto& inc = conn.getIncidence({ D - 1, D }, it->getIndex());
assert(inc.size() == 2);
auto el1 = getCell(*inc.begin());
auto el2 = getCell(*std::next(inc.begin()));
auto find = interface.find({ el1->getAttribute(), el2->getAttribute() });
if (find != interface.end())
{
setAttribute({ D - 1, it->getIndex() }, find->second);
}
else
{
find = interface.find({ el2->getAttribute(), el1->getAttribute() });
if (find != interface.end())
setAttribute({ D - 1, it->getIndex() }, find->second);
}
}
}
return *this;
Expand Down Expand Up @@ -251,6 +284,42 @@ namespace Rodin::Geometry
return m_sdim;
}

Mesh<Context::Serial>& Mesh<Context::Serial>::setPolytopeTransformation(
const std::pair<size_t, Index> p, std::unique_ptr<PolytopeTransformation> trans)
{
m_transformationIndex.track(p, std::move(trans));
return *this;
}

PolytopeTransformation*
Mesh<Context::Serial>::getDefaultPolytopeTransformation(size_t dimension, Index idx) const
{
if (dimension == 0)
{
Variational::ScalarP1Element fe(Polytope::Type::Point);
const size_t sdim = getSpaceDimension();
Math::PointMatrix pm(sdim, 1);
pm.col(0) = getVertexCoordinates(idx);
return new IsoparametricTransformation(std::move(pm), std::move(fe));
}
else
{
auto g = getGeometry(dimension, idx);
const size_t sdim = getSpaceDimension();
const size_t n = Polytope::getVertexCount(g);
Math::PointMatrix pm(sdim, n);
const auto& polytope = getConnectivity().getPolytope(dimension, idx);
assert(n == static_cast<size_t>(polytope.size()));
for (const auto& v : polytope | boost::adaptors::indexed())
{
assert(sdim == static_cast<size_t>(getVertexCoordinates(v.value()).size()));
pm.col(v.index()) = getVertexCoordinates(v.value());
}
Variational::ScalarP1Element fe(g);
return new IsoparametricTransformation(std::move(pm), std::move(fe));
}
}

const PolytopeTransformation&
Mesh<Context::Serial>::getPolytopeTransformation(size_t dimension, Index idx) const
{
Expand All @@ -262,38 +331,10 @@ namespace Rodin::Geometry
}
else
{
if (dimension == 0)
{
Variational::ScalarP1Element fe(Polytope::Type::Point);
const size_t sdim = getSpaceDimension();
Math::PointMatrix pm(sdim, 1);
pm.col(0) = getVertexCoordinates(idx);
auto trans =
std::unique_ptr<PolytopeTransformation>(
new IsoparametricTransformation(std::move(pm), std::move(fe)));
auto p = m_transformationIndex.insert(it, { dimension, idx }, std::move(trans));
return *p->second;
}
else
{
auto g = getGeometry(dimension, idx);
const size_t sdim = getSpaceDimension();
const size_t n = Polytope::getVertexCount(g);
Math::PointMatrix pm(sdim, n);
const auto& polytope = getConnectivity().getPolytope(dimension, idx);
assert(n == static_cast<size_t>(polytope.size()));
for (const auto& v : polytope | boost::adaptors::indexed())
{
assert(sdim == static_cast<size_t>(getVertexCoordinates(v.value()).size()));
pm.col(v.index()) = getVertexCoordinates(v.value());
}
Variational::ScalarP1Element fe(g);
auto trans =
std::unique_ptr<PolytopeTransformation>(
new IsoparametricTransformation(std::move(pm), std::move(fe)));
auto p = m_transformationIndex.insert(it, { dimension, idx }, std::move(trans));
return *p->second;
}
PolytopeTransformation* trans = getDefaultPolytopeTransformation(dimension, idx);
auto p = m_transformationIndex.insert(
it, { dimension, idx }, std::unique_ptr<PolytopeTransformation>(trans));
return *p->second;
}
}

Expand Down Expand Up @@ -341,7 +382,7 @@ namespace Rodin::Geometry
const FlatSet<Attribute>& attrs) const
{
FlatSet<Index> visited;
visited.reserve(getCount(d));
visited.reserve(getPolytopeCount(d));
std::deque<Index> searchQueue;
std::deque<FlatSet<Index>> res;

Expand Down Expand Up @@ -381,12 +422,12 @@ namespace Rodin::Geometry
return res;
}

size_t Mesh<Context::Serial>::getCount(size_t dimension) const
size_t Mesh<Context::Serial>::getPolytopeCount(size_t dimension) const
{
return m_connectivity.getCount(dimension);
}

size_t Mesh<Context::Serial>::getCount(Polytope::Type g) const
size_t Mesh<Context::Serial>::getPolytopeCount(Polytope::Type g) const
{
return m_connectivity.getCount(g);
}
Expand Down Expand Up @@ -437,7 +478,7 @@ namespace Rodin::Geometry

PolytopeIterator Mesh<Context::Serial>::getPolytope(size_t dimension, Index idx) const
{
return PolytopeIterator(dimension, *this, BoundedIndexGenerator(idx, getCount(dimension)));
return PolytopeIterator(dimension, *this, BoundedIndexGenerator(idx, getPolytopeCount(dimension)));
}

bool Mesh<Context::Serial>::isInterface(Index faceIdx) const
Expand Down
Loading

0 comments on commit fa5c939

Please sign in to comment.