diff --git a/CMakeLists.txt b/CMakeLists.txt
index 66cc9dd..83adc0b 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -48,10 +48,10 @@ if(${USE_POLYSCOPE})
         add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
         add_compile_definitions(POLYSCOPE_DEFINED)
         target_sources(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
-        if(${USE_POLYSCOPE})
-                message("Using polyscope here")
-                target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
-        endif()
+
+        message("Using polyscope..")
+        target_include_directories(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/include)
+        target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
 endif()
 
 ##vcglib devel branch
@@ -62,9 +62,6 @@ download_project(PROJ		    vcglib_devel
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
 add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
-
-#add_subdirectory("/home/iason/Coding/Libraries/vcglib" ${CMAKE_CURRENT_BINARY_DIR}/vcglib)
-#target_include_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/vcglib")
 target_sources(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
 
 ##matplot++ lib
@@ -90,27 +87,11 @@ download_project(PROJ                threed-beam-fea
 )
 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)
+find_package(Eigen3 3.4 REQUIRED)
 if(MSVC)
         add_compile_definitions(_HAS_STD_BYTE=0)
 endif(MSVC)
-#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
 
-#set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
 if(${MYSOURCES_STATIC_LINK})
         message("STATIC LINK MY SOURCES:" ${MYSOURCES_STATIC_LINK})
 endif()
@@ -149,33 +130,56 @@ download_project(PROJ			ENSMALLEN
                  PREFIX             	${EXTERNAL_DEPS_DIR}
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
-#        add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
-#        target_link_libraries(${PROJECT_NAME} INTERFACE ensmallen)
 target_include_directories(${PROJECT_NAME}
 PUBLIC ${ENSMALLEN_SOURCE_DIR}/include)
 add_compile_definitions(USE_ENSMALLEN)
 #endif()
 
 #Chrono
-#set(Chrono_DIR "/home/iason/Coding/Libraries/chrono-7.0.3/build/cmake")
 set(Chrono_DIR "/home/iason/Coding/Libraries/chrono-7.0.3/build/cmake")
-LIST(APPEND CMAKE_PREFIX_PATH "${CMAKE_INSTALL_PREFIX}/../Chrono/lib")
-#set(Chrono_DIR "/home/iason/Coding/build/external dependencies/CHRONO-src/cmake")
-find_package(Chrono #[[COMPONENTS FEA]]
-                 CONFIG)
-if (NOT Chrono_FOUND)
-        message("Could not find Chrono or one of its required modules")
-        return()
+#LIST(APPEND CMAKE_PREFIX_PATH "${CMAKE_INSTALL_PREFIX}/../Chrono/lib")
+#find_package(Chrono #[[COMPONENTS FEA]]
+#                 CONFIG)
+#if (NOT Chrono_FOUND)
+#        message("Could not find Chrono or one of its required modules")
+#        return()
+#endif()
+#set_target_properties(${PROJECT_NAME} PROPERTIES
+#            COMPILE_FLAGS "${CHRONO_CXX_FLAGS} ${EXTRA_COMPILE_FLAGS}"
+#            COMPILE_DEFINITIONS "CHRONO_DATA_DIR=\"${CHRONO_DATA_DIR}\""
+#            LINK_FLAGS "${CHRONO_LINKER_FLAGS}")
+#target_link_libraries(${PROJECT_NAME} PRIVATE ${CHRONO_LIBRARIES})
+target_include_directories(${PROJECT_NAME} PRIVATE ${CHRONO_INCLUDE_DIRS} #[["/home/iason/Coding/Libraries/chrono-7.0.3/src"]] 	#[["/usr/include/irrlicht"]] )
+#target_link_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/chrono-7.0.3/build/Debug/lib/")
+target_link_libraries(${PROJECT_NAME} PRIVATE "/home/iason/Coding/Libraries/chrono-7.0.3/build/Debug/lib/libChronoEngine.so" #[[ChronoEngine_irrlicht]])
+
+#add_subdirectory("/home/iason/Coding/Projects/Approximating shapes with flat patterns/Elastic rods/ElasticRods")
+#set(MESHFEM_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/meshfem)
+#download_project(PROJ			MESHFEM
+#                 GIT_REPOSITORY		https://github.com/MeshFEM/MeshFEM.git
+#                 GIT_TAG		3847ffaeb6e9c4bbe15f0f81f5c773c7c481c059
+#                 BINARY_DIR		${MESHFEM_BINARY_DIR}
+#                 PREFIX             	${EXTERNAL_DEPS_DIR}
+#                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+#)
+#add_subdirectory(${MESHFEM_SOURCE_DIR} ${MESHFEM_BINARY_DIR})
+#target_link_libraries(${PROJECT_NAME} PUBLIC MeshFEM)
+##target_include_directories(${PROJECT_NAME}
+##PUBLIC ${MESHFEM_SOURCE_DIR}/include)
+##add_compile_definitions(USE_ENSMALLEN)
+
+##TBB
+if(NOT TARGET 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})
 endif()
-set_target_properties(${PROJECT_NAME} PROPERTIES
-            COMPILE_FLAGS "${CHRONO_CXX_FLAGS} ${EXTRA_COMPILE_FLAGS}"
-            COMPILE_DEFINITIONS "CHRONO_DATA_DIR=\"${CHRONO_DATA_DIR}\""
-            LINK_FLAGS "${CHRONO_LINKER_FLAGS}")
-target_link_libraries(${PROJECT_NAME} PRIVATE ${CHRONO_LIBRARIES})
-target_include_directories(${PROJECT_NAME} PRIVATE ${CHRONO_INCLUDE_DIRS})
+link_directories(${TBB_BINARY_DIR})
 
-#find_package(OpenMP REQUIRED)
-#target_link_libraries(${PROJECT_NAME} OpenMP::OpenMP_CXX)
-#target_include_directories(${PROJECT_NAME} OpenMP::OpenMP_CXX)
-
-#message(STATUS "OpenMP_C_INCLUDE_DIRS: ${OpenMP_CXX_INCLUDE_DIRS}")
diff --git a/der_leimer.cpp b/der_leimer.cpp
index 48e293c..c43d606 100644
--- a/der_leimer.cpp
+++ b/der_leimer.cpp
@@ -2,6 +2,7 @@
 #include <unordered_set>
 
 #define MAT(mat, m, row, col) mat[m * col + row]
+//#define EXPLICIT_BC
 
 DER_leimer::DER_leimer() {}
 
@@ -255,11 +256,12 @@ SimulationResults DER_leimer::constructSimulationResults(const std::shared_ptr<S
     //    pJob->pMesh->registerForDrawing();
     //    polyscope::show();
     for (int vi = 0; vi < config->numVertices(); vi++) {
-        std::cout << vi << " pos:" << config->vertices(vi, 0) << " " << config->vertices(vi, 1)
-                  << " " << config->vertices(vi, 2) << std::endl;
+        //        std::cout << vi << " pos:" << config->vertices(vi, 0) << " " << config->vertices(vi, 1)
+        //                  << " " << config->vertices(vi, 2) << std::endl;
         for (int dofi = 0; dofi < 3; dofi++) {
             const double initialPos_dofi = pJob->pMesh->vert[vi].cP()[dofi];
-            const double displacement_dofi = config->vertices(vi, dofi) - initialPos_dofi;
+            const double deformedPos_dofi = config->vertices(vi, dofi);
+            const double displacement_dofi = deformedPos_dofi - initialPos_dofi;
             simulationResults.displacements[vi][dofi] = displacement_dofi;
         }
         for (int dofi = 3; dofi < 6; dofi++) {
@@ -353,7 +355,7 @@ SimulationResults DER_leimer::executeSimulation(const std::shared_ptr<Simulation
     //    std::vector<Anchor> anchors{a0};
     const std::shared_ptr<VCGEdgeMesh> pMesh = pJob->pMesh;
     const int numRods = pMesh->EN();
-    const int numVertsPerRod = 2;
+    constexpr int numVertsPerRod = 2;
     std::vector<int> vInd;
     vInd.reserve(pMesh->EN());
     for (int ei = 0; ei < pMesh->EN(); ei++) {
@@ -391,12 +393,12 @@ SimulationResults DER_leimer::executeSimulation(const std::shared_ptr<Simulation
         a.bary = edgeVi;
         a.position = ep->cP(edgeVi).ToEigenVector<Eigen::Vector3d>();
         a.normal = ep->cV(edgeVi)->cN().ToEigenVector<Eigen::Vector3d>();
-        a.k_pos = 1e4;
+        a.k_pos = 1e6;
         if (bc.second == fixedBC) {
             a.k_dir = 0;
         } else {
             assert(bc.second == rigidBC);
-            a.k_dir = 1e4;
+            a.k_dir = 1e6;
         }
         anchors.push_back(a);
     }
@@ -569,7 +571,7 @@ bool DER_leimer::simulateOneStepGrid(RodConfig *config)
 
     iter++;
 
-    if (iter > 50 /*&& forceResidual < 1e-5*/) {
+    if (iter > 100 /*|| solverResidual < 1e-10*/ || forceResidual < 1e-5) {
         //        for (int vi = 0; vi < config->numVertices(); vi++) {
         //            std::cout << vi << " pos:" << config->vertices(vi, 0) << " " << config->vertices(vi, 1)
         //                      << " " << config->vertices(vi, 2) << std::endl;
@@ -594,11 +596,11 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                            std::vector<Eigen::Triplet<int>> *Hu_index,
                            std::vector<double> *Hu_value)
 {
-    const bool computeHessian = [&]() {
-        if (H_index && H_value && Ju && Hu_index && Hu_value)
-            return true;
-        return false;
-    }();
+    //    const bool computeHessian = [&]() {
+    //        if (H_index && H_value && Ju && Hu_index && Hu_value)
+    //            return true;
+    //        return false;
+    //    }();
 
     const bool hasExternalForces = externalForces.rows() != 0;
 
@@ -624,8 +626,10 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
         nterms += config.rods[i]->numSegments(); // stretching
         nterms += 2 * njoints;                   // bending
         nterms += njoints;                       // twisting
-        if (params.anchorPoints)                 //??
+#ifndef EXPLICIT_BC
+        if (params.anchorPoints) //??
             nterms += 3 * config.rods[i]->numVertices();
+#endif
         //        if (hasExternalForces) {
         //            nterms += config.rods[i]->numVertices();
         //        }
@@ -635,12 +639,14 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
     }
 
     ndofs += config.numVertices() * 3; // DOF for vertices
-                                       //if (params.anchorPoints)
-    //	nterms += 3 * config.numVertices();
-    //if (params.gravityEnabled)
-    //	nterms += config.numVertices(); // floor force
+    if (params.anchorPoints)
+        nterms += 3 * config.numVertices();
+        //if (params.gravityEnabled)
+        //	nterms += config.numVertices(); // floor force
 
+#ifndef EXPLICIT_BC
     nterms += 4 * config.numAnchors();
+#endif
 
     // Barycentric Coordinates are handled as DOF's since they can change when rods are sliding
     int baryoffset = ndofs;
@@ -669,8 +675,8 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
     ndofs += 3 * config.numIntersections(); // euler angles as DOFs for intersections
 
     r.resize(nterms);
-    r.setConstant(std::numeric_limits<double>::infinity());
-    //        r.setZero();
+    //    r.setConstant(std::numeric_limits<double>::infinity());
+    r.setZero();
 
     gravityPE = 0;
     gravityForces.resize(ndofs);
@@ -864,16 +870,25 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                 //                const Eigen::Vector3d debug_forceLeimer = jacobianEntry * segr;
                 //                assert(debug_forceLeimer == debug_forcePanetta_j);
                 for (int j = 0; j < 3; j++) {
-                    // Gradients for dEs/dx_i and dEs/dx_(i+1)
-                    //                    if (!pJob->constrainedVertices.contains(m1)) {
-                    J.push_back(
-                        Eigen::Triplet<double>(roffset + si_rod, 3 * m1 + j, jacobianEntry[j]));
-                    //                    }
-                    //                    if (!pJob->constrainedVertices.contains(m2)) {
-                    J.push_back(
-                        Eigen::Triplet<double>(roffset + si_rod, 3 * m2 + j, -jacobianEntry[j]));
-                    //                    }
+// Gradients for dEs/dx_i and dEs/dx_(i+1)
+#ifdef EXPLICIT_BC
+                    if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                        J.push_back(
+                            Eigen::Triplet<double>(roffset + si_rod, 3 * m1 + j, jacobianEntry[j]));
+
+#ifdef EXPLICIT_BC
+                    }
+                    if (!pJob->constrainedVertices.contains(m2)) {
+#endif
+                        J.push_back(
+                            Eigen::Triplet<double>(roffset + si_rod, 3 * m2 + j, -jacobianEntry[j]));
+#ifdef EXPLICIT_BC
+                    }
+#endif
                 }
+
+#ifdef COMPUTE_HESSIAN
                 if (computeHessian) {
                     Eigen::Matrix3d Hlen = (1 / len) * Eigen::Matrix<double, 3, 3>::Identity()
                                            - (1 / (len * len * len)) * (v1 - v2)
@@ -945,6 +960,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                                                uistiffoffset + 4 * nsegsTotal + dofoffset + si_rod,
                                                df_dt * (len - restlen)));
                 }
+#endif
             }
         }
         roffset += rod.numSegments();
@@ -1067,18 +1083,28 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                     //                    const Eigen::Vector3d debug_force1Leimer_xip1 = dEb1_dxip1 * debug_vertr1;
                     for (int j = 0; j < 3; j++) {
                         // Gradients for dEb1/dx_(i-1), dEb1/dx_i and dEb1/dx_(i+1)
-                        //                        if (!pJob->constrainedVertices.contains(m0)) {
-                        J.push_back(
-                            Eigen::Triplet<double>(roffset + 2 * idx, 3 * m0 + j, dEb1_dxim1[j]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m1)) {
-                        J.push_back(
-                            Eigen::Triplet<double>(roffset + 2 * idx, 3 * m1 + j, dEb1_dxi[j]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m2)) {
-                        J.push_back(
-                            Eigen::Triplet<double>(roffset + 2 * idx, 3 * m2 + j, dEb1_dxip1[j]));
-                        //                        }
+#ifdef EXPLICIT_BC
+                        if (!pJob->constrainedVertices.contains(m0)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + 2 * idx,
+                                                               3 * m0 + j,
+                                                               dEb1_dxim1[j]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                            J.push_back(
+                                Eigen::Triplet<double>(roffset + 2 * idx, 3 * m1 + j, dEb1_dxi[j]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m2)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + 2 * idx,
+                                                               3 * m2 + j,
+                                                               dEb1_dxip1[j]));
+#ifdef EXPLICIT_BC
+                        }
+#endif
                     }
                     // Gradients for dk2/dx_(i-1) and dk2/dx_i
                     Eigen::Vector3d dk21 = 1.0 / (v1 - v0).norm()
@@ -1100,20 +1126,29 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                     //                    const Eigen::Vector3d debug_force2Leimer_xip1 = dEb2_dxip1 * debug_vertr2;
                     for (int j = 0; j < 3; j++) {
                         // Gradients for dEb2/dx_(i-1), dEb2/dx_i and dEb2/dx_(i+1)
-                        //                        if (!pJob->constrainedVertices.contains(m0)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
-                                                           3 * m0 + j,
-                                                           dEb2_dxim1[j]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m1)) {
-                        J.push_back(
-                            Eigen::Triplet<double>(roffset + 2 * idx + 1, 3 * m1 + j, dEb2_dxi[j]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m2)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
-                                                           3 * m2 + j,
-                                                           dEb2_dxip1[j]));
-                        //                        }
+#ifdef EXPLICIT_BC
+                        if (!pJob->constrainedVertices.contains(m0)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
+                                                               3 * m0 + j,
+                                                               dEb2_dxim1[j]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
+                                                               3 * m1 + j,
+                                                               dEb2_dxi[j]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m2)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
+                                                               3 * m2 + j,
+                                                               dEb2_dxip1[j]));
+#ifdef EXPLICIT_BC
+                        }
+#endif
                     }
 
                     // Material curvature gradients dd/dtheta1 and dd/theta2
@@ -1146,6 +1181,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                                                0.5 * sqrt(factor2 /* * width * width * width*/)
                                                    * Dd12.dot(kb)));
 
