Commit f9aabbf8 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

WIP

parent 1fb540e5
......@@ -107,6 +107,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);
......@@ -439,6 +440,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);
......@@ -477,8 +499,8 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (this->verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
x_ = newIterate;
break;
// x_ = newIterate;
// break;
}
// //////////////////////////////////////////////
......
......@@ -100,6 +100,8 @@ public:
SolutionType getSol() const {return x_;}
std::string iterateNamePrefix_;
protected:
/** \brief The grid */
......
......@@ -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> > >
......
......@@ -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> > >
......
......@@ -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> > >
......
......@@ -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> > >
......
......@@ -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> > >
......
#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
......@@ -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> > >
......
......@@ -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
......
......@@ -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
......@@ -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)
......
......@@ -92,6 +92,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);
......
#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
......@@ -144,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));
......@@ -165,6 +269,18 @@ int main (int argc, char *argv[]) try
for (int j=0; j<dim; j++)
dirichletDofs[i][j] = true;
// 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;
}
}
// //////////////////////////
// Initial iterate
// //////////////////////////
......@@ -179,14 +295,28 @@ 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"))
{
......@@ -240,10 +370,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,
......@@ -274,10 +412,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!");
......@@ -308,7 +477,8 @@ int main (int argc, char *argv[]) try
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance
baseTolerance,
1.0
);
////////////////////////////////////////////////////////
......@@ -330,9 +500,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