fix iks003 and iks007

ikarus-project · Nov 9, 2023 · 58d4d5d · 58d4d5d
1 parent 6e0f4be
commit 58d4d5d
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 11 deletions.
diff --git a/src/iks003_incompressibleLinearElasticity.cpp b/src/iks003_incompressibleLinearElasticity.cpp
@@ -192,7 +192,6 @@ int main(int argc, char **argv) {
 
   /// Create assembler
   auto sparseFlatAssembler = SparseFlatAssembler(fes, dirichletValues);
-  auto denseFlatAssembler  = DenseFlatAssembler(fes, dirichletValues);
 
   /// Create function for external forces and stiffness matrix
   double lambda = 0;
@@ -204,7 +203,7 @@ int main(int argc, char **argv) {
   auto fextFunction = [&](auto &&lambdaLocal, auto &&dLocal) -> auto & {
     req.insertGlobalSolution(Ikarus::FESolutions::displacement, dLocal)
         .insertParameter(Ikarus::FEParameter::loadfactor, lambdaLocal);
-    return denseFlatAssembler.getReducedVector(req);
+    return sparseFlatAssembler.getReducedVector(req);
   };
   auto KFunction = [&](auto &&lambdaLocal, auto &&dLocal) -> auto & {
     req.insertGlobalSolution(Ikarus::FESolutions::displacement, dLocal)
@@ -218,7 +217,7 @@ int main(int argc, char **argv) {
   ld.compute(K);
   if (ld.info() != Eigen::Success) DUNE_THROW(Dune::MathError, "Failed Compute");
 
-  d -= denseFlatAssembler.createFullVector(ld.solve(R));
+  d -= sparseFlatAssembler.createFullVector(ld.solve(R));
   if (ld.info() != Eigen::Success) DUNE_THROW(Dune::MathError, "Failed Solve");
 
   /// Postprocess

diff --git a/src/iks004_kirchhoffPlate.cpp b/src/iks004_kirchhoffPlate.cpp
@@ -285,8 +285,6 @@ int main(int argc, char **argv) {
       const auto &rule = Dune::QuadratureRules<double, 2>::rule(
           ele.type(), 2 * localView.tree().finiteElement().localBasis().order());
       for (auto gp : rule) {
-        const auto gpGlobalPos = geo.global(gp.position());
-
         const auto w_ex = localwAna(gp.position());
         const auto w_fe = localw(gp.position());
         l2_error += Dune::power(w_ex - w_fe, 2) * ele.geometry().integrationElement(gp.position()) * gp.weight();

diff --git a/src/iks007_vonMisesTruss.cpp b/src/iks007_vonMisesTruss.cpp
@@ -88,7 +88,7 @@ class Truss : public PowerBasisFE<typename Basis_::FlatBasis> {
 
     const ScalarType Egl = 0.5 * (lsquared - LRefsquared) / LRefsquared;
 
-    return 0.5 * EA / sqrt(LRefsquared) * Egl * Egl;
+    return 0.5 * EA * sqrt(LRefsquared) * Egl * Egl;
   }
 
  private:
@@ -100,7 +100,7 @@ int main(int argc, char **argv) {
   /// Construct grid
   Dune::GridFactory<Dune::FoamGrid<1, 2, double>> gridFactory;
   const double h = 1.0;
-  const double L = 1.0;
+  const double L = 2.0;
   gridFactory.insertVertex({0, 0});
   gridFactory.insertVertex({L, h});
   gridFactory.insertVertex({2 * L, 0});
@@ -200,17 +200,18 @@ int main(int argc, char **argv) {
 
   auto analyticalLoadDisplacementCurve = [&](auto &w) {
     const double Ltruss = std::sqrt(h * h + L * L);
-    return EA * Dune::power(h, 3) / Dune::power(Ltruss, 3)
+    return 2.0 * EA * Dune::power(h, 3) / Dune::power(Ltruss, 3)
            * (w / h - 1.5 * Dune::power(w / h, 2) + 0.5 * Dune::power(w / h, 3));
   };
 
   std::vector<double> x  = linspace(0.0, dVec.maxCoeff());
   std::vector<double> y1 = transform(x, [&](auto x) { return analyticalLoadDisplacementCurve(x); });
   auto p                 = plot(x, y1, dVec, lambdaVec);
-  p[0]->line_width(2);
+  p[0]->line_width(3);
   p[1]->line_width(2);
   p[1]->marker(line_spec::marker_style::asterisk);
+  // save("vonMisesTruss.png");
   // f->draw();
-  //  using namespace std::chrono_literals;
-  //  std::this_thread::sleep_for(5s);
+  // using namespace std::chrono_literals;
+  // std::this_thread::sleep_for(5s);
 }