+#ifdef COMPUTE_HESSIAN
                     if (computeHessian) {
                         double chi = 1 + t01.dot(t12);
                         double len01 = (v1 - v0).norm();
@@ -2014,6 +2050,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                             int stop = 0;
                         }
                     }
+#endif
                 }
             } else {
                 r[roffset + 2 * idx] = 0;
@@ -2083,19 +2120,27 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                     const Eigen::Vector3d debug_leimerTwistForce_pos2 = jacob_theta2
                                                                         * r[roffset + idx];
                     for (int j = 0; j < 3; j++) {
-                        //                        if (!pJob->constrainedVertices.contains(m0)) {
-                        J.push_back(
-                            Eigen::Triplet<double>(roffset + idx, 3 * m0 + j, -jacob_theta1[j]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m1)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + idx,
-                                                           3 * m1 + j,
-                                                           jacob_theta1[j] - jacob_theta2[j]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m2)) {
-                        J.push_back(
-                            Eigen::Triplet<double>(roffset + idx, 3 * m2 + j, jacob_theta2[j]));
-                        //                        }
+#ifdef EXPLICIT_BC
+                        if (!pJob->constrainedVertices.contains(m0)) {
+#endif
+                            J.push_back(
+                                Eigen::Triplet<double>(roffset + idx, 3 * m0 + j, -jacob_theta1[j]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + idx,
+                                                               3 * m1 + j,
+                                                               jacob_theta1[j] - jacob_theta2[j]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m2)) {
+#endif
+                            J.push_back(
+                                Eigen::Triplet<double>(roffset + idx, 3 * m2 + j, jacob_theta2[j]));
+#ifdef EXPLICIT_BC
+                        }
+#endif
                     }
                     // Gradients dEt/dtheta1 and dEt/dtheta2
                     const double debug_leimerTwistForce_theta1 = sqrt(debug_factor)
@@ -2111,6 +2156,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                                                        dofoffset + thetaoffset + si_rod,
                                                        sqrt(debug_factor)));
 
