diff --git a/include/GooseFEM/Iterate.h b/include/GooseFEM/Iterate.h index 04fc28c..f2e97ff 100644 --- a/include/GooseFEM/Iterate.h +++ b/include/GooseFEM/Iterate.h @@ -1,66 +1,66 @@ /** Support function for iterations. \file Iterate.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ITERATE_H #define GOOSEFEM_ITERATE_H #include "config.h" namespace GooseFEM { /** Support function for iterations in end-user programs. */ namespace Iterate { /** Class to perform a residual check based on the last "n" iterations. A typical usage is in dynamic simulations where equilibrium is checked based on a force residual. Fluctuations could however be responsible for this criterion to be triggered too early. By checking several time-steps such case can be avoided. */ class StopList { public: /** Constructor. \param n Number of consecutive iterations to consider. */ StopList(size_t n = 1); /** Reset all residuals to infinity. */ void reset(); /** Reset all residuals to infinity, and change the number of residuals to check. \param n Number of consecutive iterations to consider. */ void reset(size_t n); /** - Update list of residuals, return true if all residuals are below the tolerance. + Update list of residuals, return `true` if all residuals are below the tolerance. \param res Current residual. \param tol Tolerance below which all last "n" iterations must lie. */ bool stop(double res, double tol); private: std::vector m_res; ///< List with residuals. }; } // namespace Iterate } // namespace GooseFEM #include "Iterate.hpp" #endif diff --git a/include/GooseFEM/Iterate.hpp b/include/GooseFEM/Iterate.hpp index 5505259..11c9e11 100644 --- a/include/GooseFEM/Iterate.hpp +++ b/include/GooseFEM/Iterate.hpp @@ -1,67 +1,47 @@ /** Implementation of Iterate.h \file Iterate.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_ITERATE_HPP #define GOOSEFEM_ITERATE_HPP #include "Iterate.h" namespace GooseFEM { namespace Iterate { inline StopList::StopList(size_t n) { m_res.resize(n); reset(); } inline void StopList::reset() { - for (auto& i : m_res) { - i = std::numeric_limits::infinity(); - } + std::fill(m_res.begin(), m_res.end(), std::numeric_limits::infinity()); } inline void StopList::reset(size_t n) { m_res.resize(n); reset(); } inline bool StopList::stop(double res, double tol) { - // move residual one place back (forgetting the first) - for (size_t i = 1; i < m_res.size(); ++i) { - m_res[i - 1] = m_res[i]; - } - - // add new residual to the end - m_res[m_res.size() - 1] = res; - - // check for convergence: all residuals should be below the tolerance - for (size_t i = 0; i < m_res.size(); ++i) { - if (m_res[i] > tol) { - return false; - } + std::rotate(m_res.begin(), m_res.begin() + 1, m_res.end()); + m_res.back() = res; + if (res > tol) { + return false; } - - // check for convergence: all residuals should be decreasing - for (size_t i = 1; i < m_res.size(); ++i) { - if (m_res[i] > m_res[i - 1]) { - return false; - } - } - - // all checks passed: signal convergence - return true; + return std::is_sorted(m_res.begin(), m_res.end(), std::greater_equal()) && m_res.front() <= tol; } } // namespace Iterate } // namespace GooseFEM #endif diff --git a/python/Iterate.hpp b/python/Iterate.hpp new file mode 100644 index 0000000..c886c22 --- /dev/null +++ b/python/Iterate.hpp @@ -0,0 +1,39 @@ +/** +\file +\copyright Copyright 2017. Tom de Geus. All rights reserved. +\license This project is released under the GNU Public License (GPLv3). +*/ + +#ifndef PYGOOSEFEM_ITERATE_H +#define PYGOOSEFEM_ITERATE_H + +#include +#include + +namespace py = pybind11; + +void init_Iterate(py::module& mod) +{ + py::class_ cls(mod, "Iterate"); + + cls.def(py::init(), + "Class to perform a residual check based on the last 'n' iterations." + "See :cpp:class:`GooseFEM::Iterate::StopList`.", + py::arg("n") = 1); + + cls.def("reset", + py::overload_cast<>(&GooseFEM::Iterate::StopList::reset), + "Reset." + "See :cpp:func:`GooseFEM::Iterate::StopList::reset`."); + + cls.def("stop", + &GooseFEM::Iterate::StopList::stop, + "Update list of residuals, return `true` if all residuals are below the tolerance." + "See :cpp:func:`GooseFEM::Iterate::StopList::stop`.", + py::arg("res"), + py::arg("tol")); + + cls.def("__repr__", [](const GooseFEM::Iterate::StopList&) { return ""; }); +} + +#endif diff --git a/python/main.cpp b/python/main.cpp index 20f9473..eb8befd 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -1,137 +1,147 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #define FORCE_IMPORT_ARRAY #include #include #include #include #define GOOSEFEM_ENABLE_ASSERT #define GOOSEFEM_ENABLE_WARNING_PYTHON namespace py = pybind11; #include "version.hpp" #include "Allocate.hpp" #include "Vector.hpp" #include "VectorPartitioned.hpp" #include "VectorPartitionedTyings.hpp" #include "Matrix.hpp" #include "MatrixPartitioned.hpp" #include "MatrixPartitionedTyings.hpp" #include "MatrixDiagonal.hpp" #include "MatrixDiagonalPartitioned.hpp" #include "Element.hpp" #include "ElementQuad4.hpp" #include "ElementQuad4Planar.hpp" #include "ElementQuad4Axisymmetric.hpp" #include "ElementHex8.hpp" #include "Mesh.hpp" #include "MeshTri3.hpp" #include "MeshQuad4.hpp" #include "MeshHex8.hpp" +#include "Iterate.hpp" PYBIND11_MODULE(GooseFEM, m) { xt::import_numpy(); // -------- // GooseFEM // -------- m.doc() = "Some simple finite element meshes and operations"; init_version(m); init_Allocate(m); init_Vector(m); init_VectorPartitioned(m); init_VectorPartitionedTyings(m); init_Matrix(m); init_MatrixPartitioned(m); init_MatrixPartitionedTyings(m); init_MatrixDiagonal(m); init_MatrixDiagonalPartitioned(m); +// ---------------- +// GooseFEM.Iterate +// ---------------- + +py::module mIterate = m.def_submodule("Iterate", "Iteration support tools"); + +init_Iterate(mIterate); + + // ---------------- // GooseFEM.Element // ---------------- py::module mElement = m.def_submodule("Element", "Generic element routines"); init_Element(mElement); // ---------------------- // GooseFEM.Element.Quad4 // ---------------------- py::module mElementQuad4 = mElement.def_submodule("Quad4", "Linear quadrilateral elements (2D)"); py::module mElementQuad4Gauss = mElementQuad4.def_submodule("Gauss", "Gauss quadrature"); py::module mElementQuad4Nodal = mElementQuad4.def_submodule("Nodal", "Nodal quadrature"); py::module mElementQuad4MidPoint = mElementQuad4.def_submodule("MidPoint", "MidPoint quadrature"); init_ElementQuad4(mElementQuad4); init_ElementQuad4Planar(mElementQuad4); init_ElementQuad4Axisymmetric(mElementQuad4); init_ElementQuad4Gauss(mElementQuad4Gauss); init_ElementQuad4Nodal(mElementQuad4Nodal); init_ElementQuad4MidPoint(mElementQuad4MidPoint); // --------------------- // GooseFEM.Element.Hex8 // --------------------- py::module mElementHex8 = mElement.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)"); py::module mElementHex8Gauss = mElementHex8.def_submodule("Gauss", "Gauss quadrature"); py::module mElementHex8Nodal = mElementHex8.def_submodule("Nodal", "Nodal quadrature"); init_ElementHex8(mElementHex8); init_ElementHex8Gauss(mElementHex8Gauss); init_ElementHex8Nodal(mElementHex8Nodal); // ------------- // GooseFEM.Mesh // ------------- py::module mMesh = m.def_submodule("Mesh", "Generic mesh routines"); init_Mesh(mMesh); // ------------------ // GooseFEM.Mesh.Tri3 // ------------------ py::module mMeshTri3 = mMesh.def_submodule("Tri3", "Linear triangular elements (2D)"); init_MeshTri3(mMeshTri3); // ------------------- // GooseFEM.Mesh.Quad4 // ------------------- py::module mMeshQuad4 = mMesh.def_submodule("Quad4", "Linear quadrilateral elements (2D)"); init_MeshQuad4(mMeshQuad4); py::module mMeshQuad4Map = mMeshQuad4.def_submodule("Map", "Map mesh objects"); init_MeshQuad4Map(mMeshQuad4Map); // ------------------ // GooseFEM.Mesh.Hex8 // ------------------ py::module mMeshHex8 = mMesh.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)"); init_MeshHex8(mMeshHex8); } diff --git a/test/basic/Iterate.cpp b/test/basic/Iterate.cpp index aa57071..db6f797 100644 --- a/test/basic/Iterate.cpp +++ b/test/basic/Iterate.cpp @@ -1,28 +1,38 @@ #include #include #include #include #include #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12)); TEST_CASE("GooseFEM::Iterate", "Iterate.h") { - SECTION("StopList") + SECTION("StopList - sorted input") { GooseFEM::Iterate::StopList stop(5); - REQUIRE(stop.stop(5.e+0, 1.e-3) == false); - REQUIRE(stop.stop(5.e+1, 1.e-3) == false); - REQUIRE(stop.stop(5.e-1, 1.e-3) == false); - REQUIRE(stop.stop(5.e-2, 1.e-3) == false); - REQUIRE(stop.stop(5.e-3, 1.e-3) == false); - REQUIRE(stop.stop(5.e-4, 1.e-3) == false); - REQUIRE(stop.stop(5.e-4, 1.e-3) == false); - REQUIRE(stop.stop(5.e-4, 1.e-3) == false); - REQUIRE(stop.stop(5.e-4, 1.e-3) == false); - REQUIRE(stop.stop(5.e-4, 1.e-3) == true); + // x x x x x v v v v v + std::vector res = {5e+0, 5e+1, 5e-1, 5e-2, 5e-3, 5e-4, 4e-4, 3e-4, 2e-4, 1e-4}; + std::vector conv = {false, false, false, false, false, false, false, false, false, true}; + + for (size_t i = 0; i < res.size(); ++i) { + REQUIRE(stop.stop(res[i], 1e-3) == conv[i]); + } + } + + SECTION("StopList - unsorted input") + { + GooseFEM::Iterate::StopList stop(5); + + // x x x x x v v v v x v v v v + std::vector res = {5e+0, 5e+1, 5e-1, 5e-2, 5e-3, 5e-4, 4e-4, 3e-4, 2e-4, 3e-4, 2e-4, 1e-4, 9e-5, 8e-5}; + std::vector conv = {false, false, false, false, false, false, false, false, false, false, false, false, false, true}; + + for (size_t i = 0; i < res.size(); ++i) { + REQUIRE(stop.stop(res[i], 1e-3) == conv[i]); + } } }