Newer
Older

Lisa Julia Nebel
committed
#include <config.h>
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
#include <adolc/adouble.h>
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>

Lisa Julia Nebel
committed
#include <dune/gfe/assemblers/mixedgfeassembler.hh>

Lisa Julia Nebel
committed
#include <dune/gfe/assemblers/sumenergy.hh>
#include <dune/gfe/assemblers/localintegralenergy.hh>
#include <dune/gfe/densities/bulkcosseratdensity.hh>
#include <dune/gfe/localgeodesicfefunction.hh>

Lisa Julia Nebel
committed
#include <dune/gfe/mixedriemanniantrsolver.hh>

Lisa Julia Nebel
committed
#include <dune/gfe/neumannenergy.hh>

Lisa Julia Nebel
committed
// grid dimension
const int dim = 3;
// Order of the approximation space for the displacement
const int displacementOrder = 2;
// Order of the approximation space for the microrotations
const int rotationOrder = 1;
using namespace Dune;
int main (int argc, char *argv[])
{
MPIHelper::instance(argc, argv);
using SolutionType = TupleVector<std::vector<RealTuple<double,dim> >, std::vector<Rotation<double,dim> > >;
// solver settings
const double tolerance = 1e-6;
const int maxSolverSteps = 1000;
const double initialTrustRegionRadius = 1;
const int multigridIterations = 200;
const int baseIterations = 100;
const double mgTolerance = 1e-10;
const double baseTolerance = 1e-8;
// ///////////////////////////////////////
// Create the grid
// ///////////////////////////////////////
using GridType = UGGrid<dim>;
//Create a grid with 2 elements, one element will get refined to create different element types, the other one stays a cube

Lisa Julia Nebel
committed
std::shared_ptr<GridType> grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {2,1,1});
// Refine once
for (auto&& e : elements(grid->leafGridView())) {

Lisa Julia Nebel
committed
bool refineThisElement = false;
for (int i = 0; i < e.geometry().corners(); i++) {
refineThisElement = refineThisElement || (e.geometry().corner(i)[0] > 0.99);

Lisa Julia Nebel
committed
}
grid->mark(refineThisElement ? 1 : 0, e);
}
grid->adapt();
grid->loadBalance();
using GridView = GridType::LeafGridView;
GridView gridView = grid->leafGridView();
// ///////////////////////////////////////////
// Construct all needed function space bases
// ///////////////////////////////////////////
using namespace Dune::Indices;

Lisa Julia Nebel
committed
using namespace Functions::BasisFactory;
const int dimRotation = Rotation<double,dim>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
composite(
power<dim>(
lagrange<displacementOrder>()
),

Lisa Julia Nebel
committed
power<dimRotation>(
lagrange<rotationOrder>()
)
));

Lisa Julia Nebel
committed
using CompositeBasis = decltype(compositeBasis);
using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
using OrientationFEBasis = Functions::LagrangeBasis<GridView,rotationOrder>;
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
// ////////////////////////////////////////////
// Determine Dirichlet dofs
// ////////////////////////////////////////////
// This identityBasis is only needed to interpolate the identity for both the displacement and the rotation
auto identityBasis = makeBasis(
gridView,
composite(
power<dim>(
lagrange<displacementOrder>()
),

Lisa Julia Nebel
committed
power<dim>(
lagrange<rotationOrder>()
)
));

Lisa Julia Nebel
committed
auto deformationPowerBasis = Functions::subspaceBasis(identityBasis, _0);
auto rotationPowerBasis = Functions::subspaceBasis(identityBasis, _1);
MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >,BlockVector<FieldVector<double,dim> > > identity;
Functions::interpolate(deformationPowerBasis, identity, [&](FieldVector<double,dim> x){
return x;
});
Functions::interpolate(rotationPowerBasis, identity, [&](FieldVector<double,dim> x){
return x;
});

Lisa Julia Nebel
committed
BitSetVector<dim> deformationDirichletDofs(deformationFEBasis.size(), false);
BitSetVector<dim> orientationDirichletDofs(orientationFEBasis.size(), false);
const GridView::IndexSet& indexSet = gridView.indexSet();
// Make predicate function that computes which vertices are on the Dirichlet boundary, based on the vertex positions.
auto isDirichlet = [](FieldVector<double,dim> coordinate)
{
return coordinate[0] < 0.01;
};

