CMakeLists.txt

@@ -1,10 +1,7 @@
 cmake_minimum_required(VERSION 2.8)
 project(MySources)
-
+set(CMAKE_CXX_STANDARD 20)
-file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
-add_library(${PROJECT_NAME} ${MySourcesFiles} )
-target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20)
-
 #Download external dependencies NOTE: If the user has one of these libs it shouldn't be downloaded again.
 include(${CMAKE_CURRENT_LIST_DIR}/cmake/DownloadProject.cmake)
 if (CMAKE_VERSION VERSION_LESS 3.2)
@@ -12,114 +9,101 @@ if (CMAKE_VERSION VERSION_LESS 3.2)
 else()
     set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1")
 endif()
+
 set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies")
 if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
     set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR})
 endif()
 ##Create directory for the external libraries
 file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
+#message(${POLYSCOPE_ALREADY_COMPILED})
+
 if(${USE_POLYSCOPE})
-        set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
+download_project(PROJ               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})
+if(NOT EXISTS ${POLYSCOPE_BINARY_DIR})
-        add_compile_definitions(POLYSCOPE_DEFINED)
+add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
-        target_sources(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
+add_compile_definitions(POLYSCOPE_DEFINED)
 endif()
+endif()
+
+#dlib
+set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
+file(MAKE_DIRECTORY ${DLIB_BIN_DIR})
+download_project(PROJ               DLIB
+                 GIT_REPOSITORY     https://github.com/davisking/dlib.git
+                 GIT_TAG	    master
+                 BINARY_DIR         ${DLIB_BIN_DIR}
+                 PREFIX             ${EXTERNAL_DEPS_DIR}
+                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+)
+add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR})
 
 ##vcglib devel branch
 download_project(PROJ		    vcglib_devel
-                 GIT_REPOSITORY     https://gitea-s2i2s.isti.cnr.it/manolas/vcglib.git
+                 GIT_REPOSITORY     https://github.com/IasonManolas/vcglib.git
                  GIT_TAG	    devel
                  PREFIX             ${EXTERNAL_DEPS_DIR}
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
-add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
-target_sources(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
 
 ##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://github.com/IasonManolas/threed-beam-fea.git
+                 GIT_REPOSITORY      https://github.com/IasonManolas/threed-beam-fea.git
-                 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})
 
 ##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)
+
 if(MSVC)
         add_compile_definitions(_HAS_STD_BYTE=0)
 endif(MSVC)
-#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
 
-if(${MYSOURCES_STATIC_LINK})
+#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
-        message("Linking statically")
+file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
-        target_link_libraries(${PROJECT_NAME} PUBLIC -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a #[[tbb_static]] pthread gfortran quadmath)
+add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
-else()
+
-#    set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
+
-        target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
-        if(${USE_POLYSCOPE})
-        message("Using polyscope")
-                target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
-        endif()
-endif()
-target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
 target_include_directories(${PROJECT_NAME}
     PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
     PUBLIC ${vcglib_devel_SOURCE_DIR}
     PUBLIC ${threed-beam-fea_SOURCE_DIR}
-        PUBLIC matplot
+    PUBLIC ${MySourcesFiles}
     )
 
-if(USE_ENSMALLEN)
+if(${MYSOURCES_STATIC_LINK})
-        ##ENSMALLEN
+        target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
-        set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen)
+elseif(${USE_POLYSCOPE})
-        download_project(PROJ			ENSMALLEN
+        target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread polyscope)
-                         GIT_REPOSITORY		https://github.com/mlpack/ensmallen.git
+else()
-                         GIT_TAG		master
+        target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread)
-                         BINARY_DIR		${ENSMALLEN_BINARY_DIR}
-                         PREFIX             	${EXTERNAL_DEPS_DIR}
-                         ${UPDATE_DISCONNECTED_IF_AVAILABLE}
-        )
-        add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
-        set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo")
-        add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build)
-        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 ${ENSMALLEN_SOURCE_DIR})
 endif()
+target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)

beam.hpp

@@ -9,36 +9,15 @@ struct RectangularBeamDimensions {
   inline static std::string name{"Rectangular"};
   double b;
   double h;
-
+  RectangularBeamDimensions(const double &width, const double &height)
-  double A{0}; // cross sectional area
+      : b(width), h(height) {
-
-  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) { updateProperties(); }
+  RectangularBeamDimensions() : b(0.002), h(0.002) {}
   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 {
@@ -60,21 +39,18 @@ 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) { updateProperties(); }
+    ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) {}
     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

csvfile.hpp

@@ -21,19 +21,18 @@ class csvFile {
   const std::string escape_seq_;
   const std::string special_chars_;
 
-  public:
+public:
-  csvFile(const std::filesystem::path &filename,
+  csvFile(const std::string filename, const bool &overwrite,
-          const bool &overwrite,
           const std::string separator = ",")
-      : fs_(), is_first_(true), separator_(separator), escape_seq_("\""), special_chars_("\"")
+      : fs_(), is_first_(true), separator_(separator), escape_seq_("\""),
-  {
+        special_chars_("\"") {
     fs_.exceptions(std::ios::failbit | std::ios::badbit);
     if (filename.empty()) {
       fs_.copyfmt(std::cout);
       fs_.clear(std::cout.rdstate());
       fs_.basic_ios<char>::rdbuf(std::cout.rdbuf());
     } else {
-          if (!std::filesystem::exists(filename)) {
+        if (!std::filesystem::exists(std::filesystem::path(filename))) {
             std::ofstream outfile(filename);
             outfile.close();
         }

drmsimulationmodel.cpp

@@ -217,16 +217,28 @@ void DRMSimulationModel::runUnitTests()
 
     //  polyscope::show();
 }
-void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
+
+void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
 {
-    //#ifdef USE_ENSMALLEN
+    mSettings = settings;
-    //    this->pJob = pJob;
+    mCurrentSimulationStep = 0;
-    //#endif
+    history.clear();
+    history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
+    constrainedVertices.clear();
     pMesh.reset();
     pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
     vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
+    plotYValues.clear();
+    plotHandle.reset();
+    checkedForMaximumMoment = false;
+    externalMomentsNorm = 0;
+    Dt = settings.Dtini;
+    numOfDampings = 0;
+    shouldTemporarilyDampForces = false;
+    externalLoadStep = 1;
+    isVertexConstrained.clear();
+    minTotalResidualForcesNorm = std::numeric_limits<double>::max();
 
-    constrainedVertices.clear();
     constrainedVertices = pJob->constrainedVertices;
     if (!pJob->nodalForcedDisplacements.empty()) {
         for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : pJob->nodalForcedDisplacements) {
@@ -239,26 +251,6 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
     for (auto fixedVertex : constrainedVertices) {
         isVertexConstrained[fixedVertex.first] = true;
     }
-}
-
-void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
-{
-    mSettings = settings;
-    mCurrentSimulationStep = 0;
-    history.clear();
-    history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
-    plotYValues.clear();
-    plotHandle.reset();
-    checkedForMaximumMoment = false;
-    externalMomentsNorm = 0;
-    Dt = settings.Dtini;
-    numOfDampings = 0;
-    shouldTemporarilyDampForces = false;
-    externalLoadStep = 1;
-    isVertexConstrained.clear();
-    minTotalResidualForcesNorm = std::numeric_limits<double>::max();
-
-    reset(pJob);
 
 #ifdef POLYSCOPE_DEFINED
     if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
@@ -796,21 +788,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.dimensions.A
+        const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A
                                / (2 * element.initialLength);
         const double theta1Diff = theta1_jplus1 - theta1_j;
-        const double torsionalPE = element.material.G * element.dimensions.inertia.J
+        const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2)
-                                   * pow(theta1Diff, 2) / (2 * element.initialLength);
+                                   / (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.dimensions.inertia.I2 / (2 * element.initialLength);
+                                      * element.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.dimensions.inertia.I3
+                                       * element.material.youngsModulus * element.inertia.I3
                                        / (2 * element.initialLength);
 
         totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE;
@@ -901,8 +893,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
             //      #pragma omp parallel for schedule(static) num_threads(6)
             for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
                 const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
-                                                        && fixedVertices.at(edgeNode.vi)
+                                                    && fixedVertices.at(edgeNode.vi).contains(dofi);
-                                                               .contains(dofi);
                 if (!isDofConstrainedFor_ev) {
                     DifferentiateWithRespectTo dui{ev, dofi};
                     // Axial force computation
@@ -1073,6 +1064,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
         }
 #endif
     }
+    pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
     pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
     pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
     //    pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
@@ -1261,15 +1253,14 @@ 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.dimensions.A
+            const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
-                                           / elem.length;
             translationalSumSk += SkTranslational;
             const double lengthToThe3 = std::pow(elem.length, 3);
-            const long double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
+            const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
                                                 / lengthToThe3;
-            const long double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
+            const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
                                                 / lengthToThe3;
-            const long double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
+            const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
                                                / lengthToThe3;
             rotationalSumSk_I2 += SkRotational_I2;
             rotationalSumSk_I3 += SkRotational_I3;
@@ -1628,15 +1619,14 @@ 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.dimensions.A
+            const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
-                                           / elem.length;
             translationalSumSk += SkTranslational;
             const double lengthToThe3 = std::pow(elem.length, 3);
-            const double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
+            const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
                                            / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
-            const double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
+            const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
                                            / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
-            const double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
+            const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
                                           / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
             rotationalSumSk_I2 += SkRotational_I2;
             rotationalSumSk_I3 += SkRotational_I3;
@@ -1833,13 +1823,10 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
             results.debug_q_normal[vi] = q_normal;
             results.debug_q_nr[vi] = q_nr_nInit;
             results.rotationalDisplacementQuaternion[vi]
-                //Eigen::Quaterniond R
                 = (q_normal
                    * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal;
             //Update the displacement vector to contain the euler angles
-            const Eigen::Vector3d eulerAngles = results
+            const Eigen::Vector3d eulerAngles = results.rotationalDisplacementQuaternion[vi]
-                                                    .rotationalDisplacementQuaternion[vi]
-                                                    //            R
                                                     .toRotationMatrix()
                                                     .eulerAngles(0, 1, 2);
             results.displacements[vi][3] = eulerAngles[0];
@@ -1889,7 +1876,7 @@ void DRMSimulationModel::printCurrentState() const
 void DRMSimulationModel::printDebugInfo() const
 {
     std::cout << pMesh->elements[0].rigidity.toString() << std::endl;
-    std::cout << "Number of dampings:" << numOfDampings << std::endl;
+    std::cout << "Number of dampings:" << numOfDampings << endl;
     printCurrentState();
 }
 
@@ -2012,10 +1999,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
     std::vector<double> I3(pMesh->EN());
     for (const EdgeType &e : pMesh->edge) {
         const size_t ei = pMesh->getIndex(e);
-        A[ei] = pMesh->elements[e].dimensions.A;
+        A[ei] = pMesh->elements[e].A;
-        J[ei] = pMesh->elements[e].dimensions.inertia.J;
+        J[ei] = pMesh->elements[e].inertia.J;
-        I2[ei] = pMesh->elements[e].dimensions.inertia.I2;
+        I2[ei] = pMesh->elements[e].inertia.I2;
-        I3[ei] = pMesh->elements[e].dimensions.inertia.I3;
+        I3[ei] = pMesh->elements[e].inertia.I3;
     }
 
     polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);
@@ -2669,119 +2656,54 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
     return results;
 }
 
-//#ifdef USE_ENSMALLEN
+void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const
-
-//void DRMSimulationModel::setJob(const std::shared_ptr<SimulationJob> &pJob)
-//{
-//    reset(pJob);
-
-//    updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
-//    //    arma::mat externalForces(pMesh->VN() * NumDoF);
-//    //    for (int vi = 0; vi < pMesh->VN(); vi++) {
-//    //        const Vector6d &nodeExternalForce = pMesh.vert[vi].node.externalForce;
-//    //        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
-//    //            externalForces(vi * NumDoF + dofi) = nodeExternalForce[dofi];
-//    //        }
-//    //    }
-//}
-
-//SimulationMesh *DRMSimulationModel::getDeformedMesh(const arma::mat &x)
-//{
-//    for (int vi = 0; vi < pMesh->VN(); vi++) {
-//        Node &node = pMesh->nodes[vi];
-//        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
-//            node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
-//        }
-//    }
-
-//    applyDisplacements(constrainedVertices);
-//    if (!pJob->nodalForcedDisplacements.empty()) {
-//        applyForcedDisplacements(pJob->nodalForcedDisplacements);
-//    }
-
-//    return pMesh.get();
-//}
-
-//double DRMSimulationModel::EvaluateWithGradient(const arma::mat &x, arma::mat &g)
-//{
-//    for (int vi = 0; vi < pMesh->VN(); vi++) {
-//        Node &node = pMesh->nodes[vi];
-//        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
-//            node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
-//        }
-//    }
-
-//    applyDisplacements(constrainedVertices);
-//    if (!pJob->nodalForcedDisplacements.empty()) {
-//        applyForcedDisplacements(pJob->nodalForcedDisplacements);
-//    }
-//    updateElementalLengths();
-//    updateNormalDerivatives();
-//    updateT1Derivatives();
-//    updateRDerivatives();
-//    updateT2Derivatives();
-//    updateT3Derivatives();
-//    updateResidualForcesOnTheFly(constrainedVertices);
-//    for (int vi = 0; vi < pMesh->VN(); vi++) {
-//        Node &node = pMesh->nodes[vi];
-//        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
-//            g(vi * NumDoF + dofi) = node.force.residual[dofi];
-//        }
-//    }
-//    mCurrentSimulationStep++;
-//    const double PE = computeTotalPotentialEnergy();
-//    std::cout << "PE:" << PE << std::endl;
-//    return PE;
-//}
-
-//#endif
-
-void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePath) const
 {
+    std::filesystem::create_directories(folderPath);
     nlohmann::json json;
-    json[GET_VARIABLE_NAME(shouldDraw)] = shouldDraw;
+
-    json[GET_VARIABLE_NAME(beVerbose)] = beVerbose;
+    json[jsonLabels.gamma] = gamma;
-    json[GET_VARIABLE_NAME(shouldCreatePlots)] = shouldCreatePlots;
+    json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false";
-    json[GET_VARIABLE_NAME(Dtini)] = Dtini;
+    json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false";
-    json[GET_VARIABLE_NAME(xi)] = xi;
+    json[jsonLabels.beVerbose] = beVerbose ? "true" : "false";
-    json[GET_VARIABLE_NAME(gamma)] = gamma;
+    json[jsonLabels.Dtini] = Dtini;
-    json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold;
+    json[jsonLabels.xi] = xi;
     if (maxDRMIterations.has_value()) {
-        json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value();
+        json[jsonLabels.maxDRMIterations] = maxDRMIterations.value();
     }
     if (debugModeStep.has_value()) {
-        json[GET_VARIABLE_NAME(debugModeStep)] = debugModeStep.value();
+        json[jsonLabels.debugModeStep] = debugModeStep.value();
     }
     if (totalTranslationalKineticEnergyThreshold.has_value()) {
-        json[GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold)]
+        json[jsonLabels.totalTranslationalKineticEnergyThreshold]
             = totalTranslationalKineticEnergyThreshold.value();
     }
+
     if (averageResidualForcesCriterionThreshold.has_value()) {
-        json[GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold)]
+        json[jsonLabels.averageResidualForcesCriterionThreshold]
+            = averageResidualForcesCriterionThreshold.value();
+    }
+
+    if (averageResidualForcesCriterionThreshold.has_value()) {
+        json[jsonLabels.averageResidualForcesCriterionThreshold]
             = averageResidualForcesCriterionThreshold.value();
     }
     if (linearGuessForceScaleFactor.has_value()) {
-        json[GET_VARIABLE_NAME(linearGuessForceScaleFactor)] = linearGuessForceScaleFactor.value();
+        json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value();
     }
-    if (viscousDampingFactor.has_value()) {
+
-        json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value();
-    }
-    json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping;
     //    if (intermediateResultsSaveStep.has_value()) {
     //        json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value();
     //    }
 
-    //    if (saveIntermediateBestStates.has_value()) {
+    if (saveIntermediateBestStates.has_value()) {
-    //        json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
+        json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
-    //                                                                  ? "true"
+                                                                  ? "true"
-    //                                                                  : "false";
+                                                                  : "false";
-    //    }
+    }
 
-    std::filesystem::create_directories(jsonFilePath.parent_path());
+    const std::string jsonFilename = "drmSettings.json";
-    std::ofstream jsonFile(jsonFilePath);
+    std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
-    std::cout << "Saving DRM settings to:" << jsonFilePath << std::endl;
     jsonFile << json;
-    jsonFile.close();
 }
 
 bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath)
@@ -2805,40 +2727,41 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
     std::ifstream ifs(jsonFilePath.string());
     ifs >> json;
 
-    if (json.contains(GET_VARIABLE_NAME(shouldDraw))) {
+    if (json.contains(jsonLabels.shouldDraw)) {
-        shouldDraw = json.at(GET_VARIABLE_NAME(shouldDraw)) == "true" ? true : false;
+        shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false;
     }
-    if (json.contains(GET_VARIABLE_NAME(beVerbose))) {
+    if (json.contains(jsonLabels.beVerbose)) {
-        beVerbose = json.at(GET_VARIABLE_NAME(beVerbose)) == "true" ? true : false;
+        beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false;
     }
-    if (json.contains(GET_VARIABLE_NAME(shouldCreatePlots))) {
+    if (json.contains(jsonLabels.shouldCreatePlots)) {
-        shouldCreatePlots = json.at(GET_VARIABLE_NAME(shouldCreatePlots)) == "true" ? true : false;
+        shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false;
     }
-    if (json.contains(GET_VARIABLE_NAME(Dtini))) {
+    if (json.contains(jsonLabels.Dtini)) {
-        Dtini = json.at(GET_VARIABLE_NAME(Dtini));
+        Dtini = json.at(jsonLabels.Dtini);
     }
-    if (json.contains(GET_VARIABLE_NAME(xi))) {
+    if (json.contains(jsonLabels.xi)) {
-        xi = json.at(GET_VARIABLE_NAME(xi));
+        xi = json.at(jsonLabels.xi);
     }
-    if (json.contains(GET_VARIABLE_NAME(maxDRMIterations))) {
+    if (json.contains(jsonLabels.maxDRMIterations)) {
-        maxDRMIterations = json.at(GET_VARIABLE_NAME(maxDRMIterations));
+        maxDRMIterations = json.at(jsonLabels.maxDRMIterations);
     }
-    if (json.contains(GET_VARIABLE_NAME(debugModeStep))) {
+    if (json.contains(jsonLabels.debugModeStep)) {
-        debugModeStep = json.at(GET_VARIABLE_NAME(debugModeStep));
+        debugModeStep = json.at(jsonLabels.debugModeStep);
     }
-    if (json.contains(GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold))) {
+    if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) {
         totalTranslationalKineticEnergyThreshold = json.at(
-            GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold));
+            jsonLabels.totalTranslationalKineticEnergyThreshold);
     }
-    if (json.contains(GET_VARIABLE_NAME(gamma))) {
+    if (json.contains(jsonLabels.gamma)) {
-        gamma = json.at(GET_VARIABLE_NAME(gamma));
+        gamma = json.at(jsonLabels.gamma);
     }
-    if (json.contains(GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold))) {
+    if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) {
         averageResidualForcesCriterionThreshold = json.at(
-            GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold));
+            jsonLabels.averageResidualForcesCriterionThreshold);
     }
-    if (json.contains(GET_VARIABLE_NAME(linearGuessForceScaleFactor))) {
+
-        linearGuessForceScaleFactor = json.at(GET_VARIABLE_NAME(linearGuessForceScaleFactor));
+    if (json.contains(jsonLabels.linearGuessForceScaleFactor)) {
+        linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor);
     }
 
     //    if (json.contains(GET_VARIABLE_NAME(intermediateResultsSaveStep))) {
@@ -2851,7 +2774,6 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
                                          ? true
                                          : false;
     }
-    //    std::cout << json.dump() << std::endl;
 
     return true;
 }

drmsimulationmodel.hpp

@@ -1,21 +1,13 @@
 #ifndef BEAMFORMFINDER_HPP
 #define BEAMFORMFINDER_HPP
 
-//#ifdef USE_MATPLOT
+#include "matplot/matplot.h"
-//#include "matplot.h"
-#include <matplot/matplot.h>
-//#endif
 #include "simulationmesh.hpp"
 #include "simulationmodel.hpp"
 #include <Eigen/Dense>
 #include <filesystem>
 #include <unordered_set>
 
-#ifdef USE_ENSMALLEN
-#include <armadillo>
-#include <ensmallen.hpp>
-#endif
-
 struct SimulationJob;
 
 enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
@@ -33,6 +25,7 @@ class DRMSimulationModel : public SimulationModel
 public:
     struct Settings
     {
+        //        bool isDebugMode{false};
         bool shouldDraw{false};
         bool beVerbose{false};
         bool shouldCreatePlots{false};
@@ -45,8 +38,8 @@ public:
         int gradualForcedDisplacementSteps{50};
         //        int desiredGradualExternalLoadsSteps{1};
         double gamma{0.8};
-        double totalResidualForcesNormThreshold{1e-8};
+        double totalResidualForcesNormThreshold{1e-3};
-        double totalExternalForcesNormPercentageTermination{1e-8};
+        double totalExternalForcesNormPercentageTermination{1e-3};
         std::optional<int> maxDRMIterations;
         std::optional<int> debugModeStep;
         std::optional<double> totalTranslationalKineticEnergyThreshold;
@@ -54,11 +47,11 @@ public:
         std::optional<double> linearGuessForceScaleFactor;
         //        std::optional<int> intermediateResultsSaveStep;
         std::optional<bool> saveIntermediateBestStates;
+        Settings() {}
+        void save(const std::filesystem::path &folderPath) const;
+        bool load(const std::filesystem::path &filePath);
         std::optional<double> viscousDampingFactor;
         bool useKineticDamping{true};
-        Settings() {}
-        void save(const std::filesystem::path &jsonFilePath) const;
-        bool load(const std::filesystem::path &filePath);
         struct JsonLabels
         {
             const std::string shouldDraw{"shouldDraw"};
@@ -66,16 +59,14 @@ public:
             const std::string shouldCreatePlots{"shouldCreatePlots"};
             const std::string Dtini{"DtIni"};
             const std::string xi{"xi"};
-            const std::string gamma{"gamma"};
-            const std::string totalResidualForcesNormThreshold;
             const std::string maxDRMIterations{"maxIterations"};
             const std::string debugModeStep{"debugModeStep"};
             const std::string totalTranslationalKineticEnergyThreshold{
                 "totalTranslationaKineticEnergyThreshold"};
+            const std::string gamma{"gamma"};
             const std::string averageResidualForcesCriterionThreshold{
                 "averageResidualForcesThreshold"};
             const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
-            const std::string viscousDampingFactor{"viscousDampingFactor"};
         };
         static JsonLabels jsonLabels;
     };
@@ -103,9 +94,7 @@ private:
     //  Eigen::Tensor<double, 4> theta3Derivatives;
     //  std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
     bool shouldApplyInitialDistortion = false;
-    //#ifdef USE_ENSMALLEN
+
-    //    std::shared_ptr<SimulationJob> pJob;
-    //#endif
     void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
     void updateNodalInternalForces(
         const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
@@ -153,7 +142,7 @@ private:
         const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
 
 #ifdef POLYSCOPE_DEFINED
-    void draw(const std::string &screenshotsFolder = {});
+  void draw(const string &screenshotsFolder= {});
 #endif
   void
   updateNodalInternalForce(Vector6d &nodalInternalForce,
@@ -252,19 +241,11 @@ private:
 
   virtual SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob);
 
-  void reset(const std::shared_ptr<SimulationJob> &pJob);
-
   public:
   DRMSimulationModel();
   SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
                                       const Settings &settings,
                                       const SimulationResults &solutionGuess = SimulationResults());
-#ifdef USE_ENSMALLEN
-  double EvaluateWithGradient(const arma::mat &x, arma::mat &g);
-  void setJob(const std::shared_ptr<SimulationJob> &pJob);
-  SimulationMesh *getDeformedMesh(const arma::mat &x);
-
-#endif
 
   static void runUnitTests();
 };

