Skip to content

Commit

Permalink
gridtest work again
Browse files Browse the repository at this point in the history
  • Loading branch information
@henrij22
henrij22 committed Apr 15, 2024
1 parent 56b7956 commit d1b378a
Show file tree
Hide file tree
Showing 6 changed files with 132 additions and 48 deletions.
42 changes: 37 additions & 5 deletions dune/iga/test/gridtests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,33 @@ auto checkUniqueEdges(const auto& gridView) {
}
}

template <typename GridView> requires (GridView::dimension == 2)
auto checkUniqueVertices(const GridView& gridView) {
Dune::TestSuite t;

constexpr int gridDimensionworld = GridView::dimensionworld;
constexpr int gridDimension = GridView::dimension;

std::set<std::array<FieldVector<double, gridDimensionworld>, 1>, Compare<double, gridDimensionworld, 1>>
elementVertexPairSet;
for (int eleIndex = 0; auto&& element : elements(gridView)) {
elementVertexPairSet.clear();
for (auto vertexIdx : Dune::range(element.subEntities(gridDimension))) {
auto vertex = element.template subEntity<gridDimension>(vertexIdx);

std::array<FieldVector<double, gridDimensionworld>, 1> pair;
for (auto c : Dune::range(vertex.geometry().corners()))
pair[c] = vertex.geometry().corner(c);

bool inserted = elementVertexPairSet.insert(pair).second;
t.require(inserted) << "Duplicate vertex detected in Element " << eleIndex << " Vertex: " << pair[0];
}
++eleIndex;
}
return t;

}

