diff --git a/CMake/External/CMakeLists.txt b/CMake/External/CMakeLists.txt index a60b70369204dc2662c31a518abea367dd7e3cfd..79de0e0df0f1d897ef9fc46be9c4b5b405fcfed0 100644 --- a/CMake/External/CMakeLists.txt +++ b/CMake/External/CMakeLists.txt @@ -12,7 +12,7 @@ endif() #----------------------------------------------------------------------------- # Output Directories #----------------------------------------------------------------------------- -if(NOT DEFINED CMAKE_LIBRARY_OUTPUT_DIRECTORY) +if(NOT DEFINED CMAKE_LIBRARY_OUTPUT_DIRECTORY) set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/${${PROJECT_NAME}_LIB_DIR}) endif() if(NOT DEFINED CMAKE_ARCHIVE_OUTPUT_DIRECTORY) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5b37a885c4b43ee50a0929100a2b2086ad061694..c42f5729c80b9466d10514ea026bd2e7d92d6213 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,7 +18,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/CMake/Utilities ${CMAKE_MODULE_PATH} ) -set(${PROJECT_NAME}_CMAKE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/CMake) +set(${PROJECT_NAME}_CMAKE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/CMake) #----------------------------------------------------------------------------- # Set a default build type if none was specified diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index 274d53b038823f1011fdf7148391817637a79409..30968b0f611e3aa4138e1be5c7a97e977435b4a0 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -16,6 +16,6 @@ macro(listOfSubDir result curdir) foreach(subdir ${subDirs}) add_subdirectory(${subdir}) -endforeach() +endforeach() diff --git a/Examples/DeformableBody/CMakeLists.txt b/Examples/DeformableBody/CMakeLists.txt index 3cca576ce0959b0ecf32f941d44584e0e117d9e0..725688edcd98318500b05718b577920bca4a4eb7 100644 --- a/Examples/DeformableBody/CMakeLists.txt +++ b/Examples/DeformableBody/CMakeLists.txt @@ -27,7 +27,7 @@ add_executable(${PROJECT_NAME} DeformableBody.cpp) # Link libraries to executable #----------------------------------------------------------------------------- target_link_libraries(${PROJECT_NAME} SimulationManager) - + #----------------------------------------------------------------------------- # Add shaders #----------------------------------------------------------------------------- diff --git a/Examples/Picking/Picking.cpp b/Examples/Picking/Picking.cpp index 0be9f4255db9fec75717e107e419701310bb0682..73805a90742f960175a9f83861bd52ba3cc09b5f 100644 --- a/Examples/Picking/Picking.cpp +++ b/Examples/Picking/Picking.cpp @@ -69,6 +69,7 @@ int main() LOG(WARNING) << "Could not read mesh from file."; return 1; } + tetMesh->scale(10., Geometry::TransformType::ApplyToData); // Extract the surface mesh auto volTetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(tetMesh); if (!volTetMesh) @@ -109,10 +110,9 @@ int main() //---------------------------------------------------------- auto nlSystem = std::make_shared<NonLinearSystem>(dynaModel->getFunction(), dynaModel->getFunctionGradient()); std::vector<LinearProjectionConstraint> linProj; - for (auto id : dynaModel->getFixNodeIds()) - { - linProj.push_back(LinearProjectionConstraint(id, true)); - } + linProj.push_back(LinearProjectionConstraint(0, true)); + linProj.push_back(LinearProjectionConstraint(2, true)); + nlSystem->setUnknownVector(dynaModel->getUnknownVec()); nlSystem->setUpdateFunction(dynaModel->getUpdateFunction()); nlSystem->setUpdatePreviousStatesFunction(dynaModel->getUpdatePrevStateFunction()); @@ -141,7 +141,7 @@ int main() // Sphere0 auto sphereForPickObj = apiutils::createCollidingAnalyticalSceneObject( - Geometry::Type::Sphere, scene, "Sphere0", 1, Vec3d(0., 0., 0.)); + Geometry::Type::Sphere, scene, "Sphere0", 0.5, Vec3d(0., 0., 0.)); auto pickTrackingCtrl = std::make_shared<DeviceTracker>(client); //pickTrackingCtrl->setTranslationOffset(Vec3d(0., 0., 24.)); diff --git a/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.cpp b/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.cpp index 09b6030da2c3b9fd775a61b770b10815e0e7dad3..5986e3f74a690762a2bf884b4b65493916d4ee6b 100644 --- a/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.cpp +++ b/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.cpp @@ -41,6 +41,14 @@ LinearProjectionConstraint::setProjection(const size_t& nodeId, const Vec3d& p, m_projection = Mat3d::Identity() - p*p.transpose() - q*q.transpose(); } +void +LinearProjectionConstraint::setProjectionToLine(const size_t& nodeId, const Vec3d& p) +{ + m_nodeId = nodeId; + auto v = p/p.norm(); + m_projection = v*v.transpose(); +} + void LinearProjectionConstraint::setProjectorToDirichlet(const unsigned int& nodeId, const Vec3d z) { diff --git a/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.h b/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.h index 4fbdab7ac4f75fc033ddc3960b9fc03b7258151d..09e24f08f1d576fc8e481b7d9d2c93c7a0260c53 100644 --- a/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.h +++ b/Source/Constraint/PbdConstraints/imstkLinearProjectionConstraint.h @@ -50,6 +50,11 @@ public: /// void setProjection(const size_t& nodeId, const Vec3d& p, const Vec3d& q = Vec3d::Zero()); + /// + /// \brief Form the projection + /// + void setProjectionToLine(const size_t& nodeId, const Vec3d& p); + /// /// \brief Set the projector to simulate Dirichlet conditions /// diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp index c6ab72f4363aa4ba27841ef2becd1ca8364ca072..e5bdee5ccba294f070c91dc0eb6d9ba136aad741 100644 --- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp +++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp @@ -33,6 +33,7 @@ namespace imstk FEMDeformableBodyModel::FEMDeformableBodyModel() : DynamicalModel(DynamicalModelType::elastoDynamics) { + m_fixedNodeIds.reserve(1000); } void diff --git a/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp b/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp index ed5be3e5ec042e66f797c455ffc8c98b98b9bd17..73c2258d15fd991b887cdf78b266d8943370abb8 100644 --- a/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp +++ b/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp @@ -80,7 +80,8 @@ VTKRenderer::VTKRenderer(std::shared_ptr<Scene> scene, const bool enableVR) auto axes = vtkSmartPointer<vtkAxesActor>::New(); axes->SetShaftType(vtkAxesActor::CYLINDER_SHAFT); axes->SetAxisLabels(false); - m_debugVtkActors.push_back( axes ); + axes->SetTotalLength(40, 40, 40); + m_debugVtkActors.push_back(axes); // Camera and camera actor if (!enableVR) diff --git a/Source/SimulationManager/imstkSimulationManager.cpp b/Source/SimulationManager/imstkSimulationManager.cpp index f2c2f4ab900f63c5cd66994ca77ddf7818e5ddbb..36c2176ca1cbc84a3d35ebb6ed53dd9e6b292b5e 100644 --- a/Source/SimulationManager/imstkSimulationManager.cpp +++ b/Source/SimulationManager/imstkSimulationManager.cpp @@ -450,6 +450,11 @@ SimulationManager::runSimulation() LOG(INFO) << "Running simulation"; } + if (!m_simThreadLaunched) + { + this->launchSimulation(); + } + // Run scene m_sceneManagerMap.at(m_activeSceneName)->run(); diff --git a/Source/Solvers/imstkConjugateGradient.cpp b/Source/Solvers/imstkConjugateGradient.cpp index 39f4865b9b5665295c3b30f69c51da04416af6d7..8834f30d65f9c08f2af0a0e8feae48b55d30476f 100644 --- a/Source/Solvers/imstkConjugateGradient.cpp +++ b/Source/Solvers/imstkConjugateGradient.cpp @@ -21,6 +21,8 @@ #include "imstkConjugateGradient.h" +#include <iostream> + namespace imstk { ConjugateGradient::ConjugateGradient() @@ -38,19 +40,14 @@ ConjugateGradient::ConjugateGradient(const SparseMatrixd& A, const Vectord& rhs) void ConjugateGradient::applyLinearProjectionFilter(Vectord& x, const std::vector<LinearProjectionConstraint>& linProj, const bool setVal) { - Vec3d p; for (auto &localProjector : linProj) { const auto threeI = 3 * localProjector.getNodeId(); - p = Vec3d(x(threeI), x(threeI + 1), x(threeI + 2)); + Vec3d p = localProjector.getProjector()*Vec3d(x(threeI), x(threeI + 1), x(threeI + 2)); - if (!setVal) - { - p = localProjector.getProjector()*p; - } - else + if (setVal) { - p = (Mat3d::Identity() - localProjector.getProjector())*localProjector.getValue(); + p += (Mat3d::Identity() - localProjector.getProjector())*localProjector.getValue(); } x(threeI) = p.x(); @@ -68,7 +65,7 @@ ConjugateGradient::solve(Vectord& x) return; } - if (m_FixedLinearProjConstraints->size() == 0) + if (!(m_FixedLinearProjConstraints || m_DynamicLinearProjConstraints)) { x = m_cgSolver.solve(m_linearSystem->getRHSVector()); } @@ -86,12 +83,25 @@ ConjugateGradient::modifiedCGSolve(Vectord& x) // Set the initial guess to zero x.setZero(); - applyLinearProjectionFilter(x, *m_DynamicLinearProjConstraints, true); - applyLinearProjectionFilter(x, *m_FixedLinearProjConstraints, true); + if (m_DynamicLinearProjConstraints) + { + applyLinearProjectionFilter(x, *m_DynamicLinearProjConstraints, true); + } + + if (m_FixedLinearProjConstraints) + { + applyLinearProjectionFilter(x, *m_FixedLinearProjConstraints, true); + } auto res = b; - applyLinearProjectionFilter(res, *m_DynamicLinearProjConstraints, false); - applyLinearProjectionFilter(res, *m_FixedLinearProjConstraints, false); + if (m_DynamicLinearProjConstraints) + { + applyLinearProjectionFilter(res, *m_DynamicLinearProjConstraints, false); + } + if (m_FixedLinearProjConstraints) + { + applyLinearProjectionFilter(res, *m_FixedLinearProjConstraints, false); + } auto c = res; auto delta = res.dot(res); auto deltaPrev = delta; @@ -104,8 +114,14 @@ ConjugateGradient::modifiedCGSolve(Vectord& x) while (delta > eps) { q = A * c; - applyLinearProjectionFilter(q, *m_DynamicLinearProjConstraints, false); - applyLinearProjectionFilter(q, *m_FixedLinearProjConstraints, false); + if (m_DynamicLinearProjConstraints) + { + applyLinearProjectionFilter(q, *m_DynamicLinearProjConstraints, false); + } + if (m_FixedLinearProjConstraints) + { + applyLinearProjectionFilter(q, *m_FixedLinearProjConstraints, false); + } dotval = c.dot(q); if (dotval != 0.0) { @@ -122,12 +138,18 @@ ConjugateGradient::modifiedCGSolve(Vectord& x) delta = res.dot(res); c *= delta / deltaPrev; c += res; - applyLinearProjectionFilter(c, *m_DynamicLinearProjConstraints, false); - applyLinearProjectionFilter(c, *m_FixedLinearProjConstraints, false); + if (m_DynamicLinearProjConstraints) + { + applyLinearProjectionFilter(c, *m_DynamicLinearProjConstraints, false); + } + if (m_FixedLinearProjConstraints) + { + applyLinearProjectionFilter(c, *m_FixedLinearProjConstraints, false); + } if (++iterNum >= m_maxIterations) { - LOG(WARNING) << "ConjugateGradient::modifiedCGSolve - The solver did not converge after max. iterations"; + //LOG(WARNING) << "ConjugateGradient::modifiedCGSolve - The solver did not converge after max. iterations"; break; } } diff --git a/Source/Solvers/imstkConjugateGradient.h b/Source/Solvers/imstkConjugateGradient.h index cf16aafdbbc6013546b1e4a9f5fa2dfbaed85afc..544150ddc3df64678e772df7d6d7ba6d1bfca286 100644 --- a/Source/Solvers/imstkConjugateGradient.h +++ b/Source/Solvers/imstkConjugateGradient.h @@ -98,6 +98,7 @@ public: /// void applyLinearProjectionFilter(Vectord& x, const std::vector<LinearProjectionConstraint>& linProj, const bool setVal); + /// /// \brief Get the vector denoting the filter /// void setLinearProjectors(std::vector<LinearProjectionConstraint>* f)