edgemesh.cpp

@@ -25,13 +25,13 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
   return eigenEdgeNormals;
 }
 
-bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
+bool VCGEdgeMesh::save(const string &plyFilename)
 {
-    std::string filename = meshFilePath;
+    std::string filename = plyFilename;
     if (filename.empty()) {
         filename = std::filesystem::current_path().append(getLabel() + ".ply").string();
-    } else if (std::filesystem::is_directory(std::filesystem::path(meshFilePath))) {
+    } else if (std::filesystem::is_directory(std::filesystem::path(plyFilename))) {
-        filename = std::filesystem::path(meshFilePath).append(getLabel() + ".ply").string();
+        filename = std::filesystem::path(plyFilename).append(getLabel() + ".ply").string();
     }
     assert(std::filesystem::path(filename).extension().string() == ".ply");
     unsigned int mask = 0;
@@ -197,14 +197,14 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
   return true;
 }
 
-bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
+bool VCGEdgeMesh::load(const string &plyFilename) {
-{
+
-    std::string usedPath = meshFilePath;
+  std::string usedPath = plyFilename;
-    if (std::filesystem::path(meshFilePath).is_relative()) {
+  if (std::filesystem::path(plyFilename).is_relative()) {
-        usedPath = std::filesystem::absolute(meshFilePath).string();
+    usedPath = std::filesystem::absolute(plyFilename).string();
   }
   assert(std::filesystem::exists(usedPath));
-    Clear();
+  this->Clear();
   //  if (!loadUsingNanoply(usedPath)) {
   //    std::cerr << "Error: Unable to open " + usedPath << std::endl;
   //    return false;
@@ -217,15 +217,15 @@ bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
   getEdges(eigenEdges);
   getVertices(eigenVertices);
   vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
-    label = std::filesystem::path(meshFilePath).stem().string();
+  label = std::filesystem::path(plyFilename).stem().string();
 
   const bool printDebugInfo = false;
   if (printDebugInfo) {
-        std::cout << meshFilePath << " was loaded successfuly." << std::endl;
+    std::cout << plyFilename << " was loaded successfuly." << std::endl;
     std::cout << "Mesh has " << EN() << " edges." << std::endl;
   }
 
-    label = std::filesystem::path(meshFilePath).stem().string();
+  label=std::filesystem::path(plyFilename).stem().string();
   return true;
 }
 
@@ -248,7 +248,7 @@ bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
 
 bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename)
 {
-    Clear();
+    this->Clear();
     //  assert(plyFileHasAllRequiredFields(plyFilename));
     // Load the ply file
     int mask = 0;
@@ -337,35 +337,10 @@ bool VCGEdgeMesh::copy(VCGEdgeMesh &mesh) {
   return true;
 }
 
-void VCGEdgeMesh::set(const std::vector<double> &vertexPositions, const std::vector<int> &edges)
-{
-    Clear();
-    for (int ei = 0; ei < edges.size(); ei += 2) {
-        const int vi0 = edges[ei];
-        const int vi1 = edges[ei + 1];
-        const int vi0_offset = 3 * vi0;
-        const int vi1_offset = 3 * vi1;
-        const CoordType p0(vertexPositions[vi0_offset],
-                           vertexPositions[vi0_offset + 1],
-                           vertexPositions[vi0_offset + 2]);
-        const CoordType p1(vertexPositions[vi1_offset],
-                           vertexPositions[vi1_offset + 1],
-                           vertexPositions[vi1_offset + 2]);
-        auto eIt = vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(*this, p0, p1);
-        CoordType n(0, 0, 1);
-        eIt->cV(0)->N() = n;
-        eIt->cV(1)->N() = n;
-    }
-    removeDuplicateVertices();
-
-    updateEigenEdgeAndVertices();
-}
-
 void VCGEdgeMesh::removeDuplicateVertices()
 {
     vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.000061524);
     vcg::tri::Allocator<VCGEdgeMesh>::CompactVertexVector(*this);
-    vcg::tri::Allocator<VCGEdgeMesh>::CompactEdgeVector(*this);
     vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
 }
 
@@ -443,12 +418,8 @@ polyscope::CurveNetwork *VCGEdgeMesh::registerForDrawing(
 {
     PolyscopeInterface::init();
     const double drawingRadius = desiredRadius;
-    updateEigenEdgeAndVertices();
     polyscope::CurveNetwork *polyscopeHandle_edgeMesh
         = polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges());
-    //    std::cout << "EDGES:" << polyscopeHandle_edgeMesh->nEdges() << std::endl;
-    assert(polyscopeHandle_edgeMesh->nEdges() == getEigenEdges().rows()
-           && polyscopeHandle_edgeMesh->nNodes() == getEigenVertices().rows());
 
     polyscopeHandle_edgeMesh->setEnabled(shouldEnable);
     polyscopeHandle_edgeMesh->setRadius(drawingRadius, true);

edgemesh.hpp

@@ -58,12 +58,9 @@ public:
              * */
     bool copy(VCGEdgeMesh &mesh);
 
-    void set(const std::vector<double> &vertexPositions, const std::vector<int> &edges);
-
     void removeDuplicateVertices();
-    virtual void deleteDanglingVertices();
+    void deleteDanglingVertices();
-    virtual void deleteDanglingVertices(
+    void deleteDanglingVertices(vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
-        vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
 
     void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;
 
@@ -73,8 +70,8 @@ public:
 
     //    bool loadUsingNanoply(const std::string &plyFilename);
 
-    bool load(const std::filesystem::path &meshFilePath) override;
+    bool load(const std::string &plyFilename) override;
-    bool save(const std::filesystem::path &meshFilePath = std::filesystem::path()) override;
+    bool save(const std::string &plyFilename = std::string()) override;
 
     bool createSpanGrid(const size_t squareGridDimension);
     bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight);

hexagonremesher.hpp

@@ -11,7 +11,7 @@
 
 namespace PolygonalRemeshing {
 //std::shared_ptr<VCGPolyMesh>
-inline std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
+std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
 {
     //    const std::string tileIntoFilePath(
     //        "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply");

linearsimulationmodel.cpp

@@ -1,274 +0,0 @@
-#include "linearsimulationmodel.hpp"
-//#include "gsl/gsl"
-
-std::vector<fea::Elem> LinearSimulationModel::getFeaElements(
-    const std::shared_ptr<SimulationMesh> &pMesh)
-{
-    const int numberOfEdges = pMesh->EN();
-    std::vector<fea::Elem> elements(numberOfEdges);
-    for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
-        const SimulationMesh::CoordType &evn0 = pMesh->edge[edgeIndex].cV(0)->cN();
-        const SimulationMesh::CoordType &evn1 = pMesh->edge[edgeIndex].cV(1)->cN();
-        const std::vector<double> nAverageVector{(evn0[0] + evn1[0]) / 2,
-                                                 (evn0[1] + evn1[1]) / 2,
-                                                 (evn0[2] + evn1[2]) / 2};
-        const Element &element = pMesh->elements[edgeIndex];
-        const double E = element.material.youngsModulus;
-        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));
-        elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
-    }
-
-    return elements;
-}
-
-std::vector<fea::Node> LinearSimulationModel::getFeaNodes(
-    const std::shared_ptr<SimulationMesh> &pMesh)
-{
-    const int numberOfNodes = pMesh->VN();
-    std::vector<fea::Node> feaNodes(numberOfNodes);
-    for (int vi = 0; vi < numberOfNodes; vi++) {
-        const CoordType &p = pMesh->vert[vi].cP();
-        feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
-    }
-
-    return feaNodes;
-}
-
-std::vector<fea::BC> LinearSimulationModel::getFeaFixedNodes(
-    const std::shared_ptr<SimulationJob> &job)
-{
-    std::vector<fea::BC> boundaryConditions;
-    //        boundaryConditions.reserve(job->constrainedVertices.size() * 6
-    //                                   + job->nodalForcedDisplacements.size() * 3);
-    for (auto fixedVertex : job->constrainedVertices) {
-        const int vertexIndex = fixedVertex.first;
-        for (int dofIndex : fixedVertex.second) {
-            boundaryConditions.emplace_back(
-                fea::BC(vertexIndex, static_cast<fea::DOF>(dofIndex), 0));
-        }
-    }
-
-    for (auto forcedDisplacement : job->nodalForcedDisplacements) {
-        const int vi = forcedDisplacement.first;
-        for (int dofIndex = 0; dofIndex < 3; dofIndex++) {
-            boundaryConditions.emplace_back(
-                fea::BC(vi, static_cast<fea::DOF>(dofIndex), forcedDisplacement.second[dofIndex]));
-        }
-    }
-
-    return boundaryConditions;
-}
-
-std::vector<fea::Force> LinearSimulationModel::getFeaNodalForces(
-    const std::shared_ptr<SimulationJob> &job)
-{
-    std::vector<fea::Force> nodalForces;
-    nodalForces.reserve(job->nodalExternalForces.size() * 6);
-
-    for (auto nodalForce : job->nodalExternalForces) {
-        for (int dofIndex = 0; dofIndex < 6; dofIndex++) {
-            if (nodalForce.second[dofIndex] == 0) {
-                continue;
-            }
-            fea::Force f(nodalForce.first, dofIndex, nodalForce.second[dofIndex]);
-            nodalForces.emplace_back(f);
-        }
-    }
-
-    return nodalForces;
-}
-
-SimulationResults LinearSimulationModel::getResults(
-    const fea::Summary &feaSummary, const std::shared_ptr<SimulationJob> &simulationJob)
-{
-    SimulationResults results;
-    results.converged = feaSummary.converged;
-    if (!results.converged) {
-        return results;
-    }
-
-    results.executionTime = feaSummary.total_time_in_ms * 1000;
-    // displacements
-    results.displacements.resize(feaSummary.num_nodes);
-    results.rotationalDisplacementQuaternion.resize(feaSummary.num_nodes);
-    for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
-        results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
-
-        const Vector6d &nodalDisplacement = results.displacements[vi];
-        Eigen::Quaternion<double> q_nx;
-        q_nx = Eigen::AngleAxis<double>(nodalDisplacement[3], Eigen::Vector3d(1, 0, 0));
-        Eigen::Quaternion<double> q_ny;
-        q_ny = Eigen::AngleAxis<double>(nodalDisplacement[4], Eigen::Vector3d(0, 1, 0));
-        Eigen::Quaternion<double> q_nz;
-        q_nz = Eigen::AngleAxis<double>(nodalDisplacement[5], Eigen::Vector3d(0, 0, 1));
-        results.rotationalDisplacementQuaternion[vi] = q_nx * q_ny * q_nz;
-    }
-    results.setLabelPrefix("Linear");
-
-    //  // Convert forces
-    //  // Convert to vector of eigen matrices of the form force component-> per
-    //  // Edge
-    //  const int numDof = 6;
-    //  const size_t numberOfEdges = feaSummary.element_forces.size();
-    //  edgeForces =
-    //      std::vector<Eigen::VectorXd>(numDof, Eigen::VectorXd(2 *
-    //      numberOfEdges));
-    //  for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
-    //    for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof;
-    //         forceComponentIndex++) {
-    //      (edgeForces[forceComponentIndex])(2 * edgeIndex) =
-    //          feaSummary.element_forces[edgeIndex][forceComponentIndex];
-    //      (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) =
-    //          feaSummary.element_forces[edgeIndex][numDof +
-    //          forceComponentIndex];
-    //    }
-    //  }
-
-    for (int ei = 0; ei < feaSummary.num_elems; ei++) {
-        const std::vector<double> &elementForces = feaSummary.element_forces[ei];
-        const Element &element = simulationJob->pMesh->elements[ei];
-        //Axial
-        const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
-                                                       * element.initialLength
-                                                       / (4 * element.material.youngsModulus
-                                                          * element.dimensions.A);
-        const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
-                                                       * element.initialLength
-                                                       / (4 * element.material.youngsModulus
-                                                          * 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.dimensions.A
-                                                           * element.material.G);
-        const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
-                                                        * element.initialLength
-                                                        / (4 * element.dimensions.A
-                                                           * element.material.G);
-        const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
-                                                        * element.initialLength
-                                                        / (4 * element.dimensions.A
-                                                           * element.material.G);
-        const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
-                                                        * element.initialLength
-                                                        / (4 * element.dimensions.A
-                                                           * element.material.G);
-        const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
-                                                    + elementPotentialEnergy_shearZ_v0
-                                                    + elementPotentialEnergy_shearY_v1
-                                                    + elementPotentialEnergy_shearZ_v1;
-        //Bending
-        const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
-                                                          * element.initialLength
-                                                          / (4 * element.material.youngsModulus
-                                                             * element.dimensions.inertia.I2);
-        const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
-                                                          * element.initialLength
-                                                          / (4 * element.material.youngsModulus
-                                                             * element.dimensions.inertia.I3);
-        const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
-                                                          * element.initialLength
-                                                          / (4 * element.material.youngsModulus
-                                                             * element.dimensions.inertia.I2);
-        const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
-                                                          * element.initialLength
-                                                          / (4 * element.material.youngsModulus
-                                                             * element.dimensions.inertia.I3);
-        const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
-                                                      + elementPotentialEnergy_bendingZ_v0
-                                                      + elementPotentialEnergy_bendingY_v1
-                                                      + elementPotentialEnergy_bendingZ_v1;
-        //Torsion
-        const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
-                                                         * element.initialLength
-                                                         / (4 * element.dimensions.inertia.J
-                                                            * element.material.G);
-        const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
-                                                         * element.initialLength
-                                                         / (4 * element.dimensions.inertia.J
-                                                            * element.material.G);
-        const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
-                                                      + elementPotentialEnergy_torsion_v1;
-        const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
-                                                      + elementPotentialEnergy_shear
-                                                      + elementPotentialEnergy_bending
-                                                      + elementPotentialEnergy_torsion;
-        results.internalPotentialEnergy += elementInternalPotentialEnergy;
-    }
-    results.pJob = simulationJob;
-
-    return results;
-}
-
-SimulationResults LinearSimulationModel::executeSimulation(
-    const std::shared_ptr<SimulationJob> &pSimulationJob)
-{
-    assert(pSimulationJob->pMesh->VN() != 0);
-    const bool wasInitializedWithDifferentStructure = pSimulationJob->pMesh->EN()
-                                                          != simulator.structure.elems.size()
-                                                      || pSimulationJob->pMesh->VN()
-                                                             != simulator.structure.nodes.size();
-    if (!simulator.wasInitialized() || wasInitializedWithDifferentStructure) {
-        setStructure(pSimulationJob->pMesh);
-    }
-    //  printInfo(job);
-
-    // create the default options
-    fea::Options opts;
-    opts.save_elemental_forces = false;
-    opts.save_nodal_displacements = false;
-    opts.save_nodal_forces = false;
-    opts.save_report = false;
-    opts.save_tie_forces = false;
-    //  if (!elementalForcesOutputFilepath.empty()) {
-    //    opts.save_elemental_forces = true;
-    //    opts.elemental_forces_filename = elementalForcesOutputFilepath;
-    //  }
-    //  if (!nodalDisplacementsOutputFilepath.empty()) {
-    //    opts.save_nodal_displacements = true;
-    //    opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath;
-    //  }
-
-    // have the program output status updates
-    opts.verbose = false;
-
-    // form an empty vector of ties
-    std::vector<fea::Tie> ties;
-
-    // also create an empty list of equations
-    std::vector<fea::Equation> equations;
-    auto fixedVertices = getFeaFixedNodes(pSimulationJob);
-    auto nodalForces = getFeaNodalForces(pSimulationJob);
-    fea::Summary feaResults = simulator.solve(fixedVertices, nodalForces, ties, equations, opts);
-
-    SimulationResults results = getResults(feaResults, pSimulationJob);
-    results.pJob = pSimulationJob;
-    return results;
-}
-
-void LinearSimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
-{
-    fea::BeamStructure structure(getFeaNodes(pMesh), getFeaElements(pMesh));
-    simulator.initialize(structure);
-}
-
-void LinearSimulationModel::printInfo(const fea::BeamStructure &job)
-{
-    std::cout << "Details regarding the fea::Job:" << std::endl;
-    std::cout << "Nodes:" << std::endl;
-    for (fea::Node n : job.nodes) {
-        std::cout << n << std::endl;
-    }
-    std::cout << "Elements:" << std::endl;
-    for (Eigen::Vector2i e : job.elems) {
-        std::cout << e << std::endl;
-    }
-}

linearsimulationmodel.hpp

@@ -12,21 +12,265 @@ public:
     LinearSimulationModel(){
 
     }
-    static std::vector<fea::Elem> getFeaElements(const std::shared_ptr<SimulationMesh> &pMesh);
+    static std::vector<fea::Elem>
+    getFeaElements(const std::shared_ptr<SimulationJob> &job) {
+        const int numberOfEdges = job->pMesh->EN();
+        std::vector<fea::Elem> elements(numberOfEdges);
+        for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
+                const SimulationMesh::CoordType &evn0 =
+                        job->pMesh->edge[edgeIndex].cV(0)->cN();
+                const SimulationMesh::CoordType &evn1 = job->pMesh->edge[edgeIndex].cV(1)->cN();
+                const std::vector<double> nAverageVector{(evn0[0] + evn1[0]) / 2,
+                                                         (evn0[1] + evn1[1]) / 2,
+                                                         (evn0[2] + evn1[2]) / 2};
+                const Element &element = job->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,
+                                         nAverageVector);
+                const int vi0 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(0));
+                const int vi1 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(1));
+                elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
+            }
 
