Commit a383ae33 authored by fabianw's avatar fabianw
Browse files

added non-uniform marching cube mesh

parent 2e900236
......@@ -32,45 +32,50 @@ public:
Interpolator(ArgumentParser& p) :
m_extent(p("extent").asDouble(1.0))
{
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") _getH = [](const Real e, const int N) { return e/static_cast<Real>(N); };
else if (p("data").asString("cell") == "node") _getH = [](const Real e, const int N) { return e/static_cast<Real>(N-1); };
// load the data
if (p("load").asString("h5") == "h5") _load_h5(p);
else if (p("load").asString("h5") == "wavelet") _load_wavelet(p);
else
{
std::cerr << "ERROR: Interpolator: Unknown data attribute \"" << p("data").asString() << "\"" << std::endl;
std::cerr << "ERROR: Interpolator: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
abort();
}
// load the data
if (p("load").asString("h5") == "h5") _load_h5(p);
else if (p("load").asString("h5") == "wavelet") _load_wavelet(p);
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") _nCells = [](const int N) { return N; };
else if (p("data").asString("cell") == "node") _nCells = [](const int N) { return N-1; };
else
{
std::cerr << "ERROR: Interpolator: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
std::cerr << "ERROR: Interpolator: Unknown data attribute \"" << p("data").asString() << "\"" << std::endl;
abort();
}
const int Nmax = std::max(m_data.Nx(), std::max(m_data.Ny(), m_data.Nz()));
assert(Nmax > 1);
m_h = _getH(m_extent, Nmax);
m_h = m_extent/_nCells(Nmax);
m_invh = 1.0/m_h;
m_extentX = m_h*_nCells(m_data.Nx());
m_extentY = m_h*_nCells(m_data.Ny());
m_extentZ = m_h*_nCells(m_data.Nz());
// init coordinate transform
m_origin_off = 0.0;
_physicalDataPos = [this](const int i) { return i*this->m_h; };
m_origin_off = {0.0, 0.0, 0.0};
m_isNodal = 1;
if (p("data").asString("cell") == "cell")
{
_physicalDataPos = [this](const int i) { return (0.5+i)*this->m_h; };
m_origin_off = {0.5*m_h, 0.5*m_h, 0.5*m_h};
m_isNodal = 0;
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_origin_off;
const Real yo = y - m_origin_off;
const Real zo = z - m_origin_off;
const Real xo = x - m_origin_off[0];
const Real yo = y - m_origin_off[1];
const Real zo = z - m_origin_off[2];
const int ix0 = _idx(xo);
const int iy0 = _idx(yo);
......@@ -94,6 +99,11 @@ public:
inline Real getH() const { return m_h; }
inline Real getExtent() const { return m_extent; }
inline Real getExtentX() const { return m_extentX; }
inline Real getExtentY() const { return m_extentY; }
inline Real getExtentZ() const { return m_extentZ; }
inline Real getDataPos(const int i) const { return _physicalDataPos(i); }
inline void setOrigin(const Real o[3]) { for (int i=0; i<3; ++i) m_origin_off[i] = o[i]; }
inline int getNx() const { return m_data.Nx(); }
inline int getNy() const { return m_data.Ny(); }
inline int getNz() const { return m_data.Nz(); }
......@@ -101,7 +111,10 @@ public:
protected:
Real m_extent;
Real m_origin_off;
Real m_extentX;
Real m_extentY;
Real m_extentZ;
Real m_origin_off[3];
Real m_h, m_invh;
int m_isNodal;
......@@ -109,7 +122,8 @@ protected:
TKernel m_kernel;
// helper
std::function<Real(const Real, const int)> _getH;
std::function<int(const int)> _nCells;
std::function<Real(const int)> _physicalDataPos;
inline int _idx(const Real x) const { return std::floor(x*m_invh); }
inline Real _pos(const int i) const { return i*m_h; }
......
......@@ -8,6 +8,7 @@
#include <cassert>
#include <vector>
#include <cmath>
#include <mpi.h>
#include "Interpolator.h"
......@@ -80,23 +81,23 @@ public:
MPI_Cart_coords(cartcomm, myrank, 3, m_myPEIndex);
m_isroot = (0 == myrank);
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") this->_getH = [](const Real e, const int N) { return e/static_cast<Real>(N); };
else if (p("data").asString("cell") == "node") this->_getH = [](const Real e, const int N) { return e/static_cast<Real>(N-1); };
// load the data
if (p("load").asString("h5") == "h5") this->_load_h5_MPI(p);
else if (p("load").asString("h5") == "wavelet") this->_load_wavelet_MPI(p);
else
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Unknown data attribute \"" << p("data").asString() << "\"" << std::endl;
std::cerr << "ERROR: InterpolatorMPI: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
abort();
}
// load the data
if (p("load").asString("h5") == "h5") this->_load_h5_MPI(p);
else if (p("load").asString("h5") == "wavelet") this->_load_wavelet_MPI(p);
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") this->_nCells = [](const int N) { return N; };
else if (p("data").asString("cell") == "node") this->_nCells = [](const int N) { return N-1; };
else
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
std::cerr << "ERROR: InterpolatorMPI: Unknown data attribute \"" << p("data").asString() << "\"" << std::endl;
abort();
}
......@@ -105,17 +106,32 @@ public:
const int Nz = m_PESize[2]*this->m_data.Nz();
const int Nmax = std::max(Nx, std::max(Ny, Nz));
assert(Nmax > 1);
this->m_h = this->_getH(this->m_extent, Nmax);
this->m_h = this->m_extent/this->_nCells(Nmax);
this->m_invh = 1.0/this->m_h;
this->m_extentX = this->m_h*this->_nCells(Nx);
this->m_extentY = this->m_h*this->_nCells(Ny);
this->m_extentZ = this->m_h*this->_nCells(Nz);
// init coordinate transform
this->m_origin_off = 0.0;
this->_physicalDataPos = [this](const int i) { return i*this->m_h; };
this->m_isNodal = 1;
if (p("data").asString("cell") == "cell")
{
this->_physicalDataPos = [this](const int i) { return (0.5+i)*this->m_h; };
this->m_isNodal = 0;
this->m_origin_off += 0.5*this->m_h;
}
const Real Ox = (this->m_extentX / m_PESize[0]) * m_myPEIndex[0];
const Real Oy = (this->m_extentY / m_PESize[1]) * m_myPEIndex[1];
const Real Oz = (this->m_extentZ / m_PESize[2]) * m_myPEIndex[2];
Real origin[3];
origin[0] = this->_physicalDataPos(ceil(Ox*this->m_invh)) - Ox;
origin[1] = this->_physicalDataPos(ceil(Oy*this->m_invh)) - Oy;
origin[2] = this->_physicalDataPos(ceil(Oz*this->m_invh)) - Oz;
assert(origin[0] >= 0.0);
assert(origin[1] >= 0.0);
assert(origin[2] >= 0.0);
for (int i = 0; i < 3; ++i)
this->m_origin_off[i] = origin[i];
// fetch halos cells
_fetch_halo(cartcomm);
......@@ -174,7 +190,8 @@ void InterpolatorMPI<Tinterp>::_load_h5_MPI(ArgumentParser& parser)
const std::string filename = parser("file").asString("");
if (filename == "")
{
std::cerr << "ERROR: Interpolator: -file is not specified. No input file given" << std::endl;
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: -file is not specified. No input file given" << std::endl;
abort();
}
......@@ -201,9 +218,24 @@ void InterpolatorMPI<Tinterp>::_load_h5_MPI(ArgumentParser& parser)
MPI_Bcast(maxDim, 4, MPI_INT, 0, m_comm);
NCH = maxDim[3];
assert(maxDim[0] % m_PESize[0] == 0);
assert(maxDim[1] % m_PESize[1] == 0);
assert(maxDim[2] % m_PESize[2] == 0);
if (maxDim[0] % m_PESize[0] != 0)
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Number of processes in x-direction is not an integer multiple of available data." << std::endl;
abort();
}
if (maxDim[1] % m_PESize[1] != 0)
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Number of processes in y-direction is not an integer multiple of available data." << std::endl;
abort();
}
if (maxDim[2] % m_PESize[2] != 0)
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Number of processes in z-direction is not an integer multiple of available data." << std::endl;
abort();
}
const size_t NX = maxDim[0]/m_PESize[0];
const size_t NY = maxDim[1]/m_PESize[1];
const size_t NZ = maxDim[2]/m_PESize[2];
......
......@@ -383,12 +383,14 @@ public:
const int Nx = static_cast<Real>(m_interp.getNx()) / cscale;
const int Ny = static_cast<Real>(m_interp.getNy()) / cscale;
const int Nz = static_cast<Real>(m_interp.getNz()) / cscale;
assert(Nx > 1);
assert(Ny > 1);
assert(Nz > 1);
const Real ch = m_interp.getExtent() / (std::max(Nx, std::max(Ny, Nz)) - m_interp.isNodal());
assert(Nx > 0);
assert(Ny > 0);
assert(Nz > 0);
const Real chx = m_interp.getExtentX() / Nx;
const Real chy = m_interp.getExtentY() / Ny;
const Real chz = m_interp.getExtentZ() / Nz;
_extractIso(isoval, Nx, Ny, Nz, ch);
_extractIso(isoval, Nx, Ny, Nz, chx, chy, chz);
_writePLY(fout, m_triList);
}
......@@ -411,10 +413,23 @@ protected:
return (t - v0)/dv;
}
inline void _extractIso(const Real isoval, const int Nx, const int Ny, const int Nz, const Real ch, const Real Ox=0.0, const Real Oy=0.0, const Real Oz=0.0, const bool verbose=true);
inline void _marchCube(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz);
inline void _marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz);
inline void _evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target);
inline void _extractIso(const Real isoval,
const int Nx, const int Ny, const int Nz,
const Real chx, const Real chy, const Real chz,
const Real Ox=0.0, const Real Oy=0.0, const Real Oz=0.0,
const bool verbose=true);
inline void _marchCube(const Real isoval,
const Real x, const Real y, const Real z,
const Real hx, const Real hy, const Real hz,
const Real Ox, const Real Oy, const Real Oz);
inline void _marchTetrahedron(const Real isoval,
const Real x, const Real y, const Real z,
const Real hx, const Real hy, const Real hz,
const Real Ox, const Real Oy, const Real Oz);
inline void _evalTetrahedron(const Real coords[4][3], const Real val[4], const Real isoval);
void _writePLY(const std::string basename, const std::vector<Triangle>& tList)
{
......@@ -496,7 +511,10 @@ protected:
template <typename Tkernel, template<typename> class Tinterpolator>
void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval, const int Nx, const int Ny, const int Nz, const Real ch, const Real Ox, const Real Oy, const Real Oz, const bool verbose)
void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval,
const int Nx, const int Ny, const int Nz,
const Real chx, const Real chy, const Real chz,
const Real Ox, const Real Oy, const Real Oz, const bool verbose)
{
const int speculate = m_parser("speculate_triangle").asInt(10000);
m_triList.clear();
......@@ -507,10 +525,10 @@ void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval, const i
for (int j = 0; j < Ny; ++j)
for (int i = 0; i < Nx; ++i)
{
const Real x = i*ch;
const Real y = j*ch;
const Real z = k*ch;
_marchCube(x, y, z, ch, isoval, Ox, Oy, Oz);
const Real x = i*chx;
const Real y = j*chy;
const Real z = k*chz;
_marchCube(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz);
}
}
else if (m_parser("march").asString("cube") == "tetrahedron")
......@@ -519,10 +537,10 @@ void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval, const i
for (int j = 0; j < Ny; ++j)
for (int i = 0; i < Nx; ++i)
{
const Real x = i*ch;
const Real y = j*ch;
const Real z = k*ch;
_marchTetrahedron(x, y, z, ch, isoval, Ox, Oy, Oz);
const Real x = i*chx;
const Real y = j*chy;
const Real z = k*chz;
_marchTetrahedron(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz);
}
}
else
......@@ -535,7 +553,10 @@ void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval, const i
template <typename Tkernel, template<typename> class Tinterpolator>
void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz)
void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real isoval,
const Real x, const Real y, const Real z,
const Real hx, const Real hy, const Real hz,
const Real Ox, const Real Oy, const Real Oz)
{
Real valCubeVerts[8];
Real coordEdgeVert[12][3];
......@@ -543,15 +564,15 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real x, const Real y,
// compute cube vertex values
for (int i = 0; i < 8; ++i)
valCubeVerts[i] = m_interp(
x + h*gVertexOffset[i][0],
y + h*gVertexOffset[i][1],
z + h*gVertexOffset[i][2]);
x + hx*gVertexOffset[i][0],
y + hy*gVertexOffset[i][1],
z + hz*gVertexOffset[i][2]);
// find if there are vertices inside and outside the target value
int flagVert = 0;
for (int i = 0; i < 8; ++i)
{
if (valCubeVerts[i] <= target) flagVert |= 1 << i;
if (valCubeVerts[i] <= isoval) flagVert |= 1 << i;
}
// indentify the edge intersections
......@@ -567,10 +588,10 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real x, const Real y,
if (flagEdge & (1 << i))
{
// linear interpolation
const Real frac = _getEdgeFrac(valCubeVerts[gEdgeConnection[i][0]], valCubeVerts[gEdgeConnection[i][1]], target);
coordEdgeVert[i][0] = x + h*(gVertexOffset[gEdgeConnection[i][0]][0] + frac*gEdgeDirection[i][0]) + m_x0 + Ox;
coordEdgeVert[i][1] = y + h*(gVertexOffset[gEdgeConnection[i][0]][1] + frac*gEdgeDirection[i][1]) + m_y0 + Oy;
coordEdgeVert[i][2] = z + h*(gVertexOffset[gEdgeConnection[i][0]][2] + frac*gEdgeDirection[i][2]) + m_z0 + Oz;
const Real frac = _getEdgeFrac(valCubeVerts[gEdgeConnection[i][0]], valCubeVerts[gEdgeConnection[i][1]], isoval);
coordEdgeVert[i][0] = x + hx*(gVertexOffset[gEdgeConnection[i][0]][0] + frac*gEdgeDirection[i][0]) + m_x0 + Ox;
coordEdgeVert[i][1] = y + hy*(gVertexOffset[gEdgeConnection[i][0]][1] + frac*gEdgeDirection[i][1]) + m_y0 + Oy;
coordEdgeVert[i][2] = z + hz*(gVertexOffset[gEdgeConnection[i][0]][2] + frac*gEdgeDirection[i][2]) + m_z0 + Oz;
}
}
......@@ -593,7 +614,10 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real x, const Real y,
template <typename Tkernel, template<typename> class Tinterpolator>
void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz)
void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real isoval,
const Real x, const Real y, const Real z,
const Real hx, const Real hy, const Real hz,
const Real Ox, const Real Oy, const Real Oz)
{
Real valCubeVerts[8];
Real coordCubeVerts[8][3];
......@@ -603,9 +627,9 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real x, const
// compute cube vertex coordinates
for (int i = 0; i < 8; ++i)
{
coordCubeVerts[i][0] = x + h*gVertexOffset[i][0];
coordCubeVerts[i][1] = y + h*gVertexOffset[i][1];
coordCubeVerts[i][2] = z + h*gVertexOffset[i][2];
coordCubeVerts[i][0] = x + hx*gVertexOffset[i][0];
coordCubeVerts[i][1] = y + hy*gVertexOffset[i][1];
coordCubeVerts[i][2] = z + hz*gVertexOffset[i][2];
}
// compute cube vertex values
......@@ -630,13 +654,13 @@ void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real x, const
coordTetraVerts[j][2] = coordCubeVerts[vertID][2];
valTetraVerts[j] = valCubeVerts[vertID];
}
_evalTetrahedron(coordTetraVerts, valTetraVerts, target);
_evalTetrahedron(coordTetraVerts, valTetraVerts, isoval);
}
}
template <typename Tkernel, template<typename> class Tinterpolator>
void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target)
void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][3], const Real val[4], const Real isoval)
{
Real coordEdgeVert[6][3];
......@@ -644,7 +668,7 @@ void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][
int flagVert = 0;
for (int i = 0; i < 4; ++i)
{
if (val[i] <= target) flagVert |= 1 << i;
if (val[i] <= isoval) flagVert |= 1 << i;
}
// indentify the edge intersections
......@@ -663,7 +687,7 @@ void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][
const int vert1 = gTetrahedronEdgeConnection[i][1];
// linear interpolation
const Real frac = _getEdgeFrac(val[vert0], val[vert1], target);
const Real frac = _getEdgeFrac(val[vert0], val[vert1], isoval);
const Real ifrac= 1.0 - frac;
coordEdgeVert[i][0] = ifrac*coords[vert0][0] + frac*coords[vert1][0];
......
......@@ -21,24 +21,40 @@ public:
{
const Real cscale = this->m_parser("cube_scale").asDouble(1.0);
assert(cscale > 0.0);
const int Nx = static_cast<Real>(this->m_interp.getLocalNx()) / cscale;
const int Ny = static_cast<Real>(this->m_interp.getLocalNy()) / cscale;
const int Nz = static_cast<Real>(this->m_interp.getLocalNz()) / cscale;
assert(Nx > 1);
assert(Ny > 1);
assert(Nz > 1);
int gNx = static_cast<Real>(this->m_interp.getNx()) / cscale;
int gNy = static_cast<Real>(this->m_interp.getNy()) / cscale;
int gNz = static_cast<Real>(this->m_interp.getNz()) / cscale;
int pesize[3];
this->m_interp.getPESize(pesize);
// TODO: (fabianw@mavt.ethz.ch; Tue Nov 29 23:51:01 2016) need better
// grid def here
const Real ch = this->m_interp.getExtent() / (std::max(pesize[0]*Nx, std::max(pesize[1]*Ny, pesize[2]*Nz)) - this->m_interp.isNodal());
gNx += (gNx % pesize[0]);
gNy += (gNy % pesize[1]);
gNz += (gNz % pesize[2]);
const int Nx = gNx/pesize[0];
const int Ny = gNy/pesize[1];
const int Nz = gNz/pesize[2];
assert(Nx > 0);
assert(Ny > 0);
assert(Nz > 0);
const Real chx = this->m_interp.getExtentX() / gNx;
const Real chy = this->m_interp.getExtentY() / gNy;
const Real chz = this->m_interp.getExtentZ() / gNz;
int peidx[3];
this->m_interp.getPEIndex(peidx);
const Real Ox = ch*Nx*peidx[0];
const Real Oy = ch*Ny*peidx[1];
const Real Oz = ch*Nz*peidx[2];
this->_extractIso(isoval, Nx, Ny, Nz, ch, Ox, Oy, Oz, this->m_interp.isroot());
const Real Ox = (this->m_interp.getExtentX()/pesize[0])*peidx[0];
const Real Oy = (this->m_interp.getExtentY()/pesize[1])*peidx[1];
const Real Oz = (this->m_interp.getExtentZ()/pesize[2])*peidx[2];
const Real data_invh = 1.0 / this->m_interp.getH();
Real origin[3];
origin[0] = this->m_interp.getDataPos(ceil(Ox*data_invh)) - Ox;
origin[1] = this->m_interp.getDataPos(ceil(Oy*data_invh)) - Oy;
origin[2] = this->m_interp.getDataPos(ceil(Oz*data_invh)) - Oz;
assert(origin[0] >= 0.0);
assert(origin[1] >= 0.0);
assert(origin[2] >= 0.0);
this->m_interp.setOrigin(origin);
this->_extractIso(isoval, Nx, Ny, Nz, chx, chy, chz, Ox, Oy, Oz, this->m_interp.isroot());
_writePLY_MPI(fout, this->m_triList);
}
......
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