Skip to content

Boilerplate for different element/vertex size #87

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
vikaskurapati merged 2 commits into from
Mar 17, 2025
Merged

Boilerplate for different element/vertex size #87

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c998913
Basic different dimension/element size support
davschneller Mar 13, 2025
621da3f
Add even more higher orders
davschneller Mar 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 12 additions & 8 deletions src/input/MeshData.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,21 +19,25 @@ class MeshData {
public:
virtual ~MeshData() = default;

virtual std::size_t cellCount() = 0;
virtual std::size_t vertexCount() = 0;
virtual std::size_t cellCount() const = 0;
virtual std::size_t vertexCount() const = 0;

virtual const std::vector<uint64_t>& connectivity() = 0;
virtual const std::vector<double>& geometry() = 0;
virtual const std::vector<int>& group() = 0;
virtual const std::vector<int64_t>& boundary() = 0;

// TODO: generalize?
virtual std::size_t vertexSize() const { return 3; }
virtual std::size_t cellSize() const { return 4; }
};

class FullStorageMeshData : public MeshData {
public:
virtual ~FullStorageMeshData() = default;

virtual std::size_t cellCount() { return cellCountValue; }
virtual std::size_t vertexCount() { return vertexCountValue; }
virtual std::size_t cellCount() const { return cellCountValue; }
virtual std::size_t vertexCount() const { return vertexCountValue; }

virtual const std::vector<uint64_t>& connectivity() { return connectivityData; }
virtual const std::vector<double>& geometry() { return geometryData; }
Expand All @@ -55,7 +59,7 @@ class FullStorageMeshData : public MeshData {

void setBoundary(std::size_t cell, int face, int value) {
if (bndShift < 0) {
boundaryData[cell * 4 + face] = value;
boundaryData[cell * (vertexSize() + 1) + face] = value;
} else {
constexpr auto i32limit = std::numeric_limits<unsigned char>::max();
constexpr auto i64limit = std::numeric_limits<unsigned short>::max();
Expand All @@ -72,11 +76,11 @@ class FullStorageMeshData : public MeshData {
cellCountValue = cellCount;
vertexCountValue = vertexCount;

connectivityData.resize(cellCount * 4);
geometryData.resize(vertexCount * 3);
connectivityData.resize(cellCount * cellSize());
geometryData.resize(vertexCount * vertexSize());
groupData.resize(cellCount);
if (bndShift < 0) {
boundaryData.resize(cellCount * 4);
boundaryData.resize(cellCount * (vertexSize() + 1));
} else {
boundaryData.resize(cellCount);
}
Expand Down
9 changes: 6 additions & 3 deletions src/input/SerialMeshFile.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ template <typename T> class SerialMeshFile : public FullStorageMeshData {
public:
virtual ~SerialMeshFile() = default;

std::size_t vertexSize() const override { return T::Dim; }
std::size_t cellSize() const override { return puml::nodeCount(T::Dim, T::Order); }

private:
#ifdef PARALLEL
MPI_Comm m_comm;
Expand Down Expand Up @@ -89,11 +92,11 @@ template <typename T> class SerialMeshFile : public FullStorageMeshData {
m_meshReader.readGroups(groupData.data());

logInfo(m_rank) << "Read boundary conditions";
std::vector<int> preBoundaryData(nLocalElements * 4);
std::vector<int> preBoundaryData(nLocalElements * (vertexSize() + 1));
m_meshReader.readBoundaries(preBoundaryData.data());
for (std::size_t i = 0; i < nLocalElements; ++i) {
for (int j = 0; j < 4; ++j) {
setBoundary(i, j, preBoundaryData[4 * i + j]);
for (int j = 0; j < (vertexSize() + 1); ++j) {
setBoundary(i, j, preBoundaryData[(vertexSize() + 1) * i + j]);
}
}
}
Expand Down
150 changes: 130 additions & 20 deletions src/meshreader/GMSHBuilder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,37 +1,147 @@
#ifndef GMSHBUILDER_20201014_H
#define GMSHBUILDER_20201014_H

#include <algorithm>
#include <array>
#include <cassert>
#include <vector>

#include "utils/logger.h"

#include "third_party/GMSHMeshBuilder.h"

namespace puml {

template <std::size_t D> struct GMSHSimplexType {};
constexpr std::size_t nodeCount(std::size_t D, std::size_t Order) {
std::size_t count = 1;
for (std::size_t d = 1; d <= D; ++d) {
count *= (Order + D - d + 1);
count /= d;
}
return count;
}

template <std::size_t P, typename BaseT, template <std::size_t> typename ParamT, typename... Args>
BaseT* makePointerI(std::size_t param, Args&&... args) {
if constexpr (P < 11) {
if (param == P) {
return new ParamT<P>(std::forward<Args>(args)...);
} else {
return makePointerI<P + 1, BaseT, ParamT>(param, std::forward<Args>(args)...);
}
} else {
logError() << "Order too high:" << param;
return nullptr;
}
}

template <typename BaseT, template <std::size_t> typename ParamT, typename... Args>
BaseT* makePointer(std::size_t param, Args&&... args) {
return makePointerI<1, BaseT, ParamT>(param, std::forward<Args>(args)...);
}

template <std::size_t D, std::size_t Order> struct GMSHSimplexType {};

// clang format seems to mess up the following lines from time to time, hence we disable it
// clang-format off
template <> struct GMSHSimplexType<0u> {
template <std::size_t Order> struct GMSHSimplexType<0U, Order> {
static constexpr long type = 15;
}; // 15 -> point
template <> struct GMSHSimplexType<1u> {
};
template <> struct GMSHSimplexType<1U, 1U> {
static constexpr long type = 1;
}; // 1 -> line
template <> struct GMSHSimplexType<2u> {
};
template <> struct GMSHSimplexType<2U, 1U> {
static constexpr long type = 2;
}; // 2 -> triangle
template <> struct GMSHSimplexType<3u> {
};
template <> struct GMSHSimplexType<3U, 1U> {
static constexpr long type = 4;
}; // 4 -> tetrahedron
// clang-format on
};
template <> struct GMSHSimplexType<1U, 2U> {
static constexpr long type = 8;
};
template <> struct GMSHSimplexType<2U, 2U> {
static constexpr long type = 9;
};
template <> struct GMSHSimplexType<3U, 2U> {
static constexpr long type = 11;
};
template <> struct GMSHSimplexType<1U, 3U> {
static constexpr long type = 26;
};
template <> struct GMSHSimplexType<1U, 4U> {
static constexpr long type = 27;
};
template <> struct GMSHSimplexType<1U, 5U> {
static constexpr long type = 28;
};
template <> struct GMSHSimplexType<1U, 6U> {
static constexpr long type = 62;
};
template <> struct GMSHSimplexType<1U, 7U> {
static constexpr long type = 63;
};
template <> struct GMSHSimplexType<1U, 8U> {
static constexpr long type = 64;
};
template <> struct GMSHSimplexType<1U, 9U> {
static constexpr long type = 65;
};
template <> struct GMSHSimplexType<1U, 10U> {
static constexpr long type = 66;
};
template <> struct GMSHSimplexType<2U, 3U> {
static constexpr long type = 21;
};
template <> struct GMSHSimplexType<2U, 4U> {
static constexpr long type = 23;
};
template <> struct GMSHSimplexType<2U, 5U> {
static constexpr long type = 25;
};
template <> struct GMSHSimplexType<2U, 6U> {
static constexpr long type = 42;
};
template <> struct GMSHSimplexType<2U, 7U> {
static constexpr long type = 43;
};
template <> struct GMSHSimplexType<2U, 8U> {
static constexpr long type = 44;
};
template <> struct GMSHSimplexType<2U, 9U> {
static constexpr long type = 45;
};
template <> struct GMSHSimplexType<2U, 10U> {
static constexpr long type = 46;
};
template <> struct GMSHSimplexType<3U, 3U> {
static constexpr long type = 29;
};
template <> struct GMSHSimplexType<3U, 4U> {
static constexpr long type = 30;
};
template <> struct GMSHSimplexType<3U, 5U> {
static constexpr long type = 31;
};
template <> struct GMSHSimplexType<3U, 6U> {
static constexpr long type = 71;
};
template <> struct GMSHSimplexType<3U, 7U> {
static constexpr long type = 72;
};
template <> struct GMSHSimplexType<3U, 8U> {
static constexpr long type = 73;
};
template <> struct GMSHSimplexType<3U, 9U> {
static constexpr long type = 74;
};
template <> struct GMSHSimplexType<3U, 10U> {
static constexpr long type = 75;
};

template <std::size_t D> class GMSHBuilder : public tndm::GMSHMeshBuilder {
template <std::size_t D, std::size_t Order> class GMSHBuilder : public tndm::GMSHMeshBuilder {
public:
using vertex_t = std::array<double, D>;
using element_t = std::array<std::size_t, D + 1u>;
using facet_t = std::array<int, D>;
using element_t = std::array<std::size_t, nodeCount(D, Order)>;
using facet_t = std::array<int, nodeCount(D - 1u, Order)>;

std::vector<vertex_t> vertices;
std::vector<element_t> elements;
Expand All @@ -56,16 +166,16 @@ template <std::size_t D> class GMSHBuilder : public tndm::GMSHMeshBuilder {
bcs.reserve(numElements);
}
void addElement(long type, long tag, long* node, std::size_t numNodes) {
if (type == GMSHSimplexType<D>::type) {
assert(numNodes == D + 1u);
if (type == GMSHSimplexType<D, Order>::type) {
element_t elem;
std::copy(node, node + D + 1u, elem.begin());
assert(numNodes == elem.size());
std::copy_n(node, elem.size(), elem.begin());
elements.emplace_back(elem);
groups.emplace_back(static_cast<int>(tag));
} else if (type == GMSHSimplexType<D - 1u>::type) {
assert(numNodes == D);
} else if (type == GMSHSimplexType<D - 1u, Order>::type) {
facet_t elem;
std::copy(node, node + D, elem.begin());
assert(numNodes == elem.size());
std::copy_n(node, elem.size(), elem.begin());
facets.emplace_back(elem);
bcs.emplace_back(static_cast<int>(tag));
}
Expand Down
13 changes: 8 additions & 5 deletions src/meshreader/ParallelGMSHReader.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,15 +9,17 @@

#include <algorithm>
#include <cstddef>
#include <iostream>
#include <vector>

#include "aux/Distributor.h"
#include "aux/MPIConvenience.h"

namespace puml {

template <typename P> class ParallelGMSHReader {
template <typename P, std::size_t OrderP> class ParallelGMSHReader {
public:
constexpr static std::size_t Order = OrderP;
constexpr static std::size_t Dim = 3;
using bc_t = std::array<int, Dim + 1u>;
constexpr static std::size_t Facet2Nodes[Dim + 1][Dim] = {
Expand Down Expand Up @@ -49,11 +51,12 @@ template <typename P> class ParallelGMSHReader {
[[nodiscard]] std::size_t nVertices() const { return nVertices_; }
[[nodiscard]] std::size_t nElements() const { return nElements_; }
void readElements(std::size_t* elements) const {
static_assert(sizeof(GMSHBuilder<Dim>::element_t) == (Dim + 1) * sizeof(std::size_t));
scatter(builder_.elements.data()->data(), elements, nElements(), Dim + 1);
static_assert(sizeof(typename GMSHBuilder<Dim, Order>::element_t) ==
nodeCount(Dim, Order) * sizeof(std::size_t));
scatter(builder_.elements.data()->data(), elements, nElements(), nodeCount(Dim, Order));
}
void readVertices(double* vertices) const {
static_assert(sizeof(GMSHBuilder<Dim>::vertex_t) == Dim * sizeof(double));
static_assert(sizeof(typename GMSHBuilder<Dim, Order>::vertex_t) == Dim * sizeof(double));
scatter(builder_.vertices.data()->data(), vertices, nVertices(), Dim);
}
void readBoundaries(int* boundaries) const {
Expand Down Expand Up @@ -157,7 +160,7 @@ template <typename P> class ParallelGMSHReader {
}

MPI_Comm comm_;
GMSHBuilder<Dim> builder_;
GMSHBuilder<Dim, Order> builder_;
std::vector<bc_t> bcs_;
std::size_t nVertices_ = 0;
std::size_t nElements_ = 0;
Expand Down
2 changes: 2 additions & 0 deletions src/meshreader/ParallelMeshReader.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@ template <class R> class ParallelMeshReader {
R m_serialReader;

public:
constexpr static std::size_t Order = 1;
constexpr static std::size_t Dim = 3;
ParallelMeshReader(MPI_Comm comm = MPI_COMM_WORLD)
: m_nVertices(0), m_nElements(0), m_nBoundaries(0), m_comm(comm) {
init();
Expand Down
Loading