From 31667d23f4835caaa96c67af3a9bce7a2b4a901c Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 8 Jun 2020 09:56:34 +0200
Subject: [PATCH] WIP
---
dune/elasticity/common/trustregionsolver.cc | 2 +-
.../quasiconvexity-circle-expacosh.parset | 2 +-
src/CMakeLists.txt | 4 +-
src/quasiconvexity-test.cc | 72 +++++++++++++++++--
4 files changed, 69 insertions(+), 11 deletions(-)
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 1209d0e..6200c1e 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/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;
diff --git a/problems/quasiconvexity-circle-expacosh.parset b/problems/quasiconvexity-circle-expacosh.parset
index 26f7a26..247656d 100644
--- a/problems/quasiconvexity-circle-expacosh.parset
+++ b/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
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e8fe65f..0d66716 100644
--- a/src/CMakeLists.txt
+++ b/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)
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 0c1d7a6..07aa142 100644
--- a/src/quasiconvexity-test.cc
+++ b/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,
--
GitLab