-    static std::vector<fea::Node> getFeaNodes(const std::shared_ptr<SimulationMesh> &pMesh);
+        return elements;
-    static std::vector<fea::BC> getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job);
+    }
 
-    static std::vector<fea::Force> getFeaNodalForces(const std::shared_ptr<SimulationJob> &job);
+    static std::vector<fea::Node>
+    getFeaNodes(const std::shared_ptr<SimulationJob> &job) {
+        const int numberOfNodes = job->pMesh->VN();
+        std::vector<fea::Node> feaNodes(numberOfNodes);
+        for (int vi = 0; vi < numberOfNodes; vi++) {
+                const CoordType &p = job->pMesh->vert[vi].cP();
+                feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
+            }
+
+        return feaNodes;
+    }
+    static std::vector<fea::BC>
+    getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job) {
+        std::vector<fea::BC> boundaryConditions;
+        //        boundaryConditions.reserve(job->constrainedVertices.size() * 6
+        //                                   + job->nodalForcedDisplacements.size() * 3);
+        for (auto fixedVertex : job->constrainedVertices) {
+                const int vertexIndex = fixedVertex.first;
+                for (int dofIndex : fixedVertex.second) {
+                    boundaryConditions.emplace_back(
+                        fea::BC(vertexIndex, static_cast<fea::DOF>(dofIndex), 0));
+                }
+        }
+
+        for (auto forcedDisplacement : job->nodalForcedDisplacements) {
+            const int vi = forcedDisplacement.first;
+            for (int dofIndex = 0; dofIndex < 3; dofIndex++) {
+                boundaryConditions.emplace_back(fea::BC(vi,
+                                                        static_cast<fea::DOF>(dofIndex),
+                                                        forcedDisplacement.second[dofIndex]));
+            }
+        }
+
+        return boundaryConditions;
+    }
+
+    static std::vector<fea::Force>
+    getFeaNodalForces(const std::shared_ptr<SimulationJob> &job) {
+        std::vector<fea::Force> nodalForces;
+        nodalForces.reserve(job->nodalExternalForces.size() * 6);
+
+        for (auto nodalForce : job->nodalExternalForces) {
+                for (int dofIndex = 0; dofIndex < 6; dofIndex++) {
+                        if (nodalForce.second[dofIndex] == 0) {
+                                continue;
+                            }
+                        fea::Force f(nodalForce.first, dofIndex, nodalForce.second[dofIndex]);
+                        nodalForces.emplace_back(f);
+                    }
+            }
+
+        return nodalForces;
+    }
     static SimulationResults getResults(const fea::Summary &feaSummary,
-                                        const std::shared_ptr<SimulationJob> &simulationJob);
+                                        const std::shared_ptr<SimulationJob> &simulationJob)
+    {
+        SimulationResults results;
+        results.converged = feaSummary.converged;
+        if (!results.converged) {
+            return results;
+        }
 
-    SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob);
+        results.executionTime = feaSummary.total_time_in_ms * 1000;
-    void setStructure(const std::shared_ptr<SimulationMesh> &pMesh);
+        // displacements
+        results.displacements.resize(feaSummary.num_nodes);
+        results.rotationalDisplacementQuaternion.resize(feaSummary.num_nodes);
+        for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
+            results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
+
+            const Vector6d &nodalDisplacement = results.displacements[vi];
+            Eigen::Quaternion<double> q_nx;
+            q_nx = Eigen::AngleAxis<double>(nodalDisplacement[3], Eigen::Vector3d(1, 0, 0));
+            Eigen::Quaternion<double> q_ny;
+            q_ny = Eigen::AngleAxis<double>(nodalDisplacement[4], Eigen::Vector3d(0, 1, 0));
+            Eigen::Quaternion<double> q_nz;
+            q_nz = Eigen::AngleAxis<double>(nodalDisplacement[5], Eigen::Vector3d(0, 0, 1));
+            results.rotationalDisplacementQuaternion[vi] = q_nx * q_ny * q_nz;
+            //            results.rotationalDisplacementQuaternion[vi] = q_nz * q_ny * q_nx;
+            //            results.rotationalDisplacementQuaternion[vi] = q_nz * q_nx * q_ny;
+        }
+        results.setLabelPrefix("Linear");
+
+        //  // Convert forces
+        //  // Convert to vector of eigen matrices of the form force component-> per
+        //  // Edge
+        //  const int numDof = 6;
+        //  const size_t numberOfEdges = feaSummary.element_forces.size();
+        //  edgeForces =
+        //      std::vector<Eigen::VectorXd>(numDof, Eigen::VectorXd(2 *
+        //      numberOfEdges));
+        //  for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
+        //    for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof;
+        //         forceComponentIndex++) {
+        //      (edgeForces[forceComponentIndex])(2 * edgeIndex) =
+        //          feaSummary.element_forces[edgeIndex][forceComponentIndex];
+        //      (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) =
+        //          feaSummary.element_forces[edgeIndex][numDof +
+        //          forceComponentIndex];
+        //    }
+        //  }
+
+        for (int ei = 0; ei < feaSummary.num_elems; ei++) {
+            const std::vector<double> &elementForces = feaSummary.element_forces[ei];
+            const Element &element = simulationJob->pMesh->elements[ei];
+            //Axial
+            const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
+                                                           * element.initialLength
+                                                           / (4 * element.material.youngsModulus
+                                                              * element.A);
+            const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
+                                                           * element.initialLength
+                                                           / (4 * element.material.youngsModulus
+                                                              * element.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);
+            const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
+                                                            * element.initialLength
+                                                            / (4 * element.A * element.G);
+            const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
+                                                            * element.initialLength
+                                                            / (4 * element.A * element.G);
+            const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
+                                                            * element.initialLength
+                                                            / (4 * element.A * element.G);
+            const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
+                                                        + elementPotentialEnergy_shearZ_v0
+                                                        + elementPotentialEnergy_shearY_v1
+                                                        + elementPotentialEnergy_shearZ_v1;
+            //Bending
+            const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
+                                                              * element.initialLength
+                                                              / (4 * element.material.youngsModulus
+                                                                 * element.inertia.I2);
+            const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
+                                                              * element.initialLength
+                                                              / (4 * element.material.youngsModulus
+                                                                 * element.inertia.I3);
+            const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
+                                                              * element.initialLength
+                                                              / (4 * element.material.youngsModulus
+                                                                 * element.inertia.I2);
+            const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
+                                                              * element.initialLength
+                                                              / (4 * element.material.youngsModulus
+                                                                 * element.inertia.I3);
+            const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
+                                                          + elementPotentialEnergy_bendingZ_v0
+                                                          + elementPotentialEnergy_bendingY_v1
+                                                          + elementPotentialEnergy_bendingZ_v1;
+            //Torsion
+            const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
+                                                             * element.initialLength
+                                                             / (4 * element.inertia.J * element.G);
+            const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
+                                                             * element.initialLength
+                                                             / (4 * element.inertia.J * element.G);
+            const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
+                                                          + elementPotentialEnergy_torsion_v1;
+            const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
+                                                          + elementPotentialEnergy_shear
+                                                          + elementPotentialEnergy_bending
+                                                          + elementPotentialEnergy_torsion;
+            results.internalPotentialEnergy += elementInternalPotentialEnergy;
+        }
+        results.pJob = simulationJob;
+
+        return results;
+    }
+
+  SimulationResults
+  executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob) {
+      assert(simulationJob->pMesh->VN() != 0);
+      fea::Job job(getFeaNodes(simulationJob), getFeaElements(simulationJob));
+      //  printInfo(job);
+
+      // create the default options
+      fea::Options opts;
+      opts.save_elemental_forces = false;
+      opts.save_nodal_displacements = false;
+      opts.save_nodal_forces = false;
+      opts.save_report = false;
+      opts.save_tie_forces = false;
+      //  if (!elementalForcesOutputFilepath.empty()) {
+      //    opts.save_elemental_forces = true;
+      //    opts.elemental_forces_filename = elementalForcesOutputFilepath;
+      //  }
+      //  if (!nodalDisplacementsOutputFilepath.empty()) {
+      //    opts.save_nodal_displacements = true;
+      //    opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath;
+      //  }
+
+      // have the program output status updates
+      opts.verbose = false;
+
+      // form an empty vector of ties
+      std::vector<fea::Tie> ties;
+
+      // also create an empty list of equations
+      std::vector<fea::Equation> equations;
+      auto fixedVertices = getFeaFixedNodes(simulationJob);
+      auto nodalForces = getFeaNodalForces(simulationJob);
+      fea::Summary feaResults =
+              fea::solve(job, fixedVertices, nodalForces, ties, equations, opts);
+
+      SimulationResults results = getResults(feaResults, simulationJob);
+      results.pJob = simulationJob;
+      return results;
+  }
+  //  SimulationResults getResults() const;
+  //  void setResultsNodalDisplacementCSVFilepath(const std::string
+  //  &outputPath); void setResultsElementalForcesCSVFilepath(const std::string
+  //  &outputPath);
 
 private:
