WIP

dune/elasticity/common/trustregionsolver.cc

@@ -75,7 +75,7 @@ setup(const typename BasisType::GridView::Grid& grid,
     std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
       Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
         BasisType,
-        DuneFunctionsBasis<BasisType> > basis(assembler_.basis);
+        DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
 
     innerMultigridStep_ = std::make_unique<Dune::ParMG::Multigrid<VectorType> >();
     innerMultigridStep_->mu_ = mu;

problems/quasiconvexity-circle-expacosh.parset

@@ -12,7 +12,7 @@ gridFile = circle2ndorder.msh
 #elements = 2 2
 
 # Number of grid levels
-numLevels = 6
+numLevels = 8
 
 #############################################
 #  Solver parameters

src/CMakeLists.txt

@@ -11,8 +11,8 @@ if(ADOLC_FOUND AND IPOPT_FOUND AND PYTHONLIBS_FOUND AND dune-uggrid_FOUND)
     add_dune_ug_flags(${_program})
     add_dune_mpi_flags(${_program})
     add_dune_superlu_flags(${_program})
-    target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
-#    target_compile_options(${_program} PRIVATE "-g")
+#    target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
+    target_compile_options(${_program} PRIVATE "-g")
   endforeach()
 
   if (AMIRAMESH_FOUND)

src/quasiconvexity-test.cc

@@ -83,6 +83,35 @@ struct BurkholderDensity final
 };
 
 
+template<int dim, class field_type, class ctype=double>
+struct ExpACosHDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  ExpACosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    c1_ = parameters.get<double>("c1");
+    c2_ = parameters.get<double>("c2");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*det);
+
+    using std::acosh;
+    return exp(c1_ * acosh(K) * acosh(K)) - c2_ * log(det);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c1_;
+  double c2_;
+};
+
+
 template<int dim, class field_type, class ctype=double>
 struct KPowerDensity final
 : public Elasticity::LocalDensity<dim,field_type,ctype>
@@ -374,7 +403,7 @@ int main (int argc, char *argv[]) try
   using FEBasis = Fufem::PeriodicBasis<GridView>;
   FEBasis feBasis(gridView);
 
-#if 0
+#if 1
   // This is needed for interpolation, but PeriodicBasis
   // doesn't quite support being put into power bases yet.
   using namespace Dune::Fufem::BasisFactory;
@@ -400,9 +429,20 @@ int main (int argc, char *argv[]) try
     return FloatCmp::eq(x[0],y[0]) && FloatCmp::eq(x[1],y[1]);
   };
 
-  bool periodic = parameterSet.get<bool>("periodic");
+  // Extract all boundary vertices
+  std::vector<std::pair<std::size_t, FieldVector<double,dim> > > boundaryVertices;
+  BoundaryPatch<GridView> domainBoundary(gridView, true);
+
+  for (const auto& v1 : vertices(gridView))
+    if (domainBoundary.containsVertex(gridView.indexSet().index(v1)))
+      boundaryVertices.push_back(std::make_pair(gridView.indexSet().index(v1), v1.geometry().corner(0)));
+
+  std::cout << "Grid has " << boundaryVertices.size() << " boundary vertices" << std::endl;
+
+  bool periodic = parameterSet.get<bool>("periodic", false);
   if (periodic)
   {
+#if 0
     for (const auto& v1 : vertices(gridView))
       for (const auto& v2 : vertices(gridView))
         if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
@@ -410,6 +450,15 @@ int main (int argc, char *argv[]) try
           //std::cout << "unifying " << v1.geometry().corner(0) << " with " << v2.geometry().corner(0) << std::endl;
           feBasis.preBasis().unifyIndexPair({gridView.indexSet().index(v1)}, {gridView.indexSet().index(v2)});
         }
+#else
+    for (const auto& v1 : boundaryVertices)
+      for (const auto& v2 : boundaryVertices)
+        if (equivalent(v1.second, v2.second))
+        {
+          //std::cout << "unifying " << v1.geometry().corner(0) << " with " << v2.geometry().corner(0) << std::endl;
+          feBasis.preBasis().unifyIndexPair({v1.first}, {v2.first});
+        }
+#endif
   }
 
   feBasis.preBasis().initializeIndices();
@@ -646,11 +695,20 @@ int main (int argc, char *argv[]) try
 //                                                          FEBasis::LocalView::Tree::FiniteElement,
 //                                                          adouble> >(materialParameters);
 //
-//     if (parameterSet.get<std::string>("energy") == "expacosh")
-//       elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
-//                                                          FEBasis::LocalView::Tree::FiniteElement,
-//                                                          adouble> >(materialParameters);
-//
+#if 0
+    if (parameterSet.get<std::string>("energy") == "expacosh")
+      elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
+                                                         FEBasis::LocalView::Tree::FiniteElement,
+                                                         adouble> >(materialParameters);
+#else
+    if (parameterSet.get<std::string>("energy") == "expacosh")
+    {
+      auto density = std::make_shared<ExpACosHDensity<dim,adouble> >(F0,materialParameters);
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+                                                                       FEBasis::LocalView::Tree::FiniteElement,
+                                                                       adouble> >(density);
+    }
+#endif
 //     if (parameterSet.get<std::string>("energy") == "fung")
 //       elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView,
 //                                                   FEBasis::LocalView::Tree::FiniteElement,