diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index bdaaf12c7aa197aacc55fc96d0806e7f3f7ab5bd..1496cffd2a809270e07e25c71bb84555d24155c6 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -388,6 +388,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
std::cout << "----------------------------------------------------" << std::endl;
std::cout << " Trust-Region Step Number: " << i
<< ", radius: " << trustRegion.radius()
+ << std::setprecision(15)
<< ", energy: " << oldEnergy << std::endl;
std::cout << "----------------------------------------------------" << std::endl;
}
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 5491d419c587c5024f886ba01c2c73c8eea20dd6..7d2da5ffe68970e65c20e22dcce8c69fd31b6773 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -55,18 +55,19 @@ const int order = 1;
using namespace Dune;
-template<int dim, class field_type = double>
+template<int dim, class field_type, class ctype=double>
struct KPowerDensity final
-: public Elasticity::LocalDensity<dim,field_type>
+: public Elasticity::LocalDensity<dim,field_type,ctype>
{
- KPowerDensity(const Dune::ParameterTree& parameters)
+ KPowerDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+ : F0_(F0)
{
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
+ field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
+ auto F = derivative + F0_;
auto detF = F.determinant();
auto normSquared = F.frobenius_norm2();
@@ -75,15 +76,17 @@ struct KPowerDensity final
return std::pow(K, p_);
}
+ FieldMatrix<ctype,dim,dim> F0_;
double p_;
};
-template<int dim, class field_type = double>
+template<int dim, class field_type, class ctype=double>
struct MagicDensity final
-: public Elasticity::LocalDensity<dim,field_type>
+: public Elasticity::LocalDensity<dim,field_type,ctype>
{
- MagicDensity(const Dune::ParameterTree& parameters)
+ MagicDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+ : F0_(F0)
{
c_ = 1.0;
if (parameters.hasKey("c"))
@@ -91,8 +94,7 @@ struct MagicDensity final
}
/** \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
+ double c(const Dune::FieldVector<ctype,dim>& x) const
{
return c_;
#if 0 // Three circles
@@ -105,8 +107,9 @@ struct MagicDensity final
// 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
+ 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();
@@ -117,21 +120,24 @@ struct MagicDensity final
return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
}
+ FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
-template<int dim, class field_type = double>
+template<int dim, class field_type, class ctype=double>
struct VossDensity final
-: public Elasticity::LocalDensity<dim,field_type>
+: public Elasticity::LocalDensity<dim,field_type,ctype>
{
- VossDensity(const Dune::ParameterTree& parameters)
+ VossDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+ : F0_(F0)
{
c_ = parameters.template get<double>("c");
}
- field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+ 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();
@@ -147,6 +153,7 @@ struct VossDensity final
return h + c_ * std::log( det );
}
+ FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
@@ -210,7 +217,7 @@ int main (int argc, char *argv[]) try
upper = parameterSet.get<FieldVector<double,dim> >("upper");
std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
- grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+ grid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elements);
} else {
std::string path = parameterSet.get<std::string>("path");
@@ -297,17 +304,15 @@ int main (int argc, char *argv[]) try
blockedInterleaved()
));
- if (parameterSet.hasKey("x0"))
- {
- auto x0 = parameterSet.get<double>("x0");
+ // Set the underlying affine deformation
+ 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)};
- };
+ FieldMatrix<double,2,2> F0 = { {sqrt(x0), 0},
+ { 0, 1.0 / sqrt(x0)} };
- Dune::Functions::interpolate(powerBasis, x, initialDeformation);
- }
+ // The actual unknown is a perturbation of the affine transformation
+ // Start here with the zero function
+ std::fill(x.begin(), x.end(), FieldVector<double,2>(0));
if (parameterSet.hasKey("initialIteratePython"))
{
@@ -364,9 +369,23 @@ int main (int argc, char *argv[]) try
result.push_back(substr);
}
- x[vertex][0] = atof(result[2].c_str());
- x[vertex][1] = atof(result[3].c_str());
+ FieldVector<double,2> pos{atof(result[0].c_str()), atof(result[1].c_str())};
+ FieldVector<double,2> u{atof(result[2].c_str()), atof(result[3].c_str())};
+
+ //x[vertex] = u - (F0-Id)*pos;
+ FieldVector<double,2> phi = { u[0] - (F0[0][0]-1)*pos[0] - F0[0][1] *pos[1],
+ u[1] - F0[1][0] *pos[0] - (F0[1][1]-1)*pos[1] };
+
+ for (const auto& vert : vertices(gridView))
+ if ( (vert.geometry().corner(0)-pos).two_norm() < 1e-8)
+ {
+ x[gridView.indexSet().index(vert)] = phi;
+ std::cout << "pos: " << pos << ", vertex: " << gridView.indexSet().index(vert) << std::endl;
+ break;
+ }
+ //x[vertex] = phi;
+ //std::cout << "pos: " << pos << ", phi: " << x[vertex] << std::endl;
vertex++;
}
@@ -380,7 +399,10 @@ int main (int argc, char *argv[]) try
// Output initial iterate (of homotopy loop)
SolutionType identity;
- Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+ Dune::Functions::interpolate(powerBasis, identity, [&F0](FieldVector<double,dim> x)
+ {
+ return FieldVector<double,2>{ (F0[0][0]-1)*x[0] + F0[0][1]*x[1], F0[1][0]*x[0] + (F0[1][1]-1)*x[1] };
+ });
SolutionType displacement = x;
displacement -= identity;
@@ -388,8 +410,8 @@ int main (int argc, char *argv[]) try
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
auto localDisplacementFunction = localFunction(displacementFunction);
- // We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+ // Write the initial iterate
+ VTKWriter<GridView> vtkWriter(gridView);
vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(resultPath + "quasiconvexity-test_homotopy_0");
@@ -420,73 +442,73 @@ int main (int argc, char *argv[]) try
if (parameterSet.get<std::string>("energy") == "kpower")
{
- auto density = std::make_shared<KPowerDensity<dim,adouble> >(materialParameters);
+ auto density = std::make_shared<KPowerDensity<dim,adouble> >(F0, 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,
- adouble> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "nefflog")
- elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "nefflog2")
- elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "neffreduced")
- elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- adouble> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "cosh")
- elasticEnergy = std::make_shared<CosHEnergy<GridView,
- 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 (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") == "neff")
+// elasticEnergy = std::make_shared<NeffEnergy<GridView,
+// FEBasis::LocalView::Tree::FiniteElement,
+// adouble> >(materialParameters);
+//
+// if (parameterSet.get<std::string>("energy") == "nefflog")
+// elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
+// FEBasis::LocalView::Tree::FiniteElement,
+// adouble> >(materialParameters);
+//
+// if (parameterSet.get<std::string>("energy") == "nefflog2")
+// elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
+// FEBasis::LocalView::Tree::FiniteElement,
+// adouble> >(materialParameters);
+//
+// if (parameterSet.get<std::string>("energy") == "neffreduced")
+// elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
+// FEBasis::LocalView::Tree::FiniteElement,
+// adouble> >(materialParameters);
+//
+// if (parameterSet.get<std::string>("energy") == "cosh")
+// elasticEnergy = std::make_shared<CosHEnergy<GridView,
+// 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 (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")
{
- auto density = std::make_shared<MagicDensity<dim,adouble> >(materialParameters);
+ auto density = std::make_shared<MagicDensity<dim,adouble> >(F0,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") == "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);
+ auto density = std::make_shared<VossDensity<dim,adouble> >(F0,materialParameters);
elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(density);
@@ -654,7 +676,7 @@ int main (int argc, char *argv[]) try
for (int j=0; j<dim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
- //auto det = derivative.determinant();
+ derivative += F0;
//std::cout << derivative << std::endl;
deformationDeterminant[localP0View.index(0)] = derivative.determinant();
@@ -741,8 +763,7 @@ int main (int argc, char *argv[]) try
displacement -= identity;
auto displacementFunction2 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
- // 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.clear();
vtkWriter.addVertexData(displacementFunction2, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(resultPath + "final-result");