+#ifdef COMPUTE_HESSIAN
                     if (computeHessian) {
                         double chi = 1 + t01.dot(t12);
                         double len01 = (v1 - v0).norm();
@@ -2374,6 +2420,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                             int stop = 0;
                         }
                     }
+#endif
                 }
             } else
                 r[roffset + idx] = 0;
@@ -2384,6 +2431,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
         dofoffset += rod.numSegments();
     }
 
+#ifndef EXPLICIT_BC
     // mesh sliding energy Eslide
     // allows ribbons to slide tangentially on mesh surface they try to approximate
     // see paper Section 7.1
@@ -2491,11 +2539,11 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
             r[roffset + idx + 2] = posStiffness * disp[2];
 
             /*
-    			std::cout << "Anchor Nr. " << i << std::endl;
-    			std::cout << "\tr[" << (roffset + idx + 0) << "] = " << r[roffset + idx + 0] << std::endl;
-    			std::cout << "\tr[" << (roffset + idx + 1) << "] = " << r[roffset + idx + 1] << std::endl;
-    			std::cout << "\tr[" << (roffset + idx + 2) << r[roffset + idx + 2] << std::endl;
-    			*/
+        			std::cout << "Anchor Nr. " << i << std::endl;
+        			std::cout << "\tr[" << (roffset + idx + 0) << "] = " << r[roffset + idx + 0] << std::endl;
+        			std::cout << "\tr[" << (roffset + idx + 1) << "] = " << r[roffset + idx + 1] << std::endl;
+        			std::cout << "\tr[" << (roffset + idx + 2) << r[roffset + idx + 2] << std::endl;
+        			*/
 
             int m1 = rod->curState.centerlineIndices[a.seg] * 3;
             int m2 = rod->curState.centerlineIndices[a.seg + 1] * 3;