auto checkUniqueSurfaces(const auto& gridView) {
TestSuite t;

Expand Down Expand Up @@ -95,9 +122,14 @@ auto thoroughGridCheck(auto& grid) {
auto gvTest = [&](auto&& gv) {
TestSuite tl;

using GV = std::remove_cvref_t<decltype(gv)>;

tl.subTest(checkUniqueEdges(gv));
tl.subTest(checkUniqueSurfaces(gv));

if constexpr (GV::dimension == 2)
tl.subTest(checkUniqueVertices(gv));

auto extractGeo = std::views::transform([](const auto& ent) { return ent.geometry(); });
for (auto&& elegeo : elements(gv) | extractGeo)
checkJacobians(elegeo);
Expand All @@ -119,13 +151,13 @@ auto thoroughGridCheck(auto& grid) {
};

for (int lvl = 0; lvl <= grid.maxLevel(); ++lvl) {
auto gridView = grid.levelGridView(lvl);
t.subTest(gvTest(gridView));
t.subTest(gvTest(grid.levelGridView(lvl)));
}
t.subTest(gvTest(grid.leafGridView()));

gridcheck(grid);

try {
gridcheck(grid);
checkIntersectionIterator(grid);
checkLeafIntersections(grid);
} catch (const Dune::NotImplemented& e) {
Expand Down Expand Up @@ -865,7 +897,7 @@ auto testPlate() {
template <template <int, int, typename> typename GridFamily>
auto testGrids() {
TestSuite t("testGrids");
/*

// if constexpr (requires { testHierarchicPatch<GridFamily>(); }) {
std::cout << "testHierarchicPatch" << std::endl;
t.subTest(testHierarchicPatch<GridFamily>());
Expand Down Expand Up @@ -895,7 +927,7 @@ auto testGrids() {
// } else
// std::cout << "testNURBSGridSurface Test disabled" << std::endl;
// if constexpr (requires { testPlate<GridFamily>(); }) {
*/

std::cout << "testPlate==============================================" << std::endl;
t.subTest(testPlate<GridFamily>());
std::cout << "testPlateEND==============================================" << std::endl;
Expand Down
82 changes: 60 additions & 22 deletions dune/iga/trimmer/defaulttrimmer/createentities.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,20 @@ auto TrimmerImpl<dim, dimworld, ScalarType>::makeElementID(const HostEntity<0>&
}

template <int dim, int dimworld, typename ScalarType>
void TrimmerImpl<dim, dimworld, ScalarType>::createAndSaveElementInfo(const std::tuple<int, int, int>& indices,
auto TrimmerImpl<dim, dimworld, ScalarType>::idForTrimmedHostEdge(typename TrimmerTraits::PersistentIndexType hostEdgeId, const typename ElementTrimData::EdgeInfo& trimmedEdge) -> GlobalIdType {
auto& globalIdSetParameterSpace = parameterSpaceGrid_->globalIdSet();

// We have to check weather the edgeId is already there

GlobalIdType edgeId = {.entityIdType = GlobalIdType::EntityIdType::newId,
.id = globalIdSet_->newFreeIndex(),
.hostId = std::make_optional(hostEdgeId)};

return edgeId;
}

template <int dim, int dimworld, typename ScalarType>
void TrimmerImpl<dim, dimworld, ScalarType>::createAndSaveElementInfo(const std::tuple<unsigned int, unsigned int, int>& indices,
const HostEntity<0>& ele, bool trimmed) {
auto& globalIdSetParameterSpace = parameterSpaceGrid_->globalIdSet();

Expand All @@ -25,8 +38,8 @@ void TrimmerImpl<dim, dimworld, ScalarType>::createAndSaveElementInfo(const std:

GlobalIdType elementId = makeElementID(ele);

const int unTrimmedElementIndex = std::get<0>(indices);
const int trimmedElementIndex = std::get<1>(indices);
const unsigned int unTrimmedElementIndex = std::get<0>(indices);
const unsigned int trimmedElementIndex = std::get<1>(indices);
const int newLevel = std::get<2>(indices);
EntityInfo<0> elementInfo = {
.indexInLvlStorage = trimmedElementIndex + unTrimmedElementIndex,
Expand Down Expand Up @@ -68,14 +81,13 @@ void TrimmerImpl<dim, dimworld, ScalarType>::collectElementEdges(int level, cons
return;

auto edge = ele.template subEntity<1>(localEdgeIndex);
EntityInfo<1> edgeInfo{.indexInLvlStorage = indexSet.index(edge), .lvl = level, .stemFromTrim = false, .id = edgeId};
EntityInfo<1> edgeInfo{
.indexInLvlStorage = indexSet.index(edge), .lvl = level, .stemFromTrim = false, .id = edgeId, .hostSeed = edge.seed()};
entityContainer_.idToEdgeInfoMap.insert({edgeId, edgeInfo});

auto& edgeContainer = std::get<1>(entityContainer_.entityImps_.back());
edgeContainer.emplace_back(grid_, edge, edgeInfo);

// store vertex ids of edges, for subIndex method of indexSet
auto& edgeVertexIndices = entityContainer_.globalVertexIdOfEdgesMap_[edgeId];

if (edgeVertexIndices.size() < 2) {
const auto& cube = ReferenceElements<ctype, mydimension>::cube();
for (auto vertexLocalIndexWRTElement : cube.subEntities(localEdgeIndex, 1, 2)) {
Expand All @@ -89,13 +101,15 @@ void TrimmerImpl<dim, dimworld, ScalarType>::collectElementEdges(int level, cons
auto addTrimmedHostEdge = [&](int localEdgeIndex, const typename ElementTrimData::EdgeInfo& edgeOfTrimmedElement) {
// To get neighborhoud information we save the original hostId
auto hostEdgeId = globalIdSetParameterSpace.subId(ele, localEdgeIndex, 1);
GlobalIdType edgeId = idForTrimmedHostEdge(hostEdgeId, edgeOfTrimmedElement);

// \todo We have to somehow check if this edge is already present in the idSet
// Check for hostId, then for coordinates
// We have to get a unique .id

GlobalIdType edgeId = {.entityIdType = GlobalIdType::EntityIdType::newId,
.id = globalIdSet_->newFreeIndex(),
.hostId = std::make_optional(hostEdgeId)};
// GlobalIdType edgeId = {.entityIdType = GlobalIdType::EntityIdType::newId,
// .id = globalIdSet_->newFreeIndex(),
// .hostId = std::make_optional(hostEdgeId)};

};

Expand Down Expand Up @@ -128,26 +142,50 @@ void TrimmerImpl<dim, dimworld, ScalarType>::collectElementVertices(int level, c
auto& globalIdSetParameterSpace = parameterSpaceGrid_->globalIdSet();

if (eleTrimData.flag() == ElementTrimFlag::full) {
for (int localVertexId = 0; localVertexId < ele.subEntities(2); ++localVertexId) {
// setup all vertex indices for given element
auto hostVertexId = globalIdSetParameterSpace.subId(ele, localVertexId, 2);
IdType vertexId = {.entityIdType = GlobalIdType::EntityIdType::host, .id = hostVertexId};
// Save eleVertices in local order
for (auto localVertexIndex : Dune::range(ele.subEntities(2))) {
auto hostVertexId = globalIdSetParameterSpace.subId(ele, localVertexIndex, 2);
GlobalIdType vertexId = {.entityIdType = GlobalIdType::EntityIdType::host, .id = hostVertexId};
elementVertexIndices.emplace_back(vertexId);

if (entityContainer_.idToVertexInfoMap.back().contains(vertexId))
continue;

auto vertex = ele.template subEntity<2>(localVertexId);
auto vertex = ele.template subEntity<2>(localVertexIndex);
EntityInfo<2> vertexInfo{
.indexInLvlStorage = indexSet.index(vertex), .lvl = level, .stemFromTrim = false, .id = vertexId};
.indexInLvlStorage = indexSet.index(vertex), .lvl = level, .stemFromTrim = false, .id = vertexId, .hostSeed = vertex.seed()};
entityContainer_.idToVertexInfoMap.back().insert({vertexId, vertexInfo});

auto& vertexContainer = std::get<2>(entityContainer_.entityImps_.back());
vertexContainer.emplace_back(grid_, vertex, vertexInfo);
}

} else /* trimmed */ {

}
}

template <int dim, int dimworld, typename ScalarType>
void TrimmerImpl<dim, dimworld, ScalarType>::createSubEntities(int level) {
using EdgeEntity = typename TrimmerTraits::template Codim<1>::ParameterSpaceGridEntity;
using VertexEntity = typename TrimmerTraits::template Codim<2>::ParameterSpaceGridEntity;

auto& globalIdSetParameterSpace = parameterSpaceGrid_->globalIdSet();
auto& indexSet = parameterSpaceGrid_->levelGridView(level).indexSet();
auto& vertexMap = entityContainer_.idToVertexInfoMap.back();
auto& edgeMap = entityContainer_.idToEdgeInfoMap;

// we are resizing the containers, so we can add the entities in the indexSet order, the index is obtained from
// the infos (also this is a bit more efficient as pushing back all the time)

auto& vertexContainer = std::get<2>(entityContainer_.entityImps_.back());
vertexContainer.resize(vertexMap.size());
for (auto& [vertexId, vertexInfo] : vertexMap) {
auto vertex = parameterSpaceGrid_->entity(vertexInfo.hostSeed);
vertexContainer[vertexInfo.indexInLvlStorage] = VertexEntity{grid_, vertex, vertexInfo};
}

auto& edgeContainer = std::get<1>(entityContainer_.entityImps_.back());
edgeContainer.resize(entityContainer_.sizeOfInfos(1, level));
for (auto& [edgeId, edgeInfo] : edgeMap) {
if (edgeInfo.lvl != level)
continue;
auto edge = parameterSpaceGrid_->entity(edgeInfo.hostSeed);
edgeContainer[edgeInfo.indexInLvlStorage] = EdgeEntity{grid_, edge, edgeInfo};
}
}
} // namespace Dune::IGANEW::DefaultTrim
16 changes: 5 additions & 11 deletions dune/iga/trimmer/defaulttrimmer/createlevel.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,17 +10,13 @@ namespace Dune::IGANEW::DefaultTrim {
template <int dim, int dimworld, typename ScalarType>
void TrimmerImpl<dim, dimworld, ScalarType>::refineParameterSpaceGrid(int refCount, bool initFlag) {

using EdgeHostType = typename UntrimmedParameterSpaceGrid::template Codim<1>::Entity;
using EdgeGridType = typename GridImp::template Codim<1>::Entity;
using EleGridType = typename GridImp::template Codim<0>::Entity;
using EdgeParameterSpaceType = typename TrimmerTraits::template Codim<1>::ParameterSpaceGridEntity;
using EleParameterSpaceType = typename TrimmerTraits::template Codim<2>::ParameterSpaceGridEntity;
const int oldLevel = untrimmedParameterSpaceGrid_->maxLevel();
untrimmedParameterSpaceGrid_->globalRefine(refCount);
auto gvu = untrimmedParameterSpaceGrid_->leafGridView();
parameterSpaceGrid_->createBegin();
parameterSpaceGrid_->insertLeaf();
parameterSpaceGrid_->createEnd();

assert((initFlag and oldLevel == 0 and refCount == 0) or
!initFlag && "If we initialize the grid, the untrimmedParameterSpaceGrid_ should only have one level");
// if we init we start at 0 otherwise at 1
Expand All @@ -29,21 +25,17 @@ void TrimmerImpl<dim, dimworld, ScalarType>::refineParameterSpaceGrid(int refCou
const int newLevel = oldLevel + i;
entityContainer_.entityImps_.emplace_back();
auto& entityContainer = entityContainer_;
// auto& elementContainer = std::get<0>(entityContainer_.entityImps_.back());
// auto& edgeContainer = std::get<1>(entityContainer_.entityImps_.back());
// auto& vertexContainer = std::get<2>(entityContainer_.entityImps_.back());

auto gv = parameterSpaceGrid_->levelGridView(newLevel);
auto& globalIdSetParameterSpace = parameterSpaceGrid_->globalIdSet();
int unTrimmedElementIndex = 0;
int trimmedElementIndex = 0;
unsigned int unTrimmedElementIndex = 0;
unsigned int trimmedElementIndex = 0;

auto indices = [&]() { return std::make_tuple(unTrimmedElementIndex, trimmedElementIndex, newLevel); };

auto& indexSet = gv.indexSet();
const auto elementTrimDatas = trimElements(newLevel);

// Add a new level for vertices
entityContainer.idToVertexInfoMap.emplace_back();

for (const auto& ele : elements(gv)) {
Expand Down Expand Up @@ -76,6 +68,8 @@ void TrimmerImpl<dim, dimworld, ScalarType>::refineParameterSpaceGrid(int refCou
entityContainer.numberOfTrimmedElements.push_back(trimmedElementIndex);
entityContainer.numberOfUnTrimmedElements.push_back(unTrimmedElementIndex);

// Add a new level for vertices
createSubEntities(newLevel);
}
}
} // namespace Dune::IGANEW::DefaultTrim
17 changes: 17 additions & 0 deletions dune/iga/trimmer/defaulttrimmer/entitycontainer.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,23 @@ struct VectorEntityContainer
return 0;
}


std::size_t sizeOfInfos(int codim, int lvl) const {
if (codim == 2)
return idToVertexInfoMap[lvl].size();

auto count = [=](auto&& range) -> std::size_t {
return std::ranges::count_if(range, [=](const auto& entityInfoPair) {
return entityInfoPair.second.lvl == lvl;
});
};
if (codim == 0)
return count(idToElementInfoMap);
if (codim == 1)
return count(idToEdgeInfoMap);
__builtin_unreachable();
}

// The vector type for tuple of lists of vertices, edges, elements for each level
using EntityImps = std::decay_t<decltype(makeEntityImps_())>;

Expand Down
6 changes: 3 additions & 3 deletions dune/iga/trimmer/defaulttrimmer/patchgridentityseed.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,13 +55,13 @@ public:
/**
* @brief Check whether it is safe to create an Entity from this Seed
*/
bool isValid() const { return indexInLvlStorage_ != -1; }
bool isValid() const { return indexInLvlStorage_ != std::numeric_limits<unsigned int>::max(); }

private:
auto data() const { return std::make_pair(lvl_, indexInLvlStorage_); }

int lvl_;
int indexInLvlStorage_{-1};
int lvl_{};
unsigned int indexInLvlStorage_{std::numeric_limits<unsigned int>::max()};
};

} // namespace Dune::IGANEW::DefaultTrim
Expand Down
17 changes: 10 additions & 7 deletions dune/iga/trimmer/defaulttrimmer/trimmer.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ namespace DefaultTrim {
using HostIdType = typename GridImpl::GlobalIdSet::IdType;
using EntitySeedType = typename GridImpl::template Codim<codim>::Entity::EntitySeed;

int indexInLvlStorage{};
unsigned int indexInLvlStorage{};
int lvl{};
bool stemFromTrim{false};
IdType<HostIdType> id;
Expand All @@ -158,9 +158,9 @@ namespace DefaultTrim {
using HostIdType = typename GridImpl::GlobalIdSet::IdType;
using EntitySeedType = typename GridImpl::template Codim<0>::Entity::EntitySeed;

int indexInLvlStorage{-1};
int unTrimmedIndexInLvl{-1};
int trimmedIndexInLvl{-1};
unsigned int indexInLvlStorage{std::numeric_limits<unsigned int>::infinity()};
unsigned int unTrimmedIndexInLvl{std::numeric_limits<unsigned int>::infinity()};
unsigned int trimmedIndexInLvl{std::numeric_limits<unsigned int>::infinity()};
int lvl{};
bool stemFromTrim{false};
IdType<HostIdType> id{};
Expand Down Expand Up @@ -587,10 +587,13 @@ namespace DefaultTrim {


GlobalIdType makeElementID(const HostEntity<0>& ele);
void createAndSaveElementInfo(const std::tuple<int, int, int>& indices, const HostEntity<0>& ele, bool trimmed);
void createAndSaveElementInfo(const std::tuple<unsigned int, unsigned int, int>& indices, const HostEntity<0>& ele, bool trimmed);

void collectElementEdges(int level, const HostEntity<0>& ele, const ElementTrimData& trimData);
void collectElementVertices(int level, const HostEntity<0>& ele, const ElementTrimData& trimData);
void collectElementEdges(int level, const HostEntity<0>& ele, const ElementTrimData& eleTrimData);
void collectElementVertices(int level, const HostEntity<0>& ele, const ElementTrimData& eleTrimData);
void createSubEntities(int level);

GlobalIdType idForTrimmedHostEdge(typename TrimmerTraits::PersistentIndexType hostEdgeId, const typename ElementTrimData::EdgeInfo& trimmedEdge);

/** @brief Return maximum level defined in this grid.
*
Expand Down

0 comments on commit d1b378a

Please sign in to comment.