diff options
Diffstat (limited to 'extern/mantaflow/helper')
26 files changed, 7743 insertions, 0 deletions
diff --git a/extern/mantaflow/helper/pwrapper/manta.h b/extern/mantaflow/helper/pwrapper/manta.h new file mode 100644 index 00000000000..efbca6cc493 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/manta.h @@ -0,0 +1,31 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011-2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Include pwrapper headers + * + ******************************************************************************/ + +#ifndef _MANTA_H +#define _MANTA_H + +// Remove preprocessor keywords, so there won't infere with autocompletion etc. +#define KERNEL(...) extern int i, j, k, idx, X, Y, Z; +#define PYTHON(...) +#define returns(X) extern X; +#define alias typedef + +#include "general.h" +#include "vectorbase.h" +#include "vector4d.h" +#include "registry.h" +#include "pclass.h" +#include "pconvert.h" +#include "fluidsolver.h" + +#endif diff --git a/extern/mantaflow/helper/pwrapper/numpyWrap.cpp b/extern/mantaflow/helper/pwrapper/numpyWrap.cpp new file mode 100644 index 00000000000..d2ddb21be70 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/numpyWrap.cpp @@ -0,0 +1,132 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2017-2018 Steffen Wiewel, Moritz Becher, Rachel Chu + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Convert mantaflow grids to/from numpy arrays + * + ******************************************************************************/ + +#include "manta.h" +#include "pythonInclude.h" + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include "numpy/arrayobject.h" + +namespace Manta { + +#if PY_VERSION_HEX < 0x03000000 +PyMODINIT_FUNC initNumpy() +{ + import_array(); +} +#endif + +// ------------------------------------------------------------------------ +// Class Functions +// ------------------------------------------------------------------------ +PyArrayContainer::PyArrayContainer(void *_pParentPyArray) : pParentPyArray(_pParentPyArray) +{ + ExtractData(pParentPyArray); +} +// ------------------------------------------------------------------------ +PyArrayContainer::PyArrayContainer(const PyArrayContainer &_Other) + : pParentPyArray(_Other.pParentPyArray) +{ + ExtractData(pParentPyArray); + Py_INCREF(pParentPyArray); +} +// ------------------------------------------------------------------------ +PyArrayContainer::~PyArrayContainer() +{ + Py_DECREF(pParentPyArray); +} +// ------------------------------------------------------------------------ +PyArrayContainer &PyArrayContainer::operator=(const PyArrayContainer &_Other) +{ + if (this != &_Other) { + // DecRef the existing resource + Py_DECREF(pParentPyArray); + + // Relink new data + pParentPyArray = _Other.pParentPyArray; + ExtractData(pParentPyArray); + Py_INCREF(pParentPyArray); + } + return *this; +} +// ------------------------------------------------------------------------ +void PyArrayContainer::ExtractData(void *_pParentPyArray) +{ + PyArrayObject *pParent = reinterpret_cast<PyArrayObject *>(pParentPyArray); + + int numDims = PyArray_NDIM(pParent); + long *pDims = (long *)PyArray_DIMS(pParent); + + pData = PyArray_DATA(pParent); + TotalSize = PyArray_SIZE(pParent); + Dims = std::vector<long>(&pDims[0], &pDims[numDims]); + + int iDataType = PyArray_TYPE(pParent); + switch (iDataType) { + case NPY_FLOAT: + DataType = N_FLOAT; + break; + case NPY_DOUBLE: + DataType = N_DOUBLE; + break; + case NPY_INT: + DataType = N_INT; + break; + default: + errMsg("unknown type of Numpy array"); + break; + } +} + +// ------------------------------------------------------------------------ +// Conversion Functions +// ------------------------------------------------------------------------ + +template<> PyArrayContainer fromPy<PyArrayContainer>(PyObject *obj) +{ + if (PyArray_API == NULL) { + // python 3 uses the return value +#if PY_VERSION_HEX >= 0x03000000 + import_array(); +#else + initNumpy(); +#endif + } + + if (!PyArray_Check(obj)) { + errMsg("argument is not an numpy array"); + } + + PyArrayObject *obj_p = reinterpret_cast<PyArrayObject *>( + PyArray_CheckFromAny(obj, + NULL, + 0, + 0, + /*NPY_ARRAY_ENSURECOPY*/ NPY_ARRAY_C_CONTIGUOUS | + NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_NOTSWAPPED, + NULL)); + PyArrayContainer container = PyArrayContainer(obj_p); + + return container; +} + +// template<> PyArrayContainer* fromPyPtr<PyArrayContainer>(PyObject* obj, std::vector<void*>* tmp) +// { +// if (!tmp) throw Error("dynamic de-ref not supported for this type"); +// void* ptr = malloc(sizeof(PyArrayContainer)); +// tmp->push_back(ptr); + +// *((PyArrayContainer*) ptr) = fromPy<PyArrayContainer>(obj); +// return (PyArrayContainer*) ptr; +// } +} // namespace Manta diff --git a/extern/mantaflow/helper/pwrapper/numpyWrap.h b/extern/mantaflow/helper/pwrapper/numpyWrap.h new file mode 100644 index 00000000000..c92a2eaaa97 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/numpyWrap.h @@ -0,0 +1,86 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2017 Steffen Wiewel, Moritz Baecher, Rachel Chu + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Convert mantaflow grids to/from numpy arrays + * + ******************************************************************************/ + +#ifdef _PCONVERT_H +# ifndef _NUMPYCONVERT_H +# define _NUMPYCONVERT_H + +enum NumpyTypes { + N_BOOL = 0, + N_BYTE, + N_UBYTE, + N_SHORT, + N_USHORT, + N_INT, + N_UINT, + N_LONG, + N_ULONG, + N_LONGLONG, + N_ULONGLONG, + N_FLOAT, + N_DOUBLE, + N_LONGDOUBLE, + N_CFLOAT, + N_CDOUBLE, + N_CLONGDOUBLE, + N_OBJECT = 17, + N_STRING, + N_UNICODE, + N_VOID, + /* + * New 1.6 types appended, may be integrated + * into the above in 2.0. + */ + N_DATETIME, + N_TIMEDELTA, + N_HALF, + + N_NTYPES, + N_NOTYPE, + N_CHAR, /* special flag */ + N_USERDEF = 256, /* leave room for characters */ + + /* The number of types not including the new 1.6 types */ + N_NTYPES_ABI_COMPATIBLE = 21 +}; + +namespace Manta { +class PyArrayContainer { + public: + /// Constructors + PyArrayContainer(void *_pParentPyArray); + PyArrayContainer(const PyArrayContainer &_Other); + ~PyArrayContainer(); + /// Operators + PyArrayContainer &operator=(const PyArrayContainer &_Other); + + private: + void ExtractData(void *_pParentPyArray); + + public: + void *pData; + NumpyTypes DataType; + unsigned int TotalSize; + std::vector<long> Dims; + + private: + void *pParentPyArray; +}; + +// template<> PyArrayContainer* fromPyPtr<PyArrayContainer>(PyObject* obj, std::vector<void*>* +// tmp); +template<> PyArrayContainer fromPy<PyArrayContainer>(PyObject *obj); +} // namespace Manta + +# endif +#endif diff --git a/extern/mantaflow/helper/pwrapper/pclass.cpp b/extern/mantaflow/helper/pwrapper/pclass.cpp new file mode 100644 index 00000000000..a95254ebe11 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/pclass.cpp @@ -0,0 +1,220 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Functions for property setting/getting via python + * + ******************************************************************************/ + +#include "pythonInclude.h" +#include "structmember.h" +#include "manta.h" +#include "general.h" +#include "timing.h" + +#ifdef GUI +# include <QMutex> +#else +class QMutex { + public: + void lock(){}; + void unlock(){}; + bool tryLock() + { + return true; + }; +}; +#endif + +using namespace std; +namespace Manta { + +//****************************************************************************** +// Free functions + +void pbPreparePlugin(FluidSolver *parent, const string &name, bool doTime) +{ + if (doTime) + TimingData::instance().start(parent, name); +} + +void pbFinalizePlugin(FluidSolver *parent, const string &name, bool doTime) +{ + if (doTime) + TimingData::instance().stop(parent, name); + + // GUI update, also print name of parent if there's more than one + std::ostringstream msg; + if (name != "FluidSolver::step") { + if (parent && (parent->getNumInstances() > 0)) + msg << parent->getName() << string("."); + msg << name; + } + updateQtGui(false, 0, 0., msg.str()); + + debMsg(name << " done", 3); + // name unnamed PbClass Objects from var name + PbClass::renameObjects(); +} + +void pbSetError(const string &fn, const string &ex) +{ + debMsg("Error in " << fn, 1); + if (!ex.empty()) + PyErr_SetString(PyExc_RuntimeError, ex.c_str()); +} + +//****************************************************************************** +// Helpers + +string PbTypeVec::str() const +{ + if (T.empty()) + return ""; + string s = "<"; + for (int i = 0; i < (int)T.size(); i++) { + s += T[i].str(); + s += (i != (int)T.size() - 1) ? ',' : '>'; + } + return s; +} +string PbType::str() const +{ + if (S == "float") + return "Real"; + if (S == "manta.vec3") + return "Vec3"; + return S; +} + +//****************************************************************************** +// PbClass + +vector<PbClass *> PbClass::mInstances; + +PbClass::PbClass(FluidSolver *parent, const string &name, PyObject *obj) + : mMutex(NULL), mParent(parent), mPyObject(obj), mName(name), mHidden(false) +{ + mMutex = new QMutex(); +} + +PbClass::PbClass(const PbClass &a) + : mMutex(NULL), mParent(a.mParent), mPyObject(0), mName("_unnamed"), mHidden(false) +{ + mMutex = new QMutex(); +} + +PbClass::~PbClass() +{ + for (vector<PbClass *>::iterator it = mInstances.begin(); it != mInstances.end(); ++it) { + if (*it == this) { + mInstances.erase(it); + break; + } + } + delete mMutex; +} + +void PbClass::lock() +{ + mMutex->lock(); +} +void PbClass::unlock() +{ + mMutex->unlock(); +} +bool PbClass::tryLock() +{ + return mMutex->tryLock(); +} + +PbClass *PbClass::getInstance(int idx) +{ + if (idx < 0 || idx > (int)mInstances.size()) + errMsg("PbClass::getInstance(): invalid index"); + return mInstances[idx]; +} + +int PbClass::getNumInstances() +{ + return mInstances.size(); +} + +bool PbClass::isNullRef(PyObject *obj) +{ + return PyLong_Check(obj) && PyLong_AsDouble(obj) == 0; +} + +bool PbClass::isNoneRef(PyObject *obj) +{ + return (obj == Py_None); +} + +void PbClass::registerObject(PyObject *obj, PbArgs *args) +{ + // cross link + Pb::setReference(this, obj); + mPyObject = obj; + + mInstances.push_back(this); + + if (args) { + string _name = args->getOpt<std::string>("name", -1, ""); + if (!_name.empty()) + setName(_name); + } +} + +PbClass *PbClass::createPyObject(const string &classname, + const string &name, + PbArgs &args, + PbClass *parent) +{ + return Pb::createPy(classname, name, args, parent); +} + +void PbClass::checkParent() +{ + if (getParent() == NULL) { + errMsg("New class " + mName + ": no parent given -- specify using parent=xxx !"); + } +} +//! Assign unnamed PbClass objects their Python variable name +void PbClass::renameObjects() +{ + PyObject *sys_mod_dict = PyImport_GetModuleDict(); + PyObject *loc_mod = PyMapping_GetItemString(sys_mod_dict, (char *)"__main__"); + if (!loc_mod) + return; + PyObject *locdict = PyObject_GetAttrString(loc_mod, "__dict__"); + if (!locdict) + return; + + // iterate all PbClass instances + for (size_t i = 0; i < mInstances.size(); i++) { + PbClass *obj = mInstances[i]; + if (obj->getName().empty()) { + // empty, try to find instance in module local dictionary + + PyObject *lkey, *lvalue; + Py_ssize_t lpos = 0; + while (PyDict_Next(locdict, &lpos, &lkey, &lvalue)) { + if (lvalue == obj->mPyObject) { + string varName = fromPy<string>(PyObject_Str(lkey)); + obj->setName(varName); + // cout << "assigning variable name '" << varName << "' to unnamed instance" << endl; + break; + } + } + } + } + Py_DECREF(locdict); + Py_DECREF(loc_mod); +} + +} // namespace Manta diff --git a/extern/mantaflow/helper/pwrapper/pclass.h b/extern/mantaflow/helper/pwrapper/pclass.h new file mode 100644 index 00000000000..b34103ca9a7 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/pclass.h @@ -0,0 +1,126 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011-2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Base class for all Python-exposed classes + * + ******************************************************************************/ + +// ----------------------------------------------------------------- +// NOTE: +// Do not include this file in user code, include "manta.h" instead +// ----------------------------------------------------------------- + +#ifdef _MANTA_H +# ifndef _PTYPE_H +# define _PTYPE_H + +# include <string> +# include <vector> +# include <map> + +class QMutex; + +namespace Manta { +struct PbClassData; +class FluidSolver; +class PbArgs; + +struct PbType { + std::string S; + std::string str() const; +}; +struct PbTypeVec { + std::vector<PbType> T; + std::string str() const; +}; + +//! Base class for all classes exposed to Python +class PbClass { + public: + PbClass(FluidSolver *parent, const std::string &name = "", PyObject *obj = NULL); + PbClass(const PbClass &a); + virtual ~PbClass(); + + // basic property setter/getters + void setName(const std::string &name) + { + mName = name; + } + std::string getName() const + { + return mName; + } + PyObject *getPyObject() const + { + return mPyObject; + } + void registerObject(PyObject *obj, PbArgs *args); + FluidSolver *getParent() const + { + return mParent; + } + void setParent(FluidSolver *v) + { + mParent = v; + } + void checkParent(); + + // hidden flag for GUI, debug output + inline bool isHidden() + { + return mHidden; + } + inline void setHidden(bool v) + { + mHidden = v; + } + + void lock(); + void unlock(); + bool tryLock(); + + // PbClass instance registry + static int getNumInstances(); + static PbClass *getInstance(int index); + static void renameObjects(); + + // converters + static bool isNullRef(PyObject *o); + static bool isNoneRef(PyObject *o); + static PbClass *createPyObject(const std::string &classname, + const std::string &name, + PbArgs &args, + PbClass *parent); + inline bool canConvertTo(const std::string &classname) + { + return Pb::canConvert(mPyObject, classname); + } + + protected: + QMutex *mMutex; + FluidSolver *mParent; + PyObject *mPyObject; + std::string mName; + bool mHidden; + + static std::vector<PbClass *> mInstances; +}; + +//!\cond Register + +void pbFinalizePlugin(FluidSolver *parent, const std::string &name, bool doTime = true); +void pbPreparePlugin(FluidSolver *parent, const std::string &name, bool doTime = true); +void pbSetError(const std::string &fn, const std::string &ex); + +//!\endcond + +} // namespace Manta + +# endif +#endif diff --git a/extern/mantaflow/helper/pwrapper/pconvert.cpp b/extern/mantaflow/helper/pwrapper/pconvert.cpp new file mode 100644 index 00000000000..c8c92cbf585 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/pconvert.cpp @@ -0,0 +1,568 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Python argument wrappers and conversion tools + * + ******************************************************************************/ + +#include "pythonInclude.h" +#include <sstream> +#include <algorithm> +#include "vectorbase.h" +#include "manta.h" + +using namespace std; + +//****************************************************************************** +// Explicit definition and instantiation of python object converters + +namespace Manta { + +extern PyTypeObject PbVec3Type; +extern PyTypeObject PbVec4Type; + +struct PbVec3 { + PyObject_HEAD float data[3]; +}; + +struct PbVec4 { + PyObject_HEAD float data[4]; +}; + +PyObject *getPyNone() +{ + Py_INCREF(Py_None); + return Py_None; +} +PyObject *incref(PyObject *obj) +{ + Py_INCREF(obj); + return obj; +} + +/*template<> PyObject* toPy<PyObject*>(PyObject* obj) { + return obj; +}*/ +template<> PyObject *toPy<int>(const int &v) +{ + return PyLong_FromLong(v); +} +/*template<> PyObject* toPy<char*>(const (char*) & val) { + return PyUnicode_DecodeLatin1(val,strlen(val),"replace"); +}*/ +template<> PyObject *toPy<string>(const string &val) +{ + return PyUnicode_DecodeLatin1(val.c_str(), val.length(), "replace"); +} +template<> PyObject *toPy<float>(const float &v) +{ + return PyFloat_FromDouble(v); +} +template<> PyObject *toPy<double>(const double &v) +{ + return PyFloat_FromDouble(v); +} +template<> PyObject *toPy<bool>(const bool &v) +{ + return PyBool_FromLong(v); +} +template<> PyObject *toPy<Vec3i>(const Vec3i &v) +{ + float x = (float)v.x, y = (float)v.y, z = (float)v.z; + return PyObject_CallFunction((PyObject *)&PbVec3Type, (char *)"fff", x, y, z); +} +template<> PyObject *toPy<Vec3>(const Vec3 &v) +{ + float x = (float)v.x, y = (float)v.y, z = (float)v.z; + return PyObject_CallFunction((PyObject *)&PbVec3Type, (char *)"fff", x, y, z); +} +template<> PyObject *toPy<Vec4i>(const Vec4i &v) +{ + float x = (float)v.x, y = (float)v.y, z = (float)v.z; + return PyObject_CallFunction((PyObject *)&PbVec4Type, (char *)"ffff", x, y, z); +} +template<> PyObject *toPy<Vec4>(const Vec4 &v) +{ + float x = (float)v.x, y = (float)v.y, z = (float)v.z; + return PyObject_CallFunction((PyObject *)&PbVec4Type, (char *)"ffff", x, y, z); +} +template<> PyObject *toPy<PbClass *>(const PbClass_Ptr &obj) +{ + return obj->getPyObject(); +} + +template<> float fromPy<float>(PyObject *obj) +{ +#if PY_MAJOR_VERSION <= 2 + if (PyInt_Check(obj)) + return PyInt_AsLong(obj); +#endif + if (PyFloat_Check(obj)) + return PyFloat_AsDouble(obj); + if (PyLong_Check(obj)) + return PyLong_AsDouble(obj); + errMsg("argument is not a float"); +} +template<> double fromPy<double>(PyObject *obj) +{ +#if PY_MAJOR_VERSION <= 2 + if (PyInt_Check(obj)) + return PyInt_AsLong(obj); +#endif + if (PyFloat_Check(obj)) + return PyFloat_AsDouble(obj); + if (PyLong_Check(obj)) + return PyLong_AsDouble(obj); + errMsg("argument is not a double"); +} +template<> PyObject *fromPy<PyObject *>(PyObject *obj) +{ + return obj; +} +template<> int fromPy<int>(PyObject *obj) +{ +#if PY_MAJOR_VERSION <= 2 + if (PyInt_Check(obj)) + return PyInt_AsLong(obj); +#endif + if (PyLong_Check(obj)) + return PyLong_AsDouble(obj); + if (PyFloat_Check(obj)) { + double a = PyFloat_AsDouble(obj); + if (fabs(a - floor(a + 0.5)) > 1e-5) + errMsg("argument is not an int"); + return (int)(a + 0.5); + } + errMsg("argument is not an int"); +} +template<> string fromPy<string>(PyObject *obj) +{ + if (PyUnicode_Check(obj)) + return PyBytes_AsString(PyUnicode_AsLatin1String(obj)); +#if PY_MAJOR_VERSION <= 2 + else if (PyString_Check(obj)) + return PyString_AsString(obj); +#endif + else + errMsg("argument is not a string"); +} +template<> const char *fromPy<const char *>(PyObject *obj) +{ + if (PyUnicode_Check(obj)) + return PyBytes_AsString(PyUnicode_AsLatin1String(obj)); +#if PY_MAJOR_VERSION <= 2 + else if (PyString_Check(obj)) + return PyString_AsString(obj); +#endif + else + errMsg("argument is not a string"); +} +template<> bool fromPy<bool>(PyObject *obj) +{ + if (!PyBool_Check(obj)) + errMsg("argument is not a boolean"); + return PyLong_AsLong(obj) != 0; +} +template<> Vec3 fromPy<Vec3>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec3Type)) { + return Vec3(((PbVec3 *)obj)->data); + } + else if (PyTuple_Check(obj) && PyTuple_Size(obj) == 3) { + return Vec3(fromPy<Real>(PyTuple_GetItem(obj, 0)), + fromPy<Real>(PyTuple_GetItem(obj, 1)), + fromPy<Real>(PyTuple_GetItem(obj, 2))); + } + errMsg("argument is not a Vec3"); +} +template<> Vec3i fromPy<Vec3i>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec3Type)) { + return toVec3iChecked(((PbVec3 *)obj)->data); + } + else if (PyTuple_Check(obj) && PyTuple_Size(obj) == 3) { + return Vec3i(fromPy<int>(PyTuple_GetItem(obj, 0)), + fromPy<int>(PyTuple_GetItem(obj, 1)), + fromPy<int>(PyTuple_GetItem(obj, 2))); + } + errMsg("argument is not a Vec3i"); +} +template<> Vec4 fromPy<Vec4>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec4Type)) { + return Vec4(((PbVec4 *)obj)->data); + } + else if (PyTuple_Check(obj) && PyTuple_Size(obj) == 4) { + return Vec4(fromPy<Real>(PyTuple_GetItem(obj, 0)), + fromPy<Real>(PyTuple_GetItem(obj, 1)), + fromPy<Real>(PyTuple_GetItem(obj, 2)), + fromPy<Real>(PyTuple_GetItem(obj, 3))); + } + errMsg("argument is not a Vec4"); +} +template<> Vec4i fromPy<Vec4i>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec4Type)) { + return toVec4i(((PbVec4 *)obj)->data); + } + else if (PyTuple_Check(obj) && PyTuple_Size(obj) == 4) { + return Vec4i(fromPy<int>(PyTuple_GetItem(obj, 0)), + fromPy<int>(PyTuple_GetItem(obj, 1)), + fromPy<int>(PyTuple_GetItem(obj, 2)), + fromPy<int>(PyTuple_GetItem(obj, 3))); + } + errMsg("argument is not a Vec4i"); +} +template<> PbType fromPy<PbType>(PyObject *obj) +{ + PbType pb = {""}; + if (!PyType_Check(obj)) + return pb; + + const char *tname = ((PyTypeObject *)obj)->tp_name; + pb.S = tname; + return pb; +} +template<> PbTypeVec fromPy<PbTypeVec>(PyObject *obj) +{ + PbTypeVec vec; + if (PyType_Check(obj)) { + vec.T.push_back(fromPy<PbType>(obj)); + } + else if (PyTuple_Check(obj)) { + int sz = PyTuple_Size(obj); + for (int i = 0; i < sz; i++) + vec.T.push_back(fromPy<PbType>(PyTuple_GetItem(obj, i))); + } + else + errMsg("argument is not a type tuple"); + return vec; +} + +template<class T> T *tmpAlloc(PyObject *obj, std::vector<void *> *tmp) +{ + if (!tmp) + throw Error("dynamic de-ref not supported for this type"); + void *ptr = malloc(sizeof(T)); + tmp->push_back(ptr); + + *((T *)ptr) = fromPy<T>(obj); + return (T *)ptr; +} +template<> float *fromPyPtr<float>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<float>(obj, tmp); +} +template<> double *fromPyPtr<double>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<double>(obj, tmp); +} +template<> int *fromPyPtr<int>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<int>(obj, tmp); +} +template<> std::string *fromPyPtr<std::string>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<std::string>(obj, tmp); +} +template<> bool *fromPyPtr<bool>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<bool>(obj, tmp); +} +template<> Vec3 *fromPyPtr<Vec3>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<Vec3>(obj, tmp); +} +template<> Vec3i *fromPyPtr<Vec3i>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<Vec3i>(obj, tmp); +} +template<> Vec4 *fromPyPtr<Vec4>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<Vec4>(obj, tmp); +} +template<> Vec4i *fromPyPtr<Vec4i>(PyObject *obj, std::vector<void *> *tmp) +{ + return tmpAlloc<Vec4i>(obj, tmp); +} + +template<> bool isPy<float>(PyObject *obj) +{ +#if PY_MAJOR_VERSION <= 2 + if (PyInt_Check(obj)) + return true; +#endif + return PyFloat_Check(obj) || PyLong_Check(obj); +} +template<> bool isPy<double>(PyObject *obj) +{ +#if PY_MAJOR_VERSION <= 2 + if (PyInt_Check(obj)) + return true; +#endif + return PyFloat_Check(obj) || PyLong_Check(obj); +} +template<> bool isPy<PyObject *>(PyObject *obj) +{ + return true; +} +template<> bool isPy<int>(PyObject *obj) +{ +#if PY_MAJOR_VERSION <= 2 + if (PyInt_Check(obj)) + return true; +#endif + if (PyLong_Check(obj)) + return true; + if (PyFloat_Check(obj)) { + double a = PyFloat_AsDouble(obj); + return fabs(a - floor(a + 0.5)) < 1e-5; + } + return false; +} +template<> bool isPy<string>(PyObject *obj) +{ + if (PyUnicode_Check(obj)) + return true; +#if PY_MAJOR_VERSION <= 2 + if (PyString_Check(obj)) + return true; +#endif + return false; +} +template<> bool isPy<const char *>(PyObject *obj) +{ + if (PyUnicode_Check(obj)) + return true; +#if PY_MAJOR_VERSION <= 2 + if (PyString_Check(obj)) + return true; +#endif + return false; +} +template<> bool isPy<bool>(PyObject *obj) +{ + return PyBool_Check(obj); +} +template<> bool isPy<Vec3>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec3Type)) + return true; + if (PyTuple_Check(obj) && PyTuple_Size(obj) == 3) { + return isPy<Real>(PyTuple_GetItem(obj, 0)) && isPy<Real>(PyTuple_GetItem(obj, 1)) && + isPy<Real>(PyTuple_GetItem(obj, 2)); + } + return false; +} +template<> bool isPy<Vec3i>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec3Type)) + return true; + if (PyTuple_Check(obj) && PyTuple_Size(obj) == 3) { + return isPy<int>(PyTuple_GetItem(obj, 0)) && isPy<int>(PyTuple_GetItem(obj, 1)) && + isPy<int>(PyTuple_GetItem(obj, 2)); + } + return false; +} +template<> bool isPy<Vec4>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec4Type)) + return true; + if (PyTuple_Check(obj) && PyTuple_Size(obj) == 4) { + return isPy<Real>(PyTuple_GetItem(obj, 0)) && isPy<Real>(PyTuple_GetItem(obj, 1)) && + isPy<Real>(PyTuple_GetItem(obj, 2)) && isPy<Real>(PyTuple_GetItem(obj, 3)); + } + return false; +} +template<> bool isPy<Vec4i>(PyObject *obj) +{ + if (PyObject_IsInstance(obj, (PyObject *)&PbVec4Type)) + return true; + if (PyTuple_Check(obj) && PyTuple_Size(obj) == 4) { + return isPy<int>(PyTuple_GetItem(obj, 0)) && isPy<int>(PyTuple_GetItem(obj, 1)) && + isPy<int>(PyTuple_GetItem(obj, 2)) && isPy<int>(PyTuple_GetItem(obj, 3)); + } + return false; +} +template<> bool isPy<PbType>(PyObject *obj) +{ + return PyType_Check(obj); +} + +//****************************************************************************** +// PbArgs class defs + +PbArgs PbArgs::EMPTY(NULL, NULL); + +PbArgs::PbArgs(PyObject *linarg, PyObject *dict) : mLinArgs(0), mKwds(0) +{ + setup(linarg, dict); +} +PbArgs::~PbArgs() +{ + for (int i = 0; i < (int)mTmpStorage.size(); i++) + free(mTmpStorage[i]); + mTmpStorage.clear(); +} + +void PbArgs::copy(PbArgs &a) +{ + mKwds = a.mKwds; + mData = a.mData; + mLinData = a.mLinData; + mLinArgs = a.mLinArgs; +} +void PbArgs::clear() +{ + mLinArgs = 0; + mKwds = 0; + mData.clear(); + mLinData.clear(); +} + +PbArgs &PbArgs::operator=(const PbArgs &a) +{ + // mLinArgs = 0; + // mKwds = 0; + return *this; +} + +void PbArgs::setup(PyObject *linarg, PyObject *dict) +{ + if (dict) { + PyObject *key, *value; + Py_ssize_t pos = 0; + while (PyDict_Next(dict, &pos, &key, &value)) { + DataElement el; + el.obj = value; + el.visited = false; + mData[fromPy<string>(key)] = el; + } + mKwds = dict; + } + if (linarg) { + size_t len = PyTuple_Size(linarg); + for (size_t i = 0; i < len; i++) { + DataElement el; + el.obj = PyTuple_GetItem(linarg, i); + el.visited = false; + mLinData.push_back(el); + } + mLinArgs = linarg; + } +} + +void PbArgs::addLinArg(PyObject *obj) +{ + DataElement el = {obj, false}; + mLinData.push_back(el); +} + +void PbArgs::check() +{ + if (has("nocheck")) + return; + + for (map<string, DataElement>::iterator it = mData.begin(); it != mData.end(); it++) { + if (!it->second.visited) + errMsg("Argument '" + it->first + "' unknown"); + } + for (size_t i = 0; i < mLinData.size(); i++) { + if (!mLinData[i].visited) { + stringstream s; + s << "Function does not read argument number #" << i; + errMsg(s.str()); + } + } +} + +FluidSolver *PbArgs::obtainParent() +{ + FluidSolver *solver = getPtrOpt<FluidSolver>("solver", -1, NULL); + if (solver != 0) + return solver; + + for (map<string, DataElement>::iterator it = mData.begin(); it != mData.end(); it++) { + PbClass *obj = Pb::objFromPy(it->second.obj); + + if (obj) { + if (solver == NULL) + solver = obj->getParent(); + } + } + for (vector<DataElement>::iterator it = mLinData.begin(); it != mLinData.end(); it++) { + PbClass *obj = Pb::objFromPy(it->obj); + + if (obj) { + if (solver == NULL) + solver = obj->getParent(); + } + } + + return solver; +} + +void PbArgs::visit(int number, const string &key) +{ + if (number >= 0 && number < (int)mLinData.size()) + mLinData[number].visited = true; + map<string, DataElement>::iterator lu = mData.find(key); + if (lu != mData.end()) + lu->second.visited = true; +} + +PyObject *PbArgs::getItem(const std::string &key, bool strict, ArgLocker *lk) +{ + map<string, DataElement>::iterator lu = mData.find(key); + if (lu == mData.end()) { + if (strict) + errMsg("Argument '" + key + "' is not defined."); + return NULL; + } + PbClass *pbo = Pb::objFromPy(lu->second.obj); + // try to lock + if (pbo && lk) + lk->add(pbo); + return lu->second.obj; +} + +PyObject *PbArgs::getItem(size_t number, bool strict, ArgLocker *lk) +{ + if (number >= mLinData.size()) { + if (!strict) + return NULL; + stringstream s; + s << "Argument number #" << number << " not specified."; + errMsg(s.str()); + } + PbClass *pbo = Pb::objFromPy(mLinData[number].obj); + // try to lock + if (pbo && lk) + lk->add(pbo); + return mLinData[number].obj; +} + +//****************************************************************************** +// ArgLocker class defs + +void ArgLocker::add(PbClass *p) +{ + if (find(locks.begin(), locks.end(), p) == locks.end()) { + locks.push_back(p); + p->lock(); + } +} +ArgLocker::~ArgLocker() +{ + for (size_t i = 0; i < locks.size(); i++) + locks[i]->unlock(); + locks.clear(); +} + +} // namespace Manta diff --git a/extern/mantaflow/helper/pwrapper/pconvert.h b/extern/mantaflow/helper/pwrapper/pconvert.h new file mode 100644 index 00000000000..9c72b8b57b9 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/pconvert.h @@ -0,0 +1,251 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Python argument wrappers and conversion tools + * + ******************************************************************************/ + +// ----------------------------------------------------------------- +// NOTE: +// Do not include this file in user code, include "manta.h" instead +// ----------------------------------------------------------------- + +#ifdef _MANTA_H +# ifndef _PCONVERT_H +# define _PCONVERT_H + +# include <string> +# include <map> +# include <vector> + +namespace Manta { +template<class T> class Grid; + +//! Locks the given PbClass Arguments until ArgLocker goes out of scope +struct ArgLocker { + void add(PbClass *p); + ~ArgLocker(); + std::vector<PbClass *> locks; +}; + +PyObject *getPyNone(); + +// for PbClass-derived classes +template<class T> T *fromPyPtr(PyObject *obj, std::vector<void *> *tmp) +{ + if (PbClass::isNullRef(obj) || PbClass::isNoneRef(obj)) + return 0; + PbClass *pbo = Pb::objFromPy(obj); + const std::string &type = Namify<T>::S; + if (!pbo || !(pbo->canConvertTo(type))) + throw Error("can't convert argument to " + type + "*"); + return (T *)(pbo); +} + +template<> float *fromPyPtr<float>(PyObject *obj, std::vector<void *> *tmp); +template<> double *fromPyPtr<double>(PyObject *obj, std::vector<void *> *tmp); +template<> int *fromPyPtr<int>(PyObject *obj, std::vector<void *> *tmp); +template<> std::string *fromPyPtr<std::string>(PyObject *obj, std::vector<void *> *tmp); +template<> bool *fromPyPtr<bool>(PyObject *obj, std::vector<void *> *tmp); +template<> Vec3 *fromPyPtr<Vec3>(PyObject *obj, std::vector<void *> *tmp); +template<> Vec3i *fromPyPtr<Vec3i>(PyObject *obj, std::vector<void *> *tmp); +template<> Vec4 *fromPyPtr<Vec4>(PyObject *obj, std::vector<void *> *tmp); +template<> Vec4i *fromPyPtr<Vec4i>(PyObject *obj, std::vector<void *> *tmp); + +PyObject *incref(PyObject *obj); +template<class T> PyObject *toPy(const T &v) +{ + PyObject *obj = v.getPyObject(); + if (obj) { + return incref(obj); + } + T *co = new T(v); + const std::string &type = Namify<typename remove_pointers<T>::type>::S; + return Pb::copyObject(co, type); +} +template<class T> bool isPy(PyObject *obj) +{ + if (PbClass::isNullRef(obj) || PbClass::isNoneRef(obj)) + return false; + PbClass *pbo = Pb::objFromPy(obj); + const std::string &type = Namify<typename remove_pointers<T>::type>::S; + return pbo && pbo->canConvertTo(type); +} + +template<class T> T fromPy(PyObject *obj) +{ + throw Error( + "Unknown type conversion. Did you pass a PbClass by value? Instead always pass " + "grids/particlesystems/etc. by reference or using a pointer."); +} + +// builtin types +template<> float fromPy<float>(PyObject *obj); +template<> double fromPy<double>(PyObject *obj); +template<> int fromPy<int>(PyObject *obj); +template<> PyObject *fromPy<PyObject *>(PyObject *obj); +template<> std::string fromPy<std::string>(PyObject *obj); +template<> const char *fromPy<const char *>(PyObject *obj); +template<> bool fromPy<bool>(PyObject *obj); +template<> Vec3 fromPy<Vec3>(PyObject *obj); +template<> Vec3i fromPy<Vec3i>(PyObject *obj); +template<> Vec4 fromPy<Vec4>(PyObject *obj); +template<> Vec4i fromPy<Vec4i>(PyObject *obj); +template<> PbType fromPy<PbType>(PyObject *obj); +template<> PbTypeVec fromPy<PbTypeVec>(PyObject *obj); + +template<> PyObject *toPy<int>(const int &v); +template<> PyObject *toPy<std::string>(const std::string &val); +template<> PyObject *toPy<float>(const float &v); +template<> PyObject *toPy<double>(const double &v); +template<> PyObject *toPy<bool>(const bool &v); +template<> PyObject *toPy<Vec3i>(const Vec3i &v); +template<> PyObject *toPy<Vec3>(const Vec3 &v); +template<> PyObject *toPy<Vec4i>(const Vec4i &v); +template<> PyObject *toPy<Vec4>(const Vec4 &v); +typedef PbClass *PbClass_Ptr; +template<> PyObject *toPy<PbClass *>(const PbClass_Ptr &obj); + +template<> bool isPy<float>(PyObject *obj); +template<> bool isPy<double>(PyObject *obj); +template<> bool isPy<int>(PyObject *obj); +template<> bool isPy<PyObject *>(PyObject *obj); +template<> bool isPy<std::string>(PyObject *obj); +template<> bool isPy<const char *>(PyObject *obj); +template<> bool isPy<bool>(PyObject *obj); +template<> bool isPy<Vec3>(PyObject *obj); +template<> bool isPy<Vec3i>(PyObject *obj); +template<> bool isPy<Vec4>(PyObject *obj); +template<> bool isPy<Vec4i>(PyObject *obj); +template<> bool isPy<PbType>(PyObject *obj); + +//! Encapsulation of python arguments +class PbArgs { + public: + PbArgs(PyObject *linargs = NULL, PyObject *dict = NULL); + ~PbArgs(); + void setup(PyObject *linargs = NULL, PyObject *dict = NULL); + + void check(); + FluidSolver *obtainParent(); + + inline int numLinArgs() + { + return mLinData.size(); + } + + inline bool has(const std::string &key) + { + return getItem(key, false) != NULL; + } + inline void deleteItem(const std::string &key) + { + if (mData.find(key) != mData.end()) + mData.erase(mData.find(key)); + } + + inline PyObject *linArgs() + { + return mLinArgs; + } + inline PyObject *kwds() + { + return mKwds; + } + + void addLinArg(PyObject *obj); + + template<class T> inline void add(const std::string &key, T arg) + { + DataElement el = {toPy(arg), false}; + mData[key] = el; + } + template<class T> inline T get(const std::string &key, int number = -1, ArgLocker *lk = NULL) + { + visit(number, key); + PyObject *o = getItem(key, false, lk); + if (o) + return fromPy<T>(o); + o = getItem(number, false, lk); + if (o) + return fromPy<T>(o); + errMsg("Argument '" + key + "' is not defined."); + } + template<class T> + inline T getOpt(const std::string &key, int number, T defarg, ArgLocker *lk = NULL) + { + visit(number, key); + PyObject *o = getItem(key, false, lk); + if (o) + return fromPy<T>(o); + if (number >= 0) + o = getItem(number, false, lk); + return (o) ? fromPy<T>(o) : defarg; + } + template<class T> + inline T *getPtrOpt(const std::string &key, int number, T *defarg, ArgLocker *lk = NULL) + { + visit(number, key); + PyObject *o = getItem(key, false, lk); + if (o) + return fromPyPtr<T>(o, &mTmpStorage); + if (number >= 0) + o = getItem(number, false, lk); + return o ? fromPyPtr<T>(o, &mTmpStorage) : defarg; + } + template<class T> inline T *getPtr(const std::string &key, int number = -1, ArgLocker *lk = NULL) + { + visit(number, key); + PyObject *o = getItem(key, false, lk); + if (o) + return fromPyPtr<T>(o, &mTmpStorage); + o = getItem(number, false, lk); + if (o) + return fromPyPtr<T>(o, &mTmpStorage); + errMsg("Argument '" + key + "' is not defined."); + } + + // automatic template type deduction + template<class T> bool typeCheck(int num, const std::string &name) + { + PyObject *o = getItem(name, false, 0); + if (!o) + o = getItem(num, false, 0); + return o ? isPy<typename remove_pointers<T>::type>(o) : false; + } + + PbArgs &operator=(const PbArgs &a); // dummy + void copy(PbArgs &a); + void clear(); + void visit(int num, const std::string &key); + + static PbArgs EMPTY; + + protected: + PyObject *getItem(const std::string &key, bool strict, ArgLocker *lk = NULL); + PyObject *getItem(size_t number, bool strict, ArgLocker *lk = NULL); + + struct DataElement { + PyObject *obj; + bool visited; + }; + std::map<std::string, DataElement> mData; + std::vector<DataElement> mLinData; + PyObject *mLinArgs, *mKwds; + std::vector<void *> mTmpStorage; +}; + +} // namespace Manta + +# if NUMPY == 1 +# include "numpyWrap.h" +# endif + +# endif +#endif diff --git a/extern/mantaflow/helper/pwrapper/pvec3.cpp b/extern/mantaflow/helper/pwrapper/pvec3.cpp new file mode 100644 index 00000000000..69bde2a2ad0 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/pvec3.cpp @@ -0,0 +1,414 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Vec3 class extension for python + * + ******************************************************************************/ + +#include "pythonInclude.h" +#include <string> +#include <sstream> +#include "vectorbase.h" +#include "structmember.h" +#include "manta.h" + +using namespace std; + +namespace Manta { + +extern PyTypeObject PbVec3Type; + +struct PbVec3 { + PyObject_HEAD float data[3]; +}; + +static void PbVec3Dealloc(PbVec3 *self) +{ + Py_TYPE(self)->tp_free((PyObject *)self); +} + +static PyObject *PbVec3New(PyTypeObject *type, PyObject *args, PyObject *kwds) +{ + return type->tp_alloc(type, 0); +} + +static int PbVec3Init(PbVec3 *self, PyObject *args, PyObject *kwds) +{ + + float x1 = numeric_limits<float>::quiet_NaN(), x2 = x1, x3 = x1; + if (!PyArg_ParseTuple(args, "|fff", &x1, &x2, &x3)) + return -1; + + if (!c_isnan(x1)) { + self->data[0] = x1; + if (!c_isnan(x2) && !c_isnan(x3)) { + self->data[1] = x2; + self->data[2] = x3; + } + else { + if (!c_isnan(x2) || !c_isnan(x3)) { + errMsg("Invalid partial init of vec3"); + } + self->data[1] = x1; + self->data[2] = x1; + } + } + else { + self->data[0] = 0; + self->data[1] = 0; + self->data[2] = 0; + } + return 0; +} + +static PyObject *PbVec3Repr(PbVec3 *self) +{ + Manta::Vec3 v(self->data[0], self->data[1], self->data[2]); + return PyUnicode_FromFormat(v.toString().c_str()); +} + +static PyMemberDef PbVec3Members[] = { + {(char *)"x", T_FLOAT, offsetof(PbVec3, data), 0, (char *)"X"}, + {(char *)"y", T_FLOAT, offsetof(PbVec3, data) + sizeof(float), 0, (char *)"Y"}, + {(char *)"z", T_FLOAT, offsetof(PbVec3, data) + sizeof(float) * 2, 0, (char *)"Z"}, + {NULL} // Sentinel +}; + +static PyMethodDef PbVec3Methods[] = { + //{"name", (PyCFunction)Noddy_name, METH_NOARGS, "Return the name, combining the first and last + //name" }, + {NULL} // Sentinel +}; + +// operator overloads + +inline PyObject *PbNew(const Vec3 &a) +{ + PbVec3 *obj = (PbVec3 *)PbVec3New(&PbVec3Type, 0, 0); + obj->data[0] = a.x; + obj->data[1] = a.y; + obj->data[2] = a.z; + return (PyObject *)obj; +} + +#define CONVERTVEC(obj) \ + Vec3 v##obj; \ + if (PyObject_TypeCheck(obj, &PbVec3Type)) \ + v##obj = Vec3(&(((PbVec3 *)obj)->data[0])); \ + else if (PyFloat_Check(obj)) \ + v##obj = Vec3(PyFloat_AsDouble(obj)); \ + else if (PyLong_Check(obj)) \ + v##obj = Vec3(PyLong_AsDouble(obj)); \ + else { \ + Py_INCREF(Py_NotImplemented); \ + return Py_NotImplemented; \ + } + +#define OPHEADER \ + if (!PyObject_TypeCheck(a, &PbVec3Type) && !PyObject_TypeCheck(b, &PbVec3Type)) { \ + Py_INCREF(Py_NotImplemented); \ + return Py_NotImplemented; \ + } \ + CONVERTVEC(a) \ + CONVERTVEC(b) + +#define OPHEADER1 \ + if (!PyObject_TypeCheck(a, &PbVec3Type)) { \ + Py_INCREF(Py_NotImplemented); \ + return Py_NotImplemented; \ + } \ + CONVERTVEC(a) + +PyObject *PbVec3Add(PyObject *a, PyObject *b) +{ + OPHEADER + return PbNew(va + vb); +} + +PyObject *PbVec3Sub(PyObject *a, PyObject *b) +{ + OPHEADER + return PbNew(va - vb); +} + +PyObject *PbVec3Mult(PyObject *a, PyObject *b) +{ + OPHEADER + return PbNew(va * vb); +} + +PyObject *PbVec3Div(PyObject *a, PyObject *b) +{ + OPHEADER + return PbNew(va / vb); +} + +PyObject *PbVec3Negative(PyObject *a) +{ + OPHEADER1 + return PbNew(-va); +} + +// numbers are defined subtely different in Py3 (WTF?) +#if PY_MAJOR_VERSION >= 3 +static PyNumberMethods PbVec3NumberMethods = { + (binaryfunc)PbVec3Add, // binaryfunc nb_add; + (binaryfunc)PbVec3Sub, // binaryfunc nb_sub; + (binaryfunc)PbVec3Mult, // binaryfunc nb_mult; + 0, // binaryfunc nb_remainder; + 0, // binaryfunc nb_divmod; + 0, // ternaryfunc nb_power; + (unaryfunc)PbVec3Negative, // unaryfunc nb_negative; + 0, // unaryfunc nb_positive; + 0, // unaryfunc nb_absolute; + 0, // inquiry nb_bool; + 0, // unaryfunc nb_invert; + 0, // binaryfunc nb_lshift; + 0, // binaryfunc nb_rshift; + 0, // binaryfunc nb_and; + 0, // binaryfunc nb_xor; + 0, // binaryfunc nb_or; + 0, // unaryfunc nb_int; + 0, // void *nb_reserved; + 0, // unaryfunc nb_float; + 0, // binaryfunc nb_inplace_add; + 0, // binaryfunc nb_inplace_subtract; + 0, // binaryfunc nb_inplace_multiply; + 0, // binaryfunc nb_inplace_remainder; + 0, // ternaryfunc nb_inplace_power; + 0, // binaryfunc nb_inplace_lshift; + 0, // binaryfunc nb_inplace_rshift; + 0, // binaryfunc nb_inplace_and; + 0, // binaryfunc nb_inplace_xor; + 0, // binaryfunc nb_inplace_or; + + 0, // binaryfunc nb_floor_divide; + (binaryfunc)PbVec3Div, // binaryfunc nb_true_divide; + 0, // binaryfunc nb_inplace_floor_divide; + 0, // binaryfunc nb_inplace_true_divide; + + 0 // unaryfunc nb_index; +}; +#else +static PyNumberMethods PbVec3NumberMethods = { + (binaryfunc)PbVec3Add, // binaryfunc nb_add; + (binaryfunc)PbVec3Sub, // binaryfunc nb_sub; + (binaryfunc)PbVec3Mult, // binaryfunc nb_mult; + 0, // binaryfunc nb_divide; + 0, // binaryfunc nb_remainder; + 0, // binaryfunc nb_divmod; + 0, // ternaryfunc nb_power; + (unaryfunc)PbVec3Negative, // unaryfunc nb_negative; + 0, // unaryfunc nb_positive; + 0, // unaryfunc nb_absolute; + 0, // inquiry nb_nonzero; + 0, // unaryfunc nb_invert; + 0, // binaryfunc nb_lshift; + 0, // binaryfunc nb_rshift; + 0, // binaryfunc nb_and; + 0, // binaryfunc nb_xor; + 0, // binaryfunc nb_or; + 0, // coercion nb_coerce; + 0, // unaryfunc nb_int; + 0, // unaryfunc nb_long; + 0, // unaryfunc nb_float; + 0, // unaryfunc nb_oct; + 0, // unaryfunc nb_hex; + 0, // binaryfunc nb_inplace_add; + 0, // binaryfunc nb_inplace_subtract; + 0, // binaryfunc nb_inplace_multiply; + 0, // binaryfunc nb_inplace_divide; + 0, // binaryfunc nb_inplace_remainder; + 0, // ternaryfunc nb_inplace_power; + 0, // binaryfunc nb_inplace_lshift; + 0, // binaryfunc nb_inplace_rshift; + 0, // binaryfunc nb_inplace_and; + 0, // binaryfunc nb_inplace_xor; + 0, // binaryfunc nb_inplace_or; + 0, // binaryfunc nb_floor_divide; + (binaryfunc)PbVec3Div, // binaryfunc nb_true_divide; + 0, // binaryfunc nb_inplace_floor_divide; + 0, // binaryfunc nb_inplace_true_divide; + 0, // unaryfunc nb_index; +}; +#endif + +PyTypeObject PbVec3Type = { + PyVarObject_HEAD_INIT(NULL, 0) "manta.vec3", /* tp_name */ + sizeof(PbVec3), /* tp_basicsize */ + 0, /* tp_itemsize */ + (destructor)PbVec3Dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ + 0, /* tp_reserved */ + (reprfunc)PbVec3Repr, /* tp_repr */ + &PbVec3NumberMethods, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ +#if PY_MAJOR_VERSION >= 3 + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags */ +#else + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE | Py_TPFLAGS_CHECKTYPES, /* tp_flags */ +#endif + "float vector type", /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + PbVec3Methods, /* tp_methods */ + PbVec3Members, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + (initproc)PbVec3Init, /* tp_init */ + 0, /* tp_alloc */ + PbVec3New, /* tp_new */ +}; + +inline PyObject *castPy(PyTypeObject *p) +{ + return reinterpret_cast<PyObject *>(static_cast<void *>(p)); +} + +// 4d vector + +extern PyTypeObject PbVec4Type; + +struct PbVec4 { + PyObject_HEAD float data[4]; +}; + +static PyMethodDef PbVec4Methods[] = { + {NULL} // Sentinel +}; + +static PyMemberDef PbVec4Members[] = { + {(char *)"x", T_FLOAT, offsetof(PbVec4, data), 0, (char *)"X"}, + {(char *)"y", T_FLOAT, offsetof(PbVec4, data) + sizeof(float) * 1, 0, (char *)"Y"}, + {(char *)"z", T_FLOAT, offsetof(PbVec4, data) + sizeof(float) * 2, 0, (char *)"Z"}, + {(char *)"t", T_FLOAT, offsetof(PbVec4, data) + sizeof(float) * 3, 0, (char *)"T"}, + {NULL} // Sentinel +}; + +static void PbVec4Dealloc(PbVec4 *self) +{ + Py_TYPE(self)->tp_free((PyObject *)self); +} + +static PyObject *PbVec4New(PyTypeObject *type, PyObject *args, PyObject *kwds) +{ + return type->tp_alloc(type, 0); +} + +static int PbVec4Init(PbVec4 *self, PyObject *args, PyObject *kwds) +{ + + float x1 = numeric_limits<float>::quiet_NaN(), x2 = x1, x3 = x1, x4 = x1; + if (!PyArg_ParseTuple(args, "|ffff", &x1, &x2, &x3, &x4)) + return -1; + + if (!c_isnan(x1)) { + self->data[0] = x1; + if (!c_isnan(x2) && !c_isnan(x3) && !c_isnan(x4)) { + self->data[1] = x2; + self->data[2] = x3; + self->data[3] = x4; + } + else { + if (!c_isnan(x2) || !c_isnan(x3) || !c_isnan(x4)) { + errMsg("Invalid partial init of vec4"); + } + self->data[1] = self->data[2] = self->data[3] = x1; + } + } + else { + self->data[0] = self->data[1] = self->data[2] = self->data[3] = 0; + } + return 0; +} + +static PyObject *PbVec4Repr(PbVec4 *self) +{ + Manta::Vec4 v(self->data[0], self->data[1], self->data[2], self->data[3]); + return PyUnicode_FromFormat(v.toString().c_str()); +} + +PyTypeObject PbVec4Type = { + PyVarObject_HEAD_INIT(NULL, 0) "manta.vec4", /* tp_name */ + sizeof(PbVec4), /* tp_basicsize */ + 0, /* tp_itemsize */ + (destructor)PbVec4Dealloc, /* tp_dealloc */ + 0, /* tp_print */ + 0, /* tp_getattr */ + 0, /* tp_setattr */ + 0, /* tp_reserved */ + (reprfunc)PbVec4Repr, /* tp_repr */ + NULL, // &PbVec4NumberMethods, /* tp_as_number */ + 0, /* tp_as_sequence */ + 0, /* tp_as_mapping */ + 0, /* tp_hash */ + 0, /* tp_call */ + 0, /* tp_str */ + 0, /* tp_getattro */ + 0, /* tp_setattro */ + 0, /* tp_as_buffer */ +#if PY_MAJOR_VERSION >= 3 + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags */ +#else + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE | Py_TPFLAGS_CHECKTYPES, /* tp_flags */ +#endif + "float vector type", /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + PbVec4Methods, /* tp_methods */ + PbVec4Members, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + (initproc)PbVec4Init, /* tp_init */ + 0, /* tp_alloc */ + PbVec4New, /* tp_new */ +}; + +// register + +void PbVecInitialize(PyObject *module) +{ + if (PyType_Ready(&PbVec3Type) < 0) + errMsg("can't initialize Vec3 type"); + Py_INCREF(castPy(&PbVec3Type)); + PyModule_AddObject(module, "vec3", (PyObject *)&PbVec3Type); + + if (PyType_Ready(&PbVec4Type) < 0) + errMsg("can't initialize Vec4 type"); + Py_INCREF(castPy(&PbVec4Type)); + PyModule_AddObject(module, "vec4", (PyObject *)&PbVec4Type); +} +const static Pb::Register _REG(PbVecInitialize); + +} // namespace Manta diff --git a/extern/mantaflow/helper/pwrapper/pythonInclude.h b/extern/mantaflow/helper/pwrapper/pythonInclude.h new file mode 100644 index 00000000000..0f78c6641d2 --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/pythonInclude.h @@ -0,0 +1,48 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Base class for particle systems + * + ******************************************************************************/ + +#ifndef _PYTHONINCLUDE_H +#define _PYTHONINCLUDE_H + +#if defined(WIN32) || defined(_WIN32) + +// note - we have to include these first! +# include <string> +# include <vector> +# include <iostream> + +#endif + +// the PYTHON_DEBUG_WITH_RELEASE define enables linking with python debug libraries +#if (defined(_DEBUG) || (DEBUG == 1)) && defined(DEBUG_PYTHON_WITH_RELEASE) + +// special handling, disable linking with debug version of python libs +# undef _DEBUG +# define NDEBUG +# include <Python.h> +# if NUMPY == 1 +# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +# include "numpy/arrayobject.h" +# endif +# define _DEBUG +# undef NDEBUG + +#else +# include <Python.h> +# if NUMPY == 1 +# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +# include "numpy/arrayobject.h" +# endif +#endif + +#endif diff --git a/extern/mantaflow/helper/pwrapper/registry.cpp b/extern/mantaflow/helper/pwrapper/registry.cpp new file mode 100644 index 00000000000..332b9e158ac --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/registry.cpp @@ -0,0 +1,784 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011-2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Auto python registry + * + ******************************************************************************/ + +#include <string.h> +#include "pythonInclude.h" +#include "structmember.h" +#include "manta.h" + +using namespace std; + +const string gDefaultModuleName = "manta"; + +namespace Pb { + +//****************************************************************************** +// Custom object definition + +struct Method { + Method(const string &n, const string &d, GenericFunction f) : name(n), doc(d), func(f) + { + } + string name, doc; + GenericFunction func; + + PyMethodDef def() + { + PyMethodDef def = {&name[0], (PyCFunction)func, METH_VARARGS | METH_KEYWORDS, &doc[0]}; + return def; + } +}; +struct GetSet { + GetSet() : getter(0), setter(0) + { + } + GetSet(const string &n, const string &d, Getter g, Setter s) + : name(n), doc(d), getter(g), setter(s) + { + } + string name, doc; + Getter getter; + Setter setter; + + PyGetSetDef def() + { + PyGetSetDef def = {&name[0], getter, setter, &doc[0], NULL}; + return def; + } +}; + +struct ClassData { + string cName, pyName; + string cPureName, cTemplate; + InitFunc init; + PyTypeObject typeInfo; + PyNumberMethods numInfo; + // PySequenceMethods seqInfo; + vector<Method> methods; + map<string, GetSet> getset; + map<string, OperatorFunction> ops; + ClassData *baseclass; + string baseclassName; + Constructor constructor; + + vector<PyMethodDef> genMethods; + vector<PyGetSetDef> genGetSet; +}; + +struct PbObject { + PyObject_HEAD Manta::PbClass *instance; + ClassData *classdef; +}; + +//****************************************************** +// Internal wrapper class + +//! Registers all classes and methods exposed to Python. +/*! This class is only used internally by Pb:: framwork. + * Please use the functionality of PbClass to lookup and translate pointers. */ +class WrapperRegistry { + public: + static WrapperRegistry &instance(); + void addClass(const std::string &name, + const std::string &internalName, + const std::string &baseclass); + void addEnumEntry(const std::string &name, int value); + void addExternalInitializer(InitFunc func); + void addMethod(const std::string &classname, + const std::string &methodname, + GenericFunction method); + void addOperator(const std::string &classname, + const std::string &methodname, + OperatorFunction method); + void addConstructor(const std::string &classname, Constructor method); + void addGetSet(const std::string &classname, + const std::string &property, + Getter getfunc, + Setter setfunc); + void addPythonPath(const std::string &path); + void addPythonCode(const std::string &file, const std::string &code); + PyObject *createPyObject(const std::string &classname, + const std::string &name, + Manta::PbArgs &args, + Manta::PbClass *parent); + void construct(const std::string &scriptname, const vector<string> &args); + void cleanup(); + void renameObjects(); + void runPreInit(); + PyObject *initModule(); + ClassData *lookup(const std::string &name); + bool canConvert(ClassData *from, ClassData *to); + + private: + ClassData *getOrConstructClass(const string &name); + void registerBaseclasses(); + void registerDummyTypes(); + void registerMeta(); + void addConstants(PyObject *module); + void registerOperators(ClassData *cls); + void addParentMethods(ClassData *cls, ClassData *base); + WrapperRegistry(); + std::map<std::string, ClassData *> mClasses; + std::vector<ClassData *> mClassList; + std::vector<InitFunc> mExtInitializers; + std::vector<std::string> mPaths; + std::string mCode, mScriptName; + std::vector<std::string> args; + std::map<std::string, int> mEnumValues; +}; + +//****************************************************************************** +// Callback functions + +PyObject *cbGetClass(PbObject *self, void *cl) +{ + return Manta::toPy(self->classdef->cPureName); +} + +PyObject *cbGetTemplate(PbObject *self, void *cl) +{ + return Manta::toPy(self->classdef->cTemplate); +} + +PyObject *cbGetCName(PbObject *self, void *cl) +{ + return Manta::toPy(self->classdef->cName); +} + +void cbDealloc(PbObject *self) +{ + // cout << "dealloc " << self->instance->getName() << " " << self->classdef->cName << endl; + if (self->instance) { +#ifndef BLENDER + // don't delete top-level objects + if (self->instance->getParent() != self->instance) + delete self->instance; +#else + // in Blender we *have* to delete all objects + delete self->instance; +#endif + } + Py_TYPE(self)->tp_free((PyObject *)self); +} + +PyObject *cbNew(PyTypeObject *type, PyObject *args, PyObject *kwds) +{ + PbObject *self = (PbObject *)type->tp_alloc(type, 0); + if (self != NULL) { + // lookup and link classdef + self->classdef = WrapperRegistry::instance().lookup(type->tp_name); + self->instance = NULL; + // cout << "creating " << self->classdef->cName << endl; + } + else + errMsg("can't allocate new python class object"); + return (PyObject *)self; +} + +int cbDisableConstructor(PyObject *self, PyObject *args, PyObject *kwds) +{ + errMsg("Can't instantiate a class template without template arguments"); + return -1; +} + +PyMODINIT_FUNC PyInit_Main(void) +{ + MantaEnsureRegistration(); +#if PY_MAJOR_VERSION >= 3 + return WrapperRegistry::instance().initModule(); +#else + WrapperRegistry::instance().initModule(); +#endif +} + +//****************************************************** +// WrapperRegistry + +WrapperRegistry::WrapperRegistry() +{ + addClass("__modclass__", "__modclass__", ""); + addClass("PbClass", "PbClass", ""); +} + +ClassData *WrapperRegistry::getOrConstructClass(const string &classname) +{ + map<string, ClassData *>::iterator it = mClasses.find(classname); + + if (it != mClasses.end()) + return it->second; + ClassData *data = new ClassData; + data->cName = classname; + data->cPureName = classname; + data->cTemplate = ""; + size_t tplIdx = classname.find('<'); + if (tplIdx != string::npos) { + data->cPureName = classname.substr(0, tplIdx); + data->cTemplate = classname.substr(tplIdx + 1, classname.find('>') - tplIdx - 1); + } + data->baseclass = NULL; + data->constructor = cbDisableConstructor; + mClasses[classname] = data; + mClassList.push_back(data); + return data; +} + +void replaceAll(string &source, string const &find, string const &replace) +{ + for (string::size_type i = 0; (i = source.find(find, i)) != std::string::npos;) { + source.replace(i, find.length(), replace); + i += replace.length() - find.length() + 1; + } +} + +void WrapperRegistry::addClass(const string &pyName, + const string &internalName, + const string &baseclass) +{ + ClassData *data = getOrConstructClass(internalName); + + // regularize python name + string pythonName = pyName; + replaceAll(pythonName, "<", "_"); + replaceAll(pythonName, ">", ""); + replaceAll(pythonName, ",", "_"); + + if (data->pyName.empty()) + data->pyName = pythonName; + mClasses[pythonName] = data; + if (!baseclass.empty()) + data->baseclassName = baseclass; +} + +void WrapperRegistry::addEnumEntry(const string &name, int value) +{ + /// Gather static definitions to add them as static python objects afterwards + if (mEnumValues.insert(std::make_pair(name, value)).second == false) { + errMsg("Enum entry '" + name + "' already existing..."); + } +} + +void WrapperRegistry::addExternalInitializer(InitFunc func) +{ + mExtInitializers.push_back(func); +} + +void WrapperRegistry::addPythonPath(const string &path) +{ + mPaths.push_back(path); +} + +void WrapperRegistry::addPythonCode(const string &file, const string &code) +{ + mCode += code + "\n"; +} + +void WrapperRegistry::addGetSet(const string &classname, + const string &property, + Getter getfunc, + Setter setfunc) +{ + ClassData *classdef = getOrConstructClass(classname); + GetSet &def = classdef->getset[property]; + if (def.name.empty()) { + def.name = property; + def.doc = property; + } + if (getfunc) + def.getter = getfunc; + if (setfunc) + def.setter = setfunc; +} + +void WrapperRegistry::addMethod(const string &classname, + const string &methodname, + GenericFunction func) +{ + string aclass = classname; + if (aclass.empty()) + aclass = "__modclass__"; + + ClassData *classdef = getOrConstructClass(aclass); + for (int i = 0; i < (int)classdef->methods.size(); i++) + if (classdef->methods[i].name == methodname) + return; // avoid duplicates + classdef->methods.push_back(Method(methodname, methodname, func)); +} + +void WrapperRegistry::addOperator(const string &classname, + const string &methodname, + OperatorFunction func) +{ + if (classname.empty()) + errMsg("PYTHON operators have to be defined within classes."); + + string op = methodname.substr(8); + ClassData *classdef = getOrConstructClass(classname); + classdef->ops[op] = func; +} + +void WrapperRegistry::addConstructor(const string &classname, Constructor func) +{ + ClassData *classdef = getOrConstructClass(classname); + classdef->constructor = func; +} + +void WrapperRegistry::addParentMethods(ClassData *cur, ClassData *base) +{ + if (base == 0) + return; + + for (vector<Method>::iterator it = base->methods.begin(); it != base->methods.end(); ++it) + addMethod(cur->cName, it->name, it->func); + + for (map<string, GetSet>::iterator it = base->getset.begin(); it != base->getset.end(); ++it) + addGetSet(cur->cName, it->first, it->second.getter, it->second.setter); + + for (map<string, OperatorFunction>::iterator it = base->ops.begin(); it != base->ops.end(); ++it) + cur->ops[it->first] = it->second; + + addParentMethods(cur, base->baseclass); +} + +void WrapperRegistry::registerBaseclasses() +{ + for (int i = 0; i < (int)mClassList.size(); i++) { + string bname = mClassList[i]->baseclassName; + if (!bname.empty()) { + mClassList[i]->baseclass = lookup(bname); + if (!mClassList[i]->baseclass) + errMsg("Registering class '" + mClassList[i]->cName + "' : Base class '" + bname + + "' not found"); + } + } + + for (int i = 0; i < (int)mClassList.size(); i++) { + addParentMethods(mClassList[i], mClassList[i]->baseclass); + } +} + +void WrapperRegistry::registerMeta() +{ + for (int i = 0; i < (int)mClassList.size(); i++) { + mClassList[i]->getset["_class"] = GetSet("_class", "C class name", (Getter)cbGetClass, 0); + mClassList[i]->getset["_cname"] = GetSet("_cname", "Full C name", (Getter)cbGetCName, 0); + mClassList[i]->getset["_T"] = GetSet("_T", "C template argument", (Getter)cbGetTemplate, 0); + } +} + +void WrapperRegistry::registerOperators(ClassData *cls) +{ + PyNumberMethods &num = cls->numInfo; + for (map<string, OperatorFunction>::iterator it = cls->ops.begin(); it != cls->ops.end(); it++) { + const string &op = it->first; + OperatorFunction func = it->second; + if (op == "+=") + num.nb_inplace_add = func; + else if (op == "-=") + num.nb_inplace_subtract = func; + else if (op == "*=") + num.nb_inplace_multiply = func; + else if (op == "+") + num.nb_add = func; + else if (op == "-") + num.nb_subtract = func; + else if (op == "*") + num.nb_multiply = func; +#if PY_MAJOR_VERSION < 3 + else if (op == "/=") + num.nb_inplace_divide = func; + else if (op == "/") + num.nb_divide = func; +#else + else if (op == "/=") + num.nb_inplace_true_divide = func; + else if (op == "/") + num.nb_true_divide = func; +#endif + else + errMsg("PYTHON operator " + op + " not supported"); + } +} + +void WrapperRegistry::registerDummyTypes() +{ + vector<string> add; + for (vector<ClassData *>::iterator it = mClassList.begin(); it != mClassList.end(); ++it) { + string cName = (*it)->cName; + if (cName.find('<') != string::npos) + add.push_back(cName.substr(0, cName.find('<'))); + } + for (int i = 0; i < (int)add.size(); i++) + addClass(add[i], add[i], ""); +} + +ClassData *WrapperRegistry::lookup(const string &name) +{ + for (map<string, ClassData *>::iterator it = mClasses.begin(); it != mClasses.end(); ++it) { + if (it->first == name || it->second->cName == name) + return it->second; + } + return NULL; +} + +void WrapperRegistry::cleanup() +{ + for (vector<ClassData *>::iterator it = mClassList.begin(); it != mClassList.end(); ++it) { + delete *it; + } + mClasses.clear(); + mClassList.clear(); +} + +WrapperRegistry &WrapperRegistry::instance() +{ + static WrapperRegistry inst; + return inst; +} + +bool WrapperRegistry::canConvert(ClassData *from, ClassData *to) +{ + if (from == to) + return true; + if (from->baseclass) + return canConvert(from->baseclass, to); + return false; +} + +void WrapperRegistry::addConstants(PyObject *module) +{ + // expose arguments + PyObject *list = PyList_New(args.size()); + for (int i = 0; i < (int)args.size(); i++) + PyList_SET_ITEM(list, i, Manta::toPy(args[i])); + PyModule_AddObject(module, "args", list); + PyModule_AddObject(module, "SCENEFILE", Manta::toPy(mScriptName)); + + // expose compile flags +#ifdef DEBUG + PyModule_AddObject(module, "DEBUG", Manta::toPy<bool>(true)); +#else + PyModule_AddObject(module, "DEBUG", Manta::toPy<bool>(false)); +#endif +#ifdef MANTA_MT + PyModule_AddObject(module, "MT", Manta::toPy<bool>(true)); +#else + PyModule_AddObject(module, "MT", Manta::toPy<bool>(false)); +#endif +#ifdef GUI + PyModule_AddObject(module, "GUI", Manta::toPy<bool>(true)); +#else + PyModule_AddObject(module, "GUI", Manta::toPy<bool>(false)); +#endif +#if FLOATINGPOINT_PRECISION == 2 + PyModule_AddObject(module, "DOUBLEPRECISION", Manta::toPy<bool>(true)); +#else + PyModule_AddObject(module, "DOUBLEPRECISION", Manta::toPy<bool>(false)); +#endif + // cuda off for now + PyModule_AddObject(module, "CUDA", Manta::toPy<bool>(false)); + + // expose enum entries + std::map<std::string, int>::iterator it; + for (it = mEnumValues.begin(); it != mEnumValues.end(); it++) { + PyModule_AddObject(module, it->first.c_str(), Manta::toPy(it->second)); + // Alternative would be: + // e.g. PyModule_AddIntConstant(module, "FlagFluid", 1); + } +} + +void WrapperRegistry::runPreInit() +{ + // add python directories to path + PyObject *sys_path = PySys_GetObject((char *)"path"); + for (size_t i = 0; i < mPaths.size(); i++) { + PyObject *path = Manta::toPy(mPaths[i]); + if (sys_path == NULL || path == NULL || PyList_Append(sys_path, path) < 0) { + errMsg("unable to set python path"); + } + Py_DECREF(path); + } + if (!mCode.empty()) { + mCode = "from manta import *\n" + mCode; + PyRun_SimpleString(mCode.c_str()); + } +} + +PyObject *WrapperRegistry::createPyObject(const string &classname, + const string &name, + Manta::PbArgs &args, + Manta::PbClass *parent) +{ + ClassData *classdef = lookup(classname); + if (!classdef) + errMsg("Class " + classname + " doesn't exist."); + + // create object + PyObject *obj = cbNew(&classdef->typeInfo, NULL, NULL); + PbObject *self = (PbObject *)obj; + PyObject *nkw = 0; + + if (args.kwds()) + nkw = PyDict_Copy(args.kwds()); + else + nkw = PyDict_New(); + + PyObject *nocheck = Py_BuildValue("s", "yes"); + PyDict_SetItemString(nkw, "nocheck", nocheck); + if (parent) + PyDict_SetItemString(nkw, "parent", parent->getPyObject()); + + // create instance + if (self->classdef->constructor(obj, args.linArgs(), nkw) < 0) + errMsg("error raised in constructor"); // assume condition is already set + + Py_DECREF(nkw); + Py_DECREF(nocheck); + self->instance->setName(name); + + return obj; +} + +// prepare typeinfo and register python module +void WrapperRegistry::construct(const string &scriptname, const vector<string> &args) +{ + mScriptName = scriptname; + this->args = args; + + registerBaseclasses(); + registerMeta(); + registerDummyTypes(); + + // work around for certain gcc versions, cast to char* + PyImport_AppendInittab((char *)gDefaultModuleName.c_str(), PyInit_Main); +} + +inline PyObject *castPy(PyTypeObject *p) +{ + return reinterpret_cast<PyObject *>(static_cast<void *>(p)); +} + +PyObject *WrapperRegistry::initModule() +{ + // generate and terminate all method lists + PyMethodDef sentinelFunc = {NULL, NULL, 0, NULL}; + PyGetSetDef sentinelGetSet = {NULL, NULL, NULL, NULL, NULL}; + for (int i = 0; i < (int)mClassList.size(); i++) { + ClassData *cls = mClassList[i]; + cls->genMethods.clear(); + cls->genGetSet.clear(); + for (vector<Method>::iterator i2 = cls->methods.begin(); i2 != cls->methods.end(); ++i2) + cls->genMethods.push_back(i2->def()); + for (map<string, GetSet>::iterator i2 = cls->getset.begin(); i2 != cls->getset.end(); ++i2) + cls->genGetSet.push_back(i2->second.def()); + + cls->genMethods.push_back(sentinelFunc); + cls->genGetSet.push_back(sentinelGetSet); + } + + // prepare module info +#if PY_MAJOR_VERSION >= 3 + static PyModuleDef MainModule = {PyModuleDef_HEAD_INIT, + gDefaultModuleName.c_str(), + "Bridge module to the C++ solver", + -1, + NULL, + NULL, + NULL, + NULL, + NULL}; + // get generic methods (plugin functions) + MainModule.m_methods = &mClasses["__modclass__"]->genMethods[0]; + + // create module + PyObject *module = PyModule_Create(&MainModule); +#else + PyObject *module = Py_InitModule(gDefaultModuleName.c_str(), + &mClasses["__modclass__"]->genMethods[0]); +#endif + if (module == NULL) + return NULL; + + // load classes + for (vector<ClassData *>::iterator it = mClassList.begin(); it != mClassList.end(); ++it) { + ClassData &data = **it; + char *nameptr = (char *)data.pyName.c_str(); + + // define numeric substruct + PyNumberMethods *num = 0; + if (!data.ops.empty()) { + num = &data.numInfo; + memset(num, 0, sizeof(PyNumberMethods)); + registerOperators(&data); + } + + // define python classinfo + PyTypeObject t = { + PyVarObject_HEAD_INIT(NULL, 0)(char *) data.pyName.c_str(), // tp_name + sizeof(PbObject), // tp_basicsize + 0, // tp_itemsize + (destructor)cbDealloc, // tp_dealloc + 0, // tp_print + 0, // tp_getattr + 0, // tp_setattr + 0, // tp_reserved + 0, // tp_repr + num, // tp_as_number + 0, // tp_as_sequence + 0, // tp_as_mapping + 0, // tp_hash + 0, // tp_call + 0, // tp_str + 0, // tp_getattro + 0, // tp_setattro + 0, // tp_as_buffer + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, // tp_flags + nameptr, // tp_doc + 0, // tp_traverse + 0, // tp_clear + 0, // tp_richcompare + 0, // tp_weaklistoffset + 0, // tp_iter + 0, // tp_iternext + &data.genMethods[0], // tp_methods + 0, // tp_members + &data.genGetSet[0], // tp_getset + 0, // tp_base + 0, // tp_dict + 0, // tp_descr_get + 0, // tp_descr_set + 0, // tp_dictoffset + (initproc)(data.constructor), // tp_init + 0, // tp_alloc + cbNew // tp_new + }; + data.typeInfo = t; + + if (PyType_Ready(&data.typeInfo) < 0) + continue; + + for (map<string, ClassData *>::iterator i2 = mClasses.begin(); i2 != mClasses.end(); ++i2) { + if (*it != i2->second) + continue; + // register all aliases + Py_INCREF(castPy(&data.typeInfo)); + PyModule_AddObject(module, (char *)i2->first.c_str(), (PyObject *)&data.typeInfo); + } + } + + // externals + for (vector<InitFunc>::iterator it = mExtInitializers.begin(); it != mExtInitializers.end(); + ++it) { + (*it)(module); + } + + addConstants(module); + + return module; +} + +//****************************************************** +// Register members and exposed functions + +void setup(const std::string &filename, const std::vector<std::string> &args) +{ + WrapperRegistry::instance().construct(filename, args); + Py_Initialize(); + WrapperRegistry::instance().runPreInit(); +} + +void finalize() +{ + Py_Finalize(); + WrapperRegistry::instance().cleanup(); +} + +bool canConvert(PyObject *obj, const string &classname) +{ + ClassData *from = ((PbObject *)obj)->classdef; + ClassData *dest = WrapperRegistry::instance().lookup(classname); + if (!dest) + errMsg("Classname '" + classname + "' is not registered."); + return WrapperRegistry::instance().canConvert(from, dest); +} + +Manta::PbClass *objFromPy(PyObject *obj) +{ + if (Py_TYPE(obj)->tp_dealloc != (destructor)cbDealloc) // not a manta object + return NULL; + + return ((PbObject *)obj)->instance; +} + +PyObject *copyObject(Manta::PbClass *cls, const string &classname) +{ + ClassData *classdef = WrapperRegistry::instance().lookup(classname); + assertMsg(classdef, "python class " + classname + " does not exist."); + + // allocate new object + PbObject *obj = (PbObject *)classdef->typeInfo.tp_alloc(&(classdef->typeInfo), 0); + assertMsg(obj, "cannot allocate new python object"); + + obj->classdef = classdef; + cls->registerObject((PyObject *)obj, 0); + + return cls->getPyObject(); +} + +Manta::PbClass *createPy(const std::string &classname, + const std::string &name, + Manta::PbArgs &args, + Manta::PbClass *parent) +{ + PyObject *obj = WrapperRegistry::instance().createPyObject(classname, name, args, parent); + return ((PbObject *)obj)->instance; +} + +void setReference(Manta::PbClass *cls, PyObject *obj) +{ + ((PbObject *)obj)->instance = cls; +} + +Register::Register(const string &className, const string &funcName, GenericFunction func) +{ + WrapperRegistry::instance().addMethod(className, funcName, func); +} +Register::Register(const string &className, const string &funcName, OperatorFunction func) +{ + WrapperRegistry::instance().addOperator(className, funcName, func); +} +Register::Register(const string &className, const string &funcName, Constructor func) +{ + WrapperRegistry::instance().addConstructor(className, func); +} +Register::Register(const string &className, const string &property, Getter getter, Setter setter) +{ + WrapperRegistry::instance().addGetSet(className, property, getter, setter); +} +Register::Register(const string &className, const string &pyName, const string &baseClass) +{ + WrapperRegistry::instance().addClass(pyName, className, baseClass); +} +Register::Register(const string &name, const int value) +{ + WrapperRegistry::instance().addEnumEntry(name, value); +} +Register::Register(const string &file, const string &pythonCode) +{ + WrapperRegistry::instance().addPythonCode(file, pythonCode); +} +Register::Register(InitFunc func) +{ + WrapperRegistry::instance().addExternalInitializer(func); +} + +} // namespace Pb diff --git a/extern/mantaflow/helper/pwrapper/registry.h b/extern/mantaflow/helper/pwrapper/registry.h new file mode 100644 index 00000000000..139863df85d --- /dev/null +++ b/extern/mantaflow/helper/pwrapper/registry.h @@ -0,0 +1,106 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011-2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Auto python registry + * + ******************************************************************************/ + +#ifndef _REGISTRY_H +#define _REGISTRY_H + +#include <string> +#include <vector> + +// forward declaration to minimize Python.h includes +#ifndef PyObject_HEAD +# ifndef PyObject_Fake +struct _object; +typedef _object PyObject; +# define PyObject_Fake +# endif +#endif + +namespace Manta { +class PbClass; +class PbArgs; +} // namespace Manta + +// ************************************************** +// NOTE +// Everything in this file is intend only for internal +// use by the generated wrappers or pclass/pconvert. +// For user code, use the functionality exposed in +// pclass.h / pconvert.h instead. +// ************************************************** + +// Used to turn names into strings +namespace Manta { +template<class T> struct Namify { + static const char *S; +}; +} // namespace Manta +namespace Pb { + +// internal registry access +void setup(const std::string &filename, const std::vector<std::string> &args); +void finalize(); +bool canConvert(PyObject *obj, const std::string &to); +Manta::PbClass *objFromPy(PyObject *obj); +Manta::PbClass *createPy(const std::string &classname, + const std::string &name, + Manta::PbArgs &args, + Manta::PbClass *parent); +void setReference(Manta::PbClass *cls, PyObject *obj); +PyObject *copyObject(Manta::PbClass *cls, const std::string &classname); +void MantaEnsureRegistration(); + +#ifdef BLENDER +# ifdef PyMODINIT_FUNC +PyMODINIT_FUNC PyInit_Main(void); +# endif +#endif + +// callback type +typedef void (*InitFunc)(PyObject *); +typedef PyObject *(*GenericFunction)(PyObject *self, PyObject *args, PyObject *kwds); +typedef PyObject *(*OperatorFunction)(PyObject *self, PyObject *o); +typedef int (*Constructor)(PyObject *self, PyObject *args, PyObject *kwds); +typedef PyObject *(*Getter)(PyObject *self, void *closure); +typedef int (*Setter)(PyObject *self, PyObject *value, void *closure); + +//! Auto registry of python methods and classes +struct Register { + //! register method + Register(const std::string &className, const std::string &funcName, GenericFunction func); + //! register operator + Register(const std::string &className, const std::string &funcName, OperatorFunction func); + //! register constructor + Register(const std::string &className, const std::string &funcName, Constructor func); + //! register getter/setter + Register(const std::string &className, + const std::string &property, + Getter getter, + Setter setter); + //! register class + Register(const std::string &className, const std::string &pyName, const std::string &baseClass); + //! register enum entry + Register(const std::string &name, const int value); + //! register python code + Register(const std::string &file, const std::string &pythonCode); + //! register external code + Register(InitFunc func); +}; + +#define KEEP_UNUSED(var) \ + do { \ + (void)var; \ + } while (false); + +} // namespace Pb +#endif diff --git a/extern/mantaflow/helper/util/integrator.h b/extern/mantaflow/helper/util/integrator.h new file mode 100644 index 00000000000..5b1b02a5197 --- /dev/null +++ b/extern/mantaflow/helper/util/integrator.h @@ -0,0 +1,79 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Helper functions for simple integration + * + ******************************************************************************/ + +#ifndef _INTEGRATE_H +#define _INTEGRATE_H + +#include <vector> +#include "vectorbase.h" +#include "kernel.h" + +namespace Manta { + +enum IntegrationMode { IntEuler = 0, IntRK2, IntRK4 }; + +//! Integrate a particle set with a given velocity kernel +template<class VelKernel> void integratePointSet(VelKernel &k, int mode) +{ + typedef typename VelKernel::type0 PosType; + PosType &x = k.getArg0(); + const std::vector<Vec3> &u = k.getRet(); + const int N = x.size(); + + if (mode == IntEuler) { + for (int i = 0; i < N; i++) + x[i].pos += u[i]; + } + else if (mode == IntRK2) { + PosType x0(x); + + for (int i = 0; i < N; i++) + x[i].pos = x0[i].pos + 0.5 * u[i]; + + k.run(); + for (int i = 0; i < N; i++) + x[i].pos = x0[i].pos + u[i]; + } + else if (mode == IntRK4) { + PosType x0(x); + std::vector<Vec3> uTotal(u); + + for (int i = 0; i < N; i++) + x[i].pos = x0[i].pos + 0.5 * u[i]; + + k.run(); + for (int i = 0; i < N; i++) { + x[i].pos = x0[i].pos + 0.5 * u[i]; + uTotal[i] += 2 * u[i]; + } + + k.run(); + for (int i = 0; i < N; i++) { + x[i].pos = x0[i].pos + u[i]; + uTotal[i] += 2 * u[i]; + } + + k.run(); + for (int i = 0; i < N; i++) + x[i].pos = x0[i].pos + (Real)(1. / 6.) * (uTotal[i] + u[i]); + } + else + errMsg("unknown integration type"); + + // for(int i=0; i<N; i++) std::cout << x[i].pos.y-x[0].pos.y << std::endl; + // std::cout << "<><><>" << std::endl; +} + +} // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/interpol.h b/extern/mantaflow/helper/util/interpol.h new file mode 100644 index 00000000000..24d5d2ada06 --- /dev/null +++ b/extern/mantaflow/helper/util/interpol.h @@ -0,0 +1,324 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Helper functions for interpolation + * + ******************************************************************************/ + +#ifndef _INTERPOL_H +#define _INTERPOL_H + +#include "vectorbase.h" + +// Grid values are stored at i+0.5, j+0.5, k+0.5 +// MAC grid values are stored at i,j+0.5,k+0.5 (for x) ... + +namespace Manta { + +inline Vec3 fdTangent(const Vec3 &p0, const Vec3 &p1, const Vec3 &p2) +{ + return 0.5 * (getNormalized(p2 - p1) + getNormalized(p1 - p0)); +} + +inline Vec3 crTangent(const Vec3 &p0, const Vec3 &p1, const Vec3 &p2) +{ + return 0.5 * (p2 - p0); +} + +inline Vec3 hermiteSpline(const Vec3 &p0, const Vec3 &p1, const Vec3 &m0, const Vec3 &m1, Real t) +{ + const Real t2 = t * t, t3 = t2 * t; + return (2.0 * t3 - 3.0 * t2 + 1.0) * p0 + (t3 - 2.0 * t2 + t) * m0 + + (-2.0 * t3 + 3.0 * t2) * p1 + (t3 - t2) * m1; +} + +static inline void checkIndexInterpol(const Vec3i &size, IndexInt idx) +{ + if (idx < 0 || idx > (IndexInt)size.x * size.y * size.z) { + std::ostringstream s; + s << "Grid interpol dim " << size << " : index " << idx << " out of bound "; + errMsg(s.str()); + } +} + +// ---------------------------------------------------------------------- +// Grid interpolators +// ---------------------------------------------------------------------- + +#define BUILD_INDEX \ + Real px = pos.x - 0.5f, py = pos.y - 0.5f, pz = pos.z - 0.5f; \ + int xi = (int)px; \ + int yi = (int)py; \ + int zi = (int)pz; \ + Real s1 = px - (Real)xi, s0 = 1. - s1; \ + Real t1 = py - (Real)yi, t0 = 1. - t1; \ + Real f1 = pz - (Real)zi, f0 = 1. - f1; \ + /* clamp to border */ \ + if (px < 0.) { \ + xi = 0; \ + s0 = 1.0; \ + s1 = 0.0; \ + } \ + if (py < 0.) { \ + yi = 0; \ + t0 = 1.0; \ + t1 = 0.0; \ + } \ + if (pz < 0.) { \ + zi = 0; \ + f0 = 1.0; \ + f1 = 0.0; \ + } \ + if (xi >= size.x - 1) { \ + xi = size.x - 2; \ + s0 = 0.0; \ + s1 = 1.0; \ + } \ + if (yi >= size.y - 1) { \ + yi = size.y - 2; \ + t0 = 0.0; \ + t1 = 1.0; \ + } \ + if (size.z > 1) { \ + if (zi >= size.z - 1) { \ + zi = size.z - 2; \ + f0 = 0.0; \ + f1 = 1.0; \ + } \ + } \ + const int X = 1; \ + const int Y = size.x; + +template<class T> inline T interpol(const T *data, const Vec3i &size, const int Z, const Vec3 &pos) +{ + BUILD_INDEX + IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi; + DEBUG_ONLY(checkIndexInterpol(size, idx)); + DEBUG_ONLY(checkIndexInterpol(size, idx + X + Y + Z)); + + return ((data[idx] * t0 + data[idx + Y] * t1) * s0 + + (data[idx + X] * t0 + data[idx + X + Y] * t1) * s1) * + f0 + + ((data[idx + Z] * t0 + data[idx + Y + Z] * t1) * s0 + + (data[idx + X + Z] * t0 + data[idx + X + Y + Z] * t1) * s1) * + f1; +} + +template<int c> +inline Real interpolComponent(const Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos) +{ + BUILD_INDEX + IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi; + DEBUG_ONLY(checkIndexInterpol(size, idx)); + DEBUG_ONLY(checkIndexInterpol(size, idx + X + Y + Z)); + + return ((data[idx][c] * t0 + data[idx + Y][c] * t1) * s0 + + (data[idx + X][c] * t0 + data[idx + X + Y][c] * t1) * s1) * + f0 + + ((data[idx + Z][c] * t0 + data[idx + Y + Z][c] * t1) * s0 + + (data[idx + X + Z][c] * t0 + data[idx + X + Y + Z][c] * t1) * s1) * + f1; +} + +template<class T> +inline void setInterpol( + T *data, const Vec3i &size, const int Z, const Vec3 &pos, const T &v, Real *sumBuffer) +{ + BUILD_INDEX + IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi; + DEBUG_ONLY(checkIndexInterpol(size, idx)); + DEBUG_ONLY(checkIndexInterpol(size, idx + X + Y + Z)); + + T *ref = &data[idx]; + Real *sum = &sumBuffer[idx]; + Real s0f0 = s0 * f0, s1f0 = s1 * f0, s0f1 = s0 * f1, s1f1 = s1 * f1; + Real w0 = t0 * s0f0, wx = t0 * s1f0, wy = t1 * s0f0, wxy = t1 * s1f0; + Real wz = t0 * s0f1, wxz = t0 * s1f1, wyz = t1 * s0f1, wxyz = t1 * s1f1; + + sum[Z] += wz; + sum[X + Z] += wxz; + sum[Y + Z] += wyz; + sum[X + Y + Z] += wxyz; + ref[Z] += wz * v; + ref[X + Z] += wxz * v; + ref[Y + Z] += wyz * v; + ref[X + Y + Z] += wxyz * v; + sum[0] += w0; + sum[X] += wx; + sum[Y] += wy; + sum[X + Y] += wxy; + ref[0] += w0 * v; + ref[X] += wx * v; + ref[Y] += wy * v; + ref[X + Y] += wxy * v; +} + +#define BUILD_INDEX_SHIFT \ + BUILD_INDEX \ + /* shifted coords */ \ + int s_xi = (int)pos.x, s_yi = (int)pos.y, s_zi = (int)pos.z; \ + Real s_s1 = pos.x - (Real)s_xi, s_s0 = 1. - s_s1; \ + Real s_t1 = pos.y - (Real)s_yi, s_t0 = 1. - s_t1; \ + Real s_f1 = pos.z - (Real)s_zi, s_f0 = 1. - s_f1; \ + /* clamp to border */ \ + if (pos.x < 0) { \ + s_xi = 0; \ + s_s0 = 1.0; \ + s_s1 = 0.0; \ + } \ + if (pos.y < 0) { \ + s_yi = 0; \ + s_t0 = 1.0; \ + s_t1 = 0.0; \ + } \ + if (pos.z < 0) { \ + s_zi = 0; \ + s_f0 = 1.0; \ + s_f1 = 0.0; \ + } \ + if (s_xi >= size.x - 1) { \ + s_xi = size.x - 2; \ + s_s0 = 0.0; \ + s_s1 = 1.0; \ + } \ + if (s_yi >= size.y - 1) { \ + s_yi = size.y - 2; \ + s_t0 = 0.0; \ + s_t1 = 1.0; \ + } \ + if (size.z > 1) { \ + if (s_zi >= size.z - 1) { \ + s_zi = size.z - 2; \ + s_f0 = 0.0; \ + s_f1 = 1.0; \ + } \ + } + +inline Vec3 interpolMAC(const Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos) +{ + BUILD_INDEX_SHIFT; + DEBUG_ONLY(checkIndexInterpol(size, (zi * (IndexInt)size.y + yi) * (IndexInt)size.x + xi)); + DEBUG_ONLY(checkIndexInterpol( + size, (s_zi * (IndexInt)size.y + s_yi) * (IndexInt)size.x + s_xi + X + Y + Z)); + + // process individual components + Vec3 ret(0.); + { // X + const Vec3 *ref = &data[((zi * size.y + yi) * size.x + s_xi)]; + ret.x = f0 * ((ref[0].x * t0 + ref[Y].x * t1) * s_s0 + + (ref[X].x * t0 + ref[X + Y].x * t1) * s_s1) + + f1 * ((ref[Z].x * t0 + ref[Z + Y].x * t1) * s_s0 + + (ref[X + Z].x * t0 + ref[X + Y + Z].x * t1) * s_s1); + } + { // Y + const Vec3 *ref = &data[((zi * size.y + s_yi) * size.x + xi)]; + ret.y = f0 * ((ref[0].y * s_t0 + ref[Y].y * s_t1) * s0 + + (ref[X].y * s_t0 + ref[X + Y].y * s_t1) * s1) + + f1 * ((ref[Z].y * s_t0 + ref[Z + Y].y * s_t1) * s0 + + (ref[X + Z].y * s_t0 + ref[X + Y + Z].y * s_t1) * s1); + } + { // Z + const Vec3 *ref = &data[((s_zi * size.y + yi) * size.x + xi)]; + ret.z = s_f0 * + ((ref[0].z * t0 + ref[Y].z * t1) * s0 + (ref[X].z * t0 + ref[X + Y].z * t1) * s1) + + s_f1 * ((ref[Z].z * t0 + ref[Z + Y].z * t1) * s0 + + (ref[X + Z].z * t0 + ref[X + Y + Z].z * t1) * s1); + } + return ret; +} + +inline void setInterpolMAC( + Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos, const Vec3 &val, Vec3 *sumBuffer) +{ + BUILD_INDEX_SHIFT; + DEBUG_ONLY(checkIndexInterpol(size, (zi * (IndexInt)size.y + yi) * (IndexInt)size.x + xi)); + DEBUG_ONLY(checkIndexInterpol( + size, (s_zi * (IndexInt)size.y + s_yi) * (IndexInt)size.x + s_xi + X + Y + Z)); + + // process individual components + { // X + const IndexInt idx = (IndexInt)(zi * size.y + yi) * size.x + s_xi; + Vec3 *ref = &data[idx], *sum = &sumBuffer[idx]; + Real s0f0 = s_s0 * f0, s1f0 = s_s1 * f0, s0f1 = s_s0 * f1, s1f1 = s_s1 * f1; + Real w0 = t0 * s0f0, wx = t0 * s1f0, wy = t1 * s0f0, wxy = t1 * s1f0; + Real wz = t0 * s0f1, wxz = t0 * s1f1, wyz = t1 * s0f1, wxyz = t1 * s1f1; + + sum[Z].x += wz; + sum[X + Z].x += wxz; + sum[Y + Z].x += wyz; + sum[X + Y + Z].x += wxyz; + ref[Z].x += wz * val.x; + ref[X + Z].x += wxz * val.x; + ref[Y + Z].x += wyz * val.x; + ref[X + Y + Z].x += wxyz * val.x; + sum[0].x += w0; + sum[X].x += wx; + sum[Y].x += wy; + sum[X + Y].x += wxy; + ref[0].x += w0 * val.x; + ref[X].x += wx * val.x; + ref[Y].x += wy * val.x; + ref[X + Y].x += wxy * val.x; + } + { // Y + const IndexInt idx = (IndexInt)(zi * size.y + s_yi) * size.x + xi; + Vec3 *ref = &data[idx], *sum = &sumBuffer[idx]; + Real s0f0 = s0 * f0, s1f0 = s1 * f0, s0f1 = s0 * f1, s1f1 = s1 * f1; + Real w0 = s_t0 * s0f0, wx = s_t0 * s1f0, wy = s_t1 * s0f0, wxy = s_t1 * s1f0; + Real wz = s_t0 * s0f1, wxz = s_t0 * s1f1, wyz = s_t1 * s0f1, wxyz = s_t1 * s1f1; + + sum[Z].y += wz; + sum[X + Z].y += wxz; + sum[Y + Z].y += wyz; + sum[X + Y + Z].y += wxyz; + ref[Z].y += wz * val.y; + ref[X + Z].y += wxz * val.y; + ref[Y + Z].y += wyz * val.y; + ref[X + Y + Z].y += wxyz * val.y; + sum[0].y += w0; + sum[X].y += wx; + sum[Y].y += wy; + sum[X + Y].y += wxy; + ref[0].y += w0 * val.y; + ref[X].y += wx * val.y; + ref[Y].y += wy * val.y; + ref[X + Y].y += wxy * val.y; + } + { // Z + const IndexInt idx = (IndexInt)(s_zi * size.y + yi) * size.x + xi; + Vec3 *ref = &data[idx], *sum = &sumBuffer[idx]; + Real s0f0 = s0 * s_f0, s1f0 = s1 * s_f0, s0f1 = s0 * s_f1, s1f1 = s1 * s_f1; + Real w0 = t0 * s0f0, wx = t0 * s1f0, wy = t1 * s0f0, wxy = t1 * s1f0; + Real wz = t0 * s0f1, wxz = t0 * s1f1, wyz = t1 * s0f1, wxyz = t1 * s1f1; + + sum[0].z += w0; + sum[X].z += wx; + sum[Y].z += wy; + sum[X + Y].z += wxy; + sum[Z].z += wz; + sum[X + Z].z += wxz; + sum[Y + Z].z += wyz; + sum[X + Y + Z].z += wxyz; + ref[0].z += w0 * val.z; + ref[X].z += wx * val.z; + ref[Y].z += wy * val.z; + ref[X + Y].z += wxy * val.z; + ref[Z].z += wz * val.z; + ref[X + Z].z += wxz * val.z; + ref[Y + Z].z += wyz * val.z; + ref[X + Y + Z].z += wxyz * val.z; + } +} + +#undef BUILD_INDEX +#undef BUILD_INDEX_SHIFT + +} // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/interpolHigh.h b/extern/mantaflow/helper/util/interpolHigh.h new file mode 100644 index 00000000000..c2a77442b9c --- /dev/null +++ b/extern/mantaflow/helper/util/interpolHigh.h @@ -0,0 +1,204 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Helper functions for higher order interpolation + * + ******************************************************************************/ + +#ifndef _INTERPOLHIGH_H +#define _INTERPOLHIGH_H + +#include "vectorbase.h" + +namespace Manta { + +template<class T> inline T cubicInterp(const Real interp, const T *points) +{ + T d0 = (points[2] - points[0]) * 0.5; + T d1 = (points[3] - points[1]) * 0.5; + T deltak = (points[2] - points[1]); + + // disabled: if (deltak * d0 < 0.0) d0 = 0; + // disabled: if (deltak * d1 < 0.0) d1 = 0; + + T a0 = points[1]; + T a1 = d0; + T a2 = 3.0 * deltak - 2.0 * d0 - d1; + T a3 = -2.0 * deltak + d0 + d1; + + Real squared = interp * interp; + Real cubed = squared * interp; + return a3 * cubed + a2 * squared + a1 * interp + a0; +} + +template<class T> inline T interpolCubic2D(const T *data, const Vec3i &size, const Vec3 &pos) +{ + const Real px = pos.x - 0.5f, py = pos.y - 0.5f; + + const int x1 = (int)px; + const int x2 = x1 + 1; + const int x3 = x1 + 2; + const int x0 = x1 - 1; + + const int y1 = (int)py; + const int y2 = y1 + 1; + const int y3 = y1 + 2; + const int y0 = y1 - 1; + + if (x0 < 0 || y0 < 0 || x3 >= size[0] || y3 >= size[1]) { + return interpol(data, size, 0, pos); + } + + const Real xInterp = px - x1; + const Real yInterp = py - y1; + + const int y0x = y0 * size[0]; + const int y1x = y1 * size[0]; + const int y2x = y2 * size[0]; + const int y3x = y3 * size[0]; + + const T p0[] = {data[x0 + y0x], data[x1 + y0x], data[x2 + y0x], data[x3 + y0x]}; + const T p1[] = {data[x0 + y1x], data[x1 + y1x], data[x2 + y1x], data[x3 + y1x]}; + const T p2[] = {data[x0 + y2x], data[x1 + y2x], data[x2 + y2x], data[x3 + y2x]}; + const T p3[] = {data[x0 + y3x], data[x1 + y3x], data[x2 + y3x], data[x3 + y3x]}; + + const T finalPoints[] = {cubicInterp(xInterp, p0), + cubicInterp(xInterp, p1), + cubicInterp(xInterp, p2), + cubicInterp(xInterp, p3)}; + + return cubicInterp(yInterp, finalPoints); +} + +template<class T> +inline T interpolCubic(const T *data, const Vec3i &size, const int Z, const Vec3 &pos) +{ + if (Z == 0) + return interpolCubic2D(data, size, pos); + + const Real px = pos.x - 0.5f, py = pos.y - 0.5f, pz = pos.z - 0.5f; + + const int x1 = (int)px; + const int x2 = x1 + 1; + const int x3 = x1 + 2; + const int x0 = x1 - 1; + + const int y1 = (int)py; + const int y2 = y1 + 1; + const int y3 = y1 + 2; + const int y0 = y1 - 1; + + const int z1 = (int)pz; + const int z2 = z1 + 1; + const int z3 = z1 + 2; + const int z0 = z1 - 1; + + if (x0 < 0 || y0 < 0 || z0 < 0 || x3 >= size[0] || y3 >= size[1] || z3 >= size[2]) { + return interpol(data, size, Z, pos); + } + + const Real xInterp = px - x1; + const Real yInterp = py - y1; + const Real zInterp = pz - z1; + + const int slabsize = size[0] * size[1]; + const int z0Slab = z0 * slabsize; + const int z1Slab = z1 * slabsize; + const int z2Slab = z2 * slabsize; + const int z3Slab = z3 * slabsize; + + const int y0x = y0 * size[0]; + const int y1x = y1 * size[0]; + const int y2x = y2 * size[0]; + const int y3x = y3 * size[0]; + + const int y0z0 = y0x + z0Slab; + const int y1z0 = y1x + z0Slab; + const int y2z0 = y2x + z0Slab; + const int y3z0 = y3x + z0Slab; + + const int y0z1 = y0x + z1Slab; + const int y1z1 = y1x + z1Slab; + const int y2z1 = y2x + z1Slab; + const int y3z1 = y3x + z1Slab; + + const int y0z2 = y0x + z2Slab; + const int y1z2 = y1x + z2Slab; + const int y2z2 = y2x + z2Slab; + const int y3z2 = y3x + z2Slab; + + const int y0z3 = y0x + z3Slab; + const int y1z3 = y1x + z3Slab; + const int y2z3 = y2x + z3Slab; + const int y3z3 = y3x + z3Slab; + + // get the z0 slice + const T p0[] = {data[x0 + y0z0], data[x1 + y0z0], data[x2 + y0z0], data[x3 + y0z0]}; + const T p1[] = {data[x0 + y1z0], data[x1 + y1z0], data[x2 + y1z0], data[x3 + y1z0]}; + const T p2[] = {data[x0 + y2z0], data[x1 + y2z0], data[x2 + y2z0], data[x3 + y2z0]}; + const T p3[] = {data[x0 + y3z0], data[x1 + y3z0], data[x2 + y3z0], data[x3 + y3z0]}; + + // get the z1 slice + const T p4[] = {data[x0 + y0z1], data[x1 + y0z1], data[x2 + y0z1], data[x3 + y0z1]}; + const T p5[] = {data[x0 + y1z1], data[x1 + y1z1], data[x2 + y1z1], data[x3 + y1z1]}; + const T p6[] = {data[x0 + y2z1], data[x1 + y2z1], data[x2 + y2z1], data[x3 + y2z1]}; + const T p7[] = {data[x0 + y3z1], data[x1 + y3z1], data[x2 + y3z1], data[x3 + y3z1]}; + + // get the z2 slice + const T p8[] = {data[x0 + y0z2], data[x1 + y0z2], data[x2 + y0z2], data[x3 + y0z2]}; + const T p9[] = {data[x0 + y1z2], data[x1 + y1z2], data[x2 + y1z2], data[x3 + y1z2]}; + const T p10[] = {data[x0 + y2z2], data[x1 + y2z2], data[x2 + y2z2], data[x3 + y2z2]}; + const T p11[] = {data[x0 + y3z2], data[x1 + y3z2], data[x2 + y3z2], data[x3 + y3z2]}; + + // get the z3 slice + const T p12[] = {data[x0 + y0z3], data[x1 + y0z3], data[x2 + y0z3], data[x3 + y0z3]}; + const T p13[] = {data[x0 + y1z3], data[x1 + y1z3], data[x2 + y1z3], data[x3 + y1z3]}; + const T p14[] = {data[x0 + y2z3], data[x1 + y2z3], data[x2 + y2z3], data[x3 + y2z3]}; + const T p15[] = {data[x0 + y3z3], data[x1 + y3z3], data[x2 + y3z3], data[x3 + y3z3]}; + + // interpolate + const T z0Points[] = {cubicInterp(xInterp, p0), + cubicInterp(xInterp, p1), + cubicInterp(xInterp, p2), + cubicInterp(xInterp, p3)}; + const T z1Points[] = {cubicInterp(xInterp, p4), + cubicInterp(xInterp, p5), + cubicInterp(xInterp, p6), + cubicInterp(xInterp, p7)}; + const T z2Points[] = {cubicInterp(xInterp, p8), + cubicInterp(xInterp, p9), + cubicInterp(xInterp, p10), + cubicInterp(xInterp, p11)}; + const T z3Points[] = {cubicInterp(xInterp, p12), + cubicInterp(xInterp, p13), + cubicInterp(xInterp, p14), + cubicInterp(xInterp, p15)}; + + const T finalPoints[] = {cubicInterp(yInterp, z0Points), + cubicInterp(yInterp, z1Points), + cubicInterp(yInterp, z2Points), + cubicInterp(yInterp, z3Points)}; + + return cubicInterp(zInterp, finalPoints); +} + +inline Vec3 interpolCubicMAC(const Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos) +{ + // warning - not yet optimized... + Real vx = interpolCubic<Vec3>(data, size, Z, pos + Vec3(0.5, 0, 0))[0]; + Real vy = interpolCubic<Vec3>(data, size, Z, pos + Vec3(0, 0.5, 0))[1]; + Real vz = 0.f; + if (Z != 0) + vz = interpolCubic<Vec3>(data, size, Z, pos + Vec3(0, 0, 0.5))[2]; + return Vec3(vx, vy, vz); +} + +} // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/matrixbase.h b/extern/mantaflow/helper/util/matrixbase.h new file mode 100644 index 00000000000..9875998f0be --- /dev/null +++ b/extern/mantaflow/helper/util/matrixbase.h @@ -0,0 +1,394 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2015 Kiwon Um, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * GNU General Public License (GPL) + * http://www.gnu.org/licenses + * + * Matrix (3x3) class + * + ******************************************************************************/ + +#ifndef MATRIXBASE_H +#define MATRIXBASE_H + +#include "vectorbase.h" + +namespace Manta { + +template<typename T> class Matrix3x3 { + public: + // NOTE: default is the identity matrix! + explicit Matrix3x3(const T &p00 = 1, + const T &p01 = 0, + const T &p02 = 0, + const T &p10 = 0, + const T &p11 = 1, + const T &p12 = 0, + const T &p20 = 0, + const T &p21 = 0, + const T &p22 = 1) + { + v[0][0] = p00; + v[0][1] = p01; + v[0][2] = p02; + v[1][0] = p10; + v[1][1] = p11; + v[1][2] = p12; + v[2][0] = p20; + v[2][1] = p21; + v[2][2] = p22; + } + + explicit Matrix3x3(const Vector3D<T> &diag) + { + v[0][0] = diag.x; + v[0][1] = 0; + v[0][2] = 0; + v[1][0] = 0; + v[1][1] = diag.y; + v[1][2] = 0; + v[2][0] = 0; + v[2][1] = 0; + v[2][2] = diag.z; + } + + Matrix3x3(const Vector3D<T> &c0, const Vector3D<T> &c1, const Vector3D<T> &c2) + { + v[0][0] = c0.x; + v[0][1] = c1.x; + v[0][2] = c2.x; + v[1][0] = c0.y; + v[1][1] = c1.y; + v[1][2] = c2.y; + v[2][0] = c0.z; + v[2][1] = c1.z; + v[2][2] = c2.z; + } + + // assignment operators + Matrix3x3 &operator+=(const Matrix3x3 &m) + { + v00 += m.v00; + v01 += m.v01; + v02 += m.v02; + v10 += m.v10; + v11 += m.v11; + v12 += m.v12; + v20 += m.v20; + v21 += m.v21; + v22 += m.v22; + return *this; + } + Matrix3x3 &operator-=(const Matrix3x3 &m) + { + v00 -= m.v00; + v01 -= m.v01; + v02 -= m.v02; + v10 -= m.v10; + v11 -= m.v11; + v12 -= m.v12; + v20 -= m.v20; + v21 -= m.v21; + v22 -= m.v22; + return *this; + } + Matrix3x3 &operator*=(const T s) + { + v00 *= s; + v01 *= s; + v02 *= s; + v10 *= s; + v11 *= s; + v12 *= s; + v20 *= s; + v21 *= s; + v22 *= s; + return *this; + } + Matrix3x3 &operator/=(const T s) + { + v00 /= s; + v01 /= s; + v02 /= s; + v10 /= s; + v11 /= s; + v12 /= s; + v20 /= s; + v21 /= s; + v22 /= s; + return *this; + } + + // binary operators + Matrix3x3 operator+(const Matrix3x3 &m) const + { + return Matrix3x3(*this) += m; + } + Matrix3x3 operator-(const Matrix3x3 &m) const + { + return Matrix3x3(*this) -= m; + } + Matrix3x3 operator*(const Matrix3x3 &m) const + { + return Matrix3x3(v00 * m.v00 + v01 * m.v10 + v02 * m.v20, + v00 * m.v01 + v01 * m.v11 + v02 * m.v21, + v00 * m.v02 + v01 * m.v12 + v02 * m.v22, + + v10 * m.v00 + v11 * m.v10 + v12 * m.v20, + v10 * m.v01 + v11 * m.v11 + v12 * m.v21, + v10 * m.v02 + v11 * m.v12 + v12 * m.v22, + + v20 * m.v00 + v21 * m.v10 + v22 * m.v20, + v20 * m.v01 + v21 * m.v11 + v22 * m.v21, + v20 * m.v02 + v21 * m.v12 + v22 * m.v22); + } + Matrix3x3 operator*(const T s) const + { + return Matrix3x3(*this) *= s; + } + Vector3D<T> operator*(const Vector3D<T> &v) const + { + return Vector3D<T>(v00 * v.x + v01 * v.y + v02 * v.z, + v10 * v.x + v11 * v.y + v12 * v.z, + v20 * v.x + v21 * v.y + v22 * v.z); + } + Vector3D<T> transposedMul(const Vector3D<T> &v) const + { + // M^T*v + return Vector3D<T>(v00 * v.x + v10 * v.y + v20 * v.z, + v01 * v.x + v11 * v.y + v21 * v.z, + v02 * v.x + v12 * v.y + v22 * v.z); + } + Matrix3x3 transposedMul(const Matrix3x3 &m) const + { + // M^T*M + return Matrix3x3(v00 * m.v00 + v10 * m.v10 + v20 * m.v20, + v00 * m.v01 + v10 * m.v11 + v20 * m.v21, + v00 * m.v02 + v10 * m.v12 + v20 * m.v22, + + v01 * m.v00 + v11 * m.v10 + v21 * m.v20, + v01 * m.v01 + v11 * m.v11 + v21 * m.v21, + v01 * m.v02 + v11 * m.v12 + v21 * m.v22, + + v02 * m.v00 + v12 * m.v10 + v22 * m.v20, + v02 * m.v01 + v12 * m.v11 + v22 * m.v21, + v02 * m.v02 + v12 * m.v12 + v22 * m.v22); + } + Matrix3x3 mulTranspose(const Matrix3x3 &m) const + { + // M*m^T + return Matrix3x3(v00 * m.v00 + v01 * m.v01 + v02 * m.v02, + v00 * m.v10 + v01 * m.v11 + v02 * m.v12, + v00 * m.v20 + v01 * m.v21 + v02 * m.v22, + + v10 * m.v00 + v11 * m.v01 + v12 * m.v02, + v10 * m.v10 + v11 * m.v11 + v12 * m.v12, + v10 * m.v20 + v11 * m.v21 + v12 * m.v22, + + v20 * m.v00 + v21 * m.v01 + v22 * m.v02, + v20 * m.v10 + v21 * m.v11 + v22 * m.v12, + v20 * m.v20 + v21 * m.v21 + v22 * m.v22); + } + + bool operator==(const Matrix3x3 &m) const + { + return (v00 == m.v00 && v01 == m.v01 && v02 == m.v02 && v10 == m.v10 && v11 == m.v11 && + v12 == m.v12 && v20 == m.v20 && v21 == m.v21 && v22 == m.v22); + } + + const T &operator()(const int r, const int c) const + { + return v[r][c]; + } + T &operator()(const int r, const int c) + { + return const_cast<T &>(const_cast<const Matrix3x3 &>(*this)(r, c)); + } + + T trace() const + { + return v00 + v11 + v22; + } + T sumSqr() const + { + return (v00 * v00 + v01 * v01 + v02 * v02 + v10 * v10 + v11 * v11 + v12 * v12 + v20 * v20 + + v21 * v21 + v22 * v22); + } + + Real determinant() const + { + return (v00 * v11 * v22 - v00 * v12 * v21 + v01 * v12 * v20 - v01 * v10 * v22 + + v02 * v10 * v21 - v02 * v11 * v20); + } + Matrix3x3 &transpose() + { + return *this = transposed(); + } + Matrix3x3 transposed() const + { + return Matrix3x3(v00, v10, v20, v01, v11, v21, v02, v12, v22); + } + Matrix3x3 &invert() + { + return *this = inverse(); + } + Matrix3x3 inverse() const + { + const Real det = determinant(); // FIXME: assert(det); + const Real idet = 1e0 / det; + return Matrix3x3(idet * (v11 * v22 - v12 * v21), + idet * (v02 * v21 - v01 * v22), + idet * (v01 * v12 - v02 * v11), + idet * (v12 * v20 - v10 * v22), + idet * (v00 * v22 - v02 * v20), + idet * (v02 * v10 - v00 * v12), + idet * (v10 * v21 - v11 * v20), + idet * (v01 * v20 - v00 * v21), + idet * (v00 * v11 - v01 * v10)); + } + bool getInverse(Matrix3x3 &inv) const + { + const Real det = determinant(); + if (det == 0e0) + return false; // FIXME: is it likely to happen the floating error? + + const Real idet = 1e0 / det; + inv.v00 = idet * (v11 * v22 - v12 * v21); + inv.v01 = idet * (v02 * v21 - v01 * v22); + inv.v02 = idet * (v01 * v12 - v02 * v11); + + inv.v10 = idet * (v12 * v20 - v10 * v22); + inv.v11 = idet * (v00 * v22 - v02 * v20); + inv.v12 = idet * (v02 * v10 - v00 * v12); + + inv.v20 = idet * (v10 * v21 - v11 * v20); + inv.v21 = idet * (v01 * v20 - v00 * v21); + inv.v22 = idet * (v00 * v11 - v01 * v10); + + return true; + } + + Real normOne() const + { + // the maximum absolute column sum of the matrix + return max(std::fabs(v00) + std::fabs(v10) + std::fabs(v20), + std::fabs(v01) + std::fabs(v11) + std::fabs(v21), + std::fabs(v02) + std::fabs(v12) + std::fabs(v22)); + } + Real normInf() const + { + // the maximum absolute row sum of the matrix + return max(std::fabs(v00) + std::fabs(v01) + std::fabs(v02), + std::fabs(v10) + std::fabs(v11) + std::fabs(v12), + std::fabs(v20) + std::fabs(v21) + std::fabs(v22)); + } + + Vector3D<T> eigenvalues() const + { + Vector3D<T> eigen; + + const Real b = -v00 - v11 - v22; + const Real c = v00 * (v11 + v22) + v11 * v22 - v12 * v21 - v01 * v10 - v02 * v20; + Real d = -v00 * (v11 * v22 - v12 * v21) - v20 * (v01 * v12 - v11 * v02) - + v10 * (v02 * v21 - v22 * v01); + const Real f = (3.0 * c - b * b) / 3.0; + const Real g = (2.0 * b * b * b - 9.0 * b * c + 27.0 * d) / 27.0; + const Real h = g * g / 4.0 + f * f * f / 27.0; + + Real sign; + if (h > 0) { + Real r = -g / 2.0 + std::sqrt(h); + if (r < 0) { + r = -r; + sign = -1.0; + } + else + sign = 1.0; + Real s = sign * std::pow(r, 1.0 / 3.0); + Real t = -g / 2.0 - std::sqrt(h); + if (t < 0) { + t = -t; + sign = -1.0; + } + else + sign = 1.0; + Real u = sign * std::pow(t, 1.0 / 3.0); + eigen[0] = (s + u) - b / 3.0; + eigen[1] = eigen[2] = 0; + } + else if (h == 0) { + if (d < 0) { + d = -d; + sign = -1.0; + } + sign = 1.0; + eigen[0] = -1.0 * sign * std::pow(d, 1.0 / 3.0); + eigen[1] = eigen[2] = 0; + } + else { + const Real i = std::sqrt(g * g / 4.0 - h); + const Real j = std::pow(i, 1.0 / 3.0); + const Real k = std::acos(-g / (2.0 * i)); + const Real l = -j; + const Real m = std::cos(k / 3.0); + const Real n = std::sqrt(3.0) * std::sin(k / 3.0); + const Real p = -b / 3.0; + eigen[0] = 2e0 * j * m + p; + eigen[1] = l * (m + n) + p; + eigen[2] = l * (m - n) + p; + } + + return eigen; + } + + static Matrix3x3 I() + { + return Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1); + } + +#ifdef _WIN32 +# pragma warning(disable : 4201) +#endif + union { + struct { + T v00, v01, v02, v10, v11, v12, v20, v21, v22; + }; + T v[3][3]; + T v1[9]; + }; +#ifdef _WIN32 +# pragma warning(default : 4201) +#endif +}; + +template<typename T1, typename T> inline Matrix3x3<T> operator*(const T1 s, const Matrix3x3<T> &m) +{ + return m * static_cast<T>(s); +} + +template<typename T> inline Matrix3x3<T> crossProductMatrix(const Vector3D<T> &v) +{ + return Matrix3x3<T>(0, -v.z, v.y, v.z, 0, -v.x, -v.y, v.x, 0); +} + +template<typename T> inline Matrix3x3<T> outerProduct(const Vector3D<T> &a, const Vector3D<T> &b) +{ + return Matrix3x3<T>(a.x * b.x, + a.x * b.y, + a.x * b.z, + a.y * b.x, + a.y * b.y, + a.y * b.z, + a.z * b.x, + a.z * b.y, + a.z * b.z); +} + +typedef Matrix3x3<Real> Matrix3x3f; + +} // namespace Manta + +#endif /* MATRIXBASE_H */ diff --git a/extern/mantaflow/helper/util/mcubes.h b/extern/mantaflow/helper/util/mcubes.h new file mode 100644 index 00000000000..bd1c780e932 --- /dev/null +++ b/extern/mantaflow/helper/util/mcubes.h @@ -0,0 +1,308 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Marching cubes lookup indices + * + ******************************************************************************/ + +#ifndef _MCUBES_H_ +#define _MCUBES_H_ + +static const int mcEdges[24] = {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 int cubieOffsetX[8] = {0, 1, 1, 0, 0, 1, 1, 0}; +static const int cubieOffsetY[8] = {0, 0, 1, 1, 0, 0, 1, 1}; +static const int cubieOffsetZ[8] = {0, 0, 0, 0, 1, 1, 1, 1}; + +/* which edges are needed ? */ +/* cf. http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/ */ +static const short mcEdgeTable[256] = { + 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, + 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, + 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, + 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0xaa, + 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, + 0x569, 0x663, 0x76a, 0x66, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, + 0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, + 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c, + 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, + 0x2cf, 0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9, + 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, + 0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, + 0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, + 0x3f5, 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, + 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3, + 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, + 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, + 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, + 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, + 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}; + +/* triangles for the 256 intersection possibilities */ +/* cf. http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/ */ +static const short mcTriTable[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}}; + +#endif
\ No newline at end of file diff --git a/extern/mantaflow/helper/util/quaternion.h b/extern/mantaflow/helper/util/quaternion.h new file mode 100644 index 00000000000..c4e161baee2 --- /dev/null +++ b/extern/mantaflow/helper/util/quaternion.h @@ -0,0 +1,103 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Basic quaternion class + * + ******************************************************************************/ + +#ifndef _QUATERNION_H +#define _QUATERNION_H + +#include "vectorbase.h" + +namespace Manta { + +//! Very basic quaternion class +class Quaternion { + public: + //! default constructor + Quaternion() : x(0), y(0), z(0), w(0) + { + } + + //! copy constructor + Quaternion(const Quaternion &q) : x(q.x), y(q.y), z(q.z), w(q.w) + { + } + + //! construct a quaternion from members + Quaternion(Real _x, Real _y, Real _z, Real _w) : x(_x), y(_y), z(_z), w(_w) + { + } + + //! construct a quaternion from imag/real parts + Quaternion(Vec3 i, Real r) : x(i.x), y(i.y), z(i.z), w(r) + { + } + + //! Assign operator + inline Quaternion &operator=(const Quaternion &q) + { + x = q.x; + y = q.y; + z = q.z; + w = q.w; + return *this; + } + + //! Assign multiplication operator + inline Quaternion &operator*=(const Real a) + { + x *= a; + y *= a; + z *= a; + w *= a; + return *this; + } + + //! return inverse quaternion + inline Quaternion inverse() const + { + Real mag = 1.0 / (x * x + y * y + z * z + w * w); + return Quaternion(-x * mag, -y * mag, -z * mag, w * mag); + } + + //! imaginary part accessor + inline Vec3 imag() + { + return Vec3(x, y, z); + } + + // imaginary part + Real x; + Real y; + Real z; + + // real part + Real w; +}; + +//! Multiplication operator +inline Quaternion operator*(const Quaternion &q1, const Quaternion &q2) +{ + return Quaternion(q2.w * q1.x + q2.x * q1.w + q2.y * q1.z - q2.z * q1.y, + q2.w * q1.y + q2.y * q1.w + q2.z * q1.x - q2.x * q1.z, + q2.w * q1.z + q2.z * q1.w + q2.x * q1.y - q2.y * q1.x, + q2.w * q1.w - q2.x * q1.x - q2.y * q1.y - q2.z * q1.z); +} + +//! Multiplication operator +inline Quaternion operator*(const Quaternion &q, const Real a) +{ + return Quaternion(q.x * a, q.y * a, q.z * a, q.w * a); +} + +} // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/randomstream.h b/extern/mantaflow/helper/util/randomstream.h new file mode 100644 index 00000000000..35b9c7d8858 --- /dev/null +++ b/extern/mantaflow/helper/util/randomstream.h @@ -0,0 +1,429 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Random numbers + * + * Based on an example by Makoto Matsumoto, Takuji Nishimura, Shawn Cokus, and Richard J. Wagner + * + ******************************************************************************/ + +#ifndef _RANDOMSTREAM_H +#define _RANDOMSTREAM_H + +namespace Manta { + +#include <iostream> +#include <stdio.h> +#include <time.h> +#include "vectorbase.h" + +class MTRand { + // Data + public: + typedef unsigned long uint32; // unsigned integer type, at least 32 bits + + enum { N = 624 }; // length of state vector + enum { SAVE = N + 1 }; // length of array for save() + + protected: + enum { M = 397 }; // period parameter + + uint32 state[N]; // internal state + uint32 *pNext; // next value to get from state + int left; // number of values left before reload needed + + // Methods + public: + MTRand(const uint32 &oneSeed); // initialize with a simple uint32 + MTRand(uint32 *const bigSeed, uint32 const seedLength = N); // or an array + MTRand(); // auto-initialize with /dev/urandom or time() and clock() + + // Do NOT use for CRYPTOGRAPHY without securely hashing several returned + // values together, otherwise the generator state can be learned after + // reading 624 consecutive values. + + // Access to 32-bit random numbers + double rand(); // real number in [0,1] + double rand(const double &n); // real number in [0,n] + double randExc(); // real number in [0,1) + double randExc(const double &n); // real number in [0,n) + double randDblExc(); // real number in (0,1) + double randDblExc(const double &n); // real number in (0,n) + uint32 randInt(); // integer in [0,2^32-1] + uint32 randInt(const uint32 &n); // integer in [0,n] for n < 2^32 + double operator()() + { + return rand(); + } // same as rand() + + // Access to 53-bit random numbers (capacity of IEEE double precision) + double rand53(); // real number in [0,1) + + // Access to nonuniform random number distributions + double randNorm(const double &mean = 0.0, const double &variance = 1.0); + + // Re-seeding functions with same behavior as initializers + void seed(const uint32 oneSeed); + void seed(uint32 *const bigSeed, const uint32 seedLength = N); + void seed(); + + // Saving and loading generator state + void save(uint32 *saveArray) const; // to array of size SAVE + void load(uint32 *const loadArray); // from such array + friend std::ostream &operator<<(std::ostream &os, const MTRand &mtrand); + friend std::istream &operator>>(std::istream &is, MTRand &mtrand); + + protected: + void initialize(const uint32 oneSeed); + void reload(); + uint32 hiBit(const uint32 &u) const + { + return u & 0x80000000UL; + } + uint32 loBit(const uint32 &u) const + { + return u & 0x00000001UL; + } + uint32 loBits(const uint32 &u) const + { + return u & 0x7fffffffUL; + } + uint32 mixBits(const uint32 &u, const uint32 &v) const + { + return hiBit(u) | loBits(v); + } + uint32 twist(const uint32 &m, const uint32 &s0, const uint32 &s1) const + { + return m ^ (mixBits(s0, s1) >> 1) ^ (-loBit(s1) & 0x9908b0dfUL); + } + static uint32 hash(time_t t, clock_t c); +}; + +inline MTRand::MTRand(const uint32 &oneSeed) +{ + seed(oneSeed); +} + +inline MTRand::MTRand(uint32 *const bigSeed, const uint32 seedLength) +{ + seed(bigSeed, seedLength); +} + +inline MTRand::MTRand() +{ + seed(); +} + +inline double MTRand::rand() +{ + return double(randInt()) * (1.0 / 4294967295.0); +} + +inline double MTRand::rand(const double &n) +{ + return rand() * n; +} + +inline double MTRand::randExc() +{ + return double(randInt()) * (1.0 / 4294967296.0); +} + +inline double MTRand::randExc(const double &n) +{ + return randExc() * n; +} + +inline double MTRand::randDblExc() +{ + return (double(randInt()) + 0.5) * (1.0 / 4294967296.0); +} + +inline double MTRand::randDblExc(const double &n) +{ + return randDblExc() * n; +} + +inline double MTRand::rand53() +{ + uint32 a = randInt() >> 5, b = randInt() >> 6; + return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); // by Isaku Wada +} + +inline double MTRand::randNorm(const double &mean, const double &variance) +{ + // Return a real number from a normal (Gaussian) distribution with given + // mean and variance by Box-Muller method + double r = sqrt(-2.0 * log(1.0 - randDblExc())) * variance; + double phi = 2.0 * 3.14159265358979323846264338328 * randExc(); + return mean + r * cos(phi); +} + +inline MTRand::uint32 MTRand::randInt() +{ + // Pull a 32-bit integer from the generator state + // Every other access function simply transforms the numbers extracted here + + if (left == 0) + reload(); + --left; + + uint32 s1; + s1 = *pNext++; + s1 ^= (s1 >> 11); + s1 ^= (s1 << 7) & 0x9d2c5680UL; + s1 ^= (s1 << 15) & 0xefc60000UL; + return (s1 ^ (s1 >> 18)); +} + +inline MTRand::uint32 MTRand::randInt(const uint32 &n) +{ + // Find which bits are used in n + // Optimized by Magnus Jonsson (magnus@smartelectronix.com) + uint32 used = n; + used |= used >> 1; + used |= used >> 2; + used |= used >> 4; + used |= used >> 8; + used |= used >> 16; + + // Draw numbers until one is found in [0,n] + uint32 i; + do + i = randInt() & used; // toss unused bits to shorten search + while (i > n); + return i; +} + +inline void MTRand::seed(const uint32 oneSeed) +{ + // Seed the generator with a simple uint32 + initialize(oneSeed); + reload(); +} + +inline void MTRand::seed(uint32 *const bigSeed, const uint32 seedLength) +{ + // Seed the generator with an array of uint32's + // There are 2^19937-1 possible initial states. This function allows + // all of those to be accessed by providing at least 19937 bits (with a + // default seed length of N = 624 uint32's). Any bits above the lower 32 + // in each element are discarded. + // Just call seed() if you want to get array from /dev/urandom + initialize(19650218UL); + const unsigned int Nenum = N; + int i = 1; + uint32 j = 0; + int k = (Nenum > seedLength ? Nenum : seedLength); + for (; k; --k) { + state[i] = state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL); + state[i] += (bigSeed[j] & 0xffffffffUL) + j; + state[i] &= 0xffffffffUL; + ++i; + ++j; + if (i >= N) { + state[0] = state[N - 1]; + i = 1; + } + if (j >= seedLength) + j = 0; + } + for (k = N - 1; k; --k) { + state[i] = state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL); + state[i] -= i; + state[i] &= 0xffffffffUL; + ++i; + if (i >= N) { + state[0] = state[N - 1]; + i = 1; + } + } + state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array + reload(); +} + +inline void MTRand::seed() +{ + // Seed the generator with an array from /dev/urandom if available + // Otherwise use a hash of time() and clock() values + + // First try getting an array from /dev/urandom + FILE *urandom = fopen("/dev/urandom", "rb"); + if (urandom) { + uint32 bigSeed[N]; + uint32 *s = bigSeed; + int i = N; + bool success = true; + while (success && i--) + success = fread(s++, sizeof(uint32), 1, urandom); + fclose(urandom); + if (success) { + seed(bigSeed, N); + return; + } + } + + // Was not successful, so use time() and clock() instead + seed(hash(time(NULL), clock())); +} + +inline void MTRand::initialize(const uint32 intseed) +{ + // Initialize generator state with seed + // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier. + // In previous versions, most significant bits (MSBs) of the seed affect + // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto. + uint32 *s = state; + uint32 *r = state; + int i = 1; + *s++ = intseed & 0xffffffffUL; + for (; i < N; ++i) { + *s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL; + r++; + } +} + +inline void MTRand::reload() +{ + // Generate N new values in state + // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com) + uint32 *p = state; + int i; + for (i = N - M; i--; ++p) + *p = twist(p[M], p[0], p[1]); + for (i = M; --i; ++p) + *p = twist(p[M - N], p[0], p[1]); + *p = twist(p[M - N], p[0], state[0]); + + left = N, pNext = state; +} + +inline MTRand::uint32 MTRand::hash(time_t t, clock_t c) +{ + // Get a uint32 from t and c + // Better than uint32(x) in case x is floating point in [0,1] + // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk) + + static uint32 differ = 0; // guarantee time-based seeds will change + + uint32 h1 = 0; + unsigned char *p = (unsigned char *)&t; + for (size_t i = 0; i < sizeof(t); ++i) { + h1 *= std::numeric_limits<unsigned char>::max() + 2U; + h1 += p[i]; + } + uint32 h2 = 0; + p = (unsigned char *)&c; + for (size_t j = 0; j < sizeof(c); ++j) { + h2 *= std::numeric_limits<unsigned char>::max() + 2U; + h2 += p[j]; + } + return (h1 + differ++) ^ h2; +} + +inline void MTRand::save(uint32 *saveArray) const +{ + uint32 *sa = saveArray; + const uint32 *s = state; + int i = N; + for (; i--; *sa++ = *s++) { + } + *sa = left; +} + +inline void MTRand::load(uint32 *const loadArray) +{ + uint32 *s = state; + uint32 *la = loadArray; + int i = N; + for (; i--; *s++ = *la++) { + } + left = *la; + pNext = &state[N - left]; +} + +inline std::ostream &operator<<(std::ostream &os, const MTRand &mtrand) +{ + const MTRand::uint32 *s = mtrand.state; + int i = mtrand.N; + for (; i--; os << *s++ << "\t") { + } + return os << mtrand.left; +} + +inline std::istream &operator>>(std::istream &is, MTRand &mtrand) +{ + MTRand::uint32 *s = mtrand.state; + int i = mtrand.N; + for (; i--; is >> *s++) { + } + is >> mtrand.left; + mtrand.pNext = &mtrand.state[mtrand.N - mtrand.left]; + return is; +} + +// simple interface to mersenne twister +class RandomStream { + public: + inline RandomStream(long seed) : mtr(seed){}; + ~RandomStream() + { + } + + /*! get a random number from the stream */ + inline double getDouble(void) + { + return mtr.rand(); + }; + inline float getFloat(void) + { + return (float)mtr.rand(); + }; + + inline float getFloat(float min, float max) + { + return mtr.rand(max - min) + min; + }; + inline float getRandNorm(float mean, float var) + { + return mtr.randNorm(mean, var); + }; + +#if FLOATINGPOINT_PRECISION == 1 + inline Real getReal() + { + return getFloat(); + } + +#else + inline Real getReal() + { + return getDouble(); + } +#endif + + inline Vec3 getVec3() + { + Real a = getReal(), b = getReal(), c = getReal(); + return Vec3(a, b, c); + } + inline Vec3 getVec3Norm() + { + Vec3 a = getVec3(); + normalize(a); + return a; + } + + private: + MTRand mtr; +}; + +} // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/rcmatrix.h b/extern/mantaflow/helper/util/rcmatrix.h new file mode 100644 index 00000000000..39951cece2e --- /dev/null +++ b/extern/mantaflow/helper/util/rcmatrix.h @@ -0,0 +1,1112 @@ +// +// Helper matrix class, RCMatrix.h +// Required for PD optimizations (guiding) +// Thanks to Ryoichi Ando, and Robert Bridson +// + +#ifndef RCMATRIX3_H +#define RCMATRIX3_H + +#include <iterator> +#include <cassert> +#include <vector> +#include <fstream> + +// index type +#define int_index long long + +// link to omp & tbb for now +#if OPENMP == 1 || TBB == 1 +# define MANTA_ENABLE_PARALLEL 0 +// allow the preconditioner to be computed in parallel? (can lead to slightly non-deterministic +// results) +# define MANTA_ENABLE_PARALLEL_PC 0 +// use c++11 code? +# define MANTA_USE_CPP11 1 +#else +# define MANTA_ENABLE_PARALLEL 0 +# define MANTA_ENABLE_PARALLEL_PC 0 +# define MANTA_USE_CPP11 0 +#endif + +#if MANTA_ENABLE_PARALLEL == 1 +# include <thread> +# include <algorithm> + +static const int manta_num_threads = std::thread::hardware_concurrency(); + +// For clang +# define parallel_for(total_size) \ + { \ + int_index parallel_array_size = (total_size); \ + std::vector<std::thread> threads(manta_num_threads); \ + for (int thread_number = 0; thread_number < manta_num_threads; thread_number++) { \ + threads[thread_number] = std::thread([&](int_index parallel_array_size, int_index thread_number ) { \ + for( int_index parallel_index=thread_number; parallel_index < parallel_array_size; parallel_index += manta_num_threads ) { + +# define parallel_end \ + } \ + },parallel_array_size,thread_number); \ + } \ + for (auto &thread : threads) \ + thread.join(); \ + } + +# define parallel_block \ + { \ + std::vector<std::thread> threads; \ + { + +# define do_parallel threads.push_back( std::thread([&]() { +# define do_end \ + } ) ); + +# define block_end \ + } \ + for (auto &thread : threads) { \ + thread.join(); \ + } \ + } + +#else + +# define parallel_for(size) \ + { \ + int thread_number = 0; \ + int_index parallel_index = 0; \ + for (int_index parallel_index = 0; parallel_index < (int_index)size; parallel_index++) { +# define parallel_end \ + } \ + thread_number = parallel_index = 0; \ + } + +# define parallel_block +# define do_parallel +# define do_end +# define block_end + +#endif + +#include "vectorbase.h" + +namespace Manta { + +static const unsigned default_expected_none_zeros = 7; + +template<class N, class T> struct RCMatrix { + struct RowEntry { + std::vector<N> index; + std::vector<T> value; + }; + RCMatrix() : n(0), expected_none_zeros(default_expected_none_zeros) + { + } + RCMatrix(N size, N expected_none_zeros = default_expected_none_zeros) + : n(0), expected_none_zeros(expected_none_zeros) + { + resize(size); + } + RCMatrix(const RCMatrix &m) : n(0), expected_none_zeros(default_expected_none_zeros) + { + init(m); + } + RCMatrix &operator=(const RCMatrix &m) + { + expected_none_zeros = m.expected_none_zeros; + init(m); + return *this; + } + RCMatrix &operator=(RCMatrix &&m) + { + matrix = m.matrix; + offsets = m.offsets; + expected_none_zeros = m.expected_none_zeros; + n = m.n; + m.n = 0; + m.matrix.clear(); + m.offsets.clear(); + return *this; + } + RCMatrix(RCMatrix &&m) + : n(m.n), expected_none_zeros(m.expected_none_zeros), matrix(m.matrix), offsets(m.offsets) + { + m.n = 0; + m.matrix.clear(); + m.offsets.clear(); + } + void init(const RCMatrix &m) + { + expected_none_zeros = m.expected_none_zeros; + resize(m.n); + parallel_for(n) + { + N i = parallel_index; + if (m.matrix[i]) { + alloc_row(i); + matrix[i]->index = m.matrix[i]->index; + matrix[i]->value = m.matrix[i]->value; + } + else { + dealloc_row(i); + } + } + parallel_end + } + ~RCMatrix() + { + clear(); + } + void clear() + { + for (N i = 0; i < n; i++) { + dealloc_row(i); + matrix[i] = NULL; + if (offsets.size()) + offsets[i] = 0; + } + }; + bool empty(N i) const + { + return matrix[i] == NULL; + } + N row_nonzero_size(N i) const + { + return matrix[i] == NULL ? 0 : matrix[i]->index.size(); + } + void resize(N size, N expected_none_zeros = 0) + { + if (!expected_none_zeros) { + expected_none_zeros = this->expected_none_zeros; + } + if (n > size) { + // Shrinking + for (N i = size ? size - 1 : 0; i < n; i++) + dealloc_row(i); + matrix.resize(size); + } + else if (n < size) { + // Expanding + matrix.resize(size); + for (N i = n; i < size; i++) { + matrix[i] = NULL; + if (offsets.size()) + offsets[i] = 0; + } + } + n = size; + } + void alloc_row(N i) + { + assert(i < n); + if (!matrix[i]) { + matrix[i] = new RowEntry; + matrix[i]->index.reserve(expected_none_zeros); + matrix[i]->value.reserve(expected_none_zeros); + if (offsets.size()) + offsets[i] = 0; + } + } + void dealloc_row(N i) + { + assert(i < n); + if (matrix[i]) { + if (offsets.empty() || !offsets[i]) + delete matrix[i]; + matrix[i] = NULL; + if (offsets.size()) + offsets[i] = 0; + } + } + T operator()(N i, N j) const + { + assert(i < n); + for (Iterator it = row_begin(i); it; ++it) { + if (it.index() == j) + return it.value(); + } + return T(0.0); + } + void add_to_element_checked(N i, N j, T val) + { + if ((i < 0) || (j < 0) || (i >= n) || (j >= n)) + return; + add_to_element(i, j, val); + } + void add_to_element(N i, N j, T increment_value) + { + if (std::abs(increment_value) > VECTOR_EPSILON) { + assert(i < n); + assert(offsets.empty() || offsets[i] == 0); + alloc_row(i); + std::vector<N> &index = matrix[i]->index; + std::vector<T> &value = matrix[i]->value; + for (N k = 0; k < (N)index.size(); ++k) { + if (index[k] == j) { + value[k] += increment_value; + return; + } + else if (index[k] > j) { + index.insert(index.begin() + k, j); + value.insert(value.begin() + k, increment_value); + return; + } + } + index.push_back(j); + value.push_back(increment_value); + } + } + + void set_element(N i, N j, T v) + { + if (std::abs(v) > VECTOR_EPSILON) { + assert(i < n); + assert(offsets.empty() || offsets[i] == 0); + alloc_row(i); + std::vector<N> &index = matrix[i]->index; + std::vector<T> &value = matrix[i]->value; + for (N k = 0; k < (N)index.size(); ++k) { + if (index[k] == j) { + value[k] = v; + return; + } + else if (index[k] > j) { + index.insert(index.begin() + k, j); + value.insert(value.begin() + k, v); + return; + } + } + index.push_back(j); + value.push_back(v); + } + } + + // Make sure that j is the biggest column in the row, no duplication allowed + void fix_element(N i, N j, T v) + { + if (std::abs(v) > VECTOR_EPSILON) { + assert(i < n); + assert(offsets.empty() || offsets[i] == 0); + alloc_row(i); + std::vector<N> &index = matrix[i]->index; + std::vector<T> &value = matrix[i]->value; + index.push_back(j); + value.push_back(v); + } + } + int_index trim_zero_entries(double e = VECTOR_EPSILON) + { + std::vector<int_index> deleted_entries(n, 0); + parallel_for(n) + { + N i = parallel_index; + if (matrix[i]) { + std::vector<N> &index = matrix[i]->index; + std::vector<T> &value = matrix[i]->value; + N head = 0; + N k = 0; + for (k = 0; k < index.size(); ++k) { + if (std::abs(value[k]) > e) { + index[head] = index[k]; + value[head] = value[k]; + ++head; + } + } + if (head != k) { + index.erase(index.begin() + head, index.end()); + value.erase(value.begin() + head, value.end()); + deleted_entries[i] += k - head; + } + if (!offsets.size() && !head) { + remove_row(i); + } + } + } + parallel_end + // + int_index sum_deleted(0); + for (int_index i = 0; i < n; i++) + sum_deleted += deleted_entries[i]; + return sum_deleted; + } + void remove_reference(N i) + { + if (offsets.size() && offsets[i] && matrix[i]) { + RowEntry *save = matrix[i]; + matrix[i] = new RowEntry; + *matrix[i] = *save; + for (N &index : matrix[i]->index) + index += offsets[i]; + offsets[i] = 0; + } + } + void remove_row(N i) + { + dealloc_row(i); + } + bool is_symmetric(double e = VECTOR_EPSILON) const + { + std::vector<bool> flags(n, true); + parallel_for(n) + { + N i = parallel_index; + bool flag = true; + for (Iterator it = row_begin(i); it; ++it) { + N index = it.index(); + T value = it.value(); + if (std::abs(value) > e) { + bool found_entry = false; + for (Iterator it_i = row_begin(index); it_i; ++it_i) { + if (it_i.index() == i) { + found_entry = true; + if (std::abs(value - it_i.value()) > e) { + flag = false; + break; + } + } + } + if (!found_entry) + flag = false; + if (!flag) + break; + } + } + flags[i] = flag; + } + parallel_end for (N i = 0; i < matrix.size(); ++i) + { + if (!flags[i]) + return false; + } + return true; + } + + void expand() + { + if (offsets.empty()) + return; + for (N i = 1; i < n; i++) { + if (offsets[i]) { + RowEntry *ref = matrix[i]; + matrix[i] = new RowEntry; + *matrix[i] = *ref; + for (N j = 0; j < (N)matrix[i]->index.size(); j++) { + matrix[i]->index[j] += offsets[i]; + } + } + } + offsets.resize(0); + } + + N column(N i) const + { + return empty(i) ? 0 : row_begin(i, row_nonzero_size(i) - 1).index(); + } + N getColumnSize() const + { + N max_column(0); + auto column = [&](N i) { + N max_column(0); + for (Iterator it = row_begin(i); it; ++it) + max_column = std::max(max_column, it.index()); + return max_column + 1; + }; + for (N i = 0; i < n; i++) + max_column = std::max(max_column, column(i)); + return max_column; + } + N getNonzeroSize() const + { + N nonzeros(0); + for (N i = 0; i < n; ++i) { + nonzeros += row_nonzero_size(i); + } + return nonzeros; + } + class Iterator : std::iterator<std::input_iterator_tag, T> { + public: + Iterator(const RowEntry *rowEntry, N k, N offset) : rowEntry(rowEntry), k(k), offset(offset) + { + } + operator bool() const + { + return rowEntry != NULL && k < (N)rowEntry->index.size(); + } + Iterator &operator++() + { + ++k; + return *this; + } + T value() const + { + return rowEntry->value[k]; + } + N index() const + { + return rowEntry->index[k] + offset; + } + N index_raw() const + { + return rowEntry->index[k]; + } + N size() const + { + return rowEntry == NULL ? 0 : rowEntry->index.size(); + } + + protected: + const RowEntry *rowEntry; + N k, offset; + }; + Iterator row_begin(N n, N k = 0) const + { + return Iterator(matrix[n], k, offsets.size() ? offsets[n] : 0); + } + class DynamicIterator : public Iterator { + public: + DynamicIterator(RowEntry *rowEntry, N k, N offset) + : rowEntry(rowEntry), Iterator(rowEntry, k, offset) + { + } + void setValue(T value) + { + rowEntry->value[Iterator::k] = value; + } + void setIndex(N index) + { + rowEntry->index[Iterator::k] = index; + } + + protected: + RowEntry *rowEntry; + }; + DynamicIterator dynamic_row_begin(N n, N k = 0) + { + N offset = offsets.size() ? offsets[n] : 0; + if (offset) { + printf("---- Warning ----\n"); + printf("Dynamic iterator is not allowed for referenced rows.\n"); + printf("You should be very careful otherwise this causes some bugs.\n"); + printf( + "We encourage you that you convert this row into a raw format, then loop over it...\n"); + printf("-----------------\n"); + exit(0); + } + return DynamicIterator(matrix[n], k, offset); + } + RCMatrix transpose(N rowsize = 0, + unsigned expected_none_zeros = default_expected_none_zeros) const + { + if (!rowsize) + rowsize = getColumnSize(); + RCMatrix result(rowsize, expected_none_zeros); + for (N i = 0; i < n; i++) { + for (Iterator it = row_begin(i); it; ++it) + result.fix_element(it.index(), i, it.value()); + } +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix getKtK() const + { + RCMatrix m = transpose(); + RCMatrix result(n, expected_none_zeros); + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + for (Iterator it_A = m.row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + assert(j < n); + T a = it_A.value(); + if (std::abs(a) > VECTOR_EPSILON) { + for (Iterator it_B = row_begin(j); it_B; ++it_B) { + // result.add_to_element(i,it_B.index(),it_B.value()*a); + double value = it_B.value() * a; + if (std::abs(value) > VECTOR_EPSILON) + result.add_to_element(i, it_B.index(), value); + } + } + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix operator*(const RCMatrix &m) const + { + RCMatrix result(n, expected_none_zeros); + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + assert(j < m.n); + T a = it_A.value(); + if (std::abs(a) > VECTOR_EPSILON) { + for (Iterator it_B = m.row_begin(j); it_B; ++it_B) { + // result.add_to_element(i,it_B.index(),it_B.value()*a); + double value = it_B.value() * a; + if (std::abs(value) > VECTOR_EPSILON) + result.add_to_element(i, it_B.index(), value); + } + } + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix sqrt() const + { + RCMatrix result(n, expected_none_zeros); + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + result.set_element(i, j, std::sqrt(it_A.value())); + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix operator*(const double k) const + { + RCMatrix result(n, expected_none_zeros); + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + result.add_to_element(i, j, it_A.value() * k); + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix applyKernel(const RCMatrix &kernel, const int nx, const int ny) const + { + RCMatrix result(n, expected_none_zeros); + // find center position of kernel (half of kernel size) + int kCols = kernel.n, kRows = kernel.n, rows = nx, cols = ny; + int kCenterX = kCols / 2; + int kCenterY = kRows / 2; + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + if (i >= rows) + break; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + if (j >= cols) + break; + for (int m = 0; m < kRows; ++m) { // kernel rows + int mm = kRows - 1 - m; // row index of flipped kernel + for (int n = 0; n < kCols; ++n) { // kernel columns + int nn = kCols - 1 - n; // column index of flipped kernel + // index of input signal, used for checking boundary + int ii = i + (m - kCenterY); + int jj = j + (n - kCenterX); + // ignore input samples which are out of bound + if (ii >= 0 && ii < rows && jj >= 0 && jj < cols) + result.add_to_element(i, j, (*this)(ii, jj) * kernel(mm, nn)); + } + } + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix applyHorizontalKernel(const RCMatrix &kernel, const int nx, const int ny) const + { + RCMatrix result(n, expected_none_zeros); + // find center position of kernel (half of kernel size) + int kCols = kernel.n, kRows = 1, rows = nx, cols = ny; + int kCenterX = kCols / 2; + int kCenterY = kRows / 2; + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + if (i >= rows) + break; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + if (j >= cols) + break; + for (int m = 0; m < kRows; ++m) { // kernel rows + int mm = kRows - 1 - m; // row index of flipped kernel + for (int n = 0; n < kCols; ++n) { // kernel columns + int nn = kCols - 1 - n; // column index of flipped kernel + // index of input signal, used for checking boundary + int ii = i + (m - kCenterY); + int jj = j + (n - kCenterX); + // ignore input samples which are out of bound + if (ii >= 0 && ii < rows && jj >= 0 && jj < cols) + result.add_to_element(i, j, (*this)(ii, jj) * kernel(mm, nn)); + } + } + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix applyVerticalKernel(const RCMatrix &kernel, const int nx, const int ny) const + { + RCMatrix result(n, expected_none_zeros); + // find center position of kernel (half of kernel size) + int kCols = 1, kRows = kernel.n, rows = nx, cols = ny; + int kCenterX = kCols / 2; + int kCenterY = kRows / 2; + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + if (i >= rows) + break; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + if (j >= cols) + break; + for (int m = 0; m < kRows; ++m) { // kernel rows + int mm = kRows - 1 - m; // row index of flipped kernel + for (int n = 0; n < kCols; ++n) { // kernel columns + int nn = kCols - 1 - n; // column index of flipped kernel + // index of input signal, used for checking boundary + int ii = i + (m - kCenterY); + int jj = j + (n - kCenterX); + // ignore input samples which are out of bound + if (ii >= 0 && ii < rows && jj >= 0 && jj < cols) + result.add_to_element(i, j, (*this)(ii, jj) * kernel(mm, nn)); + } + } + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix applySeparableKernel(const RCMatrix &kernelH, + const RCMatrix &kernelV, + const int nx, + const int ny) const + { + return applyHorizontalKernel(kernelH, nx, ny).applyVerticalKernel(kernelV, nx, ny); + } + + RCMatrix applySeparableKernelTwice(const RCMatrix &kernelH, + const RCMatrix &kernelV, + const int nx, + const int ny) const + { + return applySeparableKernel(kernelH, kernelV, nx, ny) + .applySeparableKernel(kernelH, kernelV, nx, ny); + } + + std::vector<T> operator*(const std::vector<T> &rhs) const + { + std::vector<T> result(n, 0.0); + multiply(rhs, result); +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + void multiply(const std::vector<T> &rhs, std::vector<T> &result) const + { + result.resize(n); + for (N i = 0; i < n; i++) { + T new_value = 0.0; + for (Iterator it = row_begin(i); it; ++it) { + N j_index = it.index(); + assert(j_index < rhs.size()); + new_value += rhs[j_index] * it.value(); + } + result[i] = new_value; + } + } + RCMatrix operator+(const RCMatrix &m) const + { + RCMatrix A(*this); + return std::move(A.add(m)); + } + RCMatrix &add(const RCMatrix &m) + { + if (m.n > n) + resize(m.n); + parallel_for(m.n) + { + N i = parallel_index; + for (Iterator it = m.row_begin(i); it; ++it) { + add_to_element(i, it.index(), it.value()); + } + } + parallel_end return *this; + } + RCMatrix operator-(const RCMatrix &m) const + { + RCMatrix A(*this); + return std::move(A.sub(m)); + } + RCMatrix &sub(const RCMatrix &m) + { + if (m.n > n) + resize(m.n); + parallel_for(m.n) + { + N i = parallel_index; + for (Iterator it = m.row_begin(i); it; ++it) { + add_to_element(i, it.index(), -it.value()); + } + } + parallel_end return *this; + } + RCMatrix &replace(const RCMatrix &m, int rowInd, int colInd) + { + if (m.n > n) + resize(m.n); + parallel_for(m.n) + { + N i = parallel_index; + for (Iterator it = m.row_begin(i); it; ++it) { + set_element(i + rowInd, it.index() + colInd, it.value()); + } + } + parallel_end return *this; + } + Real max_residual(const std::vector<T> &lhs, const std::vector<T> &rhs) const + { + std::vector<T> r = operator*(lhs); + Real max_residual = 0.0; + for (N i = 0; i < rhs.size(); i++) { + if (!empty(i)) + max_residual = std::max(max_residual, std::abs(r[i] - rhs[i])); + } + return max_residual; + } + std::vector<T> residual_vector(const std::vector<T> &lhs, const std::vector<T> &rhs) const + { + std::vector<T> result = operator*(lhs); + assert(result.size() == rhs.size()); + for (N i = 0; i < result.size(); i++) { + result[i] = std::abs(result[i] - rhs[i]); + } +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + T norm() const + { + T result(0.0); + for (N i = 0; i < n; ++i) { + for (Iterator it = row_begin(i); it; ++it) { + result = std::max(result, std::abs(it.value())); + } + } + return result; + } + + T norm_L2_sqr() const + { + T result(0.0); + for (N i = 0; i < n; ++i) { + for (Iterator it = row_begin(i); it; ++it) { + result += (it.value()) * (it.value()); + } + } + return result; + } + + void write_matlab(std::ostream &output, + unsigned int rows, + unsigned int columns, + const char *variable_name) + { + output << variable_name << "=sparse(["; + for (N i = 0; i < n; ++i) { + if (matrix[i]) { + const std::vector<N> &index = matrix[i]->index; + for (N j = 0; j < (N)index.size(); ++j) { + output << i + 1 << " "; + } + } + } + output << "],...\n ["; + for (N i = 0; i < n; ++i) { + if (matrix[i]) { + const std::vector<N> &index = matrix[i]->index; + for (N j = 0; j < (N)index.size(); ++j) { + output << index[j] + (offsets.empty() ? 0 : offsets[i]) + 1 << " "; + } + } + } + output << "],...\n ["; + for (N i = 0; i < n; ++i) { + if (matrix[i]) { + const std::vector<T> &value = matrix[i]->value; + for (N j = 0; j < value.size(); ++j) { + output << value[j] << " "; + } + } + } + output << "], " << rows << ", " << columns << ");" << std::endl; + }; + void export_matlab(std::string filename, std::string name) + { + // Export this matrix + std::ofstream file; + file.open(filename.c_str()); + write_matlab(file, n, getColumnSize(), name.c_str()); + file.close(); + } + void print_readable(std::string name, bool printNonZero = true) + { + std::cout << name << " \n"; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (printNonZero) { + if ((*this)(i, j) == 0) { + std::cout << " ."; + continue; + } + } + else { + if ((*this)(i, j) == 0) { + continue; + } + } + + if ((*this)(i, j) >= 0) + std::cout << " "; + std::cout << " " << (*this)(i, j); + } + std::cout << " \n"; + } + } + /// + N n; + N expected_none_zeros; + std::vector<RowEntry *> matrix; + std::vector<int> offsets; +}; + +template<class N, class T> +static inline RCMatrix<N, T> operator*(const std::vector<T> &diagonal, const RCMatrix<N, T> &A) +{ + RCMatrix<N, T> result(A); + parallel_for(result.n) + { + N row(parallel_index); + for (auto it = result.dynamic_row_begin(row); it; ++it) { + it.setValue(it.value() * diagonal[row]); + } + } + parallel_end return std::move(result); +} + +template<class N, class T> struct RCFixedMatrix { + std::vector<N> rowstart; + std::vector<N> index; + std::vector<T> value; + N n; + N max_rowlength; + // + RCFixedMatrix() : n(0), max_rowlength(0) + { + } + RCFixedMatrix(const RCMatrix<N, T> &matrix) + { + n = matrix.n; + rowstart.resize(n + 1); + rowstart[0] = 0; + max_rowlength = 0; + for (N i = 0; i < n; i++) { + if (!matrix.empty(i)) { + rowstart[i + 1] = rowstart[i] + matrix.row_nonzero_size(i); + max_rowlength = std::max(max_rowlength, rowstart[i + 1] - rowstart[i]); + } + else { + rowstart[i + 1] = rowstart[i]; + } + } + value.resize(rowstart[n]); + index.resize(rowstart[n]); + N j = 0; + for (N i = 0; i < n; i++) { + for (typename RCMatrix<N, T>::Iterator it = matrix.row_begin(i); it; ++it) { + value[j] = it.value(); + index[j] = it.index(); + ++j; + } + } + } + class Iterator : std::iterator<std::input_iterator_tag, T> { + public: + Iterator(N start, N end, const std::vector<N> &index, const std::vector<T> &value) + : index_array(index), value_array(value), k(start), start(start), end(end) + { + } + operator bool() const + { + return k < end; + } + Iterator &operator++() + { + ++k; + return *this; + } + T value() const + { + return value_array[k]; + } + N index() const + { + return index_array[k]; + } + N size() const + { + return end - start; + } + + private: + const std::vector<N> &index_array; + const std::vector<T> &value_array; + N k, start, end; + }; + Iterator row_begin(N n) const + { + return Iterator(rowstart[n], rowstart[n + 1], index, value); + } + std::vector<T> operator*(const std::vector<T> &rhs) const + { + std::vector<T> result(n, 0.0); + multiply(rhs, result); +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + void multiply(const std::vector<T> &rhs, std::vector<T> &result) const + { + result.resize(n); + parallel_for(n) + { + N i = parallel_index; + T new_value = 0.0; + for (Iterator it = row_begin(i); it; ++it) { + N j_index = it.index(); + assert(j_index < rhs.size()); + new_value += rhs[j_index] * it.value(); + } + result[i] = new_value; + } + parallel_end + } + RCMatrix<N, T> operator*(const RCFixedMatrix &m) const + { + RCMatrix<N, T> result(n, max_rowlength); + // Run in parallel + parallel_for(result.n) + { + N i = parallel_index; + for (Iterator it_A = row_begin(i); it_A; ++it_A) { + N j = it_A.index(); + assert(j < m.n); + T a = it_A.value(); + if (std::abs(a) > VECTOR_EPSILON) { + for (Iterator it_B = m.row_begin(j); it_B; ++it_B) { + result.add_to_element(i, it_B.index(), it_B.value() * a); + } + } + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } + + RCMatrix<N, T> toRCMatrix() const + { + RCMatrix<N, T> result(n, 0); + parallel_for(n) + { + N i = parallel_index; + N size = rowstart[i + 1] - rowstart[i]; + result.matrix[i] = new typename RCMatrix<N, T>::RowEntry; + result.matrix[i]->index.resize(size); + result.matrix[i]->value.resize(size); + for (N j = 0; j < size; j++) { + result.matrix[i]->index[j] = index[rowstart[i] + j]; + result.matrix[i]->value[j] = value[rowstart[i] + j]; + } + } + parallel_end +#if MANTA_USE_CPP11 == 1 + return std::move(result); +#else + return result; +#endif + } +}; + +typedef RCMatrix<int, Real> Matrix; +typedef RCFixedMatrix<int, Real> FixedMatrix; + +} // namespace Manta + +#undef parallel_for +#undef parallel_end + +#undef parallel_block +#undef do_parallel +#undef do_end +#undef block_end + +#endif diff --git a/extern/mantaflow/helper/util/simpleimage.cpp b/extern/mantaflow/helper/util/simpleimage.cpp new file mode 100644 index 00000000000..9846fa5bd96 --- /dev/null +++ b/extern/mantaflow/helper/util/simpleimage.cpp @@ -0,0 +1,312 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Simple image IO + * + ******************************************************************************/ + +#include "vectorbase.h" +#include "simpleimage.h" + +namespace Manta { + +// write rectangle to ppm +bool SimpleImage::writePpm( + std::string filename, int minx, int miny, int maxx, int maxy, bool invertXY) +{ + int w = maxx - minx; + int h = maxy - miny; + + if (w <= 0 || h <= 0 || w > mSize[0] || h > mSize[1]) { + errMsg("SimpleImage::WritePPM Invalid rect: w=" + << w << ", h=" << h << ", size=" << mSize[0] << "," << mSize[1] << " min/max: " << minx + << "," << miny << " to " << maxx << "," << maxy << ", resetting... "); + minx = miny = 0; + maxx = mSize[0] - 1; + maxy = mSize[1] - 1; + w = mSize[0] - 1; + h = mSize[1] - 1; + } + + FILE *fp = fopen(filename.c_str(), "wb"); + if (fp == NULL) { + errMsg("SimpleImage::WritePPM Unable to open '" << filename << "' for writing"); + return false; + } + fprintf(fp, "P6\n%d %d\n255\n", w, h); + + int pixCnt = 0; + for (int j = maxy - 1; j >= miny; j--) + for (int i = minx; i < maxx; i++) { + unsigned char col[3]; + for (int l = 0; l < 3; l++) { + float val; + if (invertXY) + val = (float)get(j, i)[l]; + else + val = (float)get(i, j)[l]; + + val = clamp(val, (float)0., (float)1.); + col[l] = (unsigned char)(255. * val); + } + // col[1] = col[2] = col[0]; + // if (fwrite(col,1,3, fp) != 3) errMsg("SimpleImage::writePpm fwrite failed"); + fwrite(col, 1, 3, fp); + pixCnt++; + // fprintf(stderr,"%d %d %d \n",col[0],i,j); + } + + fclose(fp); + // debMsg("WritePPM Wrote '"<<filename<<"', region="<<minx<<","<<miny<<" to + // "<<maxx<<","<<maxy<<"; "<<pixCnt, 1); + + return true; +} + +bool SimpleImage::writePpm(std::string filename) +{ + return writePpm(filename, 0, 0, getSize()[0], getSize()[1]); +} + +// read in a ppm file, and init the image accordingly +bool SimpleImage::initFromPpm(std::string filename) +{ + // maximum length of a line of text + const int MAXLINE = 1024; + + int filetype = 0; + enum { PGM, PPM }; // possible file types + + FILE *fp; + char line[MAXLINE]; + int size, rowsize; + + // Read in file type + fp = fopen(filename.c_str(), "rb"); + if (!fp) { + if (mAbortOnError) + debMsg("SimpleImage Error - unable to open file '" << filename << "' for reading", 1); + return 0; + } + + // 1st line: PPM or PGM + if (fgets(line, MAXLINE, fp) == NULL) { + if (mAbortOnError) + debMsg("SimpleImage::initFromPpm fgets failed", 1); + return 0; + } + + if (line[1] == '5') + filetype = PGM; + else if (line[1] == '6') + filetype = PPM; + else { + if (mAbortOnError) + debMsg("SimpleImage Error: need PPM or PGM file as input!", 1); + return 0; + } + + // Read in width and height, & allocate space + // 2nd line: width height + if (fgets(line, MAXLINE, fp) == NULL) { + if (mAbortOnError) + errMsg("SimpleImage::initFromPpm fgets failed"); + return 0; + } + int windW = 0, windH = 0; // size of the window on the screen + int intsFound = sscanf(line, "%d %d", &windW, &windH); + if (intsFound == 1) { + // only X found, search on next line as well for Y... + if (sscanf(line, "%d", &windH) != 1) { + if (mAbortOnError) + errMsg("initFromPpm Ppm dimensions not found!" << windW << "," << windH); + return 0; + } + else { + // ok, found 2 lines + // debMsg("initFromPpm Ppm dimensions found!"<<windW<<","<<windH, 1); + } + } + else if (intsFound == 2) { + // ok! + } + else { + if (mAbortOnError) + errMsg("initFromPpm Ppm dimensions not found at all!" << windW << "," << windH); + return 0; + } + + if (filetype == PGM) { + size = windH * windW; // greymap: 1 byte per pixel + rowsize = windW; + } + else { + // filetype == PPM + size = windH * windW * 3; // pixmap: 3 bytes per pixel + rowsize = windW * 3; + } + + unsigned char *pic = new unsigned char[size]; // (GLubyte *)malloc (size); + + // Read in maximum value (ignore) , could be scanned with sscanf as well, but this should be + // 255... 3rd line + if (fgets(line, MAXLINE, fp) == NULL) { + if (mAbortOnError) + errMsg("SimpleImage::initFromPpm fgets failed"); + return 0; + } + + // Read in the pixel array row-by-row: 1st row = top scanline */ + unsigned char *ptr = NULL; + ptr = &pic[(windH - 1) * rowsize]; + for (int i = windH; i > 0; i--) { + assertMsg(fread((void *)ptr, 1, rowsize, fp) == rowsize, + "SimpleImage::initFromPpm couldn't read data"); + ptr -= rowsize; + } + + // init image + this->init(windW, windH); + if (filetype == PGM) { + // grayscale + for (int i = 0; i < windW; i++) { + for (int j = 0; j < windH; j++) { + double r = (double)pic[(j * windW + i) * 1 + 0] / 255.; + (*this)(i, j) = Vec3(r, r, r); + } + } + } + else { + // convert grid to RGB vec's + for (int i = 0; i < windW; i++) { + for (int j = 0; j < windH; j++) { + // return mpData[y*mSize[0]+x]; + double r = (double)pic[(j * windW + i) * 3 + 0] / 255.; + double g = (double)pic[(j * windW + i) * 3 + 1] / 255.; + double b = (double)pic[(j * windW + i) * 3 + 2] / 255.; + + //(*this)(i,j) = Vec3(r,g,b); + + // RGB values have to be rotated to get the right colors!? + // this might also be an artifact of photoshop export...? + (*this)(i, j) = Vec3(g, b, r); + } + } + } + + delete[] pic; + fclose(fp); + return 1; +} + +// check index is valid +bool SimpleImage::indexIsValid(int i, int j) +{ + if (i < 0) + return false; + if (j < 0) + return false; + if (i >= mSize[0]) + return false; + if (j >= mSize[1]) + return false; + return true; +} + +}; // namespace Manta + +//***************************************************************************** + +#include "grid.h" +namespace Manta { + +// simple shaded output , note requires grid functionality! +static void gridPrecompLight(const Grid<Real> &density, Grid<Real> &L, Vec3 light = Vec3(1, 1, 1)) +{ + FOR_IJK(density) + { + Vec3 n = getGradient(density, i, j, k) * -1.; + normalize(n); + + Real d = dot(light, n); + L(i, j, k) = d; + } +} + +// simple shading with pre-computed gradient +static inline void shadeCell( + Vec3 &dst, int shadeMode, Real src, Real light, int depthPos, Real depthInv) +{ + switch (shadeMode) { + + case 1: { + // surfaces + Vec3 ambient = Vec3(0.1, 0.1, 0.1); + Vec3 diffuse = Vec3(0.9, 0.9, 0.9); + Real alpha = src; + + // different color for depth? + diffuse[0] *= ((Real)depthPos * depthInv) * 0.7 + 0.3; + diffuse[1] *= ((Real)depthPos * depthInv) * 0.7 + 0.3; + + Vec3 col = ambient + diffuse * light; + + // img( 0+i, j ) = (1.-alpha) * img( 0+i, j ) + alpha * col; + dst = (1. - alpha) * dst + alpha * col; + } break; + + default: { + // volumetrics / smoke + dst += depthInv * Vec3(src, src, src); + } break; + } +} + +//! helper to project a grid intro an image (used for ppm export and GUI displauy) +void projectImg(SimpleImage &img, const Grid<Real> &val, int shadeMode = 0, Real scale = 1.) +{ + Vec3i s = val.getSize(); + Vec3 si = Vec3(1. / (Real)s[0], 1. / (Real)s[1], 1. / (Real)s[2]); + + // init image size + int imgSx = s[0]; + if (val.is3D()) + imgSx += s[2] + s[0]; // mult views in 3D + img.init(imgSx, std::max(s[0], std::max(s[1], s[2]))); + + // precompute lighting + Grid<Real> L(val); + gridPrecompLight(val, L, Vec3(1, 1, 1)); + + FOR_IJK(val) + { + Vec3i idx(i, j, k); + shadeCell(img(0 + i, j), shadeMode, val(idx), L(idx), k, si[2]); + } + + if (val.is3D()) { + + FOR_IJK(val) + { + Vec3i idx(i, j, k); + shadeCell(img(s[0] + k, j), shadeMode, val(idx), L(idx), i, si[0]); + } + + FOR_IJK(val) + { + Vec3i idx(i, j, k); + shadeCell(img(s[0] + s[2] + i, k), shadeMode, val(idx), L(idx), j, si[1]); + } + + } // 3d + + img.mapRange(1. / scale); +} + +}; // namespace Manta diff --git a/extern/mantaflow/helper/util/simpleimage.h b/extern/mantaflow/helper/util/simpleimage.h new file mode 100644 index 00000000000..d7e88b83f74 --- /dev/null +++ b/extern/mantaflow/helper/util/simpleimage.h @@ -0,0 +1,205 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2014 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Simple image IO + * + ******************************************************************************/ + +#ifndef MANTA_SIMPLEIMAGE_H +#define MANTA_SIMPLEIMAGE_H + +#include <stdio.h> +#include <string.h> +#include <stdlib.h> + +#include "manta.h" +#include "vectorbase.h" + +namespace Manta { + +//***************************************************************************** +// simple 2d image class +// template<class Scalar> +class SimpleImage { + public: + // cons/des + SimpleImage() : mSize(-1), mpData(NULL), mAbortOnError(true){}; + virtual ~SimpleImage() + { + if (mpData) + delete[] mpData; + }; + + //! set to constant + void reset(Real val = 0.) + { + const Vec3 v = Vec3(val); + for (int i = 0; i < mSize[0] * mSize[1]; i++) + mpData[i] = v; + } + //! init memory & reset to zero + void init(int x, int y) + { + mSize = Vec3i(x, y, 0); + mpData = new Vec3[x * y]; + reset(); + }; + + inline bool checkIndex(int x, int y) + { + if ((x < 0) || (y < 0) || (x > mSize[0] - 1) || (y > mSize[1] - 1)) { + errMsg("SimpleImage::operator() Invalid access to " << x << "," << y << ", size=" << mSize); + return false; + } + return true; + } + + // access element + inline Vec3 &operator()(int x, int y) + { + DEBUG_ONLY(checkIndex(x, y)); + return mpData[y * mSize[0] + x]; + }; + inline Vec3 &get(int x, int y) + { + return (*this)(x, y); + } + inline Vec3 &getMap(int x, int y, int z, int axis) + { + int i = x, j = y; + if (axis == 1) + j = z; + if (axis == 0) { + i = y; + j = z; + } + return get(i, j); + } + + // output as string, debug + std::string toString() + { + std::ostringstream out; + + for (int j = 0; j < mSize[1]; j++) { + for (int i = 0; i < mSize[0]; i++) { + // normal zyx order */ + out << (*this)(i, j); + out << " "; + } + // if (format) + out << std::endl; + } + + return out.str(); + } + + // multiply all values by f + void add(Vec3 f) + { + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + get(i, j) += f; + } + } + // multiply all values by f + void multiply(Real f) + { + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + get(i, j) *= f; + } + } + // map 0-f to 0-1 range, clamp + void mapRange(Real f) + { + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + get(i, j) /= f; + for (int c = 0; c < 3; ++c) + get(i, j)[c] = clamp(get(i, j)[c], (Real)0., (Real)1.); + } + } + + // normalize max values + void normalizeMax() + { + Real max = normSquare(get(0, 0)); + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + if (normSquare(get(i, j)) > max) + max = normSquare(get(i, j)); + } + max = sqrt(max); + Real invMax = 1. / max; + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + get(i, j) *= invMax; + } + }; + + // normalize min and max values + void normalizeMinMax() + { + Real max = normSquare(get(0, 0)); + Real min = max; + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + if (normSquare(get(i, j)) > max) + max = normSquare(get(i, j)); + if (normSquare(get(i, j)) < min) + min = normSquare(get(i, j)); + } + max = sqrt(max); + min = sqrt(min); + Real factor = 1. / (max - min); + for (int j = 0; j < mSize[1]; j++) + for (int i = 0; i < mSize[0]; i++) { + get(i, j) -= min; + get(i, j) *= factor; + } + }; + + void setAbortOnError(bool set) + { + mAbortOnError = set; + } + + // ppm in/output + + // write whole image + bool writePpm(std::string filename); + // write rectangle to ppm + bool writePpm( + std::string filename, int minx, int miny, int maxx, int maxy, bool invertXY = false); + // read in a ppm file, and init the image accordingly + bool initFromPpm(std::string filename); + + // check index is valid + bool indexIsValid(int i, int j); + + //! access + inline Vec3i getSize() const + { + return mSize; + } + + protected: + //! size + Vec3i mSize; + //! data + Vec3 *mpData; + // make errors fatal, or continue? + bool mAbortOnError; + +}; // SimpleImage + +}; // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/solvana.h b/extern/mantaflow/helper/util/solvana.h new file mode 100644 index 00000000000..9dc1ec83654 --- /dev/null +++ b/extern/mantaflow/helper/util/solvana.h @@ -0,0 +1,214 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Analytical solutions to some problems + * generated using MATLAB symbolic math ccode + * + ******************************************************************************/ + +#ifndef _SOLVANA_H +#define _SOLVANA_H + +//! solves the equation [e1 e2 e3; 1 1 1]*x = g using least squares +inline void SolveOverconstraint34(float e1x, + float e1y, + float e1z, + float e2x, + float e2y, + float e2z, + float e3x, + float e3y, + float e3z, + float g1, + float g2, + float g3, + float &x1, + float &x2, + float &x3) +{ + float e1x2 = e1x * e1x, e1y2 = e1y * e1y, e1z2 = e1z * e1z; + float e2x2 = e2x * e2x, e2y2 = e2y * e2y, e2z2 = e2z * e2z; + float e3x2 = e3x * e3x, e3y2 = e3y * e3y, e3z2 = e3z * e3z; + float e1xy = e1x * e1y, e1xz = e1x * e1z, e1yz = e1y * e1z; + float e2xy = e2x * e2y, e2xz = e2x * e2z, e2yz = e2y * e2z; + float e3xy = e3x * e3y, e3xz = e3x * e3z, e3yz = e3y * e3z; + float e12x = e1x * e2x, e12y = e1y * e2y, e12z = e1z * e2z; + float e13x = e1x * e3x, e13y = e1y * e3y, e13z = e1z * e3z; + float e23x = e2x * e3x, e23y = e2y * e3y, e23z = e2z * e3z; + float t1543 = e3y2 * e2x2; + float t1544 = e3x2 * e2y2; + float t1545 = e3z2 * e2x2; + float t1546 = e3x2 * e2z2; + float t1547 = e3z2 * e2y2; + float t1548 = e3y2 * e2z2; + float t1549 = e2y2 * e1x2; + float t1550 = e2x2 * e1y2; + float t1551 = e2z2 * e1x2; + float t1552 = e2x2 * e1z2; + float t1553 = e2z2 * e1y2; + float t1554 = e2y2 * e1z2; + float t1555 = e3y2 * e1x2; + float t1556 = e3x2 * e1y2; + float t1557 = e3z2 * e1x2; + float t1558 = e3x2 * e1z2; + float t1559 = e3z2 * e1y2; + float t1560 = e3y2 * e1z2; + float t1561 = e3z2 * e2y2 * e1x2; + float t1562 = e3y2 * e2z2 * e1x2; + float t1563 = e3z2 * e2x2 * e1y2; + float t1564 = e3x2 * e2z2 * e1y2; + float t1565 = e3y2 * e2x2 * e1z2; + float t1566 = e3x2 * e2y2 * e1z2; + float t1567 = e1xy * e2x * e3y * 2.0; + float t1568 = e1xy * e2y * e3x * 2.0; + float t1569 = e1xz * e2x * e3z * 2.0; + float t1570 = e1xz * e2z * e3x * 2.0; + float t1571 = e1yz * e2y * e3z * 2.0; + float t1572 = e1yz * e2z * e3y * 2.0; + float t1573 = e1x * e2xy * e3y * 2.0; + float t1574 = e1y * e2xy * e3x * 2.0; + float t1575 = e1x * e2xz * e3z * 2.0; + float t1576 = e1z * e2xz * e3x * 2.0; + float t1577 = e1y * e2yz * e3z * 2.0; + float t1578 = e1z * e2yz * e3y * 2.0; + float t1579 = e1x * e2y * e3xy * 2.0; + float t1580 = e1y * e2x * e3xy * 2.0; + float t1581 = e1x * e2z * e3xz * 2.0; + float t1582 = e1z * e2x * e3xz * 2.0; + float t1583 = e1y * e2z * e3yz * 2.0; + float t1584 = e1z * e2y * e3yz * 2.0; + float t1585 = e1xy * e2xz * e3yz * 2.0; + float t1586 = e1xy * e2yz * e3xz * 2.0; + float t1587 = e1xz * e2xy * e3yz * 2.0; + float t1588 = e1xz * e2yz * e3xy * 2.0; + float t1589 = e1yz * e2xy * e3xz * 2.0; + float t1590 = e1yz * e2xz * e3xy * 2.0; + float t1596 = e12x * e3y2 * 2.0; + float t1597 = e13x * e2y2 * 2.0; + float t1598 = e23x * e1y2 * 2.0; + float t1599 = e12x * e3z2 * 2.0; + float t1600 = e13x * e2z2 * 2.0; + float t1601 = e12y * e3x2 * 2.0; + float t1602 = e13y * e2x2 * 2.0; + float t1603 = e23y * e1x2 * 2.0; + float t1604 = e23x * e1z2 * 2.0; + float t1605 = e12y * e3z2 * 2.0; + float t1606 = e13y * e2z2 * 2.0; + float t1607 = e12z * e3x2 * 2.0; + float t1608 = e13z * e2x2 * 2.0; + float t1609 = e23z * e1x2 * 2.0; + float t1610 = e23y * e1z2 * 2.0; + float t1611 = e12z * e3y2 * 2.0; + float t1612 = e13z * e2y2 * 2.0; + float t1613 = e23z * e1y2 * 2.0; + float t1614 = e1xy * e2xy * 2.0; + float t1615 = e1xz * e2xz * 2.0; + float t1616 = e1yz * e2yz * 2.0; + float t1617 = e1xy * e3xy * 2.0; + float t1618 = e1xz * e3xz * 2.0; + float t1619 = e1yz * e3yz * 2.0; + float t1620 = e2xy * e3xy * 2.0; + float t1621 = e2xz * e3xz * 2.0; + float t1622 = e2yz * e3yz * 2.0; + float t1623 = e1xy * e2xy * e3z2 * 2.0; + float t1624 = e1xz * e2xz * e3y2 * 2.0; + float t1625 = e1yz * e2yz * e3x2 * 2.0; + float t1626 = e1xy * e3xy * e2z2 * 2.0; + float t1627 = e1xz * e3xz * e2y2 * 2.0; + float t1628 = e1yz * e3yz * e2x2 * 2.0; + float t1629 = e2xy * e3xy * e1z2 * 2.0; + float t1630 = e2xz * e3xz * e1y2 * 2.0; + float t1631 = e2yz * e3yz * e1x2 * 2.0; + float t1591 = t1550 + t1551 + t1560 + t1543 + t1552 + t1561 + t1570 + t1544 + t1553 + t1562 + + t1571 + t1580 + t1545 + t1554 + t1563 + t1572 + t1581 + t1590 + t1546 + t1555 + + t1564 + t1573 + t1582 + t1547 + t1556 + t1565 + t1574 + t1583 + t1548 + t1557 + + t1566 + t1575 + t1584 + t1549 + t1558 + t1567 + t1576 + t1585 + t1559 + t1568 + + t1577 + t1586 + t1569 + t1578 + t1587 - t1596 + t1579 + t1588 - t1597 + t1589 - + t1598 - t1599 - t1600 - t1601 - t1610 - t1602 - t1611 - t1620 - t1603 - t1612 - + t1621 - t1630 - t1604 - t1613 - t1622 - t1631 - t1605 - t1614 - t1623 - t1606 - + t1615 - t1624 - t1607 - t1616 - t1625 - t1608 - t1617 - t1626 - t1609 - t1618 - + t1627 - t1619 - t1628 - t1629; + float t1592 = 1.0 / t1591; + float t1635 = e13x * e2y2; + float t1636 = e13x * e2z2; + float t1637 = e13y * e2x2; + float t1638 = e13y * e2z2; + float t1639 = e13z * e2x2; + float t1640 = e13z * e2y2; + float t1653 = e23x * 2.0; + float t1654 = e23y * 2.0; + float t1655 = e23z * 2.0; + float t1641 = e3x2 + e3z2 + e3y2 + e2y2 + t1543 + e2z2 + t1544 + e2x2 + t1545 + t1546 + t1547 + + t1548 - t1620 - t1621 - t1622 - t1653 - t1654 - t1655; + float t1642 = e12x * e3y2; + float t1643 = e12x * e3z2; + float t1644 = e12y * e3x2; + float t1645 = e12y * e3z2; + float t1646 = e12z * e3x2; + float t1647 = e12z * e3y2; + float t1656 = e1x * e2y * e3xy; + float t1657 = e1y * e2x * e3xy; + float t1658 = e1x * e2z * e3xz; + float t1659 = e1z * e2x * e3xz; + float t1660 = e1y * e2z * e3yz; + float t1661 = e1z * e2y * e3yz; + float t1648 = e3x2 + e3z2 + e3y2 - e13x - e13y - e13z + e12x - e23y + e12y + t1642 - e23z - + t1660 + e12z + t1643 - t1661 + t1644 + t1645 + t1646 + t1647 - t1656 - t1657 - + e23x - t1658 - t1659; + float t1679 = e1x * e2xy * e3y; + float t1680 = e1y * e2xy * e3x; + float t1681 = e1x * e2xz * e3z; + float t1682 = e1z * e2xz * e3x; + float t1683 = e1y * e2yz * e3z; + float t1684 = e1z * e2yz * e3y; + float t1652 = e2y2 + e2z2 + e2x2 + e13x + e13y + e13z + t1640 - e12x - e23y - e12y - e23z - + e12z + t1635 - t1680 + t1636 - t1681 + t1637 - t1682 + t1638 - t1683 + t1639 - + t1684 - e23x - t1679; + float t1662 = e23x * e1y2; + float t1663 = e23y * e1x2; + float t1664 = e23x * e1z2; + float t1665 = e23z * e1x2; + float t1666 = e23y * e1z2; + float t1667 = e23z * e1y2; + float t1670 = e1xy * e2x * e3y; + float t1671 = e1xy * e2y * e3x; + float t1672 = e1xz * e2x * e3z; + float t1673 = e1xz * e2z * e3x; + float t1674 = e1yz * e2y * e3z; + float t1675 = e1yz * e2z * e3y; + float t1668 = e1x2 + e1y2 + e1z2 - e13x - e13y - e13z - e12x + e23y - e12y + e23z - e12z - + t1670 + t1662 - t1671 + t1663 - t1672 + t1664 - t1673 + t1665 - t1674 + t1666 - + t1675 + e23x + t1667; + float t1676 = e13x * 2.0; + float t1677 = e13y * 2.0; + float t1678 = e13z * 2.0; + float t1669 = e3x2 + e3z2 + e3y2 + t1560 + e1x2 + t1555 + e1y2 + t1556 + e1z2 + t1557 + t1558 + + t1559 - t1617 - t1618 - t1619 - t1676 - t1677 - t1678; + float t1686 = e12x * 2.0; + float t1687 = e12y * 2.0; + float t1688 = e12z * 2.0; + float t1685 = t1550 + t1551 + e2y2 + t1552 + e2z2 + t1553 + e2x2 + t1554 + e1x2 + e1y2 + e1z2 + + t1549 - t1614 - t1615 - t1616 - t1686 - t1687 - t1688; + x1 = -g2 * (-e1y * t1592 * t1641 + e2y * t1592 * t1648 + e3y * t1592 * t1652) - + g3 * (-e1z * t1592 * t1641 + e2z * t1592 * t1648 + e3z * t1592 * t1652) - + g1 * (-e1x * t1592 * t1641 + e2x * t1592 * t1648 + + e3x * t1592 * + (e2y2 + e2z2 + e2x2 + e13x + e13y + e13z + t1640 + t1635 + t1636 + t1637 + t1638 + + t1639 - e12x - e12y - e12z - e23x - e23y - e23z - e1x * e2xy * e3y - + e1y * e2xy * e3x - e1x * e2xz * e3z - e1z * e2xz * e3x - e1y * e2yz * e3z - + e1z * e2yz * e3y)); + x2 = -g1 * (e1x * t1592 * t1648 - e2x * t1592 * t1669 + e3x * t1592 * t1668) - + g2 * (e1y * t1592 * t1648 - e2y * t1592 * t1669 + e3y * t1592 * t1668) - + g3 * (e1z * t1592 * t1648 - e2z * t1592 * t1669 + e3z * t1592 * t1668); + x3 = -g1 * (e1x * t1592 * t1652 + e2x * t1592 * t1668 - e3x * t1592 * t1685) - + g2 * (e1y * t1592 * t1652 + e2y * t1592 * t1668 - e3y * t1592 * t1685) - + g3 * (e1z * t1592 * t1652 + e2z * t1592 * t1668 - e3z * t1592 * t1685); +} + +#endif
\ No newline at end of file diff --git a/extern/mantaflow/helper/util/vector4d.cpp b/extern/mantaflow/helper/util/vector4d.cpp new file mode 100644 index 00000000000..d342df607f5 --- /dev/null +++ b/extern/mantaflow/helper/util/vector4d.cpp @@ -0,0 +1,50 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Basic vector class + * + ******************************************************************************/ + +#include "vector4d.h" + +using namespace std; + +namespace Manta { + +template<> const Vector4D<int> Vector4D<int>::Zero(0, 0, 0, 0); +template<> const Vector4D<float> Vector4D<float>::Zero(0.f, 0.f, 0.f, 0.f); +template<> const Vector4D<double> Vector4D<double>::Zero(0., 0., 0., 0.); +template<> +const Vector4D<float> Vector4D<float>::Invalid(numeric_limits<float>::quiet_NaN(), + numeric_limits<float>::quiet_NaN(), + numeric_limits<float>::quiet_NaN(), + numeric_limits<float>::quiet_NaN()); +template<> +const Vector4D<double> Vector4D<double>::Invalid(numeric_limits<double>::quiet_NaN(), + numeric_limits<double>::quiet_NaN(), + numeric_limits<double>::quiet_NaN(), + numeric_limits<double>::quiet_NaN()); +template<> bool Vector4D<float>::isValid() const +{ + return !c_isnan(x) && !c_isnan(y) && !c_isnan(z) && !c_isnan(t); +} +template<> bool Vector4D<double>::isValid() const +{ + return !c_isnan(x) && !c_isnan(y) && !c_isnan(z) && !c_isnan(t); +} + +//! Specialization for readable ints +template<> std::string Vector4D<int>::toString() const +{ + char buf[256]; + snprintf(buf, 256, "[%d,%d,%d,%d]", (*this)[0], (*this)[1], (*this)[2], (*this)[3]); + return std::string(buf); +} + +} // namespace Manta diff --git a/extern/mantaflow/helper/util/vector4d.h b/extern/mantaflow/helper/util/vector4d.h new file mode 100644 index 00000000000..c3d72ac8aff --- /dev/null +++ b/extern/mantaflow/helper/util/vector4d.h @@ -0,0 +1,515 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * 4D vector class + * + ******************************************************************************/ + +#ifndef _VECTOR4D_H +#define _VECTOR4D_H + +#include "vectorbase.h" + +namespace Manta { + +//! Basic inlined vector class +template<class S> class Vector4D { + public: + //! Constructor + inline Vector4D() : x(0), y(0), z(0), t(0) + { + } + + //! Copy-Constructor + inline Vector4D(const Vector4D<S> &v) : x(v.x), y(v.y), z(v.z), t(v.t) + { + } + + //! Copy-Constructor + inline Vector4D(const float *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3]) + { + } + + //! Copy-Constructor + inline Vector4D(const double *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3]) + { + } + + //! Construct a vector from one S + inline Vector4D(S v) : x(v), y(v), z(v), t(v) + { + } + + //! Construct a vector from three Ss + inline Vector4D(S vx, S vy, S vz, S vw) : x(vx), y(vy), z(vz), t(vw) + { + } + + // Operators + + //! Assignment operator + inline const Vector4D<S> &operator=(const Vector4D<S> &v) + { + x = v.x; + y = v.y; + z = v.z; + t = v.t; + return *this; + } + //! Assignment operator + inline const Vector4D<S> &operator=(S s) + { + x = y = z = t = s; + return *this; + } + //! Assign and add operator + inline const Vector4D<S> &operator+=(const Vector4D<S> &v) + { + x += v.x; + y += v.y; + z += v.z; + t += v.t; + return *this; + } + //! Assign and add operator + inline const Vector4D<S> &operator+=(S s) + { + x += s; + y += s; + z += s; + t += s; + return *this; + } + //! Assign and sub operator + inline const Vector4D<S> &operator-=(const Vector4D<S> &v) + { + x -= v.x; + y -= v.y; + z -= v.z; + t -= v.t; + return *this; + } + //! Assign and sub operator + inline const Vector4D<S> &operator-=(S s) + { + x -= s; + y -= s; + z -= s; + t -= s; + return *this; + } + //! Assign and mult operator + inline const Vector4D<S> &operator*=(const Vector4D<S> &v) + { + x *= v.x; + y *= v.y; + z *= v.z; + t *= v.t; + return *this; + } + //! Assign and mult operator + inline const Vector4D<S> &operator*=(S s) + { + x *= s; + y *= s; + z *= s; + t *= s; + return *this; + } + //! Assign and div operator + inline const Vector4D<S> &operator/=(const Vector4D<S> &v) + { + x /= v.x; + y /= v.y; + z /= v.z; + t /= v.t; + return *this; + } + //! Assign and div operator + inline const Vector4D<S> &operator/=(S s) + { + x /= s; + y /= s; + z /= s; + t /= s; + return *this; + } + //! Negation operator + inline Vector4D<S> operator-() const + { + return Vector4D<S>(-x, -y, -z, -t); + } + + //! Get smallest component + // inline S min() const { return ( x<y ) ? ( ( x<z ) ? x:z ) : ( ( y<z ) ? y:z ); } + //! Get biggest component + // inline S max() const { return ( x>y ) ? ( ( x>z ) ? x:z ) : ( ( y>z ) ? y:z ); } + + //! Test if all components are zero + inline bool empty() + { + return x == 0 && y == 0 && z == 0 && t == 0; + } + + //! access operator + inline S &operator[](unsigned int i) + { + return value[i]; + } + //! constant access operator + inline const S &operator[](unsigned int i) const + { + return value[i]; + } + + //! debug output vector to a string + std::string toString() const; + + //! test if nans are present + bool isValid() const; + + //! actual values + union { + S value[4]; + struct { + S x; + S y; + S z; + S t; + }; + struct { + S X; + S Y; + S Z; + S T; + }; + }; + + // zero element + static const Vector4D<S> Zero, Invalid; + + protected: +}; + +//************************************************************************ +// Additional operators +//************************************************************************ + +//! Addition operator +template<class S> inline Vector4D<S> operator+(const Vector4D<S> &v1, const Vector4D<S> &v2) +{ + return Vector4D<S>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.t + v2.t); +} +//! Addition operator +template<class S, class S2> inline Vector4D<S> operator+(const Vector4D<S> &v, S2 s) +{ + return Vector4D<S>(v.x + s, v.y + s, v.z + s, v.t + s); +} +//! Addition operator +template<class S, class S2> inline Vector4D<S> operator+(S2 s, const Vector4D<S> &v) +{ + return Vector4D<S>(v.x + s, v.y + s, v.z + s, v.t + s); +} + +//! Subtraction operator +template<class S> inline Vector4D<S> operator-(const Vector4D<S> &v1, const Vector4D<S> &v2) +{ + return Vector4D<S>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.t - v2.t); +} +//! Subtraction operator +template<class S, class S2> inline Vector4D<S> operator-(const Vector4D<S> &v, S2 s) +{ + return Vector4D<S>(v.x - s, v.y - s, v.z - s, v.t - s); +} +//! Subtraction operator +template<class S, class S2> inline Vector4D<S> operator-(S2 s, const Vector4D<S> &v) +{ + return Vector4D<S>(s - v.x, s - v.y, s - v.z, s - v.t); +} + +//! Multiplication operator +template<class S> inline Vector4D<S> operator*(const Vector4D<S> &v1, const Vector4D<S> &v2) +{ + return Vector4D<S>(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z, v1.t * v2.t); +} +//! Multiplication operator +template<class S, class S2> inline Vector4D<S> operator*(const Vector4D<S> &v, S2 s) +{ + return Vector4D<S>(v.x * s, v.y * s, v.z * s, v.t * s); +} +//! Multiplication operator +template<class S, class S2> inline Vector4D<S> operator*(S2 s, const Vector4D<S> &v) +{ + return Vector4D<S>(s * v.x, s * v.y, s * v.z, s * v.t); +} + +//! Division operator +template<class S> inline Vector4D<S> operator/(const Vector4D<S> &v1, const Vector4D<S> &v2) +{ + return Vector4D<S>(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z, v1.t / v2.t); +} +//! Division operator +template<class S, class S2> inline Vector4D<S> operator/(const Vector4D<S> &v, S2 s) +{ + return Vector4D<S>(v.x / s, v.y / s, v.z / s, v.t / s); +} +//! Division operator +template<class S, class S2> inline Vector4D<S> operator/(S2 s, const Vector4D<S> &v) +{ + return Vector4D<S>(s / v.x, s / v.y, s / v.z, s / v.t); +} + +//! Comparison operator +template<class S> inline bool operator==(const Vector4D<S> &s1, const Vector4D<S> &s2) +{ + return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z && s1.t == s2.t; +} + +//! Comparison operator +template<class S> inline bool operator!=(const Vector4D<S> &s1, const Vector4D<S> &s2) +{ + return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z || s1.t != s2.t; +} + +//************************************************************************ +// External functions +//************************************************************************ + +//! Dot product +template<class S> inline S dot(const Vector4D<S> &t, const Vector4D<S> &v) +{ + return t.x * v.x + t.y * v.y + t.z * v.z + t.t * v.t; +} + +//! Cross product +/*template<class S> +inline Vector4D<S> cross ( const Vector4D<S> &t, const Vector4D<S> &v ) { + NYI Vector4D<S> cp ( + ( ( t.y*v.z ) - ( t.z*v.y ) ), + ( ( t.z*v.x ) - ( t.x*v.z ) ), + ( ( t.x*v.y ) - ( t.y*v.x ) ) ); + return cp; +}*/ + +//! Compute the magnitude (length) of the vector +template<class S> inline S norm(const Vector4D<S> &v) +{ + S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; + return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l); +} + +//! Compute squared magnitude +template<class S> inline S normSquare(const Vector4D<S> &v) +{ + return v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; +} + +//! Returns a normalized vector +template<class S> inline Vector4D<S> getNormalized(const Vector4D<S> &v) +{ + S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; + if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) + return v; /* normalized "enough"... */ + else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { + S fac = 1. / sqrt(l); + return Vector4D<S>(v.x * fac, v.y * fac, v.z * fac, v.t * fac); + } + else + return Vector4D<S>((S)0); +} + +//! Compute the norm of the vector and normalize it. +/*! \return The value of the norm */ +template<class S> inline S normalize(Vector4D<S> &v) +{ + S norm; + S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t; + if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) { + norm = 1.; + } + else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { + norm = sqrt(l); + v *= 1. / norm; + } + else { + v = Vector4D<S>::Zero; + norm = 0.; + } + return (S)norm; +} + +//! Outputs the object in human readable form as string +template<class S> std::string Vector4D<S>::toString() const +{ + char buf[256]; + snprintf(buf, + 256, + "[%+4.6f,%+4.6f,%+4.6f,%+4.6f]", + (double)(*this)[0], + (double)(*this)[1], + (double)(*this)[2], + (double)(*this)[3]); + // for debugging, optionally increase precision: + // snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( + // *this ) [1], ( double ) ( *this ) [2], ( double ) ( *this ) [3] ); + return std::string(buf); +} + +template<> std::string Vector4D<int>::toString() const; + +//! Outputs the object in human readable form to stream +template<class S> std::ostream &operator<<(std::ostream &os, const Vector4D<S> &i) +{ + os << i.toString(); + return os; +} + +//! Reads the contents of the object from a stream +template<class S> std::istream &operator>>(std::istream &is, Vector4D<S> &i) +{ + char c; + char dummy[4]; + is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> dummy >> i[3] >> c; + return is; +} + +/**************************************************************************/ +// Define default vector alias +/**************************************************************************/ + +//! 3D vector class of type Real (typically float) +typedef Vector4D<Real> Vec4; + +//! 3D vector class of type int +typedef Vector4D<int> Vec4i; + +//! convert to Real Vector +template<class T> inline Vec4 toVec4(T v) +{ + return Vec4(v[0], v[1], v[2], v[3]); +} +template<class T> inline Vec4i toVec4i(T v) +{ + return Vec4i(v[0], v[1], v[2], v[3]); +} + +/**************************************************************************/ +// Specializations for common math functions +/**************************************************************************/ + +template<> inline Vec4 clamp<Vec4>(const Vec4 &a, const Vec4 &b, const Vec4 &c) +{ + return Vec4( + clamp(a.x, b.x, c.x), clamp(a.y, b.y, c.y), clamp(a.z, b.z, c.z), clamp(a.t, b.t, c.t)); +} +template<> inline Vec4 safeDivide<Vec4>(const Vec4 &a, const Vec4 &b) +{ + return Vec4( + safeDivide(a.x, b.x), safeDivide(a.y, b.y), safeDivide(a.z, b.z), safeDivide(a.t, b.t)); +} +template<> inline Vec4 nmod<Vec4>(const Vec4 &a, const Vec4 &b) +{ + return Vec4(nmod(a.x, b.x), nmod(a.y, b.y), nmod(a.z, b.z), nmod(a.t, b.t)); +} + +/**************************************************************************/ +// 4d interpolation (note only 4d here, 2d/3d interpolations are in interpol.h) +/**************************************************************************/ + +#define BUILD_INDEX_4D \ + Real px = pos.x - 0.5f, py = pos.y - 0.5f, pz = pos.z - 0.5f, pt = pos.t - 0.5f; \ + int xi = (int)px; \ + int yi = (int)py; \ + int zi = (int)pz; \ + int ti = (int)pt; \ + Real s1 = px - (Real)xi, s0 = 1. - s1; \ + Real t1 = py - (Real)yi, t0 = 1. - t1; \ + Real f1 = pz - (Real)zi, f0 = 1. - f1; \ + Real g1 = pt - (Real)ti, g0 = 1. - g1; \ + /* clamp to border */ \ + if (px < 0.) { \ + xi = 0; \ + s0 = 1.0; \ + s1 = 0.0; \ + } \ + if (py < 0.) { \ + yi = 0; \ + t0 = 1.0; \ + t1 = 0.0; \ + } \ + if (pz < 0.) { \ + zi = 0; \ + f0 = 1.0; \ + f1 = 0.0; \ + } \ + if (pt < 0.) { \ + ti = 0; \ + g0 = 1.0; \ + g1 = 0.0; \ + } \ + if (xi >= size.x - 1) { \ + xi = size.x - 2; \ + s0 = 0.0; \ + s1 = 1.0; \ + } \ + if (yi >= size.y - 1) { \ + yi = size.y - 2; \ + t0 = 0.0; \ + t1 = 1.0; \ + } \ + if (zi >= size.z - 1) { \ + zi = size.z - 2; \ + f0 = 0.0; \ + f1 = 1.0; \ + } \ + if (ti >= size.t - 1) { \ + ti = size.t - 2; \ + g0 = 0.0; \ + g1 = 1.0; \ + } \ + const int sX = 1; \ + const int sY = size.x; + +static inline void checkIndexInterpol4d(const Vec4i &size, int idx) +{ + if (idx < 0 || idx > size.x * size.y * size.z * size.t) { + std::ostringstream s; + s << "Grid interpol4d dim " << size << " : index " << idx << " out of bound "; + errMsg(s.str()); + } +} + +template<class T> +inline T interpol4d( + const T *data, const Vec4i &size, const IndexInt sZ, const IndexInt sT, const Vec4 &pos) +{ + BUILD_INDEX_4D + IndexInt idx = (IndexInt)xi + sY * (IndexInt)yi + sZ * (IndexInt)zi + sT * (IndexInt)ti; + DEBUG_ONLY(checkIndexInterpol4d(size, idx)); + DEBUG_ONLY(checkIndexInterpol4d(size, idx + sX + sY + sZ + sT)); + + return (((data[idx] * t0 + data[idx + sY] * t1) * s0 + + (data[idx + sX] * t0 + data[idx + sX + sY] * t1) * s1) * + f0 + + ((data[idx + sZ] * t0 + data[idx + sY + sZ] * t1) * s0 + + (data[idx + sX + sZ] * t0 + data[idx + sX + sY + sZ] * t1) * s1) * + f1) * + g0 + + (((data[idx + sT] * t0 + data[idx + sT + sY] * t1) * s0 + + (data[idx + sT + sX] * t0 + data[idx + sT + sX + sY] * t1) * s1) * + f0 + + ((data[idx + sT + sZ] * t0 + data[idx + sT + sY + sZ] * t1) * s0 + + (data[idx + sT + sX + sZ] * t0 + data[idx + sT + sX + sY + sZ] * t1) * s1) * + f1) * + g1; +} + +}; // namespace Manta + +#endif diff --git a/extern/mantaflow/helper/util/vectorbase.cpp b/extern/mantaflow/helper/util/vectorbase.cpp new file mode 100644 index 00000000000..413ae086d1c --- /dev/null +++ b/extern/mantaflow/helper/util/vectorbase.cpp @@ -0,0 +1,49 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Basic vector class + * + ******************************************************************************/ + +#include "vectorbase.h" + +using namespace std; + +namespace Manta { + +template<> const Vector3D<int> Vector3D<int>::Zero(0, 0, 0); +template<> const Vector3D<float> Vector3D<float>::Zero(0.f, 0.f, 0.f); +template<> const Vector3D<double> Vector3D<double>::Zero(0., 0., 0.); +template<> +const Vector3D<float> Vector3D<float>::Invalid(numeric_limits<float>::quiet_NaN(), + numeric_limits<float>::quiet_NaN(), + numeric_limits<float>::quiet_NaN()); +template<> +const Vector3D<double> Vector3D<double>::Invalid(numeric_limits<double>::quiet_NaN(), + numeric_limits<double>::quiet_NaN(), + numeric_limits<double>::quiet_NaN()); + +template<> bool Vector3D<float>::isValid() const +{ + return !c_isnan(x) && !c_isnan(y) && !c_isnan(z); +} +template<> bool Vector3D<double>::isValid() const +{ + return !c_isnan(x) && !c_isnan(y) && !c_isnan(z); +} + +//! Specialization for readable ints +template<> std::string Vector3D<int>::toString() const +{ + char buf[256]; + snprintf(buf, 256, "[%d,%d,%d]", (*this)[0], (*this)[1], (*this)[2]); + return std::string(buf); +} + +} // namespace Manta diff --git a/extern/mantaflow/helper/util/vectorbase.h b/extern/mantaflow/helper/util/vectorbase.h new file mode 100644 index 00000000000..a3135431eb3 --- /dev/null +++ b/extern/mantaflow/helper/util/vectorbase.h @@ -0,0 +1,679 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011-2016 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Basic vector class + * + ******************************************************************************/ + +#ifndef _VECTORBASE_H +#define _VECTORBASE_H + +// get rid of windos min/max defines +#if defined(WIN32) || defined(_WIN32) +# define NOMINMAX +#endif + +#include <stdio.h> +#include <string> +#include <limits> +#include <iostream> +#include "general.h" + +// if min/max are still around... +#if defined(WIN32) || defined(_WIN32) +# undef min +# undef max +#endif + +// redefine usage of some windows functions +#if defined(WIN32) || defined(_WIN32) +# ifndef snprintf +# define snprintf _snprintf +# endif +#endif + +// use which fp-precision? 1=float, 2=double +#ifndef FLOATINGPOINT_PRECISION +# define FLOATINGPOINT_PRECISION 1 +#endif + +// VECTOR_EPSILON is the minimal vector length +// In order to be able to discriminate floating point values near zero, and +// to be sure not to fail a comparison because of roundoff errors, use this +// value as a threshold. +#if FLOATINGPOINT_PRECISION == 1 +typedef float Real; +# define VECTOR_EPSILON (1e-6f) +#else +typedef double Real; +# define VECTOR_EPSILON (1e-10) +#endif + +#ifndef M_PI +# define M_PI 3.1415926536 +#endif +#ifndef M_E +# define M_E 2.7182818284 +#endif + +namespace Manta { + +//! Basic inlined vector class +template<class S> class Vector3D { + public: + //! Constructor + inline Vector3D() : x(0), y(0), z(0) + { + } + + //! Copy-Constructor + inline Vector3D(const Vector3D<S> &v) : x(v.x), y(v.y), z(v.z) + { + } + + //! Copy-Constructor + inline Vector3D(const int *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) + { + } + + //! Copy-Constructor + inline Vector3D(const float *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) + { + } + + //! Copy-Constructor + inline Vector3D(const double *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) + { + } + + //! Construct a vector from one S + inline Vector3D(S v) : x(v), y(v), z(v) + { + } + + //! Construct a vector from three Ss + inline Vector3D(S vx, S vy, S vz) : x(vx), y(vy), z(vz) + { + } + + // Operators + + //! Assignment operator + inline const Vector3D<S> &operator=(const Vector3D<S> &v) + { + x = v.x; + y = v.y; + z = v.z; + return *this; + } + //! Assignment operator + inline const Vector3D<S> &operator=(S s) + { + x = y = z = s; + return *this; + } + //! Assign and add operator + inline const Vector3D<S> &operator+=(const Vector3D<S> &v) + { + x += v.x; + y += v.y; + z += v.z; + return *this; + } + //! Assign and add operator + inline const Vector3D<S> &operator+=(S s) + { + x += s; + y += s; + z += s; + return *this; + } + //! Assign and sub operator + inline const Vector3D<S> &operator-=(const Vector3D<S> &v) + { + x -= v.x; + y -= v.y; + z -= v.z; + return *this; + } + //! Assign and sub operator + inline const Vector3D<S> &operator-=(S s) + { + x -= s; + y -= s; + z -= s; + return *this; + } + //! Assign and mult operator + inline const Vector3D<S> &operator*=(const Vector3D<S> &v) + { + x *= v.x; + y *= v.y; + z *= v.z; + return *this; + } + //! Assign and mult operator + inline const Vector3D<S> &operator*=(S s) + { + x *= s; + y *= s; + z *= s; + return *this; + } + //! Assign and div operator + inline const Vector3D<S> &operator/=(const Vector3D<S> &v) + { + x /= v.x; + y /= v.y; + z /= v.z; + return *this; + } + //! Assign and div operator + inline const Vector3D<S> &operator/=(S s) + { + x /= s; + y /= s; + z /= s; + return *this; + } + //! Negation operator + inline Vector3D<S> operator-() const + { + return Vector3D<S>(-x, -y, -z); + } + + //! Get smallest component + inline S min() const + { + return (x < y) ? ((x < z) ? x : z) : ((y < z) ? y : z); + } + //! Get biggest component + inline S max() const + { + return (x > y) ? ((x > z) ? x : z) : ((y > z) ? y : z); + } + + //! Test if all components are zero + inline bool empty() + { + return x == 0 && y == 0 && z == 0; + } + + //! access operator + inline S &operator[](unsigned int i) + { + return value[i]; + } + //! constant access operator + inline const S &operator[](unsigned int i) const + { + return value[i]; + } + + //! debug output vector to a string + std::string toString() const; + + //! test if nans are present + bool isValid() const; + + //! actual values + union { + S value[3]; + struct { + S x; + S y; + S z; + }; + struct { + S X; + S Y; + S Z; + }; + }; + + //! zero element + static const Vector3D<S> Zero, Invalid; + + //! For compatibility with 4d vectors (discards 4th comp) + inline Vector3D(S vx, S vy, S vz, S vDummy) : x(vx), y(vy), z(vz) + { + } + + protected: +}; + +//! helper to check whether float/double value is non-zero +inline bool notZero(Real f) +{ + if (std::abs(f) > VECTOR_EPSILON) + return true; + return false; +} + +//************************************************************************ +// Additional operators +//************************************************************************ + +//! Addition operator +template<class S> inline Vector3D<S> operator+(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z); +} +//! Addition operator +template<class S, class S2> inline Vector3D<S> operator+(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x + s, v.y + s, v.z + s); +} +//! Addition operator +template<class S, class S2> inline Vector3D<S> operator+(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(v.x + s, v.y + s, v.z + s); +} + +//! Subtraction operator +template<class S> inline Vector3D<S> operator-(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z); +} +//! Subtraction operator +template<class S, class S2> inline Vector3D<S> operator-(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x - s, v.y - s, v.z - s); +} +//! Subtraction operator +template<class S, class S2> inline Vector3D<S> operator-(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(s - v.x, s - v.y, s - v.z); +} + +//! Multiplication operator +template<class S> inline Vector3D<S> operator*(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z); +} +//! Multiplication operator +template<class S, class S2> inline Vector3D<S> operator*(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x * s, v.y * s, v.z * s); +} +//! Multiplication operator +template<class S, class S2> inline Vector3D<S> operator*(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(s * v.x, s * v.y, s * v.z); +} + +//! Division operator +template<class S> inline Vector3D<S> operator/(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z); +} +//! Division operator +template<class S, class S2> inline Vector3D<S> operator/(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x / s, v.y / s, v.z / s); +} +//! Division operator +template<class S, class S2> inline Vector3D<S> operator/(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(s / v.x, s / v.y, s / v.z); +} + +//! Comparison operator +template<class S> inline bool operator==(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z; +} + +//! Comparison operator +template<class S> inline bool operator!=(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z; +} + +//************************************************************************ +// External functions +//************************************************************************ + +//! Min operator +template<class S> inline Vector3D<S> vmin(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::min(s1.x, s2.x), std::min(s1.y, s2.y), std::min(s1.z, s2.z)); +} + +//! Min operator +template<class S, class S2> inline Vector3D<S> vmin(const Vector3D<S> &s1, S2 s2) +{ + return Vector3D<S>(std::min(s1.x, s2), std::min(s1.y, s2), std::min(s1.z, s2)); +} + +//! Min operator +template<class S1, class S> inline Vector3D<S> vmin(S1 s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::min(s1, s2.x), std::min(s1, s2.y), std::min(s1, s2.z)); +} + +//! Max operator +template<class S> inline Vector3D<S> vmax(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::max(s1.x, s2.x), std::max(s1.y, s2.y), std::max(s1.z, s2.z)); +} + +//! Max operator +template<class S, class S2> inline Vector3D<S> vmax(const Vector3D<S> &s1, S2 s2) +{ + return Vector3D<S>(std::max(s1.x, s2), std::max(s1.y, s2), std::max(s1.z, s2)); +} + +//! Max operator +template<class S1, class S> inline Vector3D<S> vmax(S1 s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::max(s1, s2.x), std::max(s1, s2.y), std::max(s1, s2.z)); +} + +//! Dot product +template<class S> inline S dot(const Vector3D<S> &t, const Vector3D<S> &v) +{ + return t.x * v.x + t.y * v.y + t.z * v.z; +} + +//! Cross product +template<class S> inline Vector3D<S> cross(const Vector3D<S> &t, const Vector3D<S> &v) +{ + Vector3D<S> cp( + ((t.y * v.z) - (t.z * v.y)), ((t.z * v.x) - (t.x * v.z)), ((t.x * v.y) - (t.y * v.x))); + return cp; +} + +//! Project a vector into a plane, defined by its normal +/*! Projects a vector into a plane normal to the given vector, which must + have unit length. Self is modified. + \param v The vector to project + \param n The plane normal + \return The projected vector */ +template<class S> +inline const Vector3D<S> &projectNormalTo(const Vector3D<S> &v, const Vector3D<S> &n) +{ + S sprod = dot(v, n); + return v - n * dot(v, n); +} + +//! Compute the magnitude (length) of the vector +//! (clamps to 0 and 1 with VECTOR_EPSILON) +template<class S> inline S norm(const Vector3D<S> &v) +{ + S l = v.x * v.x + v.y * v.y + v.z * v.z; + if (l <= VECTOR_EPSILON * VECTOR_EPSILON) + return (0.); + return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l); +} + +//! Compute squared magnitude +template<class S> inline S normSquare(const Vector3D<S> &v) +{ + return v.x * v.x + v.y * v.y + v.z * v.z; +} + +//! compatibility, allow use of int, Real and Vec inputs with norm/normSquare +inline Real norm(const Real v) +{ + return fabs(v); +} +inline Real normSquare(const Real v) +{ + return square(v); +} +inline Real norm(const int v) +{ + return abs(v); +} +inline Real normSquare(const int v) +{ + return square(v); +} + +//! Returns a normalized vector +template<class S> inline Vector3D<S> getNormalized(const Vector3D<S> &v) +{ + S l = v.x * v.x + v.y * v.y + v.z * v.z; + if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) + return v; /* normalized "enough"... */ + else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { + S fac = 1. / sqrt(l); + return Vector3D<S>(v.x * fac, v.y * fac, v.z * fac); + } + else + return Vector3D<S>((S)0); +} + +//! Compute the norm of the vector and normalize it. +/*! \return The value of the norm */ +template<class S> inline S normalize(Vector3D<S> &v) +{ + S norm; + S l = v.x * v.x + v.y * v.y + v.z * v.z; + if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) { + norm = 1.; + } + else if (l > VECTOR_EPSILON * VECTOR_EPSILON) { + norm = sqrt(l); + v *= 1. / norm; + } + else { + v = Vector3D<S>::Zero; + norm = 0.; + } + return (S)norm; +} + +//! Obtain an orthogonal vector +/*! Compute a vector that is orthonormal to the given vector. + * Nothing else can be assumed for the direction of the new vector. + * \return The orthonormal vector */ +template<class S> Vector3D<S> getOrthogonalVector(const Vector3D<S> &v) +{ + // Determine the component with max. absolute value + int maxIndex = (fabs(v.x) > fabs(v.y)) ? 0 : 1; + maxIndex = (fabs(v[maxIndex]) > fabs(v.z)) ? maxIndex : 2; + + // Choose another axis than the one with max. component and project + // orthogonal to self + Vector3D<S> o(0.0); + o[(maxIndex + 1) % 3] = 1; + + Vector3D<S> c = cross(v, o); + normalize(c); + return c; +} + +//! Convert vector to polar coordinates +/*! Stable vector to angle conversion + *\param v vector to convert + \param phi unique angle [0,2PI] + \param theta unique angle [0,PI] + */ +template<class S> inline void vecToAngle(const Vector3D<S> &v, S &phi, S &theta) +{ + if (fabs(v.y) < VECTOR_EPSILON) + theta = M_PI / 2; + else if (fabs(v.x) < VECTOR_EPSILON && fabs(v.z) < VECTOR_EPSILON) + theta = (v.y >= 0) ? 0 : M_PI; + else + theta = atan(sqrt(v.x * v.x + v.z * v.z) / v.y); + if (theta < 0) + theta += M_PI; + + if (fabs(v.x) < VECTOR_EPSILON) + phi = M_PI / 2; + else + phi = atan(v.z / v.x); + if (phi < 0) + phi += M_PI; + if (fabs(v.z) < VECTOR_EPSILON) + phi = (v.x >= 0) ? 0 : M_PI; + else if (v.z < 0) + phi += M_PI; +} + +//! Compute vector reflected at a surface +/*! Compute a vector, that is self (as an incoming vector) + * reflected at a surface with a distinct normal vector. + * Note that the normal is reversed, if the scalar product with it is positive. + \param t The incoming vector + \param n The surface normal + \return The new reflected vector + */ +template<class S> inline Vector3D<S> reflectVector(const Vector3D<S> &t, const Vector3D<S> &n) +{ + Vector3D<S> nn = (dot(t, n) > 0.0) ? (n * -1.0) : n; + return (t - nn * (2.0 * dot(nn, t))); +} + +//! Compute vector refracted at a surface +/*! \param t The incoming vector + * \param n The surface normal + * \param nt The "inside" refraction index + * \param nair The "outside" refraction index + * \param refRefl Set to 1 on total reflection + * \return The refracted vector + */ +template<class S> +inline Vector3D<S> refractVector( + const Vector3D<S> &t, const Vector3D<S> &normal, S nt, S nair, int &refRefl) +{ + // from Glassner's book, section 5.2 (Heckberts method) + S eta = nair / nt; + S n = -dot(t, normal); + S tt = 1.0 + eta * eta * (n * n - 1.0); + if (tt < 0.0) { + // we have total reflection! + refRefl = 1; + } + else { + // normal reflection + tt = eta * n - sqrt(tt); + return (t * eta + normal * tt); + } + return t; +} + +//! Outputs the object in human readable form as string +template<class S> std::string Vector3D<S>::toString() const +{ + char buf[256]; + snprintf(buf, + 256, + "[%+4.6f,%+4.6f,%+4.6f]", + (double)(*this)[0], + (double)(*this)[1], + (double)(*this)[2]); + // for debugging, optionally increase precision: + // snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( *this ) + // [1], ( double ) ( *this ) [2] ); + return std::string(buf); +} + +template<> std::string Vector3D<int>::toString() const; + +//! Outputs the object in human readable form to stream +/*! Output format [x,y,z] */ +template<class S> std::ostream &operator<<(std::ostream &os, const Vector3D<S> &i) +{ + os << i.toString(); + return os; +} + +//! Reads the contents of the object from a stream +/*! Input format [x,y,z] */ +template<class S> std::istream &operator>>(std::istream &is, Vector3D<S> &i) +{ + char c; + char dummy[3]; + is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> c; + return is; +} + +/**************************************************************************/ +// Define default vector alias +/**************************************************************************/ + +//! 3D vector class of type Real (typically float) +typedef Vector3D<Real> Vec3; + +//! 3D vector class of type int +typedef Vector3D<int> Vec3i; + +//! convert to Real Vector +template<class T> inline Vec3 toVec3(T v) +{ + return Vec3(v[0], v[1], v[2]); +} + +//! convert to int Vector +template<class T> inline Vec3i toVec3i(T v) +{ + return Vec3i((int)v[0], (int)v[1], (int)v[2]); +} + +//! convert to int Vector +template<class T> inline Vec3i toVec3i(T v0, T v1, T v2) +{ + return Vec3i((int)v0, (int)v1, (int)v2); +} + +//! round, and convert to int Vector +template<class T> inline Vec3i toVec3iRound(T v) +{ + return Vec3i((int)round(v[0]), (int)round(v[1]), (int)round(v[2])); +} + +//! convert to int Vector if values are close enough to an int +template<class T> inline Vec3i toVec3iChecked(T v) +{ + Vec3i ret; + for (size_t i = 0; i < 3; i++) { + Real a = v[i]; + if (fabs(a - floor(a + 0.5)) > 1e-5) + errMsg("argument is not an int, cannot convert"); + ret[i] = (int)(a + 0.5); + } + return ret; +} + +//! convert to double Vector +template<class T> inline Vector3D<double> toVec3d(T v) +{ + return Vector3D<double>(v[0], v[1], v[2]); +} + +//! convert to float Vector +template<class T> inline Vector3D<float> toVec3f(T v) +{ + return Vector3D<float>(v[0], v[1], v[2]); +} + +/**************************************************************************/ +// Specializations for common math functions +/**************************************************************************/ + +template<> inline Vec3 clamp<Vec3>(const Vec3 &a, const Vec3 &b, const Vec3 &c) +{ + return Vec3(clamp(a.x, b.x, c.x), clamp(a.y, b.y, c.y), clamp(a.z, b.z, c.z)); +} +template<> inline Vec3 safeDivide<Vec3>(const Vec3 &a, const Vec3 &b) +{ + return Vec3(safeDivide(a.x, b.x), safeDivide(a.y, b.y), safeDivide(a.z, b.z)); +} +template<> inline Vec3 nmod<Vec3>(const Vec3 &a, const Vec3 &b) +{ + return Vec3(nmod(a.x, b.x), nmod(a.y, b.y), nmod(a.z, b.z)); +} + +}; // namespace Manta + +#endif |