Skip to content

Commit

Permalink
small refactoring. Mainly adding curve and loop index encoder
Browse files Browse the repository at this point in the history
  • Loading branch information
@rath3t
rath3t committed Feb 2, 2024
1 parent 4a01aa9 commit 36263ed
Show file tree
Hide file tree
Showing 9 changed files with 127 additions and 80 deletions.
2 changes: 1 addition & 1 deletion cmake/modules/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,6 @@
#
# SPDX-License-Identifier: LGPL-3.0-or-later

set(modules "DuneIgaMacros.cmake")
set(modules DuneIgaMacros.cmake AddAutoDiffFlags.cmake AddClipperLibFlags.cmake AddEarCutFlags.cmake AddEigenFlags.cmake AddMatplotppFlags.cmake AddnLohmannJsonFlags.cmake)

install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})
10 changes: 5 additions & 5 deletions dune/iga/geometrykernel/findintersection.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@

namespace Dune::IGANEW {

enum class FindIntersectionCurveAndLineResult { disjoint, intersect, parallel };
enum class IntersectionCurveAndLine { disjoint, intersect, parallel };

template <typename GeoCurve, typename ScalarType, int dim>
requires(dim == 2)
Expand All @@ -31,11 +31,11 @@ namespace Dune::IGANEW {

// Check if point is in domain
if (not geoCurve.domain()[0].checkInside(sol[0]))
return std::make_tuple(FindIntersectionCurveAndLineResult::disjoint, sol, curveP0);
return std::make_tuple(FindIntersectionCurveAndLineResult::intersect, sol, geoCurve.global(sol[0]));
return std::make_tuple(IntersectionCurveAndLine::disjoint, sol, curveP0);
return std::make_tuple(IntersectionCurveAndLine::intersect, sol, geoCurve.global(sol[0]));
}
// If system is not solvable, the lines are paralell and therefore have no intersection
return std::make_tuple(FindIntersectionCurveAndLineResult::parallel, tParameter, curveP0);
return std::make_tuple(IntersectionCurveAndLine::parallel, tParameter, curveP0);
}

