Skip to content

Commit

Permalink
Merge pull request #161 from Pressio/fix_integer_type_gradients
Browse files Browse the repository at this point in the history 
fix integer type in gradients
  • Loading branch information
@fnrizzi
fnrizzi authored Sep 17, 2023
2 parents 6ae0bd4 + 737c317 commit 7f7a555
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 26 deletions.
4 changes: 2 additions & 2 deletions include/pressiodemoapps/gradient.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,12 +110,12 @@ class GradientEvaluator
it is not the actual name of the class.
You do not need to know the actual type, just use auto.
*/
const auto & queryFace(int cellGID, FacePosition fp) const{
const auto & queryFace(typename MeshType::index_t cellGID, FacePosition fp) const{
return m_impl.queryFace(cellGID, fp);
}

private:
using face_t = impl::Face<typename MeshType::scalar_type, MaxNumDofPerCell>;
using face_t = impl::Face<typename MeshType::scalar_type, typename MeshType::index_t, MaxNumDofPerCell>;
impl::GradientEvaluatorInternal<MeshType, face_t> m_impl;
};

Expand Down
52 changes: 28 additions & 24 deletions include/pressiodemoapps/impl/gradient_2d.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ namespace pressiodemoapps{ namespace impl{
this is intentionally just impl details, the actual public API is kept separate
*/

template<class ScalarType, class CellConnectivity, class FType>
template<class IntType, class ScalarType, class CellConnectivity, class FType>
ScalarType face_normal_gradient_for_cell_centered_function_2d(const FType & f,
const CellConnectivity & cc,
char axis,
Expand All @@ -75,11 +75,11 @@ ScalarType face_normal_gradient_for_cell_centered_function_2d(const FType & f,
// fd rules taken from https://web.media.mit.edu/~crtaylor/calculator.html
// manually verified only the 2pt one

const int i_05 = cc[0]; // this always exists
const int i_p15 = (axis=='x') ? cc[3] : (axis=='y') ? cc[2] : -1;
const int i_m15 = (axis=='x') ? cc[1] : (axis=='y') ? cc[4] : -1;
const int i_p30 = (graphNumCols>=9) ? ((axis=='x') ? cc[7] : ((axis=='y') ? cc[6] : -1)) : -1;
const int i_m30 = (graphNumCols>=9) ? ((axis=='x') ? cc[5] : ((axis=='y') ? cc[8] : -1)) : -1;
const IntType i_05 = cc[0]; // this always exists
const IntType i_p15 = (axis=='x') ? cc[3] : (axis=='y') ? cc[2] : -1;
const IntType i_m15 = (axis=='x') ? cc[1] : (axis=='y') ? cc[4] : -1;
const IntType i_p30 = (graphNumCols>=9) ? ((axis=='x') ? cc[7] : ((axis=='y') ? cc[6] : -1)) : -1;
const IntType i_m30 = (graphNumCols>=9) ? ((axis=='x') ? cc[5] : ((axis=='y') ? cc[8] : -1)) : -1;

const auto f_05 = f(i_05*numDofPerCell + dofShift);
const auto f_p15 = (i_p15 != -1) ? f(i_p15*numDofPerCell + dofShift) : 0;
Expand All @@ -103,28 +103,29 @@ ScalarType face_normal_gradient_for_cell_centered_function_2d(const FType & f,
constexpr int _faceNormalX = 1;
constexpr int _faceNormalY = 2;

template<class ScT, std::size_t N>
template<class ScT, class IndexType, std::size_t N>
struct Face{
std::array<ScT, 3> centerCoordinates = {};
std::array<ScT, N> normalGradient = {};
int parentCellGraphRow = {};
IndexType parentCellGraphRow = {};
int normalDirection = {}; //1 for x, 2 for y

static constexpr std::size_t N_ = N;
};

template<class ScT>
struct Face<ScT, 1>{
template<class ScT, class IndexType>
struct Face<ScT, IndexType, 1>{
std::array<ScT, 3> centerCoordinates = {};
ScT normalGradient = {};
int parentCellGraphRow = {};
IndexType parentCellGraphRow = {};
int normalDirection = {}; //1 for x, 2 for y

static constexpr std::size_t N_ = 1;
};

template<class IndexType>
struct MyKey{
int parentCellGID;
IndexType parentCellGID;
FacePosition pos;

bool operator==(const MyKey & p) const {
Expand All @@ -133,8 +134,9 @@ struct MyKey{
};

struct HashFnc{
std::size_t operator() (const MyKey & k) const{
const std::size_t h1 = std::hash<int>()(k.parentCellGID);
template<class IndexType>
std::size_t operator() (const MyKey<IndexType> & k) const{
const std::size_t h1 = std::hash<IndexType>()(k.parentCellGID);
const std::size_t h2 = std::hash<int>()(static_cast<int>(k.pos));
return h1 ^ h2;
}
Expand All @@ -145,6 +147,8 @@ template<class MeshType, class FaceType>
class GradientEvaluatorInternal
{
using sc_t = typename MeshType::scalar_type;
using index_t = typename MeshType::index_t;
using key_t = MyKey<index_t>;

public:
// constructor for when we only want the normal grad at domain boundaries' faces
Expand All @@ -156,8 +160,8 @@ class GradientEvaluatorInternal
initializeForStoringNormalGradsAtBoundaryFaces(mesh);
}

const auto & queryFace(int cellGID, FacePosition fp) const{
const MyKey key{cellGID, fp};
const auto & queryFace(index_t cellGID, FacePosition fp) const{
const key_t key{cellGID, fp};
assert(m_data.count(key) == 1);
auto it = m_data.find(key);
return it->second;
Expand Down Expand Up @@ -216,14 +220,14 @@ class GradientEvaluatorInternal

if constexpr (FaceType::N_ == 1){
face.normalGradient =
face_normal_gradient_for_cell_centered_function_2d(field, currCellConnec, axis,
face_normal_gradient_for_cell_centered_function_2d<index_t>(field, currCellConnec, axis,
fdMode, h, numDofPerCell, 0);
}
else{
for (int j=0; j<numDofPerCell; ++j){
face.normalGradient[j] =
face_normal_gradient_for_cell_centered_function_2d(field, currCellConnec, axis,
fdMode, h, numDofPerCell, j);
face_normal_gradient_for_cell_centered_function_2d<index_t>(field, currCellConnec, axis,
fdMode, h, numDofPerCell, j);
}
}
}
Expand Down Expand Up @@ -255,25 +259,25 @@ class GradientEvaluatorInternal
if (bL){
const auto faceCX = cellX-dxHalf;
const auto faceCY = cellY;
m_data[MyKey{cellGID, FacePosition::Left}] =
m_data[key_t{cellGID, FacePosition::Left}] =
FaceType{{faceCX, faceCY, faceCZ}, {}, rowInd, _faceNormalX};
}
if (bF){
const auto faceCX = cellX;
const auto faceCY = cellY+dyHalf;
m_data[MyKey{cellGID, FacePosition::Front}] =
m_data[key_t{cellGID, FacePosition::Front}] =
FaceType{{faceCX, faceCY, faceCZ}, {}, rowInd, _faceNormalY};
}
if (bR){
const auto faceCX = cellX+dxHalf;
const auto faceCY = cellY;
m_data[MyKey{cellGID, FacePosition::Right}] =
m_data[key_t{cellGID, FacePosition::Right}] =
FaceType{{faceCX, faceCY, faceCZ}, {}, rowInd, _faceNormalX};
}
if (bB){
const auto faceCX = cellX;
const auto faceCY = cellY-dyHalf;
m_data[MyKey{cellGID, FacePosition::Back}] =
m_data[key_t{cellGID, FacePosition::Back}] =
FaceType{{faceCX, faceCY, faceCZ}, {}, rowInd, _faceNormalY};
}
}
Expand All @@ -282,7 +286,7 @@ class GradientEvaluatorInternal
private:
BoundaryFacesNormalGradientScheme m_schemeAtBoundaryFaces;
std::reference_wrapper<const MeshType> m_meshObj;
std::unordered_map<MyKey, FaceType, HashFnc> m_data = {};
std::unordered_map<key_t, FaceType, HashFnc> m_data = {};
};

}}
Expand Down

0 comments on commit 7f7a555

Please sign in to comment.