diff --git a/dune/iga/geometrykernel/nurbspatchgeometry.hh b/dune/iga/geometrykernel/nurbspatchgeometry.hh index 0e7679f3..86ce2ee3 100644 --- a/dune/iga/geometrykernel/nurbspatchgeometry.hh +++ b/dune/iga/geometrykernel/nurbspatchgeometry.hh @@ -317,6 +317,102 @@ namespace Dune::IGANEW::GeometryKernel { /* @brief Get the patch data of the NURBS patch. */ auto& patchData() { return patchData_; } + struct ConnectionResult { + bool isConnected; + int boundary; + operator bool() const { return isConnected; } + }; + + /** + * @brief Check if the patch is connected at the given boundary. + * + * @details This function checks if the patch is connected at the specified boundary. The boundary index is + * determined by the vertex numbering for lines, the edge numbering for quadrilaterals, and the face numbering for + * cubes. For more details, refer to DUNE Book Chapter 5.5. + * + * @param boundary The index of the boundary to check (vertex for lines, edge for quadrilaterals, face for cubes). + * @param tol Absolute tolerance for comparing control points. Default is 1e-8. + * + * @return A ConnectionResult object indicating whether the patch is connected at the specified boundary. + * The isConnected member is true if the patch is connected, and the boundary member is the index of the + * adjacent boundary if connected, or -1 if not connected. + * + * @tparam boundary The boundary index to check. + * @tparam tol Absolute tolerance for comparing control points. Default is 1e-8. + * + * @throws Dune::NotImplemented If the dimension is greater than 3 (higher dimensions not implemented yet). + */ + ConnectionResult isConnectedAtBoundary(int boundary, double tol = 1e-8) const { + const auto& cps = patchData_.controlPoints; + if constexpr (mydimension == 1) { + if (auto cp = cps.directGetAll(); + FloatCmp::eq, FloatCmp::CmpStyle::absolute>(cp.front().p, cp.back().p, + tol)) + return ConnectionResult(true, !boundary); + else + return ConnectionResult(false, -1); + } else if constexpr (mydimension == 2) { + auto edgeRange = [&cps](int b) { + switch (b) { + case 0: + return cps.hyperSurfFront({1}); + case 1: + return cps.hyperSurfBack({1}); + case 2: + return cps.hyperSurfFront({0}); + case 3: + return cps.hyperSurfBack({0}); + default: + assert("No valid boundary index"); + } + __builtin_unreachable(); + }; + + auto theEdgeToCheck = edgeRange(boundary); + for (auto i : Dune::range(4)) { + if (boundary == i) continue; + if (std::ranges::equal(theEdgeToCheck, edgeRange(i), [tol](auto f, auto s) { + return FloatCmp::eq, FloatCmp::CmpStyle::absolute>(f.p, s.p, tol); + })) + return ConnectionResult(true, i); + } + return ConnectionResult(false, -1); + } else if constexpr (mydimension == 3) { + auto edgeRange = [&cps](int b) { + switch (b) { + case 0: + return cps.hyperSurfFront({1, 2}); + case 1: + return cps.hyperSurfBack({1, 2}); + case 2: + return cps.hyperSurfFront({0, 1}); + case 3: + return cps.hyperSurfBack({0, 1}); + case 4: + return cps.hyperSurfFront({0, 2}); + case 5: + return cps.hyperSurfBack({0, 2}); + default: + assert("No valid boundary index"); + } + __builtin_unreachable(); + }; + + auto theEdgeToCheck = edgeRange(boundary); + for (auto i : Dune::range(6)) { + if (boundary == i) continue; + if (std::ranges::equal(theEdgeToCheck, edgeRange(i), [tol](auto f, auto s) { + return FloatCmp::eq, FloatCmp::CmpStyle::absolute>(f.p, s.p, tol); + })) + return ConnectionResult(true, i); + } + return ConnectionResult(false, -1); + } else + DUNE_THROW(Dune::NotImplemented, + "Higher dimensions not implemented yet. But should be straight-forward to do so, somebody just has " + "to come up with a resenable indexing of hyper surfaces."); + } + friend auto referenceElement(const NURBSPatch& geo) { return referenceElement(Dune::GeometryTypes::cube(mydimension)); } diff --git a/dune/iga/hierarchicpatch/patchgrid.hh b/dune/iga/hierarchicpatch/patchgrid.hh index f267920c..aba30fbe 100644 --- a/dune/iga/hierarchicpatch/patchgrid.hh +++ b/dune/iga/hierarchicpatch/patchgrid.hh @@ -456,6 +456,7 @@ namespace Dune::IGANEW { auto untrimmedElementNumbers(int lvl) const { return patchGeometries_[lvl].numberOfSpans(); } const auto& trimmer() const { return *trimmer_; } + const auto& patchGeometry(int i) const { return patchGeometries_.at(i); } private: PatchGrid() = default; diff --git a/dune/iga/test/gridtests.cpp b/dune/iga/test/gridtests.cpp index b1e6d149..cf343a24 100644 --- a/dune/iga/test/gridtests.cpp +++ b/dune/iga/test/gridtests.cpp @@ -184,6 +184,21 @@ auto testNurbsGridCylinder() { TestSuite testSuite; + auto patchGeo = grid.patchGeometry(0); + + // Since this is only a quarter of a cylinder no glueing should be present + auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(0); + testSuite.check(!glued) << "The edge " << 0 << " should be glued!"; + testSuite.check(edgeIndex == -1); + auto [glued2, edgeIndex2] = patchGeo.isConnectedAtBoundary(0); + testSuite.check(!glued2) << "The edge " << 2 << " shouldn't be glued!"; + testSuite.check(edgeIndex == -1); + + Dune::RefinementIntervals refinementIntervals1(3); + SubsamplingVTKWriter vtkWriter(grid.leafGridView(), refinementIntervals1); + // vtkWriter.write("NURBSGridTest-CurveNewFineResample"); + vtkWriter.write("NURBSGridTest-Zylinder"); + testSuite.subTest(thoroughGridCheck(grid)); return testSuite; } @@ -206,6 +221,33 @@ auto testHierarchicPatch() { Dune::IGANEW::PatchGrid<2, 3, GridFamily> patch(nurbsPatchData); + auto glueIndicator = [](int i) { + switch (i) { + case 0: + return 1; + case 1: + return 0; + case 2: + return 3; + case 3: + return 2; + default: + DUNE_THROW(RangeError, "Invalid index"); + } + }; + auto patchGeo = patch.patchGeometry(0); + for (auto i : Dune::range(4)) { + auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(i); + t.check(glued) << "The edge " << i << " should be glued!"; + t.check(edgeIndex == glueIndicator(i)) + << "The edge " << i << " should be glued to edge " << glueIndicator(i) << " but is glued to " << edgeIndex; + } + + Dune::RefinementIntervals refinementIntervals1(3); + SubsamplingVTKWriter vtkWriter(patch.leafGridView(), refinementIntervals1); + // vtkWriter.write("NURBSGridTest-CurveNewFineResample"); + vtkWriter.write("NURBSGridTest-Torus"); + t.subTest(thoroughGridCheck(patch)); return t; @@ -326,6 +368,13 @@ auto testNURBSGridCurve() { TestSuite t; IGANEW::PatchGrid grid(patchData); + + auto patchGeo = grid.patchGeometry(0); + + auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(0); + t.check(glued) << "The edge " << 0 << " should be glued!"; + t.check(edgeIndex == -1); + grid.globalRefine(3); auto gridView = grid.leafGridView(); const auto& indexSet = gridView.indexSet(); @@ -389,6 +438,13 @@ auto testNURBSGridSurface() { auto gridView = grid.leafGridView(); const auto& indexSet = gridView.indexSet(); + auto patchGeo = grid.patchGeometry(0); + for (auto i : Dune::range(4)) { + auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(i); + t.check(!glued) << "The edge " << i << " shouldn't be glued!"; + t.check(edgeIndex == -1); + } + Dune::RefinementIntervals refinementIntervals1(subSampling); SubsamplingVTKWriter vtkWriter(gridView, refinementIntervals1); vtkWriter.write("NURBSGridTest-Surface"); @@ -451,6 +507,14 @@ auto test3DGrid() { auto gridView = grid.leafGridView(); TestSuite t; + auto patchGeo = grid.patchGeometry(0); + + for (auto i : Dune::range(6)) { + auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(i); + t.check(!glued) << "The surface " << i << " shouldn't be glued!"; + t.check(edgeIndex == -1); + } + t.subTest(checkUniqueEdges(gridView)); t.subTest(checkUniqueSurfaces(gridView)); diff --git a/dune/iga/trimmer/defaulttrimmer/trimelement.hh b/dune/iga/trimmer/defaulttrimmer/trimelement.hh index e3c455f1..df2c8229 100644 --- a/dune/iga/trimmer/defaulttrimmer/trimelement.hh +++ b/dune/iga/trimmer/defaulttrimmer/trimelement.hh @@ -98,32 +98,31 @@ namespace Dune::IGANEW::DefaultTrim { = Util::callFindIntersection(patchTrimData.getCurve(currentCurveIdx), getEdgeIdx(vV2), pt2, corners); assertPoint(pt2, curvePoint); - // If currentT > tParam it can have 2 reasons. 1) the trim id only on one side of the element, and it has to connect back - // or 2) if the trimming curve consist only of one curve where it has the same front and back Controlpoint, so sometimes - // the wrong tParam (front or back) gets determined + // If currentT > tParam it can have 2 reasons. 1) the trim id only on one side of the element, and it has to + // connect back or 2) if the trimming curve consist only of one curve where it has the same front and back + // Controlpoint, so sometimes the wrong tParam (front or back) gets determined if (currentT > tParam) { bool success = false; if (currentCurveIdx.first > 0 && patchTrimData.loops()[currentCurveIdx.first].size() == 1) { auto curve = patchTrimData.getCurve(currentCurveIdx); - if (auto cp = curve.patchData().controlPoints.directGetAll(); FloatCmp::eq(cp.front().p, cp.back().p)) { + if (curve.isConnectedAtBoundary(0)) { auto elementTrimmingCurve = Util::createTrimmingCurveSlice(curve, currentT, curve.domain()[0].back()); elementTrimData.addEdgeNewNew(elementTrimmingCurve, curvePoint); - currentT = std::numeric_limits::infinity(); - success = true; + currentT = std::numeric_limits::infinity(); + success = true; } } else if (getEdgeIdx(vV1) == getEdgeIdx(vV2)) { auto trimmedEdge = Util::createHostGeometry(foundVertices.back(), curvePoint); elementTrimData.addEdgeNewNewOnHost(getEdgeIdx(vV2), trimmedEdge, curvePoint); - currentT = tParam; - success = true; + currentT = tParam; + success = true; } - if (not success) - throwGridError(); + if (not success) throwGridError(); } else { auto elementTrimmingCurve = Util::createTrimmingCurveSlice(patchTrimData.getCurve(currentCurveIdx), currentT, tParam); elementTrimData.addEdgeNewNew(elementTrimmingCurve, curvePoint); - currentT = std::numeric_limits::infinity(); + currentT = std::numeric_limits::infinity(); } foundVertices.push_back(curvePoint); diff --git a/dune/iga/utils/mdnet.hh b/dune/iga/utils/mdnet.hh index 21dc7d68..35f9432b 100644 --- a/dune/iga/utils/mdnet.hh +++ b/dune/iga/utils/mdnet.hh @@ -16,13 +16,9 @@ namespace Dune::IGANEW { namespace Impl { - template + template class HyperSurfaceIterator; - template - HyperSurfaceIterator operator+(const HyperSurfaceIterator& l, int inc); - template - HyperSurfaceIterator operator-(const HyperSurfaceIterator& l, int inc); } // namespace Impl /** @@ -366,8 +362,13 @@ namespace Dune::IGANEW { Impl::HyperSurfaceIterator hyperSurfBegin( const std::array& direction) { - return Impl::HyperSurfaceIterator(*this, direction, 0); + return Impl::HyperSurfaceIterator(*this, direction, 0); } + + auto hyperSurfFront(const std::array& direction) { + return *this->hyperSurfBegin(direction); + } + Impl::HyperSurfaceIterator hyperSurfEnd( const std::array& direction) { int directionEnd; @@ -380,7 +381,58 @@ namespace Dune::IGANEW { } else directionEnd = this->strideSizes()[0]; - return Impl::HyperSurfaceIterator(*this, direction, directionEnd); + return Impl::HyperSurfaceIterator(*this, direction, directionEnd); + } + + auto hyperSurfBack(const std::array& direction) { + int directionEnd; + if constexpr (netdim != 0) { + for (int dirI = 0, i = 0; i < netdim; ++i) { + if (dirI < direction.size() && i == direction.at(dirI++)) continue; + directionEnd = this->strideSizes()[i] - 1; + break; + } + } else + directionEnd = this->strideSizes()[0] - 1; + + return *Impl::HyperSurfaceIterator(*this, direction, directionEnd); + } + + Impl::HyperSurfaceIterator hyperSurfBegin( + const std::array& direction) const { + return Impl::HyperSurfaceIterator(*this, direction, 0); + } + + auto hyperSurfFront(const std::array& direction) const { + return *this->hyperSurfBegin(direction); + } + Impl::HyperSurfaceIterator hyperSurfEnd( + const std::array& direction) const { + int directionEnd; + if constexpr (netdim != 0) { + for (int dirI = 0, i = 0; i < netdim; ++i) { + if (dirI < direction.size() && i == direction.at(dirI++)) continue; + directionEnd = this->strideSizes()[i]; + break; + } + } else + directionEnd = this->strideSizes()[0]; + + return Impl::HyperSurfaceIterator(*this, direction, directionEnd); + } + + auto hyperSurfBack(const std::array& direction) const { + int directionBack; + if constexpr (netdim != 0) { + for (int dirI = 0, i = 0; i < netdim; ++i) { + if (dirI < direction.size() && i == direction.at(dirI++)) continue; + directionBack = this->strideSizes()[i] - 1; + break; + } + } else + directionBack = this->strideSizes()[0] - 1; + + return *Impl::HyperSurfaceIterator(*this, direction, directionBack); } private: @@ -393,11 +445,12 @@ namespace Dune::IGANEW { -> MultiDimensionalNet; namespace Impl { - template + template class HyperSurfaceIterator { public: - HyperSurfaceIterator(MultiDimensionalNet& net, const std::array& direction, - int at) + using MDNetType = std::conditional_t, + MultiDimensionalNet>; + HyperSurfaceIterator(MDNetType& net, const std::array& direction, int at) : net_{&net}, direction_{direction}, at_{at} { std::array indicesSurface; for (int i = 0; i < indicesSurface.size(); ++i) @@ -454,7 +507,7 @@ namespace Dune::IGANEW { else multiIndex[0] = at_; - auto objectExtractor = [ multiIndex, this ](auto mI) mutable -> auto& { + auto objectExtractor = [ multiIndex, *this ](auto mI) mutable -> auto& { if constexpr (netdim != 1) for (int i = 0; i < direction_.size(); ++i) multiIndex[direction_[i]] = mI[i]; @@ -466,25 +519,18 @@ namespace Dune::IGANEW { auto operator->() { return *this; } + friend HyperSurfaceIterator operator+(const HyperSurfaceIterator& l, int inc) { + return HyperSurfaceIterator(*(l.net_), l.direction_, l.at_ + inc); + } + + friend HyperSurfaceIterator operator-(const HyperSurfaceIterator& l, int inc) { return l + (-inc); } + private: - MultiDimensionalNet* net_; + MDNetType* net_; std::array direction_; int at_; - template - friend HyperSurfaceIterator operator+(const HyperSurfaceIterator& l, - int inc); - template - friend HyperSurfaceIterator operator-(const HyperSurfaceIterator& l, - int inc); }; - template - HyperSurfaceIterator operator+(const HyperSurfaceIterator& l, const int inc) { - return HyperSurfaceIterator(*(l.net_), l.direction_, l.at_ + inc); - } - template - HyperSurfaceIterator operator-(const HyperSurfaceIterator& l, const int inc) { - return l + (-inc); - } + } // namespace Impl template