Commit f4fba5b8 authored by fabianw's avatar fabianw
Browse files

minor mods in several places

parent 9c8675a7
......@@ -16,6 +16,13 @@ template <typename T, int XS=0, int XE=0, int YS=0, int YE=0, int ZS=0, int ZE=0
class Matrix3D
{
public:
static constexpr int XSTART = XS;
static constexpr int XEND = XE;
static constexpr int YSTART = YS;
static constexpr int YEND = YE;
static constexpr int ZSTART = ZS;
static constexpr int ZEND = ZE;
Matrix3D() : m_Nx(0), m_Ny(0), m_Nz(0), m_allocated(false), m_data(nullptr) { }
Matrix3D(const int Nx, const int Ny, const int Nz, const T* const data = nullptr) :
......
......@@ -6,6 +6,8 @@
#ifndef COMMON_H_ZH0IVQBA
#define COMMON_H_ZH0IVQBA
#include <cassert>
#ifdef _FLOAT_PRECISION_
typedef float Real;
#else
......@@ -25,6 +27,138 @@ typedef double Real;
#endif /* _ALIGN_ */
#include "Matrix3D.h"
template <typename T, size_t P>
struct Face
{
T* data;
const int _N;
const int _M;
const int _Nelem;
static constexpr int _P = P;
Face(const size_t N, const size_t M) :
data(new T[N*M*P]),
_N(N), _M(M), _Nelem(N*M*P)
{}
~Face() { delete [] data; }
inline T& operator()(const size_t ix, const size_t iy, const size_t iz)
{
assert(ix < _N);
assert(iy < _M);
assert(iz < P);
return data[ix + iy*_N + iz*_N*_M];
}
inline T operator()(const size_t ix, const size_t iy, const size_t iz) const
{
assert(ix < _N);
assert(iy < _M);
assert(iz < P);
return data[ix + iy*_N + iz*_N*_M];
}
};
template <typename T, size_t P, size_t Q>
struct Edge
{
T* data;
const int _N;
const int _Nelem;
static constexpr int _P = P;
static constexpr int _Q = Q;
Edge(const size_t N) :
data(new T[N*P*Q]),
_N(N), _Nelem(N*P*Q)
{}
~Edge() { delete [] data; }
inline T& operator()(const size_t ix, const size_t iy, const size_t iz)
{
assert(ix < _N);
assert(iy < P);
assert(iz < Q);
return data[ix + iy*_N + iz*_N*P];
}
inline T operator()(const size_t ix, const size_t iy, const size_t iz) const
{
assert(ix < _N);
assert(iy < P);
assert(iz < Q);
return data[ix + iy*_N + iz*_N*P];
}
};
template <typename T, size_t P, size_t Q, size_t R>
struct Corner
{
T* data;
const int _Nelem;
static constexpr int _P = P;
static constexpr int _Q = Q;
static constexpr int _R = R;
Corner() : data(new T[P*Q*R]), _Nelem(P*Q*R) {}
~Corner() { delete [] data; }
inline T& operator()(const size_t ix, const size_t iy, const size_t iz)
{
assert(ix < P);
assert(iy < Q);
assert(iz < R);
return data[ix + iy*P + iz*P*Q];
}
inline T operator()(const size_t ix, const size_t iy, const size_t iz) const
{
assert(ix < P);
assert(iy < Q);
assert(iz < R);
return data[ix + iy*P + iz*P*Q];
}
};
typedef Matrix3D<Real,-3,3,-3,3,-3,3> Matrix_t;
typedef Face<Real,3> Face_t;
typedef Edge<Real,3,3> Edge_t;
typedef Corner<Real,3,3,3> Corner_t;
struct Faces
{
// indices 0 and 1 correspond to the near and far faces relative to the
// origin.
Face_t x0, x1, y0, y1, z0, z1;
Faces(const size_t XN, const size_t XM,
const size_t YN, const size_t YM,
const size_t ZN, const size_t ZM) :
x0(XN,XM), x1(XN,XM), y0(YN,YM), y1(YN,YM), z0(ZN,ZM), z1(ZN,ZM)
{}
};
struct Edges
{
// edges are along the given coordinate direction. The counting is done
// using right hand rule for the indicated direction. The 0-edge is the
// one that goes through the origin.
Edge_t x0, x1, x2, x3, y0, y1, y2, y3, z0, z1, z2, z3;
Edges(const size_t XN, const size_t YN, const size_t ZN) :
x0(XN), x1(XN), x2(XN), x3(XN),
y0(YN), y1(YN), y2(YN), y3(YN),
z0(ZN), z1(ZN), z2(ZN), z3(ZN)
{}
};
struct Corners
{
// corners are identified by x0 or x1 (near far). The second index is
// again given by the right hand rule along the x-axis. The x{0,1}0 are
// the ones along the axis that goes through the origin.
Corner_t x00, x01, x02, x03, x10, x11, x12, x13;
};
#endif /* COMMON_H_ZH0IVQBA */
......@@ -56,25 +56,21 @@ public:
m_invh = 1.0/m_h;
// init coordinate transform
m_xorigin_off = p("xmin").asDouble(0.0);
m_yorigin_off = p("ymin").asDouble(0.0);
m_zorigin_off = p("zmin").asDouble(0.0);
m_origin_off = 0.0;
m_isNodal = 1;
if (p("data").asString("cell") == "cell")
{
m_isNodal = 0;
m_xorigin_off += 0.5*m_h;
m_yorigin_off += 0.5*m_h;
m_zorigin_off += 0.5*m_h;
m_origin_off += 0.5*m_h;
}
}
virtual ~Interpolator() {}
Real operator()(const Real x, const Real y, const Real z) const
{
const Real xo = x - m_xorigin_off;
const Real yo = y - m_yorigin_off;
const Real zo = z - m_zorigin_off;
const Real xo = x - m_origin_off;
const Real yo = y - m_origin_off;
const Real zo = z - m_origin_off;
const int ix0 = _idx(xo);
const int iy0 = _idx(yo);
......@@ -96,18 +92,16 @@ public:
return interp;
}
Real getH() const { return m_h; }
Real getExtent() const { return m_extent; }
int getNx() const { return m_data.Nx(); }
int getNy() const { return m_data.Ny(); }
int getNz() const { return m_data.Nz(); }
int isNodal() const { return m_isNodal; }
inline Real getH() const { return m_h; }
inline Real getExtent() const { return m_extent; }
inline int getNx() const { return m_data.Nx(); }
inline int getNy() const { return m_data.Ny(); }
inline int getNz() const { return m_data.Nz(); }
inline int isNodal() const { return m_isNodal; }
protected:
Real m_extent;
Real m_xorigin_off;
Real m_yorigin_off;
Real m_zorigin_off;
Real m_origin_off;
Real m_h, m_invh;
int m_isNodal;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment