commit a939059e03a9484df2330f1123a677b6eca010d2
parent 0a03d44d05e123755cf7de1a7a87897aa3db7d89
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Mon, 28 Nov 2016 15:25:12 +0100
added working version (non MPI)
Diffstat:
10 files changed, 1677 insertions(+), 0 deletions(-)
diff --git a/include/ArgumentParser.h b/include/ArgumentParser.h
@@ -0,0 +1,399 @@
+/*
+ * ArgumentParser.h
+ * Cubism
+ *
+ * This argument parser assumes that all arguments are optional ie, each of the argument names is preceded by a '-'
+ * all arguments are however NOT optional to avoid a mess with default values and returned values when not found!
+ *
+ * More converter could be required:
+ * add as needed
+ * TypeName as{TypeName}() in Value
+ *
+ * Created by Christian Conti on 6/7/10.
+ * Copyright 2010 ETH Zurich. All rights reserved.
+ *
+ */
+
+#pragma once
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <map>
+#include <vector>
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <limits>
+
+
+using namespace std;
+
+class Value
+{
+private:
+ string content;
+
+public:
+
+ Value() : content("") {}
+
+ Value(string content_) : content(content_) { /*printf("%s\n",content.c_str());*/ }
+
+ double asDouble(double def=0)
+ {
+ if (content == "")
+ {
+ ostringstream sbuf;
+ sbuf << def;
+ content = sbuf.str();
+ }
+ return (double) atof(content.c_str());
+ }
+
+ int asInt(int def=0)
+ {
+ if (content == "")
+ {
+ ostringstream sbuf;
+ sbuf << def;
+ content = sbuf.str();
+ }
+ return atoi(content.c_str());
+ }
+
+ bool asBool(bool def=false)
+ {
+ if (content == "")
+ {
+ if (def) content = "true";
+ else content = "false";
+ }
+ if (content == "0") return false;
+ if (content == "false") return false;
+
+ return true;
+ }
+
+ string asString(string def="")
+ {
+ if (content == "") content = def;
+
+ return content;
+ }
+};
+
+class CommandlineParser
+{
+private:
+ const int iArgC;
+ const char** vArgV;
+ bool bStrictMode, bVerbose;
+
+protected:
+ map<string,Value> mapArguments;
+
+public:
+
+ Value& operator()(const string arg)
+ {
+ if (bStrictMode)
+ {
+ map<string,Value>::const_iterator it = mapArguments.find(arg);
+
+ if (it == mapArguments.end())
+ {
+ printf("Runtime option NOT SPECIFIED! ABORTING! name: %s\n",arg.data());
+ abort();
+ }
+ }
+
+ if (bVerbose) printf("%s is %s\n", arg.data(), mapArguments[arg].asString().data());
+ return mapArguments[arg];
+ }
+
+ bool check(const string arg) const
+ {
+ return mapArguments.find(arg) != mapArguments.end();
+ }
+
+ CommandlineParser(const int argc, const char ** argv) : mapArguments(), iArgC(argc), vArgV(argv), bStrictMode(false), bVerbose(true)
+ {
+ for (int i=1; i<argc; i++)
+ if (argv[i][0] == '-')
+ {
+ string values = "";
+ int itemCount = 0;
+
+ for (int j=i+1; j<argc; j++)
+ {
+ const bool leadingDash = (argv[j][0] == '-');
+ const char c = argv[j][1];
+ const bool firstNumeric = ((c >= '0' && c <= '9') || c == 0) ? true : false;
+ if (leadingDash && !firstNumeric)
+ break;
+ else
+ {
+ if (strcmp(values.c_str(), ""))
+ values += ' ';
+
+ values += argv[j];
+ itemCount++;
+ }
+ }
+
+ if (itemCount == 0)
+ values = "true";
+
+ mapArguments[argv[i]] = Value(values);
+ i += itemCount;
+ }
+
+ mute();
+ //printf("found %ld arguments of %d\n",mapArguments.size(),argc);
+ }
+
+ int getargc() const { return iArgC; }
+
+ const char** getargv() const { return vArgV; }
+
+ void set_strict_mode()
+ {
+ bStrictMode = true;
+ }
+
+ void unset_strict_mode()
+ {
+ bStrictMode = false;
+ }
+
+ void mute()
+ {
+ bVerbose = false;
+ }
+
+ void loud()
+ {
+ bVerbose = true;
+ }
+
+ void save_options(string path=".")
+ {
+ string options;
+ for(map<string,Value>::iterator it=mapArguments.begin(); it!=mapArguments.end(); it++)
+ {
+ options+= it->first + " " + it->second.asString() + " ";
+ }
+ string filepath = (path + "/" + string("argumentparser.log"));
+ FILE * f = fopen(filepath.data(), "a");
+ if (f == NULL)
+ {
+ printf("impossible to write %s.\n", filepath.data());
+ return;
+ }
+ fprintf(f, "%s\n", options.data());
+ fclose(f);
+ }
+
+ void print_args()
+ {
+ for(map<string,Value>::iterator it=mapArguments.begin(); it!=mapArguments.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << it->first;
+ std::cout << ": " << it->second.asString() << std::endl;
+ }
+ }
+};
+
+
+class ArgumentParser: public CommandlineParser
+{
+ typedef std::map<std::string, Value> ArgMap;
+ typedef std::map<std::string, Value*> pArgMap;
+ typedef std::map<std::string, ArgMap* > FileMap;
+
+ const char commentStart;
+
+ // keep a reference form option origin
+ ArgMap from_commandline;
+ FileMap from_files;
+ pArgMap from_code;
+
+ // helper
+ void _ignoreComments(std::istream& stream, const char commentChar)
+ {
+ stream >> std::ws;
+ int nextchar = stream.peek();
+ while (nextchar == commentChar)
+ {
+ stream.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+ stream >> std::ws;
+ nextchar = stream.peek();
+ }
+ }
+
+ bool _existKey(std::string& key) const
+ {
+ const std::string og_key = key;
+ bool bExist = true;
+
+ // look for both possible keys (i.e. with leading "-" and without)
+ ArgMap::const_iterator it = mapArguments.find(key);
+ if (it == mapArguments.end())
+ {
+ if (key[0] == '-') key = key.erase(0, 1);
+ else key = "-" + key;
+ it = mapArguments.find(key);
+ if (it == mapArguments.end())
+ {
+ key = og_key;
+ bExist = false;
+ }
+ }
+ return bExist;
+ }
+
+ inline std::string _stripKey(std::string key) const
+ {
+ if (key[0] == '-') key = key.erase(0, 1);
+ return key;
+ }
+
+
+public:
+ ArgumentParser(const int _argc, const char ** _argv):
+ CommandlineParser(_argc, _argv), commentStart('#')
+ {
+ from_commandline = mapArguments;
+ }
+
+ virtual ~ArgumentParser()
+ {
+ for (FileMap::iterator it = from_files.begin(); it != from_files.end(); it++)
+ delete it->second;
+ }
+
+ void readFile(const std::string filepath)
+ {
+ from_files[filepath] = new ArgMap;
+ ArgMap& myFMap = *(from_files[filepath]);
+
+ std::ifstream confFile(filepath.c_str());
+ if (confFile.is_open())
+ {
+ // read (key value) pairs from input file, ignore comments
+ // beginning with "#"
+ _ignoreComments(confFile, commentStart);
+ while (!confFile.eof())
+ {
+ std::string line, key, val;
+ std::getline(confFile, line);
+ std::istringstream lineStream(line);
+ lineStream >> key;
+ lineStream >> val;
+ _ignoreComments(lineStream, commentStart);
+ while(!lineStream.eof())
+ {
+ std::string multiVal;
+ lineStream >> multiVal;
+ val += (" " + multiVal);
+ _ignoreComments(lineStream, commentStart);
+ }
+ if (_existKey(key)) continue;
+
+ std::pair<string, Value> item(key, Value(val));
+ mapArguments.insert(item); // add to parent container
+ myFMap.insert(item); // add to private container
+
+ _ignoreComments(confFile, commentStart);
+ }
+ }
+ }
+
+ Value& operator()(std::string key)
+ {
+ const bool bDefaultInCode = !_existKey(key);
+ Value& retval = CommandlineParser::operator()(key);
+ if (bDefaultInCode) from_code[key] = &retval;
+ return retval;
+ }
+
+ inline bool check(std::string key) const { return _existKey(key); }
+ inline bool exist(std::string key) const { return _existKey(key); }
+
+ void print_args()
+ {
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ std::cout << "* Summary:" << std::endl;
+ std::cout << "* Parameter read from command line: " << from_commandline.size() << std::endl;
+ size_t nFiles = 0;
+ size_t nFileParameter = 0;
+ for (FileMap::const_iterator it=from_files.begin(); it!=from_files.end(); ++it)
+ {
+ if (it->second->size() > 0)
+ {
+ ++nFiles;
+ nFileParameter += it->second->size();
+ }
+ }
+ std::cout << "* Parameter read from " << std::setw(3) << std::right << nFiles << " file(s): " << nFileParameter << std::endl;
+ std::cout << "* Parameter read from defaults in code: " << from_code.size() << std::endl;
+ std::cout << "* Total number of parameter read from all sources: " << mapArguments.size() << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+
+ // command line given arguments
+ if (!from_commandline.empty())
+ {
+ std::cout << "* Command Line:" << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ for(ArgMap::iterator it=from_commandline.begin(); it!=from_commandline.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << _stripKey(it->first);
+ std::cout << ": " << it->second.asString() << std::endl;
+ }
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ }
+
+ // options read from input files
+ if (!from_files.empty())
+ {
+ for (FileMap::iterator itFile=from_files.begin(); itFile!=from_files.end(); itFile++)
+ {
+ if (!itFile->second->empty())
+ {
+ std::cout << "* File: " << itFile->first << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ ArgMap& fileArgs = *(itFile->second);
+ for(ArgMap::iterator it=fileArgs.begin(); it!=fileArgs.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << _stripKey(it->first);
+ std::cout << ": " << it->second.asString() << std::endl;
+ }
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ }
+ }
+ }
+
+ // defaults defined in code
+ if (!from_code.empty())
+ {
+ std::cout << "* Defined in Code:" << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ for(pArgMap::iterator it=from_code.begin(); it!=from_code.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << _stripKey(it->first);
+ std::cout << ": " << it->second->asString() << std::endl;
+ }
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ }
+ }
+};
diff --git a/include/Matrix3D.h b/include/Matrix3D.h
@@ -0,0 +1,170 @@
+// File : Matrix3D.h
+// Date : Wed Nov 23 14:45:22 2016
+// Author : Fabian Wermelinger
+// Description: 3D Matrix
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef MATRIX3D_H_L9MVKQGE
+#define MATRIX3D_H_L9MVKQGE
+
+#include <cassert>
+#include <cstdlib>
+#include <cstddef>
+#include <cstring>
+#include <iostream>
+
+template <typename T, int XS=0, int XE=0, int YS=0, int YE=0, int ZS=0, int ZE=0>
+class Matrix3D
+{
+public:
+ 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) :
+ m_Nx(0), m_Ny(0), m_Nz(0), m_allocated(false), m_data(nullptr)
+ {
+ alloc(Nx, Ny, Nz, data);
+ }
+
+ Matrix3D(const Matrix3D& rhs) :
+ m_Nx(0), m_Ny(0), m_Nz(0), m_allocated(false), m_data(nullptr)
+ {
+ if (rhs.m_allocated)
+ {
+ m_data = _alloc(rhs.m_Nx, rhs.m_Ny, rhs.m_Nz);
+ _copy(rhs.m_data);
+ }
+ }
+
+ Matrix3D(Matrix3D&& rhs) :
+ m_Nx(rhs.m_Nx), m_Ny(rhs.m_Ny), m_Nz(rhs.m_Nz),
+ m_XS(rhs.m_XS), m_XE(rhs.m_XE), m_Xpitch(rhs.m_Xpitch), m_Npitched(rhs.m_Npitched),
+ m_allocated(rhs.m_allocated), m_data(std::move(rhs.m_data))
+ { }
+
+ ~Matrix3D() { if (m_allocated) free(m_data); }
+
+ Matrix3D& operator=(const Matrix3D& rhs)
+ {
+ if (this != &rhs)
+ {
+ _reset();
+ if (rhs.m_allocated)
+ {
+ m_data = _alloc(rhs.m_Nx, rhs.m_Ny, rhs.m_Nz);
+ _copy(rhs.m_data);
+ }
+ }
+ return *this;
+ }
+
+ Matrix3D& operator=(Matrix3D&& rhs)
+ {
+ if (this != &rhs)
+ {
+ _reset();
+ if (rhs.m_allocated)
+ {
+ m_Nx = rhs.m_Nx;
+ m_Ny = rhs.m_Ny;
+ m_Nz = rhs.m_Nz;
+ m_XS = rhs.m_XS;
+ m_XE = rhs.m_XE;
+ m_Xpitch = rhs.m_Xpitch;
+ m_Npitched = rhs.m_Npitched;
+ m_allocated = rhs.m_allocated;
+ m_data = std::move(rhs.m_data);
+ rhs.m_allocated = false;
+ }
+ }
+ return *this;
+ }
+
+ T operator()(const int ix, const int iy, const int iz) const
+ {
+ assert(m_allocated);
+ assert(ix >= XS && ix < m_Nx+XE);
+ assert(iy >= YS && iy < m_Ny+YE);
+ assert(iz >= ZS && iz < m_Nz+ZE);
+ return m_data[(ix-XS) + m_Xpitch*(iy-YS) + m_Xpitch*(m_Ny+YE-YS)*(iz-ZS)];
+ }
+
+ T& operator()(const int ix, const int iy, const int iz)
+ {
+ assert(m_allocated);
+ assert(ix >= XS && ix < m_Nx+XE);
+ assert(iy >= YS && iy < m_Ny+YE);
+ assert(iz >= ZS && iz < m_Nz+ZE);
+ return m_data[(ix-XS) + m_Xpitch*(iy-YS) + m_Xpitch*(m_Ny+YE-YS)*(iz-ZS)];
+ }
+
+ inline void alloc(const int Nx, const int Ny, const int Nz, const T* const data = nullptr)
+ {
+ if (m_allocated) free(m_data);
+ m_data = _alloc(Nx, Ny, Nz);
+
+ if (data)
+ {
+ for (int iz = 0; iz < m_Nz; ++iz)
+ for (int iy = 0; iy < m_Ny; ++iy)
+ {
+ const T* const src = data + m_Nx*iy + m_Nx*m_Ny*iz;
+ T* dst = &(this->operator()(0,iy,iz));
+ memcpy(dst, src, m_Nx*sizeof(T));
+ }
+ }
+ }
+
+ inline int Nx() const { return m_Nx; }
+ inline int Ny() const { return m_Ny; }
+ inline int Nz() const { return m_Nz; }
+
+private:
+ int m_Nx;
+ int m_Ny;
+ int m_Nz;
+ int m_XS, m_XE, m_Xpitch, m_Npitched;
+ bool m_allocated;
+ T* m_data;
+
+ inline T* _alloc(const int Nx, const int Ny, const int Nz)
+ {
+ m_Nx = Nx;
+ m_Ny = Ny;
+ m_Nz = Nz;
+ const int alignedElements = _ALIGN_/sizeof(T);
+ m_XS = alignedElements * (XS - (alignedElements-1))/alignedElements;
+ m_XE = alignedElements * (XE + (alignedElements-1))/alignedElements;
+ m_Xpitch = (m_Nx+m_XE) - m_XS;
+ m_Npitched = m_Xpitch * (m_Ny+YE-YS) * (m_Nz+ZE-ZS);
+
+ void* pmem;
+ if (posix_memalign(&pmem, _ALIGN_, m_Npitched*sizeof(T)))
+ {
+ std::cerr << "ERROR: Matrix3D: can not allocate memory" << std::endl;
+ abort();
+ }
+ memset(pmem, 0, m_Npitched*sizeof(T));
+ m_allocated = true;
+ return (T*)pmem;
+ }
+
+ inline void _copy(const T* const base)
+ {
+ const int Nslice = m_Xpitch * (m_Ny+YE-YS);
+ for (int iz = 0; iz < m_Nz+ZE-ZS; ++iz)
+ {
+ const T* const src = base + Nslice*iz;
+ T* dst = m_data + Nslice*iz;
+ memcpy(dst, src, Nslice*sizeof(T));
+ }
+ }
+
+ inline void _reset()
+ {
+ if (m_allocated) free(m_data);
+ m_Nx = m_Ny = m_Nz = m_XS = m_XE = m_Xpitch = m_Npitched = 0;
+ m_allocated = false;
+ m_data = nullptr;
+ }
+};
+
+#endif /* MATRIX3D_H_L9MVKQGE */
diff --git a/include/common.h b/include/common.h
@@ -0,0 +1,30 @@
+// File : common.h
+// Date : Tue Nov 22 15:51:55 2016
+// Author : Fabian Wermelinger
+// Description: common stuff
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef COMMON_H_ZH0IVQBA
+#define COMMON_H_ZH0IVQBA
+
+#ifdef _FLOAT_PRECISION_
+typedef float Real;
+#else
+typedef double Real;
+#endif
+
+#ifdef _USE_HDF_
+#ifdef _FLOAT_PRECISION_
+#define HDF_PRECISION H5T_NATIVE_FLOAT
+#else
+#define HDF_PRECISION H5T_NATIVE_DOUBLE
+#endif
+#endif /* _USE_HDF_ */
+
+#ifndef _ALIGN_
+#define _ALIGN_ 16
+#endif /* _ALIGN_ */
+
+#include "Matrix3D.h"
+typedef Matrix3D<Real,-3,3,-3,3,-3,3> Matrix_t;
+
+#endif /* COMMON_H_ZH0IVQBA */
diff --git a/src/Interpolator.h b/src/Interpolator.h
@@ -0,0 +1,262 @@
+// File : Interpolator.h
+// Date : Tue Nov 22 15:37:00 2016
+// Author : Fabian Wermelinger
+// Description: Data interpolator type
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef INTERPOLATOR_H_QCK6H0CI
+#define INTERPOLATOR_H_QCK6H0CI
+
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <cmath>
+#include <cstring>
+#include <string>
+#include <functional>
+#ifdef _USE_HDF_
+#include <hdf5.h>
+#endif /* _USE_HDF_ */
+
+// TODO: (fabianw@mavt.ethz.ch; Thu Nov 24 08:21:14 2016) testing
+#include <fstream>
+
+#include "common.h"
+#include "ArgumentParser.h"
+#ifdef _USE_CUBISMZ_
+#include "Reader_WaveletCompression.h"
+#endif /* _USE_CUBISMZ_ */
+
+
+template <typename TKernel>
+class Interpolator
+{
+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); };
+ else
+ {
+ std::cerr << "ERROR: Interpolator: Unknown data attribute \"" << p("data").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);
+ else
+ {
+ std::cerr << "ERROR: Interpolator: Undefined loader \"" << p("load").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_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_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;
+ }
+ }
+ 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 int ix0 = _idx(xo);
+ const int iy0 = _idx(yo);
+ const int iz0 = _idx(zo);
+
+ Real interp = 0.0;
+ for (int iz = TKernel::start; iz < TKernel::end; ++iz)
+ for (int iy = TKernel::start; iy < TKernel::end; ++iy)
+ for (int ix = TKernel::start; ix < TKernel::end; ++ix)
+ {
+ const Real xj = _pos(ix0+ix);
+ const Real yj = _pos(iy0+iy);
+ const Real zj = _pos(iz0+iz);
+ const Real Mx = m_kernel((xo-xj)*m_invh);
+ const Real My = m_kernel((yo-yj)*m_invh);
+ const Real Mz = m_kernel((zo-zj)*m_invh);
+ interp += m_data(ix0+ix,iy0+iy,iz0+iz)*Mx*My*Mz;
+ }
+ 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; }
+
+ // TODO: (fabianw@mavt.ethz.ch; Thu Nov 24 16:38:31 2016) trash this
+ // void dump_data(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 int iy = _idx(yo);
+ // const int iz = _idx(zo);
+
+ // ofstream idata("raw_data.dat");
+ // for (int ix = 0; ix < m_data.Nx(); ++ix)
+ // idata << _pos(ix)+m_xorigin_off << "\t" << m_data(ix,iy,iz) << std::endl;
+ // }
+
+
+private:
+ const Real m_extent;
+ Real m_xorigin_off;
+ Real m_yorigin_off;
+ Real m_zorigin_off;
+ Real m_h, m_invh;
+ int m_isNodal;
+
+ Matrix_t m_data;
+ TKernel m_kernel;
+
+ // loader
+ void _load_h5(ArgumentParser& p);
+ void _load_wavelet(ArgumentParser& p);
+
+ // helper
+ std::function<Real(const Real, const int)> _getH;
+ inline int _idx(const Real x) const { return x*m_invh; }
+ inline Real _pos(const int i) const { return i*m_h; }
+};
+
+
+template <typename TKernel>
+void Interpolator<TKernel>::_load_h5(ArgumentParser& parser)
+{
+#ifdef _USE_HDF_
+ const std::string group = parser("hdf5_group").asString("/data");
+ const std::string filename = parser("file").asString("");
+ if (filename == "")
+ {
+ std::cerr << "ERROR: Interpolator: -file is not specified. No input file given" << std::endl;
+ abort();
+ }
+
+ hid_t file_id, dataset_id, dataspace_id, file_dataspace_id;
+ hsize_t* dims;
+ hssize_t num_elem;
+ int rank, ndims, NCH;
+ int maxDim[3];
+ Real* data;
+
+ file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ dataset_id = H5Dopen(file_id, group.c_str(), H5P_DEFAULT);
+ file_dataspace_id = H5Dget_space(dataset_id);
+ rank = H5Sget_simple_extent_ndims(file_dataspace_id);
+ dims = new hsize_t[rank];
+ ndims = H5Sget_simple_extent_dims(file_dataspace_id, dims, NULL);
+ num_elem = H5Sget_simple_extent_npoints(file_dataspace_id);
+ data = new Real[num_elem];
+ maxDim[2] = dims[0];
+ maxDim[1] = dims[1];
+ maxDim[0] = dims[2];
+ NCH = dims[3];
+ dataspace_id = H5Screate_simple(rank, dims, NULL);
+ int status = H5Dread(dataset_id, HDF_PRECISION, dataspace_id, file_dataspace_id, H5P_DEFAULT, data);
+
+ /* release stuff */
+ delete [] dims;
+ status = H5Dclose(dataset_id);
+ status = H5Sclose(dataspace_id);
+ status = H5Sclose(file_dataspace_id);
+ status = H5Fclose(file_id);
+
+ if (NCH > 1)
+ {
+ const int channel = parser("channel").asInt(0); // data channel
+ const bool magnitude = parser("magnitude").asBool(false); // vector magnitude (only if NCH == 3)
+ if (magnitude && NCH == 3)
+ {
+ int k = 0;
+ for (int i = 0; i < num_elem; i += NCH)
+ {
+ const Real a = data[i+0];
+ const Real b = data[i+1];
+ const Real c = data[i+2];
+ data[k++] = std::sqrt(a*a + b*b + c*c);
+ }
+
+ }
+ else
+ {
+ assert(channel < NCH);
+ int k = 0;
+ for (int i = 0; i < num_elem; i += NCH)
+ data[k++] = data[i+channel];
+ }
+ }
+ this->m_data = std::move(Matrix_t(maxDim[0], maxDim[1], maxDim[2], data));
+ delete [] data;
+#endif /* _USE_HDF_ */
+}
+
+template <typename TKernel>
+void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p)
+{
+#ifdef _USE_CUBISMZ_
+ const bool byte_swap = p("-swap").asBool(false);
+ const int wavelet_type = p("-wtype").asInt(1);
+ const std::string filename = p("file").asString("");
+ if (filename == "")
+ {
+ std::cerr << "ERROR: Interpolator: -file is not specified. No input file given" << std::endl;
+ abort();
+ }
+
+ const string fname(filename);
+ Reader_WaveletCompression myreader(fname, byte_swap, wavelet_type);
+ myreader.load_file();
+ const int NBX = myreader.xblocks();
+ const int NBY = myreader.yblocks();
+ const int NBZ = myreader.zblocks();
+ const int maxDim[3] = {NBX*_BLOCKSIZE_, NBY*_BLOCKSIZE_, NBZ*_BLOCKSIZE_};
+ Real* const data = new Real[maxDim[0]*maxDim[1]*maxDim[2]];
+
+ Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
+ for (int iz=0; iz < NBZ; ++iz)
+ for (int iy=0; iy < NBY; ++iy)
+ for (int ix=0; ix < NBX; ++ix)
+ {
+ const double zratio = myreader.load_block2(ix, iy, iz, blockdata);
+
+ for (int z=0; z < _BLOCKSIZE_; ++z)
+ for (int y=0; y < _BLOCKSIZE_; ++y)
+ {
+ assert(iy*_BLOCKSIZE_+y < maxDim[1]);
+ assert(iz*_BLOCKSIZE_+z < maxDim[2]);
+ Real* const dst = data + _BLOCKSIZE_*(ix + NBX*(y+iy*_BLOCKSIZE_ + (z+iz*_BLOCKSIZE_)*NBY*_BLOCKSIZE_));
+ memcpy(dst, &blockdata[z][y][0], _BLOCKSIZE_*sizeof(Real));
+ }
+ }
+
+ this->m_data = std::move(Matrix_t(maxDim[0], maxDim[1], maxDim[2], data));
+ delete [] data;
+#else
+ fprintf(stderr, "WARNING: Executable was compiled without wavelet compressor support...\n");
+#endif /* _USE_CUBISMZ_ */
+}
+
+#endif /* INTERPOLATOR_H_QCK6H0CI */
diff --git a/src/IsoExtractor.h b/src/IsoExtractor.h
@@ -0,0 +1,681 @@
+// File : IsoExtractor.h
+// Date : Fri Nov 25 19:12:52 2016
+// Author : Fabian Wermelinger
+// Description: Iso surface extractor
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef ISOEXTRACTOR_H_A1WQHFFC
+#define ISOEXTRACTOR_H_A1WQHFFC
+
+#include <vector>
+#include <fstream>
+#include <string>
+#include <ctime>
+#include <cstdlib>
+#include "common.h"
+#include "ArgumentParser.h"
+#include "Interpolator.h"
+
+struct Triangle
+{
+ Real vertices[3][3];
+};
+
+static const Real gVertexOffset[8][3] =
+{
+ {0.0, 0.0, 0.0},{1.0, 0.0, 0.0},{1.0, 1.0, 0.0},{0.0, 1.0, 0.0},
+ {0.0, 0.0, 1.0},{1.0, 0.0, 1.0},{1.0, 1.0, 1.0},{0.0, 1.0, 1.0}
+};
+
+static const int gEdgeConnection[12][2] =
+{
+ {0,1}, {1,2}, {2,3}, {3,0},
+ {4,5}, {5,6}, {6,7}, {7,4},
+ {0,4}, {1,5}, {2,6}, {3,7}
+};
+
+static const Real gEdgeDirection[12][3] =
+{
+ {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},
+ {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},
+ {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0, 0.0, 1.0}
+};
+
+static const int gTetrahedronEdgeConnection[6][2] =
+{
+ {0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}
+};
+
+static const int gTetrahedronsInACube[6][4] =
+{
+ {0,5,1,6},
+ {0,1,2,6},
+ {0,2,3,6},
+ {0,3,7,6},
+ {0,7,4,6},
+ {0,4,5,6},
+};
+
+static const int gCubeEdgeFlags[256] =
+{
+ 0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
+ 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
+ 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
+ 0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
+ 0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
+ 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
+ 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
+ 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
+ 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
+ 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
+ 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
+ 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
+ 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
+ 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
+ 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
+ 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000
+};
+
+static const int gTriangleConnectionTable[256][16] =
+{
+ {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
+ {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
+ {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+ {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
+ {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
+ {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
+ {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
+ {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
+ {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
+ {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
+ {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
+ {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
+ {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
+ {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
+ {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
+ {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
+ {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
+ {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
+ {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
+ {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
+ {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
+ {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
+ {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
+ {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
+ {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
+ {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
+ {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
+ {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
+ {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
+ {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
+ {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
+ {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
+ {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
+ {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
+ {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
+ {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
+ {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
+ {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
+ {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
+ {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
+ {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
+ {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
+ {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
+ {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
+ {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
+ {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
+ {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
+ {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
+ {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
+ {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
+ {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
+ {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
+ {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
+ {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
+ {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
+ {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
+ {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
+ {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
+ {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
+ {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
+ {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
+ {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
+ {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
+ {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
+ {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
+ {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
+ {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
+ {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
+ {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
+ {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
+ {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
+ {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
+ {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
+ {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
+ {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
+ {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
+ {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
+ {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
+ {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
+ {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
+ {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
+ {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
+ {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
+ {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
+ {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
+ {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
+ {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
+ {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
+ {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
+ {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
+ {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
+ {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
+ {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
+ {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
+ {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
+ {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
+ {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
+ {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
+ {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
+ {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
+ {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
+ {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
+ {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
+ {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
+ {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
+ {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
+ {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
+ {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
+ {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
+ {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
+ {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
+ {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
+ {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
+ {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
+ {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
+ {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
+ {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
+ {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
+ {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
+ {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
+ {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
+ {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
+ {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
+ {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
+ {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
+ {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
+ {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
+ {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
+ {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
+ {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
+ {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
+ {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
+ {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
+ {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
+ {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
+ {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
+ {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
+ {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
+ {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
+ {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
+ {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
+ {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
+ {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
+ {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
+ {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
+ {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
+ {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
+ {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
+ {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
+ {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
+ {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
+ {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
+ {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
+ {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
+ {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
+ {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
+ {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
+ {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
+ {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
+ {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
+ {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
+ {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
+ {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+ {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
+ {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
+ {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
+ {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
+};
+
+static const int gTetrahedronEdgeFlags[16] =
+{
+ 0x00, 0x0d, 0x13, 0x1e, 0x26, 0x2b, 0x35, 0x38, 0x38, 0x35, 0x2b, 0x26, 0x1e, 0x13, 0x0d, 0x00
+};
+
+static const int gTetrahedronTriangles[16][7] =
+{
+ {-1, -1, -1, -1, -1, -1, -1},
+ { 0, 3, 2, -1, -1, -1, -1},
+ { 0, 1, 4, -1, -1, -1, -1},
+ { 1, 4, 2, 2, 4, 3, -1},
+
+ { 1, 2, 5, -1, -1, -1, -1},
+ { 0, 3, 5, 0, 5, 1, -1},
+ { 0, 2, 5, 0, 5, 4, -1},
+ { 5, 4, 3, -1, -1, -1, -1},
+
+ { 3, 4, 5, -1, -1, -1, -1},
+ { 4, 5, 0, 5, 2, 0, -1},
+ { 1, 5, 0, 5, 3, 0, -1},
+ { 5, 2, 1, -1, -1, -1, -1},
+
+ { 3, 4, 2, 2, 4, 1, -1},
+ { 4, 1, 0, -1, -1, -1, -1},
+ { 2, 3, 0, -1, -1, -1, -1},
+ {-1, -1, -1, -1, -1, -1, -1}
+};
+
+
+template <typename Tinterp>
+class IsoExtractor
+{
+public:
+ IsoExtractor(ArgumentParser& parser) :
+ m_parser(parser),
+ m_interp(parser),
+ m_x0(parser("xmin").asDouble(0.0)),
+ m_y0(parser("ymin").asDouble(0.0)),
+ m_z0(parser("zmin").asDouble(0.0)) {}
+ ~IsoExtractor() {}
+
+ void extract(const Real isoval, const std::string fout)
+ {
+ const Real x0 = m_parser("xmin").asDouble(0.0);
+ const Real y0 = m_parser("ymin").asDouble(0.0);
+ const Real z0 = m_parser("zmin").asDouble(0.0);
+
+ const Real cscale = m_parser("cube_scale").asDouble(1.0);
+ assert(cscale > 0.0);
+ 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());
+
+ m_triList.clear();
+
+ if (m_parser("march").asString("cube") == "cube")
+ {
+ for (int k = 0; k < Nz; ++k)
+ 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);
+ }
+ }
+ else if (m_parser("march").asString("cube") == "tetrahedron")
+ {
+ for (int k = 0; k < Nz; ++k)
+ 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);
+ }
+ }
+ else
+ {
+ std::cerr << "ERROR: IsoExtractor: Unknown algorithm " << m_parser("march").asString() << std::endl;
+ abort();
+ }
+
+ _writePLY(fout, m_triList);
+ }
+
+
+private:
+ ArgumentParser& m_parser;
+ Interpolator<Tinterp> m_interp;
+ const Real m_x0;
+ const Real m_y0;
+ const Real m_z0;
+
+ std::vector<Triangle> m_triList;
+
+ // helper
+ inline Real _getEdgeFrac(const Real v0, const Real v1, const Real t)
+ {
+ const Real dv = v1 - v0;
+ if (dv == 0.0) return 0.5;
+ return (t - v0)/dv;
+ }
+
+ inline void _marchCube(const Real x, const Real y, const Real z, const Real h, const Real target);
+ inline void _marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target);
+ inline void _evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target);
+
+ void _writePLY(const std::string basename, const std::vector<Triangle>& tList)
+ {
+ time_t rawtime;
+ std::time(&rawtime);
+ struct tm* timeinfo = std::localtime(&rawtime);
+ char buf[256];
+ std::strftime(buf, 256, "%A, %h %d %Y, %r %Z %z", timeinfo);
+
+ const unsigned int one = 1;
+ const unsigned char* const lsb = (const unsigned char*)(&one);
+ const std::string endianess((*lsb) ? "binary_little_endian" : "binary_big_endian");
+
+ // write header
+ std::ofstream ply(basename+".ply");
+ ply << "ply" << std::endl;
+ ply << "format ";
+ if (m_parser("binary").asBool(true))
+ ply << endianess;
+ else
+ ply << "ascii";
+ ply << " 1.0" << std::endl;
+ ply << "comment generated by IsoExtractor.h on " << buf << std::endl;
+ ply << "element vertex " << tList.size()*3 << std::endl;
+ ply << "property float x" << std::endl;
+ ply << "property float y" << std::endl;
+ ply << "property float z" << std::endl;
+ ply << "element face " << tList.size() << std::endl;
+ ply << "property list uchar int vertex_index" << std::endl;
+ ply << "end_header" << std::endl;
+ ply.close();
+
+ // write content
+ if (m_parser("binary").asBool(true))
+ {
+ ply.open(basename+".ply", std::ios::out | std::ios::app | std::ios::binary);
+ _writePLY_binary(ply, tList);
+ }
+ else
+ {
+ ply.open(basename+".ply", std::ios::out | std::ios::app | std::ios::binary);
+ _writePLY_ascii(ply, tList);
+ }
+ ply.close();
+ }
+
+ void _writePLY_binary(std::ostream& ply, const std::vector<Triangle>& tList)
+ {
+ for (auto t : tList)
+ {
+ float tc[3][3];
+ for (int j = 0; j < 3; ++j)
+ {
+ tc[j][0] = t.vertices[j][0];
+ tc[j][1] = t.vertices[j][1];
+ tc[j][2] = t.vertices[j][2];
+ }
+ ply.write((char*)&tc[0][0], 9*sizeof(float));
+ }
+
+ const char nVert = 3;
+ for (int i = 0; i < tList.size()*3; i += 3)
+ {
+ int tv[3] = {i, i+1, i+2};
+ ply.write(&nVert, sizeof(char));
+ ply.write((char*)&tv[0], 3*sizeof(int));
+ }
+ }
+
+ void _writePLY_ascii(std::ostream& ply, const std::vector<Triangle>& tList)
+ {
+ for (auto t : tList)
+ for (int j = 0; j < 3; ++j)
+ ply << t.vertices[j][0] << " " << t.vertices[j][1] << " " << t.vertices[j][2] << std::endl;
+ for (int i = 0; i < tList.size()*3; i += 3)
+ ply << "3 " << i << " " << i+1 << " " << i+2 << std::endl;
+ }
+};
+
+
+template <typename Tinterp>
+void IsoExtractor<Tinterp>::_marchCube(const Real x, const Real y, const Real z, const Real h, const Real target)
+{
+ Real valCubeVerts[8];
+ Real coordEdgeVert[12][3];
+
+ // 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]);
+
+ // 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;
+ }
+
+ // indentify the edge intersections
+ const int flagEdge = gCubeEdgeFlags[flagVert];
+
+ // if there are no intersections at all, bye
+ if(flagEdge == 0)
+ return;
+
+ // compute coordinates of triangle vertices on intersected edges
+ for (int i = 0; i < 12; ++i)
+ {
+ 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]);
+ coordEdgeVert[i][1] = y + h*(gVertexOffset[gEdgeConnection[i][0]][1] + frac*gEdgeDirection[i][1]);
+ coordEdgeVert[i][2] = z + h*(gVertexOffset[gEdgeConnection[i][0]][2] + frac*gEdgeDirection[i][2]);
+ }
+ }
+
+ // push the triangles found into the list
+ for (int i = 0; i < 5; ++i)
+ {
+ if (gTriangleConnectionTable[flagVert][3*i] < 0) break;
+
+ Triangle tri;
+ for (int j = 0; j < 3; ++j)
+ {
+ const int k = gTriangleConnectionTable[flagVert][3*i + j];
+ tri.vertices[j][0] = coordEdgeVert[k][0];
+ tri.vertices[j][1] = coordEdgeVert[k][1];
+ tri.vertices[j][2] = coordEdgeVert[k][2];
+ }
+ m_triList.push_back(tri);
+ }
+}
+
+
+template <typename Tinterp>
+void IsoExtractor<Tinterp>::_marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target)
+{
+ Real valCubeVerts[8];
+ Real coordCubeVerts[8][3];
+ Real valTetraVerts[4];
+ Real coordTetraVerts[4][3];
+
+ // 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];
+ }
+
+ // compute cube vertex values
+ for (int i = 0; i < 8; ++i)
+ valCubeVerts[i] = m_interp(
+ coordCubeVerts[i][0],
+ coordCubeVerts[i][1],
+ coordCubeVerts[i][2]);
+
+ for (int i = 0; i < 6; ++i)
+ {
+ for (int j = 0; j < 4; ++j)
+ {
+ const int vertID = gTetrahedronsInACube[i][j];
+ coordTetraVerts[j][0] = coordCubeVerts[vertID][0];
+ coordTetraVerts[j][1] = coordCubeVerts[vertID][1];
+ coordTetraVerts[j][2] = coordCubeVerts[vertID][2];
+ valTetraVerts[j] = valCubeVerts[vertID];
+ }
+ _evalTetrahedron(coordTetraVerts, valTetraVerts, target);
+ }
+}
+
+
+template <typename Tinterp>
+void IsoExtractor<Tinterp>::_evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target)
+{
+ Real coordEdgeVert[6][3];
+
+ // find if there are vertices inside and outside the target value
+ int flagVert = 0;
+ for (int i = 0; i < 4; ++i)
+ {
+ if (val[i] <= target) flagVert |= 1 << i;
+ }
+
+ // indentify the edge intersections
+ const int flagEdge = gTetrahedronEdgeFlags[flagVert];
+
+ // if there are no intersections at all, bye
+ if(flagEdge == 0)
+ return;
+
+ // compute coordinates of triangle vertices on intersected edges
+ for (int i = 0; i < 6; ++i)
+ {
+ if (flagEdge & (1 << i))
+ {
+ const int vert0 = gTetrahedronEdgeConnection[i][0];
+ const int vert1 = gTetrahedronEdgeConnection[i][1];
+
+ // linear interpolation
+ const Real frac = _getEdgeFrac(val[vert0], val[vert1], target);
+ const Real ifrac= 1.0 - frac;
+
+ coordEdgeVert[i][0] = ifrac*coords[vert0][0] + frac*coords[vert1][0];
+ coordEdgeVert[i][1] = ifrac*coords[vert0][1] + frac*coords[vert1][1];
+ coordEdgeVert[i][2] = ifrac*coords[vert0][2] + frac*coords[vert1][2];
+ }
+ }
+
+ // push the triangles found into the list
+ for (int i = 0; i < 2; ++i)
+ {
+ if (gTetrahedronTriangles[flagVert][3*i] < 0) break;
+
+ Triangle tri;
+ for (int j = 0; j < 3; ++j)
+ {
+ const int k = gTetrahedronTriangles[flagVert][3*i + j];
+ tri.vertices[j][0] = coordEdgeVert[k][0];
+ tri.vertices[j][1] = coordEdgeVert[k][1];
+ tri.vertices[j][2] = coordEdgeVert[k][2];
+ }
+ m_triList.push_back(tri);
+ }
+}
+
+#endif /* ISOEXTRACTOR_H_A1WQHFFC */
diff --git a/src/easy_iso.cpp b/src/easy_iso.cpp
@@ -0,0 +1,31 @@
+// File : easy_iso.cpp
+// Date : Tue Nov 22 15:00:14 2016
+// Author : Fabian Wermelinger
+// Description: Isosurface extractor
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#include <string>
+#include <sstream>
+#include "ArgumentParser.h"
+#include "interpolation/kernels.h"
+#include "IsoExtractor.h"
+
+int main(int argc, char* argv[])
+{
+ ArgumentParser parser(argc, const_cast<const char**>(argv));
+
+ using Tinterp = M2;
+ if (parser("interp_kernel").asString("M2") == "Mp4")
+ using Tinterp = Mprime4;
+ else if (parser("interp_kernel").asString("M2") == "M5")
+ using Tinterp = M5;
+
+ IsoExtractor<Tinterp> iso(parser);
+
+ const Real isoval = parser("isoval").asDouble(0.0);
+ std::ostringstream basename;
+ basename << "isosurf_" << isoval;
+ const std::string fout = parser("fout").asString(basename.str());
+ iso.extract(isoval, fout);
+
+ return 0;
+}
diff --git a/src/interpolation/M2.h b/src/interpolation/M2.h
@@ -0,0 +1,28 @@
+// File : M2.h
+// Date : Thu Nov 24 16:14:51 2016
+// Author : Fabian Wermelinger
+// Description: M2 interpolation kernel (Monaghan 19885) -- triangular function
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef M2_H_LQ69GRM8
+#define M2_H_LQ69GRM8
+
+#include <cmath>
+#include "common.h"
+
+class M2
+{
+public:
+ static constexpr int start = 0;
+ static constexpr int end = 2;
+
+ inline Real operator()(const Real x) const
+ {
+ const Real IxI = std::abs(x);
+ if (IxI <= 1.0)
+ return 1.0 - IxI;
+ else
+ return 0.0;
+ }
+};
+
+#endif /* M2_H_LQ69GRM8 */
diff --git a/src/interpolation/M5.h b/src/interpolation/M5.h
@@ -0,0 +1,33 @@
+// File : M5.h
+// Date : Thu Nov 24 09:33:24 2016
+// Author : Fabian Wermelinger
+// Description: M5 smoothing interpolation kernel (Monaghan 1985)
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef M5_H_G6FRPHTC
+#define M5_H_G6FRPHTC
+
+#include <cmath>
+#include "common.h"
+
+class M5
+{
+public:
+ static constexpr int start = -2;
+ static constexpr int end = 4;
+
+ inline Real operator()(const Real x) const
+ {
+ const Real a0 = 1.0/48.0;
+ const Real IxI = std::abs(x);
+ if (IxI < 0.5)
+ return a0*(43.125-75.0*IxI*IxI+42.0*IxI*IxI*IxI*IxI);
+ else if (IxI < 1.5)
+ return a0*(41.25+20.0*IxI-150.0*IxI*IxI+120.0*IxI*IxI*IxI-28.0*IxI*IxI*IxI*IxI);
+ else if (IxI < 2.5)
+ return a0*(IxI-2.5)*(IxI-2.5)*(IxI-2.5)*(7.0*IxI-7.5);
+ else
+ return 0.0;
+ }
+};
+
+#endif /* M5_H_G6FRPHTC */
diff --git a/src/interpolation/Mprime4.h b/src/interpolation/Mprime4.h
@@ -0,0 +1,30 @@
+// File : Mprime4.h
+// Date : Thu Nov 24 08:33:08 2016
+// Author : Fabian Wermelinger
+// Description: M^\prime_4 interpolation kernel (3rd order)
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef MPRIME4_H_KXJZRENA
+#define MPRIME4_H_KXJZRENA
+
+#include <cmath>
+#include "common.h"
+
+class Mprime4
+{
+public:
+ static constexpr int start = -1;
+ static constexpr int end = 3;
+
+ inline Real operator()(const Real x) const
+ {
+ const Real IxI = std::abs(x);
+ if (IxI < 1.0)
+ return 0.5*(IxI - 1.0)*(3.0*IxI*IxI - 2.0*IxI - 2.0);
+ else if (IxI < 2.0)
+ return -0.5*(IxI - 1.0)*(IxI - 2.0)*(IxI - 2.0);
+ else
+ return 0.0;
+ }
+};
+
+#endif /* MPRIME4_H_KXJZRENA */
diff --git a/src/interpolation/kernels.h b/src/interpolation/kernels.h
@@ -0,0 +1,13 @@
+// File : kernels.h
+// Date : Thu Nov 24 16:47:33 2016
+// Author : Fabian Wermelinger
+// Description: Interpolation kernels
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef KERNELS_H_IRWGGDVC
+#define KERNELS_H_IRWGGDVC
+
+#include "interpolation/M2.h"
+#include "interpolation/Mprime4.h"
+#include "interpolation/M5.h"
+
+#endif /* KERNELS_H_IRWGGDVC */
