diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e68454..ec08fce 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.0) +cmake_minimum_required(VERSION 2.8) project(MySources) set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) @@ -9,30 +9,42 @@ if (CMAKE_VERSION VERSION_LESS 3.2) else() set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1") endif() +<<<<<<< HEAD #set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies") if(NOT EXISTS EXTERNAL_DEPS_DIR) +======= +set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies") +if(NOT EXISTS ${EXTERNAL_DEPS_DIR}) +>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR}) message("External dependencies directory set:" ${EXTERNAL_DEPS_DIR}) endif() ##Create directory for the external libraries file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) +<<<<<<< HEAD #message(${POLYSCOPE_ALREADY_COMPILED}) if(NOT DEFINED USE_POLYSCOPE) message(FATAL_ERROR "Use polyscope was not defined") endif() +======= +>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb if(${USE_POLYSCOPE}) + set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope) download_project(PROJ POLYSCOPE GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git GIT_TAG master PREFIX ${EXTERNAL_DEPS_DIR} + BINARY_DIR ${POLYSCOPE_BINARY_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE} ) add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR}) + add_compile_definitions(POLYSCOPE_DEFINED) endif() +<<<<<<< HEAD ##dlib #set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) @@ -45,6 +57,8 @@ endif() # ${UPDATE_DISCONNECTED_IF_AVAILABLE} #) #add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR}) +======= +>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb ##vcglib devel branch download_project(PROJ vcglib_devel GIT_REPOSITORY https://github.com/IasonManolas/vcglib.git @@ -63,53 +77,80 @@ download_project(PROJ EIGEN ) add_subdirectory(${EIGEN_SOURCE_DIR} ${EIGEN_BINARY_DIR}) +##matplot++ lib +set(MATPLOTPLUSPLUS_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/matplot) +download_project(PROJ MATPLOTPLUSPLUS + GIT_REPOSITORY https://github.com/alandefreitas/matplotplusplus + GIT_TAG master + BINARY_DIR ${MATPLOTPLUSPLUS_BINARY_DIR} + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} + ) +add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR}) + ##threed-beam-fea +set(threed-beam-fea_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/threed-beam-fea) download_project(PROJ threed-beam-fea GIT_REPOSITORY https://gitea-s2i2s.isti.cnr.it/manolas/threed-beam-fea.git GIT_TAG master + BINARY_DIR ${threed-beam-fea_BINARY_DIR} PREFIX ${EXTERNAL_DEPS_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE} ) add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR}) +<<<<<<< HEAD +======= +##TBB +set(TBB_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/tbb) +download_project(PROJ TBB + GIT_REPOSITORY https://github.com/wjakob/tbb.git + GIT_TAG master + PREFIX ${EXTERNAL_DEPS_DIR} + BINARY_DIR ${TBB_BINARY_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) +add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) +link_directories(${TBB_BINARY_DIR}) +###Eigen 3 NOTE: Eigen is required on the system the code is ran +find_package(Eigen3 3.3 REQUIRED) +>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb if(MSVC) add_compile_definitions(_HAS_STD_BYTE=0) endif(MSVC) - #link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs) file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp) add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp) - -set(USE_TBB true) -if(${USE_TBB}) - ##TBB - download_project(PROJ TBB - GIT_REPOSITORY https://github.com/wjakob/tbb.git - GIT_TAG master - PREFIX ${EXTERNAL_DEPS_DIR} - ${UPDATE_DISCONNECTED_IF_AVAILABLE} - ) - option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) - option(TBB_BUILD_STATIC "Build TBB static library" ON) - add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) - link_directories(${TBB_BINARY_DIR}) - if(${MYSOURCES_STATIC_LINK}) - target_link_libraries(${PROJECT_NAME} PUBLIC -static tbb_static) - else() - target_link_libraries(${PROJECT_NAME} PUBLIC tbb) +target_include_directories(${PROJECT_NAME} + PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph + PUBLIC ${vcglib_devel_SOURCE_DIR} + PUBLIC ${threed-beam-fea_SOURCE_DIR} + ) +if(${MYSOURCES_STATIC_LINK}) + message("Linking statically") + target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread gfortran quadmath) +else() + target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot ThreedBeamFEA tbb pthread) + if(${USE_POLYSCOPE}) + target_link_libraries(${PROJECT_NAME} polyscope) endif() endif() +target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs) if(USE_ENSMALLEN) ##ENSMALLEN + set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen) download_project(PROJ ENSMALLEN GIT_REPOSITORY https://github.com/mlpack/ensmallen.git GIT_TAG master + BINARY_DIR ${ENSMALLEN_BINARY_DIR} PREFIX ${EXTERNAL_DEPS_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE} ) add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR}) +<<<<<<< HEAD download_project(PROJ ARMADILLO GIT_REPOSITORY https://gitlab.com/conradsnicta/armadillo-code.git @@ -120,14 +161,18 @@ if(USE_ENSMALLEN) add_subdirectory(${ARMADILLO_SOURCE_DIR} ${ARMADILLO_BINARY_DIR}) target_include_directories(${PROJECT_NAME} PUBLIC armadillo ) +======= + set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo") + add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build) +>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb if(${MYSOURCES_STATIC_LINK}) target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a) else() target_link_libraries(${PROJECT_NAME} PUBLIC armadillo ensmallen) endif() - target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include) endif() +<<<<<<< HEAD target_include_directories(${PROJECT_NAME} PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph @@ -183,3 +228,5 @@ endif() #target_link_libraries(${PROJECT_NAME} PUBLIC GSL) #target_include_directories(${PROJECT_NAME} PUBLIC GSL) +======= +>>>>>>> c9fc6ccd0877f92308795b99e51c8e5dfd1e31eb diff --git a/beam.hpp b/beam.hpp index cb23c66..abf9112 100755 --- a/beam.hpp +++ b/beam.hpp @@ -9,15 +9,36 @@ struct RectangularBeamDimensions { inline static std::string name{"Rectangular"}; double b; double h; - RectangularBeamDimensions(const double &width, const double &height) - : b(width), h(height) { - assert(width > 0 && height > 0); + + double A{0}; // cross sectional area + + struct MomentsOfInertia + { + double I2{0}; // second moment of inertia + double I3{0}; // third moment of inertia + double J{0}; // torsional constant (polar moment of inertia) + } inertia; + + RectangularBeamDimensions(const double &width, const double &height) : b(width), h(height) + { + assert(width > 0 && height > 0); + updateProperties(); } - RectangularBeamDimensions() : b(0.002), h(0.002) {} + RectangularBeamDimensions() : b(0.002), h(0.002) { updateProperties(); } std::string toString() const { return std::string("b=") + std::to_string(b) + std::string(" h=") + std::to_string(h); } + + void updateProperties() + { + A = b * h; + inertia.I2 = b * std::pow(h, 3) / 12; + inertia.I3 = h * std::pow(b, 3) / 12; + inertia.J = inertia.I2 + inertia.I3; + } + static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions, + MomentsOfInertia &inertia); }; struct CylindricalBeamDimensions { @@ -39,18 +60,21 @@ struct ElementMaterial { double poissonsRatio; double youngsModulus; + double G; ElementMaterial(const double &poissonsRatio, const double &youngsModulus) : poissonsRatio(poissonsRatio), youngsModulus(youngsModulus) { assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); + updateProperties(); } - ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) {} + ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) { updateProperties(); } std::string toString() const { return std::string("Material:") + std::string("\nPoisson's ratio=") + std::to_string(poissonsRatio) + std::string("\nYoung's Modulus(GPa)=") + std::to_string(youngsModulus / 1e9); } + void updateProperties() { G = youngsModulus / (2 * (1 + poissonsRatio)); } }; #endif // BEAM_HPP diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 1b335b3..4326d6e 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -796,21 +796,21 @@ double DRMSimulationModel::computeTotalInternalPotentialEnergy() const EdgeIndex ei = pMesh->getIndex(e); const double e_k = element.length - element.initialLength; - const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A + const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.dimensions.A / (2 * element.initialLength); const double theta1Diff = theta1_jplus1 - theta1_j; - const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2) - / (2 * element.initialLength); + const double torsionalPE = element.material.G * element.dimensions.inertia.J + * pow(theta1Diff, 2) / (2 * element.initialLength); const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) + 4 * theta2_j * theta2_jplus1 + 4 * pow(theta2_jplus1, 2); const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus - * element.inertia.I2 / (2 * element.initialLength); + * element.dimensions.inertia.I2 / (2 * element.initialLength); const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) + 4 * theta3_j * theta3_jplus1 + 4 * pow(theta3_jplus1, 2); const double secondBendingPE = secondBendingPE_inBracketsTerm - * element.material.youngsModulus * element.inertia.I3 + * element.material.youngsModulus * element.dimensions.inertia.I3 / (2 * element.initialLength); totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE; @@ -1261,14 +1261,15 @@ void DRMSimulationModel::updateNodalMasses() for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) { const size_t ei = pMesh->getIndex(ep); const Element &elem = pMesh->elements[ep]; - const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; + const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A + / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2 + const long double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2 / lengthToThe3; - const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3 + const long double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3 / lengthToThe3; - const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J + const long double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J / lengthToThe3; rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1627,14 +1628,15 @@ void DRMSimulationModel::updatePositionsOnTheFly( double rotationalSumSk_J = 0; for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) { const Element &elem = pMesh->elements[ep]; - const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; + const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A + / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2 + const double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3 + const double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J + const double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1831,10 +1833,13 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr I3(pMesh->EN()); for (const EdgeType &e : pMesh->edge) { const size_t ei = pMesh->getIndex(e); - A[ei] = pMesh->elements[e].A; - J[ei] = pMesh->elements[e].inertia.J; - I2[ei] = pMesh->elements[e].inertia.I2; - I3[ei] = pMesh->elements[e].inertia.I3; + A[ei] = pMesh->elements[e].dimensions.A; + J[ei] = pMesh->elements[e].dimensions.inertia.J; + I2[ei] = pMesh->elements[e].dimensions.inertia.I2; + I3[ei] = pMesh->elements[e].dimensions.inertia.I3; } polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 10c3f7d..c3e7c39 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -2,7 +2,8 @@ #define BEAMFORMFINDER_HPP //#ifdef USE_MATPLOT -#include "matplot/matplot.h" +//#include "matplot.h" +#include //#endif #include "simulationmesh.hpp" #include "simulationmodel.hpp" diff --git a/linearsimulationmodel.cpp b/linearsimulationmodel.cpp index 11e6cd4..48a6f96 100644 --- a/linearsimulationmodel.cpp +++ b/linearsimulationmodel.cpp @@ -1,5 +1,5 @@ #include "linearsimulationmodel.hpp" -#include "gsl/gsl" +//#include "gsl/gsl" std::vector LinearSimulationModel::getFeaElements( const std::shared_ptr &pMesh) @@ -14,10 +14,10 @@ std::vector LinearSimulationModel::getFeaElements( (evn0[2] + evn1[2]) / 2}; const Element &element = pMesh->elements[edgeIndex]; const double E = element.material.youngsModulus; - fea::Props feaProperties(E * element.A, - E * element.inertia.I3, - E * element.inertia.I2, - element.G * element.inertia.J, + fea::Props feaProperties(E * element.dimensions.A, + E * element.dimensions.inertia.I3, + E * element.dimensions.inertia.I2, + element.material.G * element.dimensions.inertia.J, nAverageVector); const int vi0 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(0)); const int vi1 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(1)); @@ -139,26 +139,30 @@ SimulationResults LinearSimulationModel::getResults( const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.A); + * element.dimensions.A); const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.A); + * element.dimensions.A); const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0 + elementPotentialEnergy_axial_v1; //Shear const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2) * element.initialLength - / (4 * element.A * element.G); + / (4 * element.dimensions.A + * element.material.G); const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2) * element.initialLength - / (4 * element.A * element.G); + / (4 * element.dimensions.A + * element.material.G); const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2) * element.initialLength - / (4 * element.A * element.G); + / (4 * element.dimensions.A + * element.material.G); const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2) * element.initialLength - / (4 * element.A * element.G); + / (4 * element.dimensions.A + * element.material.G); const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0 + elementPotentialEnergy_shearZ_v0 + elementPotentialEnergy_shearY_v1 @@ -167,19 +171,19 @@ SimulationResults LinearSimulationModel::getResults( const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.inertia.I2); + * element.dimensions.inertia.I2); const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.inertia.I3); + * element.dimensions.inertia.I3); const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.inertia.I2); + * element.dimensions.inertia.I2); const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.inertia.I3); + * element.dimensions.inertia.I3); const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0 + elementPotentialEnergy_bendingZ_v0 + elementPotentialEnergy_bendingY_v1 @@ -187,10 +191,12 @@ SimulationResults LinearSimulationModel::getResults( //Torsion const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2) * element.initialLength - / (4 * element.inertia.J * element.G); + / (4 * element.dimensions.inertia.J + * element.material.G); const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2) * element.initialLength - / (4 * element.inertia.J * element.G); + / (4 * element.dimensions.inertia.J + * element.material.G); const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0 + elementPotentialEnergy_torsion_v1; const double elementInternalPotentialEnergy = elementPotentialEnergy_axial diff --git a/mesh.hpp b/mesh.hpp index c08ab86..7ff46eb 100755 --- a/mesh.hpp +++ b/mesh.hpp @@ -1,7 +1,7 @@ #ifndef MESH_HPP #define MESH_HPP -#include +//#include #include class Mesh diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index d301f3f..4cc695a 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -62,140 +62,159 @@ struct xRange } }; - struct Settings - { - enum NormalizationStrategy { NonNormalized, Epsilon }; - std::vector parameterRanges; - inline static vector normalizationStrategyStrings{"NonNormalized", "Epsilon"}; - int numberOfFunctionCalls{100000}; - double solverAccuracy{1e-2}; - NormalizationStrategy normalizationStrategy{Epsilon}; - double translationNormalizationParameter{0.003}; - double rotationNormalizationParameter{3}; - bool splitGeometryMaterialOptimization{false}; - struct ObjectiveWeights - { - double translational{1}; - double rotational{1}; - } objectiveWeights; +enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters }; - struct JsonKeys - { - inline static std::string filename{"OptimizationSettings.json"}; - inline static std::string OptimizationVariables{"OptimizationVariables"}; - inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"}; - inline static std::string SolverAccuracy{"SolverAccuracy"}; - inline static std::string TranslationalObjectiveWeight{"TransObjWeight"}; - }; +struct Settings +{ + std::filesystem::path intermediateResultsDirectoryPath; + std::vector> optimizationVariables; + enum NormalizationStrategy { NonNormalized, Epsilon }; + std::vector parameterRanges; + inline static vector normalizationStrategyStrings{"NonNormalized", "Epsilon"}; + int numberOfFunctionCalls{100000}; + double solverAccuracy{1e-2}; + NormalizationStrategy normalizationStrategy{Epsilon}; + double translationNormalizationParameter{0.003}; + double rotationNormalizationParameter{3}; + struct ObjectiveWeights + { + double translational{1}; + double rotational{1}; + } objectiveWeights; - void save(const std::filesystem::path &saveToPath) - { - assert(std::filesystem::is_directory(saveToPath)); + struct JsonKeys + { + inline static std::string filename{"OptimizationSettings.json"}; + inline static std::string OptimizationVariables{"OptimizationVariables"}; + inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"}; + inline static std::string SolverAccuracy{"SolverAccuracy"}; + inline static std::string TranslationalObjectiveWeight{"TransObjWeight"}; + inline static std::string OptimizationParameters{"OptimizationParameters"}; + }; - nlohmann::json json; - //write x ranges - for (size_t xRangeIndex = 0; xRangeIndex < parameterRanges.size(); xRangeIndex++) { - const auto &xRange = parameterRanges[xRangeIndex]; - json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)] - = xRange.toString(); - } + void save(const std::filesystem::path &saveToPath) + { + assert(std::filesystem::is_directory(saveToPath)); - json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls; - json[JsonKeys::SolverAccuracy] = solverAccuracy; - json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational; + nlohmann::json json; + //write x ranges + for (size_t xRangeIndex = 0; xRangeIndex < parameterRanges.size(); xRangeIndex++) { + const auto &xRange = parameterRanges[xRangeIndex]; + json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)] + = xRange.toString(); + } - std::filesystem::path jsonFilePath( - std::filesystem::path(saveToPath).append(JsonKeys::filename)); - std::ofstream jsonFile(jsonFilePath.string()); - jsonFile << json; - } + json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls; + json[JsonKeys::SolverAccuracy] = solverAccuracy; + json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational; + json[JsonKeys::OptimizationParameters] = optimizationVariables; - bool load(const std::filesystem::path &loadFromPath) - { - assert(std::filesystem::is_directory(loadFromPath)); - //Load optimal X - nlohmann::json json; - std::filesystem::path jsonFilepath( - std::filesystem::path(loadFromPath).append(JsonKeys::filename)); - if (!std::filesystem::exists(jsonFilepath)) { - return false; - } - std::ifstream ifs(jsonFilepath); - ifs >> json; + std::filesystem::path jsonFilePath( + std::filesystem::path(saveToPath).append(JsonKeys::filename)); + std::ofstream jsonFile(jsonFilePath.string()); + jsonFile << json; + } - //read x ranges - size_t xRangeIndex = 0; - while (true) { - const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_" - + std::to_string(xRangeIndex); - if (!json.contains(jsonXRangeKey)) { - break; - } - xRange x; - x.fromString(json.at(jsonXRangeKey)); - parameterRanges.push_back(x); - xRangeIndex++; - } + bool load(const std::filesystem::path &loadFromPath) + { + assert(std::filesystem::is_directory(loadFromPath)); + //Load optimal X + nlohmann::json json; + std::filesystem::path jsonFilepath( + std::filesystem::path(loadFromPath).append(JsonKeys::filename)); + if (!std::filesystem::exists(jsonFilepath)) { + return false; + } + std::ifstream ifs(jsonFilepath); + ifs >> json; - numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls); - solverAccuracy = json.at(JsonKeys::SolverAccuracy); - objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight); - objectiveWeights.rotational = 2 - objectiveWeights.translational; - return true; - } + //read x ranges + size_t xRangeIndex = 0; + while (true) { + const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_" + + std::to_string(xRangeIndex); + if (!json.contains(jsonXRangeKey)) { + break; + } + xRange x; + x.fromString(json.at(jsonXRangeKey)); + parameterRanges.push_back(x); + xRangeIndex++; + } - std::string toString() const - { - std::string settingsString; - if (!parameterRanges.empty()) { - std::string xRangesString; - for (const xRange &range : parameterRanges) { - xRangesString += range.toString() + " "; - } - settingsString += xRangesString; - } - settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) - + " Accuracy=" + std::to_string(solverAccuracy) - + " TransWeight=" + std::to_string(objectiveWeights.translational) - + " RotWeight=" + std::to_string(objectiveWeights.rotational); + optimizationVariables + = std::vector>( + (json.at(JsonKeys::OptimizationParameters))); + numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls); + solverAccuracy = json.at(JsonKeys::SolverAccuracy); + objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight); + objectiveWeights.rotational = 2 - objectiveWeights.translational; + return true; + } - return settingsString; - } + std::string toString() const + { + std::string settingsString; + if (!parameterRanges.empty()) { + std::string xRangesString; + for (const xRange &range : parameterRanges) { + xRangesString += range.toString() + " "; + } + settingsString += xRangesString; + } + settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) + + " Accuracy=" + std::to_string(solverAccuracy) + + " TransWeight=" + std::to_string(objectiveWeights.translational) + + " RotWeight=" + std::to_string(objectiveWeights.rotational); - void writeHeaderTo(csvFile &os) const - { - if (!parameterRanges.empty()) { - for (const xRange &range : parameterRanges) { - os << range.label + " max"; - os << range.label + " min"; - } - } - os << "Function Calls"; - os << "Solution Accuracy"; - os << "Normalization strategy"; - os << "Trans weight"; - os << "Rot weight"; - os << "Splitted geo from mat opt"; - // os << std::endl; - } + return settingsString; + } - void writeSettingsTo(csvFile &os) const - { - if (!parameterRanges.empty()) { - for (const xRange &range : parameterRanges) { - os << range.max; - os << range.min; - } - } - os << numberOfFunctionCalls; - os << solverAccuracy; - os << normalizationStrategyStrings[normalizationStrategy] + "_" - + std::to_string(translationNormalizationParameter); - os << objectiveWeights.translational; - os << objectiveWeights.rotational; - os << (splitGeometryMaterialOptimization == true ? "true" : "false"); - } - }; + void writeHeaderTo(csvFile &os) const + { + if (!parameterRanges.empty()) { + for (const xRange &range : parameterRanges) { + os << range.label + " max"; + os << range.label + " min"; + } + } + os << "Function Calls"; + os << "Solution Accuracy"; + os << "Normalization strategy"; + os << "Trans weight"; + os << "Rot weight"; + os << "Optimization parameters"; + // os << std::endl; + } + + void writeSettingsTo(csvFile &os) const + { + if (!parameterRanges.empty()) { + for (const xRange &range : parameterRanges) { + os << range.max; + os << range.min; + } + } + os << numberOfFunctionCalls; + os << solverAccuracy; + os << normalizationStrategyStrings[normalizationStrategy] + "_" + + std::to_string(translationNormalizationParameter); + os << objectiveWeights.translational; + os << objectiveWeights.rotational; + + //export optimization parameters + std::vector> vv; + for (const std::vector &v : optimizationVariables) { + std::vector vi; + vi.reserve(v.size()); + for (const OptimizationParameterIndex ¶meter : v) { + vi.emplace_back(parameter); + } + vv.push_back(vi); + } + os << Utilities::toString(vv); + } +}; inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept { @@ -222,9 +241,9 @@ struct xRange Element &e = m.elements[ei]; e.setDimensions(CrossSectionType(beamWidth, beamHeight)); e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); - e.inertia.J = J; - e.inertia.I2 = I2; - e.inertia.I3 = I3; + e.dimensions.inertia.J = J; + e.dimensions.inertia.I2 = I2; + e.dimensions.inertia.I3 = I3; } CoordType center_barycentric(1, 0, 0); @@ -606,7 +625,7 @@ struct xRange const double J = optimalXVariables.at(JLabel); for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; - e.inertia.J = J; + e.dimensions.inertia.J = J; } } @@ -615,7 +634,7 @@ struct xRange const double I2 = optimalXVariables.at(I2Label); for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; - e.inertia.I2 = I2; + e.dimensions.inertia.I2 = I2; } } @@ -624,7 +643,7 @@ struct xRange const double I3 = optimalXVariables.at(I3Label); for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; - e.inertia.I3 = I3; + e.dimensions.inertia.I3 = I3; } } pTiledReducedPattern_simulationMesh->reset(); diff --git a/simulation_structs.hpp b/simulation_structs.hpp index b6e2cbe..fcc2a4e 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -539,7 +539,7 @@ struct SimulationResults const std::filesystem::path outputFolderPath = outputFolder.empty() ? std::filesystem::current_path() : std::filesystem::path(outputFolder); - std::cout << "Saving results to:" << outputFolderPath << std::endl; + // std::cout << "Saving results to:" << outputFolderPath << std::endl; std::filesystem::path simulationJobOutputFolderPath = std::filesystem::path(outputFolderPath).append("SimulationJob"); std::filesystem::create_directories(simulationJobOutputFolderPath); @@ -569,8 +569,173 @@ struct SimulationResults void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom) { - //load job pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string()); + load(loadFromPath); + } + void load(const std::filesystem::path &loadFromPath, const std::shared_ptr &pJob) + { + this->pJob = pJob; + load(loadFromPath); + } +#ifdef POLYSCOPE_DEFINED + void unregister() const + { + if (!polyscope::hasCurveNetwork(getLabel())) { + std::cerr << "No curve network registered with a name: " << getLabel() << std::endl; + std::cerr << "Nothing to remove." << std::endl; + return; + } + pJob->unregister(getLabel()); + polyscope::removeCurveNetwork(getLabel()); + } + polyscope::CurveNetwork *registerForDrawing( + const std::optional> &desiredColor = std::nullopt, + const bool &shouldEnable = true, + const double &desiredRadius = 0.001, + // const double &desiredRadius = 0.0001, + const bool &shouldShowFrames = false) const + { + PolyscopeInterface::init(); + const std::shared_ptr &mesh = pJob->pMesh; + polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh; + if (!polyscope::hasStructure(getLabel())) { + const auto verts = mesh->getEigenVertices(); + const auto edges = mesh->getEigenEdges(); + polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(), + verts, + edges); + } else { + polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel()); + } + polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable); + polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true); + if (desiredColor.has_value()) { + const glm::vec3 desiredColor_glm(desiredColor.value()[0], + desiredColor.value()[1], + desiredColor.value()[2]); + polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm); + } + Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3); + Eigen::MatrixX3d framesX(mesh->VN(), 3); + Eigen::MatrixX3d framesY(mesh->VN(), 3); + Eigen::MatrixX3d framesZ(mesh->VN(), 3); + + Eigen::MatrixX3d framesX_initial(mesh->VN(), 3); + Eigen::MatrixX3d framesY_initial(mesh->VN(), 3); + Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3); + + // std::unordered_set interfaceNodes{1, 3, 5, 7, 9, 11}; + // std::unordered_set interfaceNodes{3, 7, 11, 15, 19, 23}; + // std::unordered_set interfaceNodes{}; + for (VertexIndex vi = 0; vi < mesh->VN(); vi++) { + const Vector6d &nodalDisplacement = displacements[vi]; + nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0], + nodalDisplacement[1], + nodalDisplacement[2]); + // Eigen::Quaternion Rx(Eigen::AngleAxis(nodalDisplacement[2],Eigen::Vector3d(1, 0, 0))); + // Eigen::Quaternion Ry(Eigen::AngleAxis(nodalDisplacement[4],Eigen::Vector3d(0, 1, 0))); + // Eigen::Quaternion Rz(Eigen::AngleAxis(nodalDisplacement[5],Eigen::Vector3d(0, 0, 1))); + // Eigen::Quaternion R=Rx*Ry*Rz; + // if (interfaceNodes.contains(vi)) { + auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1); + auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0); + auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0); + framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], + deformedFrameX[1], + deformedFrameX[2]); + framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], + deformedFrameY[1], + deformedFrameY[2]); + framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], + deformedNormal[1], + deformedNormal[2]); + framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0); + framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0); + framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1); + // } else { + // framesX.row(vi) = Eigen::Vector3d(0, 0, 0); + // framesY.row(vi) = Eigen::Vector3d(0, 0, 0); + // framesZ.row(vi) = Eigen::Vector3d(0, 0, 0); + // framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + // framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + // framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + // } + } + polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices() + + nodalDisplacements); + + const double frameRadius_default = 0.035; + const double frameLength_default = 0.035; + const bool shouldEnable_default = true; + //if(showFramesOn.contains(vi)){ + auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh + ->addNodeVectorQuantity("FrameX", framesX); + polyscopeHandle_frameX->setVectorLengthScale(frameLength_default); + polyscopeHandle_frameX->setVectorRadius(frameRadius_default); + polyscopeHandle_frameX->setVectorColor( + /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)); + auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh + ->addNodeVectorQuantity("FrameY", framesY); + polyscopeHandle_frameY->setVectorLengthScale(frameLength_default); + polyscopeHandle_frameY->setVectorRadius(frameRadius_default); + polyscopeHandle_frameY->setVectorColor( + /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)); + auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh + ->addNodeVectorQuantity("FrameZ", framesZ); + polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default); + polyscopeHandle_frameZ->setVectorRadius(frameRadius_default); + polyscopeHandle_frameZ->setVectorColor( + /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); + + auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel()); + if (!polyscopeHandle_initialMesh) { + polyscopeHandle_initialMesh = mesh->registerForDrawing(); + } + + // auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh + // ->addNodeVectorQuantity("FrameX", framesX_initial); + // polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default); + // polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default); + // polyscopeHandle_frameX_initial->setVectorColor( + // /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)); + // auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh + // ->addNodeVectorQuantity("FrameY", framesY_initial); + // polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default); + // polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default); + // polyscopeHandle_frameY_initial->setVectorColor( + // /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)); + // auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh + // ->addNodeVectorQuantity("FrameZ", framesZ_initial); + // polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default); + // polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default); + // polyscopeHandle_frameZ_initial->setVectorColor( + // /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); + // //} + pJob->registerForDrawing(getLabel()); + static bool wasExecuted = false; + if (!wasExecuted) { + std::function callback = [&]() { + static bool showFrames = shouldShowFrames; + + if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) { + polyscopeHandle_frameX->setEnabled(showFrames); + polyscopeHandle_frameY->setEnabled(showFrames); + polyscopeHandle_frameZ->setEnabled(showFrames); + } + }; + + PolyscopeInterface::addUserCallback(callback); + wasExecuted = true; + } + return polyscopeHandle_deformedEdmeMesh; + } +#endif +private: + void load(const std::filesystem::path &loadFromPath) + { + converged = true; //assuming it has converged + assert(pJob != nullptr); + //load job //Use the first .eigenBin file for loading the displacements for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) { if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") { @@ -591,156 +756,6 @@ struct SimulationResults Eigen::Vector3d::UnitZ()); } } - -#ifdef POLYSCOPE_DEFINED - void unregister() const { - if (!polyscope::hasCurveNetwork(getLabel())) { - std::cerr << "No curve network registered with a name: " << getLabel() - << std::endl; - std::cerr << "Nothing to remove." << std::endl; - return; - } - pJob->unregister(getLabel()); - polyscope::removeCurveNetwork(getLabel()); - } - polyscope::CurveNetwork *registerForDrawing( - const std::optional> &desiredColor = std::nullopt, - const bool &shouldEnable = true, - const double &desiredRadius = 0.001, - // const double &desiredRadius = 0.0001, - const bool &shouldShowFrames = false) const - { - PolyscopeInterface::init(); - const std::shared_ptr &mesh = pJob->pMesh; - polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh; - if (!polyscope::hasStructure(getLabel())) { - const auto verts = mesh->getEigenVertices(); - const auto edges = mesh->getEigenEdges(); - polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(), - verts, - edges); - } else { - polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel()); - } - polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable); - polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true); - if (desiredColor.has_value()) { - const glm::vec3 desiredColor_glm(desiredColor.value()[0], - desiredColor.value()[1], - desiredColor.value()[2]); - polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm); - } - Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3); - Eigen::MatrixX3d framesX(mesh->VN(), 3); - Eigen::MatrixX3d framesY(mesh->VN(), 3); - Eigen::MatrixX3d framesZ(mesh->VN(), 3); - - Eigen::MatrixX3d framesX_initial(mesh->VN(), 3); - Eigen::MatrixX3d framesY_initial(mesh->VN(), 3); - Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3); - - // std::unordered_set interfaceNodes{1, 3, 5, 7, 9, 11}; - // std::unordered_set interfaceNodes{3, 7, 11, 15, 19, 23}; - std::unordered_set interfaceNodes{}; - for (VertexIndex vi = 0; vi < mesh->VN(); vi++) { - const Vector6d &nodalDisplacement = displacements[vi]; - nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0], - nodalDisplacement[1], - nodalDisplacement[2]); - if (interfaceNodes.contains(vi)) { - auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1); - auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0); - auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0); - framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], - deformedFrameX[1], - deformedFrameX[2]); - framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], - deformedFrameY[1], - deformedFrameY[2]); - framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], - deformedNormal[1], - deformedNormal[2]); - framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0); - framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0); - framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1); - } else { - framesX.row(vi) = Eigen::Vector3d(0, 0, 0); - framesY.row(vi) = Eigen::Vector3d(0, 0, 0); - framesZ.row(vi) = Eigen::Vector3d(0, 0, 0); - framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0); - framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0); - framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0); - } - } - polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices() - + nodalDisplacements); - - const double frameRadius_default = 0.035; - const double frameLength_default = 0.035; - const bool shouldEnable_default = true; - //if(showFramesOn.contains(vi)){ - auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh - ->addNodeVectorQuantity("FrameX", framesX); - polyscopeHandle_frameX->setVectorLengthScale(frameLength_default); - polyscopeHandle_frameX->setVectorRadius(frameRadius_default); - polyscopeHandle_frameX->setVectorColor( - /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)); - auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh - ->addNodeVectorQuantity("FrameY", framesY); - polyscopeHandle_frameY->setVectorLengthScale(frameLength_default); - polyscopeHandle_frameY->setVectorRadius(frameRadius_default); - polyscopeHandle_frameY->setVectorColor( - /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)); - auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh - ->addNodeVectorQuantity("FrameZ", framesZ); - polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default); - polyscopeHandle_frameZ->setVectorRadius(frameRadius_default); - polyscopeHandle_frameZ->setVectorColor( - /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); - - auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel()); - if (!polyscopeHandle_initialMesh) { - polyscopeHandle_initialMesh = mesh->registerForDrawing(); - } - - // auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh - // ->addNodeVectorQuantity("FrameX", framesX_initial); - // polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default); - // polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default); - // polyscopeHandle_frameX_initial->setVectorColor( - // /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)); - // auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh - // ->addNodeVectorQuantity("FrameY", framesY_initial); - // polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default); - // polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default); - // polyscopeHandle_frameY_initial->setVectorColor( - // /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)); - // auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh - // ->addNodeVectorQuantity("FrameZ", framesZ_initial); - // polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default); - // polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default); - // polyscopeHandle_frameZ_initial->setVectorColor( - // /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); - // //} - pJob->registerForDrawing(getLabel()); - static bool wasExecuted = false; - if (!wasExecuted) { - std::function callback = [&]() { - static bool showFrames = shouldShowFrames; - - if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) { - polyscopeHandle_frameX->setEnabled(showFrames); - polyscopeHandle_frameY->setEnabled(showFrames); - polyscopeHandle_frameZ->setEnabled(showFrames); - } - }; - - PolyscopeInterface::addUserCallback(callback); - wasExecuted = true; - } - return polyscopeHandle_deformedEdmeMesh; - } -#endif }; #endif // SIMULATIONHISTORY_HPP diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index b85f923..1c803f3 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -114,93 +114,81 @@ struct SimulationResultsReporter { if (YvaluesToPlot.size() < 2) { return; } - auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); - std::vector colors(x.size(), 0.5); - std::vector markerSizes(x.size(), 5); + std::vector colors(YvaluesToPlot.size(), 0.5); + std::vector markerSizes(YvaluesToPlot.size(), 5); if (!markPoints.empty()) { for (const auto pointIndex : markPoints) { colors[pointIndex] = 0.9; markerSizes[pointIndex] = 14; } } + std::vector x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo); - // matplot::figure(true); - // matplot::hold(matplot::on); - - // ->marker_indices(history.redMarks) - // ->marker_indices(truncatedRedMarks) - // .marker_color("r") - // ->marker_size(1) - // ; - // auto greenMarksY = matplot::transform(history.greenMarks, - // [&](auto x) { return YvaluesToPlot[x]; }); - // matplot::scatter(history.greenMarks, greenMarksY)->color("green").marker_size(10); - // matplot::hold(matplot::off); } void createPlots(const SimulationHistory &history, const std::string &reportFolderPath, - const std::string &graphSuffix) { - const auto graphsFolder = - std::filesystem::path(reportFolderPath).append("Graphs"); - std::filesystem::remove_all(graphsFolder); - std::filesystem::create_directory(graphsFolder.string()); + const std::string &graphSuffix) + { + const auto graphsFolder = std::filesystem::path(reportFolderPath).append("Graphs"); + std::filesystem::remove_all(graphsFolder); + std::filesystem::create_directory(graphsFolder.string()); - if (!history.kineticEnergy.empty()) { - createPlot("Number of Iterations", - "Log of Kinetic Energy log", - history.kineticEnergy, - std::filesystem::path(graphsFolder) - .append("KineticEnergyLog_" + graphSuffix + ".png") - .string(), - history.redMarks); - } + if (!history.kineticEnergy.empty()) { + createPlot("Number of Iterations", + "Log of Kinetic Energy log", + history.kineticEnergy, + std::filesystem::path(graphsFolder) + .append("KineticEnergyLog_" + graphSuffix + ".png") + .string(), + history.redMarks); + } - if (!history.logResidualForces.empty()) { - createPlot("Number of Iterations", - "Residual Forces norm log", - history.logResidualForces, - std::filesystem::path(graphsFolder) - .append("ResidualForcesLog_" + graphSuffix + ".png") - .string(), - history.redMarks); - } + if (!history.logResidualForces.empty()) { + createPlot("Number of Iterations", + "Residual Forces norm log", + history.logResidualForces, + std::filesystem::path(graphsFolder) + .append("ResidualForcesLog_" + graphSuffix + ".png") + .string(), + history.redMarks); + } - if (!history.potentialEnergies.empty()) { - createPlot("Number of Iterations", - "Potential energy", - history.potentialEnergies, - std::filesystem::path(graphsFolder) - .append("PotentialEnergy_" + graphSuffix + ".png") - .string(), - history.redMarks); - } - if (!history.residualForcesMovingAverage.empty()) { - createPlot("Number of Iterations", - "Residual forces moving average", - history.residualForcesMovingAverage, - std::filesystem::path(graphsFolder) - .append("ResidualForcesMovingAverage_" + graphSuffix + ".png") - .string(), - history.redMarks); - } - // if (!history.residualForcesMovingAverageDerivativesLog.empty()) { - // createPlot("Number of Iterations", - // "Residual forces moving average derivative log", - // history.residualForcesMovingAverageDerivativesLog, - // std::filesystem::path(graphsFolder) - // .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png") - // .string()); - // } - if (!history.sumOfNormalizedDisplacementNorms.empty()) { - createPlot("Number of Iterations", - "Sum of normalized displacement norms", - history.sumOfNormalizedDisplacementNorms, - std::filesystem::path(graphsFolder) - .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") - .string(), - history.redMarks); - } + if (!history.potentialEnergies.empty()) { + createPlot("Number of Iterations", + "Potential energy", + history.potentialEnergies, + std::filesystem::path(graphsFolder) + .append("PotentialEnergy_" + graphSuffix + ".png") + .string(), + history.redMarks); + } + if (!history.residualForcesMovingAverage.empty()) { + createPlot("Number of Iterations", + "Residual forces moving average", + history.residualForcesMovingAverage, + std::filesystem::path(graphsFolder) + .append("ResidualForcesMovingAverage_" + graphSuffix + ".png") + .string(), + history.redMarks); + } + // if (!history.residualForcesMovingAverageDerivativesLog.empty()) { + // createPlot("Number of Iterations", + // "Residual forces moving average derivative log", + // history.residualForcesMovingAverageDerivativesLog, + // std::filesystem::path(graphsFolder) + // .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png") + // .string()); + // } + if (!history.sumOfNormalizedDisplacementNorms.empty()) { + createPlot("Number of Iterations", + "Sum of normalized displacement norms", + history.sumOfNormalizedDisplacementNorms, + std::filesystem::path(graphsFolder) + .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") + .string(), + history.redMarks); + } } }; diff --git a/simulationmesh.cpp b/simulationmesh.cpp index 05afde5..446a3eb 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -198,7 +198,6 @@ void SimulationMesh::reset() { } void SimulationMesh::initializeElements() { - computeElementalProperties(); for (const EdgeType &e : edge) { Element &element = elements[e]; element.ei = getIndex(e); @@ -258,7 +257,6 @@ void SimulationMesh::setBeamCrossSection( const CrossSectionType &beamDimensions) { for (size_t ei = 0; ei < EN(); ei++) { elements[ei].dimensions = beamDimensions; - elements[ei].computeDimensionsProperties(beamDimensions); elements[ei].updateRigidity(); } } @@ -455,60 +453,59 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, return angle; } -void Element::computeMaterialProperties(const ElementMaterial &material) { - G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); -} +//void Element::computeMaterialProperties(const ElementMaterial &material) { +// G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); +//} -void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A) -{ - A = dimensions.b * dimensions.h; -} +//void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A) +//{ +// A = dimensions.b * dimensions.h; +//} -void Element::computeDimensionsProperties( - const RectangularBeamDimensions &dimensions) { - assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); - computeCrossSectionArea(dimensions, A); - computeMomentsOfInertia(dimensions, inertia); -} +//void Element::computeDimensionsProperties( +// const RectangularBeamDimensions &dimensions) { +// assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); +// computeCrossSectionArea(dimensions, A); +// computeMomentsOfInertia(dimensions, dimensions.inertia); +//} -void Element::computeDimensionsProperties( - const CylindricalBeamDimensions &dimensions) { - assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); - A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2)); - inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64; - inertia.I3 = inertia.I2; - inertia.J = inertia.I2 + inertia.I3; -} +//void Element::computeDimensionsProperties( +// const CylindricalBeamDimensions &dimensions) { +// assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); +// A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2)); +// dimensions.inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64; +// dimensions.inertia.I3 = dimensions.inertia.I2; +// dimensions.inertia.J = dimensions.inertia.I2 + dimensions.inertia.I3; +//} void Element::setDimensions(const CrossSectionType &dimensions) { this->dimensions = dimensions; - computeDimensionsProperties(dimensions); + assert(this->dimensions.A == dimensions.A); + // computeDimensionsProperties(dimensions); updateRigidity(); } -void Element::setMaterial(const ElementMaterial &material) { - this->material = material; - computeMaterialProperties(material); - updateRigidity(); +void Element::setMaterial(const ElementMaterial &material) +{ + this->material = material; + // computeMaterialProperties(material); + updateRigidity(); } double Element::getMass(const double &materialDensity) { - const double beamVolume = A * length; + const double beamVolume = dimensions.A * length; return beamVolume * materialDensity; } -void Element::computeMomentsOfInertia(const RectangularBeamDimensions &dimensions, - Element::MomentsOfInertia &inertia) -{ - inertia.I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; - inertia.I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; - inertia.J = inertia.I2 + inertia.I3; -} - void Element::updateRigidity() { - rigidity.axial = material.youngsModulus * A / initialLength; - rigidity.torsional = G * inertia.J / initialLength; - rigidity.firstBending = 2 * material.youngsModulus * inertia.I2 / initialLength; - rigidity.secondBending = 2 * material.youngsModulus * inertia.I3 / initialLength; + // assert(initialLength != 0); + rigidity.axial = material.youngsModulus * dimensions.A / initialLength; + // assert(rigidity.axial != 0); + rigidity.torsional = material.G * dimensions.inertia.J / initialLength; + // assert(rigidity.torsional != 0); + rigidity.firstBending = 2 * material.youngsModulus * dimensions.inertia.I2 / initialLength; + // assert(rigidity.firstBending != 0); + rigidity.secondBending = 2 * material.youngsModulus * dimensions.inertia.I3 / initialLength; + // assert(rigidity.secondBending != 0); } diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 35cc499..791e9ea 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -1,9 +1,9 @@ #ifndef SIMULATIONMESH_HPP #define SIMULATIONMESH_HPP -#include "Eigen/Dense" #include "edgemesh.hpp" #include "trianglepatterngeometry.hpp" +#include //extern bool drawGlobal; struct Element; @@ -62,26 +62,14 @@ struct Element { CrossSectionType dimensions; ElementMaterial material; - double G{0}; - double A{0}; // cross sectional area void computeMaterialProperties(const ElementMaterial &material); - void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); - void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); + // void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); + // void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); void setDimensions(const CrossSectionType &dimensions); void setMaterial(const ElementMaterial &material); double getMass(const double &matrialDensity); - struct MomentsOfInertia - { - double I2{0}; // second moment of inertia - double I3{0}; // third moment of inertia - double J{0}; // torsional constant (polar moment of inertia) - } inertia; - - static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions, - MomentsOfInertia &inertia); - struct LocalFrame { VectorType t1; VectorType t2; diff --git a/utilities.hpp b/utilities.hpp index 9dce841..ad42f77 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -2,12 +2,12 @@ #define UTILITIES_H #include -#include -#include -#include -#include #include #include +#include +#include +#include +#include #define GET_VARIABLE_NAME(Variable) (#Variable) @@ -140,17 +140,85 @@ struct Vector6d : public std::array { }; namespace Utilities { -inline void parseIntegers(const std::string &str, std::vector &result) { - typedef std::regex_iterator re_iterator; - typedef re_iterator::value_type re_iterated; +//inline std::vector split(std::string text, char delim) +//{ +// std::string line; +// std::vector vec; +// std::stringstream ss(text); +// while (std::getline(ss, line, delim)) { +// vec.push_back(line); +// } +// return vec; +//} - std::regex re("(\\d+)"); +inline std::string_view leftTrimSpaces(const std::string_view& str) +{ +std::string_view trimmedString=str; + const auto pos(str.find_first_not_of(" \t\n\r\f\v")); + trimmedString.remove_prefix(std::min(pos, trimmedString.length())); + return trimmedString; +} - re_iterator rit(str.begin(), str.end(), re); - re_iterator rend; +inline std::string_view rightTrimSpaces(const std::string_view& str) +{ +std::string_view trimmedString=str; + const auto pos(trimmedString.find_last_not_of(" \t\n\r\f\v")); + trimmedString.remove_suffix(std::min(trimmedString.length() - pos - 1,trimmedString.length())); + return trimmedString; +} - std::transform(rit, rend, std::back_inserter(result), - [](const re_iterated &it) { return std::stoi(it[1]); }); +inline std::string_view trimLeftAndRightSpaces(std::string_view str) +{ +std::string_view trimmedString=str; + trimmedString = leftTrimSpaces(trimmedString); + trimmedString = rightTrimSpaces(trimmedString); + return trimmedString; +} + +inline std::vector split(const std::string& text, std::string delim) +{ + std::vector vec; + size_t pos = 0, prevPos = 0; + while (1) { + pos = text.find(delim, prevPos); + if (pos == std::string::npos) { + vec.push_back(text.substr(prevPos)); + return vec; + } + + vec.push_back(text.substr(prevPos, pos - prevPos)); + prevPos = pos + delim.length(); + } +} + +inline std::string toString(const std::vector> &vv) +{ + std::string s; + s.append("{"); + for (const std::vector &v : vv) { + s.append("{"); + for (const int &i : v) { + s.append(std::to_string(i) + ","); + } + s.pop_back(); + s.append("}"); + } + s.append("}"); + return s; +} +inline void parseIntegers(const std::string &str, std::vector &result) +{ + typedef std::regex_iterator re_iterator; + typedef re_iterator::value_type re_iterated; + + std::regex re("(\\d+)"); + + re_iterator rit(str.begin(), str.end(), re); + re_iterator rend; + + std::transform(rit, rend, std::back_inserter(result), [](const re_iterated &it) { + return std::stoi(it[1]); + }); } inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { @@ -179,7 +247,6 @@ inline std::vector fromEigenMatrix(const Eigen::MatrixXd &m) return v; } - // std::string convertToLowercase(const std::string &s) { // std::string lowercase; // std::transform(s.begin(), s.end(), lowercase.begin(),