-    fea::ThreedBeamFEA simulator;
+  //  std::string nodalDisplacementsOutputFilepath{"nodal_displacement.csv"};
-    static void printInfo(const fea::BeamStructure &job);
+  //  std::string elementalForcesOutputFilepath{"elemental_forces.csv"};
+  //  SimulationResults results;
+
+  static void printInfo(const fea::Job &job) {
+      std::cout << "Details regarding the fea::Job:" << std::endl;
+      std::cout << "Nodes:" << std::endl;
+      for (fea::Node n : job.nodes) {
+              std::cout << n << std::endl;
+          }
+      std::cout << "Elements:" << std::endl;
+      for (Eigen::Vector2i e : job.elems) {
+              std::cout << e << std::endl;
+          }
+  }
 };
 
 #endif // LINEARSIMULATIONMODEL_HPP

mesh.hpp

@@ -1,7 +1,7 @@
 #ifndef MESH_HPP
 #define MESH_HPP
 
-#include <filesystem>
+#include <Eigen/Core>
 #include <string>
 
 class Mesh
@@ -11,8 +11,8 @@ protected:
 
 public:
   virtual ~Mesh() = default;
-  virtual bool load(const std::filesystem::path &meshFilePath) { return false; }
+  virtual bool load(const std::string &filePath) { return false; }
-  virtual bool save(const std::filesystem::path &meshFilePath) { return false; }
+  virtual bool save(const std::string &filePath) { return false; }
 
   std::string getLabel() const;
   void setLabel(const std::string &newLabel);

polymesh.hpp

@@ -5,7 +5,6 @@
 #include <filesystem>
 #include <wrap/io_trimesh/export_obj.h>
 #include <wrap/io_trimesh/export_ply.h>
-#include <wrap/io_trimesh/import.h>
 
 #ifdef POLYSCOPE_DEFINED
 #include <polyscope/surface_mesh.h>
@@ -65,25 +64,14 @@ class VCGPolyMesh
       public Mesh
 {
 public:
-    virtual bool load(const std::filesystem::__cxx11::path &meshFilePath) override
+    //  virtual bool load(const std::string &filename) override {
-    {
+    //    int mask;
-        int mask;
+    //    vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(filename.c_str(), mask);
-        vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(meshFilePath.c_str(), mask);
+    //    int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(
-        int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, meshFilePath.c_str(), mask);
+    //        *this, filename.c_str(), mask);
-        if (error != 0) {
+    //    if (error != 0) {
-            std::cerr << "Could not load polygonal mesh:" << meshFilePath << std::endl;
+    //      return false;
-            return false;
+    //    }
-        }
-
-        vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
-        vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
-        vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
-        vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
-        vcg::tri::Clean<VCGPolyMesh>::RemoveUnreferencedVertex(*this);
-        //finally remove valence 1 vertices on the border
-        //        vcg::PolygonalAlgorithm<PolyMeshType>::RemoveValence2Vertices(dual);
-        return true;
-    }
     //    //    vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open();
     //    //    unsigned int mask = 0;
     //    //    mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOORD;
@@ -96,6 +84,11 @@ public:
     //    //            0) {
     //    //      std::cout << "Could not load tri mesh" << std::endl;
     //    //    }
+    //    vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
+    //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
+    //    //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
+    //    vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
+    //    return true;
     //  }
     Eigen::MatrixX3d getVertices() const
     {
@@ -120,34 +113,6 @@ public:
       return faces;
   }
 
-  //  bool load(const std::filesystem::__cxx11::path &meshFilePath)
-  //  {
-  //      const std::string extension = ".ply";
-  //      std::filesystem::path filePath = meshFilePath;
-  //      assert(std::filesystem::path(filePath).extension().string() == extension);
-  //      unsigned int mask = 0;
-  //      mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
-  //      mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
-  //      mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
-  //      mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
-  //      if (vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, filePath.c_str()) != 0) {
-  //          return false;
-  //      }
-  //      label = meshFilePath.filename();
-  //      return true;
-  //  }
-
-  bool save(const std::filesystem::__cxx11::path &meshFilePath = std::filesystem::path())
-  {
-      if (meshFilePath.extension() == ".obj") {
-          return saveOBJ(meshFilePath);
-      } else if (meshFilePath.extension() == ".ply") {
-          return savePLY(meshFilePath);
-      }
-
-      return false;
-  }
-
   bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path())
   {
       const std::string extension = ".obj";
@@ -163,7 +128,7 @@ public:
       mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
       mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
       mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
-      if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask) != 0) {
+      if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.c_str(), mask) != 0) {
           return false;
       }
       return true;
@@ -184,7 +149,7 @@ public:
       mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
       mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
       mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
-      if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask, false) != 0) {
+      if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.c_str(), mask, false) != 0) {
           return false;
       }
       return true;

reducedmodelevaluator.cpp

@@ -1,317 +0,0 @@
-#include "reducedmodelevaluator.hpp"
-#include "hexagonremesher.hpp"
-#include "trianglepatterngeometry.hpp"
-#include <filesystem>
-
-using FullPatternVertexIndex = VertexIndex;
-using ReducedPatternVertexIndex = VertexIndex;
-
-ReducedModelEvaluator::ReducedModelEvaluator()
-{
-
-}
-
-std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
-		ReducedModelOptimization::Results &optimizationResults)
-{
-	//    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
-	//    std::filesystem::path tileIntoSurfaceFilePath{
-	//        "/home/iason/Coding/Projects/Approximating shapes with flat "
-	//        "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
-	//    assert(std::filesystem::exists(tileIntoSurfaceFilePath));
-	//    pTileIntoSurface->load(tileIntoSurfaceFilePath);
-
-	//Set required file paths
-	const std::filesystem::path tileInto_triMesh_filename
-		= "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
-		"instantMeshes_plane_34.ply";
-	const std::filesystem::path reducedPatternFilePath
-		= "/home/iason/Coding/Projects/Approximating shapes with flat "
-		"patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
-	const std::filesystem::path intermediateResultsDirectoryPath
-		= "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults";
-	//    const std::filesystem::path drmSettingsFilePath
-	//        = "/home/iason/Coding/Projects/Approximating shapes with flat "
-	//          "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
-	//          "defaultDRMSimulationSettings.json";
-	DRMSimulationModel::Settings drmSimulationSettings;
-	//    drmSimulationSettings.load(drmSettingsFilePath);
-	drmSimulationSettings.linearGuessForceScaleFactor = 0.7;
-	drmSimulationSettings.debugModeStep = 10000;
-	drmSimulationSettings.beVerbose = true;
-	constexpr bool shouldRerunFullPatternSimulation = false;
-	const std::vector<std::string> scenariosTestSetLabels{"22Hex_randomBending0",
-		"22Hex_randomBending1",
-		"22Hex_randomBending2",
-		"22Hex_randomBending3",
-		"22Hex_randomBending4",
-		"22Hex_randomBending5",
-		"22Hex_randomBending6",
-		"22Hex_randomBending7",
-		"22Hex_randomBending8",
-		"22Hex_randomBending9",
-		"22Hex_randomBending10",
-		"22Hex_randomBending11",
-		"22Hex_randomBending12",
-		"22Hex_randomBending13",
-		"22Hex_randomBending14",
-		"22Hex_randomBending15",
-		"22Hex_randomBending16",
-		"22Hex_randomBending17",
-		"22Hex_randomBending18",
-		"22Hex_randomBending19",
-		"22Hex_randomBending20",
-		"22Hex_bending_0.005N",
-		"22Hex_bending_0.01N",
-		"22Hex_bending_0.03N",
-		"22Hex_bending_0.05N",
-		"22Hex_pullOppositeVerts_0.05N",
-		"22Hex_pullOppositeVerts_0.1N",
-		"22Hex_pullOppositeVerts_0.3N",
-		"22Hex_shear_2N",
-		"22Hex_shear_5N",
-		"22Hex_axial_10N",
-		"22Hex_axial_20N",
-		"22Hex_cylinder_0.05N",
-		"22Hex_cylinder_0.1N",
-		"22Hex_s_0.05N",
-		"22Hex_s_0.1N"};
-
-	//Load surface
-	std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
-		VCGTriMesh tileInto_triMesh;
-		const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
-		assert(surfaceLoadSuccessfull);
-		return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
-	}();
-	const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0),
-			(optimizationResults.baseTriangle.cP(1)
-			 + optimizationResults.baseTriangle.cP(
-				 2))
-			/ 2);
-	pTileIntoSurface->moveToCenter();
-	const double scaleFactor = optimizedBaseTriangleHeight
-		/ pTileIntoSurface->getAverageFaceRadius();
-	vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
-	//    tileIntoSurface.registerForDrawing();
-	//    polyscope::show();
-
-	//Tile full pattern into surface
-	std::vector<PatternGeometry> fullPatterns(1);
-	fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
-	////        Base triangle pattern might contain dangling vertices.Remove those
-	fullPatterns[0].interfaceNodeIndex = 3;
-	fullPatterns[0].deleteDanglingVertices();
-	std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
-	std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
-	std::vector<size_t> tileIntoEdgeToTiledFullVi;
-	std::shared_ptr<PatternGeometry> pTiledFullPattern
-		= PatternGeometry::tilePattern(fullPatterns,
-				{},
-				*pTileIntoSurface,
-				perSurfaceFacePatternIndices,
-				tileIntoEdgeToTiledFullVi,
-				perPatternIndexTiledFullPatternEdgeIndices);
-	pTiledFullPattern->setLabel("Tiled_full_patterns");
-	//    pTiledFullPattern->registerForDrawing();
-	//Tile reduced pattern into surface
-	PatternGeometry reducedPattern;
-	reducedPattern.load(reducedPatternFilePath);
-	reducedPattern.deleteDanglingVertices();
-	assert(reducedPattern.interfaceNodeIndex == 1);
-	std::vector<PatternGeometry> reducedPatterns(1);
-	reducedPatterns[0].copy(reducedPattern);
-	const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
-	ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
-			optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]);
-
-	std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
-	std::vector<size_t> tileIntoEdgeToTiledReducedVi;
-	std::shared_ptr<PatternGeometry> pTiledReducedPattern
-		= PatternGeometry::tilePattern(reducedPatterns,
-				{0},
-				*pTileIntoSurface,
-				perSurfaceFacePatternIndices,
-				tileIntoEdgeToTiledReducedVi,
-				perPatternIndexTiledReducedPatternEdgeIndices);
-	pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
-	//        pTiledReducedPattern->registerForDrawing();
-
-	std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
-		fullToReducedViMap; //of only the common vertices
-	std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
-		reducedToFullViMap; //of only the common vertices
-	for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
-		fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
-	}
-	constructInverseMap(fullToReducedViMap, reducedToFullViMap);
-
-	std::vector<size_t> tilledFullPatternInterfaceVi;
-	tilledFullPatternInterfaceVi.clear();
-	tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
-	size_t viIndex = 0;
-	for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
-		tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
-	}
-
-	//Create simulation meshes
-	////Tessellated full pattern simulation mesh
-	std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
-	pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
-	//NOTE: Those should be derived from the optimization results instead of hardcoding them
-	pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
-	if (optimizationResults.fullPatternYoungsModulus == 0) {
-		std::cerr << "Full pattern's young modulus not found." << std::endl;
-		std::terminate();
-	}
-	pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
-			optimizationResults.fullPatternYoungsModulus);
-	pTiledFullPattern_simulationMesh->reset();
-
-	////Tessellated reduced pattern simulation mesh
-	std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
-	pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
-	const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
-		= perPatternIndexTiledReducedPatternEdgeIndices[0];
-	ReducedModelOptimization::Results::applyOptimizationResults_elements(
-			optimizationResults, pTiledReducedPattern_simulationMesh);
-	pTiledReducedPattern_simulationMesh->reset();
-	std::vector<double> distances_drm2reduced(scenariosTestSetLabels.size(), 0);
-	std::vector<double> distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0);
-	for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
-		const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
-		const std::filesystem::path scenariosDirectoryPath
-			= "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/";
-		const std::filesystem::path tiledReducedPatternJobFilePath
-			= std::filesystem::path(scenariosDirectoryPath)
-			.append(jobLabel)
-			.append("SimulationJobs")
-			.append("Reduced")
-			.append(SimulationJob::jsonDefaultFileName);
-		//set jobs
-		std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
-		pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
-		pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
-		pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
-		std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
-		pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
-		pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
-		pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
-		//        pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
-		//        pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
-		//    polyscope::show();
-		//Save reduced job
-		const std::filesystem::path tesellatedResultsFolderPath
-			= std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults");
-		const std::filesystem::path surfaceFolderPath = std::filesystem::path(
-				tesellatedResultsFolderPath)
-			.append(pTileIntoSurface->getLabel());
-		const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
-		const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
-			.append(scenarioLabel);
-		const std::filesystem::path reducedJobDirectoryPath
-			= std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
-		std::filesystem::create_directories(reducedJobDirectoryPath);
-		pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
-		//Run scenario
-
-		////Full
-
-		const std::string patternLabel = [&]() {
-			const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel();
-			const int numberOfOccurences = std::count_if(patternLabel.begin(),
-					patternLabel.end(),
-					[](char c) { return c == '#'; });
-			if (numberOfOccurences == 0) {
-				return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#"
-					+ optimizationResults.baseTriangleFullPattern.getLabel();
-			} else if (numberOfOccurences == 1) {
-				return optimizationResults.baseTriangleFullPattern.getLabel();
-			}
-		}();
-		const auto fullResultsFolderPath
-			= std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
-		if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
-			std::filesystem::remove_all(fullResultsFolderPath);
-		}
-
-		const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
-				scenarioDirectoryPath)
-			.append(patternLabel)
-			.append("SimulationJob");
-		SimulationResults simulationResults_tiledFullPattern_drm;
-		if (std::filesystem::exists(fullResultsFolderPath)) {
-			//Load full pattern results
-			assert(std::filesystem::exists(fullPatternJobFolderPath));
-			simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
-					fullPatternJobFolderPath);
-			simulationResults_tiledFullPattern_drm.converged = true;
-		} else {
-			//Full
-			std::cout << "Executing:" << jobLabel << std::endl;
-			DRMSimulationModel drmSimulationModel;
-			simulationResults_tiledFullPattern_drm
-				= drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
-			simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
-		}
-		std::filesystem::create_directories(fullResultsFolderPath);
-		simulationResults_tiledFullPattern_drm.save(
-				std::filesystem::path(scenarioDirectoryPath).append(patternLabel));
-		if (!simulationResults_tiledFullPattern_drm.converged) {
-			std::cerr << "Full pattern simulation failed." << std::endl;
-		}
-
-		LinearSimulationModel linearSimulationModel;
-		SimulationResults simulationResults_tiledReducedPattern
-			= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
-
-		//        simulationResults_tiledReducedPattern.registerForDrawing();
-		//        simulationResults_tiledFullPattern_drm.registerForDrawing();
-		//        polyscope::show();
-		//measure distance
-		const double distance_fullDrmToReduced
-			= simulationResults_tiledFullPattern_drm
-			.computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
-		double distance_fullSumOfAllVerts = 0;
-		for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
-			distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
-				.displacements[fullToReducedPair.first]
-				.getTranslation()
-				.norm();
-		}
-		const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
-			/ distance_fullSumOfAllVerts;
-		distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
-		distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced;
-	}
-	//#ifndef POLYSCOPE_DEFINED
-	//    return distances_drm2reduced;
-	//#else
-
-	//report distance
-	//        csvFile csv_results("", flse);
-	csvFile csv_results({}, false);
-	//    csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
-	//                            .append("results.csv")
-	//                            .string(),
-	//                        false);
-	csv_results /*<< "Job Label"*/
-		<< "drm2Reduced"
-		<< "norm_drm2Reduced";
-	csv_results << endrow;
-	double sumOfNormalizedFull2Reduced = 0;
-	for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
-		const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
-		const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex];
-		const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex];
-		csv_results /*<< jobLabel*/ << distance_fullDrmToReduced
-			<< distance_normalizedFullDrmToReduced;
-		csv_results << endrow;
-		sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
-	}
-
-	std::cout << "Average normalized error per scenario:"
-		<< sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
-	return distances_normalizedDrm2reduced;
-	//#endif
-}

