diff --git a/CMake/External/External_iMSTKData.cmake b/CMake/External/External_iMSTKData.cmake
index c3592fbab4df201fe05d05c0b236d4a5cf9d8593..d849218257774f21519b021cf13e54f8609a577f 100644
--- a/CMake/External/External_iMSTKData.cmake
+++ b/CMake/External/External_iMSTKData.cmake
@@ -15,15 +15,6 @@ set(copy_data_command
${CMAKE_INSTALL_PREFIX}/data
)
-# HS - Due to an issue where it seems that the repository does not get updated
-# we need to call fetch here to refresh the index, this way the checkout can
-# succeed for the add_external_project call
-execute_process(
- COMMAND ${GIT_EXECUTABLE} fetch --all
- WORKING_DIRECTORY ${iMSTKData_PREFIX}/src
- RESULT_VARIABLE error_code
- )
-
include(imstkAddExternalProject)
imstk_add_external_project( iMSTKData
GIT_REPOSITORY "https://gitlab.kitware.com/iMSTK/imstk-data.git"
diff --git a/CMake/Utilities/imstkAddTest.cmake b/CMake/Utilities/imstkAddTest.cmake
index 69afec7d181f8bfcd20c70f1b22d9882ca948ed5..ec40c5ebb24678eb432ae8c42a4f15ff0ddce8aa 100644
--- a/CMake/Utilities/imstkAddTest.cmake
+++ b/CMake/Utilities/imstkAddTest.cmake
@@ -38,13 +38,6 @@ function(imstk_add_test_internal target kind)
file(GLOB test_files "${CMAKE_CURRENT_SOURCE_DIR}/*Test.h"
"${CMAKE_CURRENT_SOURCE_DIR}/*Test.cpp")
- # Get all source file names
- set(test_file_names "")
- foreach(test_file ${test_files})
- get_filename_component(test_file_name ${test_file} NAME)
- list(APPEND test_file_names ${test_file_name})
- endforeach()
-
# Create test driver executable
imstk_add_executable(${test_driver_executable} ${test_files})
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9f4e1b1e6e71e31c34c9732bc30a9fcdf29ec289..1d119e5e9d044addd7f7e1d638afb48b40d8d4e7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -345,14 +345,16 @@ endif(Uncrustify_EXECUTABLE)
#----------------------------------------------------------------------
if (${PROJECT_NAME}_BUILD_TESTING OR ${PROJECT_NAME}_BUILD_EXAMPLES)
add_definitions( -DiMSTK_DATA_ROOT=\"${CMAKE_INSTALL_PREFIX}/data/\")
-endif()
-
-# Note if the target name or data name changes this has to be changed as well
-add_custom_target(CopyDataFiles ALL
+
+ # Note if the target name or data name changes this has to be changed as well
+ add_custom_target(CopyDataFiles ALL
COMMAND ${CMAKE_COMMAND} -E copy_directory
${CMAKE_BINARY_DIR}/../External/iMSTKData/src/Data
${CMAKE_INSTALL_PREFIX}/data
)
+endif()
+
+
#-----------------------------------------------------------------------------
# Testing
diff --git a/Docs/AddingDependency.md b/Docs/AddingDependency.md
deleted file mode 100644
index b994f8659a250cc79ac40c0ecdf5ddd434d82931..0000000000000000000000000000000000000000
--- a/Docs/AddingDependency.md
+++ /dev/null
@@ -1,96 +0,0 @@
-# Adding a new Dependency to iMSTK
-
-In this document the new library will be called `NewLib`. Depending on the use case some of the following steps may be skipped, but that is unlikely.
-
-### Superbuild
-
-Its the superbuilds job to download all the noted dependency, build and install them. Once that is done it will configure and build the Innerbuild.
-
-While a lot of projects at kitware use a superbuild architecture there are slight differences between various implementations, very likely most or all of the following steps are applicable to other setups but they may or may not happen in the same plce as iMSTK
-
-### Innerbuild
-
-The innerbuild is the actual `iMSTK` code, it expects all dependencies to be built, installed, and findable.
-
-## CMakeLists.txt
-
-In the `CMakeLists.txt` in the root directory. While we are using `${PROJECT_NAME}` to reflect the `iMSTK` project in this file, everywhere else `iMSTK` will be more appropriate.
-
----
-
-Add a Setting `${PROJECT_NAME}_USE_NewLib` e.g.
-
- option(${PROJECT_NAME}_USE_NewLib "Build with NewLib support" OFF)
-
-The state of the option (`ON`/`OFF`) should reflect whether the new library is should be built by default. When adding a new required dependency the option may be omitted, in that case omit the `if()` check for all the other steps
-
-This enables us to control and check for the availability of the library, any section that is only to be performed when `NewLib` is being built should be checked via `if(${PROJECT_NAME}_USE_NewLib)` or `if(iMSTK_USE_NewLIb)`
-
----
-
-Define `NewLib` as a dependency
-
- if (${PROJECT_NAME}_USE_NewLib)
- imstk_define_dependency(NewLib)
- endif()
-
-This tells the superbuild architecture that `iMSTK` depends on `NewLib` and will trigger all the other processes.
-
----
-
-Add `find_package` with the appropriate options for this library
-
- if (${PROJECT_NAME}_USE_VRPN)
- find_package( NewLib REQUIRED )
- endif()
-
-This is actually executed in the innerbuild of `iMSTK` it will enable the innerbuild to find the components of the library. Please note that `find_package` should executed is the local `FindNewLib.cmake` that will be created in a later step, and not the Cmake global script or the find scripts of the library that is being built.
-
-## CMake\External\CMakeLists.txt
-
-This file adds the project that makes up the innerbuild to the superbuild, any variables that need to be passed from the superbuild to the innerbuild need to be set here. This means you will probably want to add
-
- -D${PROJECT_NAME}_USE_NewLib:BOOL=${${PROJECT_NAME}_USE_NewLib}
-
-to the section in `ExternalProject_Add`, so that the appropriate option is visible inside the innerbuild.
-
-Additionally to enable CMake to find the library correctly _if_ the library already provides a `Config.cmake` or a `Find` pass the path to the library into the innerbuild, so use
-
- -DNewLib_DIR:PATH=${NewLib_DIR}
-
-## `CMake/External/External_NewLib.cmake`
-
-You will need to create this file, this is what describes what files to download and how to build them to support your new library to the superbuild. In general this will mean customizing `imstk_add_external_project`.
-
- include(imstkAddExternalProject)
- imstk_add_external_project(NewLib
- ....
- )
-
-and customizing the build options for the new library by passing them via the `CMAKE_CACHE_ARGS` section. e.g.
-
- -DBUILD_TESTING:BOOL=OFF
-
-There are quite a few examples for this in `iMSTK` now, best is to start simple and extend depending on the needs of the new library. `imstk_add_external_project` can be found in `CMake/Utilities/imstkAddExternalProject.cmake`
-
-## `CMake/FindNewLib.cmake`
-
-In case `NewLib` _does not_ provide a cmake `NewLibConfig.cmake` or its own find you will need to create `FindNewLib.cmake` in the `CMake` directory, this will be used during the innerbuild configuration step to initialize the include directories and library files related to the new library. The general pattern is like this. `imstkFind` can be found in `CMake/Utilities`
-
-
- include(imstkFind)
-
- # Locate the header newlib.h inside the newlib subdirectory
- imstk_find_header(NewLib newlib.h newlib)
- # Find liba and add it to the libraries for NewLib
- imstk_find_libary(NewLib liba)
- # Finish the process
- imstk_find_package(NewLib)
-
-## `CMake/iMSTKConfig.cmake.in`
-
-To expose the state of any variables that you set in the superbuild to project that depend on `iMSTK` you need to store the state inside this file. At installation time, this will be written out and is used to correctly restore the state at build time. e.g.
-
- set(iMSTK_USE_NewLib @iMSTK_USE_NewLib@)
-
-any `find_package` commands issued in the main file have to be replicated here are well, so that users of `iMSTK` can access all of `iMSTK`s dependencies.
diff --git a/Docs/Maintainance.md b/Docs/Maintainance.md
new file mode 100644
index 0000000000000000000000000000000000000000..1173e387b4de4df678c4e1943819a76e796cf687
--- /dev/null
+++ b/Docs/Maintainance.md
@@ -0,0 +1,134 @@
+# iMSTK Maintenance
+
+The following sections describe various iMSTK maintenance tasks
+
+# Adding a new iMSTK Dependency
+
+For clarity the new library will be called `NewLib`. Depending on the use case some of the following steps may be skipped, but that is unlikely.
+
+### Superbuild
+
+iMSTK's CMake-based build process proceeds in two stages. The first stage downloads, builds and installs all of iMSTK's dependencies. This stage is called 'Superbuild'. Superbuild allows developers to build iMSTK and its dependencies in one go. When a new dependency needs to be added, the dependency list needs to be updated so that the superbuild builds and installs it for iMSTK to use. This is described in sections below.
+
+### Innerbuild
+
+iMSTK's code gets build in the second stage of the build process called 'Innerbuild'. The innerbuild follows the superbuild and expects all dependencies to be built, installed, and findable.
+
+## Updating the CMakeLists.txt
+
+In the `CMakeLists.txt` in the root directory. While we are using `${PROJECT_NAME}` to reflect the `iMSTK` project in this file, everywhere else `iMSTK` will be more appropriate.
+
+
+**Step 1:** Define the external library as a iMSTK dependency
+
+---
+ if (${PROJECT_NAME}_USE_NewLib)
+ imstk_define_dependency(NewLib)
+ endif()
+
+This tells the superbuild that `iMSTK` depends on `NewLib` and will trigger all the other processes namely building, installing, finding and linking to it.
+
+
+**Step 2:** Add `find_package` with the appropriate options for this library
+
+---
+
+ if (${PROJECT_NAME}_USE_VRPN)
+ find_package( NewLib REQUIRED )
+ endif()
+
+This is executed in the innerbuild of `iMSTK`. It will enable the innerbuild to find the components of the library. Please note that `find_package` should be executed in the local `FindNewLib.cmake` that will be created in a later step, and not the CMake top-level script or the find scripts of the library that is being built.
+
+**Step 3:** Add an option to turn the building of dependency ON/OFF (Optional)
+
+---
+Add a Setting `${PROJECT_NAME}_USE_NewLib` as
+
+ option(${PROJECT_NAME}_USE_NewLib "Build with NewLib support" OFF)
+
+The state of the option (`ON`/`OFF`) should reflect whether the new library should be built by default or not. If a dependency is a required one, one may omit this. If the dependency is optional, one can conditionally execute other steps in the CMake build process by surrounding the statements with `if(${PROJECT_NAME}_USE_NewLib)` or `if(iMSTK_USE_NewLIb)`.
+
+**Step 4:** Edit CMake\External\CMakeLists.txt
+
+---
+
+This CMake script adds the project that makes up the innerbuild to the superbuild, any variables that need to be passed from the superbuild to the innerbuild need to be set here. This means you will probably want to add
+
+ -D${PROJECT_NAME}_USE_NewLib:BOOL=${${PROJECT_NAME}_USE_NewLib}
+
+to the section in `ExternalProject_Add`, so that the appropriate option is visible inside the innerbuild.
+
+Additionally to enable CMake to find the library correctly _if_ the library already provides a `Config.cmake` or a `Find` pass the path to the library into the innerbuild, so use
+
+ -DNewLib_DIR:PATH=${NewLib_DIR}
+
+**Step 5:** Add `CMake/External/External_NewLib.cmake`
+
+You will need to create this file which describes what files to download from where and how to build them to support your new library in the superbuild. In general this will mean customizing `imstk_add_external_project`.
+
+ include(imstkAddExternalProject)
+ imstk_add_external_project(NewLib
+ ....
+ )
+
+and customizing the build options for the new library by passing them via the `CMAKE_CACHE_ARGS` section. e.g.
+
+ -DBUILD_TESTING:BOOL=OFF
+
+There are quite a few examples for this in `iMSTK` now, best is to start simple and extend depending on the needs of the new library. `imstk_add_external_project` can be found in `CMake/Utilities/imstkAddExternalProject.cmake`
+
+**Step 6:** `CMake/FindNewLib.cmake`
+
+In case `NewLib` _does not_ provide a cmake `NewLibConfig.cmake` or its own find you will need to create `FindNewLib.cmake` in the `CMake` directory, this will be used during the innerbuild configuration step to initialize the include directories and library files related to the new library. The general pattern is like this. `imstkFind` can be found in `CMake/Utilities`
+
+
+ include(imstkFind)
+
+ # Locate the header newlib.h inside the newlib subdirectory
+ imstk_find_header(NewLib newlib.h newlib)
+ # Find liba and add it to the libraries for NewLib
+ imstk_find_libary(NewLib liba)
+ # Finish the process
+ imstk_find_package(NewLib)
+
+**Step 7:** Edit `CMake/iMSTKConfig.cmake.in`
+
+To expose the state of any variables that you set in the superbuild to project that depend on `iMSTK` you need to store the state inside this file. At installation time, this will be written out and is used to correctly restore the state at build time. e.g.
+
+ set(iMSTK_USE_NewLib @iMSTK_USE_NewLib@)
+
+Any `find_package` commands issued in the main file have to be replicated here are well so that users of `iMSTK` can access all of its dependencies.
+
+# Updating a Dependency in iMSTK
+
+## Updating a dependency in imstk
+
+ - Swap the url in ` /CMake/External/<DependencyLibraryName>.cmake` to a different url.
+ - Ex: `https://gitlab.kitware.com/iMSTK/assimp/-/archive/fixCompilationError/assimp-fixCompilationError.zip`. Pulling zips prevents git history with it which can be sizeable.
+ - Alternatively git shallow would be used to avoid large git histories. But it fails when the tree gets too large, thus if the HEAD of the branch pulled from moves too far, it will fail to pull. (used to cause old version of iMSTK to fail to build until we switched to zips)
+ - Update the md5 hashsum in that same file (if using zips).
+ - The hashsum is a nearly unique identifer generated from the files of the dependency. This mostly is a security measure that ensures you are getting the files you expect to be getting.
+ - To acquire the md5 hashsum one can build iMSTK with the newly updated url. It will fail on the md5 check and report both the actual and current md5.
+
+## Updating a Remote Fork
+
+Most dependencies in iMSTK are forked. This way we don't depend on the remote repository as it could change (rebased/amended/moved/deleted). A few forks also contain our own diffs in the rare case something like a CMake fix is introduced. To update a fork:
+ - Clone the fork locally: `git clone <git url of fork>` (all forks found in iMSTK group here: https://gitlab.kitware.com/iMSTK)
+ - Add upstream `git remote add upstream <url of actual repo/repo forked off>` (the description of the fork normally provides what it was forked from)
+ - Tip: Issue `git remote -v` to list all remotes
+ - Merge upstream (or rebase): `git merge upstream/<branch to update from>`
+ - Push your changes
+
+After updating your fork you can proceed with the beginning of this guide on how to pull a different source.
+
+# Adding data to the data repository
+
+iMSTKs data sits in a separate repository https://gitlab.kitware.com/iMSTK/imstk-data. It is downloaded by iMSTK's superbuild when either the testing or the examples is enabled. To add new data to the repository two steps are necessary
+
+**1. Add Data to repository**
+
+The repository is already checked out as an external dependency in your build directory (if using superbuild) as `<build_dir>\External\iMSTKData\src\Data`, data can be added to the folder here and directly commited and pushed.
+
+**2. Update SHA in `ExternalData.cmake`**
+
+The file that controls the downloading of the data is `CMake/External/External_iMSTKData.cmake`. After commiting and pushing a change in the `iMSTKData` repository, the SHA to checkout needs to be updated to the SHA matching the latest commit.
diff --git a/Docs/UpdatingDependency.md b/Docs/UpdatingDependency.md
deleted file mode 100644
index 8399e1856af512883d5bf884e5efd33e0171a3ca..0000000000000000000000000000000000000000
--- a/Docs/UpdatingDependency.md
+++ /dev/null
@@ -1,21 +0,0 @@
-# Updating a Dependency in iMSTK
-
-To update a dependency in imstk one needs to:
-
- - Swap the url in ` /CMake/External/<DependencyLibraryName>.cmake` to a different url.
- - Ex: `https://gitlab.kitware.com/iMSTK/assimp/-/archive/fixCompilationError/assimp-fixCompilationError.zip`. Pulling zips prevents git history with it which can be sizeable.
- - Alternatively git shallow would be used to avoid large git histories. But it fails when the tree gets too large, thus if the HEAD of the branch pulled from moves too far, it will fail to pull. (used to cause old version of iMSTK to fail to build until we switched to zips)
- - Update the md5 hashsum in that same file (if using zips).
- - The hashsum is a nearly unique identifer generated from the files of the dependency. This mostly is a security measure that ensures you are getting the files you expect to be getting.
- - To acquire the md5 hashsum one can build iMSTK with the newly updated url. It will fail on the md5 check and report both the actual and current md5.
-
-## Updating a Remote Fork
-
-Most dependencies in iMSTK are forked. This way we don't depend on the remote repository as it could change (rebased/amended/etc). A few forks also contain our own diffs in the rare case something like a cmake fix is introduced. To update a fork:
- - Clone the fork locally: `git clone <git url of fork>` (all forks found in iMSTK group here: https://gitlab.kitware.com/iMSTK)
- - Add upstream `git remote add upstream <url of actual repo/repo forked off>` (the description of the fork normally provides what it was forked from)
- - Tip: `git remote -v` to list all remotes
- - Merge upstream (or rebase): `git merge upstream/<branch to update from>`
- - Push your changes
-
-After updating your fork you can proceed with the beginning of this guide on how to pull a different source.
\ No newline at end of file
diff --git a/Docs/source/Dynamical_Models.rst b/Docs/source/Dynamical_Models.rst
index b260d434067c0bb32ee4ddc76986d353bf9e34d2..a0c1ebd58d5d94fba969de46f4d33671df1bc957 100644
--- a/Docs/source/Dynamical_Models.rst
+++ b/Docs/source/Dynamical_Models.rst
@@ -6,85 +6,10 @@ Overview
========
Dynamical Models, model dynamical systems which describe a function in time. Generally, this is given via some dynamical equation (usually a PDE) that can be discretized for various geometries. Commonly used geometries include tetrahedral, hexahedral, triangle meshes, pointsets, and implicit geometries.
-Position Based Dynamics (PBD)
+:doc:`Position Based Dynamics (PBD) <../PbdModel>`.
=============================
-The first is position based dynamics which is fast and stable constraint based simulation method [pbd]_. It's fast and stable because it directly solves for positions. That is, each constraint is formulated with the positions. For example, a distance constraint may try to maintain a scalar distance L between two positions. While some models may compute a velocity that would converge the position to L. PBD directly computes the position to converge to L. For a single constraint this has a very simple solution. But when multiple constraints are involved, they may begin to conflict with each other, so we reproject the constraints to quantify how satisfied a constraint is. In PBD velocities are then computed aftewards given the resulting positions. The full list of PBD constraints may be found here.
-
-xPBD
-=====
-With original PBD, constraints become very stiff as the timestep decreases. Extended PBD (xPBD) alleviates this by removing that dependency. xPBD is the default in iMSTKs PbdModel.
-
-Usage
-=====
-Position based dynamics is very flexible. If you can come up with a constraint for it, you can simulate it. This includes things such as fluids, deformables, solids. This allows:
-
-- Blood (pointsets)
-
-.. image:: media/blood.png
- :width: 400
- :alt: Alternative text
- :align: center
-
-
-- 2d/thin tissues and cloth (triangle meshes)
-
-.. image:: media/cloth1.png
- :width: 400
- :alt: Cloth simulation
- :align: center
-
-.. image:: media/cloth2.png
- :width: 400
- :alt: Cloth cutting
- :align: center
-
-- Soft tissues and organs (tetrahedral meshes)
-
-.. image:: media/heart2.png
- :width: 400
- :alt: Tetrahedral mesh heart
- :align: center
-
-.. image:: media/vess.gif
- :width: 400
- :alt: VESS
- :align: center
-
-- String (linemeshes, useful for suturing)
-
-.. image:: media/strings.png
- :width: 400
- :alt: Strings
- :align: center
-
-
-Code
-====
-To setup a PbdModel we do:
-
-::
-
- // Setup the config
- imstkNew<PBDModelConfig> pbdConfig;
-
- // Constraints
- pbdConfig->enableConstraint(PbdConstraint::Type::Distance, 1e2);
- pbdConfig->enableConstraint(PbdConstraint::Type::Dihedral, 1e1);
- pbdConfig->m_fixedNodeIds = { 0, 1 };
-
- // Other parameters
- pbdConfig->m_uniformMassValue = 1.0;
- pbdConfig->m_gravity = Vec3d(0, -9.8, 0);
- pbdConfig->m_defaultDt = 0.005;
- pbdConfig->m_iterations = 10;
-
- // Setup the model
- imstkNew<PbdModel> pbdModel;
- pbdModel->setModelGeometry(surfMesh);
- pbdModel->configure(pbdConfig);
-
-This can be given to a PbdObject for usage in the scene.
+PBD is probably the most common model used in iMSTK, it may be used for cloths, thin tissues, volumetric/tetrahedral tissues, fluids (both liquids and gasses), and threads/strings. Read more about it at the link above.
Smoothed Particle Hydrodynamics (SPH)
=====================================
@@ -337,14 +262,6 @@ Bibliography
conjugate gradient method in cloth simulation. Vis. Comput. 19, 7-8
(December 2003), 526-531.
-.. [pbd] Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John
- Ratcliff. 2007. Position based dynamics. J. Vis. Comun. Image
- Represent. 18, 2 (April 2007), 109-118.
-
-.. [xpbd] Miles Macklin, Matthias Müller, and Nuttapong Chentanez
- 2016. XPBD: position-based simulation of compliant constrained dynamics.
- In Proc. of Motion in Games. 49–54
-
.. [vrpn] Russell M. Taylor, II, Thomas C. Hudson, Adam Seeger, Hans Weber,
Jeffrey Juliano, and Aron T. Helser. 2001. VRPN: a device-independent,
network-transparent VR peripheral system. In Proceedings of the ACM
diff --git a/Docs/source/PbdModel.rst b/Docs/source/PbdModel.rst
new file mode 100644
index 0000000000000000000000000000000000000000..91f94ea60dd9474f2409e46b64822065202930bf
--- /dev/null
+++ b/Docs/source/PbdModel.rst
@@ -0,0 +1,146 @@
+Position Based Dynamics
+===================
+
+Position based dynamics [pbd]_, as used in iMSTK, is a first order, particle & constraint based dynamical model. It simulates the dynamics of objects through direct manipulation of particle positions with velocities computed afterwards. This has various pro's and cons, but a major pro is the stability of such method. Pbd in iMSTK is primarily used for deformable bodies but can also be used for fluids (both gas and liquid). Discussed here are its medical applications. The particular implementation of PBD used in iMSTK is Extended PBD [xpbd]_.
+
+The general pbd pipeline is as follows:
+
+- **Integrate Position**: Given the current velocity of every particle, compute the new/tentative position.
+- **Solve Internal Constraints**: Solve all internal constraints (directly changing positions). Resolving any violations.
+- **Update Velocities**: Compute new velocities from the displacements produced.
+
+Constraints
+--------------------------------
+
+Pbd is a constraint based model. Constraints are made for particles in pbd. This constrains the movement of a particle. Constraints are given via a constraint function q, and the gradient of the function q. To solve a constraint is to reduce the scalar, q, to 0. The gradient gives the direction to move the particle to do this. The following constraints are available in iMSTK and referred to in later sections on how to use them:
+
+There are two types of constraints in iMSTK. PbdConstraints and PbdCollisionConstraints. PbdConstraints are generally inernal constraints and work by particle index. They often exist for the entire duration of a simulation. PbdCollisionConstraints are often external, they are added at runtime upon contact, they work by particle pointer value (position and velocity). Among the PbdConstraints are:
+
+- **PbdDistanceConstraint**: Constraints two points by the initial distance between them.
+- **PbdDihedralConstraint**: Constrains two triangles by the initial angle between their planes.
+- **PbdVolumeConstraint**: Constrains all points of a tetrahedron by the initial volume.
+- **PbdFETetConstraint**:
+- **PbdBendConstraint**: Constrains 3 points by the angle between the two lines. (like a 1d Dihedral)
+- **PbdConstantDensityConstraint**: Constrains all points of a system to achieve constant density using SPH kernels (ie: repels particles in a spherical manner).
+- **PbdAreaConstraint**: Constrains 3 points of a triangle by the initial area of that triangle.
+
+To implement your own custom constraint one needs to subclass either ``PbdConstraint`` or ``PbdCollisionConstraint``, construct or initialize it appropriately, then implement/override ``computeValueAndGradient`` which should fill out the constraint value C and constraint gradient dC/dx. This constraint would then need to be added to an existing system or solved in it's own computational block somewhere (depending on when it should be solved).
+
+Deformable Membranes
+--------------------------------
+
+For thin deformable membraneous tissues (cloth like) we use surfaces made of triangles. For this simulation one may use **PbdDistanceConstraint**, **PbdDihedralConstraint**, and/or **PbdAreaConstraint**. The most common setup is just distance and dihedral constraints though. The distance constraints keep the cloth together, whilst the dihedral constraint controls how resistent to bending the membrane is.
+
+.. image:: media/PbdModel/tissue1.png
+ :width: 400
+ :alt: Cloth simulation
+ :align: center
+
+.. image:: media/PbdModel/cloth2.png
+ :width: 400
+ :alt: Cloth cutting
+ :align: center
+
+Deformable Volumetric Tissue
+--------------------------------
+
+For volumetric deformable tissues discretized with tetrahedrons may be used. With the tetrahedrons one may either use (a) **PbdVolumeConstraint** & **PbdDistanceConstraint** constraints Or (b) use **PbdFemConstraint** constraints. The FEM constraints are more accurate than the volume+distance. However, they are much slower in that one may not be able to achieve the target element count or timestep, iteration count, & stiffness. The volume constraints behave well with sign inversion, recovering well from inverted tetrahedrons.
+
+.. image:: media/PbdModel/tissue2.png
+ :width: 400
+ :alt: Cloth simulation
+ :align: center
+
+.. image:: media/PbdModel/heart2.png
+ :width: 400
+ :alt: Cloth simulation
+ :align: center
+
+.. image:: media/PbdModel/tissue3.gif
+ :width: 400
+ :alt: Cloth simulation
+ :align: center
+
+Deformable Threads
+--------------------------------
+
+Surture threads are very common in surgical scenarios. For threads one may use **PbdDistanceConstraint** & **PbdBendConstraint** constraints. The distance constraints keep the particles of the thread together, whilst the bend controls the rigidity of the the thread. The bend constraints may also be generated between multiple sets of particles to reduce iteration count.
+
+.. image:: media/PbdModel/thread1.png
+ :width: 400
+ :alt: Cloth simulation
+ :align: center
+
+Liquids
+--------------------------------
+
+Liquids can be modeled with pbd using **PbdConstantDensityConstraint**. Generally, the stiffness is kept as high as possible as liquids are incompressible. If not, you may observe "bouncey" behaviour. Liquids in iMSTK are most useful for bleeding simulation.
+
+.. image:: media/PbdModel/blood.png
+ :width: 400
+ :alt: Cloth simulation
+ :align: center
+..
+ Gasses
+ --------------------------------
+
+ The primary usage for gas is particles during electrocautery. Often these would be billboarded smoke images on particles that fade fairly quickly. There are currently no examples for gas in iMSTK. It is a fluid though, so its approach is not much different than liquids. The **PbdConstantDensityConstraint** may be used. I would suggest using a lower stiffness as liquids tend to be incompressible (constant density) whereas gasses are compressible. The other issue is the lack of proper boundary conditions. Often we are modeling a gas suspended in air. This air must be modeled too if you want accuracy. There do exist some solutions with "ghost particles" to approximate air without adding air particles, but iMSTK does not have such solutions yet. If this is for visual purposes I might suggest lowering gravity, fiddling with mass, etc to get believable behaviour without being suspended in anything.```
+
+Code
+====
+To setup a PbdModel we do:
+
+::
+
+ // Setup the config
+ imstkNew<PBDModelConfig> pbdConfig;
+
+ // Constraints
+ pbdConfig->enableConstraint(PbdConstraint::Type::Distance, 1e2);
+ pbdConfig->enableConstraint(PbdConstraint::Type::Dihedral, 1e1);
+ pbdConfig->m_fixedNodeIds = { 0, 1 };
+
+ // Other parameters
+ pbdConfig->m_uniformMassValue = 1.0;
+ pbdConfig->m_gravity = Vec3d(0, -9.8, 0);
+ pbdConfig->m_dt = 0.005;
+ pbdConfig->m_iterations = 10;
+
+ // Setup the model
+ imstkNew<PbdModel> pbdModel;
+ pbdModel->setModelGeometry(surfMesh);
+ pbdModel->configure(pbdConfig);
+
+**Constraints**: Constraints of varying types may be used via ``PBDModelConfig::enableConstraint``, internally this uses PbdConstraintFunctor's which defines how to generate constraints. If one needs more hands on with constraints you may write your own ``PbdConstraintFunctor``. Implemented by subclassing PbdConstraintFunctor and overriding the operator() function. See existing functors in ``imstkPbdConstraintFunctor.h``.
+
+::
+
+ auto myCustomFunctor = std::make_shared<MySuperCustomFunctor>();
+ myCustomFunctor->setStiffness(0.95);
+ pbdModel->addPbdConstraintFunctor(myCustomFunctor);
+
+**Fixed Node Ids**: This indicates the indices of the particles/nodes that are immovable. Immovable nodes have their inverse masses set to 0 which indicates infinite mass (hence immovable).
+
+**Uniform Mass Value**: This mass value is assigned to all particles/nodes on initialization if masses are not provided as a vertex attribute in the input mesh.
+
+**dt**: The timestep is used during integration to move the particles. Small timesteps are preferable for stability. Real time steps may be used by varying dt every update of the simulation.
+
+::
+
+ connect<Event>(sceneManager, &SceneManager::postUpdate, [&](Event*)
+ {
+ pbdConfig->m_dt = sceneManager->getDt();
+ });
+
+**Iterations**: The iterations of the solver used in the internal constraints. More iterations give changes more time to percolate through the body. For example, a really long thread with numerous segments may have a really high stiffness but if it doesn't have enough iterations it will never be able to reach maximum stiffness. In the original PBD paper stiffness varied with the number of iterations. In xPBD (default) it does not.
+
+Bibliography
+------------
+
+.. [pbd] Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John
+ Ratcliff. 2007. Position based dynamics. J. Vis. Comun. Image
+ Represent. 18, 2 (April 2007), 109-118.
+
+.. [xpbd] Miles Macklin, Matthias Müller, and Nuttapong Chentanez
+ 1. XPBD: position-based simulation of compliant constrained dynamics.
+ In Proc. of Motion in Games. 49–54
\ No newline at end of file
diff --git a/Docs/source/documentation.rst b/Docs/source/documentation.rst
index bca17516657ce7ad3b899acaec8e02d32af8d9ec..cbd44a8483c8282930d78165649d5d21768faeeb 100644
--- a/Docs/source/documentation.rst
+++ b/Docs/source/documentation.rst
@@ -140,31 +140,23 @@ The following diagram gives a high level view of iMSTK architecture. Most of the
:width: 700
:alt: iMSTK architecture
-Geometry
---------
-
+:doc:`Geometry <../Geometry>`
+-------------------
Geometries are used throughout iMSTK. The geometries themselves can be used standalone though. We have analytical, implicit, surface, volumetric, and image geometries.
-Read more about iMSTK geometry :doc:`here <../Geometry>`.
-
-
-Dynamical Models
-----------------
+:doc:`Dynamical Models <../Dynamical_Models>`
+-------------------
Geometries are great. But you're likely going to want to do something with it. The primary use case is to advance it in time via some dynamical model. In iMSTK we provide models for dynamic and static deformable/soft bodies, fluids, & rigid bodies. We include PBD, SPH, FEM, and Rigid Bodies.
-Read more about dynamical models :doc:`here <../Dynamical_Models>`.
-
-Geometric Filtering
+:doc:`Geometric Filtering <../Filtering>`
-------------------
What else can you do with geometries? Filtering! Our filtering library provides a set of geometry algorithms we have found useful for surgical simulations.
-Read more about filtering :doc:`here <../Filtering>`.
-
-Devices
---------
+:doc:`Devices <../Devices>`
+-------------------
Devices are an important part to iMSTK. This is an the interactive surgical simulation toolkit after all.
@@ -173,82 +165,54 @@ Devices are an important part to iMSTK. This is an the interactive surgical simu
- OpenHaptics: Allows one to use haptic tracking devices such as the Phantom Omni, these provide force feedback, under the right models we can stop your hand from moving when touching something, or give slight resistance.
- Coming Soon: VRPN, Bluetooth, Arduino ...
-
-Read more about iMSTK's device support :doc:`here <../Devices>`.
-
-
-Controllers
------------
+:doc:`Controllers <../Controllers>`.
+-------------------
Controllers implement the controls of a device. We provide a couple of abstract base classes such as MouseControl, KeyboardControl, TrackingDeviceControl. As well as a few subclasses such as KeyboardSceneControl and MouseSceneControl which have some default behaviors such as stopping, starting, pausing a scene. But it's heavily encouraged you to subclass your own. You may also use lambdas on the devices for fast prototyping.
-Read more about controllers :doc:`here <../Controllers>`.
-
-
-Collision Detection
+:doc:`Collision Detection <../Collision_Detection>`
-------------------
Collision detection can be standalone in iMSTK but often finds it use through Interactions, later described in Scene. Put simply, its the act of computing "CollisionData" from two geometries. Most of the time these are "contacts" such as point normal contacts which give a point, a normal, and penetration depth. Often these are then later given to constraints to be added to a dynamical model.
-Read more about iMSTK's collision detection support :doc:`here <../Collision_Detection>`.
-
-
-Collision Handling
+:doc:`Collision Handling <../Collision_Handling>`
------------------
Collision handling implements how to consume collision data. For this reason it takes input CollisionData which is generally shared with CollisionDetection. iMSTK provides a number of handling methods, generally these call upon the functions of a DynamicalModel to immediately respond (explicit solve) or add something (such as a constraint) to later implicitly solve.
-Read more about iMSTK's collision handling support :doc:`here <../Collision_Handling>`.
-
-
-Scene
---------
+:doc:`Scene <../Scene>`
+-------------------
A scene contains a flat collection of SceneObjects and can fully represent the virtual environment. These SceneObjects may be something like an OR table, a tissue, a leg, a light, or even non-visual objects. Additionally a scene contains a set of interactions via it's InteractionGraph. A number of predefined InteractionPairs are available for iMSTK physics.
-Read more about iMSTK scene :doc:`here <../Scene>`.
-
-
-Mesh IO
---------
-Geometries are great. But to fully leverage them you need to be able to import from other tools which are much better at creating them.
-
-Read more about the geometry file types supported in iMSTK. Additionally about Scene and SceneObject importing :doc:`here <../Mesh_IO>`.
-
+:doc:`Mesh IO <../Mesh_IO>`
+-------------------
+Geometries are great. But to fully leverage them you need to be able to import from other tools which are much better at creating them. Read more about the files types supported by iMSTK. Additionally about Scene and SceneObject at the link above.
-SimulationManager & Modules
----------------------------
+:doc:`SimulationManager & Modules <../SimManager_Modules>`
+-------------------
Whilst scene's define the virtual environment and how to update it. They don't define how to drive it. You can certainly just call advance on the scene in a loop. That will get you decently far, but there's a bit more too it than that.
Modules in iMSTK define something that can be init'd, update'd, and uninit'd and added to a ModuleDriver. In every iMSTK example you can simply add modules to our concrete ModuleDriver called the SimulationManager to run them. It defines a special way of updating them.
-Read more about the SimulationManager in iMSTK with code samples :doc:`here <../SimManager_Modules>`.
-
-Rendering
----------
+:doc:`Rendering <../Rendering>`
+-------------------
Rendering in iMSTK is done through delegation to support multiple backends. This means we setup delegate classes for each thing we want to render. And map what we want to render to what the backend allows us to render. Primarily we use VTK.
-Read more about Rendering in iMSTK with code samples :doc:`here <../Rendering>`.
-
Miscellaneous Topics
====================
-Parallelism
------------
-Visit :doc:`Parallelism page <../Parallelism>` which goes over loop, task, and module parallelism in iMSTK.
-
-
-Events
-------
-Visit :doc:`iMSTK events Page <../Event_System>` which goes over events. Direct and message queues.
-
-
-Computational Flow
-------------------
-Visit :doc:`Computational Flow Page <../Computational_Flow>` which goes over the flow of the advancement of a scene.
+:doc:`Parallelism <../Parallelism>`
+-------------------
+Goes over loop, task, and module parallelism in iMSTK.
+:doc:`Events <../Event_System>`
+-------------------
+Goes over events. Direct and message queues.
-
+:doc:`Computational Flow <../Computational_Flow>`
+-------------------
+Goes over the flow/advancement of a scene.
.. |image0| image:: media/logo.png
:width: 3.5in
diff --git a/Docs/source/media/blood.png b/Docs/source/media/PbdModel/blood.png
similarity index 100%
rename from Docs/source/media/blood.png
rename to Docs/source/media/PbdModel/blood.png
diff --git a/Docs/source/media/cloth.png b/Docs/source/media/PbdModel/cloth.png
similarity index 100%
rename from Docs/source/media/cloth.png
rename to Docs/source/media/PbdModel/cloth.png
diff --git a/Docs/source/media/cloth1.png b/Docs/source/media/PbdModel/cloth1.png
similarity index 100%
rename from Docs/source/media/cloth1.png
rename to Docs/source/media/PbdModel/cloth1.png
diff --git a/Docs/source/media/cloth2.png b/Docs/source/media/PbdModel/cloth2.png
similarity index 100%
rename from Docs/source/media/cloth2.png
rename to Docs/source/media/PbdModel/cloth2.png
diff --git a/Docs/source/media/heart2.png b/Docs/source/media/PbdModel/heart2.png
similarity index 100%
rename from Docs/source/media/heart2.png
rename to Docs/source/media/PbdModel/heart2.png
diff --git a/Docs/source/media/PbdModel/thread1.png b/Docs/source/media/PbdModel/thread1.png
new file mode 100644
index 0000000000000000000000000000000000000000..7008a4789d39bf55e353f6cfb39b915784dcfcb5
Binary files /dev/null and b/Docs/source/media/PbdModel/thread1.png differ
diff --git a/Docs/source/media/PbdModel/tissue1.png b/Docs/source/media/PbdModel/tissue1.png
new file mode 100644
index 0000000000000000000000000000000000000000..ee6ee00cb82aca487415199d1cfcacdfd8b0b79d
Binary files /dev/null and b/Docs/source/media/PbdModel/tissue1.png differ
diff --git a/Docs/source/media/PbdModel/tissue2.png b/Docs/source/media/PbdModel/tissue2.png
new file mode 100644
index 0000000000000000000000000000000000000000..4abd341468e0c8c5e39f9d287dd8e4f820fdf481
Binary files /dev/null and b/Docs/source/media/PbdModel/tissue2.png differ
diff --git a/Docs/source/media/vess.gif b/Docs/source/media/PbdModel/tissue3.gif
similarity index 100%
rename from Docs/source/media/vess.gif
rename to Docs/source/media/PbdModel/tissue3.gif
diff --git a/Examples/TaskGraph/Configuration/taskGraphConfigureExample.cpp b/Examples/TaskGraph/Configuration/taskGraphConfigureExample.cpp
index 977015b4cfb3d1f9bb7c748079b14b95d40ed52c..2f4f2d2f397bbbc60ccf584b8e1b19898a5db6be 100644
--- a/Examples/TaskGraph/Configuration/taskGraphConfigureExample.cpp
+++ b/Examples/TaskGraph/Configuration/taskGraphConfigureExample.cpp
@@ -188,13 +188,6 @@ main()
// This node computes displacements and sets the color to the magnitude
std::shared_ptr<TaskNode> computeVelocityScalars = std::make_shared<TaskNode>([&]()
{
- /*const StdVectorOfVec3d& initPos = clothGeometry->getInitialVertexPositions();
- const StdVectorOfVec3d& currPos = clothGeometry->getVertexPositions();
- StdVectorOfReal& scalars = *scalarsPtr;
- for (size_t i = 0; i < initPos.size(); i++)
- {
- scalars[i] = (currPos[i] - initPos[i]).norm();
- }*/
const VecDataArray<double, 3>& velocities =
*std::dynamic_pointer_cast<VecDataArray<double, 3>>(clothGeometry->getVertexAttribute("Velocities"));
DataArray<double>& scalars = *scalarsPtr;
diff --git a/Source/Animation/Particles/imstkRenderParticleEmitter.cpp b/Source/Animation/Particles/imstkRenderParticleEmitter.cpp
index 965af42874edd803352f05f2636307d0929f0678..0befcbb1072cb4fb3a74fdb5fac16962749583db 100644
--- a/Source/Animation/Particles/imstkRenderParticleEmitter.cpp
+++ b/Source/Animation/Particles/imstkRenderParticleEmitter.cpp
@@ -268,7 +268,7 @@ RenderParticleEmitter::emitParticle(std::unique_ptr<RenderParticle>& particle)
particle->m_position[2] = (float)position[2] + z;
}
- float randomRotation = getRandomNormalizedFloat() * (float)imstk::PI * 2.0f;
+ float randomRotation = getRandomNormalizedFloat() * (float)PI * 2.0f;
float randomRotationalVelocity = getRandomNormalizedFloat();
particle->m_rotation = randomRotation;
particle->m_rotationalVelocity = (randomRotationalVelocity * m_minRotationSpeed) +
diff --git a/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.cpp b/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.cpp
index e804dc09121395803187a2f46d20ff1a013318f0..77f3de9f48f1529e48f4d58081fa7335f40818d7 100644
--- a/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.cpp
+++ b/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.cpp
@@ -112,7 +112,7 @@ pointSegmentClosestDistance(const Vec3d& point, const Vec3d& x1, const Vec3d& x2
{
Vec3d dx = x2 - x1;
double m2 = dx.squaredNorm();
- if (m2 < Real(1e-20))
+ if (m2 < 1e-20)
{
return (point - x1).norm();
}
@@ -176,7 +176,7 @@ closestPointOnSegment(const Vec3d& point, const Vec3d& x1, const Vec3d& x2, int&
{
Vec3d dx = x2 - x1;
double m2 = dx.squaredNorm();
- if (m2 < Real(1e-20))
+ if (m2 < 1e-20)
{
caseType = 0;
return x1;
diff --git a/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.h b/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.h
index 29ab5173ec08dace3e8f9c1b5a628e79d1528f8a..fe8c971b63a4546d9f06174f2efc5095a1635271 100644
--- a/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.h
+++ b/Source/CollisionDetection/CollisionDetection/imstkCollisionUtils.h
@@ -104,14 +104,14 @@ bool testLineToLineAABB(const double x1, const double y1, const double z1,
/// \param prox2 Round-off precision for the test
///
inline bool
-testLineToLineAABB(const Vec3r& p1A, const Vec3r& p1B,
- const Vec3r& p2A, const Vec3r& p2B,
+testLineToLineAABB(const Vec3d& p1A, const Vec3d& p1B,
+ const Vec3d& p2A, const Vec3d& p2B,
const double prox1 = VERY_SMALL_EPSILON_D, const double prox2 = VERY_SMALL_EPSILON_D)
{
- const Real* p1Aptr = &p1A[0];
- const Real* p1Bptr = &p1B[0];
- const Real* p2Aptr = &p2A[0];
- const Real* p2Bptr = &p2B[0];
+ const double* p1Aptr = &p1A[0];
+ const double* p1Bptr = &p1B[0];
+ const double* p2Aptr = &p2A[0];
+ const double* p2Bptr = &p2B[0];
return testLineToLineAABB(p1Aptr[0], p1Aptr[1], p1Aptr[2],
p1Bptr[0], p1Bptr[1], p1Bptr[2],
p2Aptr[0], p2Aptr[1], p2Aptr[2],
@@ -153,9 +153,8 @@ testOBBToPoint(
const Vec3d& pt,
Vec3d& ptContactNormal, Vec3d& cubeContactPt, double& penetrationDepth)
{
- const Vec3d diff = (pt - cubePos);
- const Mat3d rotInv = rot.transpose();
- const Vec3d proj = rotInv * diff; // dot product on each axes
+ const Vec3d diff = (pt - cubePos);
+ const Vec3d proj = rot.transpose() * diff; // dot product on each axes
bool inside[3] =
{
@@ -172,7 +171,7 @@ testOBBToPoint(
for (int i = 0; i < 3; i++)
{
double dist = proj[i];
- const Vec3d& axes = rotInv.col(i);
+ const Vec3d& axes = rot.col(i);
if (dist < extents[i] && dist > 0.0)
{
@@ -205,7 +204,7 @@ testOBBToPoint(
for (int i = 0; i < 3; i++)
{
double dist = proj[i];
- const Vec3d& axes = rotInv.col(i);
+ const Vec3d& axes = rot.col(i);
// If distance farther than the box extents, clamp to the box
if (dist >= extents[i])
@@ -584,60 +583,41 @@ testPlaneLine(const Vec3d& p, const Vec3d& q,
}
///
-/// \brief Computes scalar triple product of three vectors
+/// \brief Compute intersection of triangle vs segment and returns triangle interpolation
+/// weights
///
-static double
-scalarTripleProduct(const Vec3d& u, const Vec3d& v, const Vec3d& w)
-{
- return u.cross(v).dot(w);
-}
-
static bool
testSegmentTriangle(
const Vec3d& p, const Vec3d& q,
const Vec3d& a, const Vec3d& b, const Vec3d& c,
Vec3d& uvw)
{
- const Vec3d pq = q - p;
- const double pqL = pq.norm();
- if (pqL < IMSTK_DOUBLE_EPS)
- {
- return false;
- }
- const Vec3d pa = a - p;
- const Vec3d pb = b - p;
- const Vec3d pc = c - p;
-
- // Test if pq is inside the edges bc, ca and ab. Done by testing
- // that the signed tetrahedral volumes, computed using scalar triple
- // products, are all positive
- uvw[0] = scalarTripleProduct(pq, pc, pb);
- if (uvw[0] < 0.0)
+ const Vec3d n = q - p;
+ const Vec3d planeNormal = (b - a).cross(c - a);
+ const double denom = n.dot(planeNormal);
+ // Plane and line are parallel
+ if (std::abs(denom) < IMSTK_DOUBLE_EPS)
{
return false;
}
- uvw[1] = scalarTripleProduct(pq, pa, pc);
- if (uvw[1] < 0.0)
+
+ const double t1 = (a - p).dot(planeNormal);
+ const double t2 = (a - q).dot(planeNormal);
+
+ // Check if p and q lie on opposites side of the plane
+ if ((t1 < 0.0 && t2 >= 0.0) || (t1 >= 0.0 && t2 < 0.0))
{
- return false;
+ // \todo: Does planeNormal need to be normalized to get valid p + t1 * n, given division by denom
+ //t1 /= denom;
+ //t2 /= denom;
+ uvw = baryCentric(p + t1 / denom * n, a, b, c);
+ // Lastly check if the point on the plane p+t1*n is inside the triangle
+ return (uvw[0] >= 0.0 && uvw[1] >= 0.0 && uvw[2] >= 0.0);
}
- uvw[2] = scalarTripleProduct(pq, pb, pa);
- if (uvw[2] < 0.0)
+ else
{
return false;
}
-
- // Compute the barycentric coordinates (u, v, w) determining the
- // intersection point r, r = u*a + v*b + w*c
- double denom = 1.0 / (uvw[0] + uvw[1] + uvw[2]);
- uvw[0] *= denom;
- uvw[1] *= denom;
- uvw[2] *= denom; // w = 1.0f - u - v;
-
- const Vec3d iPt = uvw[0] * a + uvw[1] * b + uvw[2] * c;
- const double dist = (pq / pqL).dot(iPt - p);
-
- return (dist >= 0.0 && dist <= pqL);
}
///
@@ -645,45 +625,35 @@ testSegmentTriangle(
///
static bool
testSegmentTriangle(
- const Vec3d& pA, const Vec3d& pB,
- const Vec3d& v0, const Vec3d& v1, const Vec3d& v2)
+ const Vec3d& p, const Vec3d& q,
+ const Vec3d& a, const Vec3d& b, const Vec3d& c)
{
- static const double EPSILON = 1e-8;
- Vec3d dirAB = pB - pA;
- const double lAB = dirAB.norm();
-
- if (lAB < 1e-8)
+ const Vec3d n = q - p;
+ const Vec3d planeNormal = (b - a).cross(c - a);
+ const double denom = n.dot(planeNormal);
+ // Plane and line are parallel
+ if (std::abs(denom) < IMSTK_DOUBLE_EPS)
{
return false;
}
- dirAB /= lAB;
- const Vec3d edge1 = v1 - v0;
- const Vec3d edge2 = v2 - v0;
- const Vec3d h = dirAB.cross(edge2);
- const double a = edge1.dot(h);
- if (a > -EPSILON && a < EPSILON)
- {
- return false; // This ray is parallel to this triangle.
- }
- const double f = 1.0 / a;
- const Vec3d s = pA - v0;
- const double u = f * s.dot(h);
- if (u < 0.0 || u > 1.0)
+ const double t1 = (a - p).dot(planeNormal);
+ const double t2 = (a - q).dot(planeNormal);
+
+ // Check if p and q lie on opposites side of the plane
+ if ((t1 < 0.0 && t2 >= 0.0) || (t1 >= 0.0 && t2 < 0.0))
{
- return false;
+ // \todo: Does planeNormal need to be normalized to get valid p + t1 * n, given division by denom
+ //t1 /= denom;
+ //t2 /= denom;
+ const Vec3d uvw = baryCentric(p + t1 / denom * n, a, b, c);
+ // Lastly check if the point on the plane p+t1*n is inside the triangle
+ return (uvw[0] >= 0.0 && uvw[1] >= 0.0 && uvw[2] >= 0.0);
}
- const Vec3d q = s.cross(edge1);
- const double v = f * dirAB.dot(q);
- if (v < 0.0 || u + v > 1.0)
+ else
{
return false;
}
- // At this stage we can compute t to find out where the intersection point is on the line.
- const double t = f * edge2.dot(q);
-
- // ray - triangle intersection
- return (t > EPSILON && t + EPSILON < lAB);
}
///
@@ -961,12 +931,12 @@ testRayToPlane(const Vec3d& rayOrigin, const Vec3d& rayDir,
///
/// \brief Compute closest distance from a point to a segment x1-x2
///
-Real pointSegmentClosestDistance(const Vec3d& point, const Vec3d& x1, const Vec3d& x2);
+double pointSegmentClosestDistance(const Vec3d& point, const Vec3d& x1, const Vec3d& x2);
///
/// \brief Compute closest distance from a point to a triangle x1-x2-x3
///
-Real pointTriangleClosestDistance(const Vec3d& point, const Vec3d& x1, const Vec3d& x2, const Vec3d& x3);
+double pointTriangleClosestDistance(const Vec3d& point, const Vec3d& x1, const Vec3d& x2, const Vec3d& x3);
///
/// \brief Given two triangles and their ids produce the vertex ids for edge-edge and
@@ -1038,4 +1008,4 @@ edgeToEdgeClosestPoints(
return caseType;
}
}
-}
\ No newline at end of file
+}
diff --git a/Source/CollisionDetection/CollisionDetection/imstkImplicitGeometryToPointSetCD.cpp b/Source/CollisionDetection/CollisionDetection/imstkImplicitGeometryToPointSetCD.cpp
index 63447fca3a3e2c10c8baaf406bd74c35f1a245fb..84583d2b77b713f79c88b22de7299e813fb90372 100644
--- a/Source/CollisionDetection/CollisionDetection/imstkImplicitGeometryToPointSetCD.cpp
+++ b/Source/CollisionDetection/CollisionDetection/imstkImplicitGeometryToPointSetCD.cpp
@@ -33,6 +33,7 @@ ImplicitGeometryToPointSetCD::ImplicitGeometryToPointSetCD()
{
setRequiredInputType<ImplicitGeometry>(0);
setRequiredInputType<PointSet>(1);
+ m_centralGrad.setDx(Vec3d(0.001, 0.001, 0.001));
}
void
diff --git a/Source/CollisionDetection/CollisionDetection/imstkTetraToPointSetCD.cpp b/Source/CollisionDetection/CollisionDetection/imstkTetraToPointSetCD.cpp
index 1bf274de2af362f9c6c622f4ef5f679d4fdca34b..9dd7a2aeecbbe85937e83a9228c66fed4d7c66ce 100644
--- a/Source/CollisionDetection/CollisionDetection/imstkTetraToPointSetCD.cpp
+++ b/Source/CollisionDetection/CollisionDetection/imstkTetraToPointSetCD.cpp
@@ -44,8 +44,8 @@ TetraToPointSetCD::computeCollisionDataAB(
m_hashTableB.clear();
m_hashTableB.insertPoints(*pointSet->getVertexPositions());
- constexpr const double eps = VERY_SMALL_EPSILON;
- const double eps2 = 1e-10;
+ constexpr const double eps = IMSTK_DOUBLE_EPS;
+ const double eps2 = 1.0e-10;
//const VecDataArray<double, 3>& verticesMeshA = *m_tetMesh->getVertexPositions();
std::shared_ptr<VecDataArray<double, 3>> verticesMeshBPtr = pointSet->getVertexPositions();
diff --git a/Source/CollisionHandling/imstkPBDCollisionHandling.cpp b/Source/CollisionHandling/imstkPBDCollisionHandling.cpp
index a152655a1f2310c587e61fba0255002125d77b40..358829666028fda73b67df0da5592fb554d380ce 100644
--- a/Source/CollisionHandling/imstkPBDCollisionHandling.cpp
+++ b/Source/CollisionHandling/imstkPBDCollisionHandling.cpp
@@ -68,7 +68,7 @@ getVertex(const CollisionElement& elem, const MeshSide& side)
{
if (side.m_mapPtr && side.m_mapPtr->getType() == GeometryMap::Type::OneToOne)
{
- ptId = side.m_mapPtr->getMapIdx(ptId);
+ ptId = side.m_mapPtr->getMapIdx(static_cast<size_t>(ptId));
}
results[0] = { &side.m_vertices[ptId], side.m_invMasses[ptId], &side.m_velocities[ptId] };
}
diff --git a/Source/CollisionHandling/imstkPBDPickingCH.cpp b/Source/CollisionHandling/imstkPBDPickingCH.cpp
index 77a873ff1cbc7c4c05e2a026bee486bf11ff4d3a..f8f1d271252b8679f3ed97cab539202805e131e7 100644
--- a/Source/CollisionHandling/imstkPBDPickingCH.cpp
+++ b/Source/CollisionHandling/imstkPBDPickingCH.cpp
@@ -140,7 +140,7 @@ PBDPickingCH::addPickConstraints(const std::vector<CollisionElement>& elements,
lock.unlock();
}
}
- }, 100);
+ }, elements.size() > 100);
}
void
diff --git a/Source/CollisionHandling/imstkRigidBodyCH.cpp b/Source/CollisionHandling/imstkRigidBodyCH.cpp
index f5ebfbecb3018d8da4662e97a2bbc1658be55bf2..7907cbef35a1c67ffb27082749dbfda44e05c9e3 100644
--- a/Source/CollisionHandling/imstkRigidBodyCH.cpp
+++ b/Source/CollisionHandling/imstkRigidBodyCH.cpp
@@ -134,18 +134,18 @@ RigidBodyCH::handleRbdRbdTwoWay(
const CollisionElement& colElemA = elementsA[i];
if (colElemA.m_type == CollisionElementType::PointDirection)
{
- const Vec3d& contactPt = colElemA.m_element.m_PointDirectionElement.pt;
const Vec3d& dir = colElemA.m_element.m_PointDirectionElement.dir;
const double depth = colElemA.m_element.m_PointDirectionElement.penetrationDepth;
+ const Vec3d& contactPt = colElemA.m_element.m_PointDirectionElement.pt + dir * depth;
addConstraint(rbdObjA, rbdObjB, contactPt, dir, depth);
}
else if (colElemA.m_type == CollisionElementType::PointIndexDirection)
{
// Doesn't support mapping yet
- const Vec3d& contactPt = (*geom->getVertexPositions())[colElemA.m_element.m_PointIndexDirectionElement.ptIndex];
const Vec3d& dir = colElemA.m_element.m_PointIndexDirectionElement.dir;
const double depth = colElemA.m_element.m_PointIndexDirectionElement.penetrationDepth;
+ const Vec3d& contactPt = (*geom->getVertexPositions())[colElemA.m_element.m_PointIndexDirectionElement.ptIndex] + dir * depth;
addConstraint(rbdObjA, rbdObjB, contactPt, dir, depth);
}
@@ -165,9 +165,9 @@ RigidBodyCH::handleRbdStaticOneWay(
const CollisionElement& colElem = elementsA[i];
if (colElem.m_type == CollisionElementType::PointDirection)
{
- const Vec3d& contactPt = colElem.m_element.m_PointDirectionElement.pt;
const Vec3d& dir = colElem.m_element.m_PointDirectionElement.dir;
const double depth = colElem.m_element.m_PointDirectionElement.penetrationDepth;
+ const Vec3d& contactPt = colElem.m_element.m_PointDirectionElement.pt + dir * depth;
addConstraint(rbdObj, contactPt, dir, depth);
}
@@ -175,9 +175,9 @@ RigidBodyCH::handleRbdStaticOneWay(
{
// Doesn't support mapping yet
auto geom = std::dynamic_pointer_cast<PointSet>(rbdObj->getCollidingGeometry());
- const Vec3d& contactPt = (*geom->getVertexPositions())[colElem.m_element.m_PointIndexDirectionElement.ptIndex];
const Vec3d& dir = colElem.m_element.m_PointIndexDirectionElement.dir;
const double depth = colElem.m_element.m_PointIndexDirectionElement.penetrationDepth;
+ const Vec3d& contactPt = (*geom->getVertexPositions())[colElem.m_element.m_PointIndexDirectionElement.ptIndex] + dir * depth;
addConstraint(rbdObj, contactPt, dir, depth);
}
diff --git a/Source/Common/Parallel/imstkParallelReduce.h b/Source/Common/Parallel/imstkParallelReduce.h
index 72450afa20cc4080c42fe1d5a035fd139172de24..fe661212a26b6e7689d37346acaaee0ff477fffb 100644
--- a/Source/Common/Parallel/imstkParallelReduce.h
+++ b/Source/Common/Parallel/imstkParallelReduce.h
@@ -66,8 +66,8 @@ public:
Vec2d getRange() const { return Vec2d(m_Min, m_Max); }
private:
- Real m_Min = std::numeric_limits<Real>::max();
- Real m_Max = std::numeric_limits<Real>::min();
+ double m_Min = IMSTK_DOUBLE_MAX;
+ double m_Max = IMSTK_DOUBLE_MIN;
const ContainerType& m_Data;
};
@@ -89,16 +89,16 @@ public:
{
for (size_t i = r.begin(); i != r.end(); ++i)
{
- Real mag2 = m_Data[i].squaredNorm();
+ double mag2 = m_Data[i].squaredNorm();
m_Result = m_Result > mag2 ? m_Result : mag2;
}
}
void join(MaxL2NormFunctor& pObj) { m_Result = m_Result > pObj.m_Result ? m_Result : pObj.m_Result; }
- Real getResult() const { return std::sqrt(m_Result); }
+ double getResult() const { return std::sqrt(m_Result); }
private:
- Real m_Result = 0;
+ double m_Result = 0.0;
const VecDataArray<double, 3>& m_Data;
};
@@ -147,26 +147,26 @@ public:
///
/// \brief Get the lower corner
///
- const Vec3r& getLowerCorner() const { return m_LowerCorner; }
+ const Vec3d& getLowerCorner() const { return m_LowerCorner; }
///
/// \brief Get the upper corner
///
- const Vec3r& getUpperCorner() const { return m_UpperCorner; }
+ const Vec3d& getUpperCorner() const { return m_UpperCorner; }
private:
- Vec3r m_LowerCorner = Vec3r(std::numeric_limits<Real>::max(),
- std::numeric_limits<Real>::max(),
- std::numeric_limits<Real>::max());
- Vec3r m_UpperCorner = Vec3r(-std::numeric_limits<Real>::max(),
- -std::numeric_limits<Real>::max(),
- -std::numeric_limits<Real>::max());
+ Vec3d m_LowerCorner = Vec3d(std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max());
+ Vec3d m_UpperCorner = Vec3d(-std::numeric_limits<double>::max(),
+ -std::numeric_limits<double>::max(),
+ -std::numeric_limits<double>::max());
const VecDataArray<double, 3>& m_Data;
};
///
/// \brief Find the maximum value of L2 norm from the input data array
///
-inline Real
+inline double
findMaxL2Norm(const VecDataArray<double, 3>& data)
{
MaxL2NormFunctor pObj(data);
@@ -178,7 +178,7 @@ findMaxL2Norm(const VecDataArray<double, 3>& data)
/// \brief Find the bounding box of a point set
///
inline void
-findAABB(const VecDataArray<double, 3>& points, Vec3r& lowerCorner, Vec3r& upperCorner)
+findAABB(const VecDataArray<double, 3>& points, Vec3d& lowerCorner, Vec3d& upperCorner)
{
AABBFunctor pObj(points);
tbb::parallel_reduce(tbb::blocked_range<size_t>(0, points.size()), pObj);
diff --git a/Source/Common/imstkMath.h b/Source/Common/imstkMath.h
index 65361b76a0aa237adf1eca662a958f820f8a4c52..87a8fa4c10c9aa64760aa468cd881942a0efd469 100644
--- a/Source/Common/imstkMath.h
+++ b/Source/Common/imstkMath.h
@@ -52,16 +52,6 @@ make_unique(Args&& ... args)
namespace imstk
{
-// Define Real type and dependent types
-using Real = double;
-using Vec2r = Eigen::Matrix<Real, 2, 1>;
-using Vec3r = Eigen::Matrix<Real, 3, 1>;
-using Vec4r = Eigen::Matrix<Real, 4, 1>;
-using StdVectorOfReal = std::vector<Real>;
-using StdVectorOfVec2r = std::vector<Vec2r, Eigen::aligned_allocator<Vec2r>>;
-using StdVectorOfVec3r = std::vector<Vec3r, Eigen::aligned_allocator<Vec3r>>;
-using StdVectorOfVec4r = std::vector<Vec4r, Eigen::aligned_allocator<Vec4r>>;
-
// 2D vector
// using Vec2f = Eigen::Vector2f;
// using Vec2d = Eigen::Vector2d;
@@ -152,23 +142,19 @@ using AffineTransform3d = Eigen::Affine3d;
#define WORLD_ORIGIN Vec3d::Zero()
/// Some commonly used math constants
-#define PI Real(3.14159265358979323846)
-#define PI_2 Real(1.57079632679489661923)
-#define PI_4 Real(0.785398163397448309616)
-#define INV_1_PI Real(0.318309886183790671538)
-#define INV_2_PI Real(0.636619772367581343076)
-#define TWO_OVER_SQRT_PI Real(1.12837916709551257390)
-#define SQRT2 Real(1.41421356237309504880)
-#define SQRT1_2 Real(0.707106781186547524401)
-#define NLOG_E Real(2.71828182845904523536)
-#define LOG2E Real(1.44269504088896340736)
-#define LOG10E Real(0.434294481903251827651)
-#define LN2 Real(0.693147180559945309417)
-#define LN10 Real(2.30258509299404568402)
-
-#define MAX_REAL std::numeric_limits<Real>::max()
-#define MIN_REAL std::numeric_limits<Real>::min()
-#define VERY_SMALL_EPSILON std::numeric_limits<Real>::epsilon()
+#define PI 3.14159265358979323846
+#define PI_2 1.57079632679489661923
+#define PI_4 0.785398163397448309616
+#define INV_1_PI 0.318309886183790671538
+#define INV_2_PI 0.636619772367581343076
+#define TWO_OVER_SQRT_PI 1.12837916709551257390
+#define SQRT2 1.41421356237309504880
+#define SQRT1_2 0.707106781186547524401
+#define NLOG_E 2.71828182845904523536
+#define LOG2E 1.44269504088896340736
+#define LOG10E 0.434294481903251827651
+#define LN2 0.693147180559945309417
+#define LN10 2.30258509299404568402
#define MAX_D std::numeric_limits<double>::max()
#define MIN_D std::numeric_limits<double>::min()
diff --git a/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.cpp b/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.cpp
index 93f6fbe92e2028b5d31f359e7d37ec2b3ed2cddd..c059819cba90dcf17adb63a56affea9d1c54e47b 100644
--- a/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.cpp
+++ b/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.cpp
@@ -28,13 +28,6 @@ void
PbdConstantDensityConstraint::initConstraint(const VecDataArray<double, 3>& initVertexPositions, const double)
{
const size_t numParticles = initVertexPositions.size();
-
- // constraint parameters
- // Refer: Miller, et al 2003 "Particle-Based Fluid Simulation for Interactive Applications."
- /// \todo Check if these numbers can be variable
- m_wPoly6Coeff = 315.0 / (64.0 * PI * pow(m_maxDist, 9));
- m_wSpikyCoeff = 15.0 / (PI * pow(m_maxDist, 6));
-
m_lambdas.resize(numParticles);
m_densities.resize(numParticles);
m_deltaPositions.resize(numParticles);
@@ -71,31 +64,6 @@ PbdConstantDensityConstraint::projectConstraint(const DataArray<double>& imstkNo
});
}
-double
-PbdConstantDensityConstraint::wPoly6(const Vec3d& pi, const Vec3d& pj)
-{
- double rLengthSqr = (Vec3d(pi - pj)).squaredNorm();
-
- return (rLengthSqr > m_maxDistSqr || rLengthSqr < 1e-20) ?
- 0 :
- m_wPoly6Coeff* pow(m_maxDistSqr - rLengthSqr, 3);
-}
-
-Vec3d
-PbdConstantDensityConstraint::gradSpiky(const Vec3d& pi, const Vec3d& pj)
-{
- Vec3d r = pi - pj;
- const double rLengthSqr = r.squaredNorm();
-
- if (rLengthSqr > m_maxDistSqr || rLengthSqr < 1e-20)
- {
- return Vec3d(0., 0., 0.);
- }
-
- const double rLength = std::sqrt(rLengthSqr);
- return r * (m_wSpikyCoeff * (-3.0) * (m_maxDist - rLength) * (m_maxDist - rLength));
-}
-
void
PbdConstantDensityConstraint::computeDensity(const Vec3d& pi,
const size_t index,
@@ -115,11 +83,13 @@ PbdConstantDensityConstraint::computeLambdaScalingFactor(const Vec3d& pi,
const size_t index,
const VecDataArray<double, 3>& positions)
{
- const double densityConstraint = (m_densities[index] / m_restDensity) - 1;
+ const double invRestDensity = 1.0 / m_restDensity;
+ const double densityConstraint = (m_densities[index] * invRestDensity) - 1.0;
double gradientSum = 0.0;
+
for (auto q : m_neighborList[index])
{
- gradientSum += gradSpiky(pi, positions[q]).squaredNorm() / m_restDensity;
+ gradientSum += gradSpiky(pi, positions[q]).squaredNorm() * invRestDensity;
}
m_lambdas[index] = densityConstraint / (gradientSum + m_relaxationParameter);
@@ -130,12 +100,12 @@ PbdConstantDensityConstraint::updatePositions(const Vec3d& pi,
const size_t index,
VecDataArray<double, 3>& positions)
{
- //Make sure the point is valid
- Vec3d gradientLambdaSum(0., 0., 0.);
+ // Make sure the point is valid
+ Vec3d gradientLambdaSum(0.0, 0.0, 0.0);
for (auto q : m_neighborList[index])
{
- double lambdasDiff = (m_lambdas[index] + m_lambdas[q]);
- Vec3d gradKernal = gradSpiky(pi, positions[q]);
+ const double lambdasDiff = (m_lambdas[index] + m_lambdas[q]);
+ const Vec3d gradKernal = gradSpiky(pi, positions[q]);
gradientLambdaSum += (gradKernal * lambdasDiff);
}
@@ -146,8 +116,10 @@ PbdConstantDensityConstraint::updatePositions(const Vec3d& pi,
void
PbdConstantDensityConstraint::setMaxNeighborDistance(const double dist)
{
- m_maxDist = dist;
- m_maxDistSqr = dist * dist;
+ m_maxDist = dist;
+ m_maxDistSqr = dist * dist;
+ m_wPoly6Coeff = 315.0 / (64.0 * PI * pow(m_maxDist, 9));
+ m_wSpikyCoeff = 15.0 / (PI * pow(m_maxDist, 6)) * -3.0;
if (m_NeighborSearcher)
{
m_NeighborSearcher->setSearchRadius(m_maxDist);
diff --git a/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.h b/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.h
index 80116996f335210e35a53a75aeacf4463318e12f..bd8d0e489522ed252e50a658fe6111faeed9c59c 100644
--- a/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.h
+++ b/Source/Constraint/PbdConstraints/imstkPbdConstantDensityConstraint.h
@@ -41,11 +41,7 @@ public:
///
PbdConstantDensityConstraint() : PbdConstraint()
{
- // constraint parameters
- // Refer: Miller, et al 2003 "Particle-Based Fluid Simulation for Interactive Applications."
- /// \todo Check if these numbers can be variable
- m_wPoly6Coeff = 315.0 / (64.0 * PI * pow(m_maxDist, 9));
- m_wSpikyCoeff = 15.0 / (PI * pow(m_maxDist, 6));
+ setMaxNeighborDistance(m_maxDist);
}
///
@@ -78,12 +74,36 @@ private:
///
/// \brief Smoothing kernel WPoly6 for density estimation
///
- double wPoly6(const Vec3d& pi, const Vec3d& pj);
+ inline double wPoly6(const Vec3d& pi, const Vec3d& pj) const
+ {
+ const double rLengthSqr = (Vec3d(pi - pj)).squaredNorm();
+ if (rLengthSqr > m_maxDistSqr || rLengthSqr < 1e-20)
+ {
+ return 0.0;
+ }
+ else
+ {
+ const double maxDiff = m_maxDistSqr - rLengthSqr;
+ return m_wPoly6Coeff * maxDiff * maxDiff * maxDiff;
+ }
+ }
///
/// \brief Gradient of density kernel
///
- Vec3d gradSpiky(const Vec3d& pi, const Vec3d& pj);
+ inline Vec3d gradSpiky(const Vec3d& pi, const Vec3d& pj) const
+ {
+ const Vec3d r = pi - pj;
+ const double rLengthSqr = r.squaredNorm();
+
+ if (rLengthSqr > m_maxDistSqr || rLengthSqr < 1e-20)
+ {
+ return Vec3d::Zero();
+ }
+
+ const double rLength = std::sqrt(rLengthSqr);
+ return r * (m_wSpikyCoeff * (m_maxDist - rLength) * (m_maxDist - rLength));
+ }
///
/// \brief
@@ -104,22 +124,21 @@ private:
/// \brief Set/Get rest density
///
void setDensity(const double density) { m_restDensity = density; }
- double getDensity() { return m_restDensity; }
+ double getDensity() const { return m_restDensity; }
///
/// \brief Set/Get max. neighbor distance
///
void setMaxNeighborDistance(const double dist);
- double getMaxNeighborDistance() { return m_maxDist; }
+ double getMaxNeighborDistance() const { return m_maxDist; }
///
/// \brief Set/Get neighbor search method
///
void setNeighborSearchMethod(NeighborSearch::Method method) { m_NeighborSearchMethod = method; }
- NeighborSearch::Method getNeighborSearchMethod() { return m_NeighborSearchMethod; }
+ NeighborSearch::Method getNeighborSearchMethod() const { return m_NeighborSearchMethod; }
private:
-
double m_wPoly6Coeff;
double m_wSpikyCoeff;
diff --git a/Source/Constraint/PbdConstraints/imstkPbdConstraint.cpp b/Source/Constraint/PbdConstraints/imstkPbdConstraint.cpp
index c23c0fc27acb2e5ec23db346d1ab6ecd67ce0ee5..c4374ffceda55a290d100f39020ead35803480ae 100644
--- a/Source/Constraint/PbdConstraints/imstkPbdConstraint.cpp
+++ b/Source/Constraint/PbdConstraints/imstkPbdConstraint.cpp
@@ -43,7 +43,7 @@ PbdConstraint::projectConstraint(const DataArray<double>& invMasses, const doubl
dcMidc += invMasses[m_vertexIds[i]] * m_dcdx[i].squaredNorm();
}
- if (dcMidc < VERY_SMALL_EPSILON)
+ if (dcMidc < IMSTK_DOUBLE_EPS)
{
return;
}
diff --git a/Source/Constraint/RigidBodyConstraints/imstkRbdConstraint.h b/Source/Constraint/RigidBodyConstraints/imstkRbdConstraint.h
index b91d969a47eadd831f75b7da0b11760ccd90781b..a8cfc24180d523d26687c6c0a606d55dcbaf9311 100644
--- a/Source/Constraint/RigidBodyConstraints/imstkRbdConstraint.h
+++ b/Source/Constraint/RigidBodyConstraints/imstkRbdConstraint.h
@@ -47,8 +47,6 @@ struct RigidBody
Vec3d m_prevForce = Vec3d(0.0, 0.0, 0.0);
- // Vec3d m_externalForce;
- //RigidBodyState2* m_state; // A RigidBody can only belong to one state
Vec3d* m_pos = nullptr;
Quatd* m_orientation = nullptr;
Vec3d* m_velocity = nullptr;
diff --git a/Source/Constraint/Testing/imstkPbdBendConstraintTest.cpp b/Source/Constraint/Testing/imstkPbdBendConstraintTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..93a6cb8dea175e787321cbf2c2b25022ae32e678
--- /dev/null
+++ b/Source/Constraint/Testing/imstkPbdBendConstraintTest.cpp
@@ -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.
+
+=========================================================================*/
+
+#include "gtest/gtest.h"
+
+#include "imstkPbdBendConstraint.h"
+
+using namespace imstk;
+
+///
+/// \brief Test that two connecting line segments unfold
+///
+TEST(imstkPbdBendConstraintTest, TestConvergence1)
+{
+ PbdBendConstraint constraint;
+
+ // Straight line upon initialization
+ VecDataArray<double, 3> vertices(3);
+ vertices[0] = Vec3d(0.0, 0.0, 0.0);
+ vertices[1] = Vec3d(0.5, 0.0, 0.0);
+ vertices[2] = Vec3d(1.0, 0.0, 0.0);
+ DataArray<double> invMasses(3);
+ invMasses[0] = 1.0;
+ invMasses[1] = 0.0; // Center doesn't move
+ invMasses[2] = 1.0;
+
+ constraint.initConstraint(vertices, 0, 1, 2, 1e20);
+
+ // Modify it so the line segments look like \/
+ vertices[0][1] = 0.1;
+ vertices[2][1] = 0.1;
+ for (int i = 0; i < 500; i++)
+ {
+ constraint.projectConstraint(invMasses, 0.01, PbdConstraint::SolverType::xPBD, vertices);
+ }
+
+ // Should resolve back to a flat line
+ EXPECT_NEAR(vertices[0][1], 0.0, IMSTK_DOUBLE_EPS);
+ EXPECT_NEAR(vertices[2][1], 0.0, IMSTK_DOUBLE_EPS);
+}
\ No newline at end of file
diff --git a/Source/Constraint/Testing/imstkPbdPointEdgeConstraintTest.cpp b/Source/Constraint/Testing/imstkPbdPointEdgeConstraintTest.cpp
index ee8b679fda5e46a94b5d65b5d97d6c61ef03b2a1..8affd2bc28238a9106e8e6c482e7035167c79d4d 100644
--- a/Source/Constraint/Testing/imstkPbdPointEdgeConstraintTest.cpp
+++ b/Source/Constraint/Testing/imstkPbdPointEdgeConstraintTest.cpp
@@ -25,34 +25,27 @@
using namespace imstk;
-///
-/// \brief TODO
-///
-class imstkPbdPointEdgeConstraintTest : public ::testing::Test
-{
-protected:
- PbdPointEdgeConstraint m_constraint;
-};
-
///
/// \brief Test that a point and edge meet on touching
///
-TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence1)
+TEST(imstkPbdPointEdgeConstraintTest, TestConvergence1)
{
+ PbdPointEdgeConstraint constraint;
+
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
Vec3d b = Vec3d(0.5, 0.0, 0.0);
Vec3d x = Vec3d(0.0, -1.0, 0.0);
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], IMSTK_DOUBLE_EPS);
@@ -63,22 +56,24 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence1)
///
/// \brief Test that a point and edge meet on touching
///
-TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence2)
+TEST(imstkPbdPointEdgeConstraintTest, TestConvergence2)
{
+ PbdPointEdgeConstraint constraint;
+
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
Vec3d b = Vec3d(0.5, 0.0, 0.0);
Vec3d x = Vec3d(0.0, 1.0, 0.0);
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], IMSTK_DOUBLE_EPS);
@@ -89,8 +84,10 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence2)
///
/// \brief Test that a point not within bounds of edge does not move (left)
///
-TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
+TEST(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
{
+ PbdPointEdgeConstraint constraint;
+
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
const Vec3d aInit = a;
Vec3d b = Vec3d(0.5, 0.0, 0.0);
@@ -100,14 +97,14 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
const Vec3d xInit = x;
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
EXPECT_EQ(xInit[0], x[0]);
@@ -126,8 +123,10 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
///
/// \brief Test that a point not within bounds of edge does not move (right)
///
-TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
+TEST(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
{
+ PbdPointEdgeConstraint constraint;
+
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
const Vec3d aInit = a;
Vec3d b = Vec3d(0.5, 0.0, 0.0);
@@ -137,14 +136,14 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
const Vec3d xInit = x;
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
EXPECT_EQ(xInit[0], x[0]);
@@ -158,17 +157,4 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
EXPECT_EQ(bInit[0], b[0]);
EXPECT_EQ(bInit[1], b[1]);
EXPECT_EQ(bInit[2], b[2]);
-}
-
-///
-/// \brief TODO
-///
-int
-imstkPbdPointEdgeConstraintTest(int argc, char* argv[])
-{
- // Init Google Test & Mock
- ::testing::InitGoogleTest(&argc, argv);
-
- // Run tests with gtest
- return RUN_ALL_TESTS();
-}
+}
\ No newline at end of file
diff --git a/Source/Constraint/Testing/imstkPbdPointPointConstraintTest.cpp b/Source/Constraint/Testing/imstkPbdPointPointConstraintTest.cpp
index e521198218a5c6a80f2500bca1aa498e9e323a94..17cba1daa330e1a52b6250b9e08e7af62d1a285a 100644
--- a/Source/Constraint/Testing/imstkPbdPointPointConstraintTest.cpp
+++ b/Source/Constraint/Testing/imstkPbdPointPointConstraintTest.cpp
@@ -25,44 +25,24 @@
using namespace imstk;
-///
-/// \brief TODO
-///
-class imstkPbdPointPointConstraintTest : public ::testing::Test
-{
-protected:
- PbdPointPointConstraint m_constraint;
-};
-
///
/// \brief Test that two points meet
///
-TEST_F(imstkPbdPointPointConstraintTest, TestConvergence1)
+TEST(imstkPbdPointPointConstraintTest, TestConvergence1)
{
+ PbdPointPointConstraint constraint;
+
Vec3d a = Vec3d(0.0, 0.0, 0.0);
Vec3d b = Vec3d(0.0, -1.0, 0.0);
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &a, 1.0, nullptr },
{ &b, 1.0, nullptr },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
ASSERT_EQ(a[1], b[1]);
-}
-
-///
-/// \brief TODO
-///
-int
-imstkPbdPointPointConstraintTest(int argc, char* argv[])
-{
- // Init Google Test & Mock
- ::testing::InitGoogleTest(&argc, argv);
-
- // Run tests with gtest
- return RUN_ALL_TESTS();
-}
+}
\ No newline at end of file
diff --git a/Source/Constraint/Testing/imstkPbdPointTriangleConstraintTest.cpp b/Source/Constraint/Testing/imstkPbdPointTriangleConstraintTest.cpp
index c998e1b8d4ae4be8d55c9778575b097a47ade836..2627cf165791caa416a4ba3506872eb7c25a2f4f 100644
--- a/Source/Constraint/Testing/imstkPbdPointTriangleConstraintTest.cpp
+++ b/Source/Constraint/Testing/imstkPbdPointTriangleConstraintTest.cpp
@@ -25,20 +25,13 @@
using namespace imstk;
-///
-/// \brief TODO
-///
-class imstkPbdPointTriangleConstraintTest : public ::testing::Test
-{
-protected:
- PbdPointTriangleConstraint m_constraint;
-};
-
///
/// \brief Test that a point below a triangle, and the triangle, meet on y axes
///
-TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
+TEST(imstkPbdPointTriangleConstraintTest, TestConvergence1)
{
+ PbdPointTriangleConstraint constraint;
+
Vec3d a = Vec3d(0.5, 0.0, -0.5);
Vec3d b = Vec3d(-0.5, 0.0, -0.5);
Vec3d c = Vec3d(0.0, 0.0, 0.5);
@@ -47,7 +40,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
x[1] -= 1.0;
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
@@ -55,7 +48,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], 0.00000001);
@@ -66,8 +59,10 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
///
/// \brief Test that a point above a triangle, and the triangle, meet on y axes
///
-TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
+TEST(imstkPbdPointTriangleConstraintTest, TestConvergence2)
{
+ PbdPointTriangleConstraint constraint;
+
Vec3d a = Vec3d(0.5, 0.0, -0.5);
Vec3d b = Vec3d(-0.5, 0.0, -0.5);
Vec3d c = Vec3d(0.0, 0.0, 0.5);
@@ -76,7 +71,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
x[1] += 1.0;
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
@@ -84,7 +79,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], 0.00000001);
@@ -95,8 +90,10 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
///
/// \brief Test that a point not within the triangle does not move at all
///
-TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
+TEST(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
{
+ PbdPointTriangleConstraint constraint;
+
Vec3d a = Vec3d(0.5, 0.0, -0.5);
const Vec3d aInit = a;
Vec3d b = Vec3d(-0.5, 0.0, -0.5);
@@ -120,7 +117,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
c = cInit;
Vec3d zeroVelocity = Vec3d::Zero();
- m_constraint.initConstraint(
+ constraint.initConstraint(
{ &testPts[i], 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
@@ -128,7 +125,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
1.0, 1.0);
for (int j = 0; j < 3; j++)
{
- m_constraint.solvePosition();
+ constraint.solvePosition();
}
// Test they haven't moved
@@ -148,17 +145,4 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
EXPECT_EQ(cInit[1], c[1]);
EXPECT_EQ(cInit[2], c[2]);
}
-}
-
-///
-/// \brief TODO
-///
-int
-imstkPbdPointTriangleConstraintTest(int argc, char* argv[])
-{
- // Init Google Test & Mock
- ::testing::InitGoogleTest(&argc, argv);
-
- // Run tests with gtest
- return RUN_ALL_TESTS();
-}
+}
\ No newline at end of file
diff --git a/Source/DataStructures/Testing/imstkLooseOctreeTest.cpp b/Source/DataStructures/Testing/imstkLooseOctreeTest.cpp
index aaa6740153fafcd9c714310bf80ffc579f3c5c69..d270140369250414553ebc15839420ab4f3ab966 100644
--- a/Source/DataStructures/Testing/imstkLooseOctreeTest.cpp
+++ b/Source/DataStructures/Testing/imstkLooseOctreeTest.cpp
@@ -29,9 +29,9 @@ using namespace imstk;
#define BOUND 10.0
-#define SPHERE_RADIUS Real(10)
-#define SPHERE_CENTER Vec3r(0, 0, 0)
-#define PARTICLE_RADIUS Real(2)
+#define SPHERE_RADIUS 10.0
+#define SPHERE_CENTER Vec3d::Zero()
+#define PARTICLE_RADIUS 2.0
#define ITERATIONS 10
///
@@ -40,15 +40,15 @@ using namespace imstk;
std::shared_ptr<PointSet>
generatePointSet()
{
- const Vec3r sphereCenter = SPHERE_CENTER;
+ const Vec3d sphereCenter = SPHERE_CENTER;
const auto sphereRadiusSqr = SPHERE_RADIUS * SPHERE_RADIUS;
- const auto spacing = Real(2) * PARTICLE_RADIUS;
+ const auto spacing = 2.0 * PARTICLE_RADIUS;
const int N = int(2 * SPHERE_RADIUS / spacing);
std::shared_ptr<VecDataArray<double, 3>> particles = std::make_shared<VecDataArray<double, 3>>();
VecDataArray<double, 3>& vertices = *particles;
vertices.reserve(N * N * N);
- const Vec3r corner = sphereCenter - Vec3r(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
+ const Vec3d corner = sphereCenter - Vec3d(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
for (int i = 0; i < N; ++i)
{
@@ -56,8 +56,8 @@ generatePointSet()
{
for (int k = 0; k < N; ++k)
{
- const Vec3r ppos = corner + Vec3r(spacing * Real(i), spacing * Real(j), spacing * Real(k));
- const Vec3r d = ppos - sphereCenter;
+ const Vec3d ppos = corner + Vec3d(spacing * static_cast<double>(i), spacing * static_cast<double>(j), spacing * static_cast<double>(k));
+ const Vec3d d = ppos - sphereCenter;
if (d.squaredNorm() < sphereRadiusSqr)
{
vertices.push_back(ppos);
@@ -106,10 +106,10 @@ randomizePositions(const std::shared_ptr<PointSet>& pointset)
{
for (int i = 0; i < pointset->getNumVertices(); ++i)
{
- pointset->setVertexPosition(i, Vec3r(
- (static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
- (static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
- (static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND
+ pointset->setVertexPosition(i, Vec3d(
+ (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND,
+ (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND,
+ (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND
));
}
pointset->postModified();
@@ -123,10 +123,10 @@ randomizePositions(const std::shared_ptr<SurfaceMesh>& mesh)
{
for (int i = 0; i < mesh->getNumTriangles(); ++i)
{
- const auto translation = Vec3r(
- (static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
- (static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND,
- (static_cast<Real>(rand()) / static_cast<Real>(RAND_MAX) * 2.0 - 1.0) * BOUND
+ const auto translation = Vec3d(
+ (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND,
+ (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND,
+ (static_cast<double>(rand()) / static_cast<double>(RAND_MAX) * 2.0 - 1.0) * BOUND
);
const Vec3i& face = (*mesh->getTriangleIndices())[i];
for (unsigned int j = 0; j < 3; ++j)
diff --git a/Source/DataStructures/Testing/imstkNeighborSearchTest.cpp b/Source/DataStructures/Testing/imstkNeighborSearchTest.cpp
index 90836310782ce46b9dd28fad86717823fd0af137..61cbd9b7445c81418347700366aef56797a89125 100644
--- a/Source/DataStructures/Testing/imstkNeighborSearchTest.cpp
+++ b/Source/DataStructures/Testing/imstkNeighborSearchTest.cpp
@@ -27,11 +27,11 @@
using namespace imstk;
-#define SPHERE_RADIUS Real(1)
-#define SPHERE_CENTER Vec3r(0, 0, 0)
-#define PARTICLE_RADIUS Real(0.08)
+#define SPHERE_RADIUS 1.0
+#define SPHERE_CENTER Vec3d::Zero()
+#define PARTICLE_RADIUS 0.08
#define ITERATIONS 5
-#define STEP Real(1.1)
+#define STEP 1.1
///
/// \brief Advance particle positions
@@ -41,8 +41,8 @@ advancePositions(VecDataArray<double, 3>& particles)
{
for (auto& pos: particles)
{
- Vec3r pc = pos - SPHERE_CENTER;
- Real mag = pc.norm() * STEP;
+ Vec3d pc = pos - SPHERE_CENTER;
+ double mag = pc.norm() * STEP;
pos = SPHERE_CENTER + pc.normalized() * mag;
}
}
@@ -54,8 +54,8 @@ void
neighborSearchBruteForce(VecDataArray<double, 3>& particles, std::vector<std::vector<size_t>>& neighbors)
{
neighbors.resize(particles.size());
- const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
- const Real radiusSqr = radius * radius;
+ const double radius = 4.000000000000001 * PARTICLE_RADIUS;
+ const double radiusSqr = radius * radius;
for (int p = 0; p < particles.size(); ++p)
{
@@ -71,7 +71,7 @@ neighborSearchBruteForce(VecDataArray<double, 3>& particles, std::vector<std::ve
}
const auto qpos = particles[q];
- const auto d2 = (Vec3r(ppos - qpos)).squaredNorm();
+ const auto d2 = (Vec3d(ppos - qpos)).squaredNorm();
if (d2 < radiusSqr)
{
pneighbors.push_back(q);
@@ -87,7 +87,7 @@ void
neighborSearchGridBased(VecDataArray<double, 3>& particles, std::vector<std::vector<size_t>>& neighbors)
{
neighbors.resize(particles.size());
- const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
+ const double radius = 4.000000000000001 * PARTICLE_RADIUS;
static GridBasedNeighborSearch gridSearch;
gridSearch.setSearchRadius(radius);
gridSearch.getNeighbors(neighbors, particles);
@@ -105,7 +105,7 @@ neighborSearchSpatialHashing(VecDataArray<double, 3>& particles, std::vector<std
list.resize(0);
}
- const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
+ const double radius = 4.000000000000001 * PARTICLE_RADIUS;
static SpatialHashTableSeparateChaining hashTable;
hashTable.clear();
@@ -126,8 +126,8 @@ void
neighborSearchBruteForce(VecDataArray<double, 3>& setA, VecDataArray<double, 3>& setB, std::vector<std::vector<size_t>>& neighbors)
{
neighbors.resize(setA.size());
- const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
- const Real radiusSqr = radius * radius;
+ const double radius = 4.000000000000001 * PARTICLE_RADIUS;
+ const double radiusSqr = radius * radius;
for (int p = 0; p < setA.size(); ++p)
{
@@ -138,7 +138,7 @@ neighborSearchBruteForce(VecDataArray<double, 3>& setA, VecDataArray<double, 3>&
for (int q = 0; q < setB.size(); ++q)
{
const auto qpos = setB[q];
- const auto d2 = (Vec3r(ppos - qpos)).squaredNorm();
+ const auto d2 = (Vec3d(ppos - qpos)).squaredNorm();
if (d2 < radiusSqr)
{
pneighbors.push_back(q);
@@ -154,7 +154,7 @@ void
neighborSearchGridBased(VecDataArray<double, 3>& setA, VecDataArray<double, 3>& setB, std::vector<std::vector<size_t>>& neighbors)
{
neighbors.resize(setA.size());
- const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
+ const double radius = 4.000000000000001 * PARTICLE_RADIUS;
static GridBasedNeighborSearch gridSearch;
gridSearch.setSearchRadius(radius);
gridSearch.getNeighbors(neighbors, setA, setB);
@@ -203,14 +203,14 @@ verify(std::vector<std::vector<size_t>>& neighbors1, std::vector<std::vector<siz
///
TEST(imstkNeighborSearchTest, CompareGridSearchAndSpatialHashing)
{
- const Vec3r sphereCenter = SPHERE_CENTER;
+ const Vec3d sphereCenter = SPHERE_CENTER;
const auto sphereRadiusSqr = SPHERE_RADIUS * SPHERE_RADIUS;
- const auto spacing = Real(2) * PARTICLE_RADIUS;
+ const auto spacing = 2.0 * PARTICLE_RADIUS;
const int N = int(2 * SPHERE_RADIUS / spacing);
VecDataArray<double, 3> particles;
particles.reserve(N * N * N);
- const Vec3r corner = sphereCenter - Vec3r(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
+ const Vec3d corner = sphereCenter - Vec3d(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
for (int i = 0; i < N; ++i)
{
@@ -218,8 +218,8 @@ TEST(imstkNeighborSearchTest, CompareGridSearchAndSpatialHashing)
{
for (int k = 0; k < N; ++k)
{
- const Vec3r ppos = corner + Vec3r(spacing * Real(i), spacing * Real(j), spacing * Real(k));
- const Vec3r d = ppos - sphereCenter;
+ const Vec3d ppos = corner + Vec3d(spacing * static_cast<double>(i), spacing * static_cast<double>(j), spacing * static_cast<double>(k));
+ const Vec3d d = ppos - sphereCenter;
if (d.squaredNorm() < sphereRadiusSqr)
{
particles.push_back(ppos);
@@ -249,14 +249,14 @@ TEST(imstkNeighborSearchTest, CompareGridSearchAndSpatialHashing)
///
TEST(imstkNeighborSearchTest, TestGridSearchFromDifferentPointSet)
{
- const Vec3r sphereCenter = SPHERE_CENTER;
+ const Vec3d sphereCenter = SPHERE_CENTER;
const auto sphereRadiusSqr = SPHERE_RADIUS * SPHERE_RADIUS;
- const auto spacing = Real(2) * PARTICLE_RADIUS;
+ const auto spacing = 2.0 * PARTICLE_RADIUS;
const int N = int(2 * SPHERE_RADIUS / spacing);
- StdVectorOfVec3r particles;
+ StdVectorOfVec3d particles;
particles.reserve(N * N * N);
- const Vec3r corner = sphereCenter - Vec3r(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
+ const Vec3d corner = sphereCenter - Vec3d(SPHERE_RADIUS, SPHERE_RADIUS, SPHERE_RADIUS);
for (int i = 0; i < N; ++i)
{
@@ -264,8 +264,8 @@ TEST(imstkNeighborSearchTest, TestGridSearchFromDifferentPointSet)
{
for (int k = 0; k < N; ++k)
{
- const Vec3r ppos = corner + Vec3r(spacing * Real(i), spacing * Real(j), spacing * Real(k));
- const Vec3r d = ppos - sphereCenter;
+ const Vec3d ppos = corner + Vec3d(spacing * static_cast<double>(i), spacing * static_cast<double>(j), spacing * static_cast<double>(k));
+ const Vec3d d = ppos - sphereCenter;
if (d.squaredNorm() < sphereRadiusSqr)
{
particles.push_back(ppos);
diff --git a/Source/DataStructures/imstkGridBasedNeighborSearch.cpp b/Source/DataStructures/imstkGridBasedNeighborSearch.cpp
index 527b7e85fdb4ba860fdde44c07052cd917c01523..f8851cd95963c68c73935052d0b5665c21b367d9 100644
--- a/Source/DataStructures/imstkGridBasedNeighborSearch.cpp
+++ b/Source/DataStructures/imstkGridBasedNeighborSearch.cpp
@@ -25,7 +25,7 @@
namespace imstk
{
void
-GridBasedNeighborSearch::setSearchRadius(const Real radius)
+GridBasedNeighborSearch::setSearchRadius(const double radius)
{
m_SearchRadius = radius;
m_SearchRadiusSqr = radius * radius;
@@ -48,15 +48,15 @@ GridBasedNeighborSearch::getNeighbors(std::vector<std::vector<size_t>>& result,
void
GridBasedNeighborSearch::getNeighbors(std::vector<std::vector<size_t>>& result, const VecDataArray<double, 3>& setA, const VecDataArray<double, 3>& setB)
{
- LOG_IF(FATAL, (std::abs(m_SearchRadius) < Real(1e-8))) << "Neighbor search radius is zero";
+ LOG_IF(FATAL, (std::abs(m_SearchRadius) < 1e-8)) << "Neighbor search radius is zero";
// firstly compute the bounding box of points in setB
- Vec3r lowerCorner;
- Vec3r upperCorner;
+ Vec3d lowerCorner;
+ Vec3d upperCorner;
ParallelUtils::findAABB(setB, lowerCorner, upperCorner);
// the upper corner need to be expanded a bit, to avoid round-off error during computation
- upperCorner += Vec3d(m_SearchRadius, m_SearchRadius, m_SearchRadius) * Real(0.1);
+ upperCorner += Vec3d(m_SearchRadius, m_SearchRadius, m_SearchRadius) * 0.1;
// resize grid to fit the bounding box covering setB
m_Grid.initialize(lowerCorner, upperCorner, m_SearchRadius);
@@ -117,7 +117,7 @@ GridBasedNeighborSearch::getNeighbors(std::vector<std::vector<size_t>>& result,
for (auto q : m_Grid.getCellData(cellX, cellY, cellZ).particleIndices)
{
const auto qpos = setB[q];
- const Vec3r diff = ppos - qpos;
+ const Vec3d diff = ppos - qpos;
const auto d2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
if (d2 < m_SearchRadiusSqr)
{
diff --git a/Source/DataStructures/imstkGridBasedNeighborSearch.h b/Source/DataStructures/imstkGridBasedNeighborSearch.h
index 0b3211e07a80585379468aa2ae3b4f34d8d1910b..a36e713476fff6a7461927e3c5fc1669796dfe70 100644
--- a/Source/DataStructures/imstkGridBasedNeighborSearch.h
+++ b/Source/DataStructures/imstkGridBasedNeighborSearch.h
@@ -39,18 +39,18 @@ public:
/// \brief Construct class with search radius
/// \param radius The search radius
///
- explicit GridBasedNeighborSearch(const Real radius) : m_SearchRadius(radius), m_SearchRadiusSqr(radius * radius) {}
+ explicit GridBasedNeighborSearch(const double radius) : m_SearchRadius(radius), m_SearchRadiusSqr(radius * radius) {}
///
/// \brief Set the search radius
/// \param radius The search radius
///
- void setSearchRadius(const Real radius);
+ void setSearchRadius(const double radius);
///
/// \brief Get the search radius
///
- Real getSearchRadius() const { return m_SearchRadius; }
+ double getSearchRadius() const { return m_SearchRadius; }
///
/// \brief Search neighbors for each points within the search radius
@@ -75,8 +75,8 @@ public:
void getNeighbors(std::vector<std::vector<size_t>>& result, const VecDataArray<double, 3>& setA, const VecDataArray<double, 3>& setB);
private:
- Real m_SearchRadius = 0.;
- Real m_SearchRadiusSqr = 0.;
+ double m_SearchRadius = 0.0;
+ double m_SearchRadiusSqr = 0.0;
// Data store in each grid cell
// This entire struct can be replaced by tbb::concurrent_vector<size_t>, however, with lower performance
diff --git a/Source/DataStructures/imstkLooseOctree.cpp b/Source/DataStructures/imstkLooseOctree.cpp
index 79652c1a907a0b58176c6cd4136274dff787db41..585227d073c3ee6e56853ccc236709d6ff933963 100644
--- a/Source/DataStructures/imstkLooseOctree.cpp
+++ b/Source/DataStructures/imstkLooseOctree.cpp
@@ -25,15 +25,15 @@
namespace imstk
{
-OctreeNode::OctreeNode(LooseOctree* const tree, OctreeNode* const pParent, const Vec3r& nodeCenter,
- const Real halfWidth, const uint32_t depth) :
+OctreeNode::OctreeNode(LooseOctree* const tree, OctreeNode* const pParent, const Vec3d& nodeCenter,
+ const double halfWidth, const uint32_t depth) :
m_pTree(tree),
m_pParent(pParent),
m_Center(nodeCenter),
- m_LowerBound(nodeCenter - Vec3r(halfWidth, halfWidth, halfWidth)),
- m_UpperBound(nodeCenter + Vec3r(halfWidth, halfWidth, halfWidth)),
- m_LowerExtendedBound(nodeCenter - 2.0 * Vec3r(halfWidth, halfWidth, halfWidth)),
- m_UpperExtendedBound(nodeCenter + 2.0 * Vec3r(halfWidth, halfWidth, halfWidth)),
+ m_LowerBound(nodeCenter - Vec3d(halfWidth, halfWidth, halfWidth)),
+ m_UpperBound(nodeCenter + Vec3d(halfWidth, halfWidth, halfWidth)),
+ m_LowerExtendedBound(nodeCenter - 2.0 * Vec3d(halfWidth, halfWidth, halfWidth)),
+ m_UpperExtendedBound(nodeCenter + 2.0 * Vec3d(halfWidth, halfWidth, halfWidth)),
m_HalfWidth(halfWidth),
m_Depth(depth),
m_MaxDepth(tree->m_MaxDepth),
@@ -104,7 +104,7 @@ OctreeNode::split()
m_pChildren = m_pTree->requestChildrenFromPool();
- const auto childHalfWidth = m_HalfWidth * static_cast<Real>(0.5);
+ const auto childHalfWidth = m_HalfWidth * static_cast<double>(0.5);
for (uint32_t childIdx = 0; childIdx < 8u; ++childIdx)
{
auto newCenter = m_Center;
@@ -219,7 +219,7 @@ OctreeNode::insertNonPointPrimitive(OctreePrimitive* const pPrimitive, const Oct
{
const auto lowerCorner = pPrimitive->m_LowerCorner;
const auto upperCorner = pPrimitive->m_UpperCorner;
- const Vec3r priCenter(
+ const Vec3d priCenter(
(lowerCorner[0] + upperCorner[0]) * 0.5,
(lowerCorner[1] + upperCorner[1]) * 0.5,
(lowerCorner[2] + upperCorner[2]) * 0.5);
@@ -278,14 +278,14 @@ OctreeNode::insertNonPointPrimitive(OctreePrimitive* const pPrimitive, const Oct
m_pChildren->m_Nodes[childIdx].insertNonPointPrimitive(pPrimitive, type);
}
-LooseOctree::LooseOctree(const Vec3r& center, const Real width, const Real minWidth,
- const Real minWidthRatio /*= 1.0*/, const std::string name /*= "LooseOctree"*/) :
+LooseOctree::LooseOctree(const Vec3d& center, const double width, const double minWidth,
+ const double minWidthRatio /*= 1.0*/, const std::string name /*= "LooseOctree"*/) :
m_Name(name),
m_Center(center),
m_Width(width),
m_MinWidthRatio(minWidthRatio),
m_MinWidth(minWidth),
- m_pRootNode(new OctreeNode(this, nullptr, center, width * static_cast<Real>(0.5), 1u)),
+ m_pRootNode(new OctreeNode(this, nullptr, center, width * static_cast<double>(0.5), 1u)),
m_NumAllocatedNodes(1u)
{
}
@@ -479,7 +479,7 @@ LooseOctree::build()
&& (m_vPrimitivePtrs[OctreePrimitiveType::Triangle].size() > 0
|| m_vPrimitivePtrs[OctreePrimitiveType::Analytical].size() > 0))
{
- Real minWidth = MAX_REAL;
+ double minWidth = IMSTK_DOUBLE_MAX;
for (int type = OctreePrimitiveType::Triangle; type <= OctreePrimitiveType::Analytical; ++type)
{
const auto& vPrimitivePtrs = m_vPrimitivePtrs[type];
@@ -488,14 +488,14 @@ LooseOctree::build()
continue;
}
const auto primitiveMinWidth = tbb::parallel_reduce(tbb::blocked_range<size_t>(0, vPrimitivePtrs.size()),
- MAX_REAL,
- [&](const tbb::blocked_range<size_t>& r, Real prevResult) -> Real {
+ IMSTK_DOUBLE_MAX,
+ [&](const tbb::blocked_range<size_t>& r, double prevResult) -> double {
for (auto i = r.begin(), iEnd = r.end(); i != iEnd; ++i)
{
const auto pPrimitive = vPrimitivePtrs[i];
computePrimitiveBoundingBox(pPrimitive, static_cast<OctreePrimitiveType>(type));
- Vec3r widths;
+ Vec3d widths;
for (uint32_t dim = 0; dim < 3; ++dim)
{
widths[dim] = pPrimitive->m_UpperCorner[dim] - pPrimitive->m_LowerCorner[dim];
@@ -507,7 +507,7 @@ LooseOctree::build()
}
return prevResult;
},
- [](const Real x, const Real y) -> Real {
+ [](const double x, const double y) -> double {
return x < y ? x : y;
});
@@ -528,13 +528,13 @@ LooseOctree::build()
m_MaxDepth = 1u;
uint32_t numLevelNodes = 1u;
uint32_t maxNumTreeNodes = 1u;
- Real nodeWidth = m_Width;
+ double nodeWidth = m_Width;
- while (nodeWidth * static_cast<Real>(0.5) > m_MinWidth) {
+ while (nodeWidth * 0.5 > m_MinWidth) {
++m_MaxDepth;
numLevelNodes *= 8u;
maxNumTreeNodes += numLevelNodes;
- nodeWidth *= static_cast<Real>(0.5);
+ nodeWidth *= 0.5;
}
m_pRootNode->m_MaxDepth = m_MaxDepth;
rebuild();
@@ -680,7 +680,7 @@ LooseOctree::updateBoundingBoxAndCheckValidity(const OctreePrimitiveType type)
computePrimitiveBoundingBox(pPrimitive, type);
const auto lowerCorner = pPrimitive->m_LowerCorner;
const auto upperCorner = pPrimitive->m_UpperCorner;
- const Vec3r center(
+ const Vec3d center(
(lowerCorner[0] + upperCorner[0]) * 0.5,
(lowerCorner[1] + upperCorner[1]) * 0.5,
(lowerCorner[2] + upperCorner[2]) * 0.5);
@@ -821,8 +821,8 @@ LooseOctree::computePrimitiveBoundingBox(OctreePrimitive* const pPrimitive, cons
<< "Cannot compute bounding box for point primitive";
#endif
- Vec3r lowerCorner;
- Vec3r upperCorner;
+ Vec3d lowerCorner;
+ Vec3d upperCorner;
if (type == OctreePrimitiveType::Triangle)
{
@@ -832,7 +832,7 @@ LooseOctree::computePrimitiveBoundingBox(OctreePrimitive* const pPrimitive, cons
lowerCorner = surfMesh->getVertexPosition(face[0]);
upperCorner = lowerCorner;
- const Vec3r verts12[2] = {
+ const Vec3d verts12[2] = {
surfMesh->getVertexPosition(face[1]),
surfMesh->getVertexPosition(face[2])
};
diff --git a/Source/DataStructures/imstkLooseOctree.h b/Source/DataStructures/imstkLooseOctree.h
index 885366bbdba21d6036453904445cf36336845fbd..b9803754f3e0fd99f74c838c8653ba40b7f2e09e 100644
--- a/Source/DataStructures/imstkLooseOctree.h
+++ b/Source/DataStructures/imstkLooseOctree.h
@@ -73,11 +73,11 @@ struct OctreePrimitive
union
{
- std::array<Real, 3> m_Position; ///> For a point primitive, store its position
+ std::array<double, 3> m_Position; ///> For a point primitive, store its position
struct
{
- std::array<Real, 3> m_LowerCorner; ///> For a non-point primitive, store its AABB's lower corner
- std::array<Real, 3> m_UpperCorner; ///> For a non-point primitive, store its AABB's upper corner
+ std::array<double, 3> m_LowerCorner; ///> For a non-point primitive, store its AABB's lower corner
+ std::array<double, 3> m_UpperCorner; ///> For a non-point primitive, store its AABB's upper corner
};
};
@@ -110,8 +110,8 @@ public:
OctreeNode() :
m_pTree(nullptr),
m_pParent(nullptr),
- m_Center(Vec3r(0, 0, 0)),
- m_HalfWidth(0),
+ m_Center(Vec3d::Zero()),
+ m_HalfWidth(0.0),
m_Depth(0),
m_MaxDepth(0),
m_bIsLeaf(true) {}
@@ -119,8 +119,8 @@ public:
///
/// \brief OctreeNode constructor, called during node splitting when initializing children node
///
- explicit OctreeNode(LooseOctree* const tree, OctreeNode* const pParent, const Vec3r& nodeCenter,
- const Real halfWidth, const uint32_t depth);
+ explicit OctreeNode(LooseOctree* const tree, OctreeNode* const pParent, const Vec3d& nodeCenter,
+ const double halfWidth, const uint32_t depth);
///
/// \brief Check if this node is a leaf node
@@ -204,9 +204,9 @@ public:
///
/// \brief Check if the given point is contained exactly in the node boundary (bounding box)
///
- bool contains(const Vec3r& point) { return contains(point[0], point[1], point[2]); }
- bool contains(const std::array<Real, 3>& point) { return contains(point[0], point[1], point[2]); }
- bool contains(const Real x, const Real y, const Real z)
+ bool contains(const Vec3d& point) { return contains(point[0], point[1], point[2]); }
+ bool contains(const std::array<double, 3>& point) { return contains(point[0], point[1], point[2]); }
+ bool contains(const double x, const double y, const double z)
{
return x >= m_LowerBound[0]
&& y >= m_LowerBound[1]
@@ -221,7 +221,7 @@ public:
/// \param lowerCorner The AABB's lower corner of the primitive
/// \param upperCorner The AABB's upper corner of the primitive
///
- bool contains(const std::array<Real, 3>& lowerCorner, const std::array<Real, 3>& upperCorner)
+ bool contains(const std::array<double, 3>& lowerCorner, const std::array<double, 3>& upperCorner)
{
return lowerCorner[0] >= m_LowerBound[0]
&& lowerCorner[1] >= m_LowerBound[1]
@@ -234,9 +234,9 @@ public:
///
/// \brief Check if the given point is contained in the node loose boundary (which is 2X bigger than the bounding box)
///
- bool looselyContains(const Vec3r& point) { return looselyContains(point[0], point[1], point[2]); }
- bool looselyContains(const std::array<Real, 3>& point) { return looselyContains(point[0], point[1], point[2]); }
- bool looselyContains(const Real x, const Real y, const Real z)
+ bool looselyContains(const Vec3d& point) { return looselyContains(point[0], point[1], point[2]); }
+ bool looselyContains(const std::array<double, 3>& point) { return looselyContains(point[0], point[1], point[2]); }
+ bool looselyContains(const double x, const double y, const double z)
{
return x >= m_LowerExtendedBound[0]
&& y >= m_LowerExtendedBound[1]
@@ -252,7 +252,7 @@ public:
/// \param lowerCorner The AABB's lower corner of the primitive
/// \param upperCorner The AABB's upper corner of the primitive
///
- bool looselyContains(const std::array<Real, 3>& lowerCorner, const std::array<Real, 3>& upperCorner)
+ bool looselyContains(const std::array<double, 3>& lowerCorner, const std::array<double, 3>& upperCorner)
{
return lowerCorner[0] >= m_LowerExtendedBound[0]
&& lowerCorner[1] >= m_LowerExtendedBound[1]
@@ -268,7 +268,7 @@ public:
/// \param lowerCorner The AABB's lower corner of the primitive
/// \param upperCorner The AABB's upper corner of the primitive
///
- bool looselyOverlaps(const std::array<Real, 3>& lowerCorner, const std::array<Real, 3>& upperCorner)
+ bool looselyOverlaps(const std::array<double, 3>& lowerCorner, const std::array<double, 3>& upperCorner)
{
return upperCorner[0] >= m_LowerExtendedBound[0]
&& upperCorner[1] >= m_LowerExtendedBound[1]
@@ -284,12 +284,12 @@ public:
OctreeNode* m_pParent; ///> Pointer to the parent node
OctreeNodeBlock* m_pChildren = nullptr; ///> Pointer to a memory block containing 8 children nodes
- const Vec3r m_Center; ///> Center of this node
- const Vec3r m_LowerBound; ///> The AABB's lower corner of the node
- const Vec3r m_UpperBound; ///> The AABB's upper corner of the node
- const Vec3r m_LowerExtendedBound; ///> The extended AABB's lower corner of the node, which is 2X bigger than the exact AABB
- const Vec3r m_UpperExtendedBound; ///> The extended AABB's upper corner of the node, which is 2X bigger than the exact AABB
- const Real m_HalfWidth; ///> Half width of the node AABB
+ const Vec3d m_Center; ///> Center of this node
+ const Vec3d m_LowerBound; ///> The AABB's lower corner of the node
+ const Vec3d m_UpperBound; ///> The AABB's upper corner of the node
+ const Vec3d m_LowerExtendedBound; ///> The extended AABB's lower corner of the node, which is 2X bigger than the exact AABB
+ const Vec3d m_UpperExtendedBound; ///> The extended AABB's upper corner of the node, which is 2X bigger than the exact AABB
+ const double m_HalfWidth; ///> Half width of the node AABB
const uint32_t m_Depth; ///> Depth of this node (depth > 0, depth = 1 starting at the root node)
uint32_t m_MaxDepth; ///> Cache the max depth of the tree (maximum depth level possible)
bool m_bIsLeaf = true; ///> True if this node does not have any child node (a node should have either 0 or 8 children)
@@ -338,8 +338,8 @@ public:
/// \param minWidthRatio If there is primitive that is not a point, minWidth will be recomputed as minWidth = min(width of all non-point primitives) * minWidthRatio
/// \param name Name of the octree
///
- explicit LooseOctree(const Vec3r& center, const Real width, const Real minWidth,
- const Real minWidthRatio = 1.0, const std::string name = "LooseOctree");
+ explicit LooseOctree(const Vec3d& center, const double width, const double minWidth,
+ const double minWidthRatio = 1.0, const std::string name = "LooseOctree");
///
/// Destructor, doing memory cleanup
@@ -359,17 +359,17 @@ public:
///
/// \brief Return center of the tree
///
- const Vec3r getCenter() const { return m_Center; }
+ const Vec3d getCenter() const { return m_Center; }
///
/// \brief Return width of the tree
///
- Real getWidth() const { return m_Width; }
+ double getWidth() const { return m_Width; }
///
/// \brief Return width of the lowest level tree nodes
///
- Real getMinWidth() const { return m_MinWidth; }
+ double getMinWidth() const { return m_MinWidth; }
///
/// \brief Get the maximum depth in this tree, which is computed based on the minWidth value
@@ -526,13 +526,13 @@ protected:
void deallocateMemoryPool();
const std::string m_Name; ///> Name of the tree
- const Vec3r m_Center; ///> Center of the tree
- const Real m_Width; ///> Width of the tree bounding box
+ const Vec3d m_Center; ///> Center of the tree
+ const double m_Width; ///> Width of the tree bounding box
/// If there is no point primitive, minWidth will be recomputed as minWidth = min(width of all non-point primitives) * minWidthRatio
- const Real m_MinWidthRatio;
+ const double m_MinWidthRatio;
- Real m_MinWidth; ///> Minimum width allowed for the tree nodes
+ double m_MinWidth; ///> Minimum width allowed for the tree nodes
uint32_t m_MaxDepth; ///> Max depth of the tree, which is computed based on m_MinWidth
bool m_useMaxDepth; ///> If on max depth specified by user will be used, otherwise maxdepth is based of minwidth
diff --git a/Source/DataStructures/imstkNeighborSearch.cpp b/Source/DataStructures/imstkNeighborSearch.cpp
index 55d3e1622b026d0220783f8f7c230135297f5367..a31f2243afa236ad70d097e86297b001041cd66b 100644
--- a/Source/DataStructures/imstkNeighborSearch.cpp
+++ b/Source/DataStructures/imstkNeighborSearch.cpp
@@ -27,7 +27,7 @@
namespace imstk
{
-NeighborSearch::NeighborSearch(NeighborSearch::Method searchMethod, Real searchRadius /*= 0*/) :
+NeighborSearch::NeighborSearch(NeighborSearch::Method searchMethod, double searchRadius /*= 0*/) :
m_Method(searchMethod), m_SearchRadius(searchRadius)
{
if (m_Method == Method::UniformGridBasedSearch)
@@ -43,7 +43,7 @@ NeighborSearch::NeighborSearch(NeighborSearch::Method searchMethod, Real searchR
}
void
-NeighborSearch::setSearchRadius(const Real searchRadius)
+NeighborSearch::setSearchRadius(const double searchRadius)
{
m_SearchRadius = searchRadius;
if (m_Method == Method::UniformGridBasedSearch)
diff --git a/Source/DataStructures/imstkNeighborSearch.h b/Source/DataStructures/imstkNeighborSearch.h
index 5077c0b7db6d78f4bb4a6797ab70ea33b62e15c4..1aa155218031270cae8898b69216a16155363ef4 100644
--- a/Source/DataStructures/imstkNeighborSearch.h
+++ b/Source/DataStructures/imstkNeighborSearch.h
@@ -46,17 +46,17 @@ public:
/// \brief Constructor
/// \param The selected search method
///
- explicit NeighborSearch(Method searchMethod, Real searchRadius = 0);
+ explicit NeighborSearch(Method searchMethod, double searchRadius = 0.0);
///
/// \brief Set the search radius
///
- void setSearchRadius(const Real searchRadius);
+ void setSearchRadius(const double searchRadius);
///
/// \brief Get the current search radius
///
- Real getSearchRadius() const { return m_SearchRadius; }
+ double getSearchRadius() const { return m_SearchRadius; }
///
/// \brief Search neighbors for each points within the search radius
@@ -82,7 +82,7 @@ public:
private:
Method m_Method;
- Real m_SearchRadius = 0;
+ double m_SearchRadius = 0.0;
std::shared_ptr<GridBasedNeighborSearch> m_GridBasedSearcher;
std::shared_ptr<SpatialHashTableSeparateChaining> m_SpatialHashSearcher;
diff --git a/Source/DataStructures/imstkUniformSpatialGrid.h b/Source/DataStructures/imstkUniformSpatialGrid.h
index 3ed8dbfa31861d56bf730e6f2178b7734934b25a..dc4bbe210e54b933446eeaea4c49610df281f6e3 100644
--- a/Source/DataStructures/imstkUniformSpatialGrid.h
+++ b/Source/DataStructures/imstkUniformSpatialGrid.h
@@ -40,7 +40,7 @@ public:
///
/// \brief Construct a default grid ([0, 1]^3) with cell size of 1
///
- UniformSpatialGrid() : UniformSpatialGrid(Vec3r(0, 0, 0), Vec3r(1, 1, 1), Real(1.0))
+ UniformSpatialGrid() : UniformSpatialGrid(Vec3d::Zero(), Vec3d(1.0, 1.0, 1.0), double(1.0))
{}
///
@@ -49,7 +49,7 @@ public:
/// \param upperCorner The upper corner of the grid
/// \param cellSize The length of grid cell
///
- UniformSpatialGrid(const Vec3r& lowerCorner, const Vec3r& upperCorner, Real cellSize)
+ UniformSpatialGrid(const Vec3d& lowerCorner, const Vec3d& upperCorner, double cellSize)
{
initialize(lowerCorner, upperCorner, cellSize);
}
@@ -60,7 +60,7 @@ public:
/// \param upperCorner The upper corner of the grid
/// \param cellSize the edge length of grid cell
///
- void initialize(const Vec3r& lowerCorner, const Vec3r& upperCorner, const Real cellSize)
+ void initialize(const Vec3d& lowerCorner, const Vec3d& upperCorner, const double cellSize)
{
CHECK(cellSize > 0) << "Invalid cell size";
@@ -68,7 +68,7 @@ public:
m_UpperCorner = upperCorner;
m_CellSize = cellSize;
- m_InvCellSize = Real(1.0) / m_CellSize;
+ m_InvCellSize = 1.0 / m_CellSize;
m_NTotalCells = 1u;
for (int i = 0; i < 3; ++i)
@@ -113,7 +113,7 @@ public:
/// \brief Get the 3D index (cell_x, cell_y, cell_z) of the cell containing the given positions
///
template<class IndexType>
- std::array<IndexType, 3> getCell3DIndices(const Vec3r& ppos) const
+ std::array<IndexType, 3> getCell3DIndices(const Vec3d& ppos) const
{
std::array<IndexType, 3> cellIdx;
for (int d = 0; d < 3; ++d)
@@ -137,7 +137,7 @@ public:
/// \brief Get data in a cell
/// \param A position in space
///
- CellData& getCellData(const Vec3r& ppos)
+ CellData& getCellData(const Vec3d& ppos)
{
return m_CellData[getCellLinearizedIndex<unsigned int>(ppos)];
}
@@ -146,7 +146,7 @@ public:
/// \brief Get data in a cell
/// \param A position in space
///
- const CellData& getCellData(const Vec3r& ppos) const
+ const CellData& getCellData(const Vec3d& ppos) const
{ return m_CellData[getCellLinearizedIndex<unsigned int>(ppos)]; }
///
@@ -238,7 +238,7 @@ public:
/// \param A position in space
///
template<class IndexType>
- IndexType getCellLinearizedIndex(const Vec3r& ppos) const
+ IndexType getCellLinearizedIndex(const Vec3d& ppos) const
{
auto cellIdx = getCell3DIndices<IndexType>(ppos);
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
@@ -252,10 +252,10 @@ public:
}
private:
- Vec3r m_LowerCorner; ///> Lower corner of the grid
- Vec3r m_UpperCorner; ///> Upper corner of the grid
- Real m_CellSize; ///> Length of grid cell
- Real m_InvCellSize; ///> Inverse length of grid cell
+ Vec3d m_LowerCorner; ///> Lower corner of the grid
+ Vec3d m_UpperCorner; ///> Upper corner of the grid
+ double m_CellSize; ///> Length of grid cell
+ double m_InvCellSize; ///> Inverse length of grid cell
unsigned int m_NTotalCells; ///> Number of total cells
std::array<unsigned int, 3> m_Resolution; ///> Grid resolution (number of cells in x/y/z dimensions)
diff --git a/Source/DynamicalModels/CMakeLists.txt b/Source/DynamicalModels/CMakeLists.txt
index 69025909b1a8b1d853cfba643de7eaca1ace2b0d..b2cf408807c9280b3f2e171612093172064fea6e 100644
--- a/Source/DynamicalModels/CMakeLists.txt
+++ b/Source/DynamicalModels/CMakeLists.txt
@@ -57,6 +57,6 @@ imstk_add_library(DynamicalModels
#-----------------------------------------------------------------------------
# Testing
#-----------------------------------------------------------------------------
-#if( ${PROJECT_NAME}_BUILD_TESTING )
-# add_subdirectory( Testing )
-#endif()
+if( ${PROJECT_NAME}_BUILD_TESTING )
+ add_subdirectory(Testing)
+endif()
diff --git a/Source/DynamicalModels/ObjectModels/imstkLevelSetModel.h b/Source/DynamicalModels/ObjectModels/imstkLevelSetModel.h
index cd0ee2eb99348a10580bd6dd7f56883762c3584b..358345f73d250535cd761699383e2b9c9adc53a0 100644
--- a/Source/DynamicalModels/ObjectModels/imstkLevelSetModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkLevelSetModel.h
@@ -65,7 +65,7 @@ public:
///
/// \brief Set the time step size
///
- virtual void setTimeStep(const Real timeStep) override { m_config->m_dt = timeStep; }
+ virtual void setTimeStep(const double timeStep) override { m_config->m_dt = timeStep; }
///
/// \brief Returns the time step size
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdConstraintFunctor.h b/Source/DynamicalModels/ObjectModels/imstkPbdConstraintFunctor.h
index 45a9135f7c92ca4aa819bc48cb86290f690f40cd..a4a188c904b4a88d5124c393a106439b564a640d 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdConstraintFunctor.h
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdConstraintFunctor.h
@@ -50,6 +50,7 @@ struct PbdConstraintFunctor
PbdConstraintFunctor() = default;
virtual ~PbdConstraintFunctor() = default;
+ public:
///
/// \brief Appends a set of constraint to the container given a geometry
///
@@ -58,7 +59,7 @@ struct PbdConstraintFunctor
void setGeometry(std::shared_ptr<PointSet> geom) { m_geom = geom; }
public:
- std::shared_ptr<PointSet> m_geom;
+ std::shared_ptr<PointSet> m_geom = nullptr;
};
struct PbdDistanceConstraintFunctor : public PbdConstraintFunctor
@@ -67,10 +68,12 @@ struct PbdDistanceConstraintFunctor : public PbdConstraintFunctor
PbdDistanceConstraintFunctor() = default;
~PbdDistanceConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
- const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
- auto addDistConstraint =
+ std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
+ const VecDataArray<double, 3>& vertices = *verticesPtr;
+ 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
@@ -134,10 +137,10 @@ struct PbdDistanceConstraintFunctor : public PbdConstraintFunctor
}
}
- void setStiffness(double stiffness) { m_stiffness = stiffness; }
+ void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
- double m_stiffness;
+ double m_stiffness = 0.0;
};
struct PbdFemConstraintFunctor : public PbdConstraintFunctor
@@ -146,6 +149,7 @@ struct PbdFemConstraintFunctor : public PbdConstraintFunctor
PbdFemConstraintFunctor() = default;
~PbdFemConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// Check if constraint type matches the mesh type
@@ -153,9 +157,11 @@ struct PbdFemConstraintFunctor : public PbdConstraintFunctor
<< "FEM Tetrahedral constraint should come with tetrahedral 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();
+ auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
+ std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
+ const VecDataArray<double, 3>& vertices = *verticesPtr;
+ std::shared_ptr<VecDataArray<int, 4>> elementsPtr = tetMesh->getTetrahedraIndices();
+ const VecDataArray<int, 4>& elements = *elementsPtr;
ParallelUtils::parallelFor(elements.size(),
[&](const size_t k)
@@ -168,12 +174,12 @@ struct PbdFemConstraintFunctor : public PbdConstraintFunctor
}, elements.size() > 100);
}
- void setMaterialType(PbdFEMTetConstraint::MaterialType materialType) { m_matType = materialType; }
+ void setMaterialType(const PbdFEMTetConstraint::MaterialType materialType) { m_matType = materialType; }
void setFemConfig(std::shared_ptr<PbdFEMConstraintConfig> femConfig) { m_femConfig = femConfig; }
protected:
- PbdFEMTetConstraint::MaterialType m_matType;
- std::shared_ptr<PbdFEMConstraintConfig> m_femConfig;
+ PbdFEMTetConstraint::MaterialType m_matType = PbdFEMTetConstraint::MaterialType::StVK;
+ std::shared_ptr<PbdFEMConstraintConfig> m_femConfig = nullptr;
};
struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
@@ -182,6 +188,7 @@ struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
PbdVolumeConstraintFunctor() = default;
~PbdVolumeConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// Check if constraint type matches the mesh type
@@ -189,9 +196,11 @@ struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
<< "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();
+ auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
+ std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
+ const VecDataArray<double, 3>& vertices = *verticesPtr;
+ std::shared_ptr<VecDataArray<int, 4>> elementsPtr = tetMesh->getTetrahedraIndices();
+ const VecDataArray<int, 4>& elements = *elementsPtr;
ParallelUtils::parallelFor(elements.size(),
[&](const size_t k)
@@ -204,10 +213,10 @@ struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
});
}
- void setStiffness(double stiffness) { m_stiffness = stiffness; }
+ void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
- double m_stiffness;
+ double m_stiffness = 0.0;
};
struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
@@ -216,6 +225,7 @@ struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
PbdAreaConstraintFunctor() = default;
~PbdAreaConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// check if constraint type matches the mesh type
@@ -223,9 +233,11 @@ struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
<< "Area constraint should come with a triangular mesh";
// ok, now create constraints
- auto triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
- const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
- const VecDataArray<int, 3>& elements = *triMesh->getTriangleIndices();
+ auto triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
+ std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
+ const VecDataArray<double, 3>& vertices = *verticesPtr;
+ std::shared_ptr<VecDataArray<int, 3>> elemenstPtr = triMesh->getTriangleIndices();
+ const VecDataArray<int, 3>& elements = *elemenstPtr;
ParallelUtils::parallelFor(elements.size(),
[&](const size_t k)
@@ -237,10 +249,10 @@ struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
});
}
- void setStiffness(double stiffness) { m_stiffness = stiffness; }
+ void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
- double m_stiffness;
+ double m_stiffness = 0.0;
};
struct PbdBendConstraintFunctor : public PbdConstraintFunctor
@@ -249,14 +261,17 @@ struct PbdBendConstraintFunctor : public PbdConstraintFunctor
PbdBendConstraintFunctor() = default;
~PbdBendConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
CHECK(m_geom->getTypeName() == "LineMesh")
<< "Bend constraint should come with a line mesh";
- auto lineMesh = std::dynamic_pointer_cast<LineMesh>(m_geom);
- const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
- const VecDataArray<int, 2>& elements = *lineMesh->getLinesIndices();
+ auto lineMesh = std::dynamic_pointer_cast<LineMesh>(m_geom);
+ std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
+ const VecDataArray<double, 3>& vertices = *verticesPtr;
+ /*std::shared_ptr< VecDataArray<int, 2>> indicesPtr = lineMesh->getLinesIndices();
+ const VecDataArray<int, 2>& indices = *indicesPtr*/
auto addBendConstraint =
[&](const double k, size_t i1, size_t i2, size_t i3)
@@ -277,24 +292,23 @@ struct PbdBendConstraintFunctor : public PbdConstraintFunctor
constraints.addConstraint(c);
};
- // Iterate sets of two segments
- for (int k = 0; k < elements.size() - 1; k++)
+ // Iterate sets of stride # of segments
+ for (int k = 0; k < vertices.size() - m_stride * 2; k += m_stride)
{
- auto& seg1 = elements[k];
- auto& seg2 = elements[k + 1];
- int i3 = seg2[0];
- if (i3 == seg1[0] || i3 == seg1[1])
- {
- i3 = seg2[1];
- }
- addBendConstraint(m_stiffness, seg1[0], seg1[1], i3);
+ addBendConstraint(m_stiffness, k, k + m_stride, k + 2 * m_stride);
}
}
- void setStiffness(double stiffness) { m_stiffness = stiffness; }
+ void setStiffness(const double stiffness) { m_stiffness = stiffness; }
+ void setStride(const int stride)
+ {
+ CHECK(m_stride > 1) << "Stride should be at least 1.";
+ m_stride = stride;
+ }
protected:
- double m_stiffness;
+ double m_stiffness = 0.0;
+ int m_stride = 3;
};
struct PbdDihedralConstraintFunctor : public PbdConstraintFunctor
@@ -303,17 +317,20 @@ struct PbdDihedralConstraintFunctor : public PbdConstraintFunctor
PbdDihedralConstraintFunctor() = default;
~PbdDihedralConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
CHECK(m_geom->getTypeName() == "SurfaceMesh")
<< "Dihedral constraint should come with a triangular mesh";
// Create constraints
- auto triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
- const VecDataArray<double, 3>& vertices = *triMesh->getVertexPositions();
- const VecDataArray<int, 3>& elements = *triMesh->getTriangleIndices();
- const int nV = triMesh->getNumVertices();
- std::vector<std::vector<int>> vertIdsToTriangleIds(nV);
+ auto triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
+ std::shared_ptr<VecDataArray<double, 3>> verticesPtr = triMesh->getVertexPositions();
+ const VecDataArray<double, 3>& vertices = *verticesPtr;
+ std::shared_ptr<VecDataArray<int, 3>> elementsPtr = triMesh->getTriangleIndices();
+ const VecDataArray<int, 3>& elements = *elementsPtr;
+ const int nV = triMesh->getNumVertices();
+ std::vector<std::vector<int>> vertIdsToTriangleIds(nV);
for (int k = 0; k < elements.size(); ++k)
{
@@ -389,10 +406,10 @@ struct PbdDihedralConstraintFunctor : public PbdConstraintFunctor
}
}
- void setStiffness(double stiffness) { m_stiffness = stiffness; }
+ void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
- double m_stiffness;
+ double m_stiffness = 0.0;
};
struct PbdConstantDensityConstraintFunctor : public PbdConstraintFunctor
@@ -401,6 +418,7 @@ struct PbdConstantDensityConstraintFunctor : public PbdConstraintFunctor
PbdConstantDensityConstraintFunctor() = default;
~PbdConstantDensityConstraintFunctor() override = default;
+ public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// check if constraint type matches the mesh type
@@ -412,9 +430,9 @@ struct PbdConstantDensityConstraintFunctor : public PbdConstraintFunctor
constraints.addConstraint(c);
}
- void setStiffness(double stiffness) { m_stiffness = stiffness; }
+ void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
- double m_stiffness;
+ double m_stiffness = 0.0;
};
}
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index 33133b82a3de4ab1c03c0a82bfac5c279b8481f5..7858be715a5cfc687e597fdc66ef310e475945ff 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -64,6 +64,13 @@ PBDModelConfig::enableConstraint(PbdConstraint::Type type, double stiffness)
m_regularConstraints.push_back({ type, stiffness });
}
+void
+PBDModelConfig::enableBendConstraint(const double stiffness, const int stride)
+{
+ m_regularConstraints.push_back({ PbdConstraint::Type::Bend, stiffness });
+ m_constraintStrides[m_regularConstraints.size() - 1] = stride;
+}
+
void
PBDModelConfig::enableFEMConstraint(PbdConstraint::Type type, PbdFEMConstraint::MaterialType material)
{
@@ -115,6 +122,7 @@ PbdModel::initialize()
functor->setMaterialType(constraintType.second);
m_functors.push_back(functor);
}
+ int i = 0;
for (auto constraintType : m_parameters->m_regularConstraints)
{
if (constraintType.first == PbdConstraint::Type::Area)
@@ -125,8 +133,15 @@ PbdModel::initialize()
}
else if (constraintType.first == PbdConstraint::Type::Bend)
{
+ // If bend constraint requested but no stride provided use 1
+ int stride = 1;
+ if (m_parameters->m_constraintStrides.count(i) > 0)
+ {
+ stride = m_parameters->m_constraintStrides[i];
+ }
auto functor = std::make_shared<PbdBendConstraintFunctor>();
functor->setStiffness(constraintType.second);
+ functor->setStride(stride);
m_functors.push_back(functor);
}
else if (constraintType.first == PbdConstraint::Type::ConstantDensity)
@@ -165,6 +180,7 @@ PbdModel::initialize()
functor->setStiffness(constraintType.second);
m_functors.push_back(functor);
}
+ i++;
}
for (auto functorPtr : m_functors)
@@ -177,7 +193,7 @@ PbdModel::initialize()
// Partition constraints for parallel computation
if (m_parameters->m_doPartitioning)
{
- m_constraints->partitionConstraints(m_partitionThreshold);
+ m_constraints->partitionConstraints(static_cast<int>(m_partitionThreshold));
}
else
{
@@ -286,20 +302,20 @@ PbdModel::initGraphEdges(std::shared_ptr<TaskNode> source, std::shared_ptr<TaskN
void
PbdModel::computeElasticConstants()
{
- if (std::abs(m_parameters->m_femParams->m_mu) < MIN_REAL
- && std::abs(m_parameters->m_femParams->m_lambda) < MIN_REAL)
+ if (std::abs(m_parameters->m_femParams->m_mu) < std::numeric_limits<double>::min()
+ && std::abs(m_parameters->m_femParams->m_lambda) < std::numeric_limits<double>::min())
{
- const auto E = m_parameters->m_femParams->m_YoungModulus;
- const auto nu = m_parameters->m_femParams->m_PoissonRatio;
- m_parameters->m_femParams->m_mu = E / Real(2.0) / (Real(1.0) + nu);
- m_parameters->m_femParams->m_lambda = E * nu / ((Real(1.0) + nu) * (Real(1.0) - Real(2.0) * nu));
+ const double E = m_parameters->m_femParams->m_YoungModulus;
+ const double nu = m_parameters->m_femParams->m_PoissonRatio;
+ m_parameters->m_femParams->m_mu = E / 2.0 / (1.0 + nu);
+ m_parameters->m_femParams->m_lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
}
else
{
- const auto mu = m_parameters->m_femParams->m_mu;
- const auto lambda = m_parameters->m_femParams->m_lambda;
- m_parameters->m_femParams->m_YoungModulus = mu * (Real(3.0) * lambda + Real(2.0) * mu) / (lambda + mu);
- m_parameters->m_femParams->m_PoissonRatio = lambda / Real(2.0) / (lambda + mu);
+ const double mu = m_parameters->m_femParams->m_mu;
+ const double lambda = m_parameters->m_femParams->m_lambda;
+ m_parameters->m_femParams->m_YoungModulus = mu * (3.0 * lambda + 2.0 * mu) / (lambda + mu);
+ m_parameters->m_femParams->m_PoissonRatio = lambda / 2.0 / (lambda + mu);
}
}
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
index d47b8318f58e8d0bcfe5a176b20eb6132505d295..3a07f09d60f14ac23e9bcb825dfb6a68b33e7c65 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
@@ -62,6 +62,7 @@ struct PBDModelConfig
}); ///> Info shared between the fem constraints
std::vector<std::pair<PbdConstraint::Type, double>> m_regularConstraints; ///> Constraints except FEM
+ std::unordered_map<int, int> m_constraintStrides; ///> Extra parameter for bend constraints
std::vector<std::pair<PbdConstraint::Type, PbdFEMConstraint::MaterialType>> m_FEMConstraints; ///> Constraints except FEM
///
@@ -70,6 +71,12 @@ struct PBDModelConfig
///
void enableConstraint(PbdConstraint::Type type, double stiffness);
+ ///
+ /// \brief Enables a bend constraint with given stride (stride=1 connects every neighbor vertex
+ /// so a constraint is formed every 3 verts, stride=2 forms between 5 verts)
+ ///
+ void enableBendConstraint(const double stiffness, const int stride);
+
///
/// \brief Enable a FEM constraint with mu, lambda
///
diff --git a/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.cpp b/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.cpp
index cafaa14fc771f36109153182db0bda5563e9d6a3..2a8757ee053a920d0bbb1d91e59980a4ad75fdbd 100644
--- a/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.cpp
@@ -73,17 +73,17 @@ RigidBodyModel2::initialize()
// Compute the initial state
std::shared_ptr<RigidBodyState2> state = std::make_shared<RigidBodyState2>();
state->resize(m_bodies.size());
- std::vector<bool>& isStatic = state->getIsStatic();
- StdVectorOfReal& invMasses = state->getInvMasses();
- StdVectorOfMat3d& invInteriaTensors = state->getInvIntertiaTensors();
- StdVectorOfVec3d& positions = state->getPositions();
- StdVectorOfQuatd& orientations = state->getOrientations();
- StdVectorOfVec3d& velocities = state->getVelocities();
- StdVectorOfVec3d& angularVelocities = state->getAngularVelocities();
- StdVectorOfVec3d& tentativeVelocities = state->getTentatveVelocities();
- StdVectorOfVec3d& tentativeAngularVelocities = state->getTentativeAngularVelocities();
- StdVectorOfVec3d& forces = state->getForces();
- StdVectorOfVec3d& torques = state->getTorques();
+ std::vector<bool>& isStatic = state->getIsStatic();
+ std::vector<double>& invMasses = state->getInvMasses();
+ StdVectorOfMat3d& invInteriaTensors = state->getInvIntertiaTensors();
+ StdVectorOfVec3d& positions = state->getPositions();
+ StdVectorOfQuatd& orientations = state->getOrientations();
+ StdVectorOfVec3d& velocities = state->getVelocities();
+ StdVectorOfVec3d& angularVelocities = state->getAngularVelocities();
+ StdVectorOfVec3d& tentativeVelocities = state->getTentatveVelocities();
+ StdVectorOfVec3d& tentativeAngularVelocities = state->getTentativeAngularVelocities();
+ StdVectorOfVec3d& forces = state->getForces();
+ StdVectorOfVec3d& torques = state->getTorques();
m_Minv = Eigen::SparseMatrix<double>(m_bodies.size() * 6, m_bodies.size() * 6);
std::vector<Eigen::Triplet<double>> mInvTriplets;
@@ -151,6 +151,56 @@ RigidBodyModel2::initialize()
return true;
}
+void
+RigidBodyModel2::updateMass()
+{
+ // Compute the initial state
+ std::shared_ptr<RigidBodyState2> state = m_currentState;
+
+ std::vector<double>& invMasses = state->getInvMasses();
+ StdVectorOfMat3d& invInteriaTensors = state->getInvIntertiaTensors();
+
+ m_Minv = Eigen::SparseMatrix<double>(m_bodies.size() * 6, m_bodies.size() * 6);
+ std::vector<Eigen::Triplet<double>> mInvTriplets;
+ mInvTriplets.reserve((9 + 3) * m_bodies.size());
+ for (size_t i = 0; i < m_bodies.size(); i++)
+ {
+ RigidBody& body = *m_bodies[i];
+
+ invMasses[i] = (body.m_mass == 0.0) ? 0.0 : 1.0 / body.m_mass;
+ if (body.m_intertiaTensor.determinant() == 0.0)
+ {
+ LOG(WARNING) << "Inertia tensor provided is not invertible, check that it makes sense";
+ }
+ else
+ {
+ invInteriaTensors[i] = body.m_intertiaTensor.inverse();
+ }
+
+ if (!body.m_isStatic)
+ {
+ // invMass expanded to 3x3 matrix
+ const double invMass = invMasses[i];
+ const Mat3d& invInvertia = invInteriaTensors[i];
+ int index = static_cast<int>(i * 6);
+ mInvTriplets.push_back(Eigen::Triplet<double>(index, index, invMass));
+ index++;
+ mInvTriplets.push_back(Eigen::Triplet<double>(index, index, invMass));
+ index++;
+ mInvTriplets.push_back(Eigen::Triplet<double>(index, index, invMass));
+ index++;
+ for (unsigned int c = 0; c < 3; c++)
+ {
+ for (unsigned int r = 0; r < 3; r++)
+ {
+ mInvTriplets.push_back(Eigen::Triplet<double>(index + r, index + c, invInvertia(c, r)));
+ }
+ }
+ }
+ }
+ m_Minv.setFromTriplets(mInvTriplets.begin(), mInvTriplets.end());
+}
+
void
RigidBodyModel2::configure(std::shared_ptr<RigidBodyModel2Config> config)
{
@@ -160,14 +210,14 @@ RigidBodyModel2::configure(std::shared_ptr<RigidBodyModel2Config> config)
void
RigidBodyModel2::computeTentativeVelocities()
{
- const double dt = m_config->m_dt;
- const StdVectorOfReal& invMasses = getCurrentState()->getInvMasses();
- const StdVectorOfMat3d& invInteriaTensors = getCurrentState()->getInvIntertiaTensors();
- StdVectorOfVec3d& tentativeVelocities = getCurrentState()->getTentatveVelocities();
- StdVectorOfVec3d& tentativeAngularVelocities = getCurrentState()->getTentativeAngularVelocities();
- StdVectorOfVec3d& forces = getCurrentState()->getForces();
- StdVectorOfVec3d& torques = getCurrentState()->getTorques();
- const Vec3d& fG = m_config->m_gravity;
+ const double dt = m_config->m_dt;
+ const std::vector<double>& invMasses = getCurrentState()->getInvMasses();
+ const StdVectorOfMat3d& invInteriaTensors = getCurrentState()->getInvIntertiaTensors();
+ StdVectorOfVec3d& tentativeVelocities = getCurrentState()->getTentatveVelocities();
+ StdVectorOfVec3d& tentativeAngularVelocities = getCurrentState()->getTentativeAngularVelocities();
+ StdVectorOfVec3d& forces = getCurrentState()->getForces();
+ StdVectorOfVec3d& torques = getCurrentState()->getTorques();
+ const Vec3d& fG = m_config->m_gravity;
// Sum gravity to the forces
ParallelUtils::parallelFor(forces.size(), [&forces, &fG](const size_t& i)
@@ -343,8 +393,8 @@ RigidBodyModel2::integrate()
const std::vector<bool>& isStatic = getCurrentState()->getIsStatic();
- const StdVectorOfReal& invMasses = getCurrentState()->getInvMasses();
- const StdVectorOfMat3d& invInteriaTensors = getCurrentState()->getInvIntertiaTensors();
+ const std::vector<double>& invMasses = getCurrentState()->getInvMasses();
+ const StdVectorOfMat3d& invInteriaTensors = getCurrentState()->getInvIntertiaTensors();
StdVectorOfVec3d& positions = getCurrentState()->getPositions();
StdVectorOfQuatd& orientations = getCurrentState()->getOrientations();
diff --git a/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.h b/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.h
index 359dc84c4b96511b38e3732848a46a31a886b76b..13cd35df851dd520d7cc8b60c81911a25d34661b 100644
--- a/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.h
+++ b/Source/DynamicalModels/ObjectModels/imstkRigidBodyModel2.h
@@ -71,7 +71,7 @@ public:
///
/// \brief Set the time step size
///
- virtual void setTimeStep(const Real timeStep) override { m_config->m_dt = timeStep; }
+ virtual void setTimeStep(const double timeStep) override { m_config->m_dt = timeStep; }
///
/// \brief Returns the time step size
@@ -99,10 +99,17 @@ public:
void removeRigidBody(std::shared_ptr<RigidBody> body);
///
- /// \brief Initialize the RigidBody model
+ /// \brief Initialize the RigidBody model to the initial state
///
bool initialize() override;
+ ///
+ /// \brief Updates mass and inertia matrices to those provided
+ /// by the bodies. Not often needed unless mass/inertia is changing
+ /// at runtime.
+ ///
+ void updateMass();
+
///
/// \brief Configure the model
///
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.cpp b/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.cpp
index bd65abaf4d336ef34043f14d4c142f6b7693f40f..e8ce8be97c6ee349e3cc0ba3c76fdd58c1d748d4 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.cpp
@@ -26,11 +26,11 @@ limitations under the License.
namespace imstk
{
SPHBoundaryConditions::SPHBoundaryConditions(std::pair<Vec3d, Vec3d>& inletCoords, std::vector<std::pair<Vec3d, Vec3d>>& outletCoords, std::pair<Vec3d, Vec3d>& fluidCoords,
- const Vec3d& inletNormal, const StdVectorOfVec3d&, const Real inletRadius, const Vec3d& inletCenterPt, const double inletFlowRate,
+ const Vec3d& inletNormal, const StdVectorOfVec3d&, const double inletRadius, const Vec3d& inletCenterPt, const double inletFlowRate,
StdVectorOfVec3d& mainParticlePositions, const StdVectorOfVec3d& wallParticlePositions) :
m_inletDomain(inletCoords), m_outletDomain(outletCoords),
m_fluidDomain(fluidCoords),
- m_bufferCoord(Vec3d(100, 0, 0)),
+ m_bufferCoord(Vec3d(100.0, 0.0, 0.0)),
m_inletCenterPoint(inletCenterPt),
m_inletRadius(inletRadius),
m_inletNormal(inletNormal.normalized()),
@@ -114,20 +114,20 @@ SPHBoundaryConditions::setParticleTypes(const StdVectorOfVec3d& mainParticlePosi
std::iota(std::begin(m_bufferIndices), std::end(m_bufferIndices), m_particleTypes.size() - m_numBufferParticles);
}
-Vec3r
+Vec3d
SPHBoundaryConditions::computeParabolicInletVelocity(const Vec3d& particlePosition)
{
// compute distance of point
- const Vec3d inletRegionCenterPoint = (Vec3d(1, 1, 1) + m_inletNormal).array() * m_inletCenterPoint.array() + particlePosition.dot(m_inletNormal) * m_inletNormal.array();
+ const Vec3d inletRegionCenterPoint = (Vec3d(1.0, 1.0, 1.0) + m_inletNormal).array() * m_inletCenterPoint.array() + particlePosition.dot(m_inletNormal) * m_inletNormal.array();
const double distance = (particlePosition - inletRegionCenterPoint).norm();
Vec3d inletParabolicVelocity;
if (distance > m_inletRadius)
{
- inletParabolicVelocity = Vec3d(0, 0, 0);
+ inletParabolicVelocity = Vec3d::Zero();
}
else
{
- inletParabolicVelocity = m_inletVelocity * (1 - (distance / m_inletRadius) * (distance / m_inletRadius));
+ inletParabolicVelocity = m_inletVelocity * (1.0 - (distance / m_inletRadius) * (distance / m_inletRadius));
}
return inletParabolicVelocity;
}
@@ -136,19 +136,19 @@ void
SPHBoundaryConditions::addBoundaryParticles(StdVectorOfVec3d& mainParticlePositions, const StdVectorOfVec3d& wallParticlePositions)
{
mainParticlePositions.insert(mainParticlePositions.end(), wallParticlePositions.begin(), wallParticlePositions.end());
- mainParticlePositions.insert(mainParticlePositions.end(), m_numBufferParticles, Vec3d(100, 0, 0));
+ mainParticlePositions.insert(mainParticlePositions.end(), m_numBufferParticles, Vec3d(100.0, 0.0, 0.0));
}
void
-SPHBoundaryConditions::setInletVelocity(const Real flowRate)
+SPHBoundaryConditions::setInletVelocity(const double flowRate)
{
- m_inletVelocity = -m_inletNormal * (flowRate / m_inletCrossSectionalArea * 2);
+ m_inletVelocity = -m_inletNormal * (flowRate / m_inletCrossSectionalArea * 2.0);
}
Vec3d
SPHBoundaryConditions::placeParticleAtInlet(const Vec3d& position)
{
- const Vec3d inletPosition = (Vec3d(1, 1, 1) + m_inletNormal).cwiseProduct(position) - m_inletCenterPoint.cwiseProduct(m_inletNormal);
+ const Vec3d inletPosition = (Vec3d(1.0, 1.0, 1.0) + m_inletNormal).cwiseProduct(position) - m_inletCenterPoint.cwiseProduct(m_inletNormal);
return inletPosition;
}
} // end namespace imstk
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.h b/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.h
index 70698571e381f4dec8a0413e1e6276d0d22306ec..a7501738617fde13704fc500ee1d2842fdbb655d 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.h
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHBoundaryConditions.h
@@ -43,7 +43,7 @@ public:
public:
SPHBoundaryConditions(std::pair<Vec3d, Vec3d>& inletCoords, std::vector<std::pair<Vec3d, Vec3d>>& outletCoords, std::pair<Vec3d, Vec3d>& fluidCoords,
- const Vec3d& inletNormal, const StdVectorOfVec3d& outletNormals, const Real inletRadius, const Vec3d& inletCenterPt, const double inletFlowRate,
+ const Vec3d& inletNormal, const StdVectorOfVec3d& outletNormals, const double inletRadius, const Vec3d& inletCenterPt, const double inletFlowRate,
StdVectorOfVec3d& mainParticlePositions,
const StdVectorOfVec3d& wallParticlePositions);
@@ -60,11 +60,11 @@ public:
Vec3d getBufferCoord() { return m_bufferCoord; }
- Vec3r computeParabolicInletVelocity(const Vec3d& position);
+ Vec3d computeParabolicInletVelocity(const Vec3d& position);
void addBoundaryParticles(StdVectorOfVec3d& mainParticlePositions, const StdVectorOfVec3d& wallParticlePositions);
- void setInletVelocity(const Real flowRate);
+ void setInletVelocity(const double flowRate);
Vec3d getInletCoord() { return m_inletDomain.first; }
@@ -83,13 +83,13 @@ private:
std::vector<ParticleType> m_particleTypes;
- Vec3d m_bufferCoord;
- Vec3d m_inletCenterPoint;
- Real m_inletRadius;
- Vec3r m_inletVelocity;
- Vec3d m_inletNormal;
+ Vec3d m_bufferCoord;
+ Vec3d m_inletCenterPoint;
+ double m_inletRadius;
+ Vec3d m_inletVelocity;
+ Vec3d m_inletNormal;
- Real m_inletCrossSectionalArea;
+ double m_inletCrossSectionalArea;
const size_t m_numBufferParticles = 10000;
std::vector<size_t> m_bufferIndices;
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h b/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h
index 188f58facf1ae607e7f96d9045cbc1838cd9cc13..b991c095272ba6f8cc6ccb602e84e3643be2c9de 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHKernels.h
@@ -35,7 +35,7 @@ namespace SPH
template<int N>
class Poly6Kernel
{
-using VecXr = Eigen::Matrix<Real, N, 1>;
+using VecXd = Eigen::Matrix<double, N, 1>;
public:
Poly6Kernel()
@@ -46,7 +46,7 @@ public:
///
/// \brief Set the kernel radius
///
- void setRadius(const Real radius)
+ void setRadius(const double radius)
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
@@ -60,56 +60,56 @@ public:
#pragma warning(pop)
#endif
{
- m_k = Real(4.0) / (PI * std::pow(m_radius, 8));
- m_l = -Real(24.0) / (PI * std::pow(m_radius, 8));
+ m_k = 4.0 / (PI * std::pow(m_radius, 8));
+ m_l = -24.0 / (PI * std::pow(m_radius, 8));
}
else
{
- m_k = Real(315.0) / (Real(64.0) * PI * std::pow(m_radius, 9));
- m_l = -Real(945.0) / (Real(32.0) * PI * std::pow(m_radius, 9));
+ m_k = 315.0 / (64.0 * PI * std::pow(m_radius, 9));
+ m_l = -945.0 / (32.0 * PI * std::pow(m_radius, 9));
}
m_m = m_l;
- m_W0 = W(VecXr::Zero());
+ m_W0 = W(VecXd::Zero());
}
///
/// \brief Compute weight value
/// W(r,h) = (315/(64 PI h^9))(h^2-|r|^2)^3
///
- Real W(const Real r) const
+ double W(const double r) const
{
- const auto r2 = r * r;
+ const double r2 = r * r;
const double rd = m_radius2 - r2;
- return (r2 <= m_radius2) ? rd * rd * rd * m_k : Real(0);
+ return (r2 <= m_radius2) ? rd * rd * rd * m_k : 0.0;
}
///
/// \brief Compute weight value
/// W(r,h) = (315/(64 PI h^9))(h^2-|r|^2)^3
///
- Real W(const VecXr& r) const
+ double W(const VecXd& r) const
{
- const auto r2 = r.squaredNorm();
+ const double r2 = r.squaredNorm();
const double rd = m_radius2 - r2;
- return (r2 <= m_radius2) ? rd * rd * rd * m_k : Real(0);
+ return (r2 <= m_radius2) ? rd * rd * rd * m_k : 0.0;
}
///
/// \brief Get W(0)
///
- Real W0() const { return m_W0; }
+ double W0() const { return m_W0; }
///
/// \brief Compute weight gradient
/// grad(W(r,h)) = r(-945/(32 PI h^9))(h^2-|r|^2)^2
///
- VecXr gradW(const VecXr& r) const
+ VecXd gradW(const VecXd& r) const
{
- VecXr res = VecXr::Zero();
- const auto r2 = r.squaredNorm();
- if (r2 <= m_radius2 && r2 > Real(1e-12))
+ VecXd res = VecXd::Zero();
+ const double r2 = r.squaredNorm();
+ if (r2 <= m_radius2 && r2 > 1.0e-12)
{
- Real tmp = m_radius2 - r2;
+ double tmp = m_radius2 - r2;
res = m_l * tmp * tmp * r;
}
@@ -120,14 +120,14 @@ public:
/// \brief Compute laplacian
/// laplacian(W(r,h)) = (-945/(32 PI h^9))(h^2-|r|^2)(-7|r|^2+3h^2)
///
- Real laplacian(const VecXr& r) const
+ double laplacian(const VecXd& r) const
{
- Real res = 0.;
- const auto r2 = r.squaredNorm();
+ double res = 0.;
+ const double r2 = r.squaredNorm();
if (r2 <= m_radius2)
{
- Real tmp = m_radius2 - r2;
- Real tmp2 = Real(3.0) * m_radius2 - Real(7.0) * r2;
+ double tmp = m_radius2 - r2;
+ double tmp2 = 3.0 * m_radius2 - 7.0 * r2;
res = m_m * tmp * tmp2;
}
@@ -135,12 +135,12 @@ public:
}
protected:
- Real m_radius; ///> Kernel radius
- Real m_radius2; ///> Kernel radius squared
- Real m_k; ///> Kernel coefficient for W()
- Real m_l; ///> Kernel coefficient for gradW()
- Real m_m; ///> Kernel coefficient for laplacian()
- Real m_W0; ///> Precomputed W(0)
+ double m_radius; ///> Kernel radius
+ double m_radius2; ///> Kernel radius squared
+ double m_k; ///> Kernel coefficient for W()
+ double m_l; ///> Kernel coefficient for gradW()
+ double m_m; ///> Kernel coefficient for laplacian()
+ double m_W0; ///> Precomputed W(0)
};
///
@@ -150,7 +150,7 @@ protected:
template<int N>
class SpikyKernel
{
-using VecXr = Eigen::Matrix<Real, N, 1>;
+using VecXd = Eigen::Matrix<double, N, 1>;
public:
SpikyKernel()
@@ -161,7 +161,7 @@ public:
///
/// \brief Set the kernel radius
///
- void setRadius(const Real radius)
+ void setRadius(const double radius)
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
@@ -175,58 +175,58 @@ public:
#pragma warning(pop)
#endif
{
- const auto radius5 = std::pow(m_radius, 5);
- m_k = Real(10.0) / (PI * radius5);
- m_l = -Real(30.0) / (PI * radius5);
+ const double radius5 = std::pow(m_radius, 5);
+ m_k = 10.0 / (PI * radius5);
+ m_l = -30.0 / (PI * radius5);
}
else
{
- const auto radius6 = std::pow(m_radius, 6);
- m_k = Real(15.0) / (PI * radius6);
- m_l = -Real(45.0) / (PI * radius6);
+ const double radius6 = std::pow(m_radius, 6);
+ m_k = 15.0 / (PI * radius6);
+ m_l = -45.0 / (PI * radius6);
}
- m_W0 = W(VecXr::Zero());
+ m_W0 = W(VecXd::Zero());
}
///
/// \brief Compute weight value
/// W(r,h) = 15/(PI*h^6) * (h-r)^3
///
- Real W(const Real r) const
+ double W(const double r) const
{
const double rd = m_radius - r;
- return (r <= m_radius) ? rd * rd * rd * m_k : Real(0);
+ return (r <= m_radius) ? rd * rd * rd * m_k : 0.0;
}
///
/// \brief Compute weight value
/// W(r,h) = 15/(PI*h^6) * (h-r)^3
///
- Real W(const VecXr& r) const
+ double W(const VecXd& r) const
{
const double r2 = r.squaredNorm();
const double rd = m_radius - std::sqrt(r2);
- return (r2 <= m_radius2) ? rd * rd * rd * m_k : Real(0);
+ return (r2 <= m_radius2) ? rd * rd * rd * m_k : 0.0;
}
///
/// \brief Get W(0)
///
- Real W0() const { return m_W0; }
+ double W0() const { return m_W0; }
///
/// \brief Compute weight gradient
/// grad(W(r,h)) = -r(45/(PI*h^6) * (h-r)^2)
///
- VecXr gradW(const VecXr& r) const
+ VecXd gradW(const VecXd& r) const
{
- VecXr res = VecXr::Zero();
+ VecXd res = VecXd::Zero();
const auto r2 = r.squaredNorm();
- if (r2 <= m_radius2 && r2 > Real(1e-12))
+ if (r2 <= m_radius2 && r2 > 1.0e-12)
{
- const auto rl = std::sqrt(r2);
- const auto hr = m_radius - rl;
- const auto hr2 = hr * hr;
+ const double rl = std::sqrt(r2);
+ const double hr = m_radius - rl;
+ const double hr2 = hr * hr;
res = m_l * hr2 * (r / rl);
}
@@ -234,11 +234,11 @@ public:
}
protected:
- Real m_radius; ///> Kernel radius
- Real m_radius2; ///> Kernel radius squared
- Real m_k; ///> Kernel coefficient for W()
- Real m_l; ///> Kernel coefficient for gradW()
- Real m_W0; ///> Precomputed W(0)
+ double m_radius; ///> Kernel radius
+ double m_radius2; ///> Kernel radius squared
+ double m_k; ///> Kernel coefficient for W()
+ double m_l; ///> Kernel coefficient for gradW()
+ double m_W0; ///> Precomputed W(0)
};
///
@@ -248,7 +248,7 @@ protected:
template<int N>
class CohesionKernel
{
-using VecXr = Eigen::Matrix<Real, N, 1>;
+using VecXd = Eigen::Matrix<double, N, 1>;
public:
CohesionKernel()
@@ -259,7 +259,7 @@ public:
///
/// \brief Set the kernel radius
///
- void setRadius(const Real radius)
+ void setRadius(const double radius)
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
@@ -273,24 +273,24 @@ public:
#pragma warning(pop)
#endif
- m_k = Real(32.0) / (PI * std::pow(m_radius, 9));
- m_c = std::pow(m_radius, 6) / Real(64.0);
- m_W0 = W(VecXr::Zero());
+ m_k = 32.0 / (PI * std::pow(m_radius, 9));
+ m_c = std::pow(m_radius, 6) / 64.0;
+ m_W0 = W(VecXd::Zero());
}
///
/// \brief Compute weight value
/// W(r,h) = (32/(PI h^9))(h-r)^3*r^3 if h/2 < r <= h,
/// (32/(PI h^9))(2*(h-r)^3*r^3 - h^6/64 if 0 < r <= h/2
- Real W(const Real r) const
+ double W(const double r) const
{
- Real res = 0.;
- const auto r2 = r * r;
+ double res = 0.;
+ const double r2 = r * r;
if (r2 <= m_radius2)
{
- const auto r1 = std::sqrt(r2);
- const auto r3 = r2 * r1;
- if (r1 > Real(0.5) * m_radius)
+ const double r1 = std::sqrt(r2);
+ const double r3 = r2 * r1;
+ if (r1 > 0.5 * m_radius)
{
const double rd = m_radius - r1;
res = m_k * rd * rd * rd * r3;
@@ -298,7 +298,7 @@ public:
else
{
const double rd = m_radius - r1;
- res = m_k * Real(2.0) * rd * rd * rd * r3 - m_c;
+ res = m_k * 2.0 * rd * rd * rd * r3 - m_c;
}
}
return res;
@@ -308,15 +308,15 @@ public:
/// \brief Compute weight value
/// W(r,h) = (32/(PI h^9))(h-r)^3*r^3 if h/2 < r <= h,
/// (32/(PI h^9))(2*(h-r)^3*r^3 - h^6/64 if 0 < r <= h/2
- Real W(const VecXr& r) const
+ double W(const VecXd& r) const
{
- Real res = 0.;
- const auto r2 = r.squaredNorm();
+ double res = 0.;
+ const double r2 = r.squaredNorm();
if (r2 <= m_radius2)
{
- const auto r1 = std::sqrt(r2);
- const auto r3 = r2 * r1;
- if (r1 > Real(0.5) * m_radius)
+ const double r1 = std::sqrt(r2);
+ const double r3 = r2 * r1;
+ if (r1 > 0.5 * m_radius)
{
const double rd = m_radius - r1;
res = m_k * rd * rd * rd * r3;
@@ -324,7 +324,7 @@ public:
else
{
const double rd = m_radius - r1;
- res = m_k * Real(2.0) * rd * rd * rd * r3 - m_c;
+ res = m_k * 2.0 * rd * rd * rd * r3 - m_c;
}
}
return res;
@@ -333,14 +333,14 @@ public:
///
/// \brief Get W(0)
///
- Real W0() const { return m_W0; }
+ double W0() const { return m_W0; }
protected:
- Real m_radius; ///> Kernel radius
- Real m_radius2; ///> Kernel radius squared
- Real m_k; ///> Kernel coefficient for W()
- Real m_c; ///> Kernel coefficient for W()
- Real m_W0; ///> Precomputed W(0)
+ double m_radius; ///> Kernel radius
+ double m_radius2; ///> Kernel radius squared
+ double m_k; ///> Kernel coefficient for W()
+ double m_c; ///> Kernel coefficient for W()
+ double m_W0; ///> Precomputed W(0)
};
///
@@ -350,7 +350,7 @@ protected:
template<int N>
class AdhesionKernel
{
-using VecXr = Eigen::Matrix<Real, N, 1>;
+using VecXd = Eigen::Matrix<double, N, 1>;
public:
AdhesionKernel()
@@ -361,31 +361,31 @@ public:
///
/// \brief Set the kernel radius
///
- void setRadius(const Real radius)
+ void setRadius(const double radius)
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
CHECK(N != 2) << "Unimplemented function";
- m_k = Real(0.007 / std::pow(m_radius, 3.25));
- m_W0 = W(VecXr::Zero());
+ m_k = 0.007 / std::pow(m_radius, 3.25);
+ m_W0 = W(VecXd::Zero());
}
///
/// \brief Compute weight value
/// W(r,h) = (0.007/h^3.25)(-4r^2/h + 6r -2h)^0.25 if h/2 < r <= h
///
- Real W(const Real r) const
+ double W(const double r) const
{
- Real res = 0.;
- const auto r2 = r * r;
+ double res = 0.;
+ const double r2 = r * r;
if (r2 <= m_radius2)
{
- const auto r = std::sqrt(r2);
- if (r > Real(0.5) * m_radius)
+ const double r = std::sqrt(r2);
+ if (r > 0.5 * m_radius)
{
- res = m_k * std::pow(-4.0 * r2 / m_radius + Real(6.0) * r - Real(2.0) * m_radius, 0.25);
+ res = m_k * std::pow(-4.0 * r2 / m_radius + 6.0 * r - 2.0 * m_radius, 0.25);
}
}
return res;
@@ -395,16 +395,16 @@ public:
/// \brief Compute weight value
/// W(r,h) = (0.007/h^3.25)(-4r^2/h + 6r -2h)^0.25 if h/2 < r <= h
///
- Real W(const VecXr& r) const
+ double W(const VecXd& r) const
{
- Real res = 0.;
- const auto r2 = r.squaredNorm();
+ double res = 0.;
+ const double r2 = r.squaredNorm();
if (r2 <= m_radius2)
{
- const auto r = std::sqrt(r2);
- if (r > Real(0.5) * m_radius)
+ const double r = std::sqrt(r2);
+ if (r > 0.5 * m_radius)
{
- res = m_k * std::pow(-4.0 * r2 / m_radius + Real(6.0) * r - Real(2.0) * m_radius, 0.25);
+ res = m_k * std::pow(-4.0 * r2 / m_radius + 6.0 * r - 2.0 * m_radius, 0.25);
}
}
return res;
@@ -413,13 +413,13 @@ public:
///
/// \brief Get W(0)
///
- Real W0() const { return m_W0; }
+ double W0() const { return m_W0; }
protected:
- Real m_radius; ///> Kernel radius
- Real m_radius2; ///> Kernel radius squared
- Real m_k; ///> Kernel coefficient for W()
- Real m_W0; ///> Precomputed W(0)
+ double m_radius; ///> Kernel radius
+ double m_radius2; ///> Kernel radius squared
+ double m_k; ///> Kernel coefficient for W()
+ double m_W0; ///> Precomputed W(0)
};
///
@@ -429,7 +429,7 @@ protected:
template<int N>
class ViscosityKernel
{
-using VecXr = Eigen::Matrix<Real, N, 1>;
+using VecXd = Eigen::Matrix<double, N, 1>;
public:
ViscosityKernel()
@@ -440,33 +440,33 @@ public:
///
/// \brief Set the kernel radius
///
- void setRadius(const Real radius)
+ void setRadius(const double radius)
{
m_radius = radius;
m_radius2 = radius * radius;
- m_k = Real(45.0 / PI) / (m_radius2 * m_radius2 * m_radius2);
+ m_k = (45.0 / PI) / (m_radius2 * m_radius2 * m_radius2);
}
///
/// \brief Compute laplacian
/// Laplace(r) = (45/PI/h^6) * (h - |r|)
///
- Real laplace(const VecXr& r) const
+ double laplace(const VecXd& r) const
{
- Real res = 0.;
- const auto r2 = r.squaredNorm();
+ double res = 0.;
+ const double r2 = r.squaredNorm();
if (r2 <= m_radius2)
{
- const auto d = std::sqrt(r2);
+ const double d = std::sqrt(r2);
res = m_k * (m_radius - d);
}
return res;
}
protected:
- Real m_radius; ///> Kernel radius
- Real m_radius2; ///> Kernel radius squared
- Real m_k; ///> Kernel coefficient for laplacian()
+ double m_radius; ///> Kernel radius
+ double m_radius2; ///> Kernel radius squared
+ double m_k; ///> Kernel coefficient for laplacian()
};
} // end namespace SPH
@@ -480,7 +480,7 @@ public:
///
/// \brief Initialize with kernel radius \p kernelRadius
///
- void initialize(const Real kernelRadius)
+ void initialize(const double kernelRadius)
{
m_poly6.setRadius(kernelRadius);
m_spiky.setRadius(kernelRadius);
@@ -491,27 +491,27 @@ public:
///
/// \brief Compute weight W(0) using poly6 kernel
///
- Real W0() const { return m_poly6.W0(); }
+ double W0() const { return m_poly6.W0(); }
///
/// \brief Compute weight W using poly6 kernel
///
- Real W(const Vec3r& r) const { return m_poly6.W(r); }
+ double W(const Vec3d& r) const { return m_poly6.W(r); }
///
/// \brief Compute gradW using spiky kernel
///
- Vec3r gradW(const Vec3r& r) const { return m_spiky.gradW(r); }
+ Vec3d gradW(const Vec3d& r) const { return m_spiky.gradW(r); }
///
/// \brief Compute laplacian using viscosity kernel
///
- Real laplace(const Vec3r& r) const { return m_viscosity.laplace(r); }
+ double laplace(const Vec3d& r) const { return m_viscosity.laplace(r); }
///
/// \brief Compute cohesion W using cohesion kernel
///
- Real cohesionW(const Vec3r& r) const { return m_cohesion.W(r); }
+ double cohesionW(const Vec3d& r) const { return m_cohesion.W(r); }
protected:
SPH::Poly6Kernel<3> m_poly6;
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp b/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp
index 694856143633968852e03d4116b34b6ba318221b..2dd5ca22dc3170ef3a3a67e5ba10c18a892d1b23 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHModel.cpp
@@ -27,10 +27,10 @@ limitations under the License.
namespace imstk
{
-SPHModelConfig::SPHModelConfig(const Real particleRadius)
+SPHModelConfig::SPHModelConfig(const double particleRadius)
{
// \todo Warning in all paths?
- if (std::abs(particleRadius) > Real(1.e-6))
+ if (std::abs(particleRadius) > 1.0e-6)
{
LOG_IF(WARNING, (particleRadius < 0)) << "Particle radius supplied is negative! Using absolute value of the supplied radius.";
m_particleRadius = std::abs(particleRadius);
@@ -43,9 +43,9 @@ SPHModelConfig::SPHModelConfig(const Real particleRadius)
initialize();
}
-SPHModelConfig::SPHModelConfig(const Real particleRadius, const Real speedOfSound, const Real restDensity)
+SPHModelConfig::SPHModelConfig(const double particleRadius, const double speedOfSound, const double restDensity)
{
- if (std::abs(particleRadius) > Real(1.e-6))
+ if (std::abs(particleRadius) > 1.0e-6)
{
LOG_IF(WARNING, (particleRadius < 0)) << "Particle radius supplied is negative! Using absolute value of the supplied radius.";
m_particleRadius = std::abs(particleRadius);
@@ -82,14 +82,14 @@ SPHModelConfig::initialize()
// Compute the derived quantities
m_particleRadiusSqr = m_particleRadius * m_particleRadius;
- m_particleMass = Real(std::pow(Real(2.0) * m_particleRadius, 3)) * m_restDensity * m_particleMassScale;
+ m_particleMass = std::pow(2.0 * m_particleRadius, 3) * m_restDensity * m_particleMassScale;
m_restDensitySqr = m_restDensity * m_restDensity;
- m_restDensityInv = Real(1.0) / m_restDensity;
+ m_restDensityInv = 1.0 / m_restDensity;
m_kernelRadius = m_particleRadius * m_kernelOverParticleRadiusRatio;
m_kernelRadiusSqr = m_kernelRadius * m_kernelRadius;
- m_pressureStiffness = m_restDensity * m_speedOfSound * m_speedOfSound / 7;
+ m_pressureStiffness = m_restDensity * m_speedOfSound * m_speedOfSound / 7.0;
}
SPHModel::SPHModel() : DynamicalModel<SPHState>(DynamicalModelType::SmoothedParticleHydrodynamics)
@@ -232,15 +232,15 @@ SPHModel::computeTimeStepSize()
m_dt = (this->m_timeStepSizeType == TimeSteppingType::Fixed) ? m_defaultDt : computeCFLTimeStepSize();
}
-Real
+double
SPHModel::computeCFLTimeStepSize()
{
auto maxVel = ParallelUtils::findMaxL2Norm(*getCurrentState()->getFullStepVelocities());
// dt = CFL * 2r / (speed of sound + max{|| v ||})
- Real timestep = maxVel > Real(1e-6) ?
- m_modelParameters->m_CFLFactor * (Real(2.0) * m_modelParameters->m_particleRadius / (m_modelParameters->m_speedOfSound + maxVel)) :
- m_modelParameters->m_maxTimestep;
+ double timestep = maxVel > 1.0e-6 ?
+ m_modelParameters->m_CFLFactor * (2.0 * m_modelParameters->m_particleRadius / (m_modelParameters->m_speedOfSound + maxVel)) :
+ m_modelParameters->m_maxTimestep;
// clamp the time step size to be within a given range
if (timestep > m_modelParameters->m_maxTimestep)
@@ -270,7 +270,7 @@ SPHModel::findParticleNeighbors()
void
SPHModel::computeNeighborRelativePositions()
{
- auto computeRelativePositions = [&](const Vec3r& ppos, const std::vector<size_t>& neighborList,
+ auto computeRelativePositions = [&](const Vec3d& ppos, const std::vector<size_t>& neighborList,
const VecDataArray<double, 3>& allPositions, std::vector<NeighborInfo>& neighborInfo)
{
for (const size_t q : neighborList)
@@ -366,7 +366,7 @@ SPHModel::computeDensity()
return; // the particle has no neighbor
}
- Real pdensity = 0.0;
+ double pdensity = 0.0;
for (const auto& qInfo : neighborInfo)
{
pdensity += m_kernels.W(qInfo.xpq);
@@ -418,7 +418,7 @@ SPHModel::normalizeDensity()
}
const auto& fluidNeighborList = neighborLists[p];
- Real tmp = 0;
+ double tmp = 0.0;
for (size_t i = 0; i < fluidNeighborList.size(); ++i)
{
@@ -452,7 +452,7 @@ SPHModel::computePressureAcceleration()
return;
}
- Vec3r accel(0, 0, 0);
+ Vec3d accel = Vec3d::Zero();
const std::vector<NeighborInfo>& neighborInfo = neighborInfos[p];
if (neighborInfo.size() <= 1)
{
@@ -504,20 +504,20 @@ SPHModel::computeViscosity()
const std::vector<NeighborInfo>& neighborInfo = neighborInfos[p];
if (neighborInfo.size() <= 1)
{
- neighborVelContr[p] = Vec3r(0, 0, 0);
- viscousAccels[p] = Vec3r(0, 0, 0);
+ neighborVelContr[p] = Vec3d::Zero();
+ viscousAccels[p] = Vec3d::Zero();
return;
}
- Vec3r neighborVelContributions(0, 0, 0);
- Vec3r neighborVelContributionsNumerator(0, 0, 0);
- Real neighborVelContributionsDenominator = 0;
- Vec3r particleShifts(0, 0, 0);
+ Vec3d neighborVelContributions = Vec3d::Zero();
+ Vec3d neighborVelContributionsNumerator = Vec3d::Zero();
+ double neighborVelContributionsDenominator = 0.0;
+ Vec3d particleShifts = Vec3d::Zero();
const Vec3d& pvel = halfStepVelocities[p];
const std::vector<size_t>& fluidNeighborList = neighborLists[p];
- Vec3r diffuseFluid(0, 0, 0);
+ Vec3d diffuseFluid = Vec3d::Zero();
for (size_t i = 0; i < fluidNeighborList.size(); ++i)
{
const auto q = fluidNeighborList[i];
@@ -525,7 +525,7 @@ SPHModel::computeViscosity()
const auto& qInfo = neighborInfo[i];
const auto r = qInfo.xpq;
const auto qdensity = qInfo.density;
- diffuseFluid += (Real(1.0) / qdensity) * m_kernels.laplace(r) * (qvel - pvel);
+ diffuseFluid += (1.0 / qdensity) * m_kernels.laplace(r) * (qvel - pvel);
neighborVelContributionsNumerator += (qvel - pvel) * m_kernels.W(r);
neighborVelContributionsDenominator += m_kernels.W(r);
@@ -572,7 +572,7 @@ SPHModel::computeSurfaceTension()
const auto& qInfo = neighborInfo[i];
const auto r = qInfo.xpq;
const auto qdensity = qInfo.density;
- n += (Real(1.0) / qdensity) * m_kernels.gradW(r);
+ n += (1.0 / qdensity) * m_kernels.gradW(r);
}
n *= m_modelParameters->m_kernelRadius * m_modelParameters->m_particleMass;
@@ -595,40 +595,40 @@ SPHModel::computeSurfaceTension()
return;
}
- const auto& fluidNeighborList = neighborLists[p];
+ const std::vector<size_t>& fluidNeighborList = neighborLists[p];
if (fluidNeighborList.size() <= 1)
{
return; // the particle has no neighbor
}
- const auto ni = surfaceNormals[p];
- const auto pdensity = densities[p];
+ const Vec3d& ni = surfaceNormals[p];
+ const double pdensity = densities[p];
const auto& neighborInfo = neighborInfos[p];
- Vec3r accel(0, 0, 0);
+ Vec3d accel = Vec3d::Zero();
for (size_t i = 0; i < fluidNeighborList.size(); ++i)
{
- const auto q = fluidNeighborList[i];
+ const size_t q = fluidNeighborList[i];
if (p == q)
{
continue;
}
- const auto& qInfo = neighborInfo[i];
- const auto qdensity = qInfo.density;
+ const NeighborInfo& qInfo = neighborInfo[i];
+ const double qdensity = qInfo.density;
// Correction factor
- const auto K_ij = Real(2) * m_modelParameters->m_restDensity / (pdensity + qdensity);
+ const double K_ij = 2.0 * m_modelParameters->m_restDensity / (pdensity + qdensity);
// Cohesion acc
- const auto r = qInfo.xpq;
- const auto d2 = r.squaredNorm();
- if (d2 > Real(1e-20))
+ const Vec3d& r = qInfo.xpq;
+ const double d2 = r.squaredNorm();
+ if (d2 > 1.0e-20)
{
accel -= K_ij * m_modelParameters->m_particleMass * (r / std::sqrt(d2)) * m_kernels.cohesionW(r);
}
// Curvature acc
- const auto nj = surfaceNormals[q];
+ const Vec3d& nj = surfaceNormals[q];
accel -= K_ij * (ni - nj);
}
@@ -661,7 +661,7 @@ SPHModel::sumAccels()
}
void
-SPHModel::updateVelocity(const Real timestep)
+SPHModel::updateVelocity(const double timestep)
{
VecDataArray<double, 3>& halfStepVelocities = *getCurrentState()->getHalfStepVelocities();
VecDataArray<double, 3>& fullStepVelocities = *getCurrentState()->getFullStepVelocities();
@@ -698,7 +698,7 @@ SPHModel::updateVelocity(const Real timestep)
}
void
-SPHModel::moveParticles(const Real timestep)
+SPHModel::moveParticles(const double timestep)
{
//ParallelUtils::parallelFor(getState().getNumParticles(),
// [&](const size_t p) {
@@ -718,8 +718,8 @@ SPHModel::moveParticles(const Real timestep)
continue;
}
- Vec3r oldPosition = positions[p];
- Vec3r newPosition = oldPosition + particleShift[p] * timestep + (halfStepVelocities[p] + neighborVelContr[p]) * timestep;
+ Vec3d oldPosition = positions[p];
+ Vec3d newPosition = oldPosition + particleShift[p] * timestep + (halfStepVelocities[p] + neighborVelContr[p]) * timestep;
positions[p] = newPosition;
@@ -767,15 +767,15 @@ SPHModel::moveParticles(const Real timestep)
m_timeStepCount++;
}
-Real
+double
SPHModel::particlePressure(const double density)
{
const double d = density / m_modelParameters->m_restDensity;
const double d2 = d * d;
const double d4 = d2 * d2;
- const Real error = m_modelParameters->m_pressureStiffness * (d4 * d2 * d - Real(1.0));
+ const double error = m_modelParameters->m_pressureStiffness * (d4 * d2 * d - 1.0);
// clamp pressure error to zero to maintain stability
- return error > Real(0) ? error : Real(0);
+ return error > 0.0 ? error : 0.0;
}
void
@@ -789,7 +789,7 @@ SPHModel::setInitialVelocities(const size_t numParticles, const Vec3d& initialVe
&& (m_sphBoundaryConditions->getParticleTypes()[p] == SPHBoundaryConditions::ParticleType::Buffer
|| m_sphBoundaryConditions->getParticleTypes()[p] == SPHBoundaryConditions::ParticleType::Wall))
{
- m_initialVelocities->push_back(Vec3d(0.0, 0.0, 0.0));
+ m_initialVelocities->push_back(Vec3d::Zero());
}
else
{
diff --git a/Source/DynamicalModels/ObjectModels/imstkSPHModel.h b/Source/DynamicalModels/ObjectModels/imstkSPHModel.h
index 9181d1ebc32b09b845b5c58fc1a5ed09a9c52383..a8775664c468c854e6aafa13585d903e1b116c8d 100644
--- a/Source/DynamicalModels/ObjectModels/imstkSPHModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkSPHModel.h
@@ -41,48 +41,48 @@ private:
void initialize();
public:
- explicit SPHModelConfig(const Real particleRadius);
- explicit SPHModelConfig(const Real particleRadius, const Real speedOfSound, const Real restDensity);
+ explicit SPHModelConfig(const double particleRadius);
+ explicit SPHModelConfig(const double particleRadius, const double speedOfSound, const double restDensity);
/// \todo Move this to solver or time integrator in the future
- Real m_minTimestep = Real(1e-6);
- Real m_maxTimestep = Real(1e-3);
- Real m_CFLFactor = Real(1.0);
+ double m_minTimestep = 1.0e-6;
+ double m_maxTimestep = 1.0e-3;
+ double m_CFLFactor = 1.0;
// particle parameters
- Real m_particleRadius = Real(0);
- Real m_particleRadiusSqr = Real(0); ///> \note derived quantity
+ double m_particleRadius = 0.0;
+ double m_particleRadiusSqr = 0.0; ///> \note derived quantity
// material parameters
- Real m_restDensity = Real(1000);
- Real m_restDensitySqr = Real(1000000.0); ///> \note derived quantity
- Real m_restDensityInv = Real(1.0 / 1000.0); ///> \note derived quantity
- Real m_particleMass = Real(1);
- Real m_particleMassScale = Real(1.0); ///> scale particle mass to a smaller value to maintain stability
- Real m_eta = Real(0.5); ///> proportion of position change due to neighbors velocity (XSPH method)
+ double m_restDensity = 1000.0;
+ double m_restDensitySqr = 1000000.0; ///> \note derived quantity
+ double m_restDensityInv = 1.0 / 1000.0; ///> \note derived quantity
+ double m_particleMass = 1.0;
+ double m_particleMassScale = 1.0; ///> scale particle mass to a smaller value to maintain stability
+ double m_eta = 0.5; ///> proportion of position change due to neighbors velocity (XSPH method)
bool m_bNormalizeDensity = false;
bool m_bDensityWithBoundary = false;
// pressure
- Real m_pressureStiffness = Real(50000.0);
+ double m_pressureStiffness = 50000.0;
// viscosity and surface tension/cohesion
- Real m_dynamicViscosityCoeff = Real(1e-2);
- Real m_viscosityBoundary = Real(1e-5);
- Real m_surfaceTensionStiffness = Real(1);
- Real m_frictionBoundary = Real(0.1);
+ double m_dynamicViscosityCoeff = 1.0e-2;
+ double m_viscosityBoundary = 1.0e-5;
+ double m_surfaceTensionStiffness = 1.0;
+ double m_frictionBoundary = 0.1;
// kernel properties
- Real m_kernelOverParticleRadiusRatio = Real(4.0);
- Real m_kernelRadius; ///> \note derived quantity
- Real m_kernelRadiusSqr; ///> \note derived quantity
+ double m_kernelOverParticleRadiusRatio = 4.0;
+ double m_kernelRadius; ///> \note derived quantity
+ double m_kernelRadiusSqr; ///> \note derived quantity
// gravity
- Vec3r m_gravity = Vec3r(0, -9.81, 0);
+ Vec3d m_gravity = Vec3d(0.0, -9.81, 0.0);
// sound speed
- Real m_speedOfSound = 18.7;
+ double m_speedOfSound = 18.7;
// neighbor search
NeighborSearch::Method m_NeighborSearchMethod = NeighborSearch::Method::UniformGridBasedSearch;
@@ -160,7 +160,7 @@ public:
void setBoundaryConditions(std::shared_ptr<SPHBoundaryConditions> sphBoundaryConditions) { m_sphBoundaryConditions = sphBoundaryConditions; }
std::shared_ptr<SPHBoundaryConditions> getBoundaryConditions() { return m_sphBoundaryConditions; }
- void setRestDensity(const Real restDensity) { m_modelParameters->m_restDensity = restDensity; }
+ void setRestDensity(const double restDensity) { m_modelParameters->m_restDensity = restDensity; }
std::shared_ptr<TaskNode> getFindParticleNeighborsNode() const { return m_findParticleNeighborsNode; }
std::shared_ptr<TaskNode> getComputeDensityNode() const { return m_computeDensityNode; }
@@ -188,7 +188,7 @@ private:
///
/// \brief Compute time step size based on CFL condition
///
- Real computeCFLTimeStepSize();
+ double computeCFLTimeStepSize();
///
/// \brief Find the neighbors for each particle
@@ -229,7 +229,7 @@ private:
///
/// \brief Update particle velocities due to pressure, viscous, and surface tension forces
///
- void updateVelocity(const Real timestep);
+ void updateVelocity(const double timestep);
///
/// \brief Compute viscosity
@@ -246,7 +246,7 @@ private:
///
/// \brief Move particles
///
- void moveParticles(const Real timestep);
+ void moveParticles(const double timestep);
//void computePressureOutlet();
diff --git a/Source/DynamicalModels/ObjectStates/imstkRigidBodyState2.h b/Source/DynamicalModels/ObjectStates/imstkRigidBodyState2.h
index 327f1610b06c84f615cf427abffa1a6738543468..6390757b3ed2dc23b55f3af9eed89baf0c863e39 100644
--- a/Source/DynamicalModels/ObjectStates/imstkRigidBodyState2.h
+++ b/Source/DynamicalModels/ObjectStates/imstkRigidBodyState2.h
@@ -53,8 +53,8 @@ public:
const std::vector<bool>& getIsStatic() const { return m_isStatic; }
std::vector<bool>& getIsStatic() { return m_isStatic; }
- const StdVectorOfReal& getInvMasses() const { return m_invMasses; }
- StdVectorOfReal& getInvMasses() { return m_invMasses; }
+ const std::vector<double>& getInvMasses() const { return m_invMasses; }
+ std::vector<double>& getInvMasses() { return m_invMasses; }
const StdVectorOfMat3d& getInvIntertiaTensors() const { return m_invIntertiaTensors; }
StdVectorOfMat3d& getInvIntertiaTensors() { return m_invIntertiaTensors; }
@@ -97,8 +97,8 @@ public:
}
private:
- StdVectorOfReal m_invMasses;
- StdVectorOfMat3d m_invIntertiaTensors;
+ std::vector<double> m_invMasses;
+ StdVectorOfMat3d m_invIntertiaTensors;
StdVectorOfVec3d m_positions;
StdVectorOfQuatd m_orientations;
diff --git a/Source/DynamicalModels/ObjectStates/imstkSPHState.h b/Source/DynamicalModels/ObjectStates/imstkSPHState.h
index 686a2856e6d71cc6a05b1fa56f6fbcecf0908f0a..fbb8a625caa5ccea04b2e0937b337a430d1d6a84 100644
--- a/Source/DynamicalModels/ObjectStates/imstkSPHState.h
+++ b/Source/DynamicalModels/ObjectStates/imstkSPHState.h
@@ -34,8 +34,8 @@ template<typename T, int N> class VecDataArray;
///
struct NeighborInfo
{
- Vec3r xpq; ///> relative position: xpq = x_p - x_q
- Real density; ///> density of neighbor particle q
+ Vec3d xpq; ///> relative position: xpq = x_p - x_q
+ double density; ///> density of neighbor particle q
};
///
@@ -147,6 +147,6 @@ private:
std::vector<std::vector<size_t>> m_NeighborLists; ///> store a list of neighbors for each particle, updated each time step
std::vector<std::vector<size_t>> m_BDNeighborLists; ///> store a list of boundary particle neighbors for each particle, updated each time step
- std::vector<std::vector<NeighborInfo>> m_NeighborInfo; ///> store a list of Vec4r(Vec3r(relative position), density) for neighbors, including boundary particle
+ std::vector<std::vector<NeighborInfo>> m_NeighborInfo; ///> store a list of Vec4d(Vec3d(relative position), density) for neighbors, including boundary particle
};
} // end namespace imstk
diff --git a/Source/DynamicalModels/Testing/CMakeLists.txt b/Source/DynamicalModels/Testing/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d6344deee27a7a22b30a790f5ca479ef4b8fd491
--- /dev/null
+++ b/Source/DynamicalModels/Testing/CMakeLists.txt
@@ -0,0 +1,2 @@
+include(imstkAddTest)
+imstk_add_test( DynamicalModels )
\ No newline at end of file
diff --git a/Source/DynamicalModels/Testing/imstkPbdConstraintFunctorTest.cpp b/Source/DynamicalModels/Testing/imstkPbdConstraintFunctorTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..719d348797b91de5816df7c1ba514094a3426916
--- /dev/null
+++ b/Source/DynamicalModels/Testing/imstkPbdConstraintFunctorTest.cpp
@@ -0,0 +1,316 @@
+/*=========================================================================
+
+ 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 "gtest/gtest.h"
+
+#include "imstkLineMesh.h"
+#include "imstkPbdConstraintFunctor.h"
+#include "imstkVecDataArray.h"
+
+using namespace imstk;
+
+///
+/// \brief Test that the bending correct constraint was generated
+///
+TEST(imstkPbdConstraintFunctorTest, TestBendingConstraintStride1Generation)
+{
+ // Create mesh for generation
+ auto lineMesh = std::make_shared<LineMesh>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(3);
+ (*vertices)[0] = Vec3d(-0.5, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(0.0, 0.0, 0.0);
+ (*vertices)[2] = Vec3d(0.5, 0.0, 0.0);
+ auto indices = std::make_shared<VecDataArray<int, 2>>(2);
+ (*indices)[0] = Vec2i(0, 1);
+ (*indices)[1] = Vec2i(1, 2);
+ lineMesh->initialize(vertices, indices);
+
+ // Create functor
+ PbdBendConstraintFunctor constraintFunctor;
+ constraintFunctor.setStiffness(1e20);
+ constraintFunctor.setStride(1);
+ constraintFunctor.setGeometry(lineMesh);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdBendConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::Bend);
+ EXPECT_EQ(constraint->getStiffness(), 1e20);
+ EXPECT_EQ(constraint->getVertexIds().size(), 3);
+ EXPECT_EQ(constraint->getVertexIds()[0], 0);
+ EXPECT_EQ(constraint->getVertexIds()[1], 1);
+ EXPECT_EQ(constraint->getVertexIds()[2], 2);
+}
+
+///
+/// \brief Test that the correct bending constraint was generated with differing stride
+///
+TEST(imstkPbdConstraintFunctorTest, TestBendingConstraintStride2Generation)
+{
+ // Create mesh for generation
+ auto lineMesh = std::make_shared<LineMesh>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(5);
+ (*vertices)[0] = Vec3d(-1.0, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(-0.5, 0.0, 0.0);
+ (*vertices)[2] = Vec3d(0.0, 0.0, 0.0);
+ (*vertices)[3] = Vec3d(0.5, 0.0, 0.0);
+ (*vertices)[4] = Vec3d(1.0, 0.0, 0.0);
+ auto indices = std::make_shared<VecDataArray<int, 2>>(4);
+ (*indices)[0] = Vec2i(0, 1);
+ (*indices)[1] = Vec2i(1, 2);
+ (*indices)[2] = Vec2i(2, 3);
+ (*indices)[3] = Vec2i(3, 4);
+ lineMesh->initialize(vertices, indices);
+
+ // Create functor
+ PbdBendConstraintFunctor constraintFunctor;
+ constraintFunctor.setStiffness(1e20);
+ constraintFunctor.setStride(2);
+ constraintFunctor.setGeometry(lineMesh);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdBendConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::Bend);
+ EXPECT_EQ(constraint->getStiffness(), 1e20);
+ EXPECT_EQ(constraint->getVertexIds().size(), 3);
+ EXPECT_EQ(constraint->getVertexIds()[0], 0);
+ EXPECT_EQ(constraint->getVertexIds()[1], 2);
+ EXPECT_EQ(constraint->getVertexIds()[2], 4);
+}
+
+///
+/// \brief Test that the correct distance constraint was generated
+///
+TEST(imstkPbdConstraintFunctorTest, TestDistanceConstraintGeneration)
+{
+ // Create mesh for generation
+ auto lineMesh = std::make_shared<LineMesh>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(2);
+ (*vertices)[0] = Vec3d(-0.5, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(0.0, 0.0, 0.0);
+ auto indices = std::make_shared<VecDataArray<int, 2>>(1);
+ (*indices)[0] = Vec2i(0, 1);
+ lineMesh->initialize(vertices, indices);
+
+ // Create functor
+ PbdDistanceConstraintFunctor constraintFunctor;
+ constraintFunctor.setStiffness(1.0e3);
+ constraintFunctor.setGeometry(lineMesh);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdDistanceConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::Distance);
+ EXPECT_EQ(constraint->getStiffness(), 1.0e3);
+ EXPECT_EQ(constraint->getVertexIds().size(), 2);
+ EXPECT_EQ(constraint->getVertexIds()[0], 0);
+ EXPECT_EQ(constraint->getVertexIds()[1], 1);
+}
+
+///
+/// \brief Test that the correct pbd FEM tetrahedral constraint was generated
+///
+TEST(imstkPbdConstraintFunctorTest, TestFEMTetConstraintGeneration)
+{
+ auto tetMesh = std::make_shared<TetrahedralMesh>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(4);
+ (*vertices)[0] = Vec3d(0.0, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(1.0, 0.0, 0.0);
+ (*vertices)[2] = Vec3d(0.0, 1.0, 0.0);
+ (*vertices)[3] = Vec3d(0.0, 0.0, 1.0);
+ auto indices = std::make_shared<VecDataArray<int, 4>>(1);
+ (*indices)[0] = Vec4i(0, 1, 2, 3);
+ tetMesh->initialize(vertices, indices);
+
+ // Create functor
+ PbdFemConstraintFunctor constraintFunctor;
+ constraintFunctor.setMaterialType(PbdFEMTetConstraint::MaterialType::Corotation);
+ auto FeConfig = std::make_shared<PbdFEMConstraintConfig>(0.0, 0.0, 1000.0, 0.2);
+ constraintFunctor.setFemConfig(FeConfig);
+ constraintFunctor.setGeometry(tetMesh);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdFEMTetConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::FEMTet);
+ EXPECT_EQ(constraint->m_material, PbdFEMTetConstraint::MaterialType::Corotation);
+ EXPECT_EQ(constraint->m_config->m_mu, 0.0);
+ EXPECT_EQ(constraint->m_config->m_lambda, 0.0);
+ EXPECT_EQ(constraint->m_config->m_YoungModulus, 1000.0);
+ EXPECT_EQ(constraint->m_config->m_PoissonRatio, 0.2);
+ EXPECT_EQ(constraint->getVertexIds().size(), 4);
+ EXPECT_EQ(constraint->getVertexIds()[0], 0);
+ EXPECT_EQ(constraint->getVertexIds()[1], 1);
+ EXPECT_EQ(constraint->getVertexIds()[2], 2);
+ EXPECT_EQ(constraint->getVertexIds()[3], 3);
+}
+
+///
+/// \brief Test that the correct pbd volume constraint was generated
+///
+TEST(imstkPbdConstraintFunctorTest, TestVolumeConstraintGeneration)
+{
+ auto tetMesh = std::make_shared<TetrahedralMesh>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(4);
+ (*vertices)[0] = Vec3d(0.0, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(1.0, 0.0, 0.0);
+ (*vertices)[2] = Vec3d(0.0, 1.0, 0.0);
+ (*vertices)[3] = Vec3d(0.0, 0.0, 1.0);
+ auto indices = std::make_shared<VecDataArray<int, 4>>(1);
+ (*indices)[0] = Vec4i(0, 1, 2, 3);
+ tetMesh->initialize(vertices, indices);
+
+ // Create functor
+ PbdVolumeConstraintFunctor constraintFunctor;
+ constraintFunctor.setStiffness(1.0e4);
+ constraintFunctor.setGeometry(tetMesh);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdVolumeConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::Volume);
+ EXPECT_EQ(constraint->getStiffness(), 1.0e4);
+ EXPECT_EQ(constraint->getVertexIds().size(), 4);
+ EXPECT_EQ(constraint->getVertexIds()[0], 0);
+ EXPECT_EQ(constraint->getVertexIds()[1], 1);
+ EXPECT_EQ(constraint->getVertexIds()[2], 2);
+ EXPECT_EQ(constraint->getVertexIds()[3], 3);
+}
+
+///
+/// \brief Test that the correct pbd area constraint was generated
+///
+TEST(imstkPbdConstraintFunctorTest, TestAreaConstraintGeneration)
+{
+ auto surfMesh = std::make_shared<SurfaceMesh>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(4);
+ (*vertices)[0] = Vec3d(0.0, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(1.0, 0.0, 0.0);
+ (*vertices)[2] = Vec3d(0.0, 1.0, 0.0);
+ (*vertices)[3] = Vec3d(0.0, 0.0, 1.0);
+ auto indices = std::make_shared<VecDataArray<int, 3>>(1);
+ (*indices)[0] = Vec3i(0, 1, 2);
+ surfMesh->initialize(vertices, indices);
+
+ // Create functor
+ PbdAreaConstraintFunctor constraintFunctor;
+ constraintFunctor.setStiffness(1.0e4);
+ constraintFunctor.setGeometry(surfMesh);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdAreaConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::Area);
+ EXPECT_EQ(constraint->getStiffness(), 1.0e4);
+ EXPECT_EQ(constraint->getVertexIds().size(), 3);
+ EXPECT_EQ(constraint->getVertexIds()[0], 0);
+ EXPECT_EQ(constraint->getVertexIds()[1], 1);
+ EXPECT_EQ(constraint->getVertexIds()[2], 2);
+}
+
+///
+/// \brief Test that the correct pbd constant density constraint was generated
+///
+TEST(imstkPbdConstraintFunctorTest, TestConstDensityConstraintGeneration)
+{
+ auto points = std::make_shared<PointSet>();
+ auto vertices = std::make_shared<VecDataArray<double, 3>>(4);
+ (*vertices)[0] = Vec3d(0.0, 0.0, 0.0);
+ (*vertices)[1] = Vec3d(1.0, 0.0, 0.0);
+ (*vertices)[2] = Vec3d(0.0, 1.0, 0.0);
+ (*vertices)[3] = Vec3d(0.0, 0.0, 1.0);
+ points->initialize(vertices);
+
+ // Create functor
+ PbdConstantDensityConstraintFunctor constraintFunctor;
+ constraintFunctor.setStiffness(1.0e4);
+ constraintFunctor.setGeometry(points);
+
+ // Fill container
+ PbdConstraintContainer container;
+ constraintFunctor(container);
+
+ // Check that constraint got generated
+ EXPECT_EQ(container.getConstraints().size(), 1);
+
+ // Check that correct constraint type got generated
+ auto constraint = std::dynamic_pointer_cast<PbdConstantDensityConstraint>(container.getConstraints()[0]);
+ EXPECT_NE(constraint, nullptr);
+
+ // Check constraint generated between correct elements and with correct values
+ EXPECT_EQ(constraint->getType(), PbdConstraint::Type::ConstantDensity);
+ EXPECT_EQ(constraint->getVertexIds().size(), 0);
+}
\ No newline at end of file
diff --git a/Source/Geometry/Analytic/imstkOrientedBox.cpp b/Source/Geometry/Analytic/imstkOrientedBox.cpp
index 7e4a0908df7a7efdb0a4b40e253d5e42a20f6f05..5d282d8181a90ad8fbb7bd583c288444e1d3c830 100644
--- a/Source/Geometry/Analytic/imstkOrientedBox.cpp
+++ b/Source/Geometry/Analytic/imstkOrientedBox.cpp
@@ -92,9 +92,8 @@ OrientedBox::getFunctionValue(const Vec3d& pos) const
const Mat3d rot = m_orientationPostTransform.toRotationMatrix();
const Vec3d& extents = m_extentsPostTransform;
- const Vec3d diff = (pos - m_positionPostTransform);
- const Mat3d rotInv = rot.transpose();
- const Vec3d proj = rotInv * diff; // dot product/project onto each axes
+ const Vec3d diff = (pos - m_positionPostTransform);
+ const Vec3d proj = rot * diff; // dot product/project onto each axes
bool inside[3] =
{
diff --git a/Source/Geometry/Testing/imstkGeometryTest.cpp b/Source/Geometry/Testing/imstkGeometryTest.cpp
index cdb72a67e803b9a56951fdc8592fb91e311a60e8..6cfbaf2902dab78f623d8b0fef5830c5883f5403 100644
--- a/Source/Geometry/Testing/imstkGeometryTest.cpp
+++ b/Source/Geometry/Testing/imstkGeometryTest.cpp
@@ -19,74 +19,89 @@
=========================================================================*/
-#include "gtest/gtest.h"
+#include "imstkGeometry.h"
-#include "imstkPlane.h"
+#include <gtest/gtest.h>
using namespace imstk;
-class imstkGeometryTest : public ::testing::Test
+namespace
{
-protected:
- Plane m_geometry; // Can't use imstk::Geometry since pure virtual. Should use google mock class.
+class MockGeometry : public Geometry
+{
+public:
+ MockGeometry() : Geometry("MockGeometry") { }
+
+ const std::string getTypeName() const override { return "MockGeometry"; }
};
+}
-TEST_F(imstkGeometryTest, GetSetScaling)
+TEST(imstkGeometryTest, GetSetScaling)
{
- m_geometry.setScaling(Vec3d(2.0, 2.0, 2.0));
- EXPECT_EQ(m_geometry.getScaling(), Vec3d(2.0, 2.0, 2.0));
+ MockGeometry geometry;
+
+ geometry.setScaling(Vec3d(2.0, 2.0, 2.0));
+ EXPECT_EQ(geometry.getScaling(), Vec3d(2.0, 2.0, 2.0));
- m_geometry.setScaling(Vec3d(0.003, 0.003, 0.003));
- EXPECT_EQ(m_geometry.getScaling(), Vec3d(0.003, 0.003, 0.003));
+ geometry.setScaling(Vec3d(0.003, 0.003, 0.003));
+ EXPECT_EQ(geometry.getScaling(), Vec3d(0.003, 0.003, 0.003));
- m_geometry.setScaling(Vec3d(400000000.0, 400000000.0, 400000000.0));
- EXPECT_EQ(m_geometry.getScaling(), Vec3d(400000000.0, 400000000.0, 400000000.0));
+ geometry.setScaling(Vec3d(400000000.0, 400000000.0, 400000000.0));
+ EXPECT_EQ(geometry.getScaling(), Vec3d(400000000.0, 400000000.0, 400000000.0));
}
-TEST_F(imstkGeometryTest, GetSetTranslation)
+TEST(imstkGeometryTest, GetSetTranslation)
{
+ MockGeometry geometry;
+
auto p1 = Vec3d(12.0, 0.0005, -400000.0);
auto p2 = Vec3d(-500.0, 30.0, 0.23);
- m_geometry.setTranslation(p1);
- EXPECT_EQ(m_geometry.getTranslation(), p1);
+ geometry.setTranslation(p1);
+ EXPECT_EQ(geometry.getTranslation(), p1);
- m_geometry.setTranslation(p2);
- EXPECT_EQ(m_geometry.getTranslation(), p2);
+ geometry.setTranslation(p2);
+ EXPECT_EQ(geometry.getTranslation(), p2);
- m_geometry.setTranslation(p1[0], p1[1], p1[2]);
- EXPECT_EQ(m_geometry.getTranslation(), p1);
+ geometry.setTranslation(p1[0], p1[1], p1[2]);
+ EXPECT_EQ(geometry.getTranslation(), p1);
- m_geometry.setTranslation(p2[0], p2[1], p2[2]);
- EXPECT_EQ(m_geometry.getTranslation(), p2);
+ geometry.setTranslation(p2[0], p2[1], p2[2]);
+ EXPECT_EQ(geometry.getTranslation(), p2);
}
-TEST_F(imstkGeometryTest, GetSetRotation1)
+TEST(imstkGeometryTest, GetSetRotation1)
{
+ MockGeometry geometry;
+
// NOTE: '==' not defined for Eigen::Quaternion, using 'isApprox'.
// See https://forum.kde.org/viewtopic.php?f=74&t=118598
// Rotation is normalized so comparing pre/post requires starting orientation
// and angles < 360deg. As well as normalized axes.
const Quatd q1 = Quatd(Rotd(4.1, Vec3d(12.0, 0.0, -0.5).normalized()));
- m_geometry.setRotation(q1);
- EXPECT_TRUE(Quatd(m_geometry.getRotation()).isApprox(q1));
+ geometry.setRotation(q1);
+ EXPECT_TRUE(Quatd(geometry.getRotation()).isApprox(q1));
}
-TEST_F(imstkGeometryTest, GetSetRotation2)
+TEST(imstkGeometryTest, GetSetRotation2)
{
+ MockGeometry geometry;
+
// NOTE: '==' not defined for Eigen::Quaternion, using 'isApprox'.
// See https://forum.kde.org/viewtopic.php?f=74&t=118598
// Rotation is normalized so comparing pre/post requires starting orientation
// and angles < 360deg. As well as normalized axes.
const Mat3d mat2 = Mat3d(Rotd(0.43, Vec3d(4000.0, -1.0, 0.0).normalized()));
- m_geometry.setRotation(mat2);
- EXPECT_TRUE(m_geometry.getRotation().isApprox(mat2));
+ geometry.setRotation(mat2);
+ EXPECT_TRUE(geometry.getRotation().isApprox(mat2));
}
-TEST_F(imstkGeometryTest, GetSetRotation3)
+TEST(imstkGeometryTest, GetSetRotation3)
{
+ MockGeometry geometry;
+
// NOTE: '==' not defined for Eigen::Quaternion, using 'isApprox'.
// See https://forum.kde.org/viewtopic.php?f=74&t=118598
@@ -95,6 +110,6 @@ TEST_F(imstkGeometryTest, GetSetRotation3)
const double angle = 1.57;
const Vec3d axes = Vec3d(-0.0, 100.0, 2000000.0).normalized();
const Mat3d mat3 = Mat3d(Rotd(angle, axes));
- m_geometry.setRotation(axes, angle);
- EXPECT_TRUE(m_geometry.getRotation().isApprox(mat3));
+ geometry.setRotation(axes, angle);
+ EXPECT_TRUE(geometry.getRotation().isApprox(mat3));
}
diff --git a/Source/Geometry/Testing/imstkImageDataTest.cpp b/Source/Geometry/Testing/imstkImageDataTest.cpp
index 8624544c1ab765e295a54ca3d52cf0b4a2d7f5f8..32dfe59a9e1c3a4b4c5a7ad8f401f25310886de1 100644
--- a/Source/Geometry/Testing/imstkImageDataTest.cpp
+++ b/Source/Geometry/Testing/imstkImageDataTest.cpp
@@ -19,13 +19,11 @@
=========================================================================*/
-#include "gtest/gtest.h"
-
-#include <memory>
-
#include "imstkImageData.h"
#include "imstkVecDataArray.h"
-#include "imstkMath.h"
+
+#include <gtest/gtest.h>
+#include <memory>
using namespace imstk;
namespace
diff --git a/Source/Geometry/Testing/imstkOrientedBoxTest.cpp b/Source/Geometry/Testing/imstkOrientedBoxTest.cpp
index 65e377c7e3b765fcc8bb2b663ba235471482d376..dd71f011a37eab17dd74da4cb3e32d56783271c8 100644
--- a/Source/Geometry/Testing/imstkOrientedBoxTest.cpp
+++ b/Source/Geometry/Testing/imstkOrientedBoxTest.cpp
@@ -25,41 +25,38 @@
using namespace imstk;
-class imstkOrientedBoxTest : public ::testing::Test
+TEST(imstkOrientedBoxTest, SetGetWidth)
{
-protected:
- OrientedBox m_box;
-};
-
-TEST_F(imstkOrientedBoxTest, SetGetWidth)
-{
- m_box.setExtents(1.0, 1.0, 1.0);
- const Vec3d extents = m_box.getExtents();
+ OrientedBox box;
+ box.setExtents(1.0, 1.0, 1.0);
+ const Vec3d extents = box.getExtents();
EXPECT_DOUBLE_EQ(1.0, extents[0]);
EXPECT_DOUBLE_EQ(1.0, extents[1]);
EXPECT_DOUBLE_EQ(1.0, extents[2]);
}
-TEST_F(imstkOrientedBoxTest, GetVolume)
+TEST(imstkOrientedBoxTest, GetVolume)
{
- m_box.setExtents(1.0, 1.0, 1.0);
- EXPECT_DOUBLE_EQ(8, m_box.getVolume());
+ OrientedBox box;
+ box.setExtents(1.0, 1.0, 1.0);
+ EXPECT_DOUBLE_EQ(8, box.getVolume());
}
-TEST_F(imstkOrientedBoxTest, GetFunctionValue)
+TEST(imstkOrientedBoxTest, GetFunctionValue)
{
- m_box.setExtents(1.0, 1.0, 2.0);
- m_box.updatePostTransformData();
+ OrientedBox box;
+ box.setExtents(1.0, 1.0, 2.0);
+ box.updatePostTransformData();
- EXPECT_DOUBLE_EQ(-1., m_box.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-0.5, m_box.getFunctionValue(Vec3d(0.5, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(0.0, m_box.getFunctionValue(Vec3d(1.0, 1.0, 2.0)));
- EXPECT_DOUBLE_EQ(9.0, m_box.getFunctionValue(Vec3d(0.0, -10.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-1., box.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-0.5, box.getFunctionValue(Vec3d(0.5, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(0.0, box.getFunctionValue(Vec3d(1.0, 1.0, 2.0)));
+ EXPECT_DOUBLE_EQ(9.0, box.getFunctionValue(Vec3d(0.0, -10.0, 0.0)));
- m_box.rotate(Vec3d(1.0, 0.0, 0.0), 0.5 * PI);
- m_box.updatePostTransformData();
+ box.rotate(Vec3d(1.0, 0.0, 0.0), 0.5 * PI);
+ box.updatePostTransformData();
- EXPECT_DOUBLE_EQ(-1.0, m_box.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-0.5, m_box.getFunctionValue(Vec3d(0.5, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-0.5, m_box.getFunctionValue(Vec3d(0.0, -1.5, 0.0)));
+ EXPECT_DOUBLE_EQ(-1.0, box.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-0.5, box.getFunctionValue(Vec3d(0.5, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-0.5, box.getFunctionValue(Vec3d(0.0, -1.5, 0.0)));
}
diff --git a/Source/Geometry/Testing/imstkPlaneTest.cpp b/Source/Geometry/Testing/imstkPlaneTest.cpp
index 83dc2c0a99afa2e580e4bc0a66b7a0384549c99d..7877e807ab58e8414b76b79d2f22b1828692f2ab 100644
--- a/Source/Geometry/Testing/imstkPlaneTest.cpp
+++ b/Source/Geometry/Testing/imstkPlaneTest.cpp
@@ -19,76 +19,74 @@
=========================================================================*/
-#include "gtest/gtest.h"
-
#include "imstkPlane.h"
-using namespace imstk;
+#include <gtest/gtest.h>
-class imstkPlaneTest : public ::testing::Test
-{
-protected:
- Plane m_plane;
-};
+using namespace imstk;
-TEST_F(imstkPlaneTest, SetGetWidth)
+TEST(imstkPlaneTest, SetGetWidth)
{
- m_plane.setWidth(2);
- EXPECT_DOUBLE_EQ(2, m_plane.getWidth());
+ Plane plane;
+ plane.setWidth(2);
+ EXPECT_DOUBLE_EQ(2, plane.getWidth());
- m_plane.setWidth(0.003);
- EXPECT_DOUBLE_EQ(0.003, m_plane.getWidth());
+ plane.setWidth(0.003);
+ EXPECT_DOUBLE_EQ(0.003, plane.getWidth());
- m_plane.setWidth(400000000);
- EXPECT_DOUBLE_EQ(400000000.0, m_plane.getWidth());
+ plane.setWidth(400000000);
+ EXPECT_DOUBLE_EQ(400000000.0, plane.getWidth());
- m_plane.setWidth(0.0);
- EXPECT_DOUBLE_EQ(0.0, m_plane.getWidth());
+ plane.setWidth(0.0);
+ EXPECT_DOUBLE_EQ(0.0, plane.getWidth());
- m_plane.setWidth(-5.0);
- EXPECT_LT(-5.0, m_plane.getWidth());
+ plane.setWidth(-5.0);
+ EXPECT_LT(-5.0, plane.getWidth());
}
-TEST_F(imstkPlaneTest, SetGetNormal)
+TEST(imstkPlaneTest, SetGetNormal)
{
+ Plane plane;
Vec3d n1 = Vec3d(0.2, -0.3, 0.9);
Vec3d n2 = Vec3d(0.003, -0.001, 0.002);
Vec3d n3 = Vec3d(400000000, -500000000, 600000000);
- m_plane.setNormal(n1);
- EXPECT_TRUE(m_plane.getNormal().isApprox(n1.normalized()));
+ plane.setNormal(n1);
+ EXPECT_TRUE(plane.getNormal().isApprox(n1.normalized()));
- m_plane.setNormal(n2);
- EXPECT_TRUE(m_plane.getNormal().isApprox(n2.normalized()));
+ plane.setNormal(n2);
+ EXPECT_TRUE(plane.getNormal().isApprox(n2.normalized()));
- m_plane.setNormal(n3);
- EXPECT_TRUE(m_plane.getNormal().isApprox(n3.normalized()));
+ plane.setNormal(n3);
+ EXPECT_TRUE(plane.getNormal().isApprox(n3.normalized()));
- m_plane.setNormal(Vec3d::Zero());
- EXPECT_FALSE(m_plane.getNormal().isApprox(Vec3d::Zero()));
+ plane.setNormal(Vec3d::Zero());
+ EXPECT_FALSE(plane.getNormal().isApprox(Vec3d::Zero()));
}
-TEST_F(imstkPlaneTest, GetVolume)
+TEST(imstkPlaneTest, GetVolume)
{
- EXPECT_DOUBLE_EQ(0, m_plane.getVolume());
+ Plane plane;
+ EXPECT_DOUBLE_EQ(0, plane.getVolume());
}
-TEST_F(imstkPlaneTest, GetFunctionValue)
+TEST(imstkPlaneTest, GetFunctionValue)
{
- m_plane.setNormal(Vec3d(0., 1., 0.));
- m_plane.updatePostTransformData();
-
- EXPECT_DOUBLE_EQ(0., m_plane.getFunctionValue(Vec3d(0., 0., 0.)));
- EXPECT_DOUBLE_EQ(0., m_plane.getFunctionValue(Vec3d(0.5, 0., 0.)));
- EXPECT_DOUBLE_EQ(1., m_plane.getFunctionValue(Vec3d(1., 1., 1.)));
- EXPECT_DOUBLE_EQ(-10., m_plane.getFunctionValue(Vec3d(0., -10., 0.)));
-
- m_plane.setPosition(Vec3d(1., 1., 1.));
- m_plane.setNormal(Vec3d(1., 1., 1.));
- m_plane.updatePostTransformData();
-
- EXPECT_DOUBLE_EQ(-std::sqrt(3.0), m_plane.getFunctionValue(Vec3d(0., 0., 0.)));
- EXPECT_DOUBLE_EQ(0.0, m_plane.getFunctionValue(Vec3d(1., 1., 1.)));
- EXPECT_DOUBLE_EQ(-2.0 / std::sqrt(3.0), m_plane.getFunctionValue(Vec3d(1., 0., 0.)));
- EXPECT_DOUBLE_EQ(-13.0 / std::sqrt(3.0), m_plane.getFunctionValue(Vec3d(0., -10., 0.)));
+ Plane plane;
+ plane.setNormal(Vec3d(0., 1., 0.));
+ plane.updatePostTransformData();
+
+ EXPECT_DOUBLE_EQ(0., plane.getFunctionValue(Vec3d(0., 0., 0.)));
+ EXPECT_DOUBLE_EQ(0., plane.getFunctionValue(Vec3d(0.5, 0., 0.)));
+ EXPECT_DOUBLE_EQ(1., plane.getFunctionValue(Vec3d(1., 1., 1.)));
+ EXPECT_DOUBLE_EQ(-10., plane.getFunctionValue(Vec3d(0., -10., 0.)));
+
+ plane.setPosition(Vec3d(1., 1., 1.));
+ plane.setNormal(Vec3d(1., 1., 1.));
+ plane.updatePostTransformData();
+
+ EXPECT_DOUBLE_EQ(-std::sqrt(3.0), plane.getFunctionValue(Vec3d(0., 0., 0.)));
+ EXPECT_DOUBLE_EQ(0.0, plane.getFunctionValue(Vec3d(1., 1., 1.)));
+ EXPECT_DOUBLE_EQ(-2.0 / std::sqrt(3.0), plane.getFunctionValue(Vec3d(1., 0., 0.)));
+ EXPECT_DOUBLE_EQ(-13.0 / std::sqrt(3.0), plane.getFunctionValue(Vec3d(0., -10., 0.)));
}
diff --git a/Source/Geometry/Testing/imstkPointSetTest.cpp b/Source/Geometry/Testing/imstkPointSetTest.cpp
index 7ec6049b94a939455b5aa982097ce285a0bf925c..7a015902564f59311e415ad6de8c48a0c04e013a 100644
--- a/Source/Geometry/Testing/imstkPointSetTest.cpp
+++ b/Source/Geometry/Testing/imstkPointSetTest.cpp
@@ -19,12 +19,11 @@
=========================================================================*/
-#include "gtest/gtest.h"
-#include "gmock/gmock.h"
-
#include "imstkPointSet.h"
#include "imstkVecDataArray.h"
-#include "imstkMath.h"
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
using namespace imstk;
diff --git a/Source/Geometry/Testing/imstkSphereTest.cpp b/Source/Geometry/Testing/imstkSphereTest.cpp
index faf0628091f659d504e202e29b4a6ce43d24c689..a8dd60f2b14bfe1ac798927d948e4fe932dbef95 100644
--- a/Source/Geometry/Testing/imstkSphereTest.cpp
+++ b/Source/Geometry/Testing/imstkSphereTest.cpp
@@ -19,66 +19,63 @@
=========================================================================*/
-#include "gtest/gtest.h"
-
#include "imstkSphere.h"
-using namespace imstk;
+#include <gtest/gtest.h>
-class imstkSphereTest : public ::testing::Test
-{
-protected:
- Sphere m_sphere;
-};
+using namespace imstk;
-TEST_F(imstkSphereTest, SetGetRadius)
+TEST(imstkSphereTest, SetGetRadius)
{
- m_sphere.setRadius(2);
- EXPECT_DOUBLE_EQ(2, m_sphere.getRadius());
+ Sphere sphere;
+ sphere.setRadius(2);
+ EXPECT_DOUBLE_EQ(2, sphere.getRadius());
- m_sphere.setRadius(0.003);
- EXPECT_DOUBLE_EQ(0.003, m_sphere.getRadius());
+ sphere.setRadius(0.003);
+ EXPECT_DOUBLE_EQ(0.003, sphere.getRadius());
- m_sphere.setRadius(400000000);
- EXPECT_DOUBLE_EQ(400000000, m_sphere.getRadius());
+ sphere.setRadius(400000000);
+ EXPECT_DOUBLE_EQ(400000000, sphere.getRadius());
- m_sphere.setRadius(0);
- EXPECT_LT(0, m_sphere.getRadius());
+ sphere.setRadius(0);
+ EXPECT_LT(0, sphere.getRadius());
- m_sphere.setRadius(-5);
- EXPECT_LT(0, m_sphere.getRadius());
+ sphere.setRadius(-5);
+ EXPECT_LT(0, sphere.getRadius());
}
-TEST_F(imstkSphereTest, GetVolume)
+TEST(imstkSphereTest, GetVolume)
{
- m_sphere.setRadius(2);
- EXPECT_DOUBLE_EQ(4.0 / 3.0 * 8 * PI, m_sphere.getVolume());
+ Sphere sphere;
+ sphere.setRadius(2);
+ EXPECT_DOUBLE_EQ(4.0 / 3.0 * 8 * PI, sphere.getVolume());
- m_sphere.setRadius(0.003);
- EXPECT_DOUBLE_EQ(4.0 / 3.0 * PI * 0.003 * 0.003 * 0.003, m_sphere.getVolume());
+ sphere.setRadius(0.003);
+ EXPECT_DOUBLE_EQ(4.0 / 3.0 * PI * 0.003 * 0.003 * 0.003, sphere.getVolume());
double r = 400000000;
- m_sphere.setRadius(400000000);
- EXPECT_DOUBLE_EQ(4.0 / 3.0 * PI * r * r * r, m_sphere.getVolume());
+ sphere.setRadius(400000000);
+ EXPECT_DOUBLE_EQ(4.0 / 3.0 * PI * r * r * r, sphere.getVolume());
}
-TEST_F(imstkSphereTest, GetFunctionValue)
+TEST(imstkSphereTest, GetFunctionValue)
{
- m_sphere.setRadius(20.);
- m_sphere.updatePostTransformData();
-
- EXPECT_DOUBLE_EQ(-20.0, m_sphere.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-15.0, m_sphere.getFunctionValue(Vec3d(5.0, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-20.0 + std::sqrt(3), m_sphere.getFunctionValue(Vec3d(1.0, 1.0, 1.0)));
- EXPECT_DOUBLE_EQ(0.0, m_sphere.getFunctionValue(Vec3d(0.0, 20.0, 0.0)));
- EXPECT_DOUBLE_EQ(30.0, m_sphere.getFunctionValue(Vec3d(0.0, 0.0, 50.0)));
-
- m_sphere.rotate(Vec3d(1.0, 1.0, 0.0), 0.1 * PI);
- m_sphere.updatePostTransformData();
-
- EXPECT_DOUBLE_EQ(-20.0, m_sphere.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-15.0, m_sphere.getFunctionValue(Vec3d(5.0, 0.0, 0.0)));
- EXPECT_DOUBLE_EQ(-20.0 + std::sqrt(3), m_sphere.getFunctionValue(Vec3d(1.0, 1.0, 1.0)));
- EXPECT_NEAR(0., m_sphere.getFunctionValue(Vec3d(0.0, 20., 0.0)), 1.0e-10);
- EXPECT_DOUBLE_EQ(30.0, m_sphere.getFunctionValue(Vec3d(0.0, 0.0, 50.0)));
+ Sphere sphere;
+ sphere.setRadius(20.);
+ sphere.updatePostTransformData();
+
+ EXPECT_DOUBLE_EQ(-20.0, sphere.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-15.0, sphere.getFunctionValue(Vec3d(5.0, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-20.0 + std::sqrt(3), sphere.getFunctionValue(Vec3d(1.0, 1.0, 1.0)));
+ EXPECT_DOUBLE_EQ(0.0, sphere.getFunctionValue(Vec3d(0.0, 20.0, 0.0)));
+ EXPECT_DOUBLE_EQ(30.0, sphere.getFunctionValue(Vec3d(0.0, 0.0, 50.0)));
+
+ sphere.rotate(Vec3d(1.0, 1.0, 0.0), 0.1 * PI);
+ sphere.updatePostTransformData();
+
+ EXPECT_DOUBLE_EQ(-20.0, sphere.getFunctionValue(Vec3d(0.0, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-15.0, sphere.getFunctionValue(Vec3d(5.0, 0.0, 0.0)));
+ EXPECT_DOUBLE_EQ(-20.0 + std::sqrt(3), sphere.getFunctionValue(Vec3d(1.0, 1.0, 1.0)));
+ EXPECT_NEAR(0., sphere.getFunctionValue(Vec3d(0.0, 20., 0.0)), 1.0e-10);
+ EXPECT_DOUBLE_EQ(30.0, sphere.getFunctionValue(Vec3d(0.0, 0.0, 50.0)));
}
diff --git a/Source/Geometry/Testing/imstkSurfaceMeshTest.cpp b/Source/Geometry/Testing/imstkSurfaceMeshTest.cpp
index f5dbf97a01134356c25b07756104f81dcf2a9887..28a852998fbd364e6813ae3d060dbc8fe2d6c6e0 100644
--- a/Source/Geometry/Testing/imstkSurfaceMeshTest.cpp
+++ b/Source/Geometry/Testing/imstkSurfaceMeshTest.cpp
@@ -19,15 +19,12 @@
=========================================================================*/
-#include "imstkOrientedBox.h"
#include "imstkGeometryUtilities.h"
-#include "imstkMath.h"
-#include "imstkPointSet.h"
+#include "imstkOrientedBox.h"
#include "imstkSurfaceMesh.h"
-#include "imstkVecDataArray.h"
-#include <gtest/gtest.h>
#include <gmock/gmock.h>
+
using namespace imstk;
namespace
diff --git a/Source/Geometry/Testing/imstkTetrahedralMeshTest.cpp b/Source/Geometry/Testing/imstkTetrahedralMeshTest.cpp
index 33aa1fea0f3b049cdb5b42b554c7a8fdbaeac54c..941d7e9ffe1a2168b6bfeda8a1aafaf09ffc6d37 100644
--- a/Source/Geometry/Testing/imstkTetrahedralMeshTest.cpp
+++ b/Source/Geometry/Testing/imstkTetrahedralMeshTest.cpp
@@ -27,17 +27,13 @@
using namespace imstk;
-class imstkTetrahedralMeshTest : public ::testing::Test
-{
-protected:
- TetrahedralMesh m_tetMesh;
-};
-
///
/// \brief Test the mesh extraction of a cube
///
-TEST_F(imstkTetrahedralMeshTest, ExtractSurfaceMesh)
+TEST(imstkTetrahedralMeshTest, ExtractSurfaceMesh)
{
+ TetrahedralMesh tetMesh;
+
// Setup a cube
// 0-------1
// /| /|
@@ -68,11 +64,11 @@ TEST_F(imstkTetrahedralMeshTest, ExtractSurfaceMesh)
indices[3] = Vec4i(1, 2, 0, 5);
indices[4] = Vec4i(2, 6, 7, 5);
- m_tetMesh.initialize(verticesPtr, indicesPtr);
+ tetMesh.initialize(verticesPtr, indicesPtr);
}
// Extract the surface
- std::shared_ptr<SurfaceMesh> surfMesh = m_tetMesh.extractSurfaceMesh();
+ std::shared_ptr<SurfaceMesh> surfMesh = tetMesh.extractSurfaceMesh();
std::shared_ptr<VecDataArray<double, 3>> surfVerticesPtr = surfMesh->getVertexPositions();
std::shared_ptr<VecDataArray<int, 3>> surfIndicesPtr = surfMesh->getTriangleIndices();
VecDataArray<double, 3>& surfVertices = *surfVerticesPtr;
@@ -99,8 +95,10 @@ TEST_F(imstkTetrahedralMeshTest, ExtractSurfaceMesh)
///
/// \brief Test the computation of volume
///
-TEST_F(imstkTetrahedralMeshTest, GetVolume)
+TEST(imstkTetrahedralMeshTest, GetVolume)
{
+ TetrahedralMesh tetMesh;
+
// We use a regular tetrahedron with edge lengths 2
// V = (edge length)^3/(6sqrt(2))
const double edgeLenth = 2.0;
@@ -118,6 +116,6 @@ TEST_F(imstkTetrahedralMeshTest, GetVolume)
indices[0] = Vec4i(0, 1, 2, 3);
- m_tetMesh.initialize(verticesPtr, indicesPtr);
- EXPECT_NEAR(expectedVolume, m_tetMesh.getVolume(), 0.000001);
+ tetMesh.initialize(verticesPtr, indicesPtr);
+ EXPECT_NEAR(expectedVolume, tetMesh.getVolume(), 0.000001);
}
diff --git a/Source/MeshIO/CMakeLists.txt b/Source/MeshIO/CMakeLists.txt
index df9f080ac2ec07c2fbc2fd6d2e90b22481e29970..7efb068faacc158d7af23aab89196ce4c5cf473a 100644
--- a/Source/MeshIO/CMakeLists.txt
+++ b/Source/MeshIO/CMakeLists.txt
@@ -14,7 +14,6 @@ imstk_add_library( MeshIO
#-----------------------------------------------------------------------------
# Testing
#-----------------------------------------------------------------------------
-#if( ${PROJECT_NAME}_BUILD_TESTING )
-# include(imstkAddTest)
-# imstk_add_test( Geometry )
-#endif()
+if( ${PROJECT_NAME}_BUILD_TESTING )
+ add_subdirectory("Testing")
+endif()
diff --git a/Source/MeshIO/Testing/CMakeLists.txt b/Source/MeshIO/Testing/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..afe9cf1acfe802844f7ab320937cb5d5330b88de
--- /dev/null
+++ b/Source/MeshIO/Testing/CMakeLists.txt
@@ -0,0 +1,21 @@
+###########################################################################
+#
+# 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.
+#
+###########################################################################
+
+include(imstkAddTest)
+imstk_add_test( MeshIO )
+target_link_libraries( MeshIOTests Common Geometry )
\ No newline at end of file
diff --git a/Source/MeshIO/Testing/imstkMeshIOTest.cpp b/Source/MeshIO/Testing/imstkMeshIOTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7a61534b09fb1fe743e2f2a625b6209ad6cb2bbb
--- /dev/null
+++ b/Source/MeshIO/Testing/imstkMeshIOTest.cpp
@@ -0,0 +1,33 @@
+/*=========================================================================
+
+ 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 "gtest/gtest.h"
+
+#include "imstkMeshIO.h"
+#include "imstkGeometryUtilities.h"
+#include "imstkSurfaceMesh.h"
+
+using namespace imstk;
+
+TEST(imstkMeshIODeathTest, FailOnMissingFile)
+{
+ EXPECT_DEATH(MeshIO::read<SurfaceMesh>(iMSTK_DATA_ROOT "doesntexist.obj"), "doesntexist.obj doesn't exist");
+}
\ No newline at end of file
diff --git a/Source/MeshIO/imstkMeshIO.cpp b/Source/MeshIO/imstkMeshIO.cpp
index 4c230f768bdf264ef272d570f006d95f9c92dfee..e3966184ac0ca5cdeced48f646b645a8132ba4aa 100644
--- a/Source/MeshIO/imstkMeshIO.cpp
+++ b/Source/MeshIO/imstkMeshIO.cpp
@@ -60,13 +60,10 @@ static std::unordered_map<std::string, MeshFileType> extToType =
std::shared_ptr<PointSet>
MeshIO::read(const std::string& filePath)
{
- bool isDir = false;
+ bool isDirectory = false;
+ bool exists = fileExists(filePath, isDirectory);
- if (isDir)
- {
- // Assume that the directory is a collection of DICOM files
- return VTKMeshIO::read(filePath, MeshFileType::DCM);
- }
+ CHECK(exists && !isDirectory) << "File " << filePath << " doesn't exist or is a directory.";
MeshFileType meshType = MeshIO::getFileType(filePath);
switch (meshType)
diff --git a/Source/Scene/imstkScene.h b/Source/Scene/imstkScene.h
index a5034093be10633b611465dbb619ead9b63a2948..96d99cc031bc948c2287675e11a407fd49a4f78f 100644
--- a/Source/Scene/imstkScene.h
+++ b/Source/Scene/imstkScene.h
@@ -56,7 +56,7 @@ struct SceneConfig
bool trackFPS = false;
// If off, tasks will run sequentially
- bool taskParallelizationEnabled = true;
+ bool taskParallelizationEnabled = false;
// If on, elapsed times for computational steps will be reported in map
bool taskTimingEnabled = false;
diff --git a/Source/Solvers/Testing/imstkPGSSolverTest.cpp b/Source/Solvers/Testing/imstkPGSSolverTest.cpp
index 12ef5a8bebbd68ae87f6f4092dcd95ef76d3cd33..2d6a2e48a62cb63570f69e51db5cc204fb5ff825 100644
--- a/Source/Solvers/Testing/imstkPGSSolverTest.cpp
+++ b/Source/Solvers/Testing/imstkPGSSolverTest.cpp
@@ -26,17 +26,13 @@
using namespace imstk;
-class imstkPGSSolverTest : public ::testing::Test
-{
-protected:
- ProjectedGaussSeidelSolver<double> m_solver;
-};
-
///
-/// \brief TODO
+/// \brief Tests PGS solving of a diagonal 5x5 matrix
///
-TEST_F(imstkPGSSolverTest, Solve5x5)
+TEST(imstkPGSSolverTest, Solve5x5)
{
+ ProjectedGaussSeidelSolver<double> solver;
+
// Testing Ax=b
Eigen::MatrixXd Ad(5, 5);
Ad <<
@@ -61,12 +57,12 @@ TEST_F(imstkPGSSolverTest, Solve5x5)
cu(i, 1) = IMSTK_DOUBLE_MAX;
}
- m_solver.setA(&A);
- m_solver.setMaxIterations(1000);
- m_solver.setRelaxation(0.05);
- m_solver.setEpsilon(1.0e-8);
+ solver.setA(&A);
+ solver.setMaxIterations(1000);
+ solver.setRelaxation(0.05);
+ solver.setEpsilon(1.0e-8);
- Eigen::VectorXd x = m_solver.solve(b, cu);
+ Eigen::VectorXd x = solver.solve(b, cu);
// Check that Ax now equals b
// Test this way in case multiple solutions exist, here we are only