template <typename GeoCurve, typename ScalarType, int dim>
Expand Down Expand Up @@ -93,7 +93,7 @@ namespace Dune::IGANEW {

if (not domain[0].checkInside(tParameter[0])) sucess = false;
return std::make_tuple(
sucess ? FindIntersectionCurveAndLineResult::intersect : FindIntersectionCurveAndLineResult::disjoint,
sucess ? IntersectionCurveAndLine::intersect : IntersectionCurveAndLine::disjoint,
tParameter, curvePoint);
}
} // namespace Dune::IGANEW
Expand Down
6 changes: 3 additions & 3 deletions dune/iga/test/testfindintersection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ auto test1() {
std::cout << std::setprecision(16) << std::endl;
auto [success, tParameter, curvePoint] = Dune::IGANEW::findIntersectionCurveAndLine(
curve1, Dune::FieldVector<double, 2>({0.5, 0.5}), {0.0, 5.0}, {0.5, 0.5});
t.check(success == FindIntersectionCurveAndLineResult::intersect) << "No intersection found";
t.check(success == IntersectionCurveAndLine::intersect) << "No intersection found";

t.check(Dune::FloatCmp::eq(tParameter[1], -0.3041666666666666))
<< "tLine is " << tParameter[1] << " but should be " << -0.3041666666666666;
Expand All @@ -72,7 +72,7 @@ auto test2() {
std::cout << std::setprecision(16) << std::endl;
auto [success, tParameter, curvePoint] = Dune::IGANEW::findIntersectionCurveAndLine(
curve2, Dune::FieldVector<double, 2>({0.5, 0.5}), {0.0, 5.0}, {0.5, 0.5});
t.check(success == FindIntersectionCurveAndLineResult::intersect) << "No intersection found";
t.check(success == IntersectionCurveAndLine::intersect) << "No intersection found";

t.check(Dune::FloatCmp::eq(tParameter[0], 0.4166666666666667))
<< "tCurve is " << tParameter[0] << " but should be " << 0.4166666666666667;
Expand All @@ -83,7 +83,7 @@ auto test2() {
// Parallel line
std::tie(success, tParameter, curvePoint)
= Dune::IGANEW::findIntersectionCurveAndLine(curve2, Dune::FieldVector<double, 2>({1, 1}), {1.2, -7}, {0.5, 0.5});
t.check(success == FindIntersectionCurveAndLineResult::parallel);
t.check(success == IntersectionCurveAndLine::parallel);

return t;
}
Expand Down
17 changes: 9 additions & 8 deletions dune/iga/test/trimtest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -338,7 +338,8 @@ struct ExpectedValues {
int notAffineCounter;
};

auto checkTrim(std::string filename, const ExpectedValues& expectedValues) {
template< typename ExecutionPolicy>
auto checkTrim(std::string filename, const ExpectedValues& expectedValues, ExecutionPolicy&& policy) {
// Dune::TestSuite t("", Dune::TestSuite::ThrowPolicy::ThrowOnRequired);

// Setup
Expand All @@ -355,7 +356,7 @@ auto checkTrim(std::string filename, const ExpectedValues& expectedValues) {
std::vector<int> yR(
range.begin(),
range.end()); // copy into vector sicne for_each does not work with iota_view and parallel execution
std::for_each(std::execution::par_unseq, yR.begin(), yR.end(), [&](auto refy) {
std::for_each(policy, yR.begin(), yR.end(), [&](auto refy) {
std::cout << "Thread: " << std::this_thread::get_id() << std::endl;
auto gridFactory = GridFactory();
auto brep = readJson<2>(filename);
Expand Down Expand Up @@ -425,7 +426,7 @@ auto checkTrim(std::string filename, const ExpectedValues& expectedValues) {

std::atomic<double> trimmedEdgeLengthsAccumulated{0};
std::for_each(
std::execution::par_unseq, gridView.template begin<0>(), gridView.template end<0>(), [&](const auto& ele) {
policy, gridView.template begin<0>(), gridView.template end<0>(), [&](const auto& ele) {
ElementTrimData elementTrimData = DefaultTrim::TrimmerImpl<2, 2, double>::trimElement(ele, patchTrimData);
auto [subTestEle, trimmedEdgeLength]
= elementTrimDataObstacleCourse(ele, elementTrimData, gridView, resTrimPatch);
Expand Down Expand Up @@ -465,31 +466,31 @@ int main(int argc, char** argv) try {
.trimmingCurveCurvedLength = 0.9765154639613303,
.trimmingCurveTotalLength = 3.61447516007258,
.firstLoopCurvesSize = 5,
.notAffineCounter = 1})));
.notAffineCounter = 1}),std::execution::seq));
t.subTest(checkTrim("auxiliaryfiles/element.ibra", ExpectedValues({.straightLength = 4,
.trimmingCurveCurvedLength = 0,
.trimmingCurveTotalLength = 4,
.firstLoopCurvesSize = 4,
.notAffineCounter = 0})));
.notAffineCounter = 0}),std::execution::seq));
t.subTest(
checkTrim("auxiliaryfiles/trim_2edges.ibra", ExpectedValues({.straightLength = 31.94725845850641,
.trimmingCurveCurvedLength = 6.323120188237797,
.trimmingCurveTotalLength = 38.27037864674421,
.firstLoopCurvesSize = 6,
.notAffineCounter = 2})));
.notAffineCounter = 2}),std::execution::seq));

t.subTest(
checkTrim("auxiliaryfiles/trim_multi.ibra", ExpectedValues({.straightLength = 27.06623257317474,
.trimmingCurveCurvedLength = 7.461084509983623,
.trimmingCurveTotalLength = 38.27037864674421,
.firstLoopCurvesSize = 5,
.notAffineCounter = 1})));
.notAffineCounter = 1}),std::execution::seq));
t.subTest(
checkTrim("auxiliaryfiles/element_trim_xb.ibra", ExpectedValues({.straightLength = 2.791977909806771,
.trimmingCurveCurvedLength = 1.6273832575204,
.trimmingCurveTotalLength = 4.419361167327171,
.firstLoopCurvesSize = 5,
.notAffineCounter = 1})));
.notAffineCounter = 1}),std::execution::seq));