@@ -2566,6 +2614,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                     }
                 }
 
+#ifdef COMPUTE_HESSIAN
                 if (computeHessian) {
                     // design variables - Ui anchor positions
                     for (int j = 0; j < 3; j++) {
@@ -2637,6 +2686,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                         int stop = 0;
                     }
                 }
+#endif
             }
 
             // direction constraint for anchor
@@ -2678,6 +2728,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                                                    rod1offset + a.seg,
                                                    sqrt(factor) * dtheta1.dot(d1dtheta)));
 
+#ifdef COMPUTE_HESSIAN
                 if (computeHessian) {
                     double len = (vSeg1 - vSeg0).norm();
                     Eigen::Matrix3d dt_de1 = ((1 / len) * Eigen::Matrix<double, 3, 3>::Identity()
@@ -2971,6 +3022,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                     }
                     int stop = 0;
                 }
+#endif
             }
 
             idx += 4;
@@ -2978,6 +3030,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
 
         roffset += idx;
     }
+#endif
 
     // compute energy terms for rod intersections
     int incCount = 0;
@@ -3119,16 +3172,22 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                                            * (-k1 * ttilde + st01.cross(d2tilde));
                     for (int k = 0; k < 3; k++) {
                         // Gradients for dEb1/dx_(i-1), dEb1/dx_i and dEb1/dx_(i+1)
-                        //                        if (!pJob->constrainedVertices.contains(m0)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + incCount,
-                                                           m0 * 3 + k,
-                                                           -sqrt(factor1) * dk11[k]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m1)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + incCount,
-                                                           m1 * 3 + k,
-                                                           sqrt(factor1) * dk11[k]));
-                        //                        }
+#ifdef EXPLICIT_BC
+                        if (!pJob->constrainedVertices.contains(m0)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + incCount,
+                                                               m0 * 3 + k,
+                                                               -sqrt(factor1) * dk11[k]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + incCount,
+                                                               m1 * 3 + k,
+                                                               sqrt(factor1) * dk11[k]));
+#ifdef EXPLICIT_BC
+                        }
+#endif
                     }
                     J.push_back(Eigen::Triplet<double>(roffset + incCount,
                                                        euleroffset + i * 3 + 0,
@@ -3145,16 +3204,22 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                                            * (-k2 * ttilde + st01.cross(d1tilde));
                     for (int k = 0; k < 3; k++) {
                         // Gradients for dEb2/dx_(i-1), dEb2/dx_i and dEb2/dx_(i+1)
-                        //                        if (!pJob->constrainedVertices.contains(m0)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
-                                                           m0 * 3 + k,
-                                                           -sqrt(factor2) * dk21[k]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m1)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
-                                                           m1 * 3 + k,
-                                                           sqrt(factor2) * dk21[k]));
-                        //                        }
+#ifdef EXPLICIT_BC
+                        if (!pJob->constrainedVertices.contains(m0)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
+                                                               m0 * 3 + k,
+                                                               -sqrt(factor2) * dk21[k]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
+                                                               m1 * 3 + k,
+                                                               sqrt(factor2) * dk21[k]));
+#ifdef EXPLICIT_BC
+                        }
+#endif
                     }
                     J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
                                                        euleroffset + i * 3 + 0,