reducedmodelevaluator.hpp

@@ -1,14 +0,0 @@
-#ifndef REDUCEDMODELEVALUATOR_HPP
-#define REDUCEDMODELEVALUATOR_HPP
-
-#include "reducedmodeloptimizer_structs.hpp"
-
-class ReducedModelEvaluator
-{
-public:
-    ReducedModelEvaluator();
-    static std::vector<double> evaluateReducedModel(
-        ReducedModelOptimization::Results &optimizationResults);
-};
-
-#endif // REDUCEDMODELEVALUATOR_HPP

reducedmodeloptimizer_structs.hpp

@@ -7,7 +7,7 @@
 #include "unordered_map"
 #include <string>
 
-namespace ReducedModelOptimization {
+namespace ReducedPatternOptimization {
 
 enum BaseSimulationScenario {
     Axial,
@@ -34,9 +34,9 @@ struct Colors
 
 struct xRange
 {
-    std::string label{};
+    std::string label;
-    double min{0};
+    double min;
-    double max{0};
+    double max;
 
     inline bool operator<(const xRange &other) { return label < other.label; }
     std::string toString() const
@@ -60,106 +60,32 @@ struct xRange
     {
         return label == xrange.label && min == xrange.min && max == xrange.max;
     }
-
-    std::tuple<std::string, double, double> toTuple() const
-    {
-        return std::make_tuple(label, min, max);
-    }
-
-    void set(const std::tuple<std::string, double, double> &inputTuple)
-    {
-        if (std::get<1>(inputTuple) > std::get<2>(inputTuple)) {
-            std::cerr
-                << "Invalid xRange tuple. Second argument(min) cant be smaller than the third(max)"
-                << std::endl;
-            std::terminate();
-            //            return;
-        }
-        std::tie(label, min, max) = inputTuple;
-    }
 };
 
-enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationVariables };
+  struct Settings
-
+  {
-inline int getParameterIndex(const std::string &s)
-{
-    if ("E" == s) {
-        return E;
-    } else if ("A" == s) {
-        return A;
-    } else if ("I2" == s) {
-        return I2;
-    } else if ("I3" == s) {
-        return I3;
-    } else if ("J" == s) {
-        return J;
-    } else if ("R" == s || "HexSize" == s) {
-        return R;
-    } else if ("Theta" == s || "HexAngle" == s) {
-        return Theta;
-    } else {
-        std::cerr << "Input is not recognized as a valid optimization variable index:" << s
-                  << std::endl;
-        return -1;
-    }
-}
-
-struct Settings
-{
-    inline static std::string defaultFilename{"OptimizationSettings.json"};
-    std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
-        {E, A, I2, I3, J, R, Theta}};
-    //        {A, I2, I3, J, R, Theta}};
-    std::optional<std::vector<double>>
-        optimizationVariablesGroupsWeights; //TODO:should be removed in the future if not a splitting strategy is used for the optimization
       enum NormalizationStrategy { NonNormalized, Epsilon };
-    inline static std::vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
+      std::vector<xRange> xRanges;
-    NormalizationStrategy normalizationStrategy{Epsilon};
+      inline static vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
-    std::array<xRange, NumberOfOptimizationVariables> variablesRanges{xRange{"E", 0.01, 100},
-                                                                      xRange{"A", 0.01, 100},
-                                                                      xRange{"I2", 0.01, 100},
-                                                                      xRange{"I3", 0.01, 100},
-                                                                      xRange{"J", 0.01, 100},
-                                                                      xRange{"R", 0.05, 0.95},
-                                                                      xRange{"Theta", -30, 30}};
       int numberOfFunctionCalls{100000};
-    double solverAccuracy{1e-3};
+      double solverAccuracy{1e-2};
-    double translationEpsilon{3e-4};
+      NormalizationStrategy normalizationStrategy{Epsilon};
-    //    double translationEpsilon{0};
+      double translationNormalizationParameter{0.003};
-    double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
+      double rotationNormalizationParameter{3};
+      bool splitGeometryMaterialOptimization{false};
       struct ObjectiveWeights
       {
-        double translational{1.2};
+          double translational{1};
-        double rotational{0.8};
+          double rotational{1};
-        bool operator==(const ObjectiveWeights &other) const
+      } objectiveWeights;
-        {
-            return this->translational == other.translational
-                   && this->rotational == other.rotational;
-        }
-    };
-    std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> perBaseScenarioObjectiveWeights;
-    std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
-    convertObjectiveWeightsToPairs() const
-    {
-        std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
-        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
-             baseScenario++) {
-            objectiveWeightsPairs[baseScenario]
-                = std::make_pair(perBaseScenarioObjectiveWeights[baseScenario].translational,
-                                 perBaseScenarioObjectiveWeights[baseScenario].rotational);
-        }
-        return objectiveWeightsPairs;
-    }
 
       struct JsonKeys
       {
-        inline static std::string OptimizationStrategy{"OptimizationStrategy"};
+          inline static std::string filename{"OptimizationSettings.json"};
-        inline static std::string OptimizationStrategyGroupWeights{
-            "OptimizationStrategyGroupWeights"};
           inline static std::string OptimizationVariables{"OptimizationVariables"};
           inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"};
           inline static std::string SolverAccuracy{"SolverAccuracy"};
-        inline static std::string ObjectiveWeights{"ObjectiveWeight"};
+          inline static std::string TranslationalObjectiveWeight{"TransObjWeight"};
       };
 
       void save(const std::filesystem::path &saveToPath)
@@ -167,52 +93,21 @@ struct Settings
           assert(std::filesystem::is_directory(saveToPath));
 
           nlohmann::json json;
-        json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
-        if (optimizationVariablesGroupsWeights.has_value()) {
-            json[GET_VARIABLE_NAME(ptimizationStrategyGroupWeights)]
-                = optimizationVariablesGroupsWeights.value();
-        }
           //write x ranges
-        const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
+          for (size_t xRangeIndex = 0; xRangeIndex < xRanges.size(); xRangeIndex++) {
-            xRangesAsTuples = [=]() {
+              const auto &xRange = xRanges[xRangeIndex];
-                std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
+              json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
-                    xRangesAsTuples;
+                  = xRange.toString();
-                for (int optimizationParameterIndex = E;
-                     optimizationParameterIndex != NumberOfOptimizationVariables;
-                     optimizationParameterIndex++) {
-                    xRangesAsTuples[optimizationParameterIndex]
-                        = variablesRanges[optimizationParameterIndex].toTuple();
           }
-                return xRangesAsTuples;
-            }();
-        json[JsonKeys::OptimizationVariables] = xRangesAsTuples;
-        //        for (size_t xRangeIndex = 0; xRangeIndex < variablesRanges.size(); xRangeIndex++) {
-        //            const auto &xRange = variablesRanges[xRangeIndex];
-        //            json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
-        //                = xRange.toString();
-        //        }
 
-        json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
+          json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
-        json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
+          json[JsonKeys::SolverAccuracy] = solverAccuracy;
-        //Objective weights
+          json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational;
-        std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
-        std::transform(perBaseScenarioObjectiveWeights.begin(),
-                       perBaseScenarioObjectiveWeights.end(),
-                       objectiveWeightsPairs.begin(),
-                       [](const ObjectiveWeights &objectiveWeights) {
-                           return std::make_pair(objectiveWeights.translational,
-                                                 objectiveWeights.rotational);
-                       });
-        json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
-        json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
-        json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
-            angularDistanceEpsilon);
 
           std::filesystem::path jsonFilePath(
-            std::filesystem::path(saveToPath).append(defaultFilename));
+              std::filesystem::path(saveToPath).append(JsonKeys::filename));
           std::ofstream jsonFile(jsonFilePath.string());
           jsonFile << json;
-        jsonFile.close();
       }
 
       bool load(const std::filesystem::path &loadFromPath)
@@ -221,170 +116,95 @@ struct Settings
           //Load optimal X
           nlohmann::json json;
           std::filesystem::path jsonFilepath(
-            std::filesystem::path(loadFromPath).append(defaultFilename));
+              std::filesystem::path(loadFromPath).append(JsonKeys::filename));
           if (!std::filesystem::exists(jsonFilepath)) {
-            std::cerr << "Optimization settings could not be loaded because input filepath does "
-                         "not exist:"
-                      << jsonFilepath << std::endl;
               return false;
           }
           std::ifstream ifs(jsonFilepath);
           ifs >> json;
 
-        if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
-            optimizationStrategy
-                = std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>(
-                    (json.at(GET_VARIABLE_NAME(optimizationStrategy))));
-        }
-        if (json.contains(GET_VARIABLE_NAME(optimizationStrategyGroupWeights))) {
-            optimizationVariablesGroupsWeights = std::vector<double>(
-                json[GET_VARIABLE_NAME(optimizationStrategyGroupWeights)]);
-        }
           //read x ranges
-        if (json.contains(JsonKeys::OptimizationVariables)) {
-            const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
-                xRangesAsTuples = json.at(JsonKeys::OptimizationVariables);
-            for (const auto &rangeTuple : xRangesAsTuples) {
-                variablesRanges[getParameterIndex(std::get<0>(rangeTuple))].set(rangeTuple);
-            }
-        } else { //NOTE:legacy compatibility
           size_t xRangeIndex = 0;
           while (true) {
               const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_"
-                                                  + std::to_string(xRangeIndex++);
+                                                + std::to_string(xRangeIndex);
               if (!json.contains(jsonXRangeKey)) {
                   break;
               }
               xRange x;
               x.fromString(json.at(jsonXRangeKey));
-                variablesRanges[getParameterIndex(x.label)] = x;
+              xRanges.push_back(x);
-            }
+              xRangeIndex++;
           }
 
