Merge pull request #52 from dilevin/static_solve

fix data loading for muscle example, add indefinite fix
dilevin · Mar 29, 2019 · 959db62 · 959db62
2 parents b429249 + 9d0baa6
commit 959db62
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 7 deletions.
diff --git a/src/Examples/exampleMuscle.cpp b/src/Examples/exampleMuscle.cpp
@@ -24,12 +24,12 @@ typedef PhysicalSystemFEM<double, MuscleTet> FEMLinearTets;
 typedef World<double, std::tuple<FEMLinearTets *,PhysicalSystemParticleSingle<double> *>,
                       std::tuple<ForceSpringFEMParticle<double> *>,
                       std::tuple<ConstraintFixedPoint<double> *> > MyWorld;
-typedef TimeStepperEulerImplicit<double, AssemblerEigenSparseMatrix<double>,
+typedef TimeStepperStatic<double, AssemblerEigenSparseMatrix<double>,
 AssemblerEigenVector<double> > MyTimeStepper;
 
 typedef Scene<MyWorld, MyTimeStepper> MyScene;
 FEMLinearTets *test;
-double muscleStart = 10000;
+double muscleStart = 60000;
 Eigen::MatrixXd Ufibre;
 Eigen::VectorXi imuscle;
 
@@ -57,10 +57,10 @@ int main(int argc, char **argv) {
     Eigen::MatrixXi F;
 
 
-    igl::readDMAT("../../data/simple_muscle/generated_files/tet_mesh_T.dmat", F);
-    igl::readDMAT("../../data/simple_muscle/generated_files/tet_mesh_V.dmat", V);
-    igl::readDMAT("../../data/simple_muscle/generated_files/combined_fiber_directions.dmat", Ufibre);
-    igl::readDMAT("../../data/simple_muscle/generated_files/muscle_muscle_indices.dmat", imuscle);
+    igl::readDMAT(dataDir()+"/simple_muscle/generated_files/tet_mesh_T.dmat", F);
+    igl::readDMAT(dataDir()+"/simple_muscle/generated_files/tet_mesh_V.dmat", V);
+    igl::readDMAT(dataDir()+"/simple_muscle/generated_files/combined_fiber_directions.dmat", Ufibre);
+    igl::readDMAT(dataDir()+"/simple_muscle/generated_files/muscle_muscle_indices.dmat", imuscle);
     // readTetgen(V, F, dataDir()+"/meshesTetgen/Beam/Beam.node", dataDir()+"/meshesTetgen/Beam/Beam.ele");
 
     test = new FEMLinearTets(V,F);
@@ -99,4 +99,4 @@ int main(int argc, char **argv) {
     GAUSSVIEW(scene);
 
     return app.exec();
-}
+}
diff --git a/src/FEM/include/EnergyMuscle.h b/src/FEM/include/EnergyMuscle.h
@@ -493,6 +493,21 @@ class EnergyMuscle : public virtual ShapeFunction {
             -gradZ.transpose()*ddw.block(6,6,3,3)*gradZ;
 
 
+        // hard coded for tet, need to change size for hex
+       Eigen::SelfAdjointEigenSolver<Matrix> es(-H);
+
+        Eigen::MatrixXd DiagEval = es.eigenvalues().real().asDiagonal();
+        Eigen::MatrixXd Evec = es.eigenvectors().real();
+
+        for (int i = 0; i < 12; ++i) {
+            if (es.eigenvalues()[i]<1e-6) {
+                DiagEval(i,i) = 1e-3;
+            }
+        }
+        //        saveMarket(H, "H.dat");
+        H = -Evec * DiagEval * Evec.transpose();
+        //
+
     }
 
     template<typename Matrix>