diff --git a/Examples/DeformableBody/DeformableBodyExample.cpp b/Examples/DeformableBody/DeformableBodyExample.cpp index d1059c517abf2eeaf428fefd4dc865ee75a42e84..4d3de72b0f28f2ed826846db954ecc3061947bba 100644 --- a/Examples/DeformableBody/DeformableBodyExample.cpp +++ b/Examples/DeformableBody/DeformableBodyExample.cpp @@ -23,7 +23,7 @@ #include "imstkSimulationManager.h" #include "imstkDeformableObject.h" #include "imstkBackwardEuler.h" -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" #include "imstkNewtonSolver.h" #include "imstkGaussSeidel.h" #include "imstkPlane.h" diff --git a/Examples/NonlinearSolver/CMakeLists.txt b/Examples/NonlinearSolver/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a66dd3fc570afa0725ae5830d17ba193ce86e329 --- /dev/null +++ b/Examples/NonlinearSolver/CMakeLists.txt @@ -0,0 +1,35 @@ +########################################################################### +# +# Copyright (c) Kitware, Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +########################################################################### + +project(Example-NonlinearSolver) + +#----------------------------------------------------------------------------- +# Create executable +#----------------------------------------------------------------------------- +imstk_add_executable(${PROJECT_NAME} nonlinearSolver.cpp) + +#----------------------------------------------------------------------------- +# Add the target to Examples folder +#----------------------------------------------------------------------------- +SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples) + +#----------------------------------------------------------------------------- +# Link libraries to executable +#----------------------------------------------------------------------------- +target_link_libraries(${PROJECT_NAME} Solvers Common) + diff --git a/Examples/NonlinearSolver/nonlinearSolver.cpp b/Examples/NonlinearSolver/nonlinearSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1725f74c8351bfb07f27fed3fbb7c273434782d9 --- /dev/null +++ b/Examples/NonlinearSolver/nonlinearSolver.cpp @@ -0,0 +1,91 @@ +/*========================================================================= + + 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 <ios> +#include <iostream> +#include <vector> +#include <iomanip> + +#include "imstkMath.h" +#include "imstkNewtonSolver.h" +#include "imstkDirectLinearSolver.h" +#include "imstkNonLinearSolver.h" + +using namespace imstk; + +int +main(int argc, char** argv) +{ + const int N = 2; + auto x = Vectord(N); + auto xe = Vectord(N); + auto y = Vectord(N); + auto A = Matrixd(N, N); + + x[0] = 100.0; + x[1] = 100.0; + xe[0] = 1.0; + xe[1] = 10.0; + + auto func = [&y](const Vectord& x, const bool isSemiImplicit) -> const Vectord& { + // auto y = Vectord(x.size()); + y[0] = x[0] * x[0] - 1.0; + y[1] = x[1] * x[1] - 100.0; + + return y; + }; + + auto jac = [&A](const Vectord& x) -> const Matrixd& { + // auto A = Matrixd(x.size(), x.size()); + A(0, 0) = 2 * x[0]; + A(0, 1) = 0.0; + A(1, 0) = 0.0; + A(1, 1) = 2 * x[1]; + + return A; + }; + + auto updateX = [&x](const Vectord& du, const bool isSemiImplicit) + { + x -= du; + return; + }; + + auto updateXold = [](void) {}; + + auto nlSystem = std::make_shared<NonLinearSystem<Matrixd>>(func, jac); + nlSystem->setUnknownVector(x); + nlSystem->setUpdateFunction(updateX); + nlSystem->setUpdatePreviousStatesFunction(updateXold); + + auto linSolver = std::make_shared<DirectLinearSolver<Matrixd>>(); + auto nlSolver = std::make_shared<NewtonSolver<Matrixd>>(); + nlSolver->setMaxIterations(100); + nlSolver->setRelativeTolerance(1e-8); + nlSolver->setAbsoluteTolerance(1e-10); + nlSolver->setSystem(nlSystem); + nlSolver->setLinearSolver(linSolver); + + std::cout << "init_error = " << std::setprecision(12) << std::scientific << (x-xe).norm() << std::endl; + nlSolver->solve(); + + std::cout << "final_error = " << std::setprecision(12) << std::scientific << (x-xe).norm() << std::endl; +} diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h index ec59caf229d171d1dd34abb6b78c92e942e69090..9641e0cba7dfb8409c037f238f0619cd75070c82 100644 --- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h +++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h @@ -24,7 +24,7 @@ // imstk #include "imstkDynamicalModel.h" #include "imstkInternalForceModelTypes.h" -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" // vega #include "sparseMatrix.h" diff --git a/Source/Solvers/imstkDirectLinearSolver.cpp b/Source/Solvers/imstkDirectLinearSolver.cpp index b319082a5956ede2912d8458ade6425ab5f6613b..77b4018ee1db170783bed698d08968edaaf57dce 100644 --- a/Source/Solvers/imstkDirectLinearSolver.cpp +++ b/Source/Solvers/imstkDirectLinearSolver.cpp @@ -83,7 +83,7 @@ DirectLinearSolver<Matrixd>::solve(const Vectord& rhs, Vectord& x) void DirectLinearSolver<Matrixd>::solve(Vectord& x) { - if(!m_system_set) LOG(FATAL) << "Linear system has not been set"; + if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set"; x.setZero(); auto b = m_linearSystem->getRHSVector(); diff --git a/Source/Solvers/imstkDirectLinearSolver.h b/Source/Solvers/imstkDirectLinearSolver.h index 9f1a77c5be5f79627163818df61247229abeed84..e586d811f03ec1e8fb0b21dd0cdc961d23616017 100644 --- a/Source/Solvers/imstkDirectLinearSolver.h +++ b/Source/Solvers/imstkDirectLinearSolver.h @@ -85,7 +85,6 @@ public: private: Eigen::LDLT<Matrixd> m_solver; - bool m_system_set; }; /// @@ -124,6 +123,5 @@ public: private: Eigen::SparseLU<SparseMatrixd, Eigen::COLAMDOrdering<MatrixType::StorageIndex>> m_solver;//? - bool m_system_set; }; } // imstk diff --git a/Source/Solvers/imstkLinearSystem.h b/Source/Solvers/imstkLinearSystem.h index 22e1638ef1155b9325702211292d56ddda084c95..5c0814258a80dd0e654bbfcdaf44e4e519d12b3b 100644 --- a/Source/Solvers/imstkLinearSystem.h +++ b/Source/Solvers/imstkLinearSystem.h @@ -21,7 +21,7 @@ #pragma once -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" namespace imstk { diff --git a/Source/Solvers/imstkNewtonSolver.cpp b/Source/Solvers/imstkNewtonSolver.cpp index c0097d7e0f3a15f459286c6b46cad06d21600cea..7632722119b1635cff0efc2ee94d6dd36e8faf57 100644 --- a/Source/Solvers/imstkNewtonSolver.cpp +++ b/Source/Solvers/imstkNewtonSolver.cpp @@ -114,11 +114,11 @@ NewtonSolver<SystemMatrix>::solve() size_t iterNum; const auto& u = this->m_nonLinearSystem->getUnknownVector(); Vectord du = u; // make this a class member in future + double error0 = MAX_D; double epsilon = m_relativeTolerance * m_relativeTolerance; for (iterNum = 0; iterNum < m_maxIterations; ++iterNum) { - double error0 = MAX_D; double error = updateJacobian(u); if (iterNum == 0) @@ -128,13 +128,14 @@ NewtonSolver<SystemMatrix>::solve() if (error / error0 < epsilon && iterNum > 0) { - //std::cout << "Num. of Newton Iterations: " << i << "\tError ratio: " << error/error0 << std::endl; + // std::cout << "Num. of Newton Iterations: " << iterNum << "\tError ratio: " << error/error0 << ", " << error << " " << error0 << std::endl; break; } m_linearSolver->solve(du); this->m_nonLinearSystem->m_FUpdate(du, this->m_isSemiImplicit); } + this->m_nonLinearSystem->m_FUpdatePrevState(); if (iterNum == m_maxIterations && !this->m_isSemiImplicit) @@ -167,7 +168,7 @@ NewtonSolver<SystemMatrix>::updateJacobian(const Vectord& x) //linearSystem->setLinearProjectors(this->m_nonLinearSystem->getLinearProjectors()); /// \todo Left for near future reference. Clear in future. m_linearSolver->setSystem(linearSystem); - return b.dot(b); + return std::sqrt(b.dot(b)); } template <typename SystemMatrix> diff --git a/Source/Solvers/imstkNonLinearSolver.h b/Source/Solvers/imstkNonLinearSolver.h index adff8d9de9bdbfbab67ae81a5adf908fe165a6f8..3902fb310a54ca81f1e11733ed556ce8efcc15de 100644 --- a/Source/Solvers/imstkNonLinearSolver.h +++ b/Source/Solvers/imstkNonLinearSolver.h @@ -26,7 +26,7 @@ #include <functional> // imstk includes -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" #include "imstkMath.h" #include "imstkSolverBase.h" diff --git a/Source/Solvers/imstkNonlinearSystem.cpp b/Source/Solvers/imstkNonLinearSystem.cpp similarity index 98% rename from Source/Solvers/imstkNonlinearSystem.cpp rename to Source/Solvers/imstkNonLinearSystem.cpp index a014ea8efc42dc89bbe432fe1290a1b85a942b5e..92c688c54b8a1ac4b4e0409814c7fdfb58b6bc80 100644 --- a/Source/Solvers/imstkNonlinearSystem.cpp +++ b/Source/Solvers/imstkNonLinearSystem.cpp @@ -19,7 +19,7 @@ =========================================================================*/ -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" #include "imstkMath.h" // #include <g3log/g3log.hpp> diff --git a/Source/Solvers/imstkNonlinearSystem.h b/Source/Solvers/imstkNonLinearSystem.h similarity index 100% rename from Source/Solvers/imstkNonlinearSystem.h rename to Source/Solvers/imstkNonLinearSystem.h diff --git a/Source/Solvers/imstkSOR.h b/Source/Solvers/imstkSOR.h index f310b48bc70d1dde34f439fdf74cdf76864ade3d..d71a93fd9af8aeea16b06a432ab856f6723d3c0e 100644 --- a/Source/Solvers/imstkSOR.h +++ b/Source/Solvers/imstkSOR.h @@ -22,7 +22,7 @@ #pragma once // iMSTK includes -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" #include "imstkIterativeLinearSolver.h" // Eigen includes diff --git a/Source/apiUtilities/imstkAPIUtilities.cpp b/Source/apiUtilities/imstkAPIUtilities.cpp index 59b3a71e74fcb4df88afe3dbf21e1163eb12221c..82bfcb51e895fd2093cfd7cf062fa32cf12ff387 100644 --- a/Source/apiUtilities/imstkAPIUtilities.cpp +++ b/Source/apiUtilities/imstkAPIUtilities.cpp @@ -30,7 +30,7 @@ #include "imstkColor.h" // Solvers -#include "imstkNonlinearSystem.h" +#include "imstkNonLinearSystem.h" // Geometry #include "imstkPlane.h"