diff --git a/Examples/ComputeGraph/Benchmarking/CMakeLists.txt b/Examples/ComputeGraph/Benchmarking/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c52359351e19f0e6d0986cf3181aa22b198ae7f4
--- /dev/null
+++ b/Examples/ComputeGraph/Benchmarking/CMakeLists.txt
@@ -0,0 +1,34 @@
+###########################################################################
+#
+# Copyright (c) Kitware, Inc.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+###########################################################################
+
+project(Example-ComputeGraphBenchmark)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME} computeGraphBenchmarkExample.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples/ComputeGraph)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} SimulationManager apiUtilities)
diff --git a/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp b/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8aeeca5ed7163839a2887c5e84034b93416a3f79
--- /dev/null
+++ b/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp
@@ -0,0 +1,159 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#include "imstkAPIUtilities.h"
+#include "imstkCamera.h"
+#include "imstkComputeGraphVizWriter.h"
+#include "imstkLineMesh.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkScene.h"
+#include "imstkSceneManager.h"
+#include "imstkSimulationManager.h"
+
+using namespace imstk;
+
+///
+/// \brief This examples uses the benchmarking features of the computational graph
+/// This allows one to see the elapsed time of every step of iMSTK as well as
+/// export the computational graph and show information such as the critical path
+///
+int
+main()
+{
+ auto simManager = std::make_shared<SimulationManager>();
+ auto scene = simManager->createNewScene("PBDString");
+ scene->getConfig()->benchmarkingEnabled = true;
+
+ // Setup N separate string simulations with varying bend stiffnesses
+ const unsigned int numStrings = 8;
+ const unsigned int numVerts = 30;
+ const double stringSpacing = 2.0; // How far each string is apart
+ const double stringLength = 10.0; // Total length of string
+ const Color startColor = Color::Red; // Color of first string
+ const Color endColor = Color::Green; // Color of last string
+
+ struct PbdSim
+ {
+ std::shared_ptr<LineMesh> geometry;
+ std::shared_ptr<PbdObject> object;
+ std::shared_ptr<PbdModel> model;
+ std::shared_ptr<PBDModelConfig> params;
+ std::shared_ptr<VisualModel> visualModel;
+ };
+ std::vector<PbdSim> sims(numStrings);
+
+ const double size = stringSpacing * (numStrings - 1);
+ const double vertexSpacing = stringLength / numVerts;
+ for (unsigned int i = 0; i < numStrings; i++)
+ {
+ // Setup the line mesh
+ sims[i].geometry = std::make_shared<LineMesh>();
+ StdVectorOfVec3d vertList;
+ vertList.resize(numVerts);
+ for (size_t j = 0; j < numVerts; j++)
+ {
+ vertList[j] = Vec3d(
+ static_cast<double>(i) * stringSpacing - size * 0.5,
+ stringLength * 0.5 - static_cast<double>(j) * vertexSpacing, 0.0);
+ }
+ sims[i].geometry->setInitialVertexPositions(vertList);
+ sims[i].geometry->setVertexPositions(vertList);
+
+ // Add connectivity data
+ std::vector<LineMesh::LineArray> segments;
+ for (size_t j = 0; j < numVerts - 1; j++)
+ {
+ LineMesh::LineArray seg = { j, j + 1 };
+ segments.push_back(seg);
+ }
+
+ sims[i].geometry->setLinesVertices(segments);
+
+ sims[i].object = std::make_shared<PbdObject>("String " + std::to_string(i));
+ sims[i].model = std::make_shared<PbdModel>();
+ sims[i].model->setModelGeometry(sims[i].geometry);
+
+ // Configure the parameters with bend stiffnesses varying from 0.001 to ~0.1
+ sims[i].params = std::make_shared<PBDModelConfig>();
+ sims[i].params->enableConstraint(PbdConstraint::Type::Distance, 0.001);
+ sims[i].params->enableConstraint(PbdConstraint::Type::Bend, static_cast<double>(i) * 0.1 / numStrings + 0.001);
+ sims[i].params->m_fixedNodeIds = { 0 };
+ sims[i].params->m_uniformMassValue = 5.0;
+ sims[i].params->m_gravity = Vec3d(0, -9.8, 0);
+ sims[i].params->m_defaultDt = 0.0005;
+ sims[i].params->m_iterations = 5;
+
+ // Set the parameters
+ sims[i].model->configure(sims[i].params);
+ sims[i].object->setDynamicalModel(sims[i].model);
+ sims[i].object->setPhysicsGeometry(sims[i].geometry);
+
+ sims[i].visualModel = std::make_shared<VisualModel>(sims[i].geometry);
+ std::shared_ptr<RenderMaterial> material = std::make_shared<RenderMaterial>();
+ material->setDisplayMode(RenderMaterial::DisplayMode::Wireframe);
+ material->setDebugColor(Color::lerpRgb(startColor, endColor, static_cast<double>(i) / (numStrings - 1)));
+ material->setLineWidth(2.0f);
+ sims[i].visualModel->setRenderMaterial(material);
+ sims[i].object->addVisualModel(sims[i].visualModel);
+
+ // Add in scene
+ scene->addSceneObject(sims[i].object);
+ }
+
+ // Adjust the camera
+ scene->getCamera()->setFocalPoint(0.0, 0.0, 0.0);
+ scene->getCamera()->setPosition(0.0, 0.0, 15.0);
+
+ // Move the points every frame
+ double t = 0.0;
+ const double dt = 0.0005;
+ const double radius = 1.5;
+ auto movePoints =
+ [&sims, &t, dt, radius](Module* module)
+ {
+ for (unsigned int i = 0; i < sims.size(); i++)
+ {
+ Vec3d pos = sims[i].model->getCurrentState()->getVertexPosition(0);
+ // Move in circle, derivatives of parametric eq of circle
+ sims[i].model->getCurrentState()->setVertexPosition(0, imstk::Vec3d(
+ pos.x() + -std::sin(t) * radius * dt,
+ pos.y(),
+ pos.z() + std::cos(t) * radius * dt));
+ }
+ t += dt;
+ };
+ simManager->getSceneManager(scene)->setPostUpdateCallback(movePoints);
+
+ // Start
+ simManager->setActiveScene(scene);
+ simManager->start();
+
+ // Write the graph, highlighting the critical path and putting the completion time in the name
+ ComputeGraphVizWriter writer;
+ writer.setInput(scene->getComputeGraph());
+ writer.setFileName("computeGraphBenchmarkExample.svg");
+ writer.setHighlightCriticalPath(true);
+ writer.setWriteTimes(true);
+ writer.write();
+
+ return 0;
+}
diff --git a/Examples/ComputeGraph/CMakeLists.txt b/Examples/ComputeGraph/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..964a9fe46555e89800ccc1ca30875fc8ac896dab
--- /dev/null
+++ b/Examples/ComputeGraph/CMakeLists.txt
@@ -0,0 +1,21 @@
+#-----------------------------------------------------------------------------
+# Add Example subdirectories
+#-----------------------------------------------------------------------------
+macro(listOfSubDir result curdir)
+ file(GLOB children RELATIVE ${curdir} ${curdir}/*)
+ set(dirlist "")
+ foreach(child ${children})
+ if(IS_DIRECTORY ${curdir}/${child})
+ list(APPEND dirlist ${child})
+ endif()
+ endforeach()
+ set(${result} ${dirlist})
+endmacro()
+
+listOfSubDir(subDirs ${CMAKE_CURRENT_SOURCE_DIR})
+
+foreach(subdir ${subDirs})
+ add_subdirectory(${subdir})
+endforeach()
+
+
diff --git a/Examples/ComputeGraph/Configuration/CMakeLists.txt b/Examples/ComputeGraph/Configuration/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ecf6a10ca4526d0f2fdff3122cefa3017f16692
--- /dev/null
+++ b/Examples/ComputeGraph/Configuration/CMakeLists.txt
@@ -0,0 +1,34 @@
+###########################################################################
+#
+# Copyright (c) Kitware, Inc.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+###########################################################################
+
+project(Example-ComputeGraphConfigure)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME} computeGraphConfigureExample.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples/ComputeGraph)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} SimulationManager apiUtilities)
diff --git a/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp b/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..65e03f79d75b152d024e615ba5556b7b22611aad
--- /dev/null
+++ b/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp
@@ -0,0 +1,212 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+
+#include "imstkAPIUtilities.h"
+#include "imstkCamera.h"
+#include "imstkComputeGraph.h"
+#include "imstkComputeGraphVizWriter.h"
+#include "imstkLight.h"
+#include "imstkParallelReduce.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkScene.h"
+#include "imstkSceneManager.h"
+#include "imstkSimulationManager.h"
+#include "imstkSurfaceMesh.h"
+
+#include "imstkVTKRenderer.h"
+#include "VTKRenderer.h"
+#include "vtkRenderWindowInteractor.h"
+
+using namespace imstk;
+
+static std::unique_ptr<SurfaceMesh> makeCloth(
+ const double width, const double height, const int nRows, const int nCols)
+{
+ // Create surface mesh
+ std::unique_ptr<SurfaceMesh> clothMesh = std::make_unique<SurfaceMesh>();
+ StdVectorOfVec3d vertList;
+
+ vertList.resize(nRows * nCols);
+ const double dy = width / (double)(nCols - 1);
+ const double dx = height / (double)(nRows - 1);
+ for (int i = 0; i < nRows; ++i)
+ {
+ for (int j = 0; j < nCols; j++)
+ {
+ vertList[i * nCols + j] = Vec3d((double)dx * i, 1.0, (double)dy * j);
+ }
+ }
+ clothMesh->setInitialVertexPositions(vertList);
+ clothMesh->setVertexPositions(vertList);
+
+ // Add connectivity data
+ std::vector<SurfaceMesh::TriangleArray> triangles;
+ for (std::size_t i = 0; i < nRows - 1; ++i)
+ {
+ for (std::size_t j = 0; j < nCols - 1; j++)
+ {
+ SurfaceMesh::TriangleArray tri[2];
+ const size_t index1 = i * nCols + j;
+ const size_t index2 = index1 + nCols;
+ const size_t index3 = index1 + 1;
+ const size_t index4 = index2 + 1;
+
+ // Interleave [/][\]
+ if (i % 2 ^ j % 2)
+ {
+ tri[0] = { { index1, index2, index3 } };
+ tri[1] = { { index4, index3, index2 } };
+ }
+ else
+ {
+ tri[0] = { { index2, index4, index1 } };
+ tri[1] = { { index4, index3, index1 } };
+ }
+ triangles.push_back(tri[0]);
+ triangles.push_back(tri[1]);
+ }
+ }
+
+ clothMesh->setTrianglesVertices(triangles);
+
+ return clothMesh;
+}
+
+static std::shared_ptr<PbdObject> makeClothObj(const std::string& name, double width, double height, int nRows, int nCols)
+{
+ auto clothObj = std::make_shared<PbdObject>(name);
+
+ std::shared_ptr<SurfaceMesh> clothMesh(std::move(makeCloth(width, height, nRows, nCols)));
+
+ // Setup the Parameters
+ auto pbdParams = std::make_shared<PBDModelConfig>();
+ pbdParams->m_uniformMassValue = 1.0;
+ pbdParams->m_gravity = Vec3d(0, -9.8, 0);
+ pbdParams->m_defaultDt = 0.005;
+ pbdParams->m_iterations = 5;
+ pbdParams->m_viscousDampingCoeff = 0.05;
+ pbdParams->enableConstraint(PbdConstraint::Type::Distance, 0.1);
+ pbdParams->enableConstraint(PbdConstraint::Type::Dihedral, 0.001);
+ std::vector<size_t> fixedNodes(nCols);
+ for (size_t i = 0; i < fixedNodes.size(); i++)
+ {
+ fixedNodes[i] = i;
+ }
+ pbdParams->m_fixedNodeIds = fixedNodes;
+
+ // Setup the Model
+ auto pbdModel = std::make_shared<PbdModel>();
+ pbdModel->setModelGeometry(clothMesh);
+ pbdModel->configure(pbdParams);
+
+ // Setup the VisualModel
+ auto material = std::make_shared<RenderMaterial>();
+ material->setBackFaceCulling(false);
+ material->setColor(Color::LightGray);
+ material->setDisplayMode(RenderMaterial::DisplayMode::WireframeSurface);
+
+ auto clothVisualModel = std::make_shared<VisualModel>(clothMesh);
+ clothVisualModel->setRenderMaterial(material);
+
+ // Setup the Object
+ clothObj->addVisualModel(clothVisualModel);
+ clothObj->setPhysicsGeometry(clothMesh);
+ clothObj->setDynamicalModel(pbdModel);
+
+ return clothObj;
+}
+
+///
+/// \brief This example demonstrates how to modify the computational graph, post initialization
+///
+int
+main()
+{
+ auto simManager = std::make_shared<SimulationManager>();
+ auto scene = simManager->createNewScene("PBDCloth");
+
+ const double width = 10.0;
+ const double height = 10.0;
+ const int nRows = 16;
+ const int nCols = 16;
+ std::shared_ptr<PbdObject> clothObj = makeClothObj("Cloth", width, height, nRows, nCols);
+ scene->addSceneObject(clothObj);
+
+ // Light (white)
+ auto whiteLight = std::make_shared<DirectionalLight>("whiteLight");
+ whiteLight->setFocalPoint(Vec3d(5, -8, -5));
+ whiteLight->setIntensity(7);
+ scene->addLight(whiteLight);
+
+ // Light (red)
+ auto colorLight = std::make_shared<SpotLight>("colorLight");
+ colorLight->setPosition(Vec3d(-5, -3, 5));
+ colorLight->setFocalPoint(Vec3d(0, -5, 5));
+ colorLight->setIntensity(100);
+ colorLight->setColor(Color::Red);
+ colorLight->setSpotAngle(30);
+ scene->addLight(colorLight);
+
+ // Adjust camera
+ scene->getCamera()->setFocalPoint(0, -5, 5);
+ scene->getCamera()->setPosition(-15., -5.0, 15.0);
+
+ // Adds a custom physics step to print out intermediate velocities
+ {
+ std::shared_ptr<PbdModel> pbdModel = clothObj->getPbdModel();
+ scene->setComputeGraphConfigureCallback([&](Scene* scene)
+ {
+ // Get the graph
+ std::shared_ptr<ComputeGraph> graph = scene->getComputeGraph();
+
+ // First write the graph before we make modifications, just to show the changes
+ ComputeGraphVizWriter writer;
+ writer.setInput(graph);
+ writer.setFileName("computeGraphConfigureExampleOld.svg");
+ writer.write();
+
+ std::shared_ptr<ComputeNode> printVelocities = std::make_shared<ComputeNode>([&]()
+ {
+ // Make the timestep a function of max velocity
+ const StdVectorOfVec3d& velocities = *pbdModel->getCurrentState()->getVelocities();
+ for (size_t i = 0; i < velocities.size(); i++)
+ {
+ printf("Velocity: %f, %f, %f\n", velocities[i].x(), velocities[i].y(), velocities[i].z());
+ }
+
+ }, "PrintVelocities");
+
+ // After IntegratePosition
+ graph->insertAfter(pbdModel->getIntegratePositionNode(), printVelocities);
+
+ // Write the modified graph
+ writer.setFileName("computeGraphConfigureExampleNew.svg");
+ writer.write();
+ });
+ }
+
+ // Start
+ simManager->setActiveScene(scene);
+ simManager->start(SimulationStatus::Paused);
+
+ return 0;
+}
diff --git a/Source/Common/ComputeGraph/imstkComputeGraph.cpp b/Source/Common/ComputeGraph/imstkComputeGraph.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8957da00501ccc1401be678a17830ac7ed4e2816
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeGraph.cpp
@@ -0,0 +1,830 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#include "imstkComputeGraph.h"
+#include "imstkComputeNode.h"
+#include "imstkLogger.h"
+#include <queue>
+#include <stack>
+#include <set>
+
+namespace imstk
+{
+ComputeGraph::ComputeGraph(std::string sourceName, std::string sinkName) :
+ m_source(std::make_shared<ComputeNode>()),
+ m_sink(std::make_shared<ComputeNode>())
+{
+ m_source->m_name = sourceName;
+ m_sink->m_name = sinkName;
+ addNode(m_source);
+ addNode(m_sink);
+}
+
+ComputeNodeVector::iterator
+ComputeGraph::findNode(std::string name)
+{
+ return std::find_if(m_nodes.begin(), m_nodes.end(),
+ [&name](const std::shared_ptr<ComputeNode>& x) { return x->m_name == name; });
+}
+
+ComputeNodeVector::const_iterator
+ComputeGraph::findNode(std::string name) const
+{
+ return std::find_if(m_nodes.begin(), m_nodes.end(),
+ [&name](const std::shared_ptr<ComputeNode>& x) { return x->m_name == name; });
+}
+
+ComputeNodeVector::iterator
+ComputeGraph::findNode(std::shared_ptr<ComputeNode> node)
+{
+ return std::find(m_nodes.begin(), m_nodes.end(), node);
+}
+
+ComputeNodeVector::const_iterator
+ComputeGraph::findNode(std::shared_ptr<ComputeNode> node) const
+{
+ return std::find(m_nodes.begin(), m_nodes.end(), node);
+}
+
+bool
+ComputeGraph::containsNode(std::shared_ptr<ComputeNode> node) const
+{
+ return findNode(node) != nodesEnd();
+}
+
+bool
+ComputeGraph::containsEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode)
+{
+ return (m_adjList.count(srcNode) != 0 && m_adjList[srcNode].count(destNode) != 0);
+}
+
+void
+ComputeGraph::addEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode)
+{
+ m_adjList[srcNode].insert(destNode);
+ m_invAdjList[destNode].insert(srcNode);
+}
+
+void
+ComputeGraph::nestGraph(std::shared_ptr<ComputeGraph> subgraph, std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
+{
+ // Ensure source and sink are in the graph
+ if (findNode(source) == nodesEnd())
+ {
+ LOG(WARNING) << "Tried to nest a graph using source, but source does not exist in this";
+ return;
+ }
+ if (findNode(sink) == nodesEnd())
+ {
+ LOG(WARNING) << "Tried to nest a graph using sink, but sink does not exist in this";
+ return;
+ }
+
+ // Copy all of the nodes into this graph (check duplicates)
+ for (ComputeNodeVector::iterator it = subgraph->getNodes().begin(); it != subgraph->getNodes().end(); it++)
+ {
+ addNode(*it);
+ }
+ // Copy all the edges into the graph (no need to check for duplicates)
+ const ComputeNodeAdjList& adjList = subgraph->getAdjList();
+ for (ComputeNodeAdjList::const_iterator it = adjList.begin(); it != adjList.end(); it++)
+ {
+ const ComputeNodeSet& outputNodes = it->second;
+ for (ComputeNodeSet::const_iterator jt = outputNodes.begin(); jt != outputNodes.end(); jt++)
+ {
+ addEdge(it->first, *jt);
+ }
+ }
+ // Add source and sink edges
+ addEdge(source, subgraph->getSource());
+ addEdge(subgraph->getSink(), sink);
+}
+
+void
+ComputeGraph::removeEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode)
+{
+ if (m_adjList.count(srcNode) == 0)
+ {
+ return;
+ }
+ if (m_adjList[srcNode].count(destNode) == 0)
+ {
+ return;
+ }
+ m_adjList[srcNode].erase(destNode);
+ m_invAdjList[destNode].erase(srcNode);
+ if (m_adjList[srcNode].size() == 0)
+ {
+ m_adjList.erase(srcNode);
+ }
+ if (m_invAdjList[destNode].size() == 0)
+ {
+ m_invAdjList.erase(destNode);
+ }
+}
+
+bool
+ComputeGraph::addNode(std::shared_ptr<ComputeNode> node)
+{
+ // If the node already exists return false
+ if (!containsNode(node))
+ {
+ // Put it in this graph
+ m_nodes.push_back(node);
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+std::shared_ptr<ComputeNode>
+ComputeGraph::addFunction(std::string name, std::function<void()> func)
+{
+ std::shared_ptr<ComputeNode> node = std::make_shared<ComputeNode>(func, name);
+ m_nodes.push_back(node);
+ return node;
+}
+
+bool
+ComputeGraph::removeNode(std::shared_ptr<ComputeNode> node)
+{
+ if (!containsNode(node))
+ {
+ LOG(INFO) << "Tried to remove node " + node->m_name + " from graph but it doesn't contain the node.";
+ return false;
+ }
+
+ // Erase edges
+ ComputeNodeSet inputs = m_invAdjList[node];
+ ComputeNodeSet outputs = m_adjList[node];
+ for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
+ {
+ removeEdge(*i, node);
+ }
+ for (ComputeNodeSet::iterator i = outputs.begin(); i != outputs.end(); i++)
+ {
+ removeEdge(node, *i);
+ }
+
+ // Finally remove the node from the graph
+ ComputeNodeVector::iterator it = findNode(node);
+ if (it != nodesEnd())
+ {
+ m_nodes.erase(it);
+ }
+ return true;
+}
+
+bool
+ComputeGraph::removeNodeAndFix(std::shared_ptr<ComputeNode> node)
+{
+ if (!containsNode(node))
+ {
+ LOG(INFO) << "Tried to remove node " + node->m_name + " from graph but it doesn't contain the node.";
+ return false;
+ }
+
+ // Erase edges
+ ComputeNodeSet inputs = m_invAdjList[node];
+ ComputeNodeSet outputs = m_adjList[node];
+ for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
+ {
+ removeEdge(*i, node);
+ }
+ for (ComputeNodeSet::iterator i = outputs.begin(); i != outputs.end(); i++)
+ {
+ removeEdge(node, *i);
+ }
+
+ // Fix the edges
+ for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
+ {
+ for (ComputeNodeSet::iterator j = outputs.begin(); j != outputs.end(); j++)
+ {
+ addEdge(*i, *j);
+ }
+ }
+
+ // Finally remove the node from the graph
+ ComputeNodeVector::iterator it = findNode(node);
+ if (it != nodesEnd())
+ {
+ m_nodes.erase(it);
+ }
+
+ return true;
+}
+
+void
+ComputeGraph::insertAfter(std::shared_ptr<ComputeNode> refNode, std::shared_ptr<ComputeNode> newNode)
+{
+ // Try to add to graph, if already exists, exit
+ if (!addNode(newNode))
+ {
+ return;
+ }
+
+ // Remove output edges
+ ComputeNodeSet outputs = m_adjList[refNode]; // Copy (since we are modifying)
+ for (ComputeNodeSet::iterator i = outputs.begin(); i != outputs.end(); i++)
+ {
+ removeEdge(refNode, *i);
+ }
+
+ // Add refNode->newNode
+ addEdge(refNode, newNode);
+
+ // Add newNode->outputs[i]
+ for (ComputeNodeSet::iterator i = outputs.begin(); i != outputs.end(); i++)
+ {
+ addEdge(newNode, *i);
+ }
+}
+
+void
+ComputeGraph::insertBefore(std::shared_ptr<ComputeNode> refNode, std::shared_ptr<ComputeNode> newNode)
+{
+ // Try to add to graph, if already exists, exit
+ if (!addNode(newNode))
+ {
+ return;
+ }
+
+ // Remove input edges
+ ComputeNodeSet inputs = m_invAdjList[refNode]; // Copy (since we are modifying)
+ for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
+ {
+ removeEdge(*i, refNode);
+ }
+
+ // inputs[i]->newNode
+ for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
+ {
+ addEdge(*i, newNode);
+ }
+
+ // refNode->newNode
+ addEdge(newNode, refNode);
+}
+
+bool
+ComputeGraph::reaches(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode)
+{
+ const ComputeNodeAdjList& adjList = getAdjList();
+
+ // Simple BFS
+ std::unordered_set<std::shared_ptr<ComputeNode>> visitedNodes;
+
+ // It inserts itself as well
+ std::queue<std::shared_ptr<ComputeNode>> nodeStack;
+ nodeStack.push(srcNode);
+ while (!nodeStack.empty())
+ {
+ std::shared_ptr<ComputeNode> currNode = nodeStack.front();
+ nodeStack.pop();
+
+ // If we've arrived at the desired node
+ if (destNode == currNode)
+ {
+ return true;
+ }
+
+ // Add children to stack if not yet visited
+ if (adjList.count(currNode) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(currNode);
+ for (ComputeNodeSet::const_iterator j = outputNodes.begin(); j != outputNodes.end(); j++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *j;
+ if (visitedNodes.count(childNode) == 0)
+ {
+ visitedNodes.insert(childNode);
+ nodeStack.push(childNode);
+ }
+ }
+ }
+ }
+ return false;
+}
+
+void
+ComputeGraph::clear()
+{
+ m_nodes.clear();
+ clearEdges();
+ addNode(m_source);
+ addNode(m_sink);
+}
+
+std::shared_ptr<ComputeGraph>
+ComputeGraph::sum(std::shared_ptr<ComputeGraph> graphA, std::shared_ptr<ComputeGraph> graphB)
+{
+ std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>();
+ // Get rid of the results source/sink
+ results->m_source = nullptr;
+ results->m_sink = nullptr;
+ results->m_nodes.clear();
+
+ // Sum the nodes
+ ComputeNodeVector& nodesA = graphA->getNodes();
+ for (size_t i = 0; i < nodesA.size(); i++)
+ {
+ results->addNode(nodesA[i]);
+ }
+ ComputeNodeVector& nodesB = graphB->getNodes();
+ for (size_t i = 0; i < nodesB.size(); i++)
+ {
+ results->addNode(nodesB[i]);
+ }
+
+ // Now sum the edges
+ const ComputeNodeAdjList& adjListA = graphA->getAdjList();
+ for (ComputeNodeAdjList::const_iterator it = adjListA.begin(); it != adjListA.end(); it++)
+ {
+ const ComputeNodeSet& outputNodes = it->second;
+ for (ComputeNodeSet::const_iterator jt = outputNodes.begin(); jt != outputNodes.end(); jt++)
+ {
+ results->addEdge(it->first, *jt);
+ }
+ }
+ const ComputeNodeAdjList& adjListB = graphB->getAdjList();
+ for (ComputeNodeAdjList::const_iterator it = adjListB.begin(); it != adjListB.end(); it++)
+ {
+ const ComputeNodeSet& outputNodes = it->second;
+ for (ComputeNodeSet::const_iterator jt = outputNodes.begin(); jt != outputNodes.end(); jt++)
+ {
+ results->addEdge(it->first, *jt);
+ }
+ }
+
+ return results;
+}
+
+std::shared_ptr<ComputeNodeList>
+ComputeGraph::topologicalSort(std::shared_ptr<ComputeGraph> graph)
+{
+ // Compute the number of inputs to each node (we will remove these as we go)
+ std::unordered_map<std::shared_ptr<ComputeNode>, size_t> numInputs;
+
+ const ComputeNodeAdjList& adjList = graph->getAdjList();
+ const ComputeNodeAdjList& invAdjList = graph->getInvAdjList();
+
+ for (ComputeNodeAdjList::const_iterator i = invAdjList.begin(); i != invAdjList.end(); i++)
+ {
+ numInputs[i->first] = invAdjList.size();
+ }
+
+ // Create an edge blacklist for edge removal during algorithm
+ std::unordered_map<std::shared_ptr<ComputeNode>, std::shared_ptr<ComputeNode>> removedEdges;
+
+ auto edgeHasBeenRemoved = [&](const std::shared_ptr<ComputeNode>& node1, const std::shared_ptr<ComputeNode>& node2)
+ {
+ return (removedEdges.count(node1) != 0) && (removedEdges[node1] == node2);
+ };
+
+ // Kahns algorithm (BFS/queue)
+ // iterate through all nodes (BFS or DFS) removing edges
+ // nodes are accepted when all input edges have been removed
+ std::queue<std::shared_ptr<ComputeNode>> sources;
+ sources.push(graph->m_source);
+ std::shared_ptr<ComputeNodeList> results = std::make_shared<ComputeNodeList>();
+
+ while (!sources.empty())
+ {
+ // Remove a node
+ std::shared_ptr<ComputeNode> node = sources.front();
+ sources.pop();
+
+ results->push_back(node);
+
+ // As long as there are children
+ if (adjList.count(node) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(node);
+ for (ComputeNodeSet::const_iterator it = outputNodes.begin(); it != outputNodes.end(); it++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *it;
+ // If edge is present
+ if (!edgeHasBeenRemoved(node, childNode))
+ {
+ // Mark edge as removed
+ removedEdges[node] = childNode;
+
+ // Decrement amount of inputs
+ numInputs[childNode]--;
+ if (numInputs[childNode] <= 0)
+ {
+ sources.push(childNode);
+ }
+ }
+ }
+ }
+ }
+ return results;
+}
+
+std::shared_ptr<ComputeGraph>
+ComputeGraph::resolveCriticalNodes(std::shared_ptr<ComputeGraph> graph)
+{
+ std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>(*graph);
+
+ const ComputeNodeAdjList& adjList = graph->getAdjList();
+ const ComputeNodeVector& nodes = graph->getNodes();
+
+ // Compute the levels of each node via DFS
+ std::unordered_map<std::shared_ptr<ComputeNode>, int> depths;
+ {
+ std::unordered_set<std::shared_ptr<ComputeNode>> visitedNodes;
+
+ // DFS for the dependencies
+ std::stack<std::shared_ptr<ComputeNode>> nodeStack;
+ depths[graph->getSource()] = 0;
+ nodeStack.push(graph->getSource());
+ while (!nodeStack.empty())
+ {
+ std::shared_ptr<ComputeNode> currNode = nodeStack.top();
+ int currLevel = depths[currNode];
+ nodeStack.pop();
+
+ // Add children to stack if not yet visited
+ if (adjList.count(currNode) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(currNode);
+ for (ComputeNodeSet::const_iterator i = outputNodes.begin(); i != outputNodes.end(); i++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *i;
+ if (visitedNodes.count(childNode) == 0)
+ {
+ visitedNodes.insert(childNode);
+ depths[childNode] = currLevel + 1;
+ nodeStack.push(childNode);
+ }
+ }
+ }
+ }
+ }
+
+ // Identify the set of critical nodes
+ ComputeNodeVector critNodes;
+ for (size_t i = 0; i < nodes.size(); i++)
+ {
+ if (nodes[i]->m_critical)
+ critNodes.push_back(nodes[i]);
+ }
+
+ // Compute the critical adjacency list
+ ComputeNodeAdjList critAdjList;
+ for (size_t i = 0; i < critNodes.size(); i++)
+ {
+ std::unordered_set<std::shared_ptr<ComputeNode>> visitedNodes;
+
+ // DFS for the dependencies
+ std::stack<std::shared_ptr<ComputeNode>> nodeStack;
+ nodeStack.push(critNodes[i]);
+ while (!nodeStack.empty())
+ {
+ std::shared_ptr<ComputeNode> currNode = nodeStack.top();
+ nodeStack.pop();
+
+ // If you can reach one critical node from the other then they are adjacent
+ if (currNode->m_critical)
+ {
+ critAdjList[critNodes[i]].insert(currNode);
+ }
+
+ // Add children to stack if not yet visited
+ if (adjList.count(currNode) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(currNode);
+ for (ComputeNodeSet::const_iterator j = outputNodes.begin(); j != outputNodes.end(); j++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *j;
+ if (visitedNodes.count(childNode) == 0)
+ {
+ visitedNodes.insert(childNode);
+ nodeStack.push(childNode);
+ }
+ }
+ }
+ }
+ }
+
+ // Now we know which critical nodes depend on each other (we are interested in those that aren't)
+ // For every critical pair
+ for (size_t i = 0; i < critNodes.size(); i++)
+ {
+ std::shared_ptr<ComputeNode> srcNode = critNodes[i];
+ for (size_t j = i + 1; j < critNodes.size(); j++)
+ {
+ std::shared_ptr<ComputeNode> destNode = critNodes[j];
+ // If the edge doesn't exist, either way
+ if ((critAdjList.count(srcNode) == 0 || critAdjList.at(srcNode).count(destNode) == 0) &&
+ (critAdjList.count(destNode) == 0 || critAdjList.at(destNode).count(srcNode) == 0))
+ {
+ // Add an edge between the critical nodes in the direction of levels
+ int leveli = depths[srcNode];
+ int levelj = depths[destNode];
+ // If i is below j, then j->i
+ if (leveli > levelj)
+ {
+ results->addEdge(destNode, srcNode);
+ }
+ else
+ {
+ results->addEdge(srcNode, destNode);
+ }
+ }
+ }
+ }
+
+ return results;
+}
+
+std::shared_ptr<ComputeGraph>
+ComputeGraph::transitiveReduce(std::shared_ptr<ComputeGraph> graph)
+{
+ // It's a bad idea to do this method if the graph is cyclic
+ if (isCyclic(graph))
+ return nullptr;
+
+ std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>(*graph);
+
+ // Copy the adjList
+ const ComputeNodeAdjList adjList = results->getAdjList();
+
+ // For every edge in the graph
+ for (ComputeNodeAdjList::const_iterator i = adjList.begin(); i != adjList.end(); i++)
+ {
+ std::shared_ptr<ComputeNode> inputNode = i->first;
+ const ComputeNodeSet& outputNodes = i->second;
+ for (ComputeNodeSet::const_iterator j = outputNodes.begin(); j != outputNodes.end(); j++)
+ {
+ std::shared_ptr<ComputeNode> outputNode = *j;
+ // Edge inputNode->outputNode
+
+ // Remove the edge and see if it still reaches
+ results->removeEdge(inputNode, outputNode);
+
+ // If there still exists a path inputNode->outputNode, leave it removed
+ if (!results->reaches(inputNode, outputNode))
+ {
+ results->addEdge(inputNode, outputNode);
+ }
+ }
+ }
+
+ return results;
+}
+
+std::shared_ptr<ComputeGraph>
+ComputeGraph::nonFunctionalPrune(std::shared_ptr<ComputeGraph> graph)
+{
+ std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>(*graph);
+
+ const ComputeNodeAdjList& adjList = results->getAdjList();
+ const ComputeNodeAdjList& invAdjList = results->getInvAdjList();
+ ComputeNodeVector& nodes = results->getNodes();
+
+ for (size_t i = 0; i < nodes.size(); i++)
+ {
+ std::shared_ptr<ComputeNode> node = nodes[i];
+
+ // These can't be removed (following code would break as they have no inputs/outputs)
+ if (node == graph->getSource() || node == graph->getSink())
+ {
+ continue;
+ }
+
+ // If node not functional and has single input/output it is removeable
+ if (!node->isFunctional())
+ {
+ // Get copies of the inputs and outputs of the node
+ ComputeNodeSet inputs = invAdjList.at(node);
+ ComputeNodeSet outputs = adjList.at(node);
+ if (inputs.size() == 1 && outputs.size() == 1)
+ {
+ // Remove inputs/outputs of node
+ for (ComputeNodeSet::iterator j = inputs.begin(); j != inputs.end(); j++)
+ {
+ results->removeEdge(*j, node);
+ }
+ for (ComputeNodeSet::iterator j = outputs.begin(); j != outputs.end(); j++)
+ {
+ results->removeEdge(node, *j);
+ }
+
+ // Fix the edges
+ for (ComputeNodeSet::iterator j = inputs.begin(); j != inputs.end(); j++)
+ {
+ for (ComputeNodeSet::iterator k = outputs.begin(); k != outputs.end(); k++)
+ {
+ results->addEdge(*j, *k);
+ }
+ }
+
+ // Finally remove the node from the graph
+ nodes.erase(nodes.begin() + i);
+
+ // If node was removed, don't advance
+ i--;
+ }
+ }
+ }
+ return results;
+}
+
+bool
+ComputeGraph::isCyclic(std::shared_ptr<ComputeGraph> graph)
+{
+ // Brute force, DFS every node to find it again
+ const ComputeNodeAdjList& adjList = graph->getAdjList();
+ const ComputeNodeVector& nodes = graph->getNodes();
+ for (size_t i = 0; i < nodes.size(); i++)
+ {
+ std::unordered_set<std::shared_ptr<ComputeNode>> visitedNodes;
+
+ // DFS for the dependencies (start at children instead of itself)
+ std::stack<std::shared_ptr<ComputeNode>> nodeStack;
+ // Add children to stack if not yet visited
+ if (adjList.count(nodes[i]) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(nodes[i]);
+ for (ComputeNodeSet::const_iterator j = outputNodes.begin(); j != outputNodes.end(); j++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *j;
+ if (visitedNodes.count(childNode) == 0)
+ {
+ visitedNodes.insert(childNode);
+ nodeStack.push(childNode);
+ }
+ }
+ }
+ while (!nodeStack.empty())
+ {
+ std::shared_ptr<ComputeNode> currNode = nodeStack.top();
+ nodeStack.pop();
+
+ // If we revisit a node then it must be cyclic
+ if (nodes[i] == currNode)
+ {
+ return true;
+ }
+
+ // Add children to stack if not yet visited
+ if (adjList.count(currNode) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(currNode);
+ for (ComputeNodeSet::const_iterator j = outputNodes.begin(); j != outputNodes.end(); j++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *j;
+ if (visitedNodes.count(childNode) == 0)
+ {
+ visitedNodes.insert(childNode);
+ nodeStack.push(childNode);
+ }
+ }
+ }
+ }
+ }
+
+ return false;
+}
+
+std::unordered_map<std::shared_ptr<ComputeNode>, std::string>
+ComputeGraph::getUniqueNames(std::shared_ptr<ComputeGraph> graph, bool apply)
+{
+ // Produce non colliding names
+ std::unordered_map<std::shared_ptr<ComputeNode>, std::string> nodeNames;
+ std::unordered_map<std::string, int> names;
+ const ComputeNodeVector& nodes = graph->getNodes();
+ for (size_t i = 0; i < nodes.size(); i++)
+ {
+ nodeNames[nodes[i]] = nodes[i]->m_name;
+ names[nodes[i]->m_name]++;
+ }
+ // Adjust names
+ for (std::unordered_map<std::shared_ptr<ComputeNode>, std::string>::iterator it = nodeNames.begin();
+ it != nodeNames.end(); it++)
+ {
+ int nameIter = 0;
+ std::string currName = it->second;
+ // If we can find a node with this name, increment name counter and try again
+ while (names[currName] > 1)
+ {
+ names[currName]--;
+ currName = it->second + std::to_string(nameIter);
+ names[currName]++;
+ nameIter++;
+ }
+ nodeNames[it->first] = currName;
+ }
+ if (apply)
+ {
+ for (std::shared_ptr<ComputeNode> node : nodes)
+ {
+ node->m_name = nodeNames[node];
+ }
+ }
+ return nodeNames;
+}
+
+std::unordered_map<std::shared_ptr<ComputeNode>, double>
+ComputeGraph::getTimes(std::shared_ptr<ComputeGraph> graph)
+{
+ const ComputeNodeAdjList& adjList = graph->getAdjList();
+
+ // Setup a map for total elapsed times at each node
+ std::unordered_map<std::shared_ptr<ComputeNode>, double> times;
+
+ {
+ std::unordered_set<std::shared_ptr<ComputeNode>> visitedNodes;
+
+ // BFS down the tree computing total times
+ std::stack<std::shared_ptr<ComputeNode>> nodeStack;
+ times[graph->getSource()] = 0.0;
+ nodeStack.push(graph->getSource());
+ while (!nodeStack.empty())
+ {
+ std::shared_ptr<ComputeNode> currNode = nodeStack.top();
+ nodeStack.pop();
+
+ // Add children to stack if not yet visited
+ if (adjList.count(currNode) != 0)
+ {
+ const ComputeNodeSet& outputNodes = adjList.at(currNode);
+ for (ComputeNodeSet::const_iterator i = outputNodes.begin(); i != outputNodes.end(); i++)
+ {
+ std::shared_ptr<ComputeNode> childNode = *i;
+ if (visitedNodes.count(childNode) == 0)
+ {
+ visitedNodes.insert(childNode);
+ times[childNode] = times[currNode] + childNode->m_elapsedTime;
+ nodeStack.push(childNode);
+ }
+ }
+ }
+ }
+ }
+ return times;
+}
+
+ComputeNodeList
+ComputeGraph::getCriticalPath(std::shared_ptr<ComputeGraph> graph)
+{
+ std::unordered_map<std::shared_ptr<ComputeNode>, double> times = getTimes(graph);
+
+ // Now backtrack to acquire the path of longest duration
+ const ComputeNodeAdjList& invAdjList = graph->getInvAdjList();
+ std::list<std::shared_ptr<ComputeNode>> results;
+ {
+ std::shared_ptr<ComputeNode> currNode = graph->getSink();
+ // Starting from the sink, choose the input with the longest total elapsed time
+ while (currNode != graph->getSource())
+ {
+ results.push_front(currNode);
+ std::shared_ptr<ComputeNode> longestNode = nullptr;
+ double maxTime = 0.0;
+
+ // For every parent
+ if (invAdjList.count(currNode) != 0)
+ {
+ const ComputeNodeSet& inputNodes = invAdjList.at(currNode);
+ for (ComputeNodeSet::const_iterator i = inputNodes.begin(); i != inputNodes.end(); i++)
+ {
+ std::shared_ptr<ComputeNode> parentNode = *i;
+ if (times[parentNode] >= maxTime)
+ {
+ maxTime = times[parentNode];
+ longestNode = parentNode;
+ }
+ }
+ }
+
+ // Move to the node with the longest time
+ currNode = longestNode;
+ }
+ }
+ results.push_front(graph->getSource());
+ return results;
+}
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkComputeGraph.h b/Source/Common/ComputeGraph/imstkComputeGraph.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae2422f16a821681c605afac7783040e014a6da4
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeGraph.h
@@ -0,0 +1,242 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#pragma once
+
+#include "imstkComputeNode.h"
+
+#include <list>
+#include <memory>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+namespace imstk
+{
+using ComputeNodeVector = std::vector<std::shared_ptr<ComputeNode>>;
+using ComputeNodeList = std::list<std::shared_ptr<ComputeNode>>;
+using ComputeNodeSet = std::unordered_set<std::shared_ptr<ComputeNode>>;
+using ComputeNodeAdjList = std::unordered_map<std::shared_ptr<ComputeNode>, ComputeNodeSet>;
+
+///
+/// \brief Base class for compute graphs, a collection of computational nodes
+/// edges can exist between graphs. Maintains nodes, adjList, and invAdjList
+/// TODO: move algorithms out of the container like STL
+///
+class ComputeGraph
+{
+public:
+ ComputeGraph(std::string sourceName = "Source", std::string sinkName = "Sink");
+ virtual ~ComputeGraph() = default;
+
+public:
+ std::shared_ptr<ComputeNode> getSource() { return m_source; }
+ std::shared_ptr<ComputeNode> getSink() { return m_sink; }
+
+ ///
+ /// \brief Get the nodes belonging to this graph
+ ///
+ ComputeNodeVector& getNodes() { return m_nodes; }
+
+ ///
+ /// \brief Get the edges belonging to this graph
+ ///
+ ComputeNodeAdjList& getAdjList() { return m_adjList; }
+
+ ///
+ /// \brief Get the inverse edges belonging to this graph
+ ///
+ ComputeNodeAdjList& getInvAdjList() { return m_invAdjList; }
+
+ // Node operations
+public:
+ ///
+ /// \brief Linear search for node by name within this graph
+ ///
+ ComputeNodeVector::iterator findNode(std::string name);
+ ComputeNodeVector::const_iterator findNode(std::string name) const;
+ ///
+ /// \brief Linear search for node within this graph
+ ///
+ ComputeNodeVector::iterator findNode(std::shared_ptr<ComputeNode> node);
+ ComputeNodeVector::const_iterator findNode(std::shared_ptr<ComputeNode> node) const;
+
+ ///
+ /// \brief Check if node exists in this graph
+ ///
+ bool containsNode(std::shared_ptr<ComputeNode> node) const;
+
+ ///
+ /// \brief Returns the start of the node container
+ ///
+ ComputeNodeVector::iterator nodesEnd() { return m_nodes.end(); }
+ ComputeNodeVector::const_iterator nodesEnd() const { return m_nodes.end(); }
+ ///
+ /// \brief Returns the start of the node container
+ ///
+ ComputeNodeVector::iterator nodesBegin() { return m_nodes.begin(); }
+ ComputeNodeVector::const_iterator nodesBegin() const { return m_nodes.begin(); }
+
+ ///
+ /// \brief Adds a node to the graph, returns if successful.
+ /// Returns false and fails if node already exists in graph
+ ///
+ bool addNode(std::shared_ptr<ComputeNode> node);
+
+ ///
+ /// \brief Creates a node for the function, adds it to the graph
+ ///
+ std::shared_ptr<ComputeNode> addFunction(std::string name, std::function<void()> func);
+
+ ///
+ /// \brief Removes node from the graph and all relevant edges
+ /// Returns false and fails if node is not present in graph
+ /// "Expensive" operation as no inverse adjacency list is maintained
+ /// full search done for relevant edges
+ ///
+ bool removeNode(std::shared_ptr<ComputeNode> node);
+
+ ///
+ /// \brief Removes node from the graph and all relevant edges
+ /// any incoming edges to the node are redirected to all its outputs
+ ///
+ bool removeNodeAndFix(std::shared_ptr<ComputeNode> node);
+
+ ///
+ /// \brief newNode gets placed after refNode & added to the graph. newNode takes on refNode's outputs
+ ///
+ void insertAfter(std::shared_ptr<ComputeNode> refNode, std::shared_ptr<ComputeNode> newNode);
+
+ ///
+ /// \brief newNode gets placed before refNode & added to the graph. newNode takes on refNode's inputs
+ ///
+ void insertBefore(std::shared_ptr<ComputeNode> refNode, std::shared_ptr<ComputeNode> newNode);
+
+ // Edge operations
+public:
+ ///
+ /// \brief Returns whether or not this graph contains the given directed edge
+ ///
+ bool containsEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode);
+
+ ///
+ /// \brief Adds a directed edge to the graph (doesn't check if nodes are in graph)
+ ///
+ void addEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode);
+
+ ///
+ /// \brief Attachs another ComputeGraph as a subgraph (copies nodes and edges, then connects source->subgraph::source, subgraph::sink->sink),
+ /// source and sink must exist in this graph. Also serves as a graph sum between this and subgraph
+ ///
+ void nestGraph(std::shared_ptr<ComputeGraph> subgraph, std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink);
+
+ ///
+ /// \brief Removes an edge from the graph (removes from adjList and invAdjList, cleans)
+ ///
+ void removeEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode);
+
+ ///
+ /// \brief Returns true if srcNode reaches destNode
+ ///
+ bool reaches(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode);
+
+public:
+ ///
+ /// \brief Removes all nodes & edges from the graph. Source and sink nodes maintained.
+ ///
+ void clear();
+
+ ///
+ /// \brief Removes all edges from the graph
+ ///
+ void clearEdges()
+ {
+ m_adjList.clear();
+ m_invAdjList.clear();
+ }
+
+ // Graph algorithms, todo: Move into filtering module
+public:
+ ///
+ /// \brief Graph sum, shared references are considered identical nodes, source/sink of results invalidated/nullptr
+ ///
+ static std::shared_ptr<ComputeGraph> sum(std::shared_ptr<ComputeGraph> graphA, std::shared_ptr<ComputeGraph> graphB);
+
+ ///
+ /// \brief Topological sort of all nodes within graph
+ ///
+ static std::shared_ptr<ComputeNodeList> topologicalSort(std::shared_ptr<ComputeGraph> graph);
+
+ ///
+ /// \brief Makes sure no two *critical nodes* run at the same time by establishing an edge between them.
+ ///
+ static std::shared_ptr<ComputeGraph> resolveCriticalNodes(std::shared_ptr<ComputeGraph> graph);
+
+ ///
+ /// \brief Remove redudant edges. Removal is such that all vertices are still reachable and graph goes from source->sink
+ /// returns nullptr if failed. Only fails if graph is cyclic.
+ ///
+ static std::shared_ptr<ComputeGraph> transitiveReduce(std::shared_ptr<ComputeGraph> graph);
+
+ ///
+ /// \brief Removes nullptr/nonfunctional ComputeNodes that don't split/join
+ ///
+ static std::shared_ptr<ComputeGraph> nonFunctionalPrune(std::shared_ptr<ComputeGraph> graph);
+
+ ///
+ /// \brief Simplifies the graph in a way that retains functionality
+ /// Performs Nonfunctional Prune and Transitive Reduction removing redundant edges and non functional nodes
+ ///
+ static std::shared_ptr<ComputeGraph> reduce(std::shared_ptr<ComputeGraph> graph)
+ {
+ return nonFunctionalPrune(transitiveReduce(graph));
+ }
+
+ ///
+ /// \brief Returns if Graph is cyclic or not
+ ///
+ static bool isCyclic(std::shared_ptr<ComputeGraph> graph);
+
+ ///
+ /// \brief Iterates the names with numeric postfix to generate uniques
+ ///
+ static std::unordered_map<std::shared_ptr<ComputeNode>, std::string> getUniqueNames(std::shared_ptr<ComputeGraph> graph, bool apply = false);
+
+ ///
+ /// \brief Gets the completion times of each node (source = ~0s)
+ ///
+ static std::unordered_map<std::shared_ptr<ComputeNode>, double> getTimes(std::shared_ptr<ComputeGraph> graph);
+
+ ///
+ /// \brief Computes the critical path
+ ///
+ static ComputeNodeList getCriticalPath(std::shared_ptr<ComputeGraph> graph);
+
+protected:
+ ComputeNodeVector m_nodes;
+ ComputeNodeAdjList m_adjList; ///> This gives the outputs of every node
+ ComputeNodeAdjList m_invAdjList; ///> This gives the inputs of every node
+ // inv adjlist?
+
+ std::shared_ptr<ComputeNode> m_source = nullptr;
+ std::shared_ptr<ComputeNode> m_sink = nullptr;
+};
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkComputeGraphController.h b/Source/Common/ComputeGraph/imstkComputeGraphController.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc1aab113556c6ef8c6a7aad5fcafd9d4b563922
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeGraphController.h
@@ -0,0 +1,58 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#pragma once
+#include <memory>
+
+namespace imstk
+{
+class ComputeGraph;
+
+///
+/// \brief Base class for compute graph controllers which are responsible
+/// for executing the compute graph
+///
+class ComputeGraphController
+{
+public:
+ ComputeGraphController() = default;
+ virtual ~ComputeGraphController() = default;
+
+public:
+ virtual void setComputeGraph(std::shared_ptr<ComputeGraph> graph) { this->m_graph = graph; }
+
+ std::shared_ptr<ComputeGraph> getComputeGraph() const { return m_graph; }
+
+ ///
+ /// \brief Initialization of the ComputeGraphController, good for anything the controller might need
+ /// to do before execution
+ ///
+ virtual void init() { };
+
+ ///
+ /// \brief Executes the ComputeGraph
+ ///
+ virtual void execute() = 0;
+
+protected:
+ std::shared_ptr<ComputeGraph> m_graph = nullptr;
+};
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..895416cda7d96344ef68cbacf830c213dbd1b903
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp
@@ -0,0 +1,104 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#include "imstkComputeGraphVizWriter.h"
+#include "imstkComputeGraph.h"
+#include "imstkLogger.h"
+#include <fstream>
+
+namespace imstk
+{
+void
+ComputeGraphVizWriter::write()
+{
+ if (m_inputGraph == nullptr)
+ {
+ LOG(WARNING) << "No input set, unable to write ComputeGraph";
+ return;
+ }
+
+ // This little writer fails when nodes are named the same, some more advanced
+ // method would have to be used to rename nodes that are the same with consistent
+ // numbering down the tree
+ std::ofstream file;
+ file.open(m_fileName);
+
+ if (!file.is_open() || file.fail())
+ {
+ return;
+ }
+
+ file <<
+ "digraph imstkDependency\n"
+ "{\n"
+ "style=filled;\n"
+ "color=lightgrey;\n"
+ "node [style=filled,color=cornflowerblue];\n"
+ "edge[arrowhead=vee, arrowtail=inv, arrowsize=.7, color=grey20]\n";
+
+ // Ensure unique names
+ std::unordered_map<std::shared_ptr<ComputeNode>, std::string> nodeNames = ComputeGraph::getUniqueNames(m_inputGraph);
+
+ // Compute the critical path (ie: longest path in duration)
+ ComputeNodeList critPath;
+ if (m_highlightCriticalPath)
+ {
+ critPath = ComputeGraph::getCriticalPath(m_inputGraph);
+ }
+ auto edgeExists = [&](const std::shared_ptr<ComputeNode>& a, const std::shared_ptr<ComputeNode>& b)
+ {
+ ComputeNodeList::iterator srcNode = std::find(critPath.begin(), critPath.end(), a);
+ // If srcNode was found and the next node is b
+ return (srcNode != critPath.end() && *std::next(srcNode) == b);
+ };
+
+ // Compute completion times of each node
+ if (m_writeTimes)
+ {
+ std::unordered_map<std::shared_ptr<ComputeNode>, double> times = ComputeGraph::getTimes(m_inputGraph);
+ for (auto i : nodeNames)
+ {
+ nodeNames[i.first] = i.second + "(" + std::to_string(times[i.first]) + "ms)";
+ }
+ }
+
+ // Write out all the edges
+ const ComputeNodeAdjList& adjList = m_inputGraph->getAdjList();
+ for (ComputeNodeAdjList::const_iterator it = adjList.begin(); it != adjList.end(); it++)
+ {
+ std::shared_ptr<ComputeNode> srcNode = it->first;
+ const ComputeNodeSet& outputNodes = it->second;
+ for (ComputeNodeSet::const_iterator jt = outputNodes.begin(); jt != outputNodes.end(); jt++)
+ {
+ std::shared_ptr<ComputeNode> destNode = *jt;
+ file << '\"' << nodeNames[srcNode] << '\"' << " -> " << '\"' << nodeNames[destNode] << '\"';
+
+ if (m_highlightCriticalPath && edgeExists(srcNode, destNode))
+ {
+ file << "[color=red]";
+ }
+ file << std::endl;
+ }
+ }
+ file << "}\n";
+ file.close();
+}
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h
new file mode 100644
index 0000000000000000000000000000000000000000..8d07c9692662427667ccbe7973a8c28e0002d32c
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h
@@ -0,0 +1,82 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#pragma once
+
+#include <memory>
+#include <string>
+
+namespace imstk
+{
+class ComputeGraph;
+
+///
+/// \class ComputeGraph
+///
+/// \brief Writes a ComputeGraph to an svg file. Produces unique node names from duplicates with postfix.
+/// May also highlight cyclic graph edges.
+///
+class ComputeGraphVizWriter
+{
+public:
+ ComputeGraphVizWriter() = default;
+ virtual ~ComputeGraphVizWriter() = default;
+
+public:
+ ///
+ /// \brief The graph to write
+ ///
+ void setInput(std::shared_ptr<ComputeGraph> graph) { this->m_inputGraph = graph; }
+
+ ///
+ /// \brief The filename and path to write too
+ ///
+ void setFileName(std::string fileName) { this->m_fileName = fileName; }
+
+ ///
+ /// \brief If on, will highlight the critical path in red
+ ///
+ void setHighlightCriticalPath(bool highlightCriticalPath) { this->m_highlightCriticalPath = highlightCriticalPath; }
+
+ ///
+ /// \brief If on, will write the time the node completed in the node name
+ /// This isn't exact timing (completion times of each node) as it was run in iMSTK,
+ /// this is timing purely from the timers around the execution of each step
+ ///
+ void setWriteTimes(bool writeTimes) { this->m_writeTimes = writeTimes; }
+
+ std::shared_ptr<ComputeGraph> getInput() const { return m_inputGraph; }
+ const std::string& getFileName() const { return m_fileName; }
+ bool getHighlightCriticalPath() const { return m_highlightCriticalPath; }
+ bool getWriteTimes() const { return m_writeTimes; }
+
+ ///
+ /// \brief Writes the graph to a file given the filename
+ ///
+ void write();
+
+private:
+ std::shared_ptr<ComputeGraph> m_inputGraph = nullptr;
+ std::string m_fileName = "";
+ bool m_highlightCriticalPath = false;
+ bool m_writeTimes = false;
+};
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkComputeNode.cpp b/Source/Common/ComputeGraph/imstkComputeNode.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da5d7777fb2e702ebc9d6f1b103450e2161f9cc7
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeNode.cpp
@@ -0,0 +1,50 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#include "imstkComputeNode.h"
+#include "imstkTimer.h"
+
+namespace imstk
+{
+void
+ComputeNode::execute()
+{
+ if (m_enabled && m_func != nullptr)
+ {
+ if (!m_enableBenchmarking)
+ {
+ m_func();
+ }
+ else
+ {
+ StopWatch timer;
+ timer.start();
+ m_func();
+ m_elapsedTime = timer.getTimeElapsed();
+ timer.stop();
+ }
+ }
+ else
+ {
+ m_elapsedTime = 0.0;
+ }
+}
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkComputeNode.h b/Source/Common/ComputeGraph/imstkComputeNode.h
new file mode 100644
index 0000000000000000000000000000000000000000..31c519f869a050ab9f36d971b953df1539f5b50c
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkComputeNode.h
@@ -0,0 +1,67 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#pragma once
+
+#include <functional>
+#include <string>
+
+namespace imstk
+{
+///
+/// \class ComputeNode
+///
+/// \brief Base class for compute graph nodes
+///
+class ComputeNode
+{
+public:
+ ComputeNode() = default;
+ ComputeNode(std::function<void()> func, std::string name = "none", bool critical = false) :
+ m_func(func), m_name(name), m_critical(critical), m_elapsedTime(0.0)
+ {
+ }
+
+ virtual ~ComputeNode() = default;
+
+public:
+ void setFunction(std::function<void()> func) { this->m_func = func; }
+ void setEnabled(bool enabled) { this->m_enabled = enabled; }
+
+ ///
+ /// \brief Returns true if function is nullptr
+ ///
+ bool isFunctional() const { return m_func != nullptr; }
+
+ // Executes
+ virtual void execute();
+
+public:
+ std::string m_name = "none";
+ bool m_enabled = true;
+ bool m_critical = false;
+ double m_elapsedTime = 0.0;
+ bool m_enableBenchmarking = false;
+
+protected:
+ std::function<void()> m_func = nullptr; ///> Don't allow user to call directly (must use execute)
+};
+}
\ No newline at end of file
diff --git a/Source/SceneEntities/Objects/deprecated/imstkPbdRigidObject.cpp b/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.cpp
similarity index 68%
rename from Source/SceneEntities/Objects/deprecated/imstkPbdRigidObject.cpp
rename to Source/Common/ComputeGraph/imstkSequentialComputeGraphController.cpp
index a70731e37ffd09107a38272930b3f615c45515e4..c4cd79d91110fcbed0d1b5443cf934d89ce6471d 100644
--- a/Source/SceneEntities/Objects/deprecated/imstkPbdRigidObject.cpp
+++ b/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.cpp
@@ -19,20 +19,25 @@
=========================================================================*/
-#include <memory>
-
-#include "imstkPbdRigidObject.h"
-#include "imstkGeometry.h"
-#include "imstkGeometryMap.h"
-#include "imstkPbdModel.h"
-
-#include <g3log/g3log.hpp>
+#include "imstkSequentialComputeGraphController.h"
+#include "imstkComputeGraph.h"
+#include "imstkComputeNode.h"
namespace imstk
{
void
-PbdRigidObject::updatePbdStates()
+SequentialComputeGraphController::init()
+{
+ m_executionOrderedNodes = ComputeGraph::topologicalSort(m_graph);
+}
+
+void
+SequentialComputeGraphController::execute()
{
- m_pbdModel->updatePbdStateFromPhysicsGeometry();
+ // Sequential
+ for (std::shared_ptr<ComputeNode> node : *m_executionOrderedNodes)
+ {
+ node->execute();
+ }
}
-} //imstk
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h b/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h
new file mode 100644
index 0000000000000000000000000000000000000000..c8e1f2156b32f3787fa1f8b4b9d31aae19062da6
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h
@@ -0,0 +1,50 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#pragma once
+#include "imstkComputeGraphController.h"
+#include <list>
+
+namespace imstk
+{
+class ComputeNode;
+
+///
+/// \class SequentialComputeGraphController
+///
+/// \brief This class executes a ComputeGraph by first topologically sorting them (Khans algorithm)
+/// then sequentially running them
+///
+class SequentialComputeGraphController : public ComputeGraphController
+{
+public:
+ ///
+ /// \brief Sorts the computational nodes
+ ///
+ void init() override;
+
+ void execute() override;
+
+private:
+ // The current nodes to execute, ordered
+ std::shared_ptr<std::list<std::shared_ptr<ComputeNode>>> m_executionOrderedNodes;
+};
+};
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkTbbComputeGraphController.cpp b/Source/Common/ComputeGraph/imstkTbbComputeGraphController.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..59bf013fd20714dd6f2a8045c310ed6941cd4a64
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkTbbComputeGraphController.cpp
@@ -0,0 +1,102 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+
+#include "imstkTbbComputeGraphController.h"
+#include "imstkComputeGraph.h"
+#include <tbb/tbb.h>
+
+namespace imstk
+{
+class ComputeNodeTask : public tbb::task
+{
+public:
+ ComputeNodeTask(std::shared_ptr<ComputeNode> node) : m_node(node) { }
+
+public:
+ task* execute() override
+ {
+ __TBB_ASSERT(ref_count() == 0, NULL);
+
+ m_node->execute();
+
+ for (size_t i = 0; i < successors.size(); i++)
+ {
+ if (successors[i]->decrement_ref_count() == 0)
+ {
+ spawn(*successors[i]);
+ }
+ }
+ return NULL;
+ }
+
+public:
+ std::shared_ptr<ComputeNode> m_node = nullptr;
+ std::vector<ComputeNodeTask*> successors;
+};
+
+void
+TbbComputeGraphController::execute()
+{
+ // Create a Task for every node
+ const ComputeNodeVector& nodes = m_graph->getNodes();
+ if (nodes.size() == 0)
+ {
+ return;
+ }
+
+ // Create a task for every node
+ std::unordered_map<std::shared_ptr<ComputeNode>, ComputeNodeTask*> tasks;
+ tasks.reserve(nodes.size());
+
+ for (size_t i = 0; i < nodes.size(); i++)
+ {
+ std::shared_ptr<ComputeNode> node = nodes[i];
+ tasks[node] = new (tbb::task::allocate_root())ComputeNodeTask(node);
+ }
+ // Increment successor reference counts
+ const ComputeNodeAdjList& adjList = m_graph->getAdjList();
+ // For every node in graph
+ for (size_t i = 0; i < nodes.size(); i++)
+ {
+ // If it contains outputs
+ if (adjList.count(nodes[i]) != 0)
+ {
+ // For every output
+ const ComputeNodeSet& outputNodes = adjList.at(nodes[i]);
+ for (ComputeNodeSet::const_iterator it = outputNodes.begin(); it != outputNodes.end(); it++)
+ {
+ // Lookup the task of the node
+ ComputeNodeTask* successor = tasks[*it];
+ tasks[nodes[i]]->successors.push_back(successor);
+ successor->increment_ref_count();
+ }
+ }
+ }
+
+ ComputeNodeTask* startTask = tasks[m_graph->getSource()];
+ ComputeNodeTask* finalTask = tasks[m_graph->getSink()];
+ // Extra ref count on the final task
+ finalTask->increment_ref_count();
+ finalTask->spawn_and_wait_for_all(*startTask);
+ finalTask->execute(); // Execute final task explicitly
+ tbb::task::destroy(*finalTask);
+}
+}
\ No newline at end of file
diff --git a/Source/Common/ComputeGraph/imstkTbbComputeGraphController.h b/Source/Common/ComputeGraph/imstkTbbComputeGraphController.h
new file mode 100644
index 0000000000000000000000000000000000000000..2847e5e5116cb35743830402ab94efa15a29b6b0
--- /dev/null
+++ b/Source/Common/ComputeGraph/imstkTbbComputeGraphController.h
@@ -0,0 +1,38 @@
+/*=========================================================================
+
+Library: iMSTK
+
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0.txt
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+=========================================================================*/
+
+#pragma once
+
+#include "imstkComputeGraphController.h"
+
+namespace imstk
+{
+///
+/// \class TbbComputeGraphController
+///
+/// \brief This class runs an input ComputeGraph in parallel using tbb tasks
+///
+class TbbComputeGraphController : public ComputeGraphController
+{
+public:
+ void execute() override;
+};
+};
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.cpp b/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.cpp
index 436d8d5e5d1442b235e2dbe883637100b8efd326..bf3e64f14a3c79db2cac34987c425a6e707d4ba2 100644
--- a/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.cpp
@@ -20,10 +20,16 @@
=========================================================================*/
#include "imstkAbstractDynamicalModel.h"
+#include "imstkComputeGraph.h"
#include "imstkLogger.h"
namespace imstk
{
+AbstractDynamicalModel::AbstractDynamicalModel(DynamicalModelType type) :
+ m_type(type), m_numDOF(0), m_computeGraph(std::make_shared<ComputeGraph>("AbstractDynamicalModel_Source", "AbstractDynamicalModel_Sink"))
+{
+}
+
bool
AbstractDynamicalModel::isGeometryValid(const std::shared_ptr<Geometry> geometry)
{
@@ -65,4 +71,29 @@ AbstractDynamicalModel::setModelGeometry(std::shared_ptr<Geometry> geometry)
LOG(WARNING) << "Invalid geometry for Model";
}
}
+
+void
+AbstractDynamicalModel::initGraphEdges()
+{
+ m_computeGraph->clearEdges();
+ initGraphEdges(m_computeGraph->getSource(), m_computeGraph->getSink());
+}
+
+void
+AbstractDynamicalModel::initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
+{
+ m_computeGraph->addEdge(source, sink);
+}
+
+std::shared_ptr<ComputeNode>
+AbstractDynamicalModel::addFunction(std::string name, std::function<void()> func)
+{
+ return m_computeGraph->addFunction(name, func);
+}
+
+void
+AbstractDynamicalModel::addEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode)
+{
+ m_computeGraph->addEdge(srcNode, destNode);
+}
} // imstk
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.h b/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.h
index 6d7a547292831a79e2f71e79d414da4e3cd3b5cb..2864f39f157de43d011b5d11926f8556caf21ce8 100644
--- a/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkAbstractDynamicalModel.h
@@ -26,9 +26,13 @@
#include <set>
#include <string>
+#include <unordered_map>
namespace imstk
{
+class ComputeGraph;
+class ComputeNode;
+
///
/// \brief Type of the time dependent mathematical model
///
@@ -58,7 +62,6 @@ enum class TimeSteppingType
class AbstractDynamicalModel
{
public:
-
///
/// \brief Type of the update of the state of the body
///
@@ -71,17 +74,13 @@ public:
None
};
-public:
- ///
- /// \brief Constructor
- ///
- AbstractDynamicalModel(DynamicalModelType type = DynamicalModelType::None) : m_type(type), m_numDOF(0) {}
+protected:
+ AbstractDynamicalModel(DynamicalModelType type = DynamicalModelType::None);
- ///
- /// \brief Destructor
- ///
+public:
virtual ~AbstractDynamicalModel() = default;
+public:
///
/// \brief Reset the current state to the initial state
///
@@ -93,6 +92,8 @@ public:
std::size_t getNumDegreeOfFreedom() const { return m_numDOF; }
void setNumDegreeOfFreedom(const size_t nDof) { m_numDOF = nDof; }
+ std::shared_ptr<ComputeGraph> getComputeGraph() const { return m_computeGraph; }
+
///
/// \brief Get the type of the object
///
@@ -104,9 +105,9 @@ public:
virtual void updateBodyStates(const Vectord& q, const StateUpdateType updateType = StateUpdateType::Displacement) = 0;
///
- /// \brief
+ /// \brief Update the geometry of the model
///
- virtual void updatePhysicsGeometry() {}
+ virtual void updatePhysicsGeometry() { }
///
/// \brief Set the time step size
@@ -138,20 +139,43 @@ public:
///
virtual bool initialize() = 0;
+ ///
+ /// \brief Initializes the edges of the graph
+ ///
+ void initGraphEdges();
+
///
/// \brief Set the type of approach used to update the time step size after every frame
///
virtual void setTimeStepSizeType(const TimeSteppingType type) { m_timeStepSizeType = type; }
TimeSteppingType getTimeStepSizeType() { return m_timeStepSizeType; }
+protected:
+ ///
+ /// \brief Adds a function to the computational graph of this model
+ ///
+ std::shared_ptr<ComputeNode> addFunction(std::string name, std::function<void()> func);
+
+ ///
+ /// \brief Adds an edge to the computational graph of this model
+ ///
+ void addEdge(std::shared_ptr<ComputeNode> srcNode, std::shared_ptr<ComputeNode> destNode);
+
+ ///
+ /// \brief Setup connectivity of the compute graph
+ ///
+ virtual void initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink);
+
protected:
DynamicalModelType m_type; ///> Mathematical model type
std::size_t m_numDOF; ///> Total number of degree of freedom
- std::shared_ptr<Geometry> m_geometry; ///> Physics geometry of the model
+ std::shared_ptr<Geometry> m_geometry = nullptr; ///> Physics geometry of the model
std::set<Geometry::Type> m_validGeometryTypes; ///> Valid geometry types of this model
TimeSteppingType m_timeStepSizeType = TimeSteppingType::Fixed;
+
+ std::shared_ptr<ComputeGraph> m_computeGraph = nullptr;
};
} // imstk
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
index fd2f00405633c7fa323f53093af627720964cf38..a88ea318052dbbb62a4d590536ce71cdd62cb815 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
@@ -24,17 +24,20 @@
//imstk
#include "imstkFEMDeformableBodyModel.h"
-#include "imstkStVKForceModel.h"
-#include "imstkLinearFEMForceModel.h"
+#include "imstkComputeGraph.h"
+#include "imstkComputeNode.h"
+#include "imstkConjugateGradient.h"
#include "imstkCorotationalFEMForceModel.h"
+#include "imstkInternalForceModel.h"
#include "imstkIsotropicHyperelasticFEMForceModel.h"
-#include "imstkVegaMeshIO.h"
-#include "imstkVolumetricMesh.h"
+#include "imstkLinearFEMForceModel.h"
+#include "imstkMath.h"
+#include "imstkNewtonSolver.h"
+#include "imstkSolverBase.h"
+#include "imstkStVKForceModel.h"
#include "imstkTimeIntegrator.h"
#include "imstkVegaMeshIO.h"
-#include "imstkNewtonSolver.h"
-#include "imstkInternalForceModel.h"
-#include "imstkMath.h"
+#include "imstkVolumetricMesh.h"
// vega
#include "generateMassMatrix.h"
@@ -56,6 +59,8 @@ FEMDeformableBodyModel::FEMDeformableBodyModel() :
Geometry::Type::TetrahedralMesh,
Geometry::Type::HexahedralMesh
};
+
+ m_solveNode = addFunction("FEMModel_Solve", [&]() { getSolver()->solve(); });
}
FEMDeformableBodyModel::~FEMDeformableBodyModel()
@@ -86,7 +91,6 @@ FEMDeformableBodyModel::configure(const std::string& configFileName)
vegaConfigFileOptions.addOptionOptional("deformationCompliance", &m_FEModelConfig->m_deformationCompliance, m_FEModelConfig->m_deformationCompliance);
vegaConfigFileOptions.addOptionOptional("compressionResistance", &m_FEModelConfig->m_compressionResistance, m_FEModelConfig->m_compressionResistance);
vegaConfigFileOptions.addOptionOptional("inversionThreshold", &m_FEModelConfig->m_inversionThreshold, m_FEModelConfig->m_inversionThreshold);
- vegaConfigFileOptions.addOptionOptional("numberOfThreads", &m_FEModelConfig->m_numberOfThreads, m_FEModelConfig->m_numberOfThreads);
vegaConfigFileOptions.addOptionOptional("gravity", &m_FEModelConfig->m_gravity, m_FEModelConfig->m_gravity);
// Parse the configuration file
@@ -153,47 +157,37 @@ FEMDeformableBodyModel::configure(std::shared_ptr<FEMModelConfig> config)
m_FEModelConfig = config;
}
-void
-FEMDeformableBodyModel::setForceModelConfiguration(std::shared_ptr<FEMModelConfig> fmConfig)
-{
- m_FEModelConfig = fmConfig;
-}
-
-std::shared_ptr<imstk::FEMModelConfig>
-FEMDeformableBodyModel::getForceModelConfiguration() const
+bool
+FEMDeformableBodyModel::initialize()
{
- return m_FEModelConfig;
-}
+ // prerequisite of for successfully initializing
+ CHECK(m_geometry != nullptr && m_FEModelConfig != nullptr) << "DeformableBodyModel::initialize: Physics mesh or force model configuration not set yet!";
-void
-FEMDeformableBodyModel::setInternalForceModel(std::shared_ptr<InternalForceModel> fm)
-{
- m_internalForceModel = fm;
-}
+ // Setup default solver if model doesn't yet have one
+ if (m_solver == nullptr)
+ {
+ // Create a nonlinear system
+ auto nlSystem = std::make_shared<NonLinearSystem>(getFunction(), getFunctionGradient());
-std::shared_ptr<imstk::InternalForceModel>
-FEMDeformableBodyModel::getInternalForceModel() const
-{
- return m_internalForceModel;
-}
+ nlSystem->setUnknownVector(getUnknownVec());
+ nlSystem->setUpdateFunction(getUpdateFunction());
+ nlSystem->setUpdatePreviousStatesFunction(getUpdatePrevStateFunction());
-void
-FEMDeformableBodyModel::setTimeIntegrator(std::shared_ptr<TimeIntegrator> timeIntegrator)
-{
- m_timeIntegrator = timeIntegrator;
-}
+ // Create a linear solver
+ auto linSolver = std::make_shared<ConjugateGradient>();
-std::shared_ptr<imstk::TimeIntegrator>
-FEMDeformableBodyModel::getTimeIntegrator() const
-{
- return m_timeIntegrator;
-}
+ if (linSolver->getType() == imstk::LinearSolver<imstk::SparseMatrixd>::Type::GaussSeidel
+ && isFixedBCImplemented())
+ {
+ LOG(WARNING) << "The GS solver may not be viable!";
+ }
-bool
-FEMDeformableBodyModel::initialize()
-{
- // prerequisite of for successfully initializing
- CHECK(m_geometry != nullptr && m_FEModelConfig != nullptr) << "DeformableBodyModel::initialize: Physics mesh or force model configuration not set yet!";
+ // Create a non-linear solver and add to the scene
+ auto nlSolver = std::make_shared<NewtonSolver>();
+ nlSolver->setLinearSolver(linSolver);
+ nlSolver->setSystem(nlSystem);
+ setSolver(nlSolver);
+ }
auto physicsMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(this->getModelGeometry());
m_vegaPhysicsMesh = VegaMeshIO::convertVolumetricMeshToVegaMesh(physicsMesh);
@@ -775,12 +769,6 @@ FEMDeformableBodyModel::initializeEigenMatrixFromVegaMatrix(const vega::SparseMa
eigenMatrix.makeCompressed();
}
-Vectord&
-FEMDeformableBodyModel::getContactForce()
-{
- return m_Fcontact;
-}
-
void
FEMDeformableBodyModel::setFixedSizeTimeStepping()
{
@@ -799,4 +787,12 @@ FEMDeformableBodyModel::getTimeStep() const
{
return m_timeIntegrator->getTimestepSize();
};
+
+void
+FEMDeformableBodyModel::initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
+{
+ // Setup graph connectivity
+ addEdge(source, m_solveNode);
+ addEdge(m_solveNode, sink);
+}
} // imstk
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
index ce911e5ef0d97f7573ec98b683ecdbd2ec185924..c6c066e621bba9ecd735d0c3a15f0a0bd58f309c 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
@@ -36,10 +36,11 @@ class VolumetricMesh;
namespace imstk
{
+class InternalForceModel;
+class NewtonSolver;
class TimeIntegrator;
+class SolverBase;
class VegaMeshIO;
-class NewtonSolver;
-class InternalForceModel;
struct FEMModelConfig
{
@@ -57,9 +58,6 @@ struct FEMModelConfig
double m_compressionResistance = 500.0;
double m_inversionThreshold = -std::numeric_limits<double>::max();
double m_gravity = 9.81;
-
- // \todo remove from here ?
- int m_numberOfThreads = 4;
};
///
@@ -83,6 +81,7 @@ public:
///
~FEMDeformableBodyModel();
+public:
///
/// \brief Configure the force model from external file
///
@@ -97,20 +96,20 @@ public:
///
/// \brief Set/Get force model configuration
///
- void setForceModelConfiguration(std::shared_ptr<FEMModelConfig> fmConfig);
- std::shared_ptr<FEMModelConfig> getForceModelConfiguration() const;
+ void setForceModelConfiguration(std::shared_ptr<FEMModelConfig> fmConfig) { this->m_FEModelConfig = fmConfig; }
+ std::shared_ptr<FEMModelConfig> getForceModelConfiguration() const { return m_FEModelConfig; }
///
/// \brief Set/Get internal force model
///
- void setInternalForceModel(std::shared_ptr<InternalForceModel> fm);
- std::shared_ptr<InternalForceModel> getInternalForceModel() const;
+ void setInternalForceModel(std::shared_ptr<InternalForceModel> fm) { m_internalForceModel = fm; }
+ std::shared_ptr<InternalForceModel> getInternalForceModel() const { return m_internalForceModel; }
///
/// \brief Set/Get time integrator
///
- void setTimeIntegrator(std::shared_ptr<TimeIntegrator> timeIntegrator);
- std::shared_ptr<TimeIntegrator> getTimeIntegrator() const;
+ void setTimeIntegrator(std::shared_ptr<TimeIntegrator> timeIntegrator) { this->m_timeIntegrator = timeIntegrator; }
+ std::shared_ptr<TimeIntegrator> getTimeIntegrator() const { return m_timeIntegrator; }
///
/// \brief Load the initial conditions of the deformable object
@@ -226,7 +225,7 @@ public:
///
/// \brief Get the contact force vector
///
- Vectord& getContactForce();
+ Vectord& getContactForce() { return m_Fcontact; }
///
/// \brief Returns the unknown vectors
@@ -265,39 +264,54 @@ public:
void disableFixedBC() { m_implementFixedBC = false; };
bool isFixedBCImplemented() const { return m_implementFixedBC; };
+ std::shared_ptr<ComputeNode> getSolveNode() const { return m_solveNode; }
+
+ std::shared_ptr<SolverBase> getSolver() const { return m_solver; }
+ void setSolver(std::shared_ptr<SolverBase> solver) { this->m_solver = solver; }
+
+protected:
+ ///
+ /// \brief Setup the computational graph of FEM
+ ///
+ void initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink) override;
+
protected:
- std::shared_ptr<InternalForceModel> m_internalForceModel; ///> Mathematical model for intenal forces
- std::shared_ptr<TimeIntegrator> m_timeIntegrator; ///> Time integrator
- std::shared_ptr<NonLinearSystem> m_nonLinearSystem; ///> Nonlinear system resulting from TI and force model
+ std::shared_ptr<SolverBase> m_solver = nullptr;
+ std::shared_ptr<InternalForceModel> m_internalForceModel = nullptr; ///> Mathematical model for intenal forces
+ std::shared_ptr<TimeIntegrator> m_timeIntegrator = nullptr; ///> Time integrator
+ std::shared_ptr<NonLinearSystem> m_nonLinearSystem = nullptr; ///> Nonlinear system resulting from TI and force model
- std::shared_ptr<FEMModelConfig> m_FEModelConfig;
+ std::shared_ptr<FEMModelConfig> m_FEModelConfig = nullptr;
/// Matrices typical to a elastodynamics and 2nd order analogous systems
- SparseMatrixd m_M; ///> Mass matrix
- SparseMatrixd m_C; ///> Damping coefficient matrix
- SparseMatrixd m_K; ///> Tangent (derivative of internal force w.r.t displacements) stiffness matrix
- SparseMatrixd m_Keff; ///> Effective stiffness matrix (dependent on internal force model and time integrator)
+ SparseMatrixd m_M; ///> Mass matrix
+ SparseMatrixd m_C; ///> Damping coefficient matrix
+ SparseMatrixd m_K; ///> Tangent (derivative of internal force w.r.t displacements) stiffness matrix
+ SparseMatrixd m_Keff; ///> Effective stiffness matrix (dependent on internal force model and time integrator)
- Vectord m_Finternal; ///> Vector of internal forces
- Vectord m_Feff; ///> Vector of effective forces
- Vectord m_Fcontact; ///> Vector of contact forces
- Vectord m_Fgravity; ///> Vector of gravity forces
- Vectord m_FexplicitExternal; ///> Vector of explicitly defined external forces
- Vectord m_qSol; ///> Vector to maintain solution at each iteration of nonlinear solver
+ Vectord m_Finternal; ///> Vector of internal forces
+ Vectord m_Feff; ///> Vector of effective forces
+ Vectord m_Fcontact; ///> Vector of contact forces
+ Vectord m_Fgravity; ///> Vector of gravity forces
+ Vectord m_FexplicitExternal; ///> Vector of explicitly defined external forces
+ Vectord m_qSol; ///> Vector to maintain solution at each iteration of nonlinear solver
- std::shared_ptr<vega::VolumetricMesh> m_vegaPhysicsMesh; ///> Mesh used for Physics
- std::shared_ptr<vega::SparseMatrix> m_vegaMassMatrix; ///> Vega mass matrix
- std::shared_ptr<vega::SparseMatrix> m_vegaTangentStiffnessMatrix; ///> Vega Tangent stiffness matrix
- std::shared_ptr<vega::SparseMatrix> m_vegaDampingMatrix; ///> Vega Laplacian damping matrix
+ std::shared_ptr<vega::VolumetricMesh> m_vegaPhysicsMesh = nullptr; ///> Mesh used for Physics
+ std::shared_ptr<vega::SparseMatrix> m_vegaMassMatrix = nullptr; ///> Vega mass matrix
+ std::shared_ptr<vega::SparseMatrix> m_vegaTangentStiffnessMatrix = nullptr; ///> Vega Tangent stiffness matrix
+ std::shared_ptr<vega::SparseMatrix> m_vegaDampingMatrix = nullptr; ///> Vega Laplacian damping matrix
- std::vector<std::size_t> m_fixedNodeIds; ///> Nodal IDs of the nodes that are fixed
+ std::vector<std::size_t> m_fixedNodeIds; ///> Nodal IDs of the nodes that are fixed
- StateUpdateType m_updateType = StateUpdateType::DeltaVelocity; ///> Update type of the model
+ StateUpdateType m_updateType = StateUpdateType::DeltaVelocity; ///> Update type of the model
- bool m_damped = false; ///> Viscous or structurally damped system
+ bool m_damped = false; ///> Viscous or structurally damped system
// If this is true, the tangent stiffness and force vector will be modified to
// accommodate (the rows and columns will be nullified) the fixed boundary conditions
bool m_implementFixedBC = true;
+
+private:
+ std::shared_ptr<ComputeNode> m_solveNode = nullptr;
};
} // imstk
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index 188949934a7743c086eb10660386b07ba27b1fda..c03caa4e302e36393ef67acb0db581143303448e 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -19,28 +19,50 @@
=========================================================================*/
-#include "imstkGraph.h"
#include "imstkPbdModel.h"
-#include "imstkTetrahedralMesh.h"
-#include "imstkSurfaceMesh.h"
+#include "imstkComputeGraph.h"
+#include "imstkGeometryUtilities.h"
+#include "imstkGraph.h"
#include "imstkLineMesh.h"
-#include "imstkPbdVolumeConstraint.h"
-#include "imstkPbdDistanceConstraint.h"
-#include "imstkPbdDihedralConstraint.h"
#include "imstkPbdAreaConstraint.h"
#include "imstkPbdBendConstraint.h"
-#include "imstkPbdFETetConstraint.h"
-#include "imstkPbdFEHexConstraint.h"
#include "imstkPbdConstantDensityConstraint.h"
-#include "imstkParallelUtils.h"
-#include "imstkGeometryUtilities.h"
-#include "imstkMeshIO.h"
-//#include "imstkColor.h"
+#include "imstkPbdDihedralConstraint.h"
+#include "imstkPbdDistanceConstraint.h"
+#include "imstkPbdFETetConstraint.h"
+#include "imstkPbdSolver.h"
+#include "imstkPbdVolumeConstraint.h"
+#include "imstkSurfaceMesh.h"
+#include "imstkTetrahedralMesh.h"
+
+#include "imstkTimer.h"
#include <unordered_map>
namespace imstk
{
+PbdModel::PbdModel() : DynamicalModel(DynamicalModelType::PositionBasedDynamics),
+ m_constraints(std::make_shared<PBDConstraintVector>()),
+ m_partitionedConstraints(std::make_shared<std::vector<PBDConstraintVector>>()),
+ m_mass(std::make_shared<StdVectorOfReal>()),
+ m_invMass(std::make_shared<StdVectorOfReal>()),
+ m_Parameters(std::make_shared<PBDModelConfig>())
+{
+ m_validGeometryTypes = {
+ Geometry::Type::PointSet,
+ Geometry::Type::LineMesh,
+ Geometry::Type::SurfaceMesh,
+ Geometry::Type::TetrahedralMesh,
+ Geometry::Type::HexahedralMesh
+ };
+
+ // Setup PBD compute nodes
+ m_integrationPositionNode = addFunction("PbdModel_IntegratePosition", std::bind(&PbdModel::integratePosition, this));
+ m_updateCollisionGeometryNode = addFunction("PbdModel_UpdateCollisionGeometry", std::bind(&PbdModel::updatePhysicsGeometry, this));
+ m_solveConstraintsNode = addFunction("PbdModel_SolveConstraints", [&]() { m_pbdSolver->solve(); }); // Avoids rebinding on solver swap
+ m_updateVelocityNode = addFunction("PbdModel_UpdateVelocity", std::bind(&PbdModel::updateVelocity, this));
+}
+
void
PBDModelConfig::enableConstraint(PbdConstraint::Type type, double stiffness)
{
@@ -58,7 +80,7 @@ PBDModelConfig::enableFEMConstraint(PbdConstraint::Type type, PbdFEMConstraint::
}
void
-PbdModel::configure(const std::shared_ptr<PBDModelConfig>& params)
+PbdModel::configure(std::shared_ptr<PBDModelConfig> params)
{
LOG_IF(FATAL, (!this->getModelGeometry())) << "PbdModel::configure - Set PBD Model geometry before configuration!";
@@ -76,20 +98,14 @@ PbdModel::initialize()
m_previousState = std::make_shared<PbdState>();
m_currentState = std::make_shared<PbdState>();
- bool option[3] = { 1, 0, 0 };
- m_initialState->initialize(m_mesh, option);
- m_previousState->initialize(m_mesh, option);
-
- option[1] = option[2] = 1;
- m_currentState->initialize(m_mesh, option);
-
- m_initialState->setPositions(m_mesh->getVertexPositions());
- m_currentState->setPositions(m_mesh->getVertexPositions());
+ m_initialState->initialize(m_mesh->getVertexPositions());
+ m_previousState->initialize(m_mesh->getVertexPositions());
+ m_currentState->initialize(m_mesh->getVertexPositions());
auto numParticles = m_mesh->getNumVertices();
- m_mass.resize(numParticles, 0);
- m_invMass.resize(numParticles, 0);
+ m_mass->resize(numParticles, 0);
+ m_invMass->resize(numParticles, 0);
setUniformMass(m_Parameters->m_uniformMassValue);
for (auto i : m_Parameters->m_fixedNodeIds)
@@ -99,6 +115,17 @@ PbdModel::initialize()
bool bOK = true; // Return immediately if some constraint failed to initialize
+ // Setup the default pbd solver if none exists
+ if (m_pbdSolver == nullptr)
+ {
+ m_pbdSolver = std::make_shared<PbdSolver>();
+ m_pbdSolver->setIterations(m_Parameters->m_iterations);
+ }
+ m_pbdSolver->setPositions(getCurrentState()->getPositions());
+ m_pbdSolver->setInvMasses(getInvMasses());
+ m_pbdSolver->setConstraints(getConstraints());
+ m_pbdSolver->setPartitionedConstraints(getPartitionedConstraints());
+
// Initialize FEM constraints
for (auto& constraint: m_Parameters->m_FEMConstraints)
{
@@ -155,23 +182,34 @@ PbdModel::initialize()
return bOK;
}
+void
+PbdModel::initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
+{
+ // Setup graph connectivity
+ addEdge(source, m_integrationPositionNode);
+ addEdge(m_integrationPositionNode, m_updateCollisionGeometryNode);
+ addEdge(m_updateCollisionGeometryNode, m_solveConstraintsNode);
+ addEdge(m_solveConstraintsNode, m_updateVelocityNode);
+ addEdge(m_updateVelocityNode, sink);
+}
+
void
PbdModel::computeElasticConstants()
{
- if (std::abs(m_Parameters->m_mu) < MIN_REAL
- && std::abs(m_Parameters->m_lambda) < MIN_REAL)
+ if (std::abs(m_Parameters->femParams->m_mu) < MIN_REAL
+ && std::abs(m_Parameters->femParams->m_lambda) < MIN_REAL)
{
- const auto E = m_Parameters->m_YoungModulus;
- const auto nu = m_Parameters->m_PoissonRatio;
- m_Parameters->m_mu = E / Real(2.0) / (Real(1.0) + nu);
- m_Parameters->m_lambda = E * nu / ((Real(1.0) + nu) * (Real(1.0) - Real(2.0) * nu));
+ const auto E = m_Parameters->femParams->m_YoungModulus;
+ const auto nu = m_Parameters->femParams->m_PoissonRatio;
+ m_Parameters->femParams->m_mu = E / Real(2.0) / (Real(1.0) + nu);
+ m_Parameters->femParams->m_lambda = E * nu / ((Real(1.0) + nu) * (Real(1.0) - Real(2.0) * nu));
}
else
{
- const auto mu = m_Parameters->m_mu;
- const auto lambda = m_Parameters->m_lambda;
- m_Parameters->m_YoungModulus = mu * (Real(3.0) * lambda + Real(2.0) * mu) / (lambda + mu);
- m_Parameters->m_PoissonRatio = lambda / Real(2.0) / (lambda + mu);
+ const auto mu = m_Parameters->femParams->m_mu;
+ const auto lambda = m_Parameters->femParams->m_lambda;
+ m_Parameters->femParams->m_YoungModulus = mu * (Real(3.0) * lambda + Real(2.0) * mu) / (lambda + mu);
+ m_Parameters->femParams->m_PoissonRatio = lambda / Real(2.0) / (lambda + mu);
}
}
@@ -192,9 +230,10 @@ PbdModel::initializeFEMConstraints(PbdFEMConstraint::MaterialType type)
{
auto& tet = elements[k];
auto c = std::make_shared<PbdFEMTetConstraint>(type);
- c->initConstraint(*this, tet[0], tet[1], tet[2], tet[3]);
+ c->initConstraint(*m_initialState->getPositions(),
+ tet[0], tet[1], tet[2], tet[3], m_Parameters->femParams);
lock.lock();
- m_constraints.push_back(std::move(c));
+ m_constraints->push_back(std::move(c));
lock.unlock();
});
return true;
@@ -216,9 +255,10 @@ PbdModel::initializeVolumeConstraints(const double stiffness)
{
auto& tet = elements[k];
auto c = std::make_shared<PbdVolumeConstraint>();
- c->initConstraint(*this, tet[0], tet[1], tet[2], tet[3], stiffness);
+ c->initConstraint(*m_initialState->getPositions(),
+ tet[0], tet[1], tet[2], tet[3], stiffness);
lock.lock();
- m_constraints.push_back(std::move(c));
+ m_constraints->push_back(std::move(c));
lock.unlock();
});
return true;
@@ -238,8 +278,8 @@ PbdModel::initializeDistanceConstraints(const double stiffness)
{
E[i1][i2] = 0;
auto c = std::make_shared<PbdDistanceConstraint>();
- c->initConstraint(*this, i1, i2, stiffness);
- m_constraints.push_back(std::move(c));
+ c->initConstraint(*m_initialState->getPositions(), i1, i2, stiffness);
+ m_constraints->push_back(std::move(c));
}
};
@@ -310,9 +350,9 @@ PbdModel::initializeAreaConstraints(const double stiffness)
{
auto& tri = elements[k];
auto c = std::make_shared<PbdAreaConstraint>();
- c->initConstraint(*this, tri[0], tri[1], tri[2], stiffness);
+ c->initConstraint(*m_initialState->getPositions(), tri[0], tri[1], tri[2], stiffness);
lock.lock();
- m_constraints.push_back(std::move(c));
+ m_constraints->push_back(std::move(c));
lock.unlock();
});
return true;
@@ -338,8 +378,8 @@ PbdModel::initializeBendConstraints(const double stiffness)
}
auto c = std::make_shared<PbdBendConstraint>();
- c->initConstraint(*this, i1, i2, i3, k);
- m_constraints.push_back(std::move(c));
+ c->initConstraint(*m_initialState->getPositions(), i1, i2, i3, k);
+ m_constraints->push_back(std::move(c));
};
// Create constraints
@@ -411,8 +451,8 @@ PbdModel::initializeDihedralConstraints(const double stiffness)
}
}
auto c = std::make_shared<PbdDihedralConstraint>();
- c->initConstraint(*this, tri[2], tri[idx], tri[0], tri[1], stiffness);
- m_constraints.push_back(std::move(c));
+ c->initConstraint(*m_initialState->getPositions(), tri[2], tri[idx], tri[0], tri[1], stiffness);
+ m_constraints->push_back(std::move(c));
}
}
};
@@ -448,8 +488,8 @@ PbdModel::initializeConstantDensityConstraint(const double stiffness)
<< "Constant constraint should come with a mesh!";
auto c = std::make_shared<PbdConstantDensityConstraint>();
- c->initConstraint(*this, stiffness);
- m_constraints.push_back(std::move(c));
+ c->initConstraint(*m_initialState->getPositions(), stiffness);
+ m_constraints->push_back(std::move(c));
return true;
}
@@ -458,10 +498,11 @@ void
PbdModel::partitionConstraints(const bool print)
{
// Form the map { vertex : list_of_constraints_involve_vertex }
+ PBDConstraintVector& allConstraints = *m_constraints;
std::unordered_map<size_t, std::vector<size_t>> vertexConstraints;
- for (size_t constrIdx = 0; constrIdx < m_constraints.size(); ++constrIdx)
+ for (size_t constrIdx = 0; constrIdx < allConstraints.size(); ++constrIdx)
{
- const auto& constr = m_constraints[constrIdx];
+ const auto& constr = allConstraints[constrIdx];
for (const auto& vIds : constr->getVertexIds())
{
vertexConstraints[vIds].push_back(constrIdx);
@@ -470,7 +511,7 @@ PbdModel::partitionConstraints(const bool print)
// Add edges to the constraint graph
// Each edge represent a shared vertex between two constraints
- Graph constraintGraph(m_constraints.size());
+ Graph constraintGraph(allConstraints.size());
for (const auto& kv : vertexConstraints)
{
const auto& constraints = kv.second; // the list of constraints for a vertex
@@ -488,86 +529,65 @@ PbdModel::partitionConstraints(const bool print)
const auto coloring = constraintGraph.doColoring();
const auto& partitionIndices = coloring.first;
const auto numPartitions = coloring.second;
- assert(partitionIndices.size() == m_constraints.size());
+ assert(partitionIndices.size() == allConstraints.size());
- m_partitionedConstraints.resize(0);
- m_partitionedConstraints.resize(static_cast<size_t>(numPartitions));
+ std::vector<PBDConstraintVector>& partitionedConstraints = *m_partitionedConstraints;
+ partitionedConstraints.resize(0);
+ partitionedConstraints.resize(static_cast<size_t>(numPartitions));
for (size_t constrIdx = 0; constrIdx < partitionIndices.size(); ++constrIdx)
{
const auto partitionIdx = partitionIndices[constrIdx];
- m_partitionedConstraints[partitionIdx].push_back(std::move(m_constraints[constrIdx]));
+ partitionedConstraints[partitionIdx].push_back(std::move(allConstraints[constrIdx]));
}
// If a partition has size smaller than the partition threshold, then move its constraints back
// These constraints will be processed sequentially
// Because small size partitions yield bad performance upon running in parallel
- m_constraints.resize(0);
- for (const auto& constraints : m_partitionedConstraints)
+ allConstraints.resize(0);
+ for (const auto& constraints : partitionedConstraints)
{
if (constraints.size() < m_partitionThreshold)
{
for (size_t constrIdx = 0; constrIdx < constraints.size(); ++constrIdx)
{
- m_constraints.push_back(std::move(constraints[constrIdx]));
+ allConstraints.push_back(std::move(constraints[constrIdx]));
}
}
}
// Remove all empty partitions
size_t writeIdx = 0;
- for (size_t readIdx = 0; readIdx < m_partitionedConstraints.size(); ++readIdx)
+ for (size_t readIdx = 0; readIdx < partitionedConstraints.size(); ++readIdx)
{
- if (m_partitionedConstraints[readIdx].size() >= m_partitionThreshold)
+ if (partitionedConstraints[readIdx].size() >= m_partitionThreshold)
{
- m_partitionedConstraints[writeIdx++] = std::move(m_partitionedConstraints[readIdx]);
+ partitionedConstraints[writeIdx++] = std::move(partitionedConstraints[readIdx]);
}
}
- m_partitionedConstraints.resize(writeIdx);
+ partitionedConstraints.resize(writeIdx);
// Print
if (print)
{
size_t numConstraints = 0;
int idx = 0;
- for (const auto& constraints : m_partitionedConstraints)
+ for (const auto& constraints : partitionedConstraints)
{
std::cout << "Partition # " << idx++ << " | # nodes: " << constraints.size() << std::endl;
numConstraints += constraints.size();
}
- std::cout << "Sequential processing # nodes: " << m_constraints.size() << std::endl;
- numConstraints += m_constraints.size();
+ std::cout << "Sequential processing # nodes: " << allConstraints.size() << std::endl;
+ numConstraints += allConstraints.size();
std::cout << "Total constraints: " << numConstraints << " | Graph size: "
<< constraintGraph.size() << std::endl;
}
}
-void
-PbdModel::projectConstraints()
-{
- unsigned int i = 0;
- while (++i < m_Parameters->m_maxIter)
- {
- for (auto c: m_constraints)
- {
- c->solvePositionConstraint(*this);
- }
-
- for (auto& partitionConstraints : m_partitionedConstraints)
- {
- ParallelUtils::parallelFor(partitionConstraints.size(),
- [&](const size_t idx)
- {
- partitionConstraints[idx]->solvePositionConstraint(*this);
- });
- }
- }
-}
-
void
PbdModel::updatePhysicsGeometry()
{
- m_mesh->setVertexPositions(m_currentState->getPositions());
+ m_mesh->setVertexPositions(*m_currentState->getPositions());
}
void
@@ -582,7 +602,7 @@ PbdModel::setTimeStepSizeType(const TimeSteppingType type)
m_timeStepSizeType = type;
if (type == TimeSteppingType::Fixed)
{
- m_Parameters->m_dt = m_Parameters->m_DefaultDt;
+ m_Parameters->m_dt = m_Parameters->m_defaultDt;
}
}
@@ -591,53 +611,51 @@ PbdModel::setUniformMass(const double val)
{
if (val != 0.0)
{
- std::fill(m_mass.begin(), m_mass.end(), val);
- std::fill(m_invMass.begin(), m_invMass.end(), 1 / val);
+ std::fill(m_mass->begin(), m_mass->end(), val);
+ std::fill(m_invMass->begin(), m_invMass->end(), 1.0 / val);
}
else
{
- std::fill(m_invMass.begin(), m_invMass.end(), 0.0);
- std::fill(m_mass.begin(), m_mass.end(), 0.0);
+ std::fill(m_invMass->begin(), m_invMass->end(), 0.0);
+ std::fill(m_mass->begin(), m_mass->end(), 0.0);
}
}
void
PbdModel::setParticleMass(const double val, const size_t idx)
{
+ StdVectorOfReal& masses = *m_mass;
+ StdVectorOfReal& invMasses = *m_invMass;
if (idx < m_mesh->getNumVertices())
{
- m_mass[idx] = val;
- m_invMass[idx] = 1.0 / val;
+ masses[idx] = val;
+ invMasses[idx] = 1.0 / val;
}
}
void
PbdModel::setFixedPoint(const size_t idx)
{
+ StdVectorOfReal& invMasses = *m_invMass;
if (idx < m_mesh->getNumVertices())
{
- m_invMass[idx] = 0;
+ invMasses[idx] = 0.0;
}
}
-double
-PbdModel::getInvMass(const size_t idx) const
-{
- return m_invMass[idx];
-}
-
void
PbdModel::integratePosition()
{
- const auto& accn = m_currentState->getAccelerations();
- auto& prevPos = m_previousState->getPositions();
- auto& pos = m_currentState->getPositions();
- auto& vel = m_currentState->getVelocities();
+ StdVectorOfVec3d& prevPos = *m_previousState->getPositions();
+ StdVectorOfVec3d& pos = *m_currentState->getPositions();
+ StdVectorOfVec3d& vel = *m_currentState->getVelocities();
+ const StdVectorOfVec3d& accn = *m_currentState->getAccelerations();
+ const StdVectorOfReal& invMasses = *m_invMass;
ParallelUtils::parallelFor(m_mesh->getNumVertices(),
[&](const size_t i)
{
- if (std::abs(m_invMass[i]) > MIN_REAL)
+ if (std::abs(invMasses[i]) > MIN_REAL)
{
vel[i] += (accn[i] + m_Parameters->m_gravity) * m_Parameters->m_dt;
prevPos[i] = pos[i];
@@ -649,14 +667,15 @@ PbdModel::integratePosition()
void
PbdModel::updateVelocity()
{
- const auto& prevPos = m_previousState->getPositions();
- const auto& pos = m_currentState->getPositions();
- auto& vel = m_currentState->getVelocities();
+ const StdVectorOfVec3d& prevPos = *m_previousState->getPositions();
+ const StdVectorOfVec3d& pos = *m_currentState->getPositions();
+ StdVectorOfVec3d& vel = *m_currentState->getVelocities();
+ const StdVectorOfReal& invMasses = *m_invMass;
ParallelUtils::parallelFor(m_mesh->getNumVertices(),
[&](const size_t i)
{
- if (std::abs(m_invMass[i]) > MIN_REAL && m_Parameters->m_dt > 0.)
+ if (std::abs(invMasses[i]) > MIN_REAL && m_Parameters->m_dt > 0.0)
{
vel[i] = (pos[i] - prevPos[i]) / m_Parameters->m_dt;
}
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
index a2c908755fd43012f1337ce2c4b1ff79d9210a0a..7f34a3c9250ac4c974e381b6bf6ebac1a076f9d7 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
@@ -21,18 +21,15 @@
#pragma once
-#include <vector>
-#include <Eigen/Dense>
-
-#include "imstkPbdConstraint.h"
-#include "imstkPbdFEMConstraint.h"
#include "imstkDynamicalModel.h"
+#include "imstkPbdCollisionConstraint.h"
+#include "imstkPbdFEMConstraint.h"
#include "imstkPbdState.h"
-#include "imstkMath.h"
namespace imstk
{
class PointSet;
+class PbdSolver;
///
/// \class PBDModelConfig
@@ -40,28 +37,35 @@ class PointSet;
///
struct PBDModelConfig
{
- double m_uniformMassValue = 1.0; ///> Mass properties
- double m_viscousDampingCoeff = 0.01; ///> Viscous damping coefficient [0, 1]
-
- double m_contactStiffness = 1.0; ///> Contact stiffness for collisions
- double m_proximity = 0.1; ///> Proximity for collisions
+ double m_uniformMassValue = 1.0; ///> Mass properties
+ double m_viscousDampingCoeff = 0.01; ///> Viscous damping coefficient [0, 1]
- unsigned int m_maxIter = 10; ///> Max. pbd iterations
- double m_dt = 0.0; ///> Time step size
- double m_DefaultDt = 0.01; ///> Default Time step size
+ std::shared_ptr<PbdCollisionConstraintConfig> collisionParams =
+ std::make_shared<PbdCollisionConstraintConfig>(PbdCollisionConstraintConfig
+ {
+ 0.1, // Proximity
+ 1.0 // Contact stiffness
+ }); ///> Info shared between the collision constraints
- std::vector<std::size_t> m_fixedNodeIds; ///> Nodal IDs of the nodes that are fixed
- Vec3r m_gravity = Vec3r(0, -9.81, 0); ///> Gravity
+ unsigned int m_iterations = 10; ///> Pbd iterations
+ unsigned int m_collisionIterations = 5; ///> Pbd collision iterations
+ double m_dt = 0.0; ///> Time step size
+ double m_defaultDt = 0.01; ///> Default Time step size
- double m_mu = 0; ///> Lame constant, if constraint type is FEM
- double m_lambda = 0; ///> Lame constant, if constraint type is FEM
+ std::vector<std::size_t> m_fixedNodeIds; ///> Nodal IDs of the nodes that are fixed
+ Vec3r m_gravity = Vec3r(0, -9.81, 0); ///> Gravity
- double m_YoungModulus = 1000; ///> FEM parameter, if constraint type is FEM
- double m_PoissonRatio = 0.2; ///> FEM parameter, if constraint type is FEM
+ std::shared_ptr<PbdFEMConstraintConfig> femParams =
+ std::make_shared<PbdFEMConstraintConfig>(PbdFEMConstraintConfig
+ {
+ 0.0, // Lame constant, if constraint type is FEM
+ 0.0, // Lame constant, if constraint type is FEM
+ 1000.0, // FEM parameter, if constraint type is FEM
+ 0.2 // FEM parameter, if constraint type is FEM
+ }); ///> Info shared between the fem constraints
- std::vector<std::pair<PbdConstraint::Type, double>> m_RegularConstraints; ///> Constraints except FEM
- std::vector<std::pair<PbdConstraint::Type,
- PbdFEMConstraint::MaterialType>> m_FEMConstraints; ///> FEM constraints
+ std::vector<std::pair<PbdConstraint::Type, double>> m_RegularConstraints; ///> Constraints except FEM
+ std::vector<std::pair<PbdConstraint::Type, PbdFEMConstraint::MaterialType>> m_FEMConstraints; ///> Constraints except FEM
///
/// \brief Enable a regular constraint (constraint that is not FEM constraint)
@@ -70,7 +74,7 @@ struct PBDModelConfig
void enableConstraint(PbdConstraint::Type type, double stiffness);
///
- /// \brief Enable a FEM constraint with a given FEM material
+ /// \brief Enable a FEM constraint with mu, lambda
///
void enableFEMConstraint(PbdConstraint::Type type, PbdFEMConstraint::MaterialType material);
};
@@ -86,16 +90,7 @@ public:
///
/// \brief Constructor
///
- PbdModel() : DynamicalModel(DynamicalModelType::PositionBasedDynamics)
- {
- m_validGeometryTypes = {
- Geometry::Type::PointSet,
- Geometry::Type::LineMesh,
- Geometry::Type::SurfaceMesh,
- Geometry::Type::TetrahedralMesh,
- Geometry::Type::HexahedralMesh
- };
- }
+ PbdModel();
///
/// \brief Destructor
@@ -105,7 +100,7 @@ public:
///
/// \brief Set simulation parameters
///
- void configure(const std::shared_ptr<PBDModelConfig>& params);
+ void configure(std::shared_ptr<PBDModelConfig> params);
///
/// \brief Get the simulation parameters
@@ -154,11 +149,6 @@ public:
///
bool initializeConstantDensityConstraint(const double stiffness);
- ///
- /// \brief compute delta x (position) and update position
- ///
- void projectConstraints();
-
///
/// \brief Update the model geometry from the newest PBD state
///
@@ -172,13 +162,13 @@ public:
///
/// \brief Returns true if there is at least one constraint
///
- bool hasConstraints() const { return !m_constraints.empty() || !m_partitionedConstraints.empty(); }
+ bool hasConstraints() const { return !m_constraints->empty() || !m_partitionedConstraints->empty(); }
///
/// \brief Set the time step size
///
virtual void setTimeStep(const Real timeStep) override { m_Parameters->m_dt = timeStep; }
- void setDefaultTimeStep(const Real timeStep) { m_Parameters->m_DefaultDt = static_cast<Real>(timeStep); }
+ void setDefaultTimeStep(const Real timeStep) { m_Parameters->m_defaultDt = static_cast<Real>(timeStep); }
///
/// \brief Set the time step size to fixed size
@@ -189,7 +179,17 @@ public:
/// \brief Returns the time step size
///
virtual double getTimeStep() const override { return m_Parameters->m_dt; }
- double getDefaultTimeStep() const { return m_Parameters->m_DefaultDt; }
+ double getDefaultTimeStep() const { return m_Parameters->m_defaultDt; }
+
+ ///
+ /// \brief
+ ///
+ std::shared_ptr<PBDConstraintVector> getConstraints() { return m_constraints; }
+
+ ///
+ /// \brief
+ ///
+ std::shared_ptr<std::vector<PBDConstraintVector>> getPartitionedConstraints() { return m_partitionedConstraints; }
///
/// \brief Set uniform mass to all the nodes
@@ -209,7 +209,9 @@ public:
///
/// \brief Get the inverse of mass of a certain node
///
- double getInvMass(const size_t idx) const;
+ double getInvMass(const size_t idx) const { return (*m_invMass)[idx]; }
+
+ std::shared_ptr<StdVectorOfReal> getInvMasses() { return m_invMass; }
///
/// \brief Time integrate the position
@@ -237,22 +239,46 @@ public:
///
void setConstraintPartitionThreshold(size_t threshold) { m_partitionThreshold = threshold; }
+ std::shared_ptr<PbdSolver> getSolver() const { return m_pbdSolver; }
+
+ void setSolver(std::shared_ptr<PbdSolver> solver) { this->m_pbdSolver = solver; }
+
+ std::shared_ptr<ComputeNode> getIntegratePositionNode() const { return m_integrationPositionNode; }
+
+ std::shared_ptr<ComputeNode> getUpdateCollisionGeometryNode() const { return m_updateCollisionGeometryNode; }
+
+ std::shared_ptr<ComputeNode> getSolveNode() const { return m_solveConstraintsNode; }
+
+ std::shared_ptr<ComputeNode> getUpdateVelocityNode() const { return m_updateVelocityNode; }
+
protected:
///
/// \brief Partition constraints for parallelization
///
void partitionConstraints(const bool print = false);
- size_t m_partitionThreshold = 16; ///> Threshold for constraint partitioning
+ ///
+ /// \brief Setup the computational graph of PBD
+ ///
+ void initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink) override;
+
+protected:
+ size_t m_partitionThreshold = 16; ///> Threshold for constraint partitioning
- std::shared_ptr<PointSet> m_mesh; ///> PointSet on which the pbd model operates on
- std::vector<double> m_mass; ///> Mass of nodes
- std::vector<double> m_invMass; ///> Inverse of mass of nodes
+ std::shared_ptr<PbdSolver> m_pbdSolver = nullptr; ///> PBD solver
+ std::shared_ptr<PointSet> m_mesh = nullptr; ///> PointSet on which the pbd model operates on
+ std::shared_ptr<StdVectorOfReal> m_mass = nullptr; ///> Mass of nodes
+ std::shared_ptr<StdVectorOfReal> m_invMass = nullptr; ///> Inverse of mass of nodes
- using PBDConstraintVector = std::vector<std::shared_ptr<PbdConstraint>>;
+ std::shared_ptr<PBDConstraintVector> m_constraints = nullptr; ///> List of pbd constraints
+ std::shared_ptr<std::vector<PBDConstraintVector>> m_partitionedConstraints = nullptr; ///> List of pbd constraints
+ std::shared_ptr<PBDModelConfig> m_Parameters = nullptr; ///> Model parameters, must be set before simulation
- PBDConstraintVector m_constraints; ///> List of pbd constraints
- std::vector<PBDConstraintVector> m_partitionedConstraints; ///> List of pbd constraints
- std::shared_ptr<PBDModelConfig> m_Parameters; ///> Model parameters, must be set before simulation
+protected:
+ // Computational Nodes
+ std::shared_ptr<ComputeNode> m_integrationPositionNode = nullptr;
+ std::shared_ptr<ComputeNode> m_updateCollisionGeometryNode = nullptr;
+ std::shared_ptr<ComputeNode> m_solveConstraintsNode = nullptr;
+ std::shared_ptr<ComputeNode> m_updateVelocityNode = nullptr;
};
} // imstk
diff --git a/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel.h b/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel.h
index b9d377d7bd34a1ece03e01dfd423c71ada8bbdf4..7c9e020342cfe46e1bfbcbff817a80e89916bcd7 100644
--- a/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel.h
@@ -19,8 +19,7 @@
=========================================================================*/
-#ifndef imstkRigidBodyModel_h
-#define imstkRigidBodyModel_h
+#pragma once
#include "imstkDynamicalModel.h"
#include "imstkRigidBodyState.h"
@@ -175,6 +174,4 @@ private:
void createCube();
void createMesh();
};
-} // imstk
-
-#endif
\ No newline at end of file
+} // imstk
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp b/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp
index ca91fa161179157371ec49e46d519426f85c2c32..2f0ea2549de46dbafbd77778c954d2a0f0583b13 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp
@@ -1,25 +1,28 @@
/*=========================================================================
- Library: iMSTK
+Library: iMSTK
- Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
- & Imaging in Medicine, Rensselaer Polytechnic Institute.
+Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+& Imaging in Medicine, Rensselaer Polytechnic Institute.
- Licensed under the Apache License, Version 2.0 (the "License");
- you may not use this file except in compliance with the License.
- You may obtain a copy of the License at
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
- http://www.apache.org/licenses/LICENSE-2.0.txt
+http://www.apache.org/licenses/LICENSE-2.0.txt
- Unless required by applicable law or agreed to in writing, software
- distributed under the License is distributed on an "AS IS" BASIS,
- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- See the License for the specific language governing permissions and
- limitations under the License.
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
=========================================================================*/
#include "imstkSPHModel.h"
+#include "imstkComputeGraph.h"
+#include "imstkComputeNode.h"
+#include "imstkLogger.h"
#include "imstkParallelUtils.h"
#include "imstkPointSet.h"
@@ -54,6 +57,34 @@ SPHModelConfig::initialize()
m_kernelRadiusSqr = m_kernelRadius * m_kernelRadius;
}
+SPHModel::SPHModel() : DynamicalModel<SPHKinematicState>(DynamicalModelType::SmoothedParticleHydrodynamics)
+{
+ m_validGeometryTypes = { Geometry::Type::PointSet };
+
+ m_findParticleNeighborsNode = addFunction("SPHModel_Partition", std::bind(&SPHModel::findParticleNeighbors, this));
+ m_computeDensityNode = addFunction("SPHModel_ComputeDensity", [&]()
+ {
+ computeNeighborRelativePositions();
+ computeDensity();
+ normalizeDensity();
+ collectNeighborDensity();
+ });
+ m_computePressureAccelNode =
+ addFunction("SPHModel_ComputePressureAccel", std::bind(&SPHModel::computePressureAcceleration, this));
+ m_computeSurfaceTensionNode =
+ addFunction("SPHModel_ComputeSurfaceTensionAccel", std::bind(&SPHModel::computeSurfaceTension, this));
+ m_computeTimeSetpSizeNode =
+ addFunction("SPHModel_ComputeTimestep", std::bind(&SPHModel::computeTimeStepSize, this));
+ m_integrateNode =
+ addFunction("SPHModel_Integrate", [&]()
+ {
+ sumAccels();
+ updateVelocity(getTimeStep());
+ computeViscosity();
+ moveParticles(getTimeStep());
+ });
+}
+
bool
SPHModel::initialize()
{
@@ -80,7 +111,10 @@ SPHModel::initialize()
// Initialize neighbor searcher
m_neighborSearcher = std::make_shared<NeighborSearch>(m_modelParameters->m_NeighborSearchMethod,
- m_modelParameters->m_kernelRadius);
+ m_modelParameters->m_kernelRadius);
+
+ m_pressureAccels = std::make_shared<StdVectorOfVec3d>(getState().getNumParticles());
+ m_surfaceTensionAccels = std::make_shared<StdVectorOfVec3d>(getState().getNumParticles());
return true;
}
@@ -93,19 +127,22 @@ SPHModel::updatePhysicsGeometry()
}
void
-SPHModel::advanceTimeStep()
+SPHModel::initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
{
- findParticleNeighbors();
- computeNeighborRelativePositions();
- computeDensity();
- normalizeDensity();
- collectNeighborDensity();
- computePressureAcceleration();
- computeSurfaceTension();
- computeTimeStepSize();
- updateVelocity(getTimeStep());
- computeViscosity();
- moveParticles(getTimeStep());
+ // Setup graph connectivity
+ addEdge(source, m_findParticleNeighborsNode);
+ addEdge(m_findParticleNeighborsNode, m_computeDensityNode);
+
+ // Pressure, Surface Tension, and time step size can be done in parallel
+ addEdge(m_computeDensityNode, m_computePressureAccelNode);
+ addEdge(m_computeDensityNode, m_computeSurfaceTensionNode);
+ addEdge(m_computeDensityNode, m_computeTimeSetpSizeNode);
+
+ addEdge(m_computePressureAccelNode, m_integrateNode);
+ addEdge(m_computeSurfaceTensionNode, m_integrateNode);
+ addEdge(m_computeTimeSetpSizeNode, m_integrateNode);
+
+ addEdge(m_integrateNode, sink);
}
void
@@ -282,6 +319,7 @@ SPHModel::computePressureAcceleration()
return error > Real(0) ? error : Real(0);
};
+ StdVectorOfVec3d& pressureAccels = *m_pressureAccels;
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
Vec3r accel(0, 0, 0);
@@ -307,7 +345,19 @@ SPHModel::computePressureAcceleration()
}
accel *= m_modelParameters->m_pressureStiffness * m_modelParameters->m_particleMass;
- getState().getAccelerations()[p] = accel;
+ //getState().getAccelerations()[p] = accel;
+ pressureAccels[p] = accel;
+ });
+}
+
+void
+SPHModel::sumAccels()
+{
+ const StdVectorOfVec3d& pressureAccels = *m_pressureAccels;
+ const StdVectorOfVec3d& surfaceTensionAccels = *m_surfaceTensionAccels;
+ ParallelUtils::parallelFor(getState().getNumParticles(),
+ [&](const size_t p) {
+ getState().getAccelerations()[p] = pressureAccels[p] + surfaceTensionAccels[p];
});
}
@@ -366,7 +416,7 @@ SPHModel::computeViscosity()
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
getState().getVelocities()[p] += getState().getDiffuseVelocities()[p];
- });
+ });
}
void
@@ -396,6 +446,7 @@ SPHModel::computeSurfaceTension()
});
// Second, compute surface tension acceleration
+ StdVectorOfVec3d& surfaceTensionAccels = *m_surfaceTensionAccels;
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
const auto& fluidNeighborList = getState().getFluidNeighborLists()[p];
@@ -436,7 +487,8 @@ SPHModel::computeSurfaceTension()
}
accel *= m_modelParameters->m_surfaceTensionStiffness;
- getState().getAccelerations()[p] += accel;
+ //getState().getAccelerations()[p] += accel;
+ surfaceTensionAccels[p] = accel;
});
}
@@ -446,6 +498,6 @@ SPHModel::moveParticles(Real timestep)
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
getState().getPositions()[p] += getState().getVelocities()[p] * timestep;
- });
+ });
}
} // end namespace imstk
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHModel.h b/Source/DynamicalModels/ObjectModels/imstkSPHModel.h
index d5cedcf23b8a2eabe9996d56bbaaed58cc5bae2d..d90a77786203287f9fd1c0a78b782e835b6b9cd8 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHModel.h
@@ -28,6 +28,7 @@
namespace imstk
{
+class ComputeNode;
class PointSet;
///
@@ -92,10 +93,7 @@ public:
///
/// \brief Constructor
///
- SPHModel() : DynamicalModel<SPHKinematicState>(DynamicalModelType::SmoothedParticleHydrodynamics)
- {
- m_validGeometryTypes = { Geometry::Type::PointSet };
- }
+ SPHModel();
///
/// \brief Destructor
@@ -175,10 +173,11 @@ public:
virtual double getTimeStep() const override
{ return static_cast<double>(m_dt); }
+protected:
///
- /// \brief Do one time step simulation
+ /// \brief Setup SPH compute graph connectivity
///
- void advanceTimeStep();
+ virtual void initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink) override;
private:
///
@@ -222,6 +221,11 @@ private:
///
void computePressureAcceleration();
+ ///
+ /// \brief Sum the forces computed in parallel
+ ///
+ void sumAccels();
+
///
/// \brief Update particle velocities due to pressure
///
@@ -244,6 +248,16 @@ private:
///
void moveParticles(const Real timestep);
+protected:
+ std::shared_ptr<ComputeNode> m_findParticleNeighborsNode = nullptr;
+ std::shared_ptr<ComputeNode> m_computeDensityNode = nullptr;
+ std::shared_ptr<ComputeNode> m_computePressureAccelNode = nullptr;
+ std::shared_ptr<ComputeNode> m_computeSurfaceTensionNode = nullptr;
+ std::shared_ptr<ComputeNode> m_computeTimeSetpSizeNode = nullptr;
+ std::shared_ptr<ComputeNode> m_sumAccelsNode = nullptr;
+ std::shared_ptr<ComputeNode> m_integrateNode = nullptr;
+
+private:
std::shared_ptr<PointSet> m_pointSetGeometry;
SPHSimulationState m_simulationState;
@@ -253,5 +267,8 @@ private:
SPHSimulationKernels m_kernels; ///> SPH kernels (must be initialized during model initialization)
std::shared_ptr<SPHModelConfig> m_modelParameters; ///> SPH Model parameters (must be set before simulation)
std::shared_ptr<NeighborSearch> m_neighborSearcher; ///> Neighbor Search (must be initialized during model initialization)
+
+ std::shared_ptr<StdVectorOfVec3d> m_pressureAccels = nullptr;
+ std::shared_ptr<StdVectorOfVec3d> m_surfaceTensionAccels = nullptr;
};
} // end namespace imstk
diff --git a/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp b/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp
index 4c06dc9f5a2b0df9df9bf4474a7e259c10f30f5a..13c84c079c919daf9bc67d76a1d167534e89e937 100644
--- a/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp
+++ b/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp
@@ -25,36 +25,25 @@
namespace imstk
{
void
-PbdState::initialize(const size_t numNodes, const bool (& options)[3])
+PbdState::initialize(const size_t numNodes)
{
- // CHECKBACK : m_pos could actually not be another copy
- if (options[0])
- {
- m_pos.resize(numNodes, Vec3d(0, 0, 0));
- }
-
- if (options[1])
- {
- m_vel.resize(numNodes, Vec3d(0, 0, 0));
- }
-
- if (options[2])
- {
- m_acc.resize(numNodes, Vec3d(0, 0, 0));
- }
+ m_pos->resize(numNodes, Vec3d(0.0, 0.0, 0.0));
+ m_vel->resize(numNodes, Vec3d(0.0, 0.0, 0.0));
+ m_acc->resize(numNodes, Vec3d(0.0, 0.0, 0.0));
}
void
-PbdState::initialize(const std::shared_ptr<PointSet>& m, const bool (& options)[3])
+PbdState::initialize(const StdVectorOfVec3d& vertices)
{
- this->initialize(m->getNumVertices(), options);
+ this->initialize(vertices.size());
+ setPositions(vertices);
}
void
PbdState::setState(std::shared_ptr<PbdState> rhs)
{
- m_pos = rhs->getPositions();
- m_vel = rhs->getVelocities();
- m_acc = rhs->getAccelerations();
+ *m_pos = *rhs->getPositions();
+ *m_vel = *rhs->getVelocities();
+ *m_acc = *rhs->getAccelerations();
}
} // imstk
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectStates/imstkPbdState.h b/Source/DynamicalModels/ObjectStates/imstkPbdState.h
index bcc2884edf0417daf9cf185002813c922f46a4e4..65f5019e57a89f52c6b0097059741b7a5bdf3223 100644
--- a/Source/DynamicalModels/ObjectStates/imstkPbdState.h
+++ b/Source/DynamicalModels/ObjectStates/imstkPbdState.h
@@ -43,30 +43,40 @@ public:
///
/// \brief Initialize the pbd state
///
- void initialize(const size_t numNodes, const bool (& options)[3]);
- void initialize(const std::shared_ptr<PointSet>& m, const bool (& options)[3]);
+ void initialize(const size_t numNodes);
+ void initialize(const StdVectorOfVec3d& vertices);
///
/// \brief Get/Set nodal position given the index
///
- void setVertexPosition(const size_t& idx, const Vec3d& pos) { m_pos.at(idx) = pos; }
- Vec3d& getVertexPosition(const size_t& idx) { return m_pos.at(idx); }
+ void setVertexPosition(const size_t& idx, const Vec3d& pos) { m_pos->at(idx) = pos; }
+ Vec3d& getVertexPosition(const size_t& idx) { return m_pos->at(idx); }
///
/// \brief Returns the vector of current nodal positions
///
- StdVectorOfVec3d& getPositions() { return m_pos; }
- void setPositions(const StdVectorOfVec3d& p) { m_pos = p; }
+ std::shared_ptr<StdVectorOfVec3d> getPositions() { return m_pos; }
+ void setPositions(std::shared_ptr<StdVectorOfVec3d> p)
+ {
+ m_pos->resize(p->size());
+ std::copy(p->begin(), p->end(), m_pos->begin());
+ }
+
+ void setPositions(const StdVectorOfVec3d& p)
+ {
+ m_pos->resize(p.size());
+ std::copy(p.begin(), p.end(), m_pos->begin());
+ }
///
/// \brief Returns the vector of current nodal velocities
///
- StdVectorOfVec3d& getVelocities() { return m_vel; }
+ std::shared_ptr<StdVectorOfVec3d> getVelocities() { return m_vel; }
///
/// \brief Returns the vector of current nodal accelerations
///
- StdVectorOfVec3d& getAccelerations() { return m_acc; }
+ std::shared_ptr<StdVectorOfVec3d> getAccelerations() { return m_acc; }
///
/// \brief Set the state to a given one
@@ -74,8 +84,8 @@ public:
void setState(std::shared_ptr<PbdState> rhs);
private:
- StdVectorOfVec3d m_pos; ///> Nodal positions
- StdVectorOfVec3d m_vel; ///> Nodal velocities
- StdVectorOfVec3d m_acc; ///> Nodal acelerations
+ std::shared_ptr<StdVectorOfVec3d> m_pos = std::make_shared<StdVectorOfVec3d>(); ///> Nodal positions
+ std::shared_ptr<StdVectorOfVec3d> m_vel = std::make_shared<StdVectorOfVec3d>(); ///> Nodal velocities
+ std::shared_ptr<StdVectorOfVec3d> m_acc = std::make_shared<StdVectorOfVec3d>(); ///> Nodal acelerations
};
} // imstk
diff --git a/Source/DynamicalModels/ObjectStates/imstkRigidBodyState.h b/Source/DynamicalModels/ObjectStates/imstkRigidBodyState.h
index 2f0764a6d07e0a0016081a0b14f60138a38db3ea..b67ce30c98cce95d020b5a7b7080949d4f3efad8 100644
--- a/Source/DynamicalModels/ObjectStates/imstkRigidBodyState.h
+++ b/Source/DynamicalModels/ObjectStates/imstkRigidBodyState.h
@@ -19,8 +19,7 @@
=========================================================================*/
-#ifndef imstkRigidBodyState_h
-#define imstkRigidBodyState_h
+#pragma once
#include <Eigen/Dense>
#include <memory>
@@ -77,6 +76,4 @@ private:
//Can add linear velocity and angular velocity too.
};
-} // imstk
-
-#endif // imstkPbdState_h
\ No newline at end of file
+} // imstk
\ No newline at end of file
diff --git a/Source/SceneEntities/Objects/deprecated/imstkPbdRigidObject.h b/Source/SceneEntities/Objects/deprecated/imstkPbdRigidObject.h
deleted file mode 100644
index e9bfb4de1444facae73261e35b38e69f08dac32b..0000000000000000000000000000000000000000
--- a/Source/SceneEntities/Objects/deprecated/imstkPbdRigidObject.h
+++ /dev/null
@@ -1,67 +0,0 @@
-/*=========================================================================
-
- Library: iMSTK
-
- Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
- & Imaging in Medicine, Rensselaer Polytechnic Institute.
-
- Licensed under the Apache License, Version 2.0 (the "License");
- you may not use this file except in compliance with the License.
- You may obtain a copy of the License at
-
- http://www.apache.org/licenses/LICENSE-2.0.txt
-
- Unless required by applicable law or agreed to in writing, software
- distributed under the License is distributed on an "AS IS" BASIS,
- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- See the License for the specific language governing permissions and
- limitations under the License.
-
-=========================================================================*/
-
-#pragma once
-
-#include "imstkPbdObject.h"
-
-namespace imstk
-{
-class Geometry;
-class GeometryMap;
-
-///
-/// \class PbdRigidObject
-///
-/// \brief
-///
-class PbdRigidObject : public PbdObject
-{
-public:
- ///
- /// \brief Constructor
- ///
- explicit PbdRigidObject(const std::string& name) : PbdObject(name)
- { m_type = SceneObject::Type::Pbd; }
-
- ///
- /// \brief Destructor
- ///
- ~PbdRigidObject() = default;
-
- ///
- /// \brief
- ///
- virtual void integratePosition() { return; }
-
- ///
- /// \brief
- ///
- virtual void integrateVelocity() { return; }
-
- ///
- /// \brief
- ///
- virtual void updatePbdStates();
-
-protected:
-};
-} // imstk
diff --git a/Source/SceneEntities/Objects/imstkCollidingObject.cpp b/Source/SceneEntities/Objects/imstkCollidingObject.cpp
index cb07dc45a23d0f5e97954361081d3569e1aefc7a..324215144a0e40ad5da18280047b6fd0616cffaa 100644
--- a/Source/SceneEntities/Objects/imstkCollidingObject.cpp
+++ b/Source/SceneEntities/Objects/imstkCollidingObject.cpp
@@ -20,6 +20,7 @@
=========================================================================*/
#include "imstkCollidingObject.h"
+#include "imstkComputeGraph.h"
#include "imstkGeometryMap.h"
namespace imstk
@@ -27,19 +28,17 @@ namespace imstk
bool
CollidingObject::initialize()
{
- if (SceneObject::initialize())
+ if (!SceneObject::initialize())
{
- if (m_collidingToVisualMap)
- {
- m_collidingToVisualMap->initialize();
- }
-
- return true;
+ return false;
}
- else
+
+ if (m_collidingToVisualMap)
{
- return false;
+ m_collidingToVisualMap->initialize();
}
+
+ return true;
}
const std::shared_ptr<Geometry>&
diff --git a/Source/SceneEntities/Objects/imstkCollidingObject.h b/Source/SceneEntities/Objects/imstkCollidingObject.h
index dfe6e5e14c8c4fe4d11706088d08c6a8401a3417..512f1ab55de7aa72f44107a1909a17ba42b332d4 100644
--- a/Source/SceneEntities/Objects/imstkCollidingObject.h
+++ b/Source/SceneEntities/Objects/imstkCollidingObject.h
@@ -23,6 +23,7 @@
#include "imstkSceneObject.h"
#include "imstkMath.h"
+#include <unordered_map>
namespace imstk
{
@@ -81,7 +82,6 @@ public:
virtual bool initialize() override;
protected:
-
std::shared_ptr<Geometry> m_collidingGeometry; ///> Geometry for collisions
std::shared_ptr<GeometryMap> m_collidingToVisualMap; ///> Maps transformations to visual geometry
Vec3d m_force = Vec3d::Zero();
diff --git a/Source/SceneEntities/Objects/imstkDeformableObject.cpp b/Source/SceneEntities/Objects/imstkDeformableObject.cpp
index cd760dd8a30aaf2e6aac4662701484e53a028a3b..e4d1fb73f7da0a3df888484e3fa1bc0ce7d4ed8b 100644
--- a/Source/SceneEntities/Objects/imstkDeformableObject.cpp
+++ b/Source/SceneEntities/Objects/imstkDeformableObject.cpp
@@ -25,62 +25,19 @@
namespace imstk
{
-Vectord&
-FeDeformableObject::getContactForce()
-{
- CHECK(m_defModel != nullptr) << "deformation model pointer not valid DeformableObject::getContactForce()";
-
- return m_defModel->getContactForce();
-}
-
bool
FeDeformableObject::initialize()
{
- m_defModel = std::dynamic_pointer_cast<FEMDeformableBodyModel>(m_dynamicalModel);
- if (m_defModel)
- {
- return DynamicObject::initialize();
- }
- else
+ m_femModel = std::dynamic_pointer_cast<FEMDeformableBodyModel>(m_dynamicalModel);
+ if (m_femModel == nullptr)
{
LOG(FATAL) << "Dynamics pointer cast failure in DeformableObject::initialize()";
return false;
}
-}
-
-const Vectord&
-FeDeformableObject::getDisplacements() const
-{
- return m_defModel->getCurrentState()->getQ();
-}
-
-const Vectord&
-FeDeformableObject::getPrevDisplacements() const
-{
- return m_defModel->getPreviousState()->getQ();
-}
-const Vectord&
-FeDeformableObject::getVelocities() const
-{
- return m_defModel->getCurrentState()->getQDot();
-}
+ DynamicObject::initialize();
+ m_femModel->initialize();
-const Vectord&
-FeDeformableObject::getPrevVelocities() const
-{
- return m_defModel->getPreviousState()->getQDot();
-}
-
-const Vectord&
-FeDeformableObject::getAccelerations() const
-{
- return m_defModel->getCurrentState()->getQDotDot();
-}
-
-const Vectord&
-FeDeformableObject::getPrevAccelerations() const
-{
- return m_defModel->getPreviousState()->getQDotDot();
+ return true;
}
} // imstk
diff --git a/Source/SceneEntities/Objects/imstkDeformableObject.h b/Source/SceneEntities/Objects/imstkDeformableObject.h
index 00a76e40aed650cd094cae8dbfdaa49dd9c7feca..690075e4c75d81d177b25fd555daa37a8cefad51 100644
--- a/Source/SceneEntities/Objects/imstkDeformableObject.h
+++ b/Source/SceneEntities/Objects/imstkDeformableObject.h
@@ -30,8 +30,6 @@
namespace imstk
{
-class Geometry;
-class GeometryMap;
class FEMDeformableBodyModel;
///
@@ -42,70 +40,20 @@ class FEMDeformableBodyModel;
class FeDeformableObject : public DynamicObject
{
public:
-
- ///
- /// \brief Constructor
- ///
explicit FeDeformableObject(const std::string& name) : DynamicObject(name) { m_type = Type::FEMDeformable; }
FeDeformableObject() = delete;
- ///
- /// \brief Destructor
- ///
~FeDeformableObject() = default;
+public:
///
/// \brief Initialize the deformable object
///
bool initialize() override;
- ///
- /// \brief Initialize the kinematic state of the body
- ///
- void initializeState();
- void initializeState(const Vectord& p, const Vectord& v);
-
- ///
- /// \brief Set/Get dynamical model
- ///
- //void setDynamicalModel(std::shared_ptr<DynamicalModel<VectorizedState>> dynaDefModel) override;
-
- ///
- /// \brief Get the vector that holds the contact forces
- ///
- Vectord& getContactForce();
-
- ///
- /// \brief Get the vector of current displacements
- ///
- const Vectord& getDisplacements() const;
-
- ///
- /// \brief Get the vector of displacements from previous time step
- ///
- const Vectord& getPrevDisplacements() const;
-
- ///
- /// \brief Get the vector of current velocities
- ///
- const Vectord& getVelocities() const;
-
- ///
- /// \brief Get the vector of velocities from previous time step
- ///
- const Vectord& getPrevVelocities() const;
-
- ///
- /// \brief Get the vector of current accelerations
- ///
- const Vectord& getAccelerations() const;
-
- ///
- /// \brief Get the vector of accelerations from previous time step
- ///
- const Vectord& getPrevAccelerations() const;
+ std::shared_ptr<FEMDeformableBodyModel> getFEMModel() const { return m_femModel; }
protected:
- std::shared_ptr<FEMDeformableBodyModel> m_defModel;
+ std::shared_ptr<FEMDeformableBodyModel> m_femModel = nullptr;
};
} // imstk
diff --git a/Source/SceneEntities/Objects/imstkDynamicObject.cpp b/Source/SceneEntities/Objects/imstkDynamicObject.cpp
index 50d28d7b2ed4f744a3bcbd752cc506f5242af72a..b56b3ea93cd609e4b0d5046db0c06e65f6cfc341 100644
--- a/Source/SceneEntities/Objects/imstkDynamicObject.cpp
+++ b/Source/SceneEntities/Objects/imstkDynamicObject.cpp
@@ -22,6 +22,7 @@
#include "imstkDynamicObject.h"
#include "imstkGeometryMap.h"
#include "imstkAbstractDynamicalModel.h"
+#include "imstkComputeGraph.h"
#include "imstkLogger.h"
namespace imstk
@@ -61,6 +62,17 @@ DynamicObject::updateGeometries()
}
}
+void
+DynamicObject::updatePhysicsGeometry()
+{
+ m_dynamicalModel->updatePhysicsGeometry();
+
+ if (m_physicsToCollidingGeomMap)
+ {
+ m_physicsToCollidingGeomMap->apply();
+ }
+}
+
bool
DynamicObject::initialize()
{
@@ -84,6 +96,16 @@ DynamicObject::initialize()
}
}
+void
+DynamicObject::initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
+{
+ // Copy, sum, and connect the model graph to nest within this graph
+ m_computeGraph->addEdge(source, getUpdateNode());
+ m_dynamicalModel->initGraphEdges();
+ m_computeGraph->nestGraph(m_dynamicalModel->getComputeGraph(), getUpdateNode(), getUpdateGeometryNode());
+ m_computeGraph->addEdge(getUpdateGeometryNode(), sink);
+}
+
void
DynamicObject::reset()
{
diff --git a/Source/SceneEntities/Objects/imstkDynamicObject.h b/Source/SceneEntities/Objects/imstkDynamicObject.h
index 0535d059ab62f7389f3b3d3dda621778552d056b..c562afdb12555250fb776a2854d3ab1862eb6f1a 100644
--- a/Source/SceneEntities/Objects/imstkDynamicObject.h
+++ b/Source/SceneEntities/Objects/imstkDynamicObject.h
@@ -37,7 +37,6 @@ class AbstractDynamicalModel;
class DynamicObject : public CollidingObject
{
public:
-
///
/// \brief Destructor
///
@@ -82,16 +81,27 @@ public:
///
void updateGeometries() final;
+ ///
+ /// \brief Update only the physics geometry and apply collision map
+ ///
+ void updatePhysicsGeometry();
+
///
/// \brief Initialize the scene object
///
virtual bool initialize() override;
///
- /// \brief Reset the dynamic object to its initial state
+ /// \brief Reset the dynamic object by reseting the respective DynamicalModel and Geometry
///
virtual void reset() override;
+protected:
+ ///
+ /// \brief Setup connectivity of compute graph
+ ///
+ virtual void initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink) override;
+
protected:
///
@@ -99,12 +109,12 @@ protected:
///
explicit DynamicObject(const std::string& name) : CollidingObject(name) {}
- std::shared_ptr<AbstractDynamicalModel> m_dynamicalModel; ///> Dynamical model
- std::shared_ptr<Geometry> m_physicsGeometry; ///> Geometry used for Physics
+ std::shared_ptr<AbstractDynamicalModel> m_dynamicalModel = nullptr; ///> Dynamical model
+ std::shared_ptr<Geometry> m_physicsGeometry = nullptr; ///> Geometry used for Physics
- //Maps
- std::shared_ptr<GeometryMap> m_physicsToCollidingGeomMap; ///> Maps from Physics to collision geometry
- std::shared_ptr<GeometryMap> m_physicsToVisualGeomMap; ///> Maps from Physics to visual geometry
+ // Maps
+ std::shared_ptr<GeometryMap> m_physicsToCollidingGeomMap = nullptr; ///> Maps from Physics to collision geometry
+ std::shared_ptr<GeometryMap> m_physicsToVisualGeomMap = nullptr; ///> Maps from Physics to visual geometry
bool m_updateVisualFromPhysicsGeometry = true; ///> Defines if visual is updated from colliding mapping or physics mapping
};
} // imstk
diff --git a/Source/SceneEntities/Objects/imstkPbdObject.cpp b/Source/SceneEntities/Objects/imstkPbdObject.cpp
index ae380e2f71d42a60afbbf080900b0f7d70ad8782..731aa243662478c8d175300cbc08b51910aef344 100644
--- a/Source/SceneEntities/Objects/imstkPbdObject.cpp
+++ b/Source/SceneEntities/Objects/imstkPbdObject.cpp
@@ -20,53 +20,31 @@
=========================================================================*/
#include "imstkPbdObject.h"
-#include "imstkGeometryMap.h"
-#include "imstkPbdModel.h"
#include "imstkLogger.h"
+#include "imstkPbdModel.h"
namespace imstk
{
-bool
-PbdObject::initialize()
+std::shared_ptr<PbdModel>
+PbdObject::getPbdModel()
{
m_pbdModel = std::dynamic_pointer_cast<PbdModel>(m_dynamicalModel);
+ return m_pbdModel;
+};
- CHECK(m_pbdModel != nullptr) << "Dynamics pointer cast failure in PbdObject::initialize()";
-
- return DynamicObject::initialize();
-}
-
-void
-PbdObject::integratePosition()
-{
- if (m_pbdModel && m_pbdModel->hasConstraints())
- {
- m_pbdModel->integratePosition();
- }
-}
-
-void
-PbdObject::updateVelocity()
+bool
+PbdObject::initialize()
{
- if (m_pbdModel && m_pbdModel->hasConstraints())
+ m_pbdModel = std::dynamic_pointer_cast<PbdModel>(m_dynamicalModel);
+ if (m_pbdModel == nullptr)
{
- m_pbdModel->updateVelocity();
+ LOG(FATAL) << "Dynamics pointer cast failure in PbdObject::initialize()";
+ return false;
}
-}
-void
-PbdObject::solveConstraints()
-{
- if (m_pbdModel && m_pbdModel->hasConstraints())
- {
- m_pbdModel->projectConstraints();
- }
-}
+ DynamicObject::initialize();
+ m_pbdModel->initialize();
-void
-PbdObject::reset()
-{
- DynamicObject::reset();
- this->updateVelocity();
+ return true;
}
} //imstk
diff --git a/Source/SceneEntities/Objects/imstkPbdObject.h b/Source/SceneEntities/Objects/imstkPbdObject.h
index 17feda4e8ffffc9326316bba6045e9ca91b316bd..0ee18a27ce575a4c5a19b7b5e4d9139fa2964722 100644
--- a/Source/SceneEntities/Objects/imstkPbdObject.h
+++ b/Source/SceneEntities/Objects/imstkPbdObject.h
@@ -25,22 +25,17 @@
namespace imstk
{
-class Geometry;
-class GeometryMap;
class PbdModel;
///
/// \class PbdObject
///
/// \brief Base class for scene objects that move and/or deform under position
-/// based dynamics formulation
+/// based dynamics formulation, implements the PbdModel and PbdSolver
///
class PbdObject : public DynamicObject
{
public:
- ///
- /// \brief Constructor
- ///
explicit PbdObject(const std::string& name) : DynamicObject(name)
{
m_type = SceneObject::Type::Pbd;
@@ -48,37 +43,20 @@ public:
PbdObject() = delete;
- ///
- /// \brief Destructor
- ///
~PbdObject() override = default;
+public:
///
- /// \brief Initialize the pbd scene object
- ///
- bool initialize() override;
-
- ///
- /// \brief Update the position based on Verlet time stepping rule
- ///
- virtual void integratePosition();
-
- ///
- /// \brief Update the velocity
- ///
- virtual void updateVelocity();
-
- ///
- /// \brief Solve the pbd constraints by projection
+ /// \biref Get the Pbd model of the object
///
- virtual void solveConstraints();
+ std::shared_ptr<PbdModel> getPbdModel();
///
- /// \brief Reset the PBD object to its initial state
+ /// \brief Initialize the pbd scene object
///
- void reset() override;
+ bool initialize() override;
protected:
- std::shared_ptr<PbdModel> m_pbdModel; ///> PBD mathematical model
+ std::shared_ptr<PbdModel> m_pbdModel = nullptr; ///> PBD mathematical model
};
} // imstk
diff --git a/Source/SceneEntities/Objects/imstkRigidObject.cpp b/Source/SceneEntities/Objects/imstkRigidObject.cpp
index 407520854b7fb416a8c924a1d2b758b6a99bd8c6..fa2ebba781ae404919be4cf353ecb80b111c0cd0 100644
--- a/Source/SceneEntities/Objects/imstkRigidObject.cpp
+++ b/Source/SceneEntities/Objects/imstkRigidObject.cpp
@@ -18,28 +18,32 @@
limitations under the License.
=========================================================================*/
+
#include "imstkRigidObject.h"
#include "imstkLogger.h"
+#include "imstkRigidBodyModel.h"
namespace imstk
{
bool
RigidObject::initialize()
{
- auto m_rigidBodyModel = std::dynamic_pointer_cast<RigidBodyModel>(m_dynamicalModel);
- if (m_rigidBodyModel)
- {
- /*if (!m_rigidBodyModel->getRigidBodyWorld())
- {
- LOG(WARNING) << "RigidObject::initialize() - Rigid body world not specified!";
- return false;
- }*/
- return DynamicObject::initialize();
- }
- else
+ m_rigidBodyModel = std::dynamic_pointer_cast<RigidBodyModel>(m_dynamicalModel);
+ if (m_rigidBodyModel == nullptr)
{
- LOG(WARNING) << "RigidObject::initialize() - Dynamics pointer cast failure";
+ LOG(FATAL) << "Dynamics pointer cast failure in RigidObject::initialize()";
return false;
}
+
+ DynamicObject::initialize();
+ m_rigidBodyModel->initialize();
+
+ return true;
+}
+
+void
+RigidObject::addForce(const Vec3d& force, const Vec3d& pos, bool wakeup)
+{
+ m_rigidBodyModel->addForce(force, pos, wakeup);
}
} // imstk
diff --git a/Source/SceneEntities/Objects/imstkRigidObject.h b/Source/SceneEntities/Objects/imstkRigidObject.h
index 7e0ff3a2810f87506bc782d6108e36521b841abb..eb24b27f24ab84c3a3ec4b9f9703d5ffd29efca8 100644
--- a/Source/SceneEntities/Objects/imstkRigidObject.h
+++ b/Source/SceneEntities/Objects/imstkRigidObject.h
@@ -21,12 +21,12 @@
#pragma once
-// imstk
#include "imstkDynamicObject.h"
-#include "imstkRigidBodyModel.h"
namespace imstk
{
+class RigidBodyModel;
+
///
/// \class RigidObject
///
@@ -35,20 +35,14 @@ namespace imstk
class RigidObject : public DynamicObject
{
public:
-
- ///
- /// \brief Constructor
- ///
explicit RigidObject(const std::string& name) : DynamicObject(name)
{
m_type = Type::Rigid;
}
- ///
- /// \brief Destructor
- ///
~RigidObject() = default;
+public:
///
/// \brief Initialize the rigid scene object
///
@@ -57,24 +51,11 @@ public:
///
/// \brief Add local force at a position relative to object
///
- void addForce(const Vec3d& force, const Vec3d& pos, bool wakeup = true)
- {
- getRigidBodyModel()->addForce(force, pos, wakeup);
- }
+ void addForce(const Vec3d& force, const Vec3d& pos, bool wakeup = true);
- ///
- /// \brief Get/Set rigid body model
- ///
- void setRigidBodyModel(std::shared_ptr<RigidBodyModel> rbModel) { m_dynamicalModel = rbModel; };
- std::shared_ptr<RigidBodyModel> getRigidBodyModel() const { return std::dynamic_pointer_cast<RigidBodyModel>(m_dynamicalModel); };
-
- void reset() override
- {
- m_dynamicalModel->resetToInitialState();
- this->updateGeometries();
- }
+ std::shared_ptr<RigidBodyModel> getRigidBodyModel() const { return m_rigidBodyModel; }
protected:
- //std::shared_ptr<RigidBodyModel> m_rigidBodyModel; ///> PBD mathematical model
+ std::shared_ptr<RigidBodyModel> m_rigidBodyModel;
};
} // imstk
\ No newline at end of file
diff --git a/Source/SceneEntities/Objects/imstkSPHObject.cpp b/Source/SceneEntities/Objects/imstkSPHObject.cpp
index f5896db40469425e42833a82ca5ab554a6c1645a..9be7701c33be9051e53d9ece65162713da930191 100644
--- a/Source/SceneEntities/Objects/imstkSPHObject.cpp
+++ b/Source/SceneEntities/Objects/imstkSPHObject.cpp
@@ -26,15 +26,22 @@ namespace imstk
{
SPHObject::SPHObject(const std::string& name) : DynamicObject(name)
{
- this->m_type = SceneObject::Type::SPH;
+ this->m_type = Type::SPH;
}
bool
SPHObject::initialize()
{
- m_SPHModel = std::dynamic_pointer_cast<SPHModel>(this->m_dynamicalModel);
- LOG_IF(FATAL, (!m_SPHModel)) << "Invalid SPH model";
+ m_SPHModel = std::dynamic_pointer_cast<SPHModel>(m_dynamicalModel);
+ if (m_SPHModel == nullptr)
+ {
+ LOG(FATAL) << "Dynamics pointer cast failure in SPHObject::initialize()";
+ return false;
+ }
- return m_SPHModel->initialize();
+ DynamicObject::initialize();
+ m_SPHModel->initialize();
+
+ return true;
}
} // end namespace imstk
diff --git a/Source/SceneEntities/Objects/imstkSPHObject.h b/Source/SceneEntities/Objects/imstkSPHObject.h
index c123cadb9c953c888d27de7b2e78cabf2ec9e593..d6f8936633e1e84e0043e18ae2c87b4cdc604987 100644
--- a/Source/SceneEntities/Objects/imstkSPHObject.h
+++ b/Source/SceneEntities/Objects/imstkSPHObject.h
@@ -30,22 +30,17 @@ class SPHModel;
///
/// \class SPHObject
///
-/// \brief Base class for scene objects that move and/or deform under position
-/// based dynamics formulation
+/// \brief Base class for scene objects that move and/or deform under
+/// smooth particle hydrodynamics
///
class SPHObject : public DynamicObject
{
public:
- ///
- /// \brief Constructor
- ///
explicit SPHObject(const std::string& name);
- ///
- /// \brief Destructor
- ///
virtual ~SPHObject() override = default;
+public:
///
/// \brief Initialize the SPH scene object
///
@@ -54,9 +49,9 @@ public:
///
/// \brief Get the SPH model of the object
///
- const std::shared_ptr<SPHModel>& getSPHModel() const { assert(m_SPHModel); return m_SPHModel; }
+ std::shared_ptr<SPHModel> getSPHModel() const { assert(m_SPHModel); return m_SPHModel; }
protected:
- std::shared_ptr<SPHModel> m_SPHModel;
+ std::shared_ptr<SPHModel> m_SPHModel = nullptr;
};
} // end namespace imstk
diff --git a/Source/SceneEntities/Objects/imstkSceneObject.cpp b/Source/SceneEntities/Objects/imstkSceneObject.cpp
index 65c82cc3cbef5701b0c6e0177f48a571c0d4f9ba..b28d296c354dfdc74ede2589e18f4bb2dbf8b1e2 100644
--- a/Source/SceneEntities/Objects/imstkSceneObject.cpp
+++ b/Source/SceneEntities/Objects/imstkSceneObject.cpp
@@ -20,12 +20,17 @@
=========================================================================*/
#include "imstkSceneObject.h"
+#include "imstkComputeGraph.h"
+#include "imstkComputeNode.h"
namespace imstk
{
-SceneObject::SceneObject(const std::string& name) : m_name(name), SceneEntity()
+SceneObject::SceneObject(const std::string& name) :
+ m_type(Type::Visual), m_name(name),
+ m_computeGraph(std::make_shared<ComputeGraph>("SceneObject_" + name + "_Source", "SceneObject_" + name + "_Sink")), SceneEntity()
{
- m_type = Type::Visual;
+ m_updateNode = m_computeGraph->addFunction("SceneObject_" + name + "_Update", [&]() { this->update(); });
+ m_updateGeometryNode = m_computeGraph->addFunction("SceneObject_" + name + "_UpdateGeometry", [&]() { this->updateGeometries(); });
}
std::shared_ptr<Geometry>
@@ -51,57 +56,18 @@ SceneObject::setVisualGeometry(std::shared_ptr<Geometry> geometry)
}
}
-std::shared_ptr<Geometry>
-SceneObject::getMasterGeometry() const
-{
- return this->getVisualGeometry();
-}
-
-std::shared_ptr<VisualModel>
-SceneObject::getVisualModel(unsigned int index)
-{
- return m_visualModels[index];
-}
-
void
-SceneObject::addVisualModel(std::shared_ptr<VisualModel> visualModel)
-{
- m_visualModels.push_back(visualModel);
-}
-
-const std::vector<std::shared_ptr<VisualModel>>&
-SceneObject::getVisualModels()
-{
- return m_visualModels;
-}
-
-size_t
-SceneObject::getNumVisualModels()
-{
- return m_visualModels.size();
-}
-
-const SceneObject::Type&
-SceneObject::getType() const
-{
- return m_type;
-}
-
-void
-SceneObject::setType(SceneObject::Type type)
-{
- m_type = type;
-}
-
-const std::string&
-SceneObject::getName() const
+SceneObject::initGraphEdges()
{
- return m_name;
+ m_computeGraph->clearEdges();
+ initGraphEdges(m_computeGraph->getSource(), m_computeGraph->getSink());
}
void
-SceneObject::setName(const std::string& name)
+SceneObject::initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink)
{
- m_name = name;
+ m_computeGraph->addEdge(source, m_updateNode);
+ m_computeGraph->addEdge(m_updateNode, m_updateGeometryNode);
+ m_computeGraph->addEdge(m_updateGeometryNode, sink);
}
} // imstk
diff --git a/Source/SceneEntities/Objects/imstkSceneObject.h b/Source/SceneEntities/Objects/imstkSceneObject.h
index 6d1cb9eaeca15f764fbb7e01d08be180a6794ea9..df93a23dc23d3ccd2b9111a6714bab13825358d4 100644
--- a/Source/SceneEntities/Objects/imstkSceneObject.h
+++ b/Source/SceneEntities/Objects/imstkSceneObject.h
@@ -29,8 +29,10 @@
namespace imstk
{
-class Geometry;
+class ComputeNode;
+class ComputeGraph;
class DeviceClient;
+class Geometry;
///
/// \class SceneObject
@@ -62,16 +64,22 @@ public:
///
virtual ~SceneObject() = default;
+public:
///
/// \brief Get the type of the object
///
- const Type& getType() const;
+ const Type& getType() const { return m_type; }
+
+ ///
+ /// \brief Get the computational graph
+ ///
+ std::shared_ptr<ComputeGraph> getComputeGraph() const { return m_computeGraph; }
///
/// \brief Get/Set the custom name of the scene object
///
- const std::string& getName() const;
- void setName(const std::string& name);
+ const std::string& getName() const { return m_name; }
+ void setName(const std::string& name) { m_name = name; }
///
/// \brief DEPRECATED: Get/Set geometry used for viewing
@@ -84,23 +92,38 @@ public:
///
/// \brief Get/add visual model
///
- std::shared_ptr<VisualModel> getVisualModel(unsigned int index);
- void addVisualModel(std::shared_ptr<VisualModel> visualModel);
+ std::shared_ptr<VisualModel> getVisualModel(unsigned int index) { return m_visualModels[index]; }
+ void addVisualModel(std::shared_ptr<VisualModel> visualModel) { m_visualModels.push_back(visualModel); }
///
/// \brief Get all visual models
///
- const std::vector<std::shared_ptr<VisualModel>>& getVisualModels();
+ const std::vector<std::shared_ptr<VisualModel>>& getVisualModels() { return m_visualModels; }
///
/// \brief Get number of visual models
///
- size_t getNumVisualModels();
+ size_t getNumVisualModels() { return m_visualModels.size(); }
///
/// \brief Get the master geometry
///
- virtual std::shared_ptr<Geometry> getMasterGeometry() const;
+ virtual std::shared_ptr<Geometry> getMasterGeometry() const { return this->getVisualGeometry(); }
+
+ ///
+ /// \brief Returns the computational node for updating
+ ///
+ std::shared_ptr<ComputeNode> getUpdateNode() const { return m_updateNode; }
+
+ ///
+ /// \brief Returns the computational node for updating geometry
+ ///
+ std::shared_ptr<ComputeNode> getUpdateGeometryNode() const { return m_updateGeometryNode; }
+
+ ///
+ /// \brief
+ ///
+ virtual void update() { }
///
/// \brief
@@ -112,20 +135,32 @@ public:
///
virtual bool initialize() { return true; }
+ ///
+ /// \brief Initializes the edges of the SceneObject's computational graph
+ ///
+ void initGraphEdges();
+
///
/// \brief
///
- virtual void reset() {}
+ virtual void reset() { }
protected:
///
- /// \brief Assigns the type of the object
+ /// \brief Setup connectivity of the compute graph
///
- void setType(Type type);
+ virtual void initGraphEdges(std::shared_ptr<ComputeNode> source, std::shared_ptr<ComputeNode> sink);
+protected:
Type m_type; ///> Type of the scene object
std::string m_name; ///> Custom name of the scene object
std::vector<std::shared_ptr<VisualModel>> m_visualModels; ///> Visual objects for rendering
+ std::shared_ptr<ComputeGraph> m_computeGraph = nullptr; ///> Computational Graph
+
+private:
+ // Dissallow reassignment of these in subclasses
+ std::shared_ptr<ComputeNode> m_updateNode = nullptr;
+ std::shared_ptr<ComputeNode> m_updateGeometryNode = nullptr;
};
using VisualObject = SceneObject;
diff --git a/Source/Solvers/CMakeLists.txt b/Source/Solvers/CMakeLists.txt
index 14b17c8a3dd2ecff8eb109f6dc8d1c4efdcd8a50..9f45108fb8cb52421166d8cfc28e347a680b2927 100644
--- a/Source/Solvers/CMakeLists.txt
+++ b/Source/Solvers/CMakeLists.txt
@@ -5,7 +5,6 @@ include(imstkAddLibrary)
imstk_add_library( Solvers
DEPENDS
Constraints
- SceneEntities
)
#-----------------------------------------------------------------------------
diff --git a/Source/Solvers/imstkPbdSolver.cpp b/Source/Solvers/imstkPbdSolver.cpp
index 83cbc529fbc3b47bff087936aacd81d2a1bd4bd2..16c5eb09e2b8a89a3980aaf54281bc525a23aa3d 100644
--- a/Source/Solvers/imstkPbdSolver.cpp
+++ b/Source/Solvers/imstkPbdSolver.cpp
@@ -19,39 +19,100 @@
=========================================================================*/
-#include "imstkPbdObject.h"
#include "imstkPbdSolver.h"
-#include "imstkPbdCollisionConstraint.h"
#include "imstkParallelUtils.h"
namespace imstk
{
+PbdSolver::PbdSolver() :
+ m_partitionedConstraints(std::make_shared<std::vector<PBDConstraintVector>>()),
+ m_constraints(std::make_shared<PBDConstraintVector>()),
+ m_positions(std::make_shared<StdVectorOfVec3d>()),
+ m_invMasses(std::make_shared<StdVectorOfReal>())
+{
+}
+
void
PbdSolver::solve()
{
- m_pbdObject->integratePosition();
- m_pbdObject->solveConstraints();
- resolveCollisionConstraints();
- m_pbdObject->updateVelocity();
+ // Solve the constraints and partitioned constraints
+ StdVectorOfVec3d& currPositions = *m_positions;
+ const StdVectorOfReal& invMasses = *m_invMasses;
+
+ const PBDConstraintVector& constraints = *m_constraints;
+ const std::vector<PBDConstraintVector>& partitionedConstraints = *m_partitionedConstraints;
+
+ unsigned int i = 0;
+ while (i++ < m_iterations)
+ {
+ for (size_t j = 0; j < constraints.size(); j++)
+ {
+ constraints[j]->solvePositionConstraint(currPositions, invMasses);
+ }
+
+ for (size_t j = 0; j < partitionedConstraints.size(); j++)
+ {
+ const PBDConstraintVector& constraintPartition = partitionedConstraints[j];
+ const size_t numConstraintsInPartition = constraintPartition.size();
+
+ ParallelUtils::parallelFor(constraintPartition.size(),
+ [&](const size_t idx)
+ {
+ constraintPartition[idx]->solvePositionConstraint(currPositions, invMasses);
+ });
+ // Sequential
+ /*for (size_t k = 0; k < constraintPartition.size(); k++)
+ {
+ constraintPartition[k]->solvePositionConstraint(currPositions, invMasses);
+ }*/
+ }
+ }
+}
+
+
+PbdCollisionSolver::PbdCollisionSolver() :
+ m_collisionConstraints(std::make_shared<std::list<PBDCollisionConstraintVector*>>()),
+ m_collisionConstraintsData(std::make_shared<std::list<CollisionConstraintData>>())
+{
+}
+
+void
+PbdCollisionSolver::addCollisionConstraints(PBDCollisionConstraintVector* constraints,
+ std::shared_ptr<StdVectorOfVec3d> posA, std::shared_ptr<StdVectorOfReal> invMassA,
+ std::shared_ptr<StdVectorOfVec3d> posB, std::shared_ptr<StdVectorOfReal> invMassB)
+{
+ m_collisionConstraints->push_back(constraints);
+ m_collisionConstraintsData->push_back({ posA, invMassA, posB, invMassB });
}
void
-PbdSolver::resolveCollisionConstraints()
+PbdCollisionSolver::solve()
{
- if (m_PBDConstraints.size() > 0)
+ // Solve collision constraints
+ if (m_collisionConstraints->size() > 0)
{
- uint32_t maxIter = 3u;
- uint32_t i = 0;
- while (++i < maxIter)
+ unsigned int i = 0;
+ while (i++ < m_collisionIterations)
{
- for (const auto constraintList : m_PBDConstraints)
+ std::list<CollisionConstraintData>::iterator colDataIter = m_collisionConstraintsData->begin();
+ for (auto constraintList : *m_collisionConstraints)
{
- for (size_t k = 0; k < constraintList->size(); ++k)
+ CollisionConstraintData colData = *colDataIter;
+ StdVectorOfVec3d& posA = *colData.m_posA;
+ const StdVectorOfReal& invMassA = *colData.m_invMassA;
+ StdVectorOfVec3d& posB = *colData.m_posB;
+ const StdVectorOfReal& invMassB = *colData.m_invMassB;
+ const PBDCollisionConstraintVector& constraints = *constraintList;
+ for (size_t j = 0; j < constraints.size(); j++)
{
- (*constraintList)[k]->solvePositionConstraint();
+ constraints[j]->solvePositionConstraint(posA, posB, invMassA, invMassB);
}
+ colDataIter++;
}
}
+ m_collisionConstraints->clear();
+ m_collisionConstraintsData->clear();
}
}
+
} // end namespace imstk
diff --git a/Source/Solvers/imstkPbdSolver.h b/Source/Solvers/imstkPbdSolver.h
index ef20dbdf5ce286c95b4bdf42810f46a5de65713a..87c3cd9124286aec3c691d3a0e4e42cea77222f8 100644
--- a/Source/Solvers/imstkPbdSolver.h
+++ b/Source/Solvers/imstkPbdSolver.h
@@ -21,67 +21,120 @@
#pragma once
+#include "imstkPbdCollisionConstraint.h"
+#include "imstkPbdConstraint.h"
#include "imstkSolverBase.h"
-#include "imstkGraph.h"
-
-#include <unordered_set>
namespace imstk
{
-class PbdObject;
-class PbdCollisionConstraint;
+struct CollisionConstraintData
+{
+ CollisionConstraintData(std::shared_ptr<StdVectorOfVec3d> posA,
+ std::shared_ptr<StdVectorOfReal> invMassA,
+ std::shared_ptr<StdVectorOfVec3d> posB,
+ std::shared_ptr<StdVectorOfReal> invMassB) :
+ m_posA(posA), m_invMassA(invMassA), m_posB(posB), m_invMassB(invMassB)
+ {
+ }
+
+ std::shared_ptr<StdVectorOfVec3d> m_posA = nullptr;
+ std::shared_ptr<StdVectorOfReal> m_invMassA = nullptr;
+ std::shared_ptr<StdVectorOfVec3d> m_posB = nullptr;
+ std::shared_ptr<StdVectorOfReal> m_invMassB = nullptr;
+};
+
///
/// \class PbdSolver
///
/// \brief Position Based Dynamics solver
+/// This solver can operate on both partitioned constraints (unordered_set of vector'd constraints) in parallel
+/// and sequentially on vector'd constraints. It requires a set of constraints, positions, and invMasses.
///
class PbdSolver : public SolverBase
{
-using PBDConstraintVector = std::vector<PbdCollisionConstraint*>;
-
public:
///
/// \brief Constructors/Destructor
///
- PbdSolver() = default;
+ PbdSolver();
virtual ~PbdSolver() override = default;
- PbdSolver(const PbdSolver& other) = delete;
- PbdSolver& operator=(const PbdSolver& other) = delete;
+public:
+ ///
+ /// \brief Set Iterations. The number of nonlinear iterations.
+ ///
+ void setIterations(const size_t iterations) { this->m_iterations = iterations; }
///
- /// \brief Set MaxIterations. The maximum number of nonlinear iterations.
+ /// \brief Sets the partioned constraints the solver should solve for
+ /// These will be solved in parallel
///
- void setMaxIterations(const size_t newMaxIterations) { this->m_maxIterations = newMaxIterations; }
+ void setPartitionedConstraints(std::shared_ptr<std::vector<PBDConstraintVector>> partitionedConstraints) { this->m_partitionedConstraints = partitionedConstraints; }
///
- /// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
+ /// \brief Sets the constraints the solver should solve for
+ /// These wil be solved sequentially
///
- size_t getMaxIterations() const { return this->m_maxIterations; }
+ void setConstraints(std::shared_ptr<PBDConstraintVector> constraints) { this->m_constraints = constraints; }
///
- /// \brief
+ /// \brief Sets the positions the solver should solve with
///
- void setPbdObject(const std::shared_ptr<PbdObject>& pbdObj) { m_pbdObject = pbdObj; }
+ void setPositions(std::shared_ptr<StdVectorOfVec3d> positions) { this->m_positions = positions; }
+
+ ///
+ /// \brief Sets the invMasses the solver should solve with
+ ///
+ void setInvMasses(std::shared_ptr<StdVectorOfReal> invMasses) { this->m_invMasses = invMasses; }
+
+ ///
+ /// \brief Get Iterations. Returns current nonlinear iterations.
+ ///
+ size_t getIterations() const { return this->m_iterations; }
///
/// \brief Solve the non linear system of equations G(x)=0 using Newton's method.
///
void solve() override;
+private:
+ size_t m_iterations = 20; ///> Number of NL Gauss-Seidel iterations for regular constraints
+
+ std::shared_ptr<std::vector<PBDConstraintVector>> m_partitionedConstraints = nullptr; ///> Set of vector'd/partitioned pbd constraints
+ std::shared_ptr<PBDConstraintVector> m_constraints = nullptr; ///> Vector of constraints
+
+ std::shared_ptr<StdVectorOfVec3d> m_positions = nullptr;
+ std::shared_ptr<StdVectorOfReal> m_invMasses = nullptr;
+};
+
+class PbdCollisionSolver : SolverBase
+{
+public:
+ PbdCollisionSolver();
+ virtual ~PbdCollisionSolver() override = default;
+
+public:
+ ///
+ /// \brief Get CollisionIterations
+ ///
+ size_t getCollisionIterations() const { return this->m_collisionIterations; }
+
///
/// \brief Add the global collision contraints to this solver
///
- void addCollisionConstraints(PBDConstraintVector* constraints) { m_PBDConstraints.insert(constraints); }
+ void addCollisionConstraints(PBDCollisionConstraintVector* constraints,
+ std::shared_ptr<StdVectorOfVec3d> posA, std::shared_ptr<StdVectorOfReal> invMassA,
+ std::shared_ptr<StdVectorOfVec3d> posB, std::shared_ptr<StdVectorOfReal> invMassB);
///
- /// \brief Solve the global collision contraints charged to this solver
+ /// \brief Solve the non linear system of equations G(x)=0 using Newton's method.
///
- void resolveCollisionConstraints();
+ void solve() override;
private:
- size_t m_maxIterations = 20; ///< Maximum number of NL Gauss-Seidel iterations
- std::unordered_set<PBDConstraintVector*> m_PBDConstraints; ///< Collision contraints charged to this solver
- std::shared_ptr<PbdObject> m_pbdObject;
+ size_t m_collisionIterations = 5; ///> Number of NL Gauss-Seidel iterations for collision constraints
+
+ std::shared_ptr<std::list<PBDCollisionConstraintVector*>> m_collisionConstraints = nullptr; ///< Collision contraints charged to this solver
+ std::shared_ptr<std::list<CollisionConstraintData>> m_collisionConstraintsData = nullptr;
};
} // imstk
diff --git a/Source/Solvers/imstkSPHSolver.h b/Source/Solvers/imstkSPHSolver.h
deleted file mode 100644
index 58f08edf73e468683864cbaebae71588208f4ae7..0000000000000000000000000000000000000000
--- a/Source/Solvers/imstkSPHSolver.h
+++ /dev/null
@@ -1,62 +0,0 @@
-/*=========================================================================
-
- Library: iMSTK
-
- Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
- & Imaging in Medicine, Rensselaer Polytechnic Institute.
-
- Licensed under the Apache License, Version 2.0 (the "License");
- you may not use this file except in compliance with the License.
- You may obtain a copy of the License at
-
- http://www.apache.org/licenses/LICENSE-2.0.txt
-
- Unless required by applicable law or agreed to in writing, software
- distributed under the License is distributed on an "AS IS" BASIS,
- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- See the License for the specific language governing permissions and
- limitations under the License.
-
-=========================================================================*/
-
-#pragma once
-
-#include "imstkSolverBase.h"
-#include "imstkSPHObject.h"
-#include "imstkSPHModel.h"
-
-namespace imstk
-{
-///
-/// \class SPHSolver
-/// \brief SPH solver
-///
-class SPHSolver : public SolverBase
-{
-public:
- SPHSolver() = default;
- virtual ~SPHSolver() override = default;
-
- SPHSolver(const SPHSolver& other) = delete;
- SPHSolver& operator=(const SPHSolver& other) = delete;
-
- ///
- /// \brief Set the simulation object
- ///
- void setSPHObject(const std::shared_ptr<SPHObject>& obj)
- { assert(obj); m_SPHObject = obj; }
-
- ///
- /// \brief Advance one time step
- ///
- virtual void solve() override
- {
- const auto& SPHModel = m_SPHObject->getSPHModel();
- LOG_IF(FATAL, (!SPHModel)) << "SPH model has not been initialized";
- SPHModel->advanceTimeStep();
- }
-
-private:
- std::shared_ptr<SPHObject> m_SPHObject;
-};
-} // end namespace imstk
diff --git a/Source/Solvers/imstkSolverBase.h b/Source/Solvers/imstkSolverBase.h
index 4db160b9ad3e75938bad0a36a5e7c82219b21992..6249e0b0c412196e331a2b49fc52188f4d398343 100644
--- a/Source/Solvers/imstkSolverBase.h
+++ b/Source/Solvers/imstkSolverBase.h
@@ -28,11 +28,10 @@ namespace imstk
///
class SolverBase
{
-public:
- ///
- /// \brief Default Constructor/Destructor
- ///
+protected:
SolverBase() = default;
+
+public:
virtual ~SolverBase() = default;
virtual void solve() = 0;