Skip to content

Commit

Permalink
Fix Potential
Browse files Browse the repository at this point in the history
  • Loading branch information
@cbritopacheco
cbritopacheco committed Aug 7, 2024
1 parent c651554 commit ad25629
Show file tree
Hide file tree
Showing 8 changed files with 719 additions and 895 deletions.
Binary file added doc/Images/UnitTriangle.afdesign
View file
Binary file not shown.
Binary file added doc/Images/UnitTriangle.jpg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
671 changes: 338 additions & 333 deletions examples/BoundaryOptimization/AcousticCloaking.cpp
View file

Large diffs are not rendered by default.

6 changes: 0 additions & 6 deletions examples/IntegralEquations/NewtonianPotential.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,6 @@ int main(int, char**)
std::cout << "assemblage\n";
eq.assemble();

std::cout << "save\n";
std::ofstream file("matrix.csv");
file << eq.getStiffnessOperator().format(
Eigen::IOFormat(Eigen::StreamPrecision, 0, ", ", "\n"));
file.close();

std::cout << "resolution\n";
Solver::LDLT(eq).solve();

Expand Down
211 changes: 41 additions & 170 deletions src/Rodin/Assembly/Multithreaded.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,36 +155,38 @@ namespace Rodin::Assembly
auto loop =
[&](const Index start, const Index end)
{
tl_lbfi.reset(bfi.copy());
tl_triplets.reserve(capacity / threadCount);
OperatorType triplets;
std::unique_ptr<LocalBilinearFormIntegratorBaseType> lbfi;
lbfi.reset(bfi.copy());
triplets.reserve(capacity / threadCount);
for (Index i = start; i < end; ++i)
{
if (seq.filter(i))
{
if (attrs.size() == 0 || attrs.count(mesh.getAttribute(d, i)))
{
const auto it = seq.getIterator(i);
tl_lbfi->setPolytope(*it);
lbfi->setPolytope(*it);
const auto& rows = input.getTestFES().getDOFs(d, i);
const auto& cols = input.getTrialFES().getDOFs(d, i);
for (size_t l = 0; l < static_cast<size_t>(rows.size()); l++)
{
for (size_t m = 0; m < static_cast<size_t>(cols.size()); m++)
{
const ScalarType s = tl_lbfi->integrate(m, l);
const ScalarType s = lbfi->integrate(m, l);
if (s != ScalarType(0))
tl_triplets.emplace_back(rows(l), cols(m), s);
triplets.emplace_back(rows(l), cols(m), s);
}
}
}
}
}
m_mutex.lock();
res.insert(res.end(),
std::make_move_iterator(tl_triplets.begin()),
std::make_move_iterator(tl_triplets.end()));
std::make_move_iterator(triplets.begin()),
std::make_move_iterator(triplets.end()));
m_mutex.unlock();
tl_triplets.clear();
triplets.clear();
};

if (std::holds_alternative<Threads::ThreadPool>(m_pool))
Expand All @@ -209,30 +211,32 @@ namespace Rodin::Assembly
auto loop =
[&](const Index start, const Index end)
{
tl_gbfi.reset(bfi.copy());
tl_triplets.reserve(capacity / threadCount);
OperatorType triplets;
std::unique_ptr<GlobalBilinearFormIntegratorBaseType> gbfi;
gbfi.reset(bfi.copy());
triplets.reserve(capacity / threadCount);
for (Index i = start; i < end; ++i)
{
if (testseq.filter(i))
{
if (testAttrs.size() == 0 || testAttrs.count(mesh.getAttribute(d, i)))
{
const auto teIt = testseq.getIterator(i);
Internal::SequentialIteration trialseq{ mesh, tl_gbfi->getTrialRegion() };
Internal::SequentialIteration trialseq{ mesh, gbfi->getTrialRegion() };
for (auto trIt = trialseq.getIterator(); trIt; ++trIt)
{
if (trialAttrs.size() == 0 || trialAttrs.count(trIt->getAttribute()))
{
tl_gbfi->setPolytope(*trIt, *teIt);
const auto& rows = input.getTestFES().getDOFs(d, i);
const auto& cols = input.getTrialFES().getDOFs(d, i);
gbfi->setPolytope(*trIt, *teIt);
const auto& rows = input.getTestFES().getDOFs(d, teIt->getIndex());
const auto& cols = input.getTrialFES().getDOFs(d, trIt->getIndex());
for (size_t l = 0; l < static_cast<size_t>(rows.size()); l++)
{
for (size_t m = 0; m < static_cast<size_t>(cols.size()); m++)
{
const ScalarType s = tl_gbfi->integrate(m, l);
const ScalarType s = gbfi->integrate(m, l);
if (s != ScalarType(0))
tl_triplets.emplace_back(rows(l), cols(m), s);
triplets.emplace_back(rows(l), cols(m), s);
}
}
}
Expand All @@ -242,10 +246,10 @@ namespace Rodin::Assembly
}
m_mutex.lock();
res.insert(res.end(),
std::make_move_iterator(tl_triplets.begin()),
std::make_move_iterator(tl_triplets.end()));
std::make_move_iterator(triplets.begin()),
std::make_move_iterator(triplets.end()));
m_mutex.unlock();
tl_triplets.clear();
triplets.clear();
};

if (std::holds_alternative<Threads::ThreadPool>(m_pool))
Expand Down Expand Up @@ -278,73 +282,10 @@ namespace Rodin::Assembly
}

private:
static thread_local OperatorType tl_triplets;
static thread_local std::unique_ptr<LocalBilinearFormIntegratorBaseType> tl_lbfi;
static thread_local std::unique_ptr<GlobalBilinearFormIntegratorBaseType> tl_gbfi;

mutable Threads::Mutex m_mutex;
mutable std::variant<Threads::ThreadPool, std::reference_wrapper<Threads::ThreadPool>> m_pool;
};