-        if (json.contains(GET_VARIABLE_NAME(numberOfFunctionCalls))) {
+          numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
-            numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(numberOfFunctionCalls));
+          solverAccuracy = json.at(JsonKeys::SolverAccuracy);
-        }
+          objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight);
-        if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
+          objectiveWeights.rotational = 2 - objectiveWeights.translational;
-            solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
-        }
-        //Objective weights
-        if (json.contains(JsonKeys::ObjectiveWeights)) {
-            std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
-                objectiveWeightsPairs = json.at(JsonKeys::ObjectiveWeights);
-            std::transform(objectiveWeightsPairs.begin(),
-                           objectiveWeightsPairs.end(),
-                           perBaseScenarioObjectiveWeights.begin(),
-                           [](const std::pair<double, double> &objectiveWeightsPair) {
-                               return ObjectiveWeights{objectiveWeightsPair.first,
-                                                       objectiveWeightsPair.second};
-                           });
-        }
-        if (json.contains(GET_VARIABLE_NAME(translationalNormalizationEpsilon))) {
-            translationEpsilon
-                = json[GET_VARIABLE_NAME(translationalNormalizationEpsilon)];
-        }
-
-        if (json.contains(GET_VARIABLE_NAME(angularDistanceEpsilon))) {
-            angularDistanceEpsilon = vcg::math::ToRad(
-                static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
-        }
-
-        //        perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
-        //        objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
-        //        objectiveWeights.rotational = 2 - objectiveWeights.translational;
           return true;
       }
 
       std::string toString() const
       {
           std::string settingsString;
-        if (!variablesRanges.empty()) {
+          if (!xRanges.empty()) {
               std::string xRangesString;
-            for (const xRange &range : variablesRanges) {
+              for (const xRange &range : xRanges) {
                   xRangesString += range.toString() + " ";
               }
               settingsString += xRangesString;
           }
           settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
-                          + " Accuracy=" + std::to_string(solverAccuracy);
+                            + " Accuracy=" + std::to_string(solverAccuracy)
-        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
+                            + " TransWeight=" + std::to_string(objectiveWeights.translational)
-             baseScenario++) {
+                            + " RotWeight=" + std::to_string(objectiveWeights.rotational);
-            +" W_t=" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
-                + " W_r="
-                + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational);
-        }
 
           return settingsString;
       }
 
       void writeHeaderTo(csvFile &os) const
       {
-        if (!variablesRanges.empty()) {
+          if (!xRanges.empty()) {
-            for (const xRange &range : variablesRanges) {
+              for (const xRange &range : xRanges) {
                   os << range.label + " max";
                   os << range.label + " min";
               }
           }
           os << "Function Calls";
           os << "Solution Accuracy";
-        os << "Normalization trans epsilon";
+          os << "Normalization strategy";
-        os << "Normalization rot epsilon(deg)";
+          os << "Trans weight";
-        os << JsonKeys::ObjectiveWeights;
+          os << "Rot weight";
-        os << "Optimization parameters";
+          os << "Splitted geo from mat opt";
           //        os << std::endl;
       }
 
       void writeSettingsTo(csvFile &os) const
       {
-        if (!variablesRanges.empty()) {
+          if (!xRanges.empty()) {
-            for (const xRange &range : variablesRanges) {
+              for (const xRange &range : xRanges) {
                   os << range.max;
                   os << range.min;
               }
           }
           os << numberOfFunctionCalls;
           os << solverAccuracy;
-        os << std::to_string(translationEpsilon);
+          os << normalizationStrategyStrings[normalizationStrategy] + "_"
-        os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
+                    + std::to_string(translationNormalizationParameter);
-        std::string objectiveWeightsString;
+          os << objectiveWeights.translational;
-        objectiveWeightsString += "{";
+          os << objectiveWeights.rotational;
-        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
+          os << (splitGeometryMaterialOptimization == true ? "true" : "false");
-             baseScenario++) {
-            objectiveWeightsString
-                += "{" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
-                   + "," + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational)
-                   + "}";
       }
-        objectiveWeightsString += "}";
+  };
-        os << objectiveWeightsString;
-
-        //export optimization parameters
-        std::vector<std::vector<int>> vv;
-        for (const std::vector<OptimizationParameterIndex> &v : optimizationStrategy) {
-            std::vector<int> vi;
-            vi.reserve(v.size());
-            for (const OptimizationParameterIndex &parameter : v) {
-                vi.emplace_back(parameter);
-            }
-            vv.push_back(vi);
-        }
-        os << Utilities::toString(vv);
-    }
-};
 
   inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept
   {
-      const bool haveTheSameObjectiveWeights
-          = std::mismatch(settings1.perBaseScenarioObjectiveWeights.begin(),
-                          settings1.perBaseScenarioObjectiveWeights.end(),
-                          settings2.perBaseScenarioObjectiveWeights.begin())
-                .first
-            != settings1.perBaseScenarioObjectiveWeights.end();
       return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
-             && settings1.variablesRanges == settings2.variablesRanges
+             && settings1.xRanges == settings2.xRanges
-             && settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
+             && settings1.solverAccuracy == settings2.solverAccuracy
-             && settings1.translationEpsilon
+             && settings1.objectiveWeights.translational == settings2.objectiveWeights.translational
-                    == settings2.translationEpsilon;
+             && settings1.translationNormalizationParameter
+                    == settings2.translationNormalizationParameter;
   }
 
   inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@@ -402,9 +222,9 @@ struct Settings
           Element &e = m.elements[ei];
           e.setDimensions(CrossSectionType(beamWidth, beamHeight));
           e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
-          e.dimensions.inertia.J = J;
+          e.inertia.J = J;
-          e.dimensions.inertia.I2 = I2;
+          e.inertia.I2 = I2;
-          e.dimensions.inertia.I3 = I3;
+          e.inertia.I3 = I3;
       }
 
       CoordType center_barycentric(1, 0, 0);
@@ -451,7 +271,6 @@ struct Settings
       std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
       std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
       std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
-      double fullPatternYoungsModulus{0};
 
       PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
       vcg::Triangle3<double> baseTriangle;
@@ -467,7 +286,6 @@ struct Settings
           std::vector<double> perSimulationScenario_translational;
           std::vector<double> perSimulationScenario_rotational;
           std::vector<double> perSimulationScenario_total;
-          std::vector<double> perSimulationScenario_total_unweighted;
       } objectiveValue;
       std::vector<double> perScenario_fullPatternPotentialEnergy;
       std::vector<double> objectiveValueHistory;
@@ -476,8 +294,8 @@ struct Settings
       struct CSVExportingSettings
       {
           bool exportPE{false};
-          bool exportIterationOfMinima{false};
+          bool exportIterationOfMinima{true};
-          bool exportRawObjectiveValue{true};
+          bool exportRawObjectiveValue{false};
           CSVExportingSettings() {}
       };
 
@@ -491,10 +309,9 @@ struct Settings
           inline static std::string FullPatternLabel{"FullPatternLabel"};
           inline static std::string Settings{"OptimizationSettings"};
           inline static std::string FullPatternPotentialEnergies{"PE"};
-          inline static std::string fullPatternYoungsModulus{"youngsModulus"};
       };
 
-      void save(const std::string &saveToPath, const bool shouldExportDebugFiles = false)
+      void save(const string &saveToPath, const bool shouldExportDebugFiles = false)
       {
           //clear directory
           if (std::filesystem::exists(saveToPath)) {
@@ -545,7 +362,6 @@ struct Settings
                   = perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
           }
           json_optimizationResults[JsonKeys::FullPatternPotentialEnergies] = fullPatternPE;
-          json_optimizationResults[JsonKeys::fullPatternYoungsModulus] = fullPatternYoungsModulus;
           ////Save to json file
           std::filesystem::path jsonFilePath(
               std::filesystem::path(saveToPath).append(JsonKeys::filename));
@@ -567,8 +383,7 @@ struct Settings
                       = fullPatternSimulationJobs[simulationScenarioIndex];
                   std::filesystem::path simulationJobFolderPath(
                       std::filesystem::path(simulationJobsPath)
-                          .append(std::to_string(simulationScenarioIndex) + "_"
+                          .append(pFullPatternSimulationJob->label));
-                                  + pFullPatternSimulationJob->label));
                   std::filesystem::create_directories(simulationJobFolderPath);
                   const auto fullPatternDirectoryPath
                       = std::filesystem::path(simulationJobFolderPath).append("Full");
@@ -598,13 +413,12 @@ struct Settings
 #endif
       }
 
-      bool load(const std::string &loadFromPath, const bool &shouldLoadDebugFiles = false)
+      bool load(const string &loadFromPath, const bool &shouldLoadDebugFiles = false)
       {
           assert(std::filesystem::is_directory(loadFromPath));
           std::filesystem::path jsonFilepath(
               std::filesystem::path(loadFromPath).append(JsonKeys::filename));
           if (!std::filesystem::exists(jsonFilepath)) {
-              std::cerr << "Input path does not exist:" << loadFromPath << std::endl;
               return false;
           }
           //Load optimal X
@@ -651,83 +465,61 @@ struct Settings
           baseTriangle.P2(0) = CoordType(baseTriangleVertices[6],
                                          baseTriangleVertices[7],
                                          baseTriangleVertices[8]);
-          if (json_optimizationResults.contains(JsonKeys::fullPatternYoungsModulus)) {
-              fullPatternYoungsModulus = json_optimizationResults.at(
-                  JsonKeys::fullPatternYoungsModulus);
-          } else {
-              fullPatternYoungsModulus = 1 * 1e9;
-          }
           if (shouldLoadDebugFiles) {
               const std::filesystem::path simulationJobsFolderPath(
                   std::filesystem::path(loadFromPath).append("SimulationJobs"));
-
+              for (const auto &directoryEntry :
-              std::vector<std::filesystem::path> sortedByName;
+                   filesystem::directory_iterator(simulationJobsFolderPath)) {
-              for (auto &entry : std::filesystem::directory_iterator(simulationJobsFolderPath))
+                  const auto simulationScenarioPath = directoryEntry.path();
-                  sortedByName.push_back(entry.path());
-
-              std::sort(sortedByName.begin(), sortedByName.end(), &Utilities::compareNat);
-
-              for (const auto &simulationScenarioPath : sortedByName) {
                   if (!std::filesystem::is_directory(simulationScenarioPath)) {
                       continue;
                   }
                   // Load full pattern files
-                  for (const auto &fileEntry : std::filesystem::directory_iterator(
+                  for (const auto &fileEntry : filesystem::directory_iterator(
                            std::filesystem::path(simulationScenarioPath).append("Full"))) {
                       const auto filepath = fileEntry.path();
                       if (filepath.extension() == ".json") {
                           SimulationJob job;
+
                           job.load(filepath.string());
-                          job.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
                           fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(job));
                       }
                   }
 
-                  const auto fullJobFilepath = Utilities::getFilepathWithExtension(
+                  // Load reduced pattern files and apply optimal parameters
-                      std::filesystem::path(simulationScenarioPath).append("Full"), ".json");
+                  for (const auto &fileEntry : filesystem::directory_iterator(
-                  SimulationJob fullJob;
+                           std::filesystem::path(simulationScenarioPath).append("Reduced"))) {
-                  fullJob.load(fullJobFilepath.string());
+                      const auto filepath = fileEntry.path();
-                  fullJob.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
+                      if (filepath.extension() == ".json") {
-                  fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(fullJob));
+                          SimulationJob job;
-
+                          job.load(filepath.string());
-                  const auto reducedJobFilepath = Utilities::getFilepathWithExtension(
-                      std::filesystem::path(simulationScenarioPath).append("Reduced"), ".json");
-                  SimulationJob reducedJob;
-                  reducedJob.load(reducedJobFilepath.string());
                           //                          applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
-                  applyOptimizationResults_elements(*this, reducedJob.pMesh);
+                          applyOptimizationResults_elements(*this, job.pMesh);
                           reducedPatternSimulationJobs.push_back(
-                      std::make_shared<SimulationJob>(reducedJob));
+                              std::make_shared<SimulationJob>(job));
                       }
                   }
-          settings.load(std::filesystem::path(loadFromPath).append(Settings::defaultFilename));
+              }
+          }
+          settings.load(loadFromPath);
 
           return true;
       }
 
       template<typename MeshType>
       static void applyOptimizationResults_innerHexagon(
-          const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
+          const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
           const vcg::Triangle3<double> &patternBaseTriangle,
           MeshType &reducedPattern)
       {
           std::unordered_map<std::string, double>
               optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
                                 reducedPattern_optimizationResults.optimalXNameValuePairs.end());
-          assert(
+          assert(optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle"));
-              (optimalXVariables.contains("R") && optimalXVariables.contains("Theta"))
-              || (optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")));
-          if (optimalXVariables.contains("HexSize")) {
           applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
                                                 optimalXVariables.at("HexAngle"),
                                                 patternBaseTriangle,
                                                 reducedPattern);
-              return;
-          }
-          applyOptimizationResults_innerHexagon(optimalXVariables.at("R"),
-                                                optimalXVariables.at("Theta"),
-                                                patternBaseTriangle,
-                                                reducedPattern);
       }
 
       template<typename MeshType>
@@ -768,8 +560,8 @@ struct Settings
       }
 
       static void applyOptimizationResults_elements(
-          const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
+          const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
-          const std::shared_ptr<SimulationMesh> &pReducedPattern_simulationMesh)
+          const shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
       {
           assert(CrossSectionType::name == RectangularBeamDimensions::name);
           // Set reduced parameters fron the optimization results
@@ -783,8 +575,8 @@ struct Settings
               const double beamWidth = std::sqrt(A);
               const double beamHeight = beamWidth;
               CrossSectionType elementDimensions(beamWidth, beamHeight);
-              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
                   e.setDimensions(elementDimensions);
               }
           }
@@ -794,8 +586,8 @@ struct Settings
           if (optimalXVariables.contains(ymLabel)) {
               const double E = optimalXVariables.at(ymLabel);
               const ElementMaterial elementMaterial(poissonsRatio, E);
-              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
                   e.setMaterial(elementMaterial);
               }
           }
@@ -803,30 +595,30 @@ struct Settings
           const std::string JLabel = "J";
           if (optimalXVariables.contains(JLabel)) {
               const double J = optimalXVariables.at(JLabel);
-              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.dimensions.inertia.J = J;
+                  e.inertia.J = J;
               }
           }
 
           const std::string I2Label = "I2";
           if (optimalXVariables.contains(I2Label)) {
               const double I2 = optimalXVariables.at(I2Label);
-              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.dimensions.inertia.I2 = I2;
+                  e.inertia.I2 = I2;
               }
           }
 
           const std::string I3Label = "I3";
           if (optimalXVariables.contains(I3Label)) {
               const double I3 = optimalXVariables.at(I3Label);
-              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.dimensions.inertia.I3 = I3;
+                  e.inertia.I3 = I3;
               }
           }
-          pReducedPattern_simulationMesh->reset();
+          pTiledReducedPattern_simulationMesh->reset();
       }
 
 #if POLYSCOPE_DEFINED
@@ -923,7 +715,7 @@ struct Settings
       }
     }
 