@@ -3183,16 +3248,22 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                     Eigen::Vector3d dpsi_dx = 0.5 * kb / (cv1 - cv0).norm();
                     for (int k = 0; k < 3; k++) {
                         // Gradients dEt/dx_(i-1), dEt/dx_i and dEt/dx_(i+1)
-                        //                        if (!pJob->constrainedVertices.contains(m0)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
-                                                           m0 * 3 + k,
-                                                           -sqrt(factor3) * dpsi_dx[k]));
-                        //                        }
-                        //                        if (!pJob->constrainedVertices.contains(m1)) {
-                        J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
-                                                           m1 * 3 + k,
-                                                           sqrt(factor3) * dpsi_dx[k]));
-                        //                        }
+#ifdef EXPLICIT_BC
+                        if (!pJob->constrainedVertices.contains(m0)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
+                                                               m0 * 3 + k,
+                                                               -sqrt(factor3) * dpsi_dx[k]));
+#ifdef EXPLICIT_BC
+                        }
+                        if (!pJob->constrainedVertices.contains(m1)) {
+#endif
+                            J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
+                                                               m1 * 3 + k,
+                                                               sqrt(factor3) * dpsi_dx[k]));
+#ifdef EXPLICIT_BC
+                        }
+#endif
                     }
                     // Gradients dEt/dtheta1 and dEt/dtheta2
                     J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