return t.exit();
} catch (Dune::Exception& e) {
Expand Down
110 changes: 76 additions & 34 deletions dune/iga/trimmer/defaulttrimmer/patchtrimdata.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,8 @@

namespace Dune::IGANEW::DefaultTrim {

template <typename GridImp>
struct PatchTrimDataImpl {
using TrimmingCurve = typename GridImp::GridFamily::TrimmerTraits::TrimmingCurve;
using ParameterType = typename GridImp::GridFamily::Trimmer::ParameterType;
using ctype = typename GridImp::GridFamily::ctype;

namespace Impl {
template<typename TrimmingCurve>
struct BoundaryLoop {
void insertTrimCurve(const TrimmingCurve& curve) { curves_.push_back(curve); }

Expand All @@ -21,58 +17,99 @@ namespace Dune::IGANEW::DefaultTrim {
std::vector<TrimmingCurve> curves_;
};

template <typename ctype>
struct PointInPatch {
FieldVector<ctype, 2> pt;
size_t curveIdxI;
size_t curveIdxJ;
};



class CurveLoopIndexEncoder {
static constexpr int loopBits = 4;
static constexpr u_int64_t loopMask{(1ULL << loopBits) - 1};
public:

static int64_t encode(int curveIndex, int loopIndex) {
// the left hand side shifts the curve index value away from loop values and then inserts the curve index
// example curveBits = 4 and loopBits = 3;
// curveIndex = 3 and loopIndex = 2 = 00010_b
// curveIndex = 3 = 00011_b
// curveIndex<< loopBits = 11000_b
// curveIndex<< loopBits = 11000_b
// (curveIndex<< loopBits) | loopIndex = 11000_b | 00010_b = 11010_b = 26
return (static_cast<u_int64_t>(curveIndex) << loopBits) | loopIndex;
}

static std::pair<int, int> decode(int64_t singleIndex) {
// example curveBits = 4 and loopBits = 3;
// loopMask = 00111
// example singleIndex = 11010_b = 26
// curveIndex = 3 and loopIndex = 2 = 00010_b
// curveIndex = 3 = 00011_b
// loopMask =
// (singleIndex >> loopBits) = 00011_b = curveIndex
int curveIndex = (singleIndex >> loopBits);
// since the loop index is stored in the right-most bits, we only have to apply the bitwise and, which returns simply the loop index
// 11010_b AND
// 00111_b
// 00010_b = 2
int loopIndex = singleIndex & loopMask;
return {loopIndex,curveIndex};
}

};
}

template <typename GridImp>
struct PatchTrimDataImpl {
using TrimmingCurve = typename GridImp::GridFamily::TrimmerTraits::TrimmingCurve;
using ParameterType = typename GridImp::GridFamily::Trimmer::ParameterType;
using ctype = typename GridImp::GridFamily::ctype;



class CurveManager {
friend PatchTrimDataImpl;

public:
using idx_t = u_int64_t;
static constexpr idx_t indexOffSet = 4; //Since the first 4 values are reserved for the corners of the untrimmed square all indices start at 4

void addLoop(const BoundaryLoop& loop) {
void addLoop(const Impl::BoundaryLoop<TrimmingCurve>& loop) {
Clipper2Lib::PathD path;
for (idx_t i = loops_.empty() ? splitter_ : loops_.back().back().z + 1; const auto& curve : loop.curves()) {
if (curve.degree().front() == 1 and loops_.empty()) {
// since the first
for (int curveIndex=0; const auto& curve : loop.curves()) {
auto index= getZValue(curveIndex,loops_.size());
if (curve.affine() == 1 and loops_.empty()) {
auto p1 = curve.corner(0);
auto p2 = curve.corner(1);
path.emplace_back(p1[0], p1[1], i);
path.emplace_back(p2[0], p2[1], i + splitter_ - 1);
i += splitter_;
continue;
}
for (const auto v : Utilities::linspace(curve.domain()[0], splitter_)) {
auto fV = curve.global({v});
path.emplace_back(fV[0], fV[1], i++);
path.emplace_back(p1[0], p1[1], index);
path.emplace_back(p2[0], p2[1], index );
}else {
for (const auto v : Utilities::linspace(curve.domain()[0], splitter_)) {
auto fV = curve.global({v});
path.emplace_back(fV[0], fV[1], index);
}
}
++curveIndex;
}
loops_.push_back(path);
loopIndices_.emplace_back(loops_.back().back().z, loop.curves().size());
}

[[nodiscard]] auto getIndices(const idx_t val) const -> std::pair<size_t, size_t> {
size_t loopIdx = 0;
if (loops_.size() > 1) {
auto it = std::ranges::find_if(loopIndices_, [&](const std::pair<idx_t, size_t>& t) { return val > t.first; });
if (it == loopIndices_.end()) DUNE_THROW(Dune::IOError, "Invalid z-Value");
loopIdx = std::ranges::distance(loopIndices_.begin(), it) + 1;
}

auto curveIdx = static_cast<size_t>(std::floor(val / splitter_) - 1);
for (const auto& [_, sizeOfLoop] : std::ranges::take_view(loopIndices_, static_cast<long>(loopIdx))) {
curveIdx -= sizeOfLoop;
}
return Impl::CurveLoopIndexEncoder::decode(val-indexOffSet);
}

return std::make_pair(loopIdx, curveIdx);
[[nodiscard]] auto getZValue(int edgeIndex, int loopIndex )const {
return Impl::CurveLoopIndexEncoder::encode(edgeIndex, loopIndex)+indexOffSet;
}

private:
Clipper2Lib::PathsD loops_;
idx_t splitter_{};
std::vector<std::pair<idx_t, size_t>> loopIndices_{};

};

void addLoop() { loops_.push_back({}); }
Expand All @@ -91,6 +128,11 @@ namespace Dune::IGANEW::DefaultTrim {
auto getIndices(const typename CurveManager::idx_t val) const -> std::pair<size_t, size_t> {
return manager_.getIndices(val);
}

auto getZValue(int edgeIndex, int loopIndex ) const {
return manager_.getZValue(edgeIndex,loopIndex);
}

auto getCurve(const typename CurveManager::idx_t val) const -> const TrimmingCurve& {
auto [loopIdx, curveIdx] = manager_.getIndices(val);
return loops_[loopIdx].curves()[curveIdx];
Expand All @@ -100,7 +142,7 @@ namespace Dune::IGANEW::DefaultTrim {
}
auto getSplitter() const -> typename CurveManager::idx_t { return manager_.splitter_; }

auto getPointsInPatch(size_t loopIndex) const -> const std::vector<PointInPatch>& {
auto getPointsInPatch(size_t loopIndex) const -> const std::vector<Impl::PointInPatch<ctype>>& {
return pointsInPatch_.at(loopIndex);
}

Expand Down Expand Up @@ -133,9 +175,9 @@ namespace Dune::IGANEW::DefaultTrim {
}

private:
std::vector<std::vector<PointInPatch>> pointsInPatch_;
std::vector<std::vector<Impl::PointInPatch<ctype>>> pointsInPatch_;
bool finished_ = false;
std::vector<BoundaryLoop> loops_;
std::vector<Impl::BoundaryLoop<TrimmingCurve>> loops_;
CurveManager manager_;
};

Expand Down
8 changes: 4 additions & 4 deletions dune/iga/trimmer/defaulttrimmer/trimelement.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,10 @@ namespace Dune::IGANEW::DefaultTrim {
auto TrimmerImpl<dim, dimworld, ScalarType>::trimElement(const auto& element, const PatchTrimData& patchTrimData) {
// std::cout << "START " << std::endl;
auto geo = element.geometry();
std::array<FieldVector<double, 2>, 4> corners; // see dune book page 127 Figure 5.12
for (auto i : std::views::iota(0, geo.corners())) {
corners[i] = geo.corner(Util::vIdxMapping[i]);
}
static constexpr int numberOfCorners = 4;
std::array<FieldVector<double, 2>, numberOfCorners> corners; // see dune book page 127 Figure 5.12
for (auto i : Dune::range(numberOfCorners))
corners[i] = geo.corner(Util::vertexIndexMapping[i]);

auto [flag, result] = Util::clipElementRectangle(element, patchTrimData);

Expand Down
Loading

0 comments on commit 36263ed

Please sign in to comment.