template <class TrialFES, class TestFES>
thread_local
std::vector<
Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>
Multithreaded<
std::vector<
Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>,
Variational::BilinearForm<TrialFES, TestFES,
std::vector<
Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>>>::tl_triplets;

template <class TrialFES, class TestFES>
thread_local
std::unique_ptr<
Variational::LocalBilinearFormIntegratorBase<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>
Multithreaded<
std::vector<
Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>,
Variational::BilinearForm<TrialFES, TestFES,
std::vector<
Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>>>::tl_lbfi;

template <class TrialFES, class TestFES>
thread_local
std::unique_ptr<
Variational::GlobalBilinearFormIntegratorBase<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>
Multithreaded<
std::vector<
Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>,
Variational::BilinearForm<TrialFES, TestFES,
std::vector<Eigen::Triplet<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>>>::tl_gbfi;

/**
* @brief Multithreaded assembly of the Math::SparseMatrix associated to a
* BilinearFormBase object.
Expand Down Expand Up @@ -541,27 +482,29 @@ namespace Rodin::Assembly
auto loop =
[&](const Index start, const Index end)
{
tl_lbfi.reset(bfi.copy());
tl_res.resize(input.getTestFES().getSize(), input.getTrialFES().getSize());
tl_res.setZero();
OperatorType pres;
std::unique_ptr<LocalBilinearFormIntegratorBaseType> lbfi;
lbfi.reset(bfi.copy());
pres.resize(input.getTestFES().getSize(), input.getTrialFES().getSize());
pres.setZero();
for (Index i = start; i < end; ++i)
{
if (seq.filter(i))
{
if (attrs.size() == 0 || attrs.count(mesh.getAttribute(d, i)))
{
const auto it = seq.getIterator(i);
tl_lbfi->setPolytope(*it);
lbfi->setPolytope(*it);
const auto& rows = input.getTestFES().getDOFs(d, i);
const auto& cols = input.getTrialFES().getDOFs(d, i);
for (size_t l = 0; l < static_cast<size_t>(rows.size()); l++)
for (size_t m = 0; m < static_cast<size_t>(cols.size()); m++)
tl_res(rows(l), cols(m)) += tl_lbfi->integrate(m, l);
pres(rows(l), cols(m)) += lbfi->integrate(m, l);
}
}
}
m_mutex.lock();
res += tl_res;
res += pres;
m_mutex.unlock();
};

Expand All @@ -587,7 +530,9 @@ namespace Rodin::Assembly
const auto loop =
[&](const Index start, const Index end)
{
tl_gbfi.reset(bfi.copy());
OperatorType tl_res;
std::unique_ptr<GlobalBilinearFormIntegratorBaseType> gbfi;
gbfi.reset(bfi.copy());
tl_res.resize(input.getTestFES().getSize(), input.getTrialFES().getSize());
tl_res.setZero();
for (Index i = start; i < end; ++i)
Expand All @@ -597,17 +542,17 @@ namespace Rodin::Assembly
if (testAttrs.size() == 0 || testAttrs.count(mesh.getAttribute(d, i)))
{
const auto teIt = testseq.getIterator(i);
Internal::SequentialIteration trialseq{ mesh, tl_gbfi->getTrialRegion() };
Internal::SequentialIteration trialseq{ mesh, gbfi->getTrialRegion() };
for (auto trIt = trialseq.getIterator(); trIt; ++trIt)
{
if (trialAttrs.size() == 0 || trialAttrs.count(trIt->getAttribute()))
{
tl_gbfi->setPolytope(*trIt, *teIt);
const auto& rows = input.getTestFES().getDOFs(d, i);
const auto& cols = input.getTrialFES().getDOFs(d, i);
gbfi->setPolytope(*trIt, *teIt);
const auto& rows = input.getTestFES().getDOFs(d, teIt->getIndex());
const auto& cols = input.getTrialFES().getDOFs(d, trIt->getIndex());
for (size_t l = 0; l < static_cast<size_t>(rows.size()); l++)
for (size_t m = 0; m < static_cast<size_t>(cols.size()); m++)
tl_res(rows(l), cols(m)) += tl_gbfi->integrate(m, l);
tl_res(rows(l), cols(m)) += gbfi->integrate(m, l);
}
}
}
Expand Down Expand Up @@ -648,67 +593,10 @@ namespace Rodin::Assembly
}

private:
static thread_local OperatorType tl_res;
static thread_local std::unique_ptr<LocalBilinearFormIntegratorBaseType> tl_lbfi;
static thread_local std::unique_ptr<GlobalBilinearFormIntegratorBaseType> tl_gbfi;

mutable Threads::Mutex m_mutex;
mutable std::variant<Threads::ThreadPool, std::reference_wrapper<Threads::ThreadPool>> m_pool;
};

template <class TrialFES, class TestFES>
thread_local
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>
Multithreaded<
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>,
Variational::BilinearForm<TrialFES, TestFES,
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>>::tl_res;

template <class TrialFES, class TestFES>
thread_local
std::unique_ptr<
Variational::LocalBilinearFormIntegratorBase<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>
Multithreaded<
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>,
Variational::BilinearForm<TrialFES, TestFES,
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>>::tl_lbfi;

template <class TrialFES, class TestFES>
thread_local
std::unique_ptr<
Variational::GlobalBilinearFormIntegratorBase<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>
Multithreaded<
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>,
Variational::BilinearForm<TrialFES, TestFES,
Math::Matrix<
typename FormLanguage::Dot<
typename FormLanguage::Traits<TrialFES>::ScalarType,
typename FormLanguage::Traits<TestFES>::ScalarType>::Type>>>::tl_gbfi;

/**
* @brief %Multithreaded assembly of the Math::Vector associated to a
* LinearForm object.
Expand Down Expand Up @@ -791,6 +679,8 @@ namespace Rodin::Assembly
const auto loop =
[&](const Index start, const Index end)
{
VectorType tl_res;
std::unique_ptr<Variational::LinearFormIntegratorBase<ScalarType>> tl_lfi;
tl_lfi.reset(lfi.copy());
tl_res.resize(input.getFES().getSize());
tl_res.setZero();
Expand Down Expand Up @@ -843,28 +733,9 @@ namespace Rodin::Assembly
}

private:
static thread_local VectorType tl_res;
static thread_local std::unique_ptr<Variational::LinearFormIntegratorBase<ScalarType>> tl_lfi;

mutable Threads::Mutex m_mutex;
mutable std::variant<Threads::ThreadPool, std::reference_wrapper<Threads::ThreadPool>> m_pool;
};

template <class FES>
thread_local
Math::Vector<typename FormLanguage::Traits<FES>::ScalarType>
Multithreaded<
Math::Vector<typename FormLanguage::Traits<FES>::ScalarType>,
Variational::LinearForm<FES, Math::Vector<typename FormLanguage::Traits<FES>::ScalarType>>>
::tl_res;

template <class FES>
thread_local
std::unique_ptr<Variational::LinearFormIntegratorBase<typename FormLanguage::Traits<FES>::ScalarType>>
Multithreaded<
Math::Vector<typename FormLanguage::Traits<FES>::ScalarType>,
Variational::LinearForm<FES, Math::Vector<typename FormLanguage::Traits<FES>::ScalarType>>>
::tl_lfi;
}

#endif
Expand Down
17 changes: 17 additions & 0 deletions src/Rodin/Variational/BilinearFormIntegrator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,6 +197,23 @@ namespace Rodin::Variational

using Parent::Parent;

template <class OtherNumber>
GlobalBilinearFormIntegratorBase(const GlobalBilinearFormIntegratorBase<OtherNumber>& other)
: Parent(other),
m_trialAttrs(other.m_trialAttrs),
m_testAttrs(other.m_testAttrs)
{}

/**
* @brief Move constructor.
*/
template <class OtherNumber>
GlobalBilinearFormIntegratorBase(GlobalBilinearFormIntegratorBase<OtherNumber>&& other)
: Parent(std::move(other)),
m_trialAttrs(std::move(other.m_trialAttrs)),
m_testAttrs(std::move(other.m_testAttrs))
{}

/**
* @brief Gets the attributes of the elements being integrated.
*/
Expand Down
Loading

0 comments on commit ad25629

Please sign in to comment.