Skip to content
Snippets Groups Projects
Commit 75fb1e6e authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Use perturbation from an affine displacement as the unknown 

Previously, the unknown was the displacement.  In this setting
it is difficult to implement periodic boundary conditions,
because the displacement itself is never periodic.
We therefore change the code to use the deviation from a
given affine transformation as the new variable.
That is surprisingly simple.  It is also what the ETH group does.
parent 09e51252
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/elasticity/common/trustregionsolver.cc
+ 1
− 0
View file @ 75fb1e6e
...@@ -388,6 +388,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -388,6 +388,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
std::cout << "----------------------------------------------------" << std::endl; std::cout << "----------------------------------------------------" << std::endl;
std::cout << " Trust-Region Step Number: " << i std::cout << " Trust-Region Step Number: " << i
<< ", radius: " << trustRegion.radius() << ", radius: " << trustRegion.radius()
<< std::setprecision(15)
<< ", energy: " << oldEnergy << std::endl; << ", energy: " << oldEnergy << std::endl;
std::cout << "----------------------------------------------------" << std::endl; std::cout << "----------------------------------------------------" << std::endl;
} }
......
src/quasiconvexity-test.cc
+ 106
− 85
View file @ 75fb1e6e
...@@ -55,18 +55,19 @@ const int order = 1; ...@@ -55,18 +55,19 @@ const int order = 1;
using namespace Dune; using namespace Dune;
template<int dim, class field_type = double> template<int dim, class field_type, class ctype=double>
struct KPowerDensity final 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"); p_ = parameters.template get<double>("p");
} }
#warning The type of 'x' should by ctype, not field_type field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
{ {
auto F = derivative + F0_;
auto detF = F.determinant(); auto detF = F.determinant();
auto normSquared = F.frobenius_norm2(); auto normSquared = F.frobenius_norm2();
...@@ -75,15 +76,17 @@ struct KPowerDensity final ...@@ -75,15 +76,17 @@ struct KPowerDensity final
return std::pow(K, p_); return std::pow(K, p_);
} }
FieldMatrix<ctype,dim,dim> F0_;
double p_; double p_;
}; };
template<int dim, class field_type = double> template<int dim, class field_type, class ctype=double>
struct MagicDensity final 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; c_ = 1.0;
if (parameters.hasKey("c")) if (parameters.hasKey("c"))
...@@ -91,8 +94,7 @@ struct MagicDensity final ...@@ -91,8 +94,7 @@ struct MagicDensity final
} }
/** \brief Implement space-dependent coefficients */ /** \brief Implement space-dependent coefficients */
#warning The parameter should actually use ctype, not field_type double c(const Dune::FieldVector<ctype,dim>& x) const
double c(const Dune::FieldVector<field_type,dim>& x) const
{ {
return c_; return c_;
#if 0 // Three circles #if 0 // Three circles
...@@ -105,8 +107,9 @@ struct MagicDensity final ...@@ -105,8 +107,9 @@ struct MagicDensity final
// return (x.two_norm() < 0.9) ? 1 : c_; // 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 det = F.determinant();
auto normSquared = F.frobenius_norm2(); auto normSquared = F.frobenius_norm2();
...@@ -117,21 +120,24 @@ struct MagicDensity final ...@@ -117,21 +120,24 @@ struct MagicDensity final
return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det); return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
} }
FieldMatrix<ctype,dim,dim> F0_;
double c_; double c_;
}; };
template<int dim, class field_type = double> template<int dim, class field_type, class ctype=double>
struct VossDensity final 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"); 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 det = F.determinant();
auto normSquared = F.frobenius_norm2(); auto normSquared = F.frobenius_norm2();
...@@ -147,6 +153,7 @@ struct VossDensity final ...@@ -147,6 +153,7 @@ struct VossDensity final
return h + c_ * std::log( det ); return h + c_ * std::log( det );
} }
FieldMatrix<ctype,dim,dim> F0_;
double c_; double c_;
}; };
...@@ -210,7 +217,7 @@ int main (int argc, char *argv[]) try ...@@ -210,7 +217,7 @@ int main (int argc, char *argv[]) try
upper = parameterSet.get<FieldVector<double,dim> >("upper"); upper = parameterSet.get<FieldVector<double,dim> >("upper");
std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements"); 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 { } else {
std::string path = parameterSet.get<std::string>("path"); std::string path = parameterSet.get<std::string>("path");
...@@ -297,17 +304,15 @@ int main (int argc, char *argv[]) try ...@@ -297,17 +304,15 @@ int main (int argc, char *argv[]) try
blockedInterleaved() blockedInterleaved()
)); ));
if (parameterSet.hasKey("x0")) // Set the underlying affine deformation
{ auto x0 = parameterSet.get<double>("x0");
auto x0 = parameterSet.get<double>("x0");
auto initialDeformation = [&x0](FieldVector<double,dim> x) -> FieldVector<double,dim> FieldMatrix<double,2,2> F0 = { {sqrt(x0), 0},
{ { 0, 1.0 / sqrt(x0)} };
return {x[0]*sqrt(x0), x[1] / 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")) if (parameterSet.hasKey("initialIteratePython"))
{ {
...@@ -364,9 +369,23 @@ int main (int argc, char *argv[]) try ...@@ -364,9 +369,23 @@ int main (int argc, char *argv[]) try
result.push_back(substr); result.push_back(substr);
} }
x[vertex][0] = atof(result[2].c_str()); FieldVector<double,2> pos{atof(result[0].c_str()), atof(result[1].c_str())};
x[vertex][1] = atof(result[3].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++; vertex++;
} }
...@@ -380,7 +399,10 @@ int main (int argc, char *argv[]) try ...@@ -380,7 +399,10 @@ int main (int argc, char *argv[]) try
// Output initial iterate (of homotopy loop) // Output initial iterate (of homotopy loop)
SolutionType identity; 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; SolutionType displacement = x;
displacement -= identity; displacement -= identity;
...@@ -388,8 +410,8 @@ int main (int argc, char *argv[]) try ...@@ -388,8 +410,8 @@ int main (int argc, char *argv[]) try
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement); auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
auto localDisplacementFunction = localFunction(displacementFunction); auto localDisplacementFunction = localFunction(displacementFunction);
// We need to subsample, because VTK cannot natively display real second-order functions // Write the initial iterate
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0)); VTKWriter<GridView> vtkWriter(gridView);
vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim)); vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(resultPath + "quasiconvexity-test_homotopy_0"); vtkWriter.write(resultPath + "quasiconvexity-test_homotopy_0");
...@@ -420,73 +442,73 @@ int main (int argc, char *argv[]) try ...@@ -420,73 +442,73 @@ int main (int argc, char *argv[]) try
if (parameterSet.get<std::string>("energy") == "kpower") 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, elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, FEBasis::LocalView::Tree::FiniteElement,
adouble> >(density); adouble> >(density);
} }
if (parameterSet.get<std::string>("energy") == "neff") // if (parameterSet.get<std::string>("energy") == "neff")
elasticEnergy = std::make_shared<NeffEnergy<GridView, // elasticEnergy = std::make_shared<NeffEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "nefflog") // if (parameterSet.get<std::string>("energy") == "nefflog")
elasticEnergy = std::make_shared<NeffLogEnergy<GridView, // elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "nefflog2") // if (parameterSet.get<std::string>("energy") == "nefflog2")
elasticEnergy = std::make_shared<NeffLogEnergy2<GridView, // elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "neffreduced") // if (parameterSet.get<std::string>("energy") == "neffreduced")
elasticEnergy = std::make_shared<NeffReducedEnergy<GridView, // elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "cosh") // if (parameterSet.get<std::string>("energy") == "cosh")
elasticEnergy = std::make_shared<CosHEnergy<GridView, // elasticEnergy = std::make_shared<CosHEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "expacosh") // if (parameterSet.get<std::string>("energy") == "expacosh")
elasticEnergy = std::make_shared<ExpACosHEnergy<GridView, // elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "fung") // if (parameterSet.get<std::string>("energy") == "fung")
elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView, // elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
if (parameterSet.get<std::string>("energy") == "magic") 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, elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, FEBasis::LocalView::Tree::FiniteElement,
adouble> >(density); adouble> >(density);
} }
if (parameterSet.get<std::string>("energy") == "burkholder") // if (parameterSet.get<std::string>("energy") == "burkholder")
elasticEnergy = std::make_shared<Elasticity::BurkholderEnergy<GridView, // elasticEnergy = std::make_shared<Elasticity::BurkholderEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "w2") // if (parameterSet.get<std::string>("energy") == "w2")
elasticEnergy = std::make_shared<Elasticity::NeffW2Energy<GridView, // elasticEnergy = std::make_shared<Elasticity::NeffW2Energy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "dacorogna") // if (parameterSet.get<std::string>("energy") == "dacorogna")
elasticEnergy = std::make_shared<Elasticity::DacorognaEnergy<GridView, // elasticEnergy = std::make_shared<Elasticity::DacorognaEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, // FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters); // adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "voss") 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, elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement, FEBasis::LocalView::Tree::FiniteElement,
adouble> >(density); adouble> >(density);
...@@ -654,7 +676,7 @@ int main (int argc, char *argv[]) try ...@@ -654,7 +676,7 @@ int main (int argc, char *argv[]) try
for (int j=0; j<dim; j++) for (int j=0; j<dim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]); derivative[j].axpy(localConfiguration[i][j], gradients[i]);
//auto det = derivative.determinant(); derivative += F0;
//std::cout << derivative << std::endl; //std::cout << derivative << std::endl;
deformationDeterminant[localP0View.index(0)] = derivative.determinant(); deformationDeterminant[localP0View.index(0)] = derivative.determinant();
...@@ -741,8 +763,7 @@ int main (int argc, char *argv[]) try ...@@ -741,8 +763,7 @@ int main (int argc, char *argv[]) try
displacement -= identity; displacement -= identity;
auto displacementFunction2 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement); auto displacementFunction2 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions VTKWriter<GridView> vtkWriter(gridView);
//SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
vtkWriter.clear(); vtkWriter.clear();
vtkWriter.addVertexData(displacementFunction2, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim)); vtkWriter.addVertexData(displacementFunction2, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(resultPath + "final-result"); vtkWriter.write(resultPath + "final-result");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment