diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 92e45badc309eeecd4e8d128103f6fbb6913dede..bdaaf12c7aa197aacc55fc96d0806e7f3f7ab5bd 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -120,6 +120,7 @@ setup(const typename BasisType::GridView::Grid& grid,
auto mmgStep = std::make_shared< MonotoneMGStep<MatrixType, CorrectionType> >();
+ mmgStep->setVerbosity(NumProc::REDUCED);
mmgStep->setMGType(mu, nu1, nu2);
mmgStep->ignoreNodes_ = &dirichletNodes;
mmgStep->setBaseSolver(baseSolver);
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
index ca65bde1b302f6d6942bdcaadbc7e286494329ae..838f2a18f12f4021314a79a140002feff3062563 100644
--- a/dune/elasticity/common/trustregionsolver.hh
+++ b/dune/elasticity/common/trustregionsolver.hh
@@ -113,6 +113,8 @@ public:
SolutionType getSol() const {return x_;}
+ std::string iterateNamePrefix_;
+
protected:
/** \brief The grid */
diff --git a/dune/elasticity/materials/coshenergy.hh b/dune/elasticity/materials/coshenergy.hh
index cbe6943082580d5e1c519d039a04c8e89e3346e5..f3d1070b00bccc74660531888830bb82e7c58ed1 100644
--- a/dune/elasticity/materials/coshenergy.hh
+++ b/dune/elasticity/materials/coshenergy.hh
@@ -11,7 +11,7 @@
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class CosHEnergy
+class CosHEnergy final
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
diff --git a/dune/elasticity/materials/expacoshenergy.hh b/dune/elasticity/materials/expacoshenergy.hh
index d6aa6ace12e2315f4dc2c46c8ee046f11ee30d05..b4baf69ab3f2540388ea60b2db5d1be05ddfdc58 100644
--- a/dune/elasticity/materials/expacoshenergy.hh
+++ b/dune/elasticity/materials/expacoshenergy.hh
@@ -11,7 +11,7 @@
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class ExpACosHEnergy
+class ExpACosHEnergy final
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
diff --git a/dune/elasticity/materials/neffenergy.hh b/dune/elasticity/materials/neffenergy.hh
index b4d8960c9950949ca7266b2bd73d7f0acaff8cc5..f1de30e75a35d28019feaa8addec3e2e17a0597e 100644
--- a/dune/elasticity/materials/neffenergy.hh
+++ b/dune/elasticity/materials/neffenergy.hh
@@ -11,7 +11,7 @@
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffEnergy
+class NeffEnergy final
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
diff --git a/dune/elasticity/materials/nefflogenergy.hh b/dune/elasticity/materials/nefflogenergy.hh
index 8b98a906b7bd59adb0b392e528c12cc60d8930ac..2e9acc0633e15d8b01aea37aac8e5b088f638197 100644
--- a/dune/elasticity/materials/nefflogenergy.hh
+++ b/dune/elasticity/materials/nefflogenergy.hh
@@ -11,7 +11,7 @@
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffLogEnergy
+class NeffLogEnergy final
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
diff --git a/dune/elasticity/materials/nefflogenergy2.hh b/dune/elasticity/materials/nefflogenergy2.hh
index cff60e4b1f76c5e9735401394c728072dbebc01f..e7449b8ec2073d07e28fbe46bd82ecc73b1b7a75 100644
--- a/dune/elasticity/materials/nefflogenergy2.hh
+++ b/dune/elasticity/materials/nefflogenergy2.hh
@@ -11,7 +11,7 @@
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffLogEnergy2
+class NeffLogEnergy2 final
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
diff --git a/dune/elasticity/materials/neffmagicenergy.hh b/dune/elasticity/materials/neffmagicenergy.hh
deleted file mode 100644
index ec6f4e0a2af726e5c8440cdffa0c027bf5c96807..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/neffmagicenergy.hh
+++ /dev/null
@@ -1,96 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune::Elasticity {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffMagicEnergy
- : public Elasticity::LocalEnergy<GridView,
- LocalFiniteElement,
- std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
- // grid types
- typedef typename GridView::Grid::ctype DT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- enum {gridDim=GridView::dimension};
- enum {dim=GridView::dimension};
-
-public:
-
- /** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
- */
- NeffMagicEnergy(const Dune::ParameterTree& parameters)
- {}
-
- /** \brief Assemble the energy for a single element */
- field_type energy (const Entity& e,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffMagicEnergy<GridView, LocalFiniteElement, field_type>::
-energy(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
-{
- assert(element.type() == localFiniteElement.type());
-
- field_type energy = 0;
-
- // store gradients of shape functions and base functions
- std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
- std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
-
- int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
- : localFiniteElement.localBasis().order() * gridDim;
-
- const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
- for (const auto& qp : quad)
- {
- const DT integrationElement = element.geometry().integrationElement(qp.position());
-
- const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
-
- // get gradients of shape functions
- localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
-
- // compute gradients of base functions
- for (size_t i=0; i<gradients.size(); ++i)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
- Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
- for (size_t i=0; i<gradients.size(); i++)
- for (int j=0; j<gridDim; j++)
- derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
- auto det = derivative.determinant();
- auto normSquared = derivative.frobenius_norm2();
-
- auto K = normSquared / (2*det);
-
- using std::acosh;
-
- auto energyDensity = K + std::sqrt(K*K-1) - acosh(K) + log(det);
-
- energy += qp.weight() * integrationElement * energyDensity;
- }
-
- return energy;
-}
-
-} // namespace Dune::Elasticity
-
-#endif //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
diff --git a/dune/elasticity/materials/neffreducedenergy.hh b/dune/elasticity/materials/neffreducedenergy.hh
index 0241ab9f84cacfcd71abefd4ca404b5d9685e371..85cfd0afc4460b0ce3196127b2b397ab39f05fe1 100644
--- a/dune/elasticity/materials/neffreducedenergy.hh
+++ b/dune/elasticity/materials/neffreducedenergy.hh
@@ -11,7 +11,7 @@
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffReducedEnergy
+class NeffReducedEnergy final
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
diff --git a/grids/circle2ndorder.msh b/grids/circle2ndorder.msh
index 74c3a3a79287816a5692684b0473b526a358a74f..5379fa4edc6ad6e8482b783a861f45984b36dad3 100644
--- a/grids/circle2ndorder.msh
+++ b/grids/circle2ndorder.msh
@@ -6,35 +6,35 @@ $Nodes
1 -1 0 0
2 1 0 0
3 0 0 0
-4 -0.5000000000004422 -0.8660254037841834 0
-5 0.5000000000016158 -0.8660254037835057 0
-6 -0.8660254037845663 -0.4999999999997789 0
-7 6.776189039135044e-13 -1 0
-8 0.8660254037849052 -0.4999999999991919 0
-9 0.5000000000004422 0.8660254037841834 0
-10 -0.5000000000016158 0.8660254037835057 0
-11 0.8660254037845663 0.4999999999997789 0
-12 -6.776189039135044e-13 1 0
-13 -0.8660254037849052 0.4999999999991919 0
-14 -0.02222222222229404 0.03849001794593358 0
-15 -0.261111111111955 0.4522577108647196 0
-16 0.2388888888890741 0.4522577108650584 0
-17 0.488888888888853 0.01924500897296679 0
-18 -0.2611111111113681 -0.4137676929191249 0
-19 -0.511111111111147 0.01924500897296679 0
-20 0.2388888888896609 -0.4137676929187861 0
+4 -0.5 -0.8660254037844386 0
+5 0.5 -0.8660254037844386 0
+6 -0.8660254037844386 -0.5 0
+7 0 -1 0
+8 0.8660254037844386 -0.5 0
+9 0.5 0.8660254037844386 0
+10 -0.5 0.8660254037844386 0
+11 0.8660254037844386 0.5 0
+12 0 1 0
+13 -0.8660254037844386 0.5 0
+14 0 0 0
+15 -0.25 0.4330127018922193 0
+16 0.25 0.4330127018922193 0
+17 0.5 0 0
+18 -0.25 -0.4330127018922193 0
+19 -0.5 0 0
+20 0.25 -0.4330127018922193 0
$EndNodes
$Elements
15
1 15 2 0 1 1
2 15 2 0 2 2
3 15 2 0 3 3
-4 8 3 0 1 0 1 4 6
-5 8 3 0 1 0 4 5 7
-6 8 3 0 1 0 5 2 8
-7 8 3 0 2 0 2 9 11
-8 8 3 0 2 0 9 10 12
-9 8 3 0 2 0 10 1 13
+4 8 3 0 1 0 1 4 6
+5 8 3 0 1 0 4 5 7
+6 8 3 0 1 0 5 2 8
+7 8 3 0 2 0 2 9 11
+8 8 3 0 2 0 9 10 12
+9 8 3 0 2 0 10 1 13
10 9 3 0 200 0 9 10 14 12 15 16
11 9 3 0 200 0 14 2 9 17 11 16
12 9 3 0 200 0 4 14 1 18 19 6
diff --git a/problems/quasiconvexity-circle-expacosh.parset b/problems/quasiconvexity-circle-expacosh.parset
index 3d680683759ae16cc369cfe5398e71424c716134..26f7a26e885667ee40a85f3d9bf3a7d40156ffce 100644
--- a/problems/quasiconvexity-circle-expacosh.parset
+++ b/problems/quasiconvexity-circle-expacosh.parset
@@ -6,8 +6,13 @@ structuredGrid = false
path = /home/sander/dune/dune-elasticity/grids/
gridFile = circle2ndorder.msh
+#structuredGrid = true
+#lower = 0 0
+#upper = 1 1
+#elements = 2 2
+
# Number of grid levels
-numLevels = 7
+numLevels = 6
#############################################
# Solver parameters
@@ -72,7 +77,9 @@ dirichletVerticesPredicate = "True"
dirichletValues = affine
# Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
-x0 = 4.38
+#x0 = 4.38
+x0 = 12.0186
+
# Setting initialIterateFile overrides the x0 parameter
-initialIterateFile = experiment-6.data
-oldX0 = 10
+#initialIterateFile = experiment-6.data
+#oldX0 = 10
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fffab1adba3802f898568f791793d7ca1ca1e229..e8fe65fd23a9601153418081d25eb5eda7ae885b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -10,7 +10,9 @@ if(ADOLC_FOUND AND IPOPT_FOUND AND PYTHONLIBS_FOUND AND dune-uggrid_FOUND)
add_dune_pythonlibs_flags(${_program})
add_dune_ug_flags(${_program})
add_dune_mpi_flags(${_program})
-# target_compile_options(${_program} PRIVATE "-fpermissive")
+ add_dune_superlu_flags(${_program})
+ target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
+# target_compile_options(${_program} PRIVATE "-g")
endforeach()
if (AMIRAMESH_FOUND)
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 411f9050013cf614d1b24c51f910be8eb24db30e..feff529f6a6aa342b1282da32de8ab77a2aeea1c 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -73,6 +73,8 @@ int main (int argc, char *argv[]) try
// parse data file
ParameterTree parameterSet;
+ if (argc < 2)
+ DUNE_THROW(Exception, "Usage: ./finite-strain-elasticity <parameter file>");
ParameterTreeParser::readINITree(argv[1], parameterSet);
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 91459e70a18b196353b18419bf8a7b376a866e24..5491d419c587c5024f886ba01c2c73c8eea20dd6 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -1,5 +1,7 @@
#include <config.h>
+#undef HAVE_DUNE_PARMG
+
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
#include <adolc/adouble.h>
@@ -35,12 +37,16 @@
#include <dune/elasticity/assemblers/localadolcstiffness.hh>
#include <dune/elasticity/assemblers/feassembler.hh>
#include <dune/elasticity/materials/neffenergy.hh>
+#include <dune/elasticity/materials/burkholderenergy.hh>
#include <dune/elasticity/materials/coshenergy.hh>
+#include <dune/elasticity/materials/dacorognaenergy.hh>
#include <dune/elasticity/materials/expacoshenergy.hh>
+#include <dune/elasticity/materials/fungenergy.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
#include <dune/elasticity/materials/nefflogenergy.hh>
#include <dune/elasticity/materials/nefflogenergy2.hh>
-#include <dune/elasticity/materials/neffmagicenergy.hh>
#include <dune/elasticity/materials/neffreducedenergy.hh>
+#include <dune/elasticity/materials/neffw2energy.hh>
// grid dimension
const int dim = 2;
@@ -49,6 +55,103 @@ const int order = 1;
using namespace Dune;
+template<int dim, class field_type = double>
+struct KPowerDensity final
+: public Elasticity::LocalDensity<dim,field_type>
+{
+ KPowerDensity(const Dune::ParameterTree& parameters)
+ {
+ p_ = parameters.template get<double>("p");
+ }
+
+#warning The type of 'x' should by ctype, not field_type
+ field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+ {
+ auto detF = F.determinant();
+ auto normSquared = F.frobenius_norm2();
+
+ auto K = normSquared / (2*detF);
+
+ return std::pow(K, p_);
+ }
+
+ double p_;
+};
+
+
+template<int dim, class field_type = double>
+struct MagicDensity final
+: public Elasticity::LocalDensity<dim,field_type>
+{
+ MagicDensity(const Dune::ParameterTree& parameters)
+ {
+ c_ = 1.0;
+ if (parameters.hasKey("c"))
+ c_ = parameters.template get<double>("c");
+ }
+
+ /** \brief Implement space-dependent coefficients */
+#warning The parameter should actually use ctype, not field_type
+ double c(const Dune::FieldVector<field_type,dim>& x) const
+ {
+ return c_;
+#if 0 // Three circles
+ bool inCircle0 = (x-FieldVector<double,dim>{-0.5,0}).two_norm() <0.2;
+ bool inCircle1 = (x-FieldVector<double,dim>{0.3535,0.3535}).two_norm() <0.2;
+ bool inCircle2 = (x-FieldVector<double,dim>{0.3535,-0.3535}).two_norm() <0.2;
+
+ return (inCircle0 or inCircle1 or inCircle2) ? c_ : 1;
+#endif
+// return (x.two_norm() < 0.9) ? 1 : c_;
+ }
+
+ field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+ {
+ auto det = F.determinant();
+ auto normSquared = F.frobenius_norm2();
+
+ auto K = normSquared / (2*det);
+
+ using std::acosh;
+
+ return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
+ }
+
+ double c_;
+};
+
+
+template<int dim, class field_type = double>
+struct VossDensity final
+: public Elasticity::LocalDensity<dim,field_type>
+{
+ VossDensity(const Dune::ParameterTree& parameters)
+ {
+ c_ = parameters.template get<double>("c");
+ }
+
+ field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+ {
+ auto det = F.determinant();
+ auto normSquared = F.frobenius_norm2();
+
+ auto KBold = normSquared / (2*det);
+
+ using std::acosh;
+
+ auto K = KBold + std::sqrt(KBold*KBold - 1);
+
+ auto h = (1.0/(2*K)) * (1 + K*K + (1+K)*std::sqrt(1 + 6*K + K*K)) - std::log(K)
+ + std::log(1 + 4*K + K*K + (1+K)*std::sqrt(1 + 6*K + K*K));
+
+ return h + c_ * std::log( det );
+ }
+
+ double c_;
+};
+
+
+
int main (int argc, char *argv[]) try
{
// initialize MPI, finalize is done automatically on exit
@@ -70,6 +173,8 @@ int main (int argc, char *argv[]) try
// parse data file
ParameterTree parameterSet;
+ if (argc < 2)
+ DUNE_THROW(Exception, "Usage: ./quasiconvexity-text <parameter file>");
ParameterTreeParser::readINITree(argv[1], parameterSet);
@@ -142,6 +247,7 @@ int main (int argc, char *argv[]) try
// Make Python function that computes which vertices are on the Dirichlet boundary,
// based on the vertex positions.
+
std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
@@ -163,6 +269,20 @@ int main (int argc, char *argv[]) try
for (int j=0; j<dim; j++)
dirichletDofs[i][j] = true;
+#if 0
+ // HACK!
+#warning Teichmueller hack!
+ for (auto&& v : vertices(gridView))
+ {
+ if (v.geometry().corner(0).two_norm() < 1e-6)
+ {
+ auto idx = indexSet.index(v);
+ for (int j=0; j<dim; j++)
+ dirichletDofs[idx][j] = true;
+ }
+ }
+#endif
+
// //////////////////////////
// Initial iterate
// //////////////////////////
@@ -177,25 +297,81 @@ int main (int argc, char *argv[]) try
blockedInterleaved()
));
- auto x0 = parameterSet.get<double>("x0");
+ if (parameterSet.hasKey("x0"))
+ {
+ auto x0 = parameterSet.get<double>("x0");
+
+ auto initialDeformation = [&x0](FieldVector<double,dim> x) -> FieldVector<double,dim>
+ {
+ return {x[0]*sqrt(x0), x[1] / sqrt(x0)};
+ };
+
+ Dune::Functions::interpolate(powerBasis, x, initialDeformation);
+ }
- auto initialDeformation = [&x0](FieldVector<double,dim> x) -> FieldVector<double,dim>
+ if (parameterSet.hasKey("initialIteratePython"))
{
- return {x[0]*sqrt(x0), x[1] / sqrt(x0)};
- };
+ // The python class that implements the Dirichlet values
+ Python::Module initialIterateModule = Python::import(parameterSet.get<std::string>("initialIteratePython"));
- Dune::Functions::interpolate(powerBasis, x, initialDeformation);
+ auto f = Python::Callable(initialIterateModule.get("initialIterate"));
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > initialIteratePython(f);
+
+ Dune::Functions::interpolate(powerBasis, x, initialIteratePython);
+ }
if (parameterSet.hasKey("initialIterateFile"))
{
- std::ifstream inFile(parameterSet.get<std::string>("initialIterateFile"), std::ios_base::binary);
- if (not inFile)
- DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("initialIterateFile") << " could not be opened.");
- MatrixVector::Generic::readBinary(inFile, x);
- inFile.peek(); // try to advance beyond the end of the file
- if (not inFile.eof())
- DUNE_THROW(IOError, "File '" << parameterSet.get<std::string>("initialIterateFile") << "' does not have the correct size!");
- inFile.close();
+ std::string filename = parameterSet.get<std::string>("initialIterateFile");
+
+ // Guess the file format by looking at the file name suffix
+ auto dotPos = filename.rfind('.');
+ if (dotPos == std::string::npos)
+ DUNE_THROW(IOError, "Could not determine grid input file format");
+ std::string suffix = filename.substr(dotPos, filename.length()-dotPos);
+
+ if (suffix == ".data")
+ {
+ std::ifstream inFile(filename, std::ios_base::binary);
+ if (not inFile)
+ DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("initialIterateFile") << " could not be opened.");
+ MatrixVector::Generic::readBinary(inFile, x);
+ inFile.peek(); // try to advance beyond the end of the file
+ if (not inFile.eof())
+ DUNE_THROW(IOError, "File '" << parameterSet.get<std::string>("initialIterateFile") << "' does not have the correct size!");
+ inFile.close();
+ }
+
+ // Read CSV files as provided by Sid Kumar
+ if (suffix == ".csv")
+ {
+ std::ifstream inFile(filename, std::ios_base::binary);
+ if (not inFile)
+ DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("initialIterateFile") << " could not be opened.");
+
+ int vertex = 0;
+ for (std::string line; std::getline(inFile, line);)
+ {
+ // The first line starts with 'x' -- skip it
+ if (line[0]=='x')
+ continue;
+
+ std::vector<std::string> result;
+ std::stringstream s_stream(line); //create string stream from the string
+ while(s_stream.good()) {
+ std::string substr;
+ getline(s_stream, substr, ','); //get first string delimited by comma
+ result.push_back(substr);
+ }
+
+ x[vertex][0] = atof(result[2].c_str());
+ x[vertex][1] = atof(result[3].c_str());
+
+ vertex++;
+ }
+
+ inFile.close();
+ }
}
////////////////////////////////////////////////////////
@@ -238,10 +414,18 @@ int main (int argc, char *argv[]) try
// Assembler using ADOL-C
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
- std::shared_ptr<LocalFEStiffness<GridView,
+ std::shared_ptr<Elasticity::LocalEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
std::vector<Dune::FieldVector<adouble, dim> > > > elasticEnergy;
+ if (parameterSet.get<std::string>("energy") == "kpower")
+ {
+ auto density = std::make_shared<KPowerDensity<dim,adouble> >(materialParameters);
+ elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(density);
+ }
+
if (parameterSet.get<std::string>("energy") == "neff")
elasticEnergy = std::make_shared<NeffEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
@@ -272,10 +456,41 @@ int main (int argc, char *argv[]) try
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "fung")
+ elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(materialParameters);
+
if (parameterSet.get<std::string>("energy") == "magic")
- elasticEnergy = std::make_shared<Elasticity::NeffMagicEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(materialParameters);
+ {
+ auto density = std::make_shared<MagicDensity<dim,adouble> >(materialParameters);
+ elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(density);
+ }
+
+ if (parameterSet.get<std::string>("energy") == "burkholder")
+ elasticEnergy = std::make_shared<Elasticity::BurkholderEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(materialParameters);
+
+ if (parameterSet.get<std::string>("energy") == "w2")
+ elasticEnergy = std::make_shared<Elasticity::NeffW2Energy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(materialParameters);
+
+ if (parameterSet.get<std::string>("energy") == "dacorogna")
+ elasticEnergy = std::make_shared<Elasticity::DacorognaEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(materialParameters);
+
+ if (parameterSet.get<std::string>("energy") == "voss")
+ {
+ auto density = std::make_shared<VossDensity<dim,adouble> >(materialParameters);
+ elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ adouble> >(density);
+ }
if(!elasticEnergy)
DUNE_THROW(Exception, "Error: Selected energy not available!");
@@ -306,7 +521,8 @@ int main (int argc, char *argv[]) try
mgTolerance,
mu, nu1, nu2,
baseIterations,
- baseTolerance
+ baseTolerance,
+ 1.0
);
////////////////////////////////////////////////////////
@@ -328,9 +544,9 @@ int main (int argc, char *argv[]) try
#else
// lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
// PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > dirichletValues(Python::evaluate(lambda));
- Python::Module module = Python::import(parameterSet.get<std::string>("dirichletValues"));
- auto dirichletValues = Python::makeFunction<FieldVector<double,dim>(const FieldVector<double,dim>&)>(module.get("f"));
- Dune::Functions::interpolate(powerBasis, x, dirichletValues, dirichletDofs);
+// Python::Module module = Python::import(parameterSet.get<std::string>("dirichletValues"));
+// auto dirichletValues = Python::makeFunction<FieldVector<double,dim>(const FieldVector<double,dim>&)>(module.get("f"));
+// Dune::Functions::interpolate(powerBasis, x, dirichletValues, dirichletDofs);
#endif
//::Functions::interpolate(fufemFEBasis, x, dirichletValues, dirichletDofs);
@@ -339,6 +555,7 @@ int main (int argc, char *argv[]) try
// Solve!
// /////////////////////////////////////////////////////
+ solver.iterateNamePrefix_ = "stage1-";
solver.setInitialIterate(x);
solver.solve();
@@ -348,6 +565,7 @@ int main (int argc, char *argv[]) try
// Output result
/////////////////////////////////
+#if 0
std::cout << "Total energy: " << assembler.computeEnergy(x) << std::endl;
elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
@@ -363,7 +581,7 @@ int main (int argc, char *argv[]) try
localADOLCStiffness.localEnergy_ = elasticEnergy.get();
std::cout << "Energy part 2: " << assembler.computeEnergy(x) << std::endl;
-
+#endif
// Compute the displacement
auto displacement = x;
displacement -= identity;
@@ -445,11 +663,14 @@ int main (int argc, char *argv[]) try
}
#endif
+ std::cout << "p: " << parameterSet.sub("materialParameters")["p"] << ", "
+ << "K range: " << (*std::min_element(K.begin(), K.end()))
+ << ", " << (*std::max_element(K.begin(), K.end())) << std::endl;
+
auto deformationDeterminantFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, deformationDeterminant);
auto KFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, K);
- // We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+ VTKWriter<GridView> vtkWriter(gridView);
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addCellData(deformationDeterminantFunction,
VTK::FieldInfo("determinant", VTK::FieldInfo::Type::scalar, 1));
@@ -463,22 +684,24 @@ int main (int argc, char *argv[]) try
MatrixVector::Generic::writeBinary(outFile, x);
outFile.close();
+ exit(0);
+#if 0
//////////////////////////////////////////////////////////////////////////////////
// Use the solution as initial iterate for the 'real' energy
//////////////////////////////////////////////////////////////////////////////////
ParameterTree neffMaterialParameters;
- neffMaterialParameters["scaling"] = "1.15";
+ neffMaterialParameters["c"] = "1";
// Careful: The following code computes the initial iterate given by the Python file,
// and *assumes* that that file describes the homogeneous deformation!
SolutionType homogeneous(feBasis.size());
Dune::Functions::interpolate(powerBasis, homogeneous, initialDeformation);
- std::shared_ptr<LocalFEStiffness<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- std::vector<Dune::FieldVector<adouble, dim> > > >
- energy = std::make_shared<NeffEnergy<GridView,
+ std::shared_ptr<Elasticity::LocalEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ std::vector<Dune::FieldVector<adouble, dim> > > >
+ energy = std::make_shared<Elasticity::NeffW2Energy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(neffMaterialParameters);
@@ -510,6 +733,7 @@ int main (int argc, char *argv[]) try
baseTolerance
);
+ solver.iterateNamePrefix_ = "stage2-";
solver.setInitialIterate(x);
solver.solve();
@@ -522,7 +746,7 @@ int main (int argc, char *argv[]) try
vtkWriter.clear();
vtkWriter.addVertexData(displacementFunction2, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(resultPath + "final-result");
-
+#endif
} catch (Exception& e) {
std::cout << e.what() << std::endl;
}