mrk it up!
代码语言：javascript复制#include "SimulatorBase.h"
#include "SPlisHSPlasH/Utilities/SceneLoader.h"
#include "Utilities/FileSystem.h"
#include "SPlisHSPlasH/TimeManager.h"
#include "Utilities/PartioReaderWriter.h"
#include "SPlisHSPlasH/Emitter.h"
#include "SPlisHSPlasH/EmitterSystem.h"
#include "SPlisHSPlasH/Simulation.h"
#include "SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.h"
#include "NumericParameter.h"
#include "Utilities/Logger.h"
#include "Utilities/Timing.h"
#include "Utilities/Counting.h"
#include "Utilities/Version.h"
#include "Utilities/SystemInfo.h"
#include "extern/partio/src/lib/Partio.h"
#include "SPlisHSPlasH/Utilities/GaussQuadrature.h"
#include "SPlisHSPlasH/Utilities/SimpleQuadrature.h"
#include "extern/cxxopts/cxxopts.hpp"
#include "Simulator/GUI/Simulator_GUI_Base.h"
#include "SPlisHSPlasH/Utilities/VolumeSampling.h"
#include "Utilities/OBJLoader.h"
#include "Utilities/BinaryFileReaderWriter.h"
#include "PositionBasedDynamicsWrapper/PBDBoundarySimulator.h"
#include "StaticBoundarySimulator.h"
INIT_LOGGING
INIT_TIMING
INIT_COUNTING
using namespace SPH;
using namespace std;
using namespace GenParam;
using namespace Utilities;
int SimulatorBase::PAUSE = -1;
int SimulatorBase::PAUSE_AT = -1;
int SimulatorBase::STOP_AT = -1;
int SimulatorBase::NUM_STEPS_PER_RENDER = -1;
int SimulatorBase::PARTIO_EXPORT = -1;
int SimulatorBase::VTK_EXPORT = -1;
int SimulatorBase::RB_VTK_EXPORT = -1;
int SimulatorBase::RB_EXPORT = -1;
int SimulatorBase::DATA_EXPORT_FPS = -1;
int SimulatorBase::PARTICLE_EXPORT_ATTRIBUTES = -1;
int SimulatorBase::STATE_EXPORT = -1;
int SimulatorBase::STATE_EXPORT_FPS = -1;
int SimulatorBase::RENDER_WALLS = -1;
int SimulatorBase::ENUM_WALLS_NONE = -1;
int SimulatorBase::ENUM_WALLS_PARTICLES_ALL = -1;
int SimulatorBase::ENUM_WALLS_PARTICLES_NO_WALLS = -1;
int SimulatorBase::ENUM_WALLS_GEOMETRY_ALL = -1;
int SimulatorBase::ENUM_WALLS_GEOMETRY_NO_WALLS = -1;
SimulatorBase::SimulatorBase()
{
Utilities::logger.addSink(unique_ptr<Utilities::ConsoleSink>(new Utilities::ConsoleSink(Utilities::LogLevel::INFO)));
m_boundarySimulator = nullptr;
m_gui = nullptr;
m_isStaticScene = true;
m_numberOfStepsPerRenderUpdate = 4;
m_sceneFile = "";
m_renderWalls = 4;
m_doPause = true;
m_pauseAt = -1.0;
m_stopAt = -1.0;
m_useParticleCaching = true;
m_useGUI = true;
m_enablePartioExport = false;
m_enableVTKExport = false;
m_enableRigidBodyVTKExport = false;
m_enableRigidBodyExport = false;
m_enableStateExport = false;
m_framesPerSecond = 25;
m_framesPerSecondState = 1;
m_nextFrameTime = 0.0;
m_nextFrameTimeState = 0.0;
m_frameCounter = 1;
m_isFirstFrame = true;
m_isFirstFrameVTK = true;
m_firstState = true;
m_colorField.resize(1, "velocity");
m_colorMapType.resize(1, 0);
m_renderMinValue.resize(1, 0.0);
m_renderMaxValue.resize(1, 5.0);
m_particleAttributes = "velocity";
m_timeStepCB = nullptr;
#ifdef DL_OUTPUT
m_nextTiming = 1.0;
#endif
}
SimulatorBase::~SimulatorBase()
{
Utilities::Timing::printAverageTimes();
Utilities::Timing::printTimeSums();
Utilities::Counting::printAverageCounts();
Utilities::Counting::printCounterSums();
delete m_boundarySimulator;
}
void SimulatorBase::initParameters()
{
ParameterObject::initParameters();
PAUSE = createBoolParameter("pause", "Pause", &m_doPause);
setGroup(PAUSE, "General");
setDescription(PAUSE, "Pause simulation.");
setHotKey(PAUSE, "space");
PAUSE_AT = createNumericParameter("pauseAt", "Pause simulation at", &m_pauseAt);
setGroup(PAUSE_AT, "General");
setDescription(PAUSE_AT, "Pause simulation at the given time. When the value is negative, the simulation is not paused.");
STOP_AT = createNumericParameter("stopAt", "Stop simulation at", &m_stopAt);
setGroup(STOP_AT, "General");
setDescription(STOP_AT, "Stop simulation at the given time. When the value is negative, the simulation is not stopped.");
NUM_STEPS_PER_RENDER = createNumericParameter("numberOfStepsPerRenderUpdate", "# time steps / update", &m_numberOfStepsPerRenderUpdate);
setGroup(NUM_STEPS_PER_RENDER, "Visualization");
setDescription(NUM_STEPS_PER_RENDER, "Number of simulation steps per rendered frame.");
static_cast<NumericParameter<unsigned int>*>(getParameter(NUM_STEPS_PER_RENDER))->setMinValue(1);
RENDER_WALLS = createEnumParameter("renderWalls", "Render walls", &m_renderWalls);
setGroup(RENDER_WALLS, "Visualization");
setDescription(RENDER_WALLS, "Make walls visible/invisible.");
EnumParameter *enumParam = static_cast<EnumParameter*>(getParameter(RENDER_WALLS));
enumParam->addEnumValue("None", ENUM_WALLS_NONE);
enumParam->addEnumValue("Particles (all)", ENUM_WALLS_PARTICLES_ALL);
enumParam->addEnumValue("Particles (no walls)", ENUM_WALLS_PARTICLES_NO_WALLS);
enumParam->addEnumValue("Geometry (all)", ENUM_WALLS_GEOMETRY_ALL);
enumParam->addEnumValue("Geometry (no walls)", ENUM_WALLS_GEOMETRY_NO_WALLS);
PARTIO_EXPORT = createBoolParameter("enablePartioExport", "Partio export", &m_enablePartioExport);
setGroup(PARTIO_EXPORT, "Export");
setDescription(PARTIO_EXPORT, "Enable/disable partio export.");
RB_EXPORT = createBoolParameter("enableRigidBodyExport", "Rigid body export", &m_enableRigidBodyExport);
setGroup(RB_EXPORT, "Export");
setDescription(RB_EXPORT, "Enable/disable rigid body export.");
VTK_EXPORT = createBoolParameter("enableVTKExport", "VTK export", &m_enableVTKExport);
setGroup(VTK_EXPORT, "Export");
setDescription(VTK_EXPORT, "Enable/disable VTK export.");
RB_VTK_EXPORT = createBoolParameter("enableRigidBodyVTKExport", "Rigid body VTK export", &m_enableRigidBodyVTKExport);
setGroup(RB_VTK_EXPORT, "Export");
setDescription(RB_VTK_EXPORT, "Enable/disable rigid body VTK export.");
DATA_EXPORT_FPS = createNumericParameter("dataExportFPS", "Export FPS", &m_framesPerSecond);
setGroup(DATA_EXPORT_FPS, "Export");
setDescription(DATA_EXPORT_FPS, "Frame rate of partio, vtk and rigid body export.");
STATE_EXPORT = createBoolParameter("enableStateExport", "Simulation state export", &m_enableStateExport);
setGroup(STATE_EXPORT, "Export");
setDescription(STATE_EXPORT, "Enable/disable export of complete simulation state.");
STATE_EXPORT_FPS = createNumericParameter("stateExportFPS", "State export FPS", &m_framesPerSecondState);
setGroup(STATE_EXPORT_FPS, "Export");
setDescription(STATE_EXPORT_FPS, "Frame rate of simulation state export.");
PARTICLE_EXPORT_ATTRIBUTES = createStringParameter("particleAttributes", "Export attributes", &m_particleAttributes);
getParameter(PARTICLE_EXPORT_ATTRIBUTES)->setReadOnly(true);
setGroup(PARTICLE_EXPORT_ATTRIBUTES, "Export");
setDescription(PARTICLE_EXPORT_ATTRIBUTES, "Attributes that are exported in the partio files (except id and position).");
}
void SimulatorBase::run()
{
initSimulation();
runSimulation();
cleanup();
}
void SimulatorBase::init(std::vector<std::string> argv, const std::string &windowName){
m_argc = static_cast<int>(argv.size());
m_argv_vec.clear();
m_argv_vec.reserve(argv.size());
for (auto & a : argv)
m_argv_vec.push_back(&a[0]);
m_argv = m_argv_vec.data();
init(m_argc, m_argv, windowName);
}
void SimulatorBase::init(int argc, char **argv, const std::string &windowName)
{
m_argc = argc;
m_argv = argv;
m_windowName = windowName;
initParameters();
m_exePath = FileSystem::getProgramPath();
setUseParticleCaching(true);
try
{
cxxopts::Options options(argv[0], "SPlisHSPlasH - An open-source library for the physically-based simulation of fluids.");
options
.positional_help("[scene file]")
.show_positional_help();
options.add_options()
("h,help", "Print help")
("v,version", "Print version")
("no-cache", "Disable caching of boundary samples/maps.")
("state-file", "State file (state_<time>.bin) that should be loaded.", cxxopts::value<std::string>())
("output-dir", "Output directory for log file and partio files.", cxxopts::value<std::string>())
("no-initial-pause", "Disable initial pause when starting the simulation.")
("no-gui", "Disable GUI.")
("stopAt", "Sets or overwrites the stopAt parameter of the scene.", cxxopts::value<Real>())
("param", "Sets or overwrites a parameter of the scene.nn"
"- Setting a fluid parameter:nt<fluid-id>:<parameter-name>:<value>n"
"- Example: --param Fluid:viscosity:0.01nn"
"- Setting a configuration parameter:nt<parameter-name>:<value>n"
"- Example: --param cflMethod:1n"
, cxxopts::value<std::string>())
;
options.add_options("invisible")
("scene-file", "Scene file", cxxopts::value<std::string>());
options.parse_positional("scene-file");
auto result = options.parse(argc, argv);
if (result.count("help"))
{
std::cout << options.help({ "" }) << std::endl;
exit(0);
}
if (result.count("version"))
{
std::cout << SPLISHSPLASH_VERSION << std::endl;
exit(0);
}
if (result.count("no-cache"))
{
setUseParticleCaching(false);
}
if (result.count("no-gui"))
{
setUseGUI(false);
}
if (result.count("stopAt"))
{
m_stopAt = result["stopAt"].as<Real>();
}
if (result.count("param"))
{
const string paramStr = result["param"].as<std::string>();
Utilities::StringTools::tokenize(paramStr, m_paramTokens, ":");
// add third element to get a unified method for all parameters
if (m_paramTokens.size() == 2)
m_paramTokens.insert(m_paramTokens.begin(), "config");
if (m_paramTokens.size() != 3)
{
LOG_ERR << "--param has wrong syntax!";
LOG_ERR << "Example 1: --param Fluid:viscosity:0.01";
LOG_ERR << "Example 2: --param cflMethod:1";
exit(1);
}
}
m_sceneFile = "";
if (result.count("scene-file"))
{
m_sceneFile = result["scene-file"].as<std::string>();
if (FileSystem::isRelativePath(m_sceneFile))
m_sceneFile = FileSystem::normalizePath(m_exePath "/" m_sceneFile);
}
#ifdef WIN32
else
{
std::string scenePath = FileSystem::normalizePath(m_exePath "/../data/Scenes/");
if (!FileSystem::isDirectory(scenePath))
scenePath = m_exePath;
std::replace(scenePath.begin(), scenePath.end(), '/', '\');
m_sceneFile = FileSystem::fileDialog(0, scenePath.c_str(), "*.json");
if (m_sceneFile == "")
exit(0);
}
#endif
m_outputPath = FileSystem::normalizePath(getExePath() "/output/" FileSystem::getFileName(m_sceneFile));
if (result.count("output-dir"))
{
m_outputPath = result["output-dir"].as<std::string>();
}
if (result.count("no-initial-pause"))
{
m_doPause = false;
}
setStateFile("");
if (result.count("state-file"))
{
setStateFile(result["state-file"].as<std::string>());
if (FileSystem::isRelativePath(getStateFile()))
setStateFile(FileSystem::normalizePath(m_exePath "/" getStateFile()));
}
}
catch (const cxxopts::OptionException& e)
{
std::cout << "error parsing options: " << e.what() << std::endl;
exit(1);
}
std::string logPath = FileSystem::normalizePath(m_outputPath "/log");
FileSystem::makeDirs(logPath);
Utilities::logger.addSink(unique_ptr<Utilities::FileSink>(new Utilities::FileSink(Utilities::LogLevel::DEBUG, logPath "/SPH_log.txt")));
LOG_INFO << "SPlisHSPlasH version: " << SPLISHSPLASH_VERSION;
LOG_DEBUG << "Git refspec: " << GIT_REFSPEC;
LOG_DEBUG << "Git SHA1: " << GIT_SHA1;
LOG_DEBUG << "Git status: " << GIT_LOCAL_STATUS;
LOG_DEBUG << "Host name: " << SystemInfo::getHostName();
if (!getUseParticleCaching())
LOG_INFO << "Boundary cache disabled.";
LOG_INFO << "Output directory: " << m_outputPath;
//////////////////////////////////////////////////////////////////////////
// read scene
//////////////////////////////////////////////////////////////////////////
m_sceneLoader = std::unique_ptr<SceneLoader>(new SceneLoader());
if (m_sceneFile != "")
m_sceneLoader->readScene(m_sceneFile.c_str(), m_scene);
else
return;
//////////////////////////////////////////////////////////////////////////
// init boundary simulation
//////////////////////////////////////////////////////////////////////////
m_isStaticScene = true;
for (unsigned int i = 0; i < m_scene.boundaryModels.size(); i )
{
if (m_scene.boundaryModels[i]->dynamic)
{
m_isStaticScene = false;
break;
}
}
if (m_isStaticScene)
{
LOG_INFO << "Initialize static boundary simulation";
m_boundarySimulator = new StaticBoundarySimulator(this);
}
else
{
LOG_INFO << "Initialize dynamic boundary simulation";
m_boundarySimulator = new PBDBoundarySimulator(this);
}
}
void SimulatorBase::initSimulation()
{
#ifdef DL_OUTPUT
// copy scene files in output so that the simulation can be reproduced
std::string sceneFilePath = FileSystem::normalizePath(m_outputPath "/scene");
FileSystem::makeDirs(sceneFilePath);
FileSystem::copyFile(m_sceneFile, sceneFilePath "/" FileSystem::getFileNameWithExt(m_sceneFile));
std::string modelsFilePath = FileSystem::normalizePath(m_outputPath "/models");
FileSystem::makeDirs(modelsFilePath);
for (unsigned int i = 0; i < m_scene.boundaryModels.size(); i )
{
std::string meshFileName = m_scene.boundaryModels[i]->meshFile;
if (meshFileName != "")
{
if (FileSystem::isRelativePath(meshFileName))
meshFileName = FileSystem::normalizePath(FileSystem::getFilePath(m_sceneFile) "/" meshFileName);
FileSystem::copyFile(meshFileName, modelsFilePath "/" FileSystem::getFileNameWithExt(meshFileName));
}
std::string mapFileName = m_scene.boundaryModels[i]->mapFile;
if (mapFileName != "")
{
if (FileSystem::isRelativePath(mapFileName))
mapFileName = FileSystem::normalizePath(FileSystem::getFilePath(m_sceneFile) "/" mapFileName);
FileSystem::copyFile(m_scene.boundaryModels[i]->mapFile, modelsFilePath "/" FileSystem::getFileNameWithExt(m_scene.boundaryModels[i]->mapFile));
}
}
for (unsigned int i = 0; i < m_scene.fluidModels.size(); i )
{
std::string samplesFileName = m_scene.fluidModels[i]->samplesFile;
if (samplesFileName != "")
{
if (FileSystem::isRelativePath(samplesFileName))
samplesFileName = FileSystem::normalizePath(FileSystem::getFilePath(m_sceneFile) "/" samplesFileName);
FileSystem::copyFile(samplesFileName, modelsFilePath "/" FileSystem::getFileNameWithExt(samplesFileName));
}
}
std::string progFilePath = FileSystem::normalizePath(m_outputPath "/program");
FileSystem::makeDirs(progFilePath);
FileSystem::copyFile(m_argv[0], progFilePath "/" FileSystem::getFileNameWithExt(m_argv[0]));
#endif
Simulation *sim = Simulation::getCurrent();
sim->init(getScene().particleRadius, getScene().sim2D);
buildModel();
if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
{
unsigned int nBoundaryParticles = 0;
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
nBoundaryParticles = static_cast<BoundaryModel_Akinci2012*>(sim->getBoundaryModel(i))->numberOfParticles();
LOG_INFO << "Number of boundary particles: " << nBoundaryParticles;
}
if (m_useGUI)
m_gui->init(m_argc, m_argv, m_windowName.c_str());
m_boundarySimulator->init();
if (m_useGUI)
{
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
const std::string &key = model->getId();
model->setDragMethodChangedCallback([this, model]() { m_gui->initSimulationParameterGUI(); getSceneLoader()->readMaterialParameterObject(model->getId(), (ParameterObject*)model->getDragBase()); });
model->setSurfaceMethodChangedCallback([this, model]() { m_gui->initSimulationParameterGUI(); getSceneLoader()->readMaterialParameterObject(model->getId(), (ParameterObject*)model->getSurfaceTensionBase()); });
model->setViscosityMethodChangedCallback([this, model]() { m_gui->initSimulationParameterGUI(); getSceneLoader()->readMaterialParameterObject(model->getId(), (ParameterObject*)model->getViscosityBase()); });
model->setVorticityMethodChangedCallback([this, model]() { m_gui->initSimulationParameterGUI(); getSceneLoader()->readMaterialParameterObject(model->getId(), (ParameterObject*)model->getVorticityBase()); });
model->setElasticityMethodChangedCallback([this, model]() { reset(); m_gui->initSimulationParameterGUI(); getSceneLoader()->readMaterialParameterObject(model->getId(), (ParameterObject*)model->getElasticityBase()); });
}
m_gui->initSimulationParameterGUI();
Simulation::getCurrent()->setSimulationMethodChangedCallback([this]() { reset(); m_gui->initSimulationParameterGUI(); getSceneLoader()->readParameterObject("Configuration", Simulation::getCurrent()->getTimeStep()); });
}
readParameters();
setCommandLineParameter();
}
void SimulatorBase::runSimulation()
{
m_boundarySimulator->initBoundaryData();
if (getStateFile() != "")
loadState(getStateFile());
if (!m_useGUI)
{
const Real stopAt = getValue<Real>(SimulatorBase::STOP_AT);
if (stopAt < 0.0)
{
LOG_ERR << "StopAt parameter must be set when starting without GUI.";
exit(1);
}
while (true)
{
if (!timeStepNoGUI())
break;
}
}
if (m_useGUI)
m_gui->run();
}
void SimulatorBase::cleanup()
{
for (unsigned int i = 0; i < m_scene.boundaryModels.size(); i )
delete m_scene.boundaryModels[i];
m_scene.boundaryModels.clear();
for (unsigned int i = 0; i < m_scene.fluidModels.size(); i )
delete m_scene.fluidModels[i];
m_scene.fluidModels.clear();
for (unsigned int i = 0; i < m_scene.fluidBlocks.size(); i )
delete m_scene.fluidBlocks[i];
m_scene.fluidBlocks.clear();
for (unsigned int i = 0; i < m_scene.emitters.size(); i )
delete m_scene.emitters[i];
m_scene.emitters.clear();
for (unsigned int i = 0; i < m_scene.animatedFields.size(); i )
delete m_scene.animatedFields[i];
m_scene.animatedFields.clear();
if (m_useGUI)
m_gui->cleanup();
delete Simulation::getCurrent();
}
void SimulatorBase::readParameters()
{
m_sceneLoader->readParameterObject("Configuration", this);
m_sceneLoader->readParameterObject("Configuration", Simulation::getCurrent());
m_sceneLoader->readParameterObject("Configuration", Simulation::getCurrent()->getTimeStep());
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
const std::string &key = model->getId();
m_sceneLoader->readMaterialParameterObject(key, model);
m_sceneLoader->readMaterialParameterObject(key, (ParameterObject*) model->getDragBase());
m_sceneLoader->readMaterialParameterObject(key, (ParameterObject*) model->getSurfaceTensionBase());
m_sceneLoader->readMaterialParameterObject(key, (ParameterObject*) model->getViscosityBase());
m_sceneLoader->readMaterialParameterObject(key, (ParameterObject*) model->getVorticityBase());
m_sceneLoader->readMaterialParameterObject(key, (ParameterObject*) model->getElasticityBase());
for (auto material : m_scene.materials)
{
if (material->id == key)
{
setColorField(i, material->colorField);
setColorMapType(i, material->colorMapType);
setRenderMinValue(i, material->minVal);
setRenderMaxValue(i, material->maxVal);
}
}
}
}
void SimulatorBase::setCommandLineParameter()
{
Simulation *sim = Simulation::getCurrent();
if (m_paramTokens.size() != 3)
return;
setCommandLineParameter((ParameterObject*)sim);
setCommandLineParameter((ParameterObject*)sim->getTimeStep());
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
const std::string &key = model->getId();
if (m_paramTokens[0] == key)
{
setCommandLineParameter((ParameterObject*)model);
setCommandLineParameter((ParameterObject*)model->getDragBase());
setCommandLineParameter((ParameterObject*)model->getSurfaceTensionBase());
setCommandLineParameter((ParameterObject*)model->getViscosityBase());
setCommandLineParameter((ParameterObject*)model->getVorticityBase());
setCommandLineParameter((ParameterObject*)model->getElasticityBase());
}
}
}
void SimulatorBase::setCommandLineParameter(GenParam::ParameterObject *paramObj)
{
if (paramObj == nullptr)
return;
const unsigned int numParams = paramObj->numParameters();
for (unsigned int j = 0; j < numParams; j )
{
ParameterBase *paramBase = paramObj->getParameter(j);
if (m_paramTokens[1] == paramBase->getName())
{
if (paramBase->getType() == RealParameterType)
{
const Real val = stof(m_paramTokens[2]);
static_cast<NumericParameter<Real>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::UINT32)
{
const unsigned int val = stoi(m_paramTokens[2]);
static_cast<NumericParameter<unsigned int>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::UINT16)
{
const unsigned short val = stoi(m_paramTokens[2]);
static_cast<NumericParameter<unsigned short>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::UINT8)
{
const unsigned char val = stoi(m_paramTokens[2]);
static_cast<NumericParameter<unsigned char>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::INT32)
{
const int val = stoi(m_paramTokens[2]);
static_cast<NumericParameter<int>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::INT16)
{
const short val = stoi(m_paramTokens[2]);
static_cast<NumericParameter<short>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::INT8)
{
const char val = stoi(m_paramTokens[2]);
static_cast<NumericParameter<char>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::ENUM)
{
const int val = stoi(m_paramTokens[2]);
static_cast<EnumParameter*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::BOOL)
{
const bool val = stoi(m_paramTokens[2]);
static_cast<BoolParameter*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == RealVectorParameterType)
{
if (static_cast<VectorParameter<Real>*>(paramBase)->getDim() == 3)
{
vector<string> tokens;
Utilities::StringTools::tokenize(m_paramTokens[2], tokens, ",");
if (tokens.size() == 3)
{
Vector3r val(stof(tokens[0]), stof(tokens[1]), stof(tokens[2]));
static_cast<VectorParameter<Real>*>(paramBase)->setValue(val.data());
}
}
}
else if (paramBase->getType() == ParameterBase::STRING)
{
const std::string val = m_paramTokens[2];
static_cast<StringParameter*>(paramBase)->setValue(val);
}
}
}
}
void SimulatorBase::buildModel()
{
TimeManager::getCurrent()->setTimeStepSize(m_scene.timeStepSize);
initFluidData();
createEmitters();
createAnimationFields();
Simulation *sim = Simulation::getCurrent();
if (sim->getTimeStep())
sim->getTimeStep()->resize();
if (!sim->is2DSimulation())
{
sim->setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_PRECOMPUTED_CUBIC);
sim->setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_PRECOMPUTED_CUBIC);
}
else
{
sim->setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_CUBIC_2D);
sim->setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_CUBIC_2D);
}
}
void SimulatorBase::reset()
{
Utilities::Timing::printAverageTimes();
Utilities::Timing::reset();
Utilities::Counting::printAverageCounts();
Utilities::Counting::reset();
Simulation::getCurrent()->reset();
m_boundarySimulator->reset();
if (m_gui)
m_gui->reset();
if (Simulation::getCurrent()->getValue<int>(Simulation::CFL_METHOD) != Simulation::ENUM_CFL_NONE)
TimeManager::getCurrent()->setTimeStepSize(m_scene.timeStepSize);
m_nextFrameTime = 0.0;
m_nextFrameTimeState = 0.0;
m_frameCounter = 1;
m_isFirstFrame = true;
m_isFirstFrameVTK = true;
#ifdef DL_OUTPUT
m_nextTiming = 1.0;
#endif
}
void SimulatorBase::timeStep()
{
const Real stopAt = getValue<Real>(SimulatorBase::STOP_AT);
if (m_gui && (stopAt > 0.0) && (stopAt < TimeManager::getCurrent()->getTime()))
m_gui->stop();
const Real pauseAt = getValue<Real>(SimulatorBase::PAUSE_AT);
if ((pauseAt > 0.0) && (pauseAt < TimeManager::getCurrent()->getTime()))
setValue(SimulatorBase::PAUSE, true);
if (getValue<bool>(SimulatorBase::PAUSE))
return;
// Simulation code
Simulation *sim = Simulation::getCurrent();
const bool sim2D = sim->is2DSimulation();
const unsigned int numSteps = getValue<unsigned int>(SimulatorBase::NUM_STEPS_PER_RENDER);
for (unsigned int i = 0; i < numSteps; i )
{
START_TIMING("SimStep");
Simulation::getCurrent()->getTimeStep()->step();
STOP_TIMING_AVG;
m_boundarySimulator->timeStep();
step();
INCREASE_COUNTER("Time step size", TimeManager::getCurrent()->getTimeStepSize());
// Make sure that particles stay in xy-plane in a 2D simulation
if (sim2D)
{
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
for (unsigned int i = 0; i < model->numActiveParticles(); i )
{
model->getPosition(i)[2] = 0.0;
model->getVelocity(i)[2] = 0.0;
}
}
}
if (m_timeStepCB)
m_timeStepCB();
}
}
bool SimulatorBase::timeStepNoGUI()
{
const Real stopAt = getValue<Real>(SimulatorBase::STOP_AT);
if ((stopAt > 0.0) && (stopAt < TimeManager::getCurrent()->getTime()))
return false;
// Simulation code
Simulation *sim = Simulation::getCurrent();
const bool sim2D = sim->is2DSimulation();
START_TIMING("SimStep");
Simulation::getCurrent()->getTimeStep()->step();
STOP_TIMING_AVG;
m_boundarySimulator->timeStep();
step();
INCREASE_COUNTER("Time step size", TimeManager::getCurrent()->getTimeStepSize());
// Make sure that particles stay in xy-plane in a 2D simulation
if (sim2D)
{
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
for (unsigned int i = 0; i < model->numActiveParticles(); i )
{
model->getPosition(i)[2] = 0.0;
model->getVelocity(i)[2] = 0.0;
}
}
}
if (m_timeStepCB)
m_timeStepCB();
return true;
}
void SimulatorBase::loadObj(const std::string &filename, TriangleMesh &mesh, const Vector3r &scale)
{
std::vector<OBJLoader::Vec3f> x;
std::vector<OBJLoader::Vec3f> normals;
std::vector<MeshFaceIndices> faces;
OBJLoader::Vec3f s = { (float)scale[0], (float)scale[1], (float)scale[2] };
OBJLoader::loadObj(filename, &x, &faces, &normals, nullptr, s);
mesh.release();
const unsigned int nPoints = (unsigned int)x.size();
const unsigned int nFaces = (unsigned int)faces.size();
mesh.initMesh(nPoints, nFaces);
for (unsigned int i = 0; i < nPoints; i )
{
mesh.addVertex(Vector3r(x[i][0], x[i][1], x[i][2]));
}
for (unsigned int i = 0; i < nFaces; i )
{
// Reduce the indices by one
int posIndices[3];
for (int j = 0; j < 3; j )
{
posIndices[j] = faces[i].posIndices[j] - 1;
}
mesh.addFace(&posIndices[0]);
}
LOG_INFO << "Number of triangles: " << nFaces;
LOG_INFO << "Number of vertices: " << nPoints;
}
void SimulatorBase::initFluidData()
{
LOG_INFO << "Initialize fluid particles";
Simulation *sim = Simulation::getCurrent();
//////////////////////////////////////////////////////////////////////////
// Determine number of different fluid IDs
//////////////////////////////////////////////////////////////////////////
std::map<std::string, unsigned int> fluidIDs;
unsigned int index = 0;
for (unsigned int i = 0; i < m_scene.fluidBlocks.size(); i )
{
if (fluidIDs.find(m_scene.fluidBlocks[i]->id) == fluidIDs.end())
fluidIDs[m_scene.fluidBlocks[i]->id] = index ;
}
for (unsigned int i = 0; i < m_scene.fluidModels.size(); i )
{
if (fluidIDs.find(m_scene.fluidModels[i]->id) == fluidIDs.end())
fluidIDs[m_scene.fluidModels[i]->id] = index ;
}
for (unsigned int i = 0; i < m_scene.emitters.size(); i )
{
if (fluidIDs.find(m_scene.emitters[i]->id) == fluidIDs.end())
fluidIDs[m_scene.emitters[i]->id] = index ;
}
const unsigned int numberOfFluidModels = static_cast<unsigned int>(fluidIDs.size());
std::vector<std::vector<Vector3r>> fluidParticles;
std::vector<std::vector<Vector3r>> fluidVelocities;
fluidParticles.resize(numberOfFluidModels);
fluidVelocities.resize(numberOfFluidModels);
createFluidBlocks(fluidIDs, fluidParticles, fluidVelocities);
std::string base_path = FileSystem::getFilePath(m_sceneFile);
const bool useCache = getUseParticleCaching();
std::string scene_path = FileSystem::getFilePath(getSceneFile());
string cachePath = scene_path "/Cache";
unsigned int startIndex = 0;
unsigned int endIndex = 0;
for (unsigned int i = 0; i < m_scene.fluidModels.size(); i )
{
const unsigned int fluidIndex = fluidIDs[m_scene.fluidModels[i]->id];
std::string fileName;
if (FileSystem::isRelativePath(m_scene.fluidModels[i]->samplesFile))
fileName = base_path "/" m_scene.fluidModels[i]->samplesFile;
else
fileName = m_scene.fluidModels[i]->samplesFile;
string ext = FileSystem::getFileExt(fileName);
transform(ext.begin(), ext.end(), ext.begin(), ::toupper);
if (ext == "OBJ")
{
// check if mesh file has changed
std::string md5FileName = FileSystem::normalizePath(cachePath "/" FileSystem::getFileNameWithExt(fileName) "_fluid.md5");
bool md5 = false;
if (useCache)
{
string md5Str = FileSystem::getFileMD5(fileName);
if (FileSystem::fileExists(md5FileName))
md5 = FileSystem::checkMD5(md5Str, md5FileName);
}
// Cache sampling
std::string mesh_base_path = FileSystem::getFilePath(fileName);
std::string mesh_file_name = FileSystem::getFileName(fileName);
std::string mode = to_string(m_scene.fluidModels[i]->mode);
const string scaleStr = real2String(m_scene.fluidModels[i]->scale[0]) "_" real2String(m_scene.fluidModels[i]->scale[1]) "_" real2String(m_scene.fluidModels[i]->scale[2]);
const string resStr = to_string(m_scene.fluidModels[i]->resolutionSDF[0]) "_" to_string(m_scene.fluidModels[i]->resolutionSDF[1]) "_" to_string(m_scene.fluidModels[i]->resolutionSDF[2]);
const string particleFileName = FileSystem::normalizePath(cachePath "/" mesh_file_name "_fluid_" real2String(m_scene.particleRadius) "_m" mode "_s" scaleStr "_r" resStr ".bgeo");
// check MD5 if cache file is available
bool foundCacheFile = false;
if (useCache)
foundCacheFile = FileSystem::fileExists(particleFileName);
if (useCache && foundCacheFile && md5)
{
PartioReaderWriter::readParticles(particleFileName, m_scene.fluidModels[i]->translation, m_scene.fluidModels[i]->rotation, 1.0, fluidParticles[fluidIndex], fluidVelocities[fluidIndex]);
LOG_INFO << "Loaded cached fluid sampling: " << particleFileName;
}
if (!useCache || !foundCacheFile || !md5)
{
LOG_INFO << "Volume sampling of " << fileName;
TriangleMesh mesh;
loadObj(fileName, mesh, m_scene.fluidModels[i]->scale);
bool invert = m_scene.fluidModels[i]->invert;
int mode = m_scene.fluidModels[i]->mode;
std::array<unsigned int, 3> resolutionSDF = m_scene.fluidModels[i]->resolutionSDF;
LOG_INFO << "SDF resolution: " << resolutionSDF[0] << ", " << resolutionSDF[1] << ", " << resolutionSDF[2];
START_TIMING("Volume sampling");
Utilities::VolumeSampling::sampleMesh(mesh.numVertices(), mesh.getVertices().data(), mesh.numFaces(), mesh.getFaces().data(),
m_scene.particleRadius, nullptr, resolutionSDF, invert, mode, fluidParticles[fluidIndex]);
STOP_TIMING_AVG;
fluidVelocities[fluidIndex].resize(fluidParticles[fluidIndex].size(), m_scene.fluidModels[i]->initialVelocity);
// Cache sampling
if (useCache && (FileSystem::makeDir(cachePath) == 0))
{
LOG_INFO << "Save particle sampling: " << particleFileName;
PartioReaderWriter::writeParticles(particleFileName, (unsigned int)fluidParticles[fluidIndex].size(), fluidParticles[fluidIndex].data(), fluidVelocities[fluidIndex].data(), 0.0);
FileSystem::writeMD5File(fileName, md5FileName);
}
// transform particles
for (unsigned int j = 0; j < (unsigned int)fluidParticles[fluidIndex].size(); j )
fluidParticles[fluidIndex][j] = m_scene.fluidModels[i]->rotation * fluidParticles[fluidIndex][j] m_scene.fluidModels[i]->translation;
}
}
else
{
if (!PartioReaderWriter::readParticles(fileName, m_scene.fluidModels[i]->translation, m_scene.fluidModels[i]->rotation, m_scene.fluidModels[i]->scale[0], fluidParticles[fluidIndex], fluidVelocities[fluidIndex]))
LOG_ERR << "File not found: " << fileName;
}
Simulation::getCurrent()->setValue(Simulation::PARTICLE_RADIUS, m_scene.particleRadius);
}
unsigned int nParticles = 0;
for (auto it = fluidIDs.begin(); it != fluidIDs.end(); it )
{
const unsigned int index = it->second;
unsigned int maxEmitterParticles = 10000;
for (auto material : m_scene.materials)
{
if (material->id == it->first)
maxEmitterParticles = material->maxEmitterParticles;
}
sim->addFluidModel(it->first, (unsigned int)fluidParticles[index].size(), fluidParticles[index].data(), fluidVelocities[index].data(), maxEmitterParticles);
nParticles = (unsigned int)fluidParticles[index].size();
}
m_colorField.resize(sim->numberOfFluidModels(), "velocity");
m_colorMapType.resize(sim->numberOfFluidModels(), 0);
m_renderMinValue.resize(sim->numberOfFluidModels(), 0.0);
m_renderMaxValue.resize(sim->numberOfFluidModels(), 5.0);
LOG_INFO << "Number of fluid particles: " << nParticles;
}
void SimulatorBase::createEmitters()
{
Simulation *sim = Simulation::getCurrent();
//////////////////////////////////////////////////////////////////////////
// emitters
//////////////////////////////////////////////////////////////////////////
for (unsigned int i = 0; i < m_scene.emitters.size(); i )
{
SceneLoader::EmitterData *ed = m_scene.emitters[i];
FluidModel *model = nullptr;
unsigned int j;
for (j = 0; j < sim->numberOfFluidModels(); j )
{
model = sim->getFluidModel(j);
if (model->getId() == ed->id)
break;
}
if (j < sim->numberOfFluidModels())
{
if (sim->is2DSimulation())
{
// in 2D convert all emitters to box emitters
// and set width to 1
if (ed->type == 1)
ed->height = ed->width;
ed->width = 1;
ed->type = 0;
}
model->getEmitterSystem()->addEmitter(
ed->width, ed->height,
ed->x, ed->rotation,
ed->velocity,
ed->type);
// Generate boundary geometry around emitters
Emitter *emitter = model->getEmitterSystem()->getEmitters().back();
SceneLoader::BoundaryData *emitterBoundary = new SceneLoader::BoundaryData();
emitterBoundary->dynamic = false;
emitterBoundary->isWall = false;
emitterBoundary->color = { 0.2f, 0.2f, 0.2f, 1.0f };
emitterBoundary->rotation = ed->rotation;
const Real supportRadius = sim->getSupportRadius();
const Vector3r & emitDir = ed->rotation.col(0);
emitterBoundary->scale = Emitter::getSize(static_cast<Real>(ed->width), static_cast<Real>(ed->height), ed->type);
const Vector3r pos = ed->x;
emitterBoundary->translation = pos;
emitterBoundary->samplesFile = "";
emitterBoundary->mapInvert = false;
emitterBoundary->mapResolution = Eigen::Matrix<unsigned int, 3, 1, Eigen::DontAlign>(20, 20, 20);
emitterBoundary->mapThickness = 0.0;
if (sim->is2DSimulation())
emitterBoundary->scale[2] = 2 * supportRadius;
if (ed->type == 0)
emitterBoundary->meshFile = FileSystem::normalizePath(getExePath() "/resources/emitter_boundary/EmitterBox.obj");
else if (ed->type == 1)
emitterBoundary->meshFile = FileSystem::normalizePath(getExePath() "/resources/emitter_boundary/EmitterCylinder.obj");
m_scene.boundaryModels.push_back(emitterBoundary);
// reuse particles if they are outside of a bounding box
bool emitterReuseParticles = false;
Vector3r emitterBoxMin(-1.0, -1.0, -1.0);
Vector3r emitterBoxMax(1.0, 1.0, 1.0);
for (auto material : m_scene.materials)
{
if (material->id == model->getId())
{
emitterReuseParticles = material->emitterReuseParticles;
emitterBoxMin = material->emitterBoxMin;
emitterBoxMax = material->emitterBoxMax;
}
}
if (emitterReuseParticles)
model->getEmitterSystem()->enableReuseParticles(emitterBoxMin, emitterBoxMax);
emitter->setEmitStartTime(ed->emitStartTime);
emitter->setEmitEndTime(ed->emitEndTime);
}
}
}
void SimulatorBase::createAnimationFields()
{
Simulation *sim = Simulation::getCurrent();
//////////////////////////////////////////////////////////////////////////
// animation fields
//////////////////////////////////////////////////////////////////////////
for (unsigned int i = 0; i < m_scene.animatedFields.size(); i )
{
SceneLoader::AnimationFieldData *data = m_scene.animatedFields[i];
sim->getAnimationFieldSystem()->addAnimationField(
data->particleFieldName,
data->x, data->rotation, data->scale,
data->expression,
data->shapeType);
AnimationField *field = sim->getAnimationFieldSystem()->getAnimationFields().back();
field->setStartTime(data->startTime);
field->setEndTime(data->endTime);
}
}
void SimulatorBase::createFluidBlocks(std::map<std::string, unsigned int> &fluidIDs, std::vector<std::vector<Vector3r>> &fluidParticles, std::vector<std::vector<Vector3r>> &fluidVelocities)
{
for (unsigned int i = 0; i < m_scene.fluidBlocks.size(); i )
{
const unsigned int fluidIndex = fluidIDs[m_scene.fluidBlocks[i]->id];
const Real diam = static_cast<Real>(2.0)*m_scene.particleRadius;
Real xshift = diam;
Real yshift = diam;
const Real eps = static_cast<Real>(1.0e-9);
if (m_scene.fluidBlocks[i]->mode == 1)
yshift = sqrt(static_cast<Real>(3.0)) * m_scene.particleRadius eps;
else if (m_scene.fluidBlocks[i]->mode == 2)
{
xshift = sqrt(static_cast<Real>(6.0)) * diam / static_cast<Real>(3.0) eps;
yshift = sqrt(static_cast<Real>(3.0)) * m_scene.particleRadius eps;
}
Vector3r diff = m_scene.fluidBlocks[i]->box.m_maxX - m_scene.fluidBlocks[i]->box.m_minX;
if (m_scene.fluidBlocks[i]->mode == 1)
{
diff[0] -= diam;
diff[2] -= diam;
}
else if (m_scene.fluidBlocks[i]->mode == 2)
{
diff[0] -= xshift;
diff[2] -= diam;
}
const int stepsX = (int)round(diff[0] / xshift) - 1;
const int stepsY = (int)round(diff[1] / yshift) - 1;
int stepsZ = (int)round(diff[2] / diam) - 1;
Vector3r start = m_scene.fluidBlocks[i]->box.m_minX static_cast<Real>(2.0)*m_scene.particleRadius*Vector3r::Ones();
fluidParticles[fluidIndex].reserve(fluidParticles[fluidIndex].size() stepsX*stepsY*stepsZ);
fluidVelocities[fluidIndex].resize(fluidVelocities[fluidIndex].size() stepsX*stepsY*stepsZ, m_scene.fluidBlocks[i]->initialVelocity);
if (Simulation::getCurrent()->is2DSimulation())
{
stepsZ = 1;
start[2] = 0.0;
}
for (int j = 0; j < stepsX; j )
{
for (int k = 0; k < stepsY; k )
{
for (int l = 0; l < stepsZ; l )
{
Vector3r currPos = Vector3r(j*xshift, k*yshift, l*diam) start;
if (m_scene.fluidBlocks[i]->mode == 1)
{
if (k % 2 == 0)
currPos = Vector3r(0, 0, m_scene.particleRadius);
else
currPos = Vector3r(m_scene.particleRadius, 0, 0);
}
else if (m_scene.fluidBlocks[i]->mode == 2)
{
currPos = Vector3r(0, 0, m_scene.particleRadius);
Vector3r shift_vec(0, 0, 0);
if ((j % 2) && !Simulation::getCurrent()->is2DSimulation())
{
shift_vec[2] = diam / (static_cast<Real>(2.0) * (k % 2 ? -1 : 1));
}
if (k % 2 == 0)
{
shift_vec[0] = xshift / static_cast<Real>(2.0);
}
currPos = shift_vec;
}
fluidParticles[fluidIndex].push_back(currPos);
}
}
}
}
}
void SimulatorBase::particleInfo(std::vector<std::vector<unsigned int>> &particles)
{
Simulation *sim = Simulation::getCurrent();
const int maxWidth = 25;
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
if (particles[i].size() > 0)
{
LOG_INFO << "---------------------------------------------------------------------------";
LOG_INFO << model->getId();
LOG_INFO << "---------------------------------------------------------------------------";
}
for (unsigned int j = 0; j < particles[i].size(); j )
{
unsigned int index = particles[i][j];
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << "Index:" << index;
for (unsigned int k = 0; k < model->numberOfFields(); k )
{
const FieldDescription &field = model->getField(k);
if (field.type == Scalar)
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name ":" << *((Real*) field.getFct(index));
else if (field.type == UInt)
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name ":" << *((unsigned int*)field.getFct(index));
else if (field.type == Vector3)
{
Eigen::Map<Vector3r> vec((Real*) field.getFct(index));
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name ":" << vec.transpose();
}
else if (field.type == Vector6)
{
Eigen::Map<Vector6r> vec((Real*)field.getFct(index));
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name ":" << vec.transpose();
}
else if (field.type == Matrix3)
{
Eigen::Map<Matrix3r> mat((Real*)field.getFct(index));
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name ":" << mat.row(0);
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << " " << mat.row(1);
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << " " << mat.row(2);
}
else if (field.type == Matrix6)
{
Eigen::Map<Matrix6r> mat((Real*) field.getFct(index));
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name ":" << mat.row(0);
for (unsigned int k = 1; k < 6; k )
LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << " " << mat.row(k);
}
}
LOG_INFO << "---------------------------------------------------------------------------n";
}
}
}
void SimulatorBase::rigidBodyExport()
{
std::string exportPath = FileSystem::normalizePath(m_outputPath "/rigid_bodies");
std::string exportPathVTK = FileSystem::normalizePath(m_outputPath "/vtk");
if (m_enableRigidBodyExport)
{
FileSystem::makeDirs(exportPath);
writeRigidBodiesBIN(exportPath);
}
if (m_enableRigidBodyVTKExport)
{
FileSystem::makeDirs(exportPathVTK);
writeRigidBodiesVTK(exportPathVTK);
}
}
void SimulatorBase::particleExport()
{
std::string partioExportPath = FileSystem::normalizePath(m_outputPath "/partio");
std::string vtkExportPath = FileSystem::normalizePath(m_outputPath "/vtk");
if (m_enablePartioExport)
FileSystem::makeDirs(partioExportPath);
if (m_enableVTKExport)
FileSystem::makeDirs(vtkExportPath);
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
std::string fileName = "ParticleData";
fileName = fileName "_" model->getId() "_" std::to_string(m_frameCounter);
if (m_enablePartioExport)
{
std::string exportFileName = FileSystem::normalizePath(partioExportPath "/" fileName);
writeParticlesPartio(exportFileName ".bgeo", model);
}
if (m_enableVTKExport)
{
std::string exportFileName = FileSystem::normalizePath(vtkExportPath "/" fileName);
writeParticlesVTK(exportFileName ".vtk", model);
}
}
}
void SimulatorBase::writeParticlesPartio(const std::string &fileName, FluidModel *model)
{
Partio::ParticlesDataMutable& particleData = *Partio::create();
Partio::ParticleAttribute posAttr = particleData.addAttribute("position", Partio::VECTOR, 3);
Partio::ParticleAttribute idAttr = particleData.addAttribute("id", Partio::INT, 1);
// add attributes
std::vector<std::string> attributes;
StringTools::tokenize(m_particleAttributes, attributes, ";");
std::map<unsigned int, int> attrMap;
std::map<unsigned int, Partio::ParticleAttribute> partioAttrMap;
for (unsigned int i = 0; i < attributes.size(); i )
{
// position is exported anyway
if (attributes[i] == "position")
{
attrMap[i] = -1;
continue;
}
bool found = false;
for (unsigned int j = 0; j < model->numberOfFields(); j )
{
const FieldDescription &field = model->getField(j);
if (field.name == attributes[i])
{
found = true;
if (field.type == Scalar)
{
attrMap[i] = j;
partioAttrMap[i] = particleData.addAttribute(attributes[i].c_str(), Partio::FLOAT, 1);
}
else if (field.type == UInt)
{
attrMap[i] = j;
partioAttrMap[i] = particleData.addAttribute(attributes[i].c_str(), Partio::INT, 1);
}
else if (field.type == Vector3)
{
attrMap[i] = j;
partioAttrMap[i] = particleData.addAttribute(attributes[i].c_str(), Partio::VECTOR, 3);
}
else
{
attrMap[i] = -1;
LOG_WARN << "Only scalar and vector fields are currently supported by the partio exporter.";
}
break;
}
}
if (!found)
{
attrMap[i] = -1;
LOG_WARN << "Unknown field cannot be exported in partio file: " << attributes[i];
}
}
const unsigned int numParticles = model->numActiveParticles();
for (unsigned int i = 0; i < numParticles; i )
{
Partio::ParticleIndex index = particleData.addParticle();
float* pos = particleData.dataWrite<float>(posAttr, index);
int* id = particleData.dataWrite<int>(idAttr, index);
const Vector3r &x = model->getPosition(i);
pos[0] = (float)x[0];
pos[1] = (float)x[1];
pos[2] = (float)x[2];
id[0] = model->getParticleId(i);
for (unsigned int j = 0; j < attributes.size(); j )
{
const int fieldIndex = attrMap[j];
if (fieldIndex != -1)
{
const FieldDescription &field = model->getField(fieldIndex);
if (field.type == FieldType::Scalar)
{
float* val = particleData.dataWrite<float>(partioAttrMap[j], index);
*val = (float) *((Real*) field.getFct(i));
}
else if (field.type == FieldType::UInt)
{
int* val = particleData.dataWrite<int>(partioAttrMap[j], index);
*val = (int) *((unsigned int*)field.getFct(i));
}
else if (field.type == FieldType::Vector3)
{
float* val = particleData.dataWrite<float>(partioAttrMap[j], index);
Eigen::Map<Vector3r> vec((Real*) field.getFct(i));
val[0] = (float)vec[0];
val[1] = (float)vec[1];
val[2] = (float)vec[2];
}
}
}
}
Partio::write(fileName.c_str(), particleData, true);
particleData.release();
}
void SPH::SimulatorBase::writeParticlesVTK(const std::string &fileName, FluidModel *model)
{
const unsigned int numParticles = model->numActiveParticles();
if (0 == numParticles)
return;
#ifdef USE_DOUBLE
const char * real_str = " doublen";
#else
const char * real_str = " floatn";
#endif
// Open the file
std::ofstream outfile{ fileName, std::ios::binary };
if (!outfile.is_open())
{
LOG_WARN << "Cannot open a file to save VTK particles.";
return;
}
outfile << "# vtk DataFile Version 4.1n";
outfile << "SPlisHSPlasH particle datan"; // title of the data set, (any string up to 256 characters n)
outfile << "BINARYn";
outfile << "DATASET UNSTRUCTURED_GRIDn";
// add attributes
std::vector<std::string> attributes;
StringTools::tokenize(m_particleAttributes, attributes, ";");
//////////////////////////////////////////////////////////////////////////
// positions and ids exported anyways
attributes.erase(
std::remove_if(attributes.begin(), attributes.end(), [](const std::string&s) { return (s == "position" || s == "id"); }),
attributes.end());
//////////////////////////////////////////////////////////////////////////
// export position attribute as POINTS
{
std::vector<Vector3r> positions;
positions.reserve(numParticles);
for (unsigned int i = 0u; i < numParticles; i )
positions.emplace_back(model->getPosition(i));
// swap endianess
for (unsigned int i = 0; i < numParticles; i )
for (unsigned int c = 0; c < 3; c )
swapByteOrder(&positions[i][c]);
// export to vtk
outfile << "POINTS " << numParticles << real_str;
outfile.write(reinterpret_cast<char*>(positions[0].data()), 3 * numParticles * sizeof(Real));
outfile << "n";
}
//////////////////////////////////////////////////////////////////////////
// export particle IDs as CELLS
{
std::vector<Eigen::Vector2i> cells;
cells.reserve(numParticles);
unsigned int nodes_per_cell_swapped = 1;
swapByteOrder(&nodes_per_cell_swapped);
for (unsigned int i = 0u; i < numParticles; i )
{
unsigned int idSwapped = model->getParticleId(i);
swapByteOrder(&idSwapped);
cells.emplace_back(nodes_per_cell_swapped, idSwapped);
}
// particles are cells with one element and the index of the particle
outfile << "CELLS " << numParticles << " " << 2 * numParticles << "n";
outfile.write(reinterpret_cast<char*>(cells[0].data()), 2 * numParticles * sizeof(unsigned int));
outfile << "n";
}
//////////////////////////////////////////////////////////////////////////
// export cell types
{
// the type of a particle cell is always 1
std::vector<int> cellTypes;
int cellTypeSwapped = 1;
swapByteOrder(&cellTypeSwapped);
cellTypes.resize(numParticles, cellTypeSwapped);
outfile << "CELL_TYPES " << numParticles << "n";
outfile.write(reinterpret_cast<char*>(cellTypes.data()), numParticles * sizeof(int));
outfile << "n";
}
//////////////////////////////////////////////////////////////////////////
// write additional attributes as per-particle data
{
outfile << "POINT_DATA " << numParticles << "n";
// write IDs
outfile << "SCALARS id unsigned_int 1n";
outfile << "LOOKUP_TABLE id_tablen";
// copy data
std::vector<unsigned int> attrData;
attrData.reserve(numParticles);
for (unsigned int i = 0u; i < numParticles; i )
attrData.emplace_back(model->getParticleId(i));
// swap endianess
for (unsigned int i = 0; i < numParticles; i )
swapByteOrder(&attrData[i]);
// export to vtk
outfile.write(reinterpret_cast<char*>(attrData.data()), numParticles * sizeof(unsigned int));
outfile << "n";
}
//////////////////////////////////////////////////////////////////////////
// per point fields (all attributes except for positions)
const auto numFields = attributes.size();
outfile << "FIELD FieldData " << std::to_string(numFields) << "n";
// iterate over attributes
for (const std::string & a : attributes)
{
const FieldDescription & field = model->getField(a);
std::string attrNameVTK;
std::regex_replace(std::back_inserter(attrNameVTK), a.begin(), a.end(), std::regex("\s "), "_");
if (field.type == FieldType::Scalar)
{
// write header information
outfile << attrNameVTK << " 1 " << numParticles << real_str;
// copy data
std::vector<Real> attrData;
attrData.reserve(numParticles);
for (unsigned int i = 0u; i < numParticles; i )
attrData.emplace_back(*((Real*) field.getFct(i)));
// swap endianess
for (unsigned int i = 0; i < numParticles; i )
swapByteOrder(&attrData[i]);
// export to vtk
outfile.write(reinterpret_cast<char*>(attrData.data()), numParticles * sizeof(Real));
}
else if (field.type == FieldType::Vector3)
{
// write header information
outfile << attrNameVTK << " 3 " << numParticles << real_str;
// copy from partio data
std::vector<Vector3r> attrData;
attrData.reserve(numParticles);
for (unsigned int i = 0u; i < numParticles; i )
attrData.emplace_back((Real*) field.getFct(i));
// swap endianess
for (unsigned int i = 0; i < numParticles; i )
for (unsigned int c = 0; c < 3; c )
swapByteOrder(&attrData[i][c]);
// export to vtk
outfile.write(reinterpret_cast<char*>(attrData[0].data()), 3 * numParticles * sizeof(Real));
}
// TODO support other field types
else
{
LOG_WARN << "Skipping attribute " << a << ", because it is of unsupported typen";
continue;
}
// end of block
outfile << "n";
}
outfile.close();
}
void SimulatorBase::writeRigidBodiesBIN(const std::string &exportPath)
{
std::string fileName = "rb_data_";
fileName = fileName std::to_string(m_frameCounter) ".bin";
std::string exportFileName = FileSystem::normalizePath(exportPath "/" fileName);
Simulation *sim = Simulation::getCurrent();
const unsigned int nBoundaryModels = sim->numberOfBoundaryModels();
std::string scene_path = FileSystem::getFilePath(getSceneFile());
// check if we have a static model
bool isStatic = true;
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *bm = sim->getBoundaryModel(i);
if (bm->getRigidBodyObject()->isDynamic())
{
isStatic = false;
break;
}
}
BinaryFileWriter binWriter;
if (m_isFirstFrame || !isStatic)
binWriter.openFile(exportFileName.c_str());
if (m_isFirstFrame)
{
binWriter.write(nBoundaryModels);
for (unsigned int i = 0; i < m_scene.boundaryModels.size(); i )
{
std::string meshFileName = m_scene.boundaryModels[i]->meshFile;
if (FileSystem::isRelativePath(meshFileName))
meshFileName = FileSystem::normalizePath(scene_path "/" meshFileName);
const string fileName = Utilities::FileSystem::getFileNameWithExt(meshFileName);
binWriter.write(fileName);
Eigen::Vector3f s = m_scene.boundaryModels[i]->scale.template cast<float>();
binWriter.writeMatrix(s);
std::string targetFilePath = exportPath "/" fileName;
if (!Utilities::FileSystem::fileExists(targetFilePath))
{
Utilities::FileSystem::copyFile(meshFileName, targetFilePath);
}
binWriter.write((char)m_scene.boundaryModels[i]->isWall);
binWriter.writeMatrix(m_scene.boundaryModels[i]->color);
}
}
if (m_isFirstFrame || !isStatic)
{
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *bm = sim->getBoundaryModel(i);
const Vector3r &x = bm->getRigidBodyObject()->getWorldSpacePosition();
const Eigen::Vector3f x_f = x.template cast<float>();
binWriter.writeMatrix(x_f);
const Matrix3r &R = bm->getRigidBodyObject()->getWorldSpaceRotation();
const Eigen::Matrix3f R_f = R.template cast<float>();
//const Eigen::Matrix3f RT = R.transpose().template cast<float>();
binWriter.writeMatrix(R_f);
}
binWriter.closeFile();
}
m_isFirstFrame = false;
}
void SimulatorBase::writeRigidBodiesVTK(const std::string &exportPath)
{
Simulation *sim = Simulation::getCurrent();
const unsigned int nBoundaryModels = sim->numberOfBoundaryModels();
// check if we have a static model
bool isStatic = true;
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *bm = sim->getBoundaryModel(i);
if (bm->getRigidBodyObject()->isDynamic())
{
isStatic = false;
break;
}
}
#ifdef USE_DOUBLE
const char * real_str = " doublen";
#else
const char * real_str = " floatn";
#endif
if (m_isFirstFrameVTK || !isStatic)
{
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
std::string fileName = "rb_data_";
fileName = fileName std::to_string(i) "_" std::to_string(m_frameCounter) ".vtk";
std::string exportFileName = FileSystem::normalizePath(exportPath "/" fileName);
// Open the file
std::ofstream outfile(exportFileName, std::ios::binary);
if (!outfile)
{
LOG_WARN << "Cannot open a file to save VTK mesh.";
return;
}
// Header
outfile << "# vtk DataFile Version 4.2n";
outfile << "SPlisHSPlasH mesh datan";
outfile << "BINARYn";
outfile << "DATASET UNSTRUCTURED_GRIDn";
BoundaryModel *bm = sim->getBoundaryModel(i);
const std::vector<Vector3r> &vertices = bm->getRigidBodyObject()->getVertices();
const std::vector<unsigned int> &faces = bm->getRigidBodyObject()->getFaces();
int n_vertices = (int)vertices.size();
int n_triangles = (int)faces.size() / 3;
// Vertices
{
std::vector<Vector3r> positions;
positions.reserve(n_vertices);
for (int j = 0u; j < n_vertices; j )
{
Vector3r x = vertices[j];
swapByteOrder(&x[0]);
swapByteOrder(&x[1]);
swapByteOrder(&x[2]);
positions.emplace_back(x);
}
// export to vtk
outfile << "POINTS " << n_vertices << real_str;
outfile.write(reinterpret_cast<char*>(positions[0].data()), 3 * n_vertices * sizeof(Real));
outfile << "n";
}
// Connectivity
{
std::vector<int> connectivity_to_write;
connectivity_to_write.reserve(4 * n_triangles);
for (int tri_i = 0; tri_i < n_triangles; tri_i )
{
int val = 3;
swapByteOrder(&val);
connectivity_to_write.push_back(val);
val = faces[3 * tri_i 0];
swapByteOrder(&val);
connectivity_to_write.push_back(val);
val = faces[3 * tri_i 1];
swapByteOrder(&val);
connectivity_to_write.push_back(val);
val = faces[3 * tri_i 2];
swapByteOrder(&val);
connectivity_to_write.push_back(val);
}
// export to vtk
outfile << "CELLS " << n_triangles << " " << 4 * n_triangles << "n";
outfile.write(reinterpret_cast<char*>(&connectivity_to_write[0]), connectivity_to_write.size() * sizeof(int));
outfile << "n";
}
// Cell types
{
outfile << "CELL_TYPES " << n_triangles << "n";
int cell_type_swapped = 5;
swapByteOrder(&cell_type_swapped);
std::vector<int> cell_type_arr(n_triangles, cell_type_swapped);
outfile.write(reinterpret_cast<char*>(&cell_type_arr[0]), cell_type_arr.size() * sizeof(int));
outfile << "n";
}
outfile.close();
}
}
m_isFirstFrameVTK = false;
}
void SimulatorBase::step()
{
if (TimeManager::getCurrent()->getTime() >= m_nextFrameTime)
{
m_nextFrameTime = static_cast<Real>(1.0) / m_framesPerSecond;
if (m_enablePartioExport || m_enableVTKExport)
particleExport();
if (m_enableRigidBodyExport || m_enableRigidBodyVTKExport)
rigidBodyExport();
m_frameCounter ;
}
if (TimeManager::getCurrent()->getTime() >= m_nextFrameTimeState)
{
m_nextFrameTimeState = static_cast<Real>(1.0) / m_framesPerSecondState;
if (m_enableStateExport)
saveState();
}
#ifdef DL_OUTPUT
if (TimeManager::getCurrent()->getTime() >= m_nextTiming)
{
LOG_INFO << "---------------------------------------------------------------------------";
LOG_INFO << "Time: " << TimeManager::getCurrent()->getTime();
Timing::printAverageTimes();
Timing::printTimeSums();
Counting::printAverageCounts();
Counting::printCounterSums();
m_nextTiming = 1.0;
}
#endif
}
void SimulatorBase::updateBoundaryParticles(const bool forceUpdate = false)
{
Simulation *sim = Simulation::getCurrent();
SceneLoader::Scene &scene = getScene();
const unsigned int nObjects = sim->numberOfBoundaryModels();
for (unsigned int i = 0; i < nObjects; i )
{
BoundaryModel_Akinci2012 *bm = static_cast<BoundaryModel_Akinci2012*>(sim->getBoundaryModel(i));
RigidBodyObject *rbo = bm->getRigidBodyObject();
if (rbo->isDynamic() || forceUpdate)
{
#pragma omp parallel default(shared)
{
#pragma omp for schedule(static)
for (int j = 0; j < (int)bm->numberOfParticles(); j )
{
bm->getPosition(j) = rbo->getRotation() * bm->getPosition0(j) rbo->getPosition();
if (rbo->isDynamic())
bm->getVelocity(j) = rbo->getAngularVelocity().cross(bm->getPosition(j) - rbo->getPosition()) rbo->getVelocity();
else
bm->getVelocity(j).setZero();
}
}
#ifdef GPU_NEIGHBORHOOD_SEARCH
// copy the particle data to the GPU
if (forceUpdate)
sim->getNeighborhoodSearch()->update_point_sets();
#endif
}
}
}
void SPH::SimulatorBase::updateDMVelocity()
{
Simulation *sim = Simulation::getCurrent();
SceneLoader::Scene &scene = getScene();
const unsigned int nObjects = sim->numberOfBoundaryModels();
for (unsigned int i = 0; i < nObjects; i )
{
BoundaryModel_Koschier2017 *bm = static_cast<BoundaryModel_Koschier2017*>(sim->getBoundaryModel(i));
RigidBodyObject *rbo = bm->getRigidBodyObject();
if (rbo->isDynamic())
{
const Real maxDist = bm->getMaxDist();
const Vector3r x(maxDist, 0.0, 0.0);
const Vector3r vel = rbo->getAngularVelocity().cross(x) rbo->getVelocity();
bm->setMaxVel(vel.norm());
}
}
}
void SPH::SimulatorBase::updateVMVelocity()
{
Simulation *sim = Simulation::getCurrent();
SceneLoader::Scene &scene = getScene();
const unsigned int nObjects = sim->numberOfBoundaryModels();
for (unsigned int i = 0; i < nObjects; i )
{
BoundaryModel_Bender2019 *bm = static_cast<BoundaryModel_Bender2019*>(sim->getBoundaryModel(i));
RigidBodyObject *rbo = bm->getRigidBodyObject();
if (rbo->isDynamic())
{
const Real maxDist = bm->getMaxDist();
const Vector3r x(maxDist, 0.0, 0.0);
const Vector3r vel = rbo->getAngularVelocity().cross(x) rbo->getVelocity();
bm->setMaxVel(vel.norm());
}
}
}
std::string SimulatorBase::real2String(const Real r)
{
string str = to_string(r);
str.erase(str.find_last_not_of('0') 1, std::string::npos);
str.erase(str.find_last_not_of('.') 1, std::string::npos);
return str;
}
void SPH::SimulatorBase::saveState()
{
std::string stateFilePath = FileSystem::normalizePath(m_outputPath "/state");
FileSystem::makeDirs(stateFilePath);
string md5Str = FileSystem::getFileMD5(m_sceneFile);
Simulation *sim = Simulation::getCurrent();
const Real time = TimeManager::getCurrent()->getTime();
const std::string timeStr = real2String(time);
// Save additional data
BinaryFileWriter binWriter;
std::string exportFileName = FileSystem::normalizePath(stateFilePath "/state_" timeStr);
binWriter.openFile(exportFileName ".bin");
binWriter.write(md5Str);
binWriter.write(m_nextFrameTime);
binWriter.write(m_nextFrameTimeState);
binWriter.write(m_frameCounter);
binWriter.write(m_isFirstFrame);
writeParameterState(binWriter);
TimeManager::getCurrent()->saveState(binWriter);
Simulation::getCurrent()->saveState(binWriter);
// fluid models
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
std::string fileName = "particle";
fileName = fileName "_" model->getId(); // "_" std::to_string(m_frameCounter);
// Save particle data
std::string expFileName = FileSystem::normalizePath(exportFileName "_" fileName);
writeFluidParticlesState(expFileName ".bgeo", model);
}
// boundary models
if (m_firstState)
{
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *model = sim->getBoundaryModel(i);
std::string fileName = "state_boundary";
fileName = fileName "_" to_string(i); // "_" std::to_string(m_frameCounter);
// Save particle data
std::string expFileName = FileSystem::normalizePath(stateFilePath "/" fileName);
writeBoundaryState(expFileName ".bgeo", model);
}
m_firstState = false;
}
// dynamic bodies
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *bm = sim->getBoundaryModel(i);
if (bm->getRigidBodyObject()->isDynamic())
{
binWriter.writeMatrix(bm->getRigidBodyObject()->getPosition());
binWriter.writeMatrix(bm->getRigidBodyObject()->getRotation());
binWriter.writeMatrix(bm->getRigidBodyObject()->getVelocity());
binWriter.writeMatrix(bm->getRigidBodyObject()->getAngularVelocity());
}
}
binWriter.closeFile();
LOG_INFO << "Saved state: " << exportFileName ".bin";
}
#ifdef WIN32
void SPH::SimulatorBase::loadStateDialog()
{
const std::string stateFilePath = FileSystem::normalizePath(m_outputPath "/state");
const std::string stateFileName = FileSystem::fileDialog(0, stateFilePath, "*.bin");
if (stateFileName == "")
return;
loadState(stateFileName);
}
#endif
void SPH::SimulatorBase::loadState(const std::string &stateFile)
{
Simulation *sim = Simulation::getCurrent();
string md5Str = FileSystem::getFileMD5(m_sceneFile);
// Load additional data
BinaryFileReader binReader;
if (!binReader.openFile(stateFile))
return;
// check if scene file has changed
std::string md5StrState;
binReader.read(md5StrState);
binReader.read(m_nextFrameTime);
binReader.read(m_nextFrameTimeState);
binReader.read(m_frameCounter);
binReader.read(m_isFirstFrame);
readParameterState(binReader);
if (md5Str != md5StrState)
LOG_WARN << "State was stored for another scene file.";
TimeManager::getCurrent()->loadState(binReader);
Simulation::getCurrent()->loadState(binReader);
const std::string importFilePath = FileSystem::getFilePath(stateFile);
const std::string importFileName = FileSystem::getFileName(stateFile);
// fluid models
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
std::string fileName = "particle";
fileName = fileName "_" model->getId(); // "_" std::to_string(m_frameCounter);
std::string impFileName = FileSystem::normalizePath(importFilePath "/" importFileName "_" fileName);
readFluidParticlesState(impFileName ".bgeo", model);
}
// boundary models
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *model = sim->getBoundaryModel(i);
std::string fileName = "state_boundary";
fileName = fileName "_" to_string(i); // "_" std::to_string(m_frameCounter);
// Save particle data
std::string impFileName = FileSystem::normalizePath(importFilePath "/" fileName);
readBoundaryState(impFileName ".bgeo", model);
}
// dynamic bodies
bool dynamic = false;
for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i )
{
BoundaryModel *bm = sim->getBoundaryModel(i);
if (bm->getRigidBodyObject()->isDynamic())
{
Vector3r x;
binReader.readMatrix(x);
bm->getRigidBodyObject()->setPosition(x);
Matrix3r R;
binReader.readMatrix(R);
bm->getRigidBodyObject()->setRotation(R);
Vector3r v;
binReader.readMatrix(v);
bm->getRigidBodyObject()->setVelocity(v);
binReader.readMatrix(v);
bm->getRigidBodyObject()->setAngularVelocity(v);
dynamic = true;
}
}
if (dynamic)
{
if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
updateBoundaryParticles(true);
else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
updateDMVelocity();
else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
updateVMVelocity();
}
binReader.closeFile();
sim->performNeighborhoodSearchSort();
}
void SimulatorBase::writeParameterState(BinaryFileWriter &binWriter)
{
writeParameterObjectState(binWriter, this);
writeParameterObjectState(binWriter, Simulation::getCurrent());
writeParameterObjectState(binWriter, Simulation::getCurrent()->getTimeStep());
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
writeParameterObjectState(binWriter, model);
writeParameterObjectState(binWriter, (ParameterObject*)model->getDragBase());
writeParameterObjectState(binWriter, (ParameterObject*)model->getSurfaceTensionBase());
writeParameterObjectState(binWriter, (ParameterObject*)model->getViscosityBase());
writeParameterObjectState(binWriter, (ParameterObject*)model->getVorticityBase());
writeParameterObjectState(binWriter, (ParameterObject*)model->getElasticityBase());
binWriter.write(getColorField(model->getPointSetIndex()));
binWriter.write(getColorMapType(model->getPointSetIndex()));
binWriter.write(getRenderMinValue(model->getPointSetIndex()));
binWriter.write(getRenderMaxValue(model->getPointSetIndex()));
}
}
void SimulatorBase::writeParameterObjectState(BinaryFileWriter &binWriter, GenParam::ParameterObject *paramObj)
{
if (paramObj == nullptr)
return;
const unsigned int numParams = paramObj->numParameters();
for (unsigned int i = 0; i < numParams; i )
{
ParameterBase *paramBase = paramObj->getParameter(i);
if (paramBase->getType() == RealParameterType)
binWriter.write(static_cast<NumericParameter<Real>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::UINT32)
binWriter.write(static_cast<NumericParameter<unsigned int>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::UINT16)
binWriter.write(static_cast<NumericParameter<unsigned short>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::UINT8)
binWriter.write(static_cast<NumericParameter<unsigned char>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::INT32)
binWriter.write(static_cast<NumericParameter<int>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::INT16)
binWriter.write(static_cast<NumericParameter<short>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::INT8)
binWriter.write(static_cast<NumericParameter<char>*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::ENUM)
binWriter.write(static_cast<EnumParameter*>(paramBase)->getValue());
else if (paramBase->getType() == ParameterBase::BOOL)
binWriter.write(static_cast<BoolParameter*>(paramBase)->getValue());
else if (paramBase->getType() == RealVectorParameterType)
{
VectorParameter<Real> *vec = static_cast<VectorParameter<Real>*>(paramBase);
binWriter.writeBuffer((char*)vec->getValue(), vec->getDim() * sizeof(Real));;
}
else if (paramBase->getType() == ParameterBase::STRING)
binWriter.write(static_cast<StringParameter*>(paramBase)->getValue());
}
}
void SimulatorBase::readParameterState(BinaryFileReader &binReader)
{
readParameterObjectState(binReader, this);
readParameterObjectState(binReader, Simulation::getCurrent());
readParameterObjectState(binReader, Simulation::getCurrent()->getTimeStep());
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i )
{
FluidModel *model = sim->getFluidModel(i);
readParameterObjectState(binReader, model);
readParameterObjectState(binReader, (ParameterObject*)model->getDragBase());
readParameterObjectState(binReader, (ParameterObject*)model->getSurfaceTensionBase());
readParameterObjectState(binReader, (ParameterObject*)model->getViscosityBase());
readParameterObjectState(binReader, (ParameterObject*)model->getVorticityBase());
readParameterObjectState(binReader, (ParameterObject*)model->getElasticityBase());
std::string field;
binReader.read(field);
setColorField(model->getPointSetIndex(), field);
int type;
binReader.read(type);
setColorMapType(model->getPointSetIndex(), type);
Real v;
binReader.read(v);
setRenderMinValue(model->getPointSetIndex(), v);
binReader.read(v);
setRenderMaxValue(model->getPointSetIndex(), v);
}
}
void SimulatorBase::readParameterObjectState(BinaryFileReader &binReader, GenParam::ParameterObject *paramObj)
{
if (paramObj == nullptr)
return;
const unsigned int numParams = paramObj->numParameters();
for (unsigned int i = 0; i < numParams; i )
{
ParameterBase *paramBase = paramObj->getParameter(i);
if (paramBase->getType() == RealParameterType)
{
Real val;
binReader.read(val);
static_cast<NumericParameter<Real>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::UINT32)
{
unsigned int val;
binReader.read(val);
static_cast<NumericParameter<unsigned int>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::UINT16)
{
unsigned short val;
binReader.read(val);
static_cast<NumericParameter<unsigned short>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::UINT8)
{
unsigned char val;
binReader.read(val);
static_cast<NumericParameter<unsigned char>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::INT32)
{
int val;
binReader.read(val);
static_cast<NumericParameter<int>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::INT16)
{
short val;
binReader.read(val);
static_cast<NumericParameter<short>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::INT8)
{
char val;
binReader.read(val);
static_cast<NumericParameter<char>*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::ENUM)
{
int val;
binReader.read(val);
static_cast<EnumParameter*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == ParameterBase::BOOL)
{
bool val;
binReader.read(val);
static_cast<BoolParameter*>(paramBase)->setValue(val);
}
else if (paramBase->getType() == RealVectorParameterType)
{
VectorParameter<Real> *vec = static_cast<VectorParameter<Real>*>(paramBase);
binReader.readBuffer((char*)vec->getValue(), vec->getDim() * sizeof(Real));;
}
else if (paramBase->getType() == ParameterBase::STRING)
{
std::string val;
binReader.read(val);
static_cast<StringParameter*>(paramBase)->setValue(val);
}
}
}
void SimulatorBase::writeFluidParticlesState(const std::string &fileName, FluidModel *model)
{
Partio::ParticlesDataMutable& particleData = *Partio::create();
std::vector<std::pair<unsigned int, Partio::ParticleAttribute>> partioAttr;
for (unsigned int j = 0; j < model->numberOfFields(); j )
{
const FieldDescription &field = model->getField(j);
if (field.storeData)
{
// LOG_INFO << "Store field: " << field.name;
if (field.type == Scalar)
{
partioAttr.push_back({ j, particleData.addAttribute(field.name.c_str(), Partio::FLOAT, 1) });
}
else if (field.type == UInt)
{
partioAttr.push_back({ j, particleData.addAttribute(field.name.c_str(), Partio::INT, 1) });
}
else if (field.type == Vector3)
{
partioAttr.push_back({ j, particleData.addAttribute(field.name.c_str(), Partio::VECTOR, 3) });
}
else
{
LOG_WARN << "Only scalar and vector fields are currently supported by the partio exporter.";
}
}
}
const unsigned int numParticles = model->numActiveParticles();
for (unsigned int i = 0; i < numParticles; i )
{
Partio::ParticleIndex index = particleData.addParticle();
for (unsigned int j = 0; j < partioAttr.size(); j )
{
const unsigned int fieldIndex = partioAttr[j].first;
const Partio::ParticleAttribute &attr = partioAttr[j].second;
const FieldDescription &field = model->getField(fieldIndex);
if (field.type == FieldType::Scalar)
{
float* val = particleData.dataWrite<float>(attr, index);
*val = (float)* ((Real*) field.getFct(i));
}
else if (field.type == FieldType::UInt)
{
int* val = particleData.dataWrite<int>(attr, index);
*val = (int)* ((unsigned int*)field.getFct(i));
}
else if (field.type == FieldType::Vector3)
{
float* val = particleData.dataWrite<float>(attr, index);
Eigen::Map<Vector3r> vec((Real*) field.getFct(i));
val[0] = (float)vec[0];
val[1] = (float)vec[1];
val[2] = (float)vec[2];
}
}
}
Partio::write(fileName.c_str(), particleData, true);
particleData.release();
}
void SimulatorBase::readFluidParticlesState(const std::string &fileName, FluidModel *model)
{
if (!FileSystem::fileExists(fileName))
{
LOG_WARN << "File " << fileName << " does not exist.";
return;
}
Partio::ParticlesDataMutable* data = Partio::read(fileName.c_str());
if (!data)
{
LOG_WARN << "Partio file " << fileName << " not readable.";
return;
}
std::vector<std::pair<unsigned int, Partio::ParticleAttribute>> partioAttr;
for (int i = 0; i < data->numAttributes(); i )
{
Partio::ParticleAttribute attr;
data->attributeInfo(i, attr);
for (unsigned int j = 0; j < model->numberOfFields(); j )
{
const FieldDescription &field = model->getField(j);
if (field.name == attr.name)
{
//LOG_INFO << "Read field: " << field.name;
partioAttr.push_back({ j, attr});
}
}
}
for (unsigned int j = 0; j < partioAttr.size(); j )
{
const unsigned int fieldIndex = partioAttr[j].first;
const Partio::ParticleAttribute &attr = partioAttr[j].second;
const FieldDescription &field = model->getField(fieldIndex);
for (int i = 0; i < data->numParticles(); i )
{
if (field.type == FieldType::Scalar)
{
const float *value = data->data<float>(attr, i);
*((Real*) field.getFct(i)) = value[0];
}
else if (field.type == FieldType::UInt)
{
const int *value = data->data<int>(attr, i);
*((unsigned int*)field.getFct(i)) = value[0];
}
else if (field.type == FieldType::Vector3)
{
const float *value = data->data<float>(attr, i);
Eigen::Map<Vector3r> vec((Real*) field.getFct(i));
vec[0] = value[0];
vec[1] = value[1];
vec[2] = value[2];
}
}
}
data->release();
}
void SimulatorBase::writeBoundaryState(const std::string &fileName, BoundaryModel *bm)
{
Simulation *sim = Simulation::getCurrent();
if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
{
BoundaryModel_Akinci2012 *model = static_cast<BoundaryModel_Akinci2012*>(bm);
Partio::ParticlesDataMutable& particleData = *Partio::create();
const Partio::ParticleAttribute &attrX0 = particleData.addAttribute("position0", Partio::VECTOR, 3);
const Partio::ParticleAttribute &attrX = particleData.addAttribute("position", Partio::VECTOR, 3);
const Partio::ParticleAttribute &attrVel = particleData.addAttribute("velocity", Partio::VECTOR, 3);
const Partio::ParticleAttribute &attrVol = particleData.addAttribute("volume", Partio::FLOAT, 1);
const unsigned int numParticles = model->numberOfParticles();
for (unsigned int i = 0; i < numParticles; i )
{
Partio::ParticleIndex index = particleData.addParticle();
float* val = particleData.dataWrite<float>(attrX0, index);
const Vector3r &x0 = model->getPosition0(i);
val[0] = (float)x0[0];
val[1] = (float)x0[1];
val[2] = (float)x0[2];
val = particleData.dataWrite<float>(attrX, index);
const Vector3r &x = model->getPosition(i);
val[0] = (float)x[0];
val[1] = (float)x[1];
val[2] = (float)x[2];
val = particleData.dataWrite<float>(attrVel, index);
const Vector3r &v = model->getVelocity(i);
val[0] = (float)v[0];
val[1] = (float)v[1];
val[2] = (float)v[2];
val = particleData.dataWrite<float>(attrVol, index);
val[0] = (float)model->getVolume(i);
}
Partio::write(fileName.c_str(), particleData, true);
particleData.release();
}
}
void SimulatorBase::readBoundaryState(const std::string &fileName, BoundaryModel *bm)
{
Simulation *sim = Simulation::getCurrent();
if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
{
if (!FileSystem::fileExists(fileName))
{
LOG_WARN << "File " << fileName << " does not exist.";
return;
}
BoundaryModel_Akinci2012 *model = static_cast<BoundaryModel_Akinci2012*>(bm);
Partio::ParticlesDataMutable* data = Partio::read(fileName.c_str());
if (!data)
{
LOG_WARN << "Partio file " << fileName << " not readable.";
return;
}
unsigned int pos0Index = 0xffffffff;
unsigned int posIndex = 0xffffffff;
unsigned int velIndex = 0xffffffff;
unsigned int volIndex = 0xffffffff;
for (int i = 0; i < data->numAttributes(); i )
{
Partio::ParticleAttribute attr;
data->attributeInfo(i, attr);
if (attr.name == "position0")
pos0Index = i;
else if (attr.name == "position")
posIndex = i;
else if (attr.name == "velocity")
velIndex = i;
else if (attr.name == "volume")
volIndex = i;
}
if ((pos0Index == 0xffffffff) ||
(posIndex == 0xffffffff) ||
(velIndex == 0xffffffff) ||
(volIndex == 0xffffffff))
{
LOG_WARN << "File " << fileName << " does not has the correct attributes.";
return;
}
Partio::ParticleAttribute attrX0;
Partio::ParticleAttribute attrX;
Partio::ParticleAttribute attrVel;
Partio::ParticleAttribute attrVol;
data->attributeInfo(pos0Index, attrX0);
data->attributeInfo(posIndex, attrX);
data->attributeInfo(velIndex, attrVel);
data->attributeInfo(volIndex, attrVol);
model->resize(data->numParticles());
for (int i = 0; i < data->numParticles(); i )
{
const float *pos0 = data->data<float>(attrX0, i);
model->getPosition0(i) = Vector3r(pos0[0], pos0[1], pos0[2]);
const float *pos = data->data<float>(attrX, i);
model->getPosition(i) = Vector3r(pos[0], pos[1], pos[2]);
const float *vel = data->data<float>(attrVel, i);
model->getVelocity(i) = Vector3r(vel[0], vel[1], vel[2]);
const float *vol = data->data<float>(attrVol, i);
model->setVolume(i, vol[0]);
}
data->release();
NeighborhoodSearch *neighborhoodSearch = Simulation::getCurrent()->getNeighborhoodSearch();
neighborhoodSearch->update_point_sets();
neighborhoodSearch->resize_point_set(model->getPointSetIndex(), &model->getPosition(0)[0], model->numberOfParticles());
}
}
void SimulatorBase::initDensityMap(std::vector<Vector3r> &x, std::vector<unsigned int> &faces, const Utilities::SceneLoader::BoundaryData *boundaryData, const bool md5, const bool isDynamic, BoundaryModel_Koschier2017 *boundaryModel)
{
Simulation *sim = Simulation::getCurrent();
const Real supportRadius = sim->getSupportRadius();
std::string scene_path = FileSystem::getFilePath(getSceneFile());
std::string scene_file_name = FileSystem::getFileName(getSceneFile());
SceneLoader::Scene &scene = getScene();
const bool useCache = getUseParticleCaching();
Discregrid::CubicLagrangeDiscreteGrid *densityMap;
// if a map file is given, use this one
if (boundaryData->mapFile != "")
{
std::string mapFileName = boundaryData->mapFile;
if (FileSystem::isRelativePath(mapFileName))
mapFileName = FileSystem::normalizePath(scene_path "/" mapFileName);
densityMap = new Discregrid::CubicLagrangeDiscreteGrid(mapFileName);
boundaryModel->setMap(densityMap);
LOG_INFO << "Loaded density map: " << mapFileName;
return;
}
string cachePath = scene_path "/Cache";
// Cache map
std::string mesh_base_path = FileSystem::getFilePath(boundaryData->meshFile);
std::string mesh_file_name = FileSystem::getFileName(boundaryData->meshFile);
Eigen::Matrix<unsigned int, 3, 1> resolutionSDF = boundaryData->mapResolution;
const string scaleStr = "s" real2String(boundaryData->scale[0]) "_" real2String(boundaryData->scale[1]) "_" real2String(boundaryData->scale[2]);
const string resStr = "r" to_string(resolutionSDF[0]) "_" to_string(resolutionSDF[1]) "_" to_string(resolutionSDF[2]);
const string invertStr = "i" to_string((int)boundaryData->mapInvert);
const string thicknessStr = "t" real2String(boundaryData->mapThickness);
const string kernelStr = "k" to_string(sim->getKernel());
string densityMapFileName = "";
if (isDynamic)
densityMapFileName = FileSystem::normalizePath(cachePath "/" mesh_file_name "_db_dm_" real2String(scene.particleRadius) "_" scaleStr "_" resStr "_" invertStr "_" thicknessStr "_" kernelStr ".cdm");
else
densityMapFileName = FileSystem::normalizePath(cachePath "/" mesh_file_name "_sb_dm_" real2String(scene.particleRadius) "_" scaleStr "_" resStr "_" invertStr "_" thicknessStr "_" kernelStr ".cdm");
// check MD5 if cache file is available
bool foundCacheFile = false;
if (useCache)
foundCacheFile = FileSystem::fileExists(densityMapFileName);
if (useCache && foundCacheFile && md5)
{
densityMap = new Discregrid::CubicLagrangeDiscreteGrid(densityMapFileName);
boundaryModel->setMap(densityMap);
LOG_INFO << "Loaded cached density map: " << densityMapFileName;
return;
}
if (!useCache || !foundCacheFile || !md5)
{
//////////////////////////////////////////////////////////////////////////
// Generate distance field of object using Discregrid
//////////////////////////////////////////////////////////////////////////
#ifdef USE_DOUBLE
Discregrid::TriangleMesh sdfMesh(&x[0][0], faces.data(), x.size(), faces.size() / 3);
#else
// if type is float, copy vector to double vector
std::vector<double> doubleVec;
doubleVec.resize(3 * x.size());
for (unsigned int i = 0; i < x.size(); i )
for (unsigned int j = 0; j < 3; j )
doubleVec[3 * i j] = x[i][j];
Discregrid::TriangleMesh sdfMesh(&doubleVec[0], faces.data(), x.size(), faces.size() / 3);
#endif
Discregrid::MeshDistance md(sdfMesh);
Eigen::AlignedBox3d domain;
for (auto const& x_ : x)
{
domain.extend(x_.cast<double>());
}
const Real tolerance = boundaryData->mapThickness;
domain.max() = (4.0*supportRadius tolerance) * Eigen::Vector3d::Ones();
domain.min() -= (4.0*supportRadius tolerance) * Eigen::Vector3d::Ones();
LOG_INFO << "Domain - min: " << domain.min()[0] << ", " << domain.min()[1] << ", " << domain.min()[2];
LOG_INFO << "Domain - max: " << domain.max()[0] << ", " << domain.max()[1] << ", " << domain.max()[2];
LOG_INFO << "Set SDF resolution: " << resolutionSDF[0] << ", " << resolutionSDF[1] << ", " << resolutionSDF[2];
densityMap = new Discregrid::CubicLagrangeDiscreteGrid(domain, std::array<unsigned int, 3>({ resolutionSDF[0], resolutionSDF[1], resolutionSDF[2] }));
auto func = Discregrid::DiscreteGrid::ContinuousFunction{};
Real sign = 1.0;
if (boundaryData->mapInvert)
sign = -1.0;
func = [&md, &sign, &tolerance](Eigen::Vector3d const& xi) {return sign * (md.signedDistanceCached(xi) - tolerance); };
LOG_INFO << "Generate SDF";
START_TIMING("SDF Construction");
densityMap->addFunction(func, false);
STOP_TIMING_AVG
const bool sim2D = sim->is2DSimulation();
//////////////////////////////////////////////////////////////////////////
// Generate density map of object using Discregrid
//////////////////////////////////////////////////////////////////////////
if (sim2D)
SimpleQuadrature::determineSamplePointsInCircle(supportRadius, 30);
auto int_domain = Eigen::AlignedBox3d(Eigen::Vector3d::Constant(-supportRadius), Eigen::Vector3d::Constant(supportRadius));
Real factor = 5.0;
if (sim2D)
factor = 1.75;
auto density_func = [&](Eigen::Vector3d const& x)
{
auto d = densityMap->interpolate(0u, x);
if (d > (1.0 1.0 / factor) * supportRadius)
{
return 0.0;
}
auto integrand = [&](Eigen::Vector3d const& xi)
{
if (xi.squaredNorm() > supportRadius*supportRadius)
return 0.0;
auto dist = densityMap->interpolate(0u, x xi);
// Linear function gamma
if (dist > 1.0 / factor * supportRadius)
return 0.0;
return static_cast<double>((1.0 - factor * dist / supportRadius) * sim->W(xi.cast<Real>()));
};
double res = 0.0;
if (sim2D)
res = 0.8 * SimpleQuadrature::integrate(integrand);
else
res = 0.8 * GaussQuadrature::integrate(integrand, int_domain, 50);
return res;
};
auto cell_diag = densityMap->cellSize().norm();
std::cout << "Generate density map..." << std::endl;
const bool no_reduction = true;
START_TIMING("Density Map Construction");
densityMap->addFunction(density_func, false, [&](Eigen::Vector3d const& x_)
{
if (no_reduction)
{
return true;
}
auto x = x_.cwiseMax(densityMap->domain().min()).cwiseMin(densityMap->domain().max());
auto dist = densityMap->interpolate(0u, x);
if (dist == std::numeric_limits<double>::max())
{
return false;
}
return fabs(dist) < 2.5 * supportRadius;
});
STOP_TIMING_PRINT;
// reduction
if (!no_reduction)
{
std::cout << "Reduce discrete fields...";
densityMap->reduceField(0u, [&](const Eigen::Vector3d &, double v)
{
return fabs(v) < 2.5 * supportRadius;
});
densityMap->reduceField(1u, [&](const Eigen::Vector3d &, double v)->double
{
if (v == std::numeric_limits<double>::max())
return false;
return true;
});
std::cout << "DONE" << std::endl;
}
boundaryModel->setMap(densityMap);
// Store cache file
if (useCache && (FileSystem::makeDir(cachePath) == 0))
{
LOG_INFO << "Save density map: " << densityMapFileName;
densityMap->save(densityMapFileName);
}
}
}
void SimulatorBase::initVolumeMap(std::vector<Vector3r> &x, std::vector<unsigned int> &faces, const Utilities::SceneLoader::BoundaryData *boundaryData, const bool md5, const bool isDynamic, BoundaryModel_Bender2019 *boundaryModel)
{
Simulation *sim = Simulation::getCurrent();
const Real supportRadius = sim->getSupportRadius();
std::string scene_path = FileSystem::getFilePath(getSceneFile());
std::string scene_file_name = FileSystem::getFileName(getSceneFile());
SceneLoader::Scene &scene = getScene();
const bool useCache = getUseParticleCaching();
Discregrid::CubicLagrangeDiscreteGrid *volumeMap;
// if a map file is given, use this one
if (boundaryData->mapFile != "")
{
std::string mapFileName = boundaryData->mapFile;
if (FileSystem::isRelativePath(mapFileName))
mapFileName = FileSystem::normalizePath(scene_path "/" mapFileName);
volumeMap = new Discregrid::CubicLagrangeDiscreteGrid(mapFileName);
boundaryModel->setMap(volumeMap);
LOG_INFO << "Loaded volume map: " << mapFileName;
return;
}
string cachePath = scene_path "/Cache";
// Cache map
std::string mesh_base_path = FileSystem::getFilePath(boundaryData->meshFile);
std::string mesh_file_name = FileSystem::getFileName(boundaryData->meshFile);
Eigen::Matrix<unsigned int, 3, 1> resolutionSDF = boundaryData->mapResolution;
const string scaleStr = "s" real2String(boundaryData->scale[0]) "_" real2String(boundaryData->scale[1]) "_" real2String(boundaryData->scale[2]);
const string resStr = "r" to_string(resolutionSDF[0]) "_" to_string(resolutionSDF[1]) "_" to_string(resolutionSDF[2]);
const string invertStr = "i" to_string((int)boundaryData->mapInvert);
const string thicknessStr = "t" real2String(boundaryData->mapThickness);
string volumeMapFileName = "";
if (isDynamic)
volumeMapFileName = FileSystem::normalizePath(cachePath "/" mesh_file_name "_db_vm_" real2String(scene.particleRadius) "_" scaleStr "_" resStr "_" invertStr "_" thicknessStr ".cdm");
else
volumeMapFileName = FileSystem::normalizePath(cachePath "/" mesh_file_name "_sb_vm_" real2String(scene.particleRadius) "_" scaleStr "_" resStr "_" invertStr "_" thicknessStr ".cdm");
// check MD5 if cache file is available
bool foundCacheFile = false;
if (useCache)
foundCacheFile = FileSystem::fileExists(volumeMapFileName);
if (useCache && foundCacheFile && md5)
{
volumeMap = new Discregrid::CubicLagrangeDiscreteGrid(volumeMapFileName);
boundaryModel->setMap(volumeMap);
LOG_INFO << "Loaded cached volume map: " << volumeMapFileName;
return;
}
if (!useCache || !foundCacheFile || !md5)
{
//////////////////////////////////////////////////////////////////////////
// Generate distance field of object using Discregrid
//////////////////////////////////////////////////////////////////////////
#ifdef USE_DOUBLE
Discregrid::TriangleMesh sdfMesh(&x[0][0], faces.data(), x.size(), faces.size() / 3);
#else
// if type is float, copy vector to double vector
std::vector<double> doubleVec;
doubleVec.resize(3 * x.size());
for (unsigned int i = 0; i < x.size(); i )
for (unsigned int j = 0; j < 3; j )
doubleVec[3 * i j] = x[i][j];
Discregrid::TriangleMesh sdfMesh(&doubleVec[0], faces.data(), x.size(), faces.size() / 3);
#endif
Discregrid::MeshDistance md(sdfMesh);
Eigen::AlignedBox3d domain;
for (auto const& x_ : x)
{
domain.extend(x_.cast<double>());
}
const Real tolerance = boundaryData->mapThickness;
domain.max() = (4.0*supportRadius tolerance) * Eigen::Vector3d::Ones();
domain.min() -= (4.0*supportRadius tolerance) * Eigen::Vector3d::Ones();
LOG_INFO << "Domain - min: " << domain.min()[0] << ", " << domain.min()[1] << ", " << domain.min()[2];
LOG_INFO << "Domain - max: " << domain.max()[0] << ", " << domain.max()[1] << ", " << domain.max()[2];
LOG_INFO << "Set SDF resolution: " << resolutionSDF[0] << ", " << resolutionSDF[1] << ", " << resolutionSDF[2];
volumeMap = new Discregrid::CubicLagrangeDiscreteGrid(domain, std::array<unsigned int, 3>({ resolutionSDF[0], resolutionSDF[1], resolutionSDF[2] }));
auto func = Discregrid::DiscreteGrid::ContinuousFunction{};
//volumeMap->setErrorTolerance(0.001);
Real sign = 1.0;
if (boundaryData->mapInvert)
sign = -1.0;
const Real particleRadius = sim->getParticleRadius();
// subtract 0.5 * particle radius to prevent penetration of particles and the boundary
func = [&md, &sign, &tolerance, &particleRadius](Eigen::Vector3d const& xi) {return sign * (md.signedDistanceCached(xi) - tolerance - 0.5*particleRadius); };
LOG_INFO << "Generate SDF";
START_TIMING("SDF Construction");
volumeMap->addFunction(func, false);
STOP_TIMING_PRINT
//////////////////////////////////////////////////////////////////////////
// Generate volume map of object using Discregrid
//////////////////////////////////////////////////////////////////////////
Simulation *sim = Simulation::getCurrent();
const bool sim2D = sim->is2DSimulation();
if (sim2D)
SimpleQuadrature::determineSamplePointsInCircle(supportRadius, 30);
auto int_domain = Eigen::AlignedBox3d(Eigen::Vector3d::Constant(-supportRadius), Eigen::Vector3d::Constant(supportRadius));
Real factor = 1.0;
if (sim2D)
factor = 1.75;
auto volume_func = [&](Eigen::Vector3d const& x)
{
auto dist = volumeMap->interpolate(0u, x);
if (dist > (1.0 1.0 /*/ factor*/) * supportRadius)
{
return 0.0;
}
auto integrand = [&volumeMap, &x, &supportRadius, &factor, &sim, &sim2D](Eigen::Vector3d const& xi) -> double
{
if (xi.squaredNorm() > supportRadius*supportRadius)
return 0.0;
auto dist = volumeMap->interpolate(0u, x xi);
if (dist <= 0.0)
return 1.0 - 0.1 * dist / supportRadius;
if (dist < 1.0 / factor * supportRadius)
return static_cast<double>(CubicKernel::W(factor * static_cast<Real>(dist)) / CubicKernel::W_zero());
return 0.0;
};
double res = 0.0;
if (sim2D)
res = 0.8 * SimpleQuadrature::integrate(integrand);
else
res = 0.8 * GaussQuadrature::integrate(integrand, int_domain, 30);
return res;
};
auto cell_diag = volumeMap->cellSize().norm();
std::cout << "Generate volume map..." << std::endl;
const bool no_reduction = true;
START_TIMING("Volume Map Construction");
volumeMap->addFunction(volume_func, false, [&](Eigen::Vector3d const& x_)
{
if (no_reduction)
{
return true;
}
auto x = x_.cwiseMax(volumeMap->domain().min()).cwiseMin(volumeMap->domain().max());
auto dist = volumeMap->interpolate(0u, x);
if (dist == std::numeric_limits<double>::max())
{
return false;
}
return fabs(dist) < 4.0 * supportRadius;
});
STOP_TIMING_PRINT;
// reduction
if (!no_reduction)
{
std::cout << "Reduce discrete fields...";
volumeMap->reduceField(0u, [&](const Eigen::Vector3d &, double v)
{
return fabs(v) < 4.0 * supportRadius;
});
volumeMap->reduceField(1u, [&](const Eigen::Vector3d &, double v)->double
{
if (v == std::numeric_limits<double>::max())
return false;
return true;
});
std::cout << "DONE" << std::endl;
}
boundaryModel->setMap(volumeMap);
// Store cache file
if (useCache && (FileSystem::makeDir(cachePath) == 0))
{
LOG_INFO << "Save volume map: " << volumeMapFileName;
volumeMap->save(volumeMapFileName);
}
}
// store maximal distance of a point to center of mass for CFL
if (boundaryData->dynamic)
{
// determine center of mass
Vector3r com;
com.setZero();
for (unsigned int i = 0; i < x.size(); i )
{
com = x[i];
}
com /= static_cast<Real>(x.size());
// determine point with maximal distance to center of mass
Real maxDist = 0.0;
for (unsigned int i = 0; i < x.size(); i )
{
const Vector3r diff = x[i] - com;
const Real dist = diff.norm();
if (dist > maxDist)
{
maxDist = dist;
}
}
boundaryModel->setMaxDist(maxDist);
}
}
