diff --git a/src/static_sim/StaticSim.cpp b/src/static_sim/StaticSim.cpp
index 51bb387..d6cbad8 100644
--- a/src/static_sim/StaticSim.cpp
+++ b/src/static_sim/StaticSim.cpp
@@ -420,61 +420,6 @@ namespace ssim {
 
     }
 
-    void StaticSim::setBC() {
-        spdlog::info("set Boundary Conditions");
-
-        // DBC
-        int nDBC = 0;
-        DBC_nI.resize(nN);
-        isDBC.setZero(nN);
-        int DBCNum = (int) DirichletBCs.size();
-        for (int nI = 0; nI < nN; ++nI) {
-            Eigen::Vector3d p = TV.row(nI);
-            for (int _i = 0; _i < DBCNum; ++_i) {
-                if (DirichletBCs[_i].inDBC(p)) {
-                    DBC_nI(nDBC) = nI;
-                    isDBC(nI) = 1;
-                    ++nDBC;
-
-                    break;
-                }
-            }
-        }
-        DBC_nI.conservativeResize(nDBC);
-        // Utils::writeOBJ(outputPath + "DBCV.obj", TV(DBC_nI, Eigen::all),
-        //                 Eigen::VectorXi::LinSpaced(nDBC, 0, nDBC-1));
-        // NBC
-        load.resize(0);
-        load.setZero(nDof);
-        int nNBC = 0;
-        Eigen::VectorXi NBC_nI(nN);
-        int NBCNum = (int) NeumannBCs.size();
-        for (int nI = 0; nI < nN; ++nI) {
-            Eigen::Vector3d p = TV.row(nI);
-            for (int _i = 0; _i < NBCNum; ++_i) {
-                if (NeumannBCs[_i].inNBC(p)) {
-                    load.segment<DIM_>(nI * DIM_) = NeumannBCs[_i].force;
-                    NBC_nI(nNBC) = nI;
-                    ++nNBC;
-
-                    break;
-                }
-            }
-        }
-        NBC_nI.conservativeResize(nNBC);
-        // Utils::writeOBJ(outputPath + "NBCV.obj", TV(NBC_nI, Eigen::all),
-        //                 Eigen::VectorXi::LinSpaced(nNBC, 0, nNBC-1));
-
-        spdlog::info("#DBC nodes: {}, #NBC particles: {}", nDBC, nNBC);
-
-        // ensure (DBC intersect NBC) = (empty)
-        for (int i_ = 0; i_ < DBC_nI.size(); ++i_) {
-            int nI = DBC_nI(i_);
-            load.segment<DIM_>(nI * DIM_).setZero();
-        }
-    }
-
-    // BC only apply on surface nodes
     // void StaticSim::setBC() {
     //     spdlog::info("set Boundary Conditions");
     //
@@ -483,13 +428,12 @@ namespace ssim {
     //     DBC_nI.resize(nN);
     //     isDBC.setZero(nN);
     //     int DBCNum = (int) DirichletBCs.size();
-    //     for (int svI = 0; svI < SVI.size(); ++svI) {
-    //         int vI = SVI(svI);
-    //         Eigen::Vector3d p = TV.row(vI);
+    //     for (int nI = 0; nI < nN; ++nI) {
+    //         Eigen::Vector3d p = TV.row(nI);
     //         for (int _i = 0; _i < DBCNum; ++_i) {
     //             if (DirichletBCs[_i].inDBC(p)) {
-    //                 DBC_nI(nDBC) = vI;
-    //                 isDBC(vI) = 1;
+    //                 DBC_nI(nDBC) = nI;
+    //                 isDBC(nI) = 1;
     //                 ++nDBC;
     //
     //                 break;
@@ -505,13 +449,12 @@ namespace ssim {
     //     int nNBC = 0;
     //     Eigen::VectorXi NBC_nI(nN);
     //     int NBCNum = (int) NeumannBCs.size();
-    //     for (int svI = 0; svI < SVI.size(); ++svI) {
-    //         int vI = SVI(svI);
-    //         Eigen::Vector3d p = TV.row(vI);
+    //     for (int nI = 0; nI < nN; ++nI) {
+    //         Eigen::Vector3d p = TV.row(nI);
     //         for (int _i = 0; _i < NBCNum; ++_i) {
     //             if (NeumannBCs[_i].inNBC(p)) {
-    //                 load.segment<DIM_>(vI * DIM_) = NeumannBCs[_i].force;
-    //                 NBC_nI(nNBC) = vI;
+    //                 load.segment<DIM_>(nI * DIM_) = NeumannBCs[_i].force;
+    //                 NBC_nI(nNBC) = nI;
     //                 ++nNBC;
     //
     //                 break;
@@ -531,6 +474,63 @@ namespace ssim {
     //     }
     // }
 
+    // BC only apply on surface nodes
+    void StaticSim::setBC() {
+        spdlog::info("set Boundary Conditions");
+
+        // DBC
+        int nDBC = 0;
+        DBC_nI.resize(nN);
+        isDBC.setZero(nN);
+        int DBCNum = (int) DirichletBCs.size();
+        for (int svI = 0; svI < SVI.size(); ++svI) {
+            int vI = SVI(svI);
+            Eigen::Vector3d p = TV.row(vI);
+            for (int _i = 0; _i < DBCNum; ++_i) {
+                if (DirichletBCs[_i].inDBC(p)) {
+                    DBC_nI(nDBC) = vI;
+                    isDBC(vI) = 1;
+                    ++nDBC;
+
+                    break;
+                }
+            }
+        }
+        DBC_nI.conservativeResize(nDBC);
+        // Utils::writeOBJ(outputPath + "DBCV.obj", TV(DBC_nI, Eigen::all),
+        //                 Eigen::VectorXi::LinSpaced(nDBC, 0, nDBC-1));
+        // NBC
+        load.resize(0);
+        load.setZero(nDof);
+        int nNBC = 0;
+        Eigen::VectorXi NBC_nI(nN);
+        int NBCNum = (int) NeumannBCs.size();
+        for (int svI = 0; svI < SVI.size(); ++svI) {
+            int vI = SVI(svI);
+            Eigen::Vector3d p = TV.row(vI);
+            for (int _i = 0; _i < NBCNum; ++_i) {
+                if (NeumannBCs[_i].inNBC(p)) {
+                    load.segment<DIM_>(vI * DIM_) = NeumannBCs[_i].force;
+                    NBC_nI(nNBC) = vI;
+                    ++nNBC;
+
+                    break;
+                }
+            }
+        }
+        NBC_nI.conservativeResize(nNBC);
+        // Utils::writeOBJ(outputPath + "NBCV.obj", TV(NBC_nI, Eigen::all),
+        //                 Eigen::VectorXi::LinSpaced(nNBC, 0, nNBC-1));
+
+        spdlog::info("#DBC nodes: {}, #NBC particles: {}", nDBC, nNBC);
+
+        // ensure (DBC intersect NBC) = (empty)
+        for (int i_ = 0; i_ < DBC_nI.size(); ++i_) {
+            int nI = DBC_nI(i_);
+            load.segment<DIM_>(nI * DIM_).setZero();
+        }
+    }
+
     Model StaticSim::get_mesh() {
         if (mesh_.NumVertex() == 0) {
             // fill mesh_