-    void writeHeaderTo(csvFile &os) const
+    void writeHeaderTo(csvFile &os)
     {
         if (exportSettings.exportRawObjectiveValue) {
             os << "Total raw Obj value";
@@ -959,20 +751,13 @@ struct Settings
         //        os << "notes";
     }
 
-    void writeHeaderTo(std::vector<csvFile *> &vectorOfPointersToOutputStreams) const
-    {
-        for (csvFile *outputStream : vectorOfPointersToOutputStreams) {
-            writeHeaderTo(*outputStream);
-        }
-    }
-
     void writeResultsTo(csvFile &os) const
     {
         if (exportSettings.exportRawObjectiveValue) {
             os << objectiveValue.totalRaw;
         }
         os << objectiveValue.total;
-        if (exportSettings.exportIterationOfMinima && !objectiveValueHistory_iteration.empty()) {
+        if (exportSettings.exportIterationOfMinima) {
             os << objectiveValueHistory_iteration.back();
         }
         for (int simulationScenarioIndex = 0;
@@ -1002,13 +787,6 @@ struct Settings
         //        os << notes;
     }
 
-    void writeResultsTo(std::vector<csvFile *> &vectorOfPointersToOutputStreams) const
-    {
-        for (csvFile *&outputStream : vectorOfPointersToOutputStreams) {
-            writeResultsTo(*outputStream);
-        }
-    }
-
     void writeMinimaInfoTo(csvFile &outputCsv)
     {
         outputCsv << "Iteration";

simulation_structs.hpp

@@ -1,6 +1,5 @@
 #ifndef SIMULATIONSTRUCTS_HPP
 #define SIMULATIONSTRUCTS_HPP
-#include <fstream>
 
 namespace Eigen {
 template <class Matrix>
@@ -24,13 +23,13 @@ void read_binary(const std::string &filename, Matrix &matrix) {
   in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
   in.close();
 }
-//const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
+const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
-//template<class Matrix>
+template<class Matrix>
-//void writeToCSV(const std::string &filename, Matrix &matrix)
+void writeToCSV(const std::string &filename, Matrix &matrix)
-//{
+{
-//    ofstream file(filename.c_str());
+    ofstream file(filename.c_str());
-//    file << matrix.format(CSVFormat);
+    file << matrix.format(CSVFormat);
-//}
+}
 } // namespace Eigen
 
 #include "simulationmesh.hpp"
@@ -113,20 +112,18 @@ class SimulationJob {
   } jsonLabels;
 
 public:
-    inline static std::string jsonDefaultFileName = "SimulationJob.json";
   std::shared_ptr<SimulationMesh> pMesh;
   std::string label{"empty_job"};
   std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
   std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
   std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
-    SimulationJob(const std::shared_ptr<SimulationMesh> &m,
+  SimulationJob(
-                  const std::string &label,
+      const std::shared_ptr<SimulationMesh> &m, const std::string &label,
       const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
       const std::unordered_map<VertexIndex, Vector6d> &ef = {},
       const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
-        : pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef),
+      : pMesh(m), label(label), constrainedVertices(cv),
-          nodalForcedDisplacements(fd)
+        nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
-    {}
   SimulationJob() {}
   SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
 
@@ -135,32 +132,32 @@ public:
       return constrainedVertices.empty() && nodalExternalForces.empty()
              && nodalForcedDisplacements.empty() && pMesh == nullptr;
   }
-  void remap(const std::unordered_map<size_t, size_t> &sourceToDestinationViMap,
+  void remap(const std::unordered_map<size_t, size_t> &thisToOtherViMap,
-             SimulationJob &destination_simulationJob) const
+             SimulationJob &simulationJobMapped) const
   {
-      std::unordered_map<VertexIndex, std::unordered_set<int>> destination_fixedVertices;
+      assert(simulationJobMapped.pMesh->VN() != 0);
-      for (const auto &source_fixedVertex : this->constrainedVertices) {
+      std::unordered_map<VertexIndex, std::unordered_set<int>> reducedModelFixedVertices;
-          destination_fixedVertices[sourceToDestinationViMap.at(source_fixedVertex.first)]
+      for (auto fullModelFixedVertex : this->constrainedVertices) {
-              = source_fixedVertex.second;
+          reducedModelFixedVertices[thisToOtherViMap.at(fullModelFixedVertex.first)]
+              = fullModelFixedVertex.second;
       }
 
-      std::unordered_map<VertexIndex, Vector6d> destination_nodalForces;
+      std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
-      for (const auto &source_nodalForces : this->nodalExternalForces) {
+      for (auto fullModelNodalForce : this->nodalExternalForces) {
-          destination_nodalForces[sourceToDestinationViMap.at(source_nodalForces.first)]
+          reducedModelNodalForces[thisToOtherViMap.at(fullModelNodalForce.first)]
-              = source_nodalForces.second;
+              = fullModelNodalForce.second;
       }
 
-      std::unordered_map<VertexIndex, Eigen::Vector3d> destination_forcedDisplacements;
+      std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
-      for (const auto &source_forcedDisplacements : this->nodalForcedDisplacements) {
+      for (auto fullForcedDisplacement : this->nodalForcedDisplacements) {
-          destination_forcedDisplacements[sourceToDestinationViMap.at(
+          reducedNodalForcedDisplacements[thisToOtherViMap.at(fullForcedDisplacement.first)]
-              source_forcedDisplacements.first)]
+              = fullForcedDisplacement.second;
-              = source_forcedDisplacements.second;
       }
 
-      destination_simulationJob.constrainedVertices = destination_fixedVertices;
+      simulationJobMapped.constrainedVertices = reducedModelFixedVertices;
-      destination_simulationJob.nodalExternalForces = destination_nodalForces;
+      simulationJobMapped.nodalExternalForces = reducedModelNodalForces;
-      destination_simulationJob.label = this->getLabel();
+      simulationJobMapped.label = this->getLabel();
-      destination_simulationJob.nodalForcedDisplacements = destination_forcedDisplacements;
+      simulationJobMapped.nodalForcedDisplacements = reducedNodalForcedDisplacements;
   }
   SimulationJob getCopy() const {
     SimulationJob jobCopy;
@@ -379,7 +376,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
           return;
       }
       std::vector<std::array<double, 3>> nodeColors(pMesh->VN());
-      for (const auto &fixedVertex : constrainedVertices) {
+      for (auto fixedVertex : constrainedVertices) {
           const bool hasRotationalDoFConstrained = fixedVertex.second.contains(3)
                                                    || fixedVertex.second.contains(4)
                                                    || fixedVertex.second.contains(5);
@@ -395,8 +392,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
           }
       }
       if (!nodalForcedDisplacements.empty()) {
-          for (const std::pair<VertexIndex, Eigen::Vector3d> &viDisplPair :
+          for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : nodalForcedDisplacements) {
-               nodalForcedDisplacements) {
               const VertexIndex vi = viDisplPair.first;
               nodeColors[vi][0] += 1;
               nodeColors[vi][0] /= 2;
@@ -491,7 +487,6 @@ struct SimulationResults
     /*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
  * vector contains the same info using euler angles
 */
-    inline const static std::string defaultJsonFilename{"SimulationResults.json"};
     bool converged{false};
     std::shared_ptr<SimulationJob> pJob;
     SimulationHistory history;
@@ -543,7 +538,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);
@@ -556,16 +551,6 @@ struct SimulationResults
         std::filesystem::create_directories(resultsFolderPath);
         Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m);
 
-        nlohmann::json json;
-
-        json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
-
-        std::filesystem::path jsonFilePath(
-            std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
-        std::ofstream jsonFile(jsonFilePath.string());
-        jsonFile << json;
-        jsonFile.close();
-
         saveDeformedModel(resultsFolderPath.string());
     }
 
@@ -583,41 +568,34 @@ 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);
+        //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") {
+                Eigen::MatrixXd displacements_eigen;
+                Eigen::read_binary(entry.path().string(), displacements_eigen);
+                displacements = Utilities::fromEigenMatrix(displacements_eigen);
+                break;
             }
-    void load(const std::filesystem::path &loadFromPath, const std::shared_ptr<SimulationJob> &pJob)
-    {
-        this->pJob = pJob;
-        load(loadFromPath);
-    }
-    static double computeDistance(
-        const SimulationResults &resultsA,
-        const SimulationResults &resultsB,
-        const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap)
-    {
-        double distance = 0;
-        for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
-            const double vertexToVertexDistance
-                = (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
-                   - resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
-                      .norm();
-            distance += vertexToVertexDistance;
-        }
-        return distance;
         }
 
-    double computeDistance(const SimulationResults &other,
+        rotationalDisplacementQuaternion.resize(pJob->pMesh->VN());
-                           const std::unordered_map<VertexIndex, VertexIndex> &thisToOtherViMap) const
+        for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
-    {
+            rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3],
-        return computeDistance(*this, other, thisToOtherViMap);
+                                                                     Eigen::Vector3d::UnitX())
+                                                   * Eigen::AngleAxisd(displacements[vi][4],
+                                                                       Eigen::Vector3d::UnitY())
+                                                   * Eigen::AngleAxisd(displacements[vi][5],
+                                                                       Eigen::Vector3d::UnitZ());
+        }
     }
 
 #ifdef POLYSCOPE_DEFINED
-    void unregister() const
+  void unregister() const {
-    {
     if (!polyscope::hasCurveNetwork(getLabel())) {
-            std::cerr << "No curve network registered with a name: " << getLabel() << std::endl;
+      std::cerr << "No curve network registered with a name: " << getLabel()
+                << std::endl;
       std::cerr << "Nothing to remove." << std::endl;
       return;
     }
@@ -635,13 +613,9 @@ struct SimulationResults
       const std::shared_ptr<SimulationMesh> &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,
+                                                                             mesh->getEigenVertices(),
-                                                                               edges);
+                                                                             mesh->getEigenEdges());
-        } else {
-            polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel());
       }
       polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
       polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
@@ -661,18 +635,13 @@ struct SimulationResults
       Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);
 
       //      std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
-        //                  std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
+      std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
-        //      std::unordered_set<int> 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<double> Rx(Eigen::AngleAxis(nodalDisplacement[2],Eigen::Vector3d(1, 0, 0)));
+          if (interfaceNodes.contains(vi)) {
-            //           Eigen::Quaternion<double> Ry(Eigen::AngleAxis(nodalDisplacement[4],Eigen::Vector3d(0, 1, 0)));
-            //           Eigen::Quaternion<double> Rz(Eigen::AngleAxis(nodalDisplacement[5],Eigen::Vector3d(0, 0, 1)));
-            //           Eigen::Quaternion<double> 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);
@@ -688,14 +657,14 @@ struct SimulationResults
               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 {
+          } else {
-            //              framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
+              framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
-            //              framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
+              framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
-            //              framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
+              framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
-            //              framesX_initial.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);
+              framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
-            //              framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
+              framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
-            //          }
+          }
       }
       polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
                                                             + nodalDisplacements);
@@ -723,10 +692,10 @@ struct SimulationResults
       polyscopeHandle_frameZ->setVectorColor(
           /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
 
-        auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
+      //      auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
-        if (!polyscopeHandle_initialMesh) {
+      //      if (!polyscopeHandle_initialMesh) {
-            polyscopeHandle_initialMesh = mesh->registerForDrawing();
+      //          polyscopeHandle_initialMesh = mesh->registerForDrawing();
-        }
+      //      }
 
       //      auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
       //                                                ->addNodeVectorQuantity("FrameX", framesX_initial);
@@ -748,66 +717,24 @@ struct SimulationResults
       //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
       //      //}
       pJob->registerForDrawing(getLabel());
-        //        static bool wasExecuted =false;
+      static bool wasExecuted = false;
-        //        if (!wasExecuted) {
+      if (!wasExecuted) {
-        //            std::function<void()> callback = [&]() {
+          std::function<void()> callback = [&]() {
-        //                static bool showFrames = shouldShowFrames;
+              static bool showFrames = shouldShowFrames;
 
-        //                if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
+              if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
-        //                    polyscopeHandle_frameX->setEnabled(showFrames);
+                  polyscopeHandle_frameX->setEnabled(showFrames);
-        //                    polyscopeHandle_frameY->setEnabled(showFrames);
+                  polyscopeHandle_frameY->setEnabled(showFrames);
-        //                    polyscopeHandle_frameZ->setEnabled(showFrames);
+                  polyscopeHandle_frameZ->setEnabled(showFrames);
-        //                }
+              }
-        //            };
+          };
 
-        //            PolyscopeInterface::addUserCallback(callback);
+          PolyscopeInterface::addUserCallback(callback);
-        //            wasExecuted = true;
+          wasExecuted = true;
-        //        }
+      }
       return polyscopeHandle_deformedEdmeMesh;
   }
 #endif
-private:
-    bool 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 (std::filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
-                Eigen::MatrixXd displacements_eigen;
-                Eigen::read_binary(entry.path().string(), displacements_eigen);
-                displacements = Utilities::fromEigenMatrix(displacements_eigen);
-                break;
-            }
-        }
-
-        rotationalDisplacementQuaternion.resize(pJob->pMesh->VN());
-        for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
-            rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3],
-                                                                     Eigen::Vector3d::UnitX())
-                                                   * Eigen::AngleAxisd(displacements[vi][4],
-                                                                       Eigen::Vector3d::UnitY())
-                                                   * Eigen::AngleAxisd(displacements[vi][5],
-                                                                       Eigen::Vector3d::UnitZ());
-        }
-
-        const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
-                                                       .append(defaultJsonFilename);
-        if (!std::filesystem::exists(jsonFilepath)) {
-            std::cerr << "Simulation results could not be loaded because filepath does "
-                         "not exist:"
-                      << jsonFilepath << std::endl;
-            return false;
-        }
-        std::ifstream ifs(jsonFilepath);
-        nlohmann::json json;
-        ifs >> json;
-        if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
-            internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
-        }
-        return true;
-    }
 };
 
 #endif // SIMULATIONHISTORY_HPP

simulationhistoryplotter.hpp

@@ -1,8 +1,8 @@
 #ifndef SIMULATIONHISTORYPLOTTER_HPP
 #define SIMULATIONHISTORYPLOTTER_HPP
 
-#include "simulation_structs.hpp"
 #include "simulationmesh.hpp"
+#include "simulation_structs.hpp"
 #include "utilities.hpp"
 #include <algorithm>
 #include <matplot/matplot.h>
@@ -15,8 +15,10 @@ struct SimulationResultsReporter {
 
   void writeStatistics(const SimulationResults &results,
                        const std::string &reportFolderPath) {
-      std::ofstream file;
+
-      file.open(std::filesystem::path(reportFolderPath).append("results.txt").string());
+    ofstream file;
+    file.open(
+        std::filesystem::path(reportFolderPath).append("results.txt").string());
 
     const size_t numberOfSteps = results.history.numberOfSteps;
     file << "Number of steps " << numberOfSteps << "\n";
@@ -112,23 +114,35 @@ struct SimulationResultsReporter {
       if (YvaluesToPlot.size() < 2) {
           return;
       }
-      std::vector<double> colors(YvaluesToPlot.size(), 0.5);
+      auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
-      std::vector<double> markerSizes(YvaluesToPlot.size(), 5);
+      std::vector<double> colors(x.size(), 0.5);
+      std::vector<double> markerSizes(x.size(), 5);
       if (!markPoints.empty()) {
           for (const auto pointIndex : markPoints) {
               colors[pointIndex] = 0.9;
               markerSizes[pointIndex] = 14;
           }
       }
-      std::vector<double> 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 std::string &graphSuffix) {
-  {
+    const auto graphsFolder =
-      const auto graphsFolder = std::filesystem::path(reportFolderPath).append("Graphs");
+        std::filesystem::path(reportFolderPath).append("Graphs");
     std::filesystem::remove_all(graphsFolder);
     std::filesystem::create_directory(graphsFolder.string());
 

simulationmesh.cpp

@@ -15,8 +15,8 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
 
   elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
       *this, std::string("Elements"));
-  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(*this,
+  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
-                                                                           std::string("Nodes"));
+      *this, std::string("Nodes"));
   initializeNodes();
   initializeElements();
 }
@@ -198,6 +198,7 @@ void SimulationMesh::reset() {
 }
 
 void SimulationMesh::initializeElements() {
+  computeElementalProperties();
   for (const EdgeType &e : edge) {
     Element &element = elements[e];
     element.ei = getIndex(e);
@@ -242,6 +243,8 @@ void SimulationMesh::updateElementalLengths() {
     const vcg::Segment3<double> s(p0, p1);
     const double elementLength = s.Length();
     elements[e].length = elementLength;
+    int i = 0;
+    i++;
   }
 }
 
@@ -257,6 +260,7 @@ 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();
   }
 }
@@ -284,8 +288,7 @@ std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
   return beamMaterial;
 }
 
-bool SimulationMesh::load(const std::string &plyFilename)
+bool SimulationMesh::load(const string &plyFilename) {
-{
   this->Clear();
   //  assert(plyFileHasAllRequiredFields(plyFilename));
   // Load the ply file
@@ -325,7 +328,7 @@ bool SimulationMesh::load(const std::string &plyFilename)
   //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
   //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
   //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
-    if (!VCGEdgeMesh::load(plyFilename)) {
+  if (!load(plyFilename)) {
       return false;
   }
 
@@ -336,7 +339,6 @@ bool SimulationMesh::load(const std::string &plyFilename)
   vcg::tri::UpdateTopology<SimulationMesh>::VertexEdge(*this);
   initializeNodes();
   initializeElements();
-    setBeamMaterial(0.3, 1 * 1e9);
   updateEigenEdgeAndVertices();
 
   //  if (!handleBeamProperties._handle->data.empty()) {
@@ -356,7 +358,8 @@ bool SimulationMesh::load(const std::string &plyFilename)
   if (normalsAreAbsent) {
     CoordType normalVector(0, 0, 1);
     std::cout << "Warning: Normals are missing from " << plyFilename
-                  << ". Added normal vector:" << toString(normalVector) << std::endl;
+              << ". Added normal vector:" << toString(normalVector)
+              << std::endl;
     for (auto &v : vert) {
       v.N() = normalVector;
     }
@@ -395,7 +398,7 @@ bool SimulationMesh::save(const std::string &plyFilename)
     //        != 1) {
     //        return false;
     //    }
-    if (!VCGEdgeMesh::save(plyFilename)) {
+    if (!save(plyFilename)) {
         return false;
     }
     return true;
@@ -453,59 +456,60 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
   return angle;
 }
 
-//void Element::computeMaterialProperties(const ElementMaterial &material) {
+void Element::computeMaterialProperties(const ElementMaterial &material) {
-//    G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
+    G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
-//}
+}
 
-//void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
+void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
-//{
+{
-//    A = dimensions.b * dimensions.h;
+    A = dimensions.b * dimensions.h;
-//}
+}
 
-//void Element::computeDimensionsProperties(
+void Element::computeDimensionsProperties(
-//    const RectangularBeamDimensions &dimensions) {
+    const RectangularBeamDimensions &dimensions) {
-//  assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
+  assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
-//  computeCrossSectionArea(dimensions, A);
+  computeCrossSectionArea(dimensions, A);
-//  computeMomentsOfInertia(dimensions, dimensions.inertia);
+  computeMomentsOfInertia(dimensions, inertia);
-//}
+}
 
-//void Element::computeDimensionsProperties(
+void Element::computeDimensionsProperties(
-//    const CylindricalBeamDimensions &dimensions) {
+    const CylindricalBeamDimensions &dimensions) {
-//  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
+  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
-//  A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
+  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;
+  inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
-//  dimensions.inertia.I3 = dimensions.inertia.I2;
+  inertia.I3 = inertia.I2;
-//  dimensions.inertia.J = dimensions.inertia.I2 + dimensions.inertia.I3;
+  inertia.J = inertia.I2 + inertia.I3;
-//}
+}
 
 void Element::setDimensions(const CrossSectionType &dimensions) {
   this->dimensions = dimensions;
-  assert(this->dimensions.A == dimensions.A);
+  computeDimensionsProperties(dimensions);
-  //  computeDimensionsProperties(dimensions);
   updateRigidity();
 }
 
-void Element::setMaterial(const ElementMaterial &material)
+void Element::setMaterial(const ElementMaterial &material) {
-{
   this->material = material;
-    //  computeMaterialProperties(material);
+  computeMaterialProperties(material);
   updateRigidity();
 }
 
 double Element::getMass(const double &materialDensity)
 {
-    const double beamVolume = dimensions.A * length;
+    const double beamVolume = A * length;
     return beamVolume * materialDensity;
 }
 
-void Element::updateRigidity() {
+void Element::computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
-    //    assert(initialLength != 0);
+                                      Element::MomentsOfInertia &inertia)
-    rigidity.axial = material.youngsModulus * dimensions.A / initialLength;
+{
-    //    assert(rigidity.axial != 0);
+    inertia.I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
-    rigidity.torsional = material.G * dimensions.inertia.J / initialLength;
+    inertia.I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
-    //    assert(rigidity.torsional != 0);
+    inertia.J = inertia.I2 + inertia.I3;
-    rigidity.firstBending = 2 * material.youngsModulus * dimensions.inertia.I2 / initialLength;
+}
-    //    assert(rigidity.firstBending != 0);
+
-    rigidity.secondBending = 2 * material.youngsModulus * dimensions.inertia.I3 / initialLength;
+void Element::updateRigidity() {
-    //    assert(rigidity.secondBending != 0);
+  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;
 }

simulationmesh.hpp

@@ -1,9 +1,9 @@
 #ifndef SIMULATIONMESH_HPP
 #define SIMULATIONMESH_HPP
 
+#include "Eigen/Dense"
 #include "edgemesh.hpp"
 #include "trianglepatterngeometry.hpp"
-#include <Eigen/Dense>
 
 //extern bool drawGlobal;
 struct Element;
@@ -34,7 +34,7 @@ public:
 
   std::vector<VCGEdgeMesh::EdgePointer>
   getIncidentElements(const VCGEdgeMesh::VertexType &v);
-  virtual bool load(const std::string &plyFilename);
+  virtual bool load(const string &plyFilename);
   std::vector<CrossSectionType> getBeamDimensions();
   std::vector<ElementMaterial> getBeamMaterial();
 
@@ -62,14 +62,26 @@ 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 RectangularBeamDimensions &dimensions);
-  //  void computeDimensionsProperties(const CylindricalBeamDimensions &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;

trianglepatterngeometry.cpp

@@ -86,8 +86,8 @@ std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
   return verts;
 }
 
-PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
+PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
-{
+
   const size_t fanSize = PatternGeometry().getFanSize();
   PatternGeometry fan(createFan(pattern));
   PatternGeometry tile(fan);
@@ -99,15 +99,11 @@ PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
     vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
   }
   vcg::Matrix44d T;
-    //    const double centerAngle = 2 * M_PI / fanSize;
+  const double centerAngle = 2 * M_PI / fanSize;
-    //    const double triangleHeight = std::sin((M_PI - centerAngle) / 2)
+  const double triangleHeight =
-    //                                  * PatternGeometry().triangleEdgeSize;
+      std::sin((M_PI - centerAngle) / 2) * PatternGeometry().triangleEdgeSize;
-    vcg::tri::UpdateBounding<PatternGeometry>::Box(fan);
-    const double triangleHeight = fan.bbox.DimY() / 2;
   T.SetTranslate(0, -2 * triangleHeight, 0);
   vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
-    //    fan.registerForDrawing();
-    //    polyscope::show();
 
   PatternGeometry fanOfFan = createFan(fan);
   vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
@@ -149,9 +145,9 @@ PatternGeometry::PatternGeometry(PatternGeometry &other) {
   vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
 }
 
-bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
+bool PatternGeometry::load(const string &plyFilename)
 {
-    if (!VCGEdgeMesh::load(meshFilePath)) {
+    if (!VCGEdgeMesh::load(plyFilename)) {
         return false;
     }
     vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
@@ -781,22 +777,6 @@ double PatternGeometry::computeBaseTriangleHeight() const
     return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
 }
 
-void PatternGeometry::deleteDanglingVertices()
-{
-    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu;
-    VCGEdgeMesh::deleteDanglingVertices(pu);
-    if (!pu.remap.empty()) {
-        interfaceNodeIndex
-            = pu.remap[interfaceNodeIndex]; //TODO:Could this be automatically be determined?
-    }
-}
-
-void PatternGeometry::deleteDanglingVertices(
-    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu)
-{
-    VCGEdgeMesh::deleteDanglingVertices(pu);
-}
-
 vcg::Triangle3<double> PatternGeometry::getBaseTriangle() const
 {
     return baseTriangle;
@@ -834,7 +814,6 @@ PatternGeometry::PatternGeometry(
   updateEigenEdgeAndVertices();
 }
 
-//TODO: refactor such that it takes as an input a single pattern and tiles into the desired faces without requiring the each face will contain a single pattern
 std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
     std::vector<ConstPatternGeometry> &patterns,
     const std::vector<int> &connectToNeighborsVi,

trianglepatterngeometry.hpp

@@ -46,7 +46,7 @@ private:
    * vcglib interface.
    * */
   PatternGeometry(PatternGeometry &other);
-  bool load(const std::filesystem::path &meshFilePath) override;
+  bool load(const std::string &plyFilename) override;
   void add(const std::vector<vcg::Point3d> &vertices);
   void add(const std::vector<vcg::Point2i> &edges);
   void add(const std::vector<vcg::Point3d> &vertices, const std::vector<vcg::Point2i> &edges);
@@ -124,9 +124,6 @@ private:
       std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
   std::unordered_set<VertexIndex> getInterfaceNodes(const std::vector<size_t> &numberOfNodesPerSlot);
   bool isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes);
-  void deleteDanglingVertices() override;
-  void deleteDanglingVertices(
-      vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu) override;
 };
 
 #endif // FLATPATTERNGEOMETRY_HPP

utilities.hpp

@@ -2,14 +2,12 @@
 #define UTILITIES_H
 
 #include <Eigen/Dense>
-#include <algorithm>
-#include <array>
-#include <chrono>
 #include <filesystem>
 #include <fstream>
-#include <iterator>
-#include <numeric>
 #include <regex>
+#include <iterator>
+#include <algorithm>
+#include <array>
 
 #define GET_VARIABLE_NAME(Variable) (#Variable)
 
@@ -142,110 +140,7 @@ struct Vector6d : public std::array<double, 6> {
 };
 
 namespace Utilities {
-inline bool compareNat(const std::string &a, const std::string &b)
+inline void parseIntegers(const std::string &str, std::vector<size_t> &result) {
-{
-    if (a.empty())
-        return true;
-    if (b.empty())
-        return false;
-    if (std::isdigit(a[0]) && !std::isdigit(b[0]))
-        return true;
-    if (!std::isdigit(a[0]) && std::isdigit(b[0]))
-        return false;
-    if (!std::isdigit(a[0]) && !std::isdigit(b[0])) {
-        if (std::toupper(a[0]) == std::toupper(b[0]))
-            return compareNat(a.substr(1), b.substr(1));
-        return (std::toupper(a[0]) < std::toupper(b[0]));
-    }
-
-    // Both strings begin with digit --> parse both numbers
-    std::istringstream issa(a);
-    std::istringstream issb(b);
-    int ia, ib;
-    issa >> ia;
-    issb >> ib;
-    if (ia != ib)
-        return ia < ib;
-
-    // Numbers are the same --> remove numbers and recurse
-    std::string anew, bnew;
-    std::getline(issa, anew);
-    std::getline(issb, bnew);
-    return (compareNat(anew, bnew));
-}
-
-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;
-}
-
-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;
-}
-
-inline std::string_view trimLeftAndRightSpaces(std::string_view str)
-{
-std::string_view trimmedString=str;
-    trimmedString = leftTrimSpaces(trimmedString);
-    trimmedString = rightTrimSpaces(trimmedString);
-    return trimmedString;
-}
-
-template<typename InputIt>
-inline void normalize(InputIt itBegin, InputIt itEnd)
-{
-    const auto squaredSumOfElements = std::accumulate(itBegin,
-                                                      itEnd,
-                                                      0.0,
-                                                      [](const auto &sum, const auto &el) {
-                                                          return sum + el * el;
-                                                      });
-    assert(squaredSumOfElements != 0);
-    std::transform(itBegin, itEnd, itBegin, [&](auto &element) {
-        return element / std::sqrt(squaredSumOfElements);
-    });
-}
-
-inline std::vector<std::string> split(const std::string& text, std::string delim)
-{
-    std::vector<std::string> 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<std::vector<int>> &vv)
-{
-    std::string s;
-    s.append("{");
-    for (const std::vector<int> &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<size_t> &result)
-{
   typedef std::regex_iterator<std::string::const_iterator> re_iterator;
   typedef re_iterator::value_type re_iterated;
 
@@ -254,9 +149,8 @@ inline void parseIntegers(const std::string &str, std::vector<size_t> &result)
   re_iterator rit(str.begin(), str.end(), re);
   re_iterator rend;
 
-    std::transform(rit, rend, std::back_inserter(result), [](const re_iterated &it) {
+  std::transform(rit, rend, std::back_inserter(result),
-        return std::stoi(it[1]);
+                 [](const re_iterated &it) { return std::stoi(it[1]); });
-    });
 }
 
 inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
@@ -285,6 +179,7 @@ inline std::vector<Vector6d> 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(),
@@ -311,19 +206,6 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
 //  return true;
 //}
 
-inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
-                                                      const std::string &extension)
-{
-    for (const std::filesystem::directory_entry &dirEntry :
-         std::filesystem::directory_iterator(folderPath)) {
-        if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
-            return std::filesystem::path(dirEntry);
-        }
-    }
-
-    return "";
-}
-
 } // namespace Utilities
 
 #ifdef POLYSCOPE_DEFINED
@@ -344,7 +226,7 @@ inline void mainCallback()
                                // instead of full width. Must have
                                // matching PopItemWidth() below.
 
-    for (std::function<void()> &userCallback : globalPolyscopeData.userCallbacks) {
+    for (std::function<void()> userCallback : globalPolyscopeData.userCallbacks) {
         userCallback();
     }
 
@@ -375,8 +257,7 @@ inline void init()
     polyscope::view::upDir = polyscope::view::UpDir::ZUp;
 
     polyscope::state::userCallback = &mainCallback;
-    polyscope::options::autocenterStructures = false;
+    //    polyscope::options::autocenterStructures = true;
-    polyscope::options::autoscaleStructures = false;
 }
 using PolyscopeLabel = std::string;
 inline std::pair<PolyscopeLabel, size_t> getSelection()
@@ -394,12 +275,9 @@ inline void registerWorldAxes()
 
     Eigen::MatrixX3d axesPositions(4, 3);
     axesPositions.row(0) = Eigen::Vector3d(0, 0, 0);
-    //    axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
+    axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
-    //    axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
+    axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
-    //    axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
+    axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
-    axesPositions.row(1) = Eigen::Vector3d(1, 0, 0);
-    axesPositions.row(2) = Eigen::Vector3d(0, 1, 0);
-    axesPositions.row(3) = Eigen::Vector3d(0, 0, 1);
 
     Eigen::MatrixX2i axesEdges(3, 2);
     axesEdges.row(0) = Eigen::Vector2i(0, 1);
@@ -469,4 +347,17 @@ inline size_t computeHashOrdered(const std::vector<int> &v)
     return std::hash<std::string>{}(elementsString);
 }
 
+//inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
+//                                                      const std::string &extension)
+//{
+//    for (const std::filesystem::directory_entry &dirEntry :
+//         std::filesystem::directory_iterator(folderPath)) {
+//        if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
+//            return std::filesystem::path(dirEntry);
+//        }
+//    }
+
+//    return "";
+//}
+
 #endif // UTILITIES_H

vcgtrimesh.cpp

@@ -7,8 +7,8 @@
 #include <wrap/io_trimesh/import.h>
 //#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
-bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
+bool VCGTriMesh::load(const std::string &filename) {
-    assert(std::filesystem::exists(meshFilePath));
+    assert(std::filesystem::exists(filename));
     unsigned int mask = 0;
     mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
     mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
@@ -20,7 +20,7 @@ bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
     //            std::filesystem::absolute(filename).string().c_str(), *this, mask)
     //        != 0) {
     if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
-                                                 std::filesystem::absolute(meshFilePath).string().c_str())
+                                                 std::filesystem::absolute(filename).string().c_str())
         != 0) {
         std::cout << "Could not load tri mesh" << std::endl;
         return false;
@@ -32,7 +32,7 @@ bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
     vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
     vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
 
-    label = std::filesystem::path(meshFilePath).stem().string();
+    label = std::filesystem::path(filename).stem().string();
     return true;
 }
 

vcgtrimesh.hpp

@@ -39,7 +39,7 @@ class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
 public:
   VCGTriMesh();
   VCGTriMesh(const std::string &filename);
-  bool load(const std::filesystem::path &meshFilePath) override;
+  bool load(const std::string &filename) override;
   Eigen::MatrixX3d getVertices() const;
   Eigen::MatrixX3i getFaces() const;
   bool save(const std::string plyFilename);