From f466f865308a0b0784caadd41d9c338dfc58019b Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 30 Aug 2019 12:52:19 +0200
Subject: [PATCH] WIP
---
dune/elasticity/assemblers/feassembler.hh | 2 +-
dune/elasticity/common/trustregionsolver.cc | 70 ++++++++++++
dune/elasticity/common/trustregionsolver.hh | 2 +
dune/elasticity/materials/burkholderenergy.hh | 99 +++++++++++++++++
.../materials/burkholderstarenergy.hh | 102 ++++++++++++++++++
dune/elasticity/materials/cosh2energy.hh | 96 +++++++++++++++++
dune/elasticity/materials/dacorognaenergy.hh | 97 +++++++++++++++++
dune/elasticity/materials/neffmagicenergy.hh | 27 ++++-
dune/elasticity/materials/neffw2energy.hh | 102 ++++++++++++++++++
.../quasiconvexity-circle-burkholder.parset | 81 ++++++++++++++
problems/quasiconvexity-circle-cosh2.parset | 75 +++++++++++++
.../quasiconvexity-circle-dacorogna.parset | 84 +++++++++++++++
.../quasiconvexity-circle-expacosh.parset | 15 ++-
problems/quasiconvexity-circle-magic.parset | 84 +++++++++++++++
problems/quasiconvexity-circle-w2.parset | 85 +++++++++++++++
src/CMakeLists.txt | 2 +
src/finite-strain-elasticity.cc | 2 +
src/quasiconvexity-test.cc | 90 ++++++++++++++--
18 files changed, 1100 insertions(+), 15 deletions(-)
create mode 100644 dune/elasticity/materials/burkholderenergy.hh
create mode 100644 dune/elasticity/materials/burkholderstarenergy.hh
create mode 100644 dune/elasticity/materials/cosh2energy.hh
create mode 100644 dune/elasticity/materials/dacorognaenergy.hh
create mode 100644 dune/elasticity/materials/neffw2energy.hh
create mode 100644 problems/quasiconvexity-circle-burkholder.parset
create mode 100644 problems/quasiconvexity-circle-cosh2.parset
create mode 100644 problems/quasiconvexity-circle-dacorogna.parset
create mode 100644 problems/quasiconvexity-circle-magic.parset
create mode 100644 problems/quasiconvexity-circle-w2.parset
diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh
index 0486478..0ce8dbb 100644
--- a/dune/elasticity/assemblers/feassembler.hh
+++ b/dune/elasticity/assemblers/feassembler.hh
@@ -70,7 +70,7 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
nb.resize(n, n);
- for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
+ for (const auto& element : elements(basis_.getGridView()/*, Dune::Partitions::interior*/))
{
const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(element);
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 5d16bf1..7b39adb 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -98,6 +98,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);
@@ -288,6 +289,54 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
i==0 // assemble occupation pattern only for the first call
);
+
+
+ ////////////////////////////////////////////////////////////////////
+ // Try to determine the largest negative eigen value
+ ////////////////////////////////////////////////////////////////////
+#if 0
+ double upperBoundForEV = 1e4;
+ auto A = *hessianMatrix_;
+ for (std::size_t j=0; j<A.N(); j++)
+ {
+ A[i][i][0][0] -= upperBoundForEV;
+ A[i][i][1][1] -= upperBoundForEV;
+ }
+
+ // If the upper bound has been large enough, A is negative definite
+ CorrectionType eigenVector(rhs.size());
+ std::fill(eigenVector.begin(), eigenVector.end(), 1.0);
+
+ for (int j=0; j<1e6; j++)
+ {
+ if ( (j%1000)==0)
+ std::cout << "power iteration: " <<j << std::endl;
+ CorrectionType nextEV(eigenVector.size());
+
+ A.mv(eigenVector,nextEV);
+ nextEV /= nextEV.two_norm();
+
+ eigenVector = nextEV;
+
+ // Compute the currect approximation to the eigen value
+ CorrectionType tmp(eigenVector.size());
+ A.mv(eigenVector,tmp);
+
+ double eigenValue = (eigenVector * tmp) / eigenVector.two_norm2();
+
+ //std::cout << "eigen value: " << (-eigenValue + upperBoundForEV) << std::endl;
+// std::cout << "eigen value: " << eigenValue << std::endl;
+ }
+
+ CorrectionType tmp(eigenVector.size());
+ hessianMatrix_->mv(eigenVector,tmp);
+
+ std::cout << "real eigen value: " << ((eigenVector * tmp) / eigenVector.two_norm2()) << std::endl;
+#endif
+
+ ///////////////////////////////////////////////////////////////////
+ // Continue with the actual computation
+ ///////////////////////////////////////////////////////////////////
rhs *= -1; // The right hand side is the _negative_ gradient
// Compute gradient norm to monitor convergence
@@ -482,6 +531,27 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
for (size_t j=0; j<newIterate.size(); j++)
newIterate[j] += corr[j];
+ // Write newIterate to a file, so we can see whether it is worth waiting for the end of
+ // the simulation
+ using namespace Dune::Functions::BasisFactory;
+ auto powerBasis = makeBasis(
+ grid_->leafGridView(),
+ power<2>(
+ lagrange<1>(),
+ blockedInterleaved()
+ ));
+
+ SolutionType identity;
+ Dune::Functions::interpolate(powerBasis, identity, [](Dune::FieldVector<double,2> x){ return x; });
+
+
+ auto displacement = newIterate;
+ displacement -= identity;
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,2>>(basis, displacement);
+ Dune::VTKWriter<typename BasisType::GridView> vtkWriter(grid_->leafGridView());
+ vtkWriter.addVertexData(displacementFunction, Dune::VTK::FieldInfo("displacement", Dune::VTK::FieldInfo::Type::vector, 2));
+ vtkWriter.write(iterateNamePrefix_ + "tr-iterate-" + std::to_string(i));
+
double energy = assembler_->computeEnergy(newIterate);
energy = grid_->comm().sum(energy);
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
index 1b87165..2039d9f 100644
--- a/dune/elasticity/common/trustregionsolver.hh
+++ b/dune/elasticity/common/trustregionsolver.hh
@@ -99,6 +99,8 @@ public:
SolutionType getSol() const {return x_;}
+ std::string iterateNamePrefix_;
+
protected:
/** \brief The grid */
diff --git a/dune/elasticity/materials/burkholderenergy.hh b/dune/elasticity/materials/burkholderenergy.hh
new file mode 100644
index 0000000..1f63ffc
--- /dev/null
+++ b/dune/elasticity/materials/burkholderenergy.hh
@@ -0,0 +1,99 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_BURKHOLDERENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_BURKHOLDERENERGY_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 BurkholderEnergy
+ : public LocalFEStiffness<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
+ */
+ BurkholderEnergy(const Dune::ParameterTree& parameters)
+ {
+ p_ = parameters.template get<double>("p");
+ }
+
+ /** \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;
+ double p_;
+
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+BurkholderEnergy<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 lambdaMaxSquared = 0.5 * (normSquared + std::sqrt(normSquared*normSquared - 4*det*det));
+ auto lambdaMax = std::sqrt(lambdaMaxSquared);
+
+ auto energyDensity = 0.5 * p_ * det * std::pow(lambdaMax,p_-2) + 0.5*(2-p_) * std::pow(lambdaMax,p_);
+
+ energy -= qp.weight() * integrationElement * energyDensity;
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_BURKHOLDERENERGY_HH
diff --git a/dune/elasticity/materials/burkholderstarenergy.hh b/dune/elasticity/materials/burkholderstarenergy.hh
new file mode 100644
index 0000000..efa2cc9
--- /dev/null
+++ b/dune/elasticity/materials/burkholderstarenergy.hh
@@ -0,0 +1,102 @@
+#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 LocalFEStiffness<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)
+ {
+ c_ = 1.0;
+ if (parameters.hasKey("c"))
+ c_ = parameters.template get<double>("c");
+ }
+
+ /** \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;
+ double c_;
+
+};
+
+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) + c_*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/cosh2energy.hh b/dune/elasticity/materials/cosh2energy.hh
new file mode 100644
index 0000000..43e17fb
--- /dev/null
+++ b/dune/elasticity/materials/cosh2energy.hh
@@ -0,0 +1,96 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_COSH2ENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_COSH2ENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class CosH2Energy
+ : public LocalFEStiffness<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
+ */
+ CosH2Energy(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
+CosH2Energy<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::cosh;
+
+ auto energyDensity = cosh(K-2)-1;
+
+ energy += qp.weight() * integrationElement * energyDensity;
+ }
+
+ return energy;
+}
+
+} // namespace Dune
+
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_COSH2ENERGY_HH
diff --git a/dune/elasticity/materials/dacorognaenergy.hh b/dune/elasticity/materials/dacorognaenergy.hh
new file mode 100644
index 0000000..b1632e8
--- /dev/null
+++ b/dune/elasticity/materials/dacorognaenergy.hh
@@ -0,0 +1,97 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_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 DacorognaEnergy
+ : public LocalFEStiffness<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
+ */
+ DacorognaEnergy(const Dune::ParameterTree& parameters)
+ {
+ gamma_ = parameters.template get<double>("gamma");
+ }
+
+ /** \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;
+
+ double gamma_;
+
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+DacorognaEnergy<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 energyDensity = normSquared * ( normSquared - gamma_ * det);
+
+ energy += qp.weight() * integrationElement * energyDensity;
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
diff --git a/dune/elasticity/materials/neffmagicenergy.hh b/dune/elasticity/materials/neffmagicenergy.hh
index 4c09fa2..a06450e 100644
--- a/dune/elasticity/materials/neffmagicenergy.hh
+++ b/dune/elasticity/materials/neffmagicenergy.hh
@@ -30,12 +30,31 @@ public:
* \param parameters The material parameters
*/
NeffMagicEnergy(const Dune::ParameterTree& parameters)
- {}
+ {
+ c_ = 1.0;
+ if (parameters.hasKey("c"))
+ c_ = parameters.template get<double>("c");
+ }
/** \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;
+
+
+ double c(const Dune::FieldVector<double,dim>& x) const
+ {
+#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_;
+ }
+ double c_;
+
};
template <class GridView, class LocalFiniteElement, class field_type>
@@ -83,7 +102,9 @@ energy(const Entity& element,
using std::acosh;
- auto energyDensity = K + std::sqrt(K*K-1) - acosh(K) + log(det);
+ auto worldPos = element.geometry().global(qp.position());
+
+ auto energyDensity = K + std::sqrt(K*K-1) - acosh(K) + c(worldPos)*log(det);
energy += qp.weight() * integrationElement * energyDensity;
}
@@ -93,4 +114,4 @@ energy(const Entity& element,
} // namespace Dune::Elasticity
-#endif //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
diff --git a/dune/elasticity/materials/neffw2energy.hh b/dune/elasticity/materials/neffw2energy.hh
new file mode 100644
index 0000000..5b5f850
--- /dev/null
+++ b/dune/elasticity/materials/neffw2energy.hh
@@ -0,0 +1,102 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_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 NeffW2Energy
+ : public LocalFEStiffness<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
+ */
+ NeffW2Energy(const Dune::ParameterTree& parameters)
+ {
+ c_ = 1.0;
+ if (parameters.hasKey("c"))
+ c_ = parameters.template get<double>("c");
+ }
+
+ /** \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;
+ double c_;
+
+};
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+NeffW2Energy<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 + acosh(K) - c_*log(det);
+
+ energy += qp.weight() * integrationElement * energyDensity;
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
diff --git a/problems/quasiconvexity-circle-burkholder.parset b/problems/quasiconvexity-circle-burkholder.parset
new file mode 100644
index 0000000..6050295
--- /dev/null
+++ b/problems/quasiconvexity-circle-burkholder.parset
@@ -0,0 +1,81 @@
+#############################################
+# Grid parameters
+#############################################
+
+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 = 8
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.01
+
+# Number of multigrid iterations per trust-region step
+numIt = 200
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+# Material parameters
+############################
+
+energy = burkholder
+
+[materialParameters]
+
+p = 3
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = neff-quasiconvexity
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "True"
+dirichletValues = affine
+
+# Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
+x0 = 3
+
+# Setting initialIterateFile overrides the x0 parameter
+#initialIterateFile = quasiconvexity-test-result-1.075-x0_1.1_8levels.data
diff --git a/problems/quasiconvexity-circle-cosh2.parset b/problems/quasiconvexity-circle-cosh2.parset
new file mode 100644
index 0000000..124a58a
--- /dev/null
+++ b/problems/quasiconvexity-circle-cosh2.parset
@@ -0,0 +1,75 @@
+#############################################
+# Grid parameters
+#############################################
+
+structuredGrid = false
+sphereGrid = true
+#path = /home/sander/dune/dune-elasticity/grids/
+#gridFile = circle2ndorder.msh
+
+# Number of grid levels
+numLevels = 9
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 10000
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.01
+
+# Number of multigrid iterations per trust-region step
+numIt = 500
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+# Material parameters
+############################
+
+energy = cosh2
+
+[materialParameters]
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = neff-quasiconvexity
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "True"
+#dirichletValues = affine
+
+# Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
+x0 = 1
+# Setting initialIterateFile overrides the x0 parameter
+#initialIterateFile = experiment-6.data
+#oldX0 = 10
diff --git a/problems/quasiconvexity-circle-dacorogna.parset b/problems/quasiconvexity-circle-dacorogna.parset
new file mode 100644
index 0000000..27abad9
--- /dev/null
+++ b/problems/quasiconvexity-circle-dacorogna.parset
@@ -0,0 +1,84 @@
+#############################################
+# Grid parameters
+#############################################
+
+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 = 6
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.01
+
+# Number of multigrid iterations per trust-region step
+numIt = 200
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+# Material parameters
+############################
+
+energy = dacorogna
+
+[materialParameters]
+
+gamma = 2.31
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = neff-quasiconvexity
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "True"
+dirichletValues = affine
+
+# Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
+x0 = 1.75
+
+# Initial iterate given as a Python function
+#initialIteratePython = one-circle
+
+# Setting initialIterateFile overrides the x0 parameter
+#initialIterateFile = quasiconvexity-test-result-1.075-x0_1.1_8levels.data
diff --git a/problems/quasiconvexity-circle-expacosh.parset b/problems/quasiconvexity-circle-expacosh.parset
index 3d68068..26f7a26 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/problems/quasiconvexity-circle-magic.parset b/problems/quasiconvexity-circle-magic.parset
new file mode 100644
index 0000000..6e45ab9
--- /dev/null
+++ b/problems/quasiconvexity-circle-magic.parset
@@ -0,0 +1,84 @@
+#############################################
+# Grid parameters
+#############################################
+
+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 = 6
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.01
+
+# Number of multigrid iterations per trust-region step
+numIt = 200
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+# Material parameters
+############################
+
+energy = magic
+
+[materialParameters]
+
+c = 2
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = neff-quasiconvexity
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "True"
+dirichletValues = affine
+
+# Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
+x0 = 2
+
+# Initial iterate given as a Python function
+#initialIteratePython = one-circle
+
+# Setting initialIterateFile overrides the x0 parameter
+#initialIterateFile = quasiconvexity-test-result-1.075-x0_1.1_8levels.data
diff --git a/problems/quasiconvexity-circle-w2.parset b/problems/quasiconvexity-circle-w2.parset
new file mode 100644
index 0000000..232b3f4
--- /dev/null
+++ b/problems/quasiconvexity-circle-w2.parset
@@ -0,0 +1,85 @@
+#############################################
+# Grid parameters
+#############################################
+
+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 = 6
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.01
+
+# Number of multigrid iterations per trust-region step
+numIt = 200
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+# Material parameters
+############################
+
+energy = w2
+
+[materialParameters]
+
+#c = 1.075
+c = 0.5
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = neff-quasiconvexity
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "True"
+dirichletValues = affine
+
+# Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
+x0 = 10
+
+# Initial iterate given as a Python function
+#initialIteratePython = one-circle
+
+# Setting initialIterateFile overrides the x0 parameter
+#initialIterateFile = quasiconvexity-test-result-1.075-x0_1.1_8levels.data
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fffab1a..9b5395f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -10,6 +10,8 @@ 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})
+ add_dune_superlu_flags(${_program})
+ target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
# target_compile_options(${_program} PRIVATE "-fpermissive")
endforeach()
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 519a98d..696d2f7 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -91,6 +91,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 91459e7..300a0ce 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>
@@ -26,6 +28,7 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
+#include <dune/fufem/functionspacebases/periodicbasis.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -35,12 +38,15 @@
#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/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;
@@ -70,6 +76,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);
@@ -129,9 +137,46 @@ int main (int argc, char *argv[]) try
#endif
// FE basis spanning the FE space that we are working in
- typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
+// typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
+// FEBasis feBasis(gridView);
+ using FEBasis = Fufem::PeriodicBasis<GridView>;
FEBasis feBasis(gridView);
+ bool periodicY = false;
+ if (periodicY)
+ {
+ if (order!=1)
+ DUNE_THROW(NotImplemented, "Periodic boundary conditions are only implemented for order 1!");
+
+ std::map<double,std::size_t> bottomMap, topMap;
+
+ for (const auto& vertex : vertices(gridView))
+ {
+ auto p = vertex.geometry().corner(0);
+ auto idx = gridView.indexSet().index(vertex);
+
+ if (p[1] < 1e-9)
+ bottomMap.insert(std::make_pair(p[0], idx));
+ if (p[1] > 1 - 1e-9)
+ topMap.insert(std::make_pair(p[0], idx));
+
+ }
+
+ if (topMap.size() != bottomMap.size())
+ DUNE_THROW(Exception, "Bottom and top vertex lists do not have the same size");
+
+ auto bottomIt = bottomMap.begin();
+ auto topIt = topMap.begin();
+
+ while (bottomIt != bottomMap.end())
+ {
+ feBasis.preBasis().unifyIndexPair({bottomIt->second}, {topIt->second});
+ ++bottomIt;
+ ++topIt;
+ }
+ feBasis.preBasis().initializeIndices();
+ }
+
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
@@ -186,6 +231,17 @@ int main (int argc, char *argv[]) try
Dune::Functions::interpolate(powerBasis, x, initialDeformation);
+ if (parameterSet.hasKey("initialIteratePython"))
+ {
+ // The python class that implements the Dirichlet values
+ Python::Module initialIterateModule = Python::import(parameterSet.get<std::string>("initialIteratePython"));
+
+ 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);
@@ -277,6 +333,21 @@ int main (int argc, char *argv[]) try
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters);
+ 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(!elasticEnergy)
DUNE_THROW(Exception, "Error: Selected energy not available!");
@@ -328,9 +399,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 +410,7 @@ int main (int argc, char *argv[]) try
// Solve!
// /////////////////////////////////////////////////////
+ solver.iterateNamePrefix_ = "stage1-";
solver.setInitialIterate(x);
solver.solve();
@@ -348,6 +420,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 +436,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;
@@ -463,12 +536,14 @@ int main (int argc, char *argv[]) try
MatrixVector::Generic::writeBinary(outFile, x);
outFile.close();
+ exit(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!
@@ -478,7 +553,7 @@ int main (int argc, char *argv[]) try
std::shared_ptr<LocalFEStiffness<GridView,
FEBasis::LocalView::Tree::FiniteElement,
std::vector<Dune::FieldVector<adouble, dim> > > >
- energy = std::make_shared<NeffEnergy<GridView,
+ energy = std::make_shared<Elasticity::NeffW2Energy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(neffMaterialParameters);
@@ -510,6 +585,7 @@ int main (int argc, char *argv[]) try
baseTolerance
);
+ solver.iterateNamePrefix_ = "stage2-";
solver.setInitialIterate(x);
solver.solve();
--
GitLab