Lisa Julia Nebel
committed
for (size_t i=0; i<deformationFEBasis.size(); i++) {
bool isDirichletDeformation = isDirichlet(identity[_0][i]);
for (size_t j=0; j<dim; j++)
deformationDirichletDofs[i][j] = isDirichletDeformation;

Lisa Julia Nebel
committed
}
for (size_t i=0; i<orientationFEBasis.size(); i++) {
bool isDirichletOrientation = isDirichlet(identity[_1][i]);
for (size_t j=0; j<dim; j++)
orientationDirichletDofs[i][j] = isDirichletOrientation;

Lisa Julia Nebel
committed
}
// /////////////////////////////////////////
// Determine Neumann dofs and values
// /////////////////////////////////////////

Lisa Julia Nebel
committed
std::function<bool(FieldVector<double,dim>)> isNeumann = [](FieldVector<double,dim> coordinate) {
return coordinate[0] > 0.99;
};

Lisa Julia Nebel
committed
BitSetVector<1> neumannVertices(gridView.size(dim), false);

Lisa Julia Nebel
committed
for (auto&& vertex : vertices(gridView))
neumannVertices[indexSet.index(vertex)] = isNeumann(vertex.geometry().corner(0));

Lisa Julia Nebel
committed
auto neumannBoundary = std::make_shared<BoundaryPatch<GridType::LeafGridView> >(gridView, neumannVertices);

Lisa Julia Nebel
committed
FieldVector<double,dim> values_ = {0,0,0};
auto neumannFunction = [&](FieldVector<double, dim>){
return values_;
};

Lisa Julia Nebel
committed
// //////////////////////////
// Initial iterate
// //////////////////////////
SolutionType x;
x[_0].resize(compositeBasis.size({0}));
x[_1].resize(compositeBasis.size({1}));
// ////////////////////////////
// Parameters for the problem
// ////////////////////////////
ParameterTree parameters;
parameters["mu"] = "1";
parameters["lambda"] = "1";
parameters["mu_c"] = "1";
parameters["L_c"] = "0.1";
parameters["curvatureType"] = "wryness";

Lisa Julia Nebel
committed
parameters["q"] = "2";
parameters["b1"] = "1";
parameters["b2"] = "1";
parameters["b3"] = "1";
// ////////////////////////////
// Create an assembler
// ////////////////////////////
using RigidBodyMotion = GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> >;
using ARigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>,Rotation<adouble,dim> >;
auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,dim>,Rotation<adouble,dim> > >();
auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
// Select which type of geometric interpolation to use
using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
auto bulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<FieldVector<double,dim>,adouble> >(parameters);
auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ARigidBodyMotion> >(bulkCosseratDensity);
sumEnergy->addLocalEnergy(bulkCosseratEnergy);
sumEnergy->addLocalEnergy(neumannEnergy);

Lisa Julia Nebel
committed
LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> localGFEADOLCStiffness(sumEnergy);
MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssembler(compositeBasis, localGFEADOLCStiffness);

Lisa Julia Nebel
committed
MixedRiemannianTrustRegionSolver<GridType,
CompositeBasis,
DeformationFEBasis, RealTuple<double,dim>,
OrientationFEBasis, Rotation<double,dim> > solver;

Lisa Julia Nebel
committed
solver.setup(*grid,
&mixedAssembler,
deformationFEBasis,
orientationFEBasis,
x,
deformationDirichletDofs,
orientationDirichletDofs,
tolerance,
maxSolverSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
3, 3, 1, // Multigrid V-cycle
baseIterations,
baseTolerance,
false);

Lisa Julia Nebel
committed
solver.setInitialIterate(x);
solver.solve();
// //////////////////////////////////////////////
// Check if solver returns the expected values
// //////////////////////////////////////////////
size_t expectedFinalIteration = 14;
double expectedEnergy = 0.67188585;
if (solver.getStatistics().finalIteration != expectedFinalIteration)
{
std::cerr << "Trust-region solver did " << solver.getStatistics().finalIteration+1
<< " iterations, instead of the expected '" << expectedFinalIteration+1 << "'!" << std::endl;
return 1;

Lisa Julia Nebel
committed
}
if ( std::abs(solver.getStatistics().finalEnergy - expectedEnergy) > 1e-7)
{
std::cerr << std::setprecision(9);
std::cerr << "Final energy is " << solver.getStatistics().finalEnergy
<< " but '" << expectedEnergy << "' was expected!" << std::endl;
return 1;

Lisa Julia Nebel
committed
}
return 0;
}