@@ -3283,6 +3354,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
                         }
                     }
 
+#ifdef COMPUTE_HESSIAN
                     if (computeHessian) {
                         double chi = 1 + ct01.dot(st01);
                         double len01 = (cv1 - cv0).norm();
@@ -4078,6 +4150,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
 						H_value->push_back(-sqrt(factor3*width) * dpsi22);
 						*/
                     }
+#endif
                 }
                 incCount += 3;
             }
@@ -4089,12 +4162,31 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
         exit(-1);
 
     if (Jr) {
+        //        std::vector<Eigen::Triplet<double>> J_nonConstrained;
+        //        std::unordered_set<int> constrainedDof_global;
+        //        for (const auto &constrainedViDof : pJob->constrainedVertices) {
+        //            const int vi = constrainedViDof.first;
+        //            std::unordered_set<int> &constrainedDof_local = constrainedViDof.second;
+
+        //        }
+        //        for (auto &triplet : J) {
+        //        }
         Jr->setFromTriplets(J.begin(), J.end());
     }
+#ifdef COMPUTE_HESSIAN
     if (computeHessian) {
         Ju->resize(nterms, udofs);
         Ju->setFromTriplets(Jui.begin(), Jui.end());
     }
+#endif
+}
+
+Eigen::Vector3d DER_leimer::getPerpendicularVector(const Eigen::Vector3d &t) const
+{
+    Eigen::Vector3d candidate1 = Eigen::Vector3d(1, 0, 0).cross(t),
+                    candidate2 = Eigen::Vector3d(0, 1, 0).cross(t);
+    return (candidate1.norm() > candidate2.norm()) ? candidate1.normalized()
+                                                   : candidate2.normalized();
 }
 
 DER_leimer::RodConfig *DER_leimer::readRodGrid(
@@ -4203,29 +4295,8 @@ DER_leimer::RodConfig *DER_leimer::readRodGrid(
             const Eigen::Vector3d v1 = rodState.centerlinePositions.row(
                 (si_rod + 1) % numVertices_rod); //the % numVertices_rod is not needed
             const Eigen::Vector3d e = (v1 - v0).normalized();
-            Eigen::Vector3d segmentDirector = Eigen::Vector3d(0, 0, 1);
-            /*= ((pMesh->edge[rodIndex].cV(0)->cN() + pMesh->edge[rodIndex].cV(1)->cN()) / 2)
-                      .ToEigenVector<Eigen::Vector3d>()
-                      .normalized();*/
-            //                = (e.cross(Eigen::Vector3d(0, 0, 1).normalized())).normalized();
-            //            switch (rodIndex) {
-            //            case 0:
-            //                segmentDirector = Eigen::Vector3d(1, 0, 0);
-            //                break;
-
-            //            case 1:
-            //                segmentDirector = Eigen::Vector3d(0, -1, 0);
-            //                break;
-
-            //            case 2:
-            //                segmentDirector = Eigen::Vector3d(-1, 0, 0);
-            //                segmentDirector = Eigen::Vector3d(-1, 0, 0);
-            //                break;
-            //            default:
-            //                assert(false);
-            //            }
-
-            rodState.directors.row(si_rod) = segmentDirector;
+            //            rodState.directors.row(si_rod) = getPerpendicularVector(e);
+            rodState.directors.row(si_rod) = Eigen::Vector3d(0, 0, 1);
 
             //            rodState.directors.row(si_rod) = binorms[si_rod];
             //            for (int k = 0; k < 3; k++) {
@@ -4745,3 +4816,10 @@ int DER_leimer::RodConfig::getNumDoF() const
     ndofs += 3 * numIntersections(); // euler angles as DOFs for intersections
     return ndofs;
 }
+
+void DER_leimer::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
+{
+    std::cout << "This function is currently not implemented" << std::endl;
+    assert(false);
+    std::terminate();
+}
diff --git a/der_leimer.hpp b/der_leimer.hpp
index 17b0199..015c383 100644
--- a/der_leimer.hpp
+++ b/der_leimer.hpp
@@ -139,6 +139,8 @@ public:
         int getNumDoF() const;
     };
 
+    void setStructure(const std::shared_ptr<SimulationMesh> &pMesh) override;
+
 private:
     bool simulateOneStepGrid(RodConfig *config);
 
@@ -181,6 +183,7 @@ private:
     SimulationResults constructSimulationResults(const std::shared_ptr<SimulationJob> &pJob,
                                                  RodConfig *config) const;
     std::shared_ptr<SimulationJob> pJob;
+    Eigen::Vector3d getPerpendicularVector(const Eigen::Vector3d &t) const;
 };
 
 Eigen::Matrix3d eulerRotation(double e0, double e1, double e2, int deriv = -1, int deriv2 = -1);