diff --git a/Examples/PBD/PBDInjection/CMakeLists.txt b/Examples/PBD/PBDInjection/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6c514bbce02afa50b945347831de47c94b4045c0
--- /dev/null
+++ b/Examples/PBD/PBDInjection/CMakeLists.txt
@@ -0,0 +1,37 @@
+###########################################################################
+#
+# 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-PBDInjection)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME}
+ PBDInjectExample.cpp InflatableObject.h InflatableObject.cpp
+ imstkPbdInflatableVolumeConstraint.h imstkPbdInflatableVolumeConstraint.cpp
+ imstkPbdInflatableDistanceConstraint.h imstkPbdInflatableDistanceConstraint.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples/PBD)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} SimulationManager)
diff --git a/Examples/PBD/PBDInjection/InflatableObject.cpp b/Examples/PBD/PBDInjection/InflatableObject.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f0aa716cac778fd8757d7a85eb924fda4e8c491
--- /dev/null
+++ b/Examples/PBD/PBDInjection/InflatableObject.cpp
@@ -0,0 +1,395 @@
+/*=========================================================================
+
+ 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 "InflatableObject.h"
+#include "imstkNew.h"
+#include "imstkVecDataArray.h"
+#include "imstkTetrahedralMesh.h"
+#include "imstkSurfaceMesh.h"
+#include "imstkCollisionUtils.h"
+#include "imstkImageData.h"
+#include "imstkLogger.h"
+#include "imstkMeshIO.h"
+#include "imstkOneToOneMap.h"
+#include "imstkParallelFor.h"
+#include "imstkPbdConstraint.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkRenderMaterial.h"
+#include "imstkTexture.h"
+#include "imstkVisualModel.h"
+#include "imstkCollisionUtils.h"
+
+#include "imstkPbdInflatableDistanceConstraint.h"
+#include "imstkPbdInflatableVolumeConstraint.h"
+
+InflatableObject::InflatableObject(const std::string& name, const Vec3d& tissueSize, const Vec3i& tissueDim, const Vec3d& tissueCenter) : PbdObject(name)
+{
+ // Setup the Geometry
+ makeTetGrid(tissueSize, tissueDim, tissueCenter);
+ m_objectSurfMesh = m_objectTetMesh->extractSurfaceMesh();
+ m_objectSurfMesh->flipNormals();
+ //setXZPlaneTexCoords(4.0);
+ setSphereTexCoords(4.0);
+
+ // Setup the Parameters
+ imstkNew<PBDModelConfig> pbdParams;
+ pbdParams->m_doPartitioning = false;
+ pbdParams->m_uniformMassValue = 0.1;
+ pbdParams->m_gravity = Vec3d(0.0, 0.0, 0.0);
+ pbdParams->m_dt = 0.05;
+ pbdParams->m_iterations = 2;
+ pbdParams->m_viscousDampingCoeff = 0.05;
+
+ // Fix the borders
+ for (int z = 0; z < tissueDim[2]; z++)
+ {
+ for (int y = 0; y < tissueDim[1]; y++)
+ {
+ for (int x = 0; x < tissueDim[0]; x++)
+ {
+ if (x == 0 || z == 0 || x == tissueDim[0] - 1 || z == tissueDim[2] - 1 || y == 0)
+ {
+ pbdParams->m_fixedNodeIds.push_back(x + tissueDim[0] * (y + tissueDim[1] * z));
+ }
+ }
+ }
+ }
+
+ // Setup the Model
+ imstkNew<PbdModel> pbdModel;
+ pbdModel->setModelGeometry(m_objectTetMesh);
+ pbdModel->configure(pbdParams);
+
+ auto distanceFunctor = std::make_shared<PbdInflatableDistanceConstraintFunctor>();
+ distanceFunctor->setStiffness(0.95);
+ auto volumeFunctor = std::make_shared<PbdInflatableVolumeConstraintFunctor>();
+ volumeFunctor->setStiffness(0.9);
+ pbdModel->addPbdConstraintFunctor(distanceFunctor);
+ pbdModel->addPbdConstraintFunctor(volumeFunctor);
+
+ // Setup the material
+ imstkNew<RenderMaterial> material;
+ material->setBackFaceCulling(false);
+ material->setDisplayMode(RenderMaterial::DisplayMode::WireframeSurface);
+ material->setShadingModel(RenderMaterial::ShadingModel::PBR);
+ auto diffuseTex = MeshIO::read<ImageData>(iMSTK_DATA_ROOT "/textures/fleshDiffuse.jpg");
+ material->addTexture(std::make_shared<Texture>(diffuseTex, Texture::Type::Diffuse));
+ auto normalTex = MeshIO::read<ImageData>(iMSTK_DATA_ROOT "/textures/fleshNormal.jpg");
+ material->addTexture(std::make_shared<Texture>(normalTex, Texture::Type::Normal));
+ auto ormTex = MeshIO::read<ImageData>(iMSTK_DATA_ROOT "/textures/fleshORM.jpg");
+ material->addTexture(std::make_shared<Texture>(ormTex, Texture::Type::ORM));
+
+ // Add a visual model to render the surface of the tet mesh
+ imstkNew<VisualModel> visualModel(m_objectSurfMesh);
+ visualModel->setRenderMaterial(material);
+ addVisualModel(visualModel);
+
+ // Add a visual model to render the normals of the surface
+ imstkNew<VisualModel> normalsVisualModel(m_objectSurfMesh);
+ normalsVisualModel->getRenderMaterial()->setDisplayMode(RenderMaterial::DisplayMode::Surface);
+ normalsVisualModel->getRenderMaterial()->setPointSize(0.5);
+ addVisualModel(normalsVisualModel);
+
+ // Setup the Object
+ setPhysicsGeometry(m_objectTetMesh);
+ setCollidingGeometry(m_objectSurfMesh);
+ setPhysicsToCollidingMap(std::make_shared<OneToOneMap>(m_objectTetMesh, m_objectSurfMesh));
+ setDynamicalModel(pbdModel);
+}
+
+bool
+InflatableObject::initialize()
+{
+ PbdObject::initialize();
+
+ m_constraintContainer = m_pbdModel->getConstraints();
+
+ return true;
+}
+
+void
+InflatableObject::makeTetGrid(const Vec3d& size, const Vec3i& dim, const Vec3d& center)
+{
+ imstkNew<TetrahedralMesh> tissueMesh;
+
+ imstkNew<VecDataArray<double, 3>> verticesPtr(dim[0] * dim[1] * dim[2]);
+ VecDataArray<double, 3>& vertices = *verticesPtr.get();
+
+ const Vec3d dx = size.cwiseQuotient((dim - Vec3i(1, 1, 1)).cast<double>());
+ for (int z = 0; z < dim[2]; z++)
+ {
+ for (int y = 0; y < dim[1]; y++)
+ {
+ for (int x = 0; x < dim[0]; x++)
+ {
+ vertices[x + dim[0] * (y + dim[1] * z)] = Vec3i(x, y, z).cast<double>().cwiseProduct(dx) - size * 0.5 + center;
+ }
+ }
+ }
+
+ // Add connectivity data
+ imstkNew<VecDataArray<int, 4>> indicesPtr;
+ VecDataArray<int, 4>& indices = *indicesPtr.get();
+ for (int z = 0; z < dim[2] - 1; z++)
+ {
+ for (int y = 0; y < dim[1] - 1; y++)
+ {
+ for (int x = 0; x < dim[0] - 1; x++)
+ {
+ int cubeIndices[8] =
+ {
+ x + dim[0] * (y + dim[1] * z),
+ (x + 1) + dim[0] * (y + dim[1] * z),
+ (x + 1) + dim[0] * (y + dim[1] * (z + 1)),
+ x + dim[0] * (y + dim[1] * (z + 1)),
+ x + dim[0] * ((y + 1) + dim[1] * z),
+ (x + 1) + dim[0] * ((y + 1) + dim[1] * z),
+ (x + 1) + dim[0] * ((y + 1) + dim[1] * (z + 1)),
+ x + dim[0] * ((y + 1) + dim[1] * (z + 1))
+ };
+
+ // Alternate the pattern so the edges line up on the sides of each voxel
+ if ((z % 2 ^ x % 2) ^ y % 2)
+ {
+ indices.push_back(Vec4i(cubeIndices[0], cubeIndices[7], cubeIndices[5], cubeIndices[4]));
+ indices.push_back(Vec4i(cubeIndices[3], cubeIndices[7], cubeIndices[2], cubeIndices[0]));
+ indices.push_back(Vec4i(cubeIndices[2], cubeIndices[7], cubeIndices[5], cubeIndices[0]));
+ indices.push_back(Vec4i(cubeIndices[1], cubeIndices[2], cubeIndices[0], cubeIndices[5]));
+ indices.push_back(Vec4i(cubeIndices[2], cubeIndices[6], cubeIndices[7], cubeIndices[5]));
+ }
+ else
+ {
+ indices.push_back(Vec4i(cubeIndices[3], cubeIndices[7], cubeIndices[6], cubeIndices[4]));
+ indices.push_back(Vec4i(cubeIndices[1], cubeIndices[3], cubeIndices[6], cubeIndices[4]));
+ indices.push_back(Vec4i(cubeIndices[3], cubeIndices[6], cubeIndices[2], cubeIndices[1]));
+ indices.push_back(Vec4i(cubeIndices[1], cubeIndices[6], cubeIndices[5], cubeIndices[4]));
+ indices.push_back(Vec4i(cubeIndices[0], cubeIndices[3], cubeIndices[1], cubeIndices[4]));
+ }
+ }
+ }
+ }
+
+ imstkNew<VecDataArray<float, 2>> uvCoordsPtr(dim[0] * dim[1] * dim[2]);
+ VecDataArray<float, 2>& uvCoords = *uvCoordsPtr.get();
+ for (int z = 0; z < dim[2]; z++)
+ {
+ for (int y = 0; y < dim[1]; y++)
+ {
+ for (int x = 0; x < dim[0]; x++)
+ {
+ uvCoords[x + dim[0] * (y + dim[1] * z)] = Vec2f(static_cast<float>(x) / dim[0], static_cast<float>(z) / dim[2]) * 3.0;
+ }
+ }
+ }
+
+ // Ensure correct windings
+ for (int i = 0; i < indices.size(); i++)
+ {
+ if (tetVolume(vertices[indices[i][0]], vertices[indices[i][1]], vertices[indices[i][2]], vertices[indices[i][3]]) < 0.0)
+ {
+ std::swap(indices[i][0], indices[i][2]);
+ }
+ }
+
+ tissueMesh->initialize(verticesPtr, indicesPtr);
+ tissueMesh->setVertexTCoords("uvs", uvCoordsPtr);
+
+ m_objectTetMesh = tissueMesh;
+}
+
+void
+InflatableObject::setXZPlaneTexCoords(const double uvScale)
+{
+ Vec3d min, max;
+ m_objectSurfMesh->computeBoundingBox(min, max);
+ const Vec3d size = max - min;
+
+ imstkNew<VecDataArray<float, 2>> uvCoordsPtr(m_objectSurfMesh->getNumVertices());
+ VecDataArray<float, 2>& uvCoords = *uvCoordsPtr.get();
+ for (int i = 0; i < m_objectSurfMesh->getNumVertices(); i++)
+ {
+ const Vec3d& vertex = m_objectSurfMesh->getVertexPosition(i);
+ uvCoords[i] = Vec2f((vertex[0] - min[0]) / size[0], (vertex[2] - min[2]) / size[2]) * uvScale;
+ }
+ m_objectSurfMesh->setVertexTCoords("tcoords", uvCoordsPtr);
+}
+
+void
+InflatableObject::setSphereTexCoords(const double uvScale)
+{
+ Vec3d min, max;
+ m_objectSurfMesh->computeBoundingBox(min, max);
+ const Vec3d size = max - min;
+ const Vec3d center = (max + min) * 0.5;
+
+ const double radius = (size * 0.5).norm();
+
+ imstkNew<VecDataArray<float, 2>> uvCoordsPtr(m_objectSurfMesh->getNumVertices());
+ VecDataArray<float, 2>& uvCoords = *uvCoordsPtr.get();
+ for (int i = 0; i < m_objectSurfMesh->getNumVertices(); i++)
+ {
+ Vec3d vertex = m_objectSurfMesh->getVertexPosition(i) - center;
+
+ // Compute phi and theta on the sphere
+ const double theta = asin(vertex[0] / radius);
+ const double phi = atan2(vertex[1], vertex[2]);
+ uvCoords[i] = Vec2f(phi / (PI * 2.0) + 0.5, theta / (PI * 2.0) + 0.5) * uvScale;
+ }
+ m_objectSurfMesh->setVertexTCoords("tcoords", uvCoordsPtr);
+}
+
+void
+InflatableObject::findAffectedConstraint(const Vec3d& toolTip, const double radius)
+{
+ m_constraintIDandWeight.clear();
+
+ Vec3d min, max;
+ m_objectTetMesh->computeBoundingBox(min, max);
+ if (!CollisionUtils::testAABBToAABB(toolTip[0], toolTip[0], toolTip[1], toolTip[1], toolTip[2], toolTip[2],
+ min[0], max[0], min[1], max[1], min[2], max[2]))
+ {
+ return;
+ }
+
+ const auto& vertices = m_objectTetMesh->getVertexPositions();
+
+ std::vector<int> constraintIDs;
+ std::vector<double> weights;
+ int counter = 0;
+ double minDistance = LONG_MAX;
+ for (auto& c : m_constraintContainer->getConstraints())
+ {
+ auto& ids = c->getVertexIds();
+
+ Vec3d center(0.0, 0.0, 0.0);
+ for (auto i : ids)
+ {
+ center += (*vertices)[i];
+ }
+
+ double distance = (center / ids.size() - toolTip).norm();
+
+ if (distance < radius)
+ {
+ constraintIDs.push_back(counter);
+ weights.push_back(computeGaussianWeight(distance));
+ }
+
+ if (distance < minDistance)
+ {
+ minDistance = distance;
+ }
+
+ counter++;
+ }
+
+ if (minDistance > 0.5)
+ {
+ return;
+ }
+ else
+ {
+ for (int i = 0; i < constraintIDs.size(); i++)
+ {
+ m_constraintIDandWeight.push_back(std::make_pair(constraintIDs[i], weights[i]));
+ }
+
+ m_affectedAreaUpdated = true;
+ }
+}
+
+void
+InflatableObject::inject(const Vec3d& toolTip, const double radius, double rate)
+{
+ if (!m_affectedAreaUpdated)
+ {
+ findAffectedConstraint(toolTip, radius);
+ }
+
+ if (m_inlationTpye == InflationType::Exponential)
+ {
+ m_inflationRatio *= (1 + rate);
+ }
+ else if (m_inlationTpye == InflationType::Linear)
+ {
+ m_inflationRatio += rate;
+ }
+
+ for (auto& id : m_constraintIDandWeight)
+ {
+ auto& c = (m_constraintContainer->getConstraints())[id.first];
+ if (c->getType() == PbdConstraint::Type::Volume)
+ {
+ const auto& constraint = std::dynamic_pointer_cast<PbdInflatableVolumeConstraint>(c);
+ constraint->multiplyInitRestVolumeBy(m_inflationRatio * id.second);
+ }
+ else if (c->getType() == PbdConstraint::Type::Distance)
+ {
+ //const auto& constraint = std::dynamic_pointer_cast<PbdInflatableDistanceConstraint>(c);
+ //constraint->multiplyInitRestLengthBy(0.1 * m_inflationRatio * id.second);
+ }
+ }
+}
+
+void
+InflatableObject::switchInflationType()
+{
+ if (m_inlationTpye == InflationType::Linear)
+ {
+ m_inlationTpye = InflationType::Exponential;
+ std::cout << "Inflation Type: Exponential." << std::endl;
+ }
+ else if (m_inlationTpye == InflationType::Exponential)
+ {
+ m_inlationTpye = InflationType::Linear;
+ std::cout << "Inflation Type: Linear." << std::endl;
+ }
+}
+
+void
+InflatableObject::setUpdateAffectedConstraint()
+{
+ m_affectedAreaUpdated = false;
+ m_inflationRatio = 1.0;
+}
+
+void
+InflatableObject::reset()
+{
+ PbdObject::reset();
+
+ m_inflationRatio = 1.0;
+
+ for (auto& c : m_constraintContainer->getConstraints())
+ {
+ if (c->getType() == PbdConstraint::Type::Volume)
+ {
+ const auto& constraint = std::dynamic_pointer_cast<PbdInflatableVolumeConstraint>(c);
+ constraint->resetRestVolume();
+ }
+ else if (c->getType() == PbdConstraint::Type::Distance)
+ {
+ const auto& constraint = std::dynamic_pointer_cast<PbdInflatableDistanceConstraint>(c);
+ constraint->resetRestLength();
+ }
+ }
+}
diff --git a/Examples/PBD/PBDInjection/InflatableObject.h b/Examples/PBD/PBDInjection/InflatableObject.h
new file mode 100644
index 0000000000000000000000000000000000000000..a9c26fdf1bd1c91f4ac7ecf8ac9415250c08a1ab
--- /dev/null
+++ b/Examples/PBD/PBDInjection/InflatableObject.h
@@ -0,0 +1,118 @@
+/*=========================================================================
+
+ 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"
+#include "imstkPbdConstraintContainer.h"
+
+namespace imstk
+{
+class ImageData;
+class SurfaceMesh;
+class TetrahedralMesh;
+class Texture;
+};
+
+using namespace imstk;
+
+///
+/// \class InflatableObject
+///
+/// \brief Inflatable object based on PBD, with inflatable volume and distance constraints
+///
+class InflatableObject : public PbdObject
+{
+public:
+ enum class InflationType
+ {
+ Linear,
+ Exponential,
+ None
+ };
+
+ InflatableObject(const std::string& name, const Vec3d& tissueSize, const Vec3i& tissueDim, const Vec3d& tissueCenter);
+
+public:
+ const std::string getTypeName() const override { return "InflatableObject"; }
+
+ bool initialize() override;
+
+ ///
+ /// \brief Perform injection on the tissue given tool tip position
+ ///
+ void inject(const Vec3d& toolTip, const double radius, double rate);
+
+ ///
+ /// \brief Switch between linear and exponential inflation type
+ ///
+ void switchInflationType();
+
+ void setInflationRatio(double ratio) { m_inflationRatio = ratio; }
+ void setInflationSize(double sigma) { m_sigma = sigma; }
+
+ ///
+ /// \brief set update affected constraints flag
+ ///
+ void setUpdateAffectedConstraint();
+
+ void reset() override;
+
+protected:
+ ///
+ /// \brief Creates a tetraheral grid
+ ///
+ void makeTetGrid(const Vec3d& size, const Vec3i& dim, const Vec3d& center);
+
+ ///
+ /// \brief Compute the texture coordinates of the tissue on a plane
+ ///
+ void setXZPlaneTexCoords(const double uvScale);
+
+ ///
+ /// \brief Spherically project the texture coordinates
+ ///
+ void setSphereTexCoords(const double uvScale);
+
+ ///
+ /// \brief find affected constraints id and distance in the injection area
+ ///
+ void findAffectedConstraint(const Vec3d& toolTip, const double radius);
+
+ ///
+ /// \brief Compute weight Gaussian distribution
+ ///
+ inline double computeGaussianWeight(double x) { return 10 * exp(-0.5 * x * x / m_sigma / m_sigma) / m_sigma; }
+
+protected:
+ std::shared_ptr<TetrahedralMesh> m_objectTetMesh;
+ std::shared_ptr<SurfaceMesh> m_objectSurfMesh;
+
+ std::shared_ptr<PbdConstraintContainer> m_constraintContainer;
+ std::vector<std::pair<int, double>> m_constraintIDandWeight;
+
+ bool m_affectedAreaUpdated = false;
+
+ double m_inflationRatio = 1.0;
+ double m_sigma = 1.0;
+
+ InflationType m_inlationTpye = InflationType::Linear;
+};
\ No newline at end of file
diff --git a/Examples/PBD/PBDInjection/PBDInjectExample.cpp b/Examples/PBD/PBDInjection/PBDInjectExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..37de20a2fe9e76c7c02a4310be583ac9047cac69
--- /dev/null
+++ b/Examples/PBD/PBDInjection/PBDInjectExample.cpp
@@ -0,0 +1,249 @@
+/*=========================================================================
+
+ 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 "imstkCamera.h"
+#include "imstkCollisionGraph.h"
+#include "imstkImageData.h"
+#include "imstkKeyboardDeviceClient.h"
+#include "imstkKeyboardSceneControl.h"
+#include "imstkLight.h"
+#include "imstkMouseSceneControl.h"
+#include "imstkOneToOneMap.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkRbdConstraint.h"
+#include "imstkRigidBodyModel2.h"
+#include "imstkRigidObject2.h"
+#include "imstkScene.h"
+#include "imstkSceneManager.h"
+#include "imstkSimulationManager.h"
+#include "imstkTetrahedralMesh.h"
+#include "imstkVTKViewer.h"
+#include "imstkLineMesh.h"
+#include "imstkVisualModel.h"
+#include "imstkRenderMaterial.h"
+#include "imstkLogger.h"
+
+#include "imstkMeshToMeshBruteForceCD.h"
+#include "imstkNew.h"
+#include "InflatableObject.h"
+
+#ifdef EXAMPLE_USE_HAPTICS
+#include "imstkHapticDeviceManager.h"
+#include "imstkHapticDeviceClient.h"
+#include "imstkRigidObjectController.h"
+#endif
+
+using namespace imstk;
+
+static std::shared_ptr<RigidObject2>
+makeToolObj(const std::string& name)
+{
+ // Setup a Sphere tool
+ //imstkNew<Sphere> toolGeometry(Vec3d(0.0, 0.0, 0.0), 1.0);
+
+ imstkNew<LineMesh> toolGeometry;
+ imstkNew<VecDataArray<double, 3>> verticesPtr(2);
+ (*verticesPtr)[0] = Vec3d(0.0, 0.0, 0.0);
+ (*verticesPtr)[1] = Vec3d(0.0, 2.0, 0.0);
+ imstkNew<VecDataArray<int, 2>> indicesPtr(1);
+ (*indicesPtr)[0] = Vec2i(0, 1);
+ toolGeometry->initialize(verticesPtr, indicesPtr);
+
+ imstkNew<RigidObject2> toolObj(name);
+ toolObj->setVisualGeometry(toolGeometry);
+ toolObj->setCollidingGeometry(toolGeometry);
+ toolObj->setPhysicsGeometry(toolGeometry);
+ toolObj->getVisualModel(0)->getRenderMaterial()->setColor(Color::Blue);
+ toolObj->getVisualModel(0)->getRenderMaterial()->setDisplayMode(RenderMaterial::DisplayMode::Wireframe);
+ toolObj->getVisualModel(0)->getRenderMaterial()->setBackFaceCulling(false);
+ toolObj->getVisualModel(0)->getRenderMaterial()->setLineWidth(10.0);
+
+ std::shared_ptr<RigidBodyModel2> rbdModel = std::make_shared<RigidBodyModel2>();
+ rbdModel->getConfig()->m_gravity = Vec3d::Zero();
+ rbdModel->getConfig()->m_maxNumIterations = 2;
+ toolObj->setDynamicalModel(rbdModel);
+
+ toolObj->getRigidBody()->m_mass = 1000.0;
+ toolObj->getRigidBody()->m_intertiaTensor = Mat3d::Identity() * 10000.0;
+ toolObj->getRigidBody()->m_initPos = Vec3d(0.0, 0.8, 0.0);
+ toolObj->getRigidBody()->m_isStatic = true;
+
+ return toolObj;
+}
+
+///
+/// \brief This example demonstrates the inflatable distance and volume Constraint
+/// using Position based dynamics
+///
+int
+main()
+{
+ // Setup logger (write to file and stdout)
+ Logger::startLogger();
+
+ // Setup the scene
+ imstkNew<Scene> scene("PBD inflatable object example");
+ scene->getActiveCamera()->setPosition(0.12, 4.51, 16.51);
+ scene->getActiveCamera()->setFocalPoint(0.0, 0.0, 0.0);
+ scene->getActiveCamera()->setViewUp(0.0, 0.96, -0.28);
+
+ // Setup a tissue
+ Vec3d tissueSize = Vec3d(10.0, 3.0, 10.0);
+ Vec3i tissueDim = Vec3i(20, 5, 20);
+ Vec3d tissueCenter = Vec3d(0.1, -1.0, 0.0);
+ double radius = tissueSize[0] / 5.0;
+ imstkNew<InflatableObject> tissueObj("PBDTissue", tissueSize, tissueDim, tissueCenter);
+ scene->addSceneObject(tissueObj);
+
+ // Setup a tool
+ Vec3d toolTip = tissueCenter + Vec3d(0.0, tissueSize[1] / 2.0, 0.0);
+ std::shared_ptr<RigidObject2> toolObj = makeToolObj("RBDTool");
+ scene->addSceneObject(toolObj);
+
+ //auto interaction = std::make_shared<PbdRigidObjectCollision>(tissueObj, toolObj, "MeshToMeshBruteForceCD");
+ //scene->getCollisionGraph()->addInteraction(interaction);
+
+ // Light
+ imstkNew<DirectionalLight> light("Light");
+ light->setFocalPoint(Vec3d(5, -8, -5));
+ light->setIntensity(1);
+ scene->addLight(light);
+
+ // Run the simulation
+ {
+ // Setup a viewer to render
+ imstkNew<VTKViewer> viewer("Viewer");
+ viewer->setActiveScene(scene);
+ viewer->setVtkLoggerMode(VTKViewer::VTKLoggerMode::MUTE);
+ viewer->setDebugAxesLength(0.1, 0.1, 0.1);
+
+ // Setup a scene manager to advance the scene
+ imstkNew<SceneManager> sceneManager("Scene Manager");
+ sceneManager->setActiveScene(scene);
+ sceneManager->setExecutionType(Module::ExecutionType::ADAPTIVE);
+ sceneManager->pause(); // Start simulation paused
+
+ imstkNew<SimulationManager> driver;
+ driver->addModule(viewer);
+ driver->addModule(sceneManager);
+ driver->setDesiredDt(0.01);
+
+ // Add mouse and keyboard controls to the viewer
+ imstkNew<MouseSceneControl> mouseControl(viewer->getMouseDevice());
+ mouseControl->setSceneManager(sceneManager);
+ viewer->addControl(mouseControl);
+ imstkNew<KeyboardSceneControl> keyControl(viewer->getKeyboardDevice());
+ keyControl->setSceneManager(sceneManager);
+ keyControl->setModuleDriver(driver);
+ viewer->addControl(keyControl);
+
+#ifdef EXAMPLE_USE_HAPTICS
+ imstkNew<HapticDeviceManager> hapticManager;
+ hapticManager->setSleepDelay(1.0); // Delay for 1ms (haptics thread is limited to max 1000hz)
+ std::shared_ptr<HapticDeviceClient> hapticDeviceClient = hapticManager->makeDeviceClient();
+ driver->addModule(hapticManager);
+
+ imstkNew<RigidObjectController> controller(toolObj, hapticDeviceClient);
+ controller->setTranslationScaling(0.05);
+ controller->setLinearKs(1000.0);
+ controller->setLinearKd(50.0);
+ controller->setAngularKs(10000000.0);
+ controller->setAngularKd(500000.0);
+ controller->setForceScaling(0.005);
+ controller->setSmoothingKernelSize(25);
+ controller->setUseForceSmoothening(true);
+ scene->addController(controller);
+#else
+ // Use keyboard controls
+ connect<Event>(sceneManager, SceneManager::preUpdate, [&](Event*)
+ {
+ if (viewer->getKeyboardDevice()->getButton('k') == KEY_PRESS)
+ {
+ (*toolObj->getRigidBody()->m_pos) += Vec3d(0.0, -0.01, 0.0);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('i') == KEY_PRESS)
+ {
+ (*toolObj->getRigidBody()->m_pos) += Vec3d(0.0, 0.01, 0.0);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('j') == KEY_PRESS)
+ {
+ (*toolObj->getRigidBody()->m_pos) += Vec3d(-0.01, 0.0, 0.0);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('l') == KEY_PRESS)
+ {
+ (*toolObj->getRigidBody()->m_pos) += Vec3d(0.01, 0.0, 0.0);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('u') == KEY_PRESS)
+ {
+ (*toolObj->getRigidBody()->m_pos) += Vec3d(0.0, 0.0, -0.01);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('o') == KEY_PRESS)
+ {
+ (*toolObj->getRigidBody()->m_pos) += Vec3d(0.0, 0.0, 0.01);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('s') == KEY_PRESS)
+ {
+ // The LineMesh used for collision with the PBD tissue
+ std::shared_ptr<LineMesh> lineMesh = std::dynamic_pointer_cast<LineMesh>(toolObj->getCollidingGeometry());
+ Vec3d vertice = lineMesh->getVertexPosition(0);
+
+ if ((toolTip - vertice).norm() > 0.01)
+ {
+ toolTip = vertice;
+ tissueObj->setUpdateAffectedConstraint();
+ }
+
+ tissueObj->inject(toolTip, radius, 0.01);
+ }
+ else if (viewer->getKeyboardDevice()->getButton('a') == KEY_PRESS)
+ {
+ // The LineMesh used for collision with the PBD tissue
+ std::shared_ptr<LineMesh> lineMesh = std::dynamic_pointer_cast<LineMesh>(toolObj->getCollidingGeometry());
+ Vec3d vertice = lineMesh->getVertexPosition(0);
+
+ if ((toolTip - vertice).norm() > 0.01)
+ {
+ toolTip = vertice;
+ tissueObj->setUpdateAffectedConstraint();
+ }
+
+ tissueObj->inject(toolTip, radius, -0.01);
+ }
+ });
+#endif
+
+ connect<Event>(sceneManager, &SceneManager::postUpdate, [&](Event*)
+ {
+ // Keep the tool moving in real time
+ toolObj->getRigidBodyModel2()->getConfig()->m_dt = sceneManager->getDt();
+ //tissueObj->getPbdModel()->getParameters()->m_dt = sceneManager->getDt();
+ });
+
+ std::cout << "================================================" << std::endl;
+ std::cout << "Key s : injection \n" << "Key a : deflation \n";
+ std::cout << "================================================" << std::endl << std::endl;
+
+ driver->start();
+ }
+
+ return 0;
+}
\ No newline at end of file
diff --git a/Examples/PBD/PBDInjection/imstkPbdInflatableDistanceConstraint.cpp b/Examples/PBD/PBDInjection/imstkPbdInflatableDistanceConstraint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..db72f94c2e24e4999505e3ec808152d00683d011
--- /dev/null
+++ b/Examples/PBD/PBDInjection/imstkPbdInflatableDistanceConstraint.cpp
@@ -0,0 +1,53 @@
+/*=========================================================================
+
+ 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 "imstkPbdInflatableDistanceConstraint.h"
+
+namespace imstk
+{
+void
+PbdInflatableDistanceConstraint::initConstraint(const VecDataArray<double, 3>& initVertexPositions,
+ const size_t& pIdx0,
+ const size_t& pIdx1,
+ const double k)
+{
+ PbdDistanceConstraint::initConstraint(initVertexPositions, pIdx0, pIdx1, k);
+ m_initialRestLength = m_restLength;
+}
+
+void
+PbdInflatableDistanceConstraint::multiplyRestLengthBy(const double ratio)
+{
+ // Multiplied ratio must be positive
+ if (ratio < DBL_EPSILON)
+ {
+ return;
+ }
+
+ m_restLength *= ratio;
+
+ // Modified RestLenght can't be less than initial RestLenght
+ if (m_restLength < m_initialRestLength)
+ {
+ m_restLength = m_initialRestLength;
+ }
+}
+}
diff --git a/Examples/PBD/PBDInjection/imstkPbdInflatableDistanceConstraint.h b/Examples/PBD/PBDInjection/imstkPbdInflatableDistanceConstraint.h
new file mode 100644
index 0000000000000000000000000000000000000000..9efad07af4875db891d7494f5406af82d8e0c877
--- /dev/null
+++ b/Examples/PBD/PBDInjection/imstkPbdInflatableDistanceConstraint.h
@@ -0,0 +1,139 @@
+/*=========================================================================
+
+ 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 "imstkPbdDistanceConstraint.h"
+#include "imstkPbdConstraintFunctor.h"
+
+namespace imstk
+{
+///
+/// \class PbdInflatableDistanceConstraint
+///
+class PbdInflatableDistanceConstraint : public PbdDistanceConstraint
+{
+public:
+ ///
+ /// \brief Constructor
+ ///
+ PbdInflatableDistanceConstraint() : PbdDistanceConstraint() {}
+
+ ///
+ /// \brief Initializes the inflatable distance constraint
+ ///
+ void initConstraint(const VecDataArray<double, 3>& initVertexPositions,
+ const size_t& pIdx0,
+ const size_t& pIdx1,
+ const double k = 1e5);
+
+ ///
+ /// \brief Soften or strengthen constraint by multiplying current restLength by a ratio
+ ///
+ void multiplyRestLengthBy(const double ratio);
+
+ ///
+ /// \brief Soften or strengthen constraint by multiplying initial restLength by a ratio
+ ///
+ void multiplyInitRestLengthBy(const double ratio) { m_restLength = ratio * m_initialRestLength; }
+
+ ///
+ /// \brief Reset constraint rest measurement
+ ///
+ void resetRestLength() { m_restLength = m_initialRestLength; }
+
+public:
+ double m_initialRestLength = 0.0; ///> Rest measurement(length, area, volume, etc.)
+};
+
+struct PbdInflatableDistanceConstraintFunctor : public PbdDistanceConstraintFunctor
+{
+ PbdInflatableDistanceConstraintFunctor() : PbdDistanceConstraintFunctor() {}
+ ~PbdInflatableDistanceConstraintFunctor() override = default;
+
+ virtual void operator()(PbdConstraintContainer& constraints) override
+ {
+ const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
+ auto addDistConstraint =
+ [&](std::vector<std::vector<bool>>& E, size_t i1, size_t i2)
+ {
+ if (i1 > i2) // Make sure i1 is always smaller than i2
+ {
+ std::swap(i1, i2);
+ }
+ if (E[i1][i2])
+ {
+ E[i1][i2] = 0;
+ auto c = std::make_shared<PbdInflatableDistanceConstraint>();
+ c->initConstraint(vertices, i1, i2, m_stiffness);
+ constraints.addConstraint(c);
+ }
+ };
+
+ if (m_geom->getTypeName() == "TetrahedralMesh")
+ {
+ const auto& tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
+ const VecDataArray<int, 4>& elements = *tetMesh->getTetrahedraIndices();
+ const auto nV = tetMesh->getNumVertices();
+ std::vector<std::vector<bool>> E(nV, std::vector<bool>(nV, 1));
+
+ for (int k = 0; k < elements.size(); ++k)
+ {
+ auto& tet = elements[k];
+ addDistConstraint(E, tet[0], tet[1]);
+ addDistConstraint(E, tet[0], tet[2]);
+ addDistConstraint(E, tet[0], tet[3]);
+ addDistConstraint(E, tet[1], tet[2]);
+ addDistConstraint(E, tet[1], tet[3]);
+ addDistConstraint(E, tet[2], tet[3]);
+ }
+ }
+ else if (m_geom->getTypeName() == "SurfaceMesh")
+ {
+ const auto& triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
+ const VecDataArray<int, 3>& elements = *triMesh->getTriangleIndices();
+ const auto nV = triMesh->getNumVertices();
+ std::vector<std::vector<bool>> E(nV, std::vector<bool>(nV, 1));
+
+ for (int k = 0; k < elements.size(); ++k)
+ {
+ auto& tri = elements[k];
+ addDistConstraint(E, tri[0], tri[1]);
+ addDistConstraint(E, tri[0], tri[2]);
+ addDistConstraint(E, tri[1], tri[2]);
+ }
+ }
+ else if (m_geom->getTypeName() == "LineMesh")
+ {
+ const auto& lineMesh = std::dynamic_pointer_cast<LineMesh>(m_geom);
+ const VecDataArray<int, 2>& elements = *lineMesh->getLinesIndices();
+ const auto& nV = lineMesh->getNumVertices();
+ std::vector<std::vector<bool>> E(nV, std::vector<bool>(nV, 1));
+
+ for (int k = 0; k < elements.size(); k++)
+ {
+ auto& seg = elements[k];
+ addDistConstraint(E, seg[0], seg[1]);
+ }
+ }
+ }
+};
+} // imstk
diff --git a/Examples/PBD/PBDInjection/imstkPbdInflatableVolumeConstraint.cpp b/Examples/PBD/PBDInjection/imstkPbdInflatableVolumeConstraint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bb4a6f4453eec774423cabbe502f4eb143801745
--- /dev/null
+++ b/Examples/PBD/PBDInjection/imstkPbdInflatableVolumeConstraint.cpp
@@ -0,0 +1,53 @@
+/*=========================================================================
+
+ 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 "imstkPbdInflatableVolumeConstraint.h"
+
+namespace imstk
+{
+void
+PbdInflatableVolumeConstraint::initConstraint(const VecDataArray<double, 3>& initVertexPositions,
+ const size_t& pIdx0, const size_t& pIdx1,
+ const size_t& pIdx2, const size_t& pIdx3,
+ const double k)
+{
+ PbdVolumeConstraint::initConstraint(initVertexPositions, pIdx0, pIdx1, pIdx2, pIdx3, k);
+ m_initialRestVolume = m_restVolume;
+}
+
+void
+PbdInflatableVolumeConstraint::multiplyRestVolumeBy(const double ratio)
+{
+ // Multiplied ratio must be positive
+ if (ratio < DBL_EPSILON)
+ {
+ return;
+ }
+
+ m_restVolume *= ratio;
+
+ // Modified RestVolume can't be less than initial RestVolume
+ if (abs(m_restVolume) < abs(m_initialRestVolume))
+ {
+ m_restVolume = m_initialRestVolume;
+ }
+}
+} // imstk
diff --git a/Examples/PBD/PBDInjection/imstkPbdInflatableVolumeConstraint.h b/Examples/PBD/PBDInjection/imstkPbdInflatableVolumeConstraint.h
new file mode 100644
index 0000000000000000000000000000000000000000..c5db102a893a6f06e6eaef863ea33b2d00f866a7
--- /dev/null
+++ b/Examples/PBD/PBDInjection/imstkPbdInflatableVolumeConstraint.h
@@ -0,0 +1,94 @@
+/*=========================================================================
+
+ 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 "imstkPbdVolumeConstraint.h"
+#include "imstkPbdConstraintFunctor.h"
+
+namespace imstk
+{
+///
+/// \class PbdInflatableVolumeConstraint
+///
+class PbdInflatableVolumeConstraint : public PbdVolumeConstraint
+{
+public:
+ ///
+ /// \brief constructor
+ ///
+ PbdInflatableVolumeConstraint() : PbdVolumeConstraint() {}
+
+ ///
+ /// \brief Initializes the inflatable volume constraint
+ ///
+ void initConstraint(const VecDataArray<double, 3>& initVertexPositions,
+ const size_t& pIdx1, const size_t& pIdx2,
+ const size_t& pIdx3, const size_t& pIdx4,
+ const double k = 2.0);
+
+ ///
+ /// \brief Soften or strengthen constraint by multiplying current restVolume by a ratio
+ ///
+ void multiplyRestVolumeBy(const double ratio);
+
+ ///
+ /// \brief Soften or strengthen constraint by multiplying initial restVolume by a ratio
+ ///
+ void multiplyInitRestVolumeBy(const double ratio) { m_restVolume = ratio * m_initialRestVolume; }
+
+ ///
+ /// \brief Reset constraint rest volume
+ ///
+ void resetRestVolume() { m_restVolume = m_initialRestVolume; }
+
+protected:
+ double m_initialRestVolume = 0.0; ///> Rest measurement(length, area, volume, etc.)
+};
+
+struct PbdInflatableVolumeConstraintFunctor : public PbdVolumeConstraintFunctor
+{
+ PbdInflatableVolumeConstraintFunctor() : PbdVolumeConstraintFunctor() {}
+ ~PbdInflatableVolumeConstraintFunctor() override = default;
+
+ virtual void operator()(PbdConstraintContainer& constraints) override
+ {
+ // Check if constraint type matches the mesh type
+ CHECK(m_geom->getTypeName() == "TetrahedralMesh")
+ << "Volume constraint should come with volumetric mesh";
+
+ // Create constraints
+ auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
+ const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
+ const VecDataArray<int, 4>& elements = *tetMesh->getTetrahedraIndices();
+
+ ParallelUtils::parallelFor(elements.size(),
+ [&](const size_t k)
+ {
+ auto& tet = elements[k];
+ auto c = std::make_shared<PbdInflatableVolumeConstraint>();
+ c->initConstraint(vertices,
+ tet[0], tet[1], tet[2], tet[3], m_stiffness);
+ constraints.addConstraint(c);
+ });
+ }
+};
+} // imstk