Newer
Older
#ifndef LFE_ORDER
#define LFE_ORDER 2
#endif
#ifndef GFE_ORDER
#define GFE_ORDER 1
#endif
#ifndef MIXED_SPACE
#if LFE_ORDER != GFE_ORDER
#define MIXED_SPACE 1
#else
#define MIXED_SPACE 0
#endif
#endif
//#define PROJECTED_INTERPOLATION
#include <config.h>
#include <fenv.h>
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
#include <adolc/adouble.h>
#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
#include <adolc/taping.h>

Oliver Sander
committed
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/tuplevector.hh>
#include <dune/common/version.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#if HAVE_DUNE_FOAMGRID
#include <dune/foamgrid/foamgrid.hh>
#endif
#include <dune/istl/multitypeblockvector.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
#include <dune/gfe/assemblers/nonplanarcosseratshellenergy.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/cosseratvtkreader.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#if MIXED_SPACE
#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
#include <dune/gfe/riemannianpnsolver.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
#if HAVE_DUNE_VTK
#include <dune/vtk/vtkreader.hh>
#endif
#include <dune/gmsh4/gmsh4reader.hh>
#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>

Lisa Julia Nebel
committed
// grid dimension
const int dim = GRID_DIM;
const int dimworld = WORLD_DIM;
// Order of the approximation space for the displacement
const int displacementOrder = LFE_ORDER;
// Order of the approximation space for the microrotations
const int rotationOrder = GFE_ORDER;
#if !MIXED_SPACE
static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
#endif
// Image space of the geodesic fe functions
using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
int main (int argc, char *argv[]) try
{
// initialize MPI, finalize is done automatically on exit
Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
if (mpiHelper.rank()==0) {
std::cout << "DISPLACEMENTORDER = " << displacementOrder << ", ROTATIONORDER = " << rotationOrder << std::endl;
std::cout << "dim = " << dim << ", dimworld = " << dimworld << std::endl;
#if MIXED_SPACE
std::cout << "MIXED_SPACE = 1" << std::endl;
std::cout << "MIXED_SPACE = 0" << std::endl;
// Check for appropriate number of command line arguments
if (argc < 3)
DUNE_THROW(Exception, "Usage: ./cosserat-continuum <python path> <parameter file>");
// Start Python interpreter
Python::start();
auto pyMain = Python::main();
Python::run("import math");
//feenableexcept(FE_INVALID);
Python::runStream()
<< std::endl << "import sys"
<< std::endl << "sys.path.append('" << argv[1] << "')"
<< std::endl;
// parse data file
auto pyModule = pyMain.import(argv[2]);
// Get main parameter set
ParameterTree parameterSet;
pyModule.get("parameterSet").toC(parameterSet);
ParameterTreeParser::readOptions(argc, argv, parameterSet);
// read solver settings
const int numLevels = parameterSet.get<int>("numLevels");
int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
const double tolerance = parameterSet.get<double>("tolerance");
const int maxSolverSteps = parameterSet.get<int>("maxSolverSteps");
const double initialRegularization = parameterSet.get<double>("initialRegularization", 1000);
const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
const int multigridIterations = parameterSet.get<int>("numIt");
const int nu1 = parameterSet.get<int>("nu1");
const int nu2 = parameterSet.get<int>("nu2");
const int mu = parameterSet.get<int>("mu");
const int baseIterations = parameterSet.get<int>("baseIt");
const double mgTolerance = parameterSet.get<double>("mgTolerance");
const double baseTolerance = parameterSet.get<double>("baseTolerance");
const bool instrumented = parameterSet.get<bool>("instrumented");
const bool adolcScalarMode = parameterSet.get<bool>("adolcScalarMode", false);
const std::string resultPath = parameterSet.get("resultPath", "");
using namespace Dune::Indices;
using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
std::vector<Rotation<double,3> > >;
// ///////////////////////////////////////
// Create the grid
// ///////////////////////////////////////
typedef std::conditional<dim==dimworld,UGGrid<dim>, FoamGrid<dim,dimworld> >::type GridType;
static_assert(dim==dimworld, "FoamGrid needs to be installed to allow problems with dim != dimworld.");
typedef UGGrid<dim> GridType;
std::shared_ptr<GridType> grid;
GridFactory<GridType> factory;
Gmsh4::LagrangeGridCreator creator{factory};
FieldVector<double,dimworld> lower(0), upper(1);
std::string structuredGridType = parameterSet["structuredGrid"];
if (structuredGridType == "cube") {
if (dim!=dimworld)
DUNE_THROW(GridError, "Please use FoamGrid and read in a grid for problems with dim != dimworld.");
lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
upper = parameterSet.get<FieldVector<double,dimworld> >("upper");

Oliver Sander
committed
auto elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
} else {
std::string path = parameterSet.get<std::string>("path", "");
std::string gridFile = parameterSet.get<std::string>("gridFile");
// Guess the grid file format by looking at the file name suffix
auto dotPos = gridFile.rfind('.');
if (dotPos == std::string::npos)
DUNE_THROW(IOError, "Could not determine grid input file format");
std::string suffix = gridFile.substr(dotPos, gridFile.length()-dotPos);
if (suffix == ".msh") {
Gmsh4Reader reader{creator};
reader.read(path + "/" + gridFile);
grid = factory.createGrid();
} else if (suffix == ".vtu" or suffix == ".vtp")
grid = VtkReader<GridType>::createGridFromFile(path + "/" + gridFile);
DUNE_THROW(NotImplemented, "Please install dune-vtk for VTK reading support!");
grid->globalRefine(numLevels-1);
grid->loadBalance();
if (mpiHelper.rank()==0)
std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
using namespace Dune::Functions::BasisFactory;
const int dimRotation = Rotation<double,3>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
composite(
power<3>(
lagrange<displacementOrder>()
power<dimRotation>(
lagrange<rotationOrder>()
));
using CompositeBasis = decltype(compositeBasis);
auto deformationPowerBasis = makeBasis(
gridView,
power<3>(
lagrange<displacementOrder>()
));
auto orientationPowerBasis = makeBasis(
gridView,
power<3>(
power<3>(
lagrange<rotationOrder>()

Lisa Julia Nebel
committed
)

Lisa Julia Nebel
committed
BlockVector<FieldVector<double,dimworld> > identityDeformation(compositeBasis.size({0}));
auto identityDeformationBasis = makeBasis(
gridView,
power<dimworld>(
lagrange<displacementOrder>()
));
Dune::Functions::interpolate(identityDeformationBasis, identityDeformation, [&](FieldVector<double,dimworld> x){
return x;
});
BlockVector<FieldVector<double,dimworld> > identityRotation(compositeBasis.size({1}));
auto identityRotationBasis = makeBasis(
gridView,
power<dimworld>(
lagrange<rotationOrder>()
));
Dune::Functions::interpolate(identityRotationBasis, identityRotation, [&](FieldVector<double,dimworld> x){
return x;
});
typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
BitSetVector<1> neumannVertices(gridView.size(dim), false);
BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
const GridView::IndexSet& indexSet = gridView.indexSet();
// 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(")");
auto pythonDirichletVertices = Python::make_function<FieldVector<bool,3> >(Python::evaluate(lambda));
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletRotationVerticesPredicate") + std::string(")");
auto pythonOrientationDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
// Same for the Neumann boundary
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));

Oliver Sander
committed
for (auto&& vertex : vertices(gridView))
{
bool isNeumann = pythonNeumannVertices(vertex.geometry().corner(0));
neumannVertices[indexSet.index(vertex)] = isNeumann;
}
BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces()
<< " faces and " << neumannVertices.count() << " degrees of freedom.\n";
BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
#if DUNE_VERSION_GTE(DUNE_FUFEM, 2, 10)
Fufem::markBoundaryPatchDofs(neumannBoundary,deformationFEBasis,neumannNodes);
#else
constructBoundaryDofs(neumannBoundary,deformationFEBasis,neumannNodes);
for (size_t i=0; i<deformationFEBasis.size(); i++) {
FieldVector<bool,3> isDirichlet;
isDirichlet = pythonDirichletVertices(identityDeformation[i]);
for (size_t j=0; j<3; j++)
deformationDirichletDofs[i][j] = isDirichlet[j];
}
for (size_t i=0; i<orientationFEBasis.size(); i++) {
bool isDirichletOrientation = pythonOrientationDirichletVertices(identityRotation[i]);
for (size_t j=0; j<3; j++)
orientationDirichletDofs[i][j] = isDirichletOrientation;
}
std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << deformationDirichletDofs.count() << " degrees of freedom.\n";
std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary (orientation) has " << orientationDirichletDofs.count() << " degrees of freedom.\n";
// //////////////////////////
// Initial iterate
// //////////////////////////
SolutionType x;
x[_0].resize(compositeBasis.size({0}));
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
std::vector<FieldVector<double,3> > v;
Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
for (size_t i=0; i<x[_0].size(); i++)
x[_0][i] = v[i];
x[_1].resize(compositeBasis.size({1}));
#if !MIXED_SPACE
if (parameterSet.hasKey("startFromFile"))
{
std::shared_ptr<GridType> initialIterateGrid;
if (parameterSet.get<bool>("structuredGrid"))
std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
initialIterateGrid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
}
else
{
std::string path = parameterSet.get<std::string>("path");
std::string initialIterateGridFilename = parameterSet.get<std::string>("initialIterateGridFilename");
initialIterateGrid = std::shared_ptr<GridType>(Gmsh4Reader<GridType>::createGridFromFile(path + "/" + initialIterateGridFilename));
}
std::vector<TargetSpace> initialIterate;
GFE::CosseratVTKReader::read(initialIterate, parameterSet.get<std::string>("initialIterateFilename"));
typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, 2> InitialBasis;
InitialBasis initialBasis(initialIterateGrid->leafGridView());
#ifdef PROJECTED_INTERPOLATION
using LocalInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
using LocalInterpolationRule = LocalGeodesicFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
#endif
GFE::EmbeddedGlobalGFEFunction<InitialBasis,LocalInterpolationRule,TargetSpace> initialFunction(initialBasis,initialIterate);
auto powerBasis = makeBasis(
gridView,
power<7>(
lagrange<displacementOrder>()
));
std::vector<FieldVector<double,7> > v;
//TODO: Interpolate does not work with an GFE:EmbeddedGlobalGFEFunction
//Dune::Functions::interpolate(powerBasis,v,initialFunction);
for (size_t i=0; i<x.size(); i++) {
auto vTargetSpace = TargetSpace(v[i]);
x[_0][i] = vTargetSpace[_0];
x[_1][i] = vTargetSpace[_1];
////////////////////////////////////////////////////////
// Main homotopy loop
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
if (dim == 2 && dimworld == 2) {
CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
orientationFEBasis,x[_1],
resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
} else if (dim == 2 && dimworld == 3) {
#if MIXED_SPACE
CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x[_0], resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
#endif
} else if (dim == 3 && dimworld == 3) {
BlockVector<FieldVector<double,dimworld> > displacement(x[_0].size());
for (size_t i=0; i<x[_0].size(); i++) {
for (int j = 0; j < 3; j ++)
displacement[i][j] = x[_0][i][j] - identityDeformation[i][j];
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(deformationPowerBasis, displacement);
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write(resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
}
for (int i=0; i<numHomotopySteps; i++) {
double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
if (mpiHelper.rank()==0)
std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
// ////////////////////////////////////////////////////////////
// Create an assembler for the energy functional
// ////////////////////////////////////////////////////////////
const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
FieldVector<double,3> neumannValues {0,0,0};
if (parameterSet.hasKey("neumannValues"))
neumannValues = parameterSet.get<FieldVector<double,3> >("neumannValues");

Lisa Julia Nebel
committed
auto neumannFunction = [&]( FieldVector<double,dimworld> ) {
auto nV = neumannValues;
nV *= (-homotopyParameter);
return nV;
};

Lisa Julia Nebel
committed
Python::Reference volumeLoadClass = Python::import(parameterSet.get<std::string>("volumeLoadPythonFunction", "zero-volume-load"));
Python::Callable volumeLoadCallable = volumeLoadClass.get("VolumeLoad");
// Call a constructor
Python::Reference volumeLoadPythonObject = volumeLoadCallable(homotopyParameter);
// Extract object member functions as Dune functions
auto volumeLoad = Python::make_function<FieldVector<double,3> > (volumeLoadPythonObject.get("volumeLoad"));
std::cout << "Material parameters:" << std::endl;
materialParameters.report();
////////////////////////////////////////////////////////
// Set Dirichlet values
////////////////////////////////////////////////////////
// Dirichlet boundary value function, depending on the homotopy parameter
Python::Callable dirichletValuesPythonClass = pyModule.get("DirichletValues");
// Construct with a particular value of the homotopy parameter
Python::Reference dirichletValuesPythonObject = dirichletValuesPythonClass(homotopyParameter);
// Extract object member functions as Dune functions
auto deformationDirichletValues = Python::make_function<FieldVector<double,3> > (dirichletValuesPythonObject.get("deformation"));
auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> > (dirichletValuesPythonObject.get("orientation"));
BlockVector<FieldVector<double,3> > ddV;
Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues);
BlockVector<FieldMatrix<double,3,3> > dOV;
Dune::Functions::interpolate(orientationPowerBasis, dOV, orientationDirichletValues);
for (std::size_t i = 0; i < compositeBasis.size({0}); i++) {
FieldVector<double,3> x0i = x[_0][i].globalCoordinates();
for (int j=0; j<3; j++) {
if (deformationDirichletDofs[i][j])
x0i[j] = ddV[i][j];
}
x[_0][i] = x0i;
}
for (std::size_t i = 0; i < compositeBasis.size({1}); i++)
if (orientationDirichletDofs[i][0])
x[_1][i].set(dOV[i]);
if (dim==dimworld) {
auto localCosseratEnergy = std::make_shared<CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> > (materialParameters,
&neumannBoundary,
neumannFunction,
volumeLoad);
LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(localCosseratEnergy,
adolcScalarMode);
MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
#if MIXED_SPACE
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion")
{
MixedRiemannianTrustRegionSolver<GridType,
CompositeBasis,
DeformationFEBasis, RealTuple<double,3>,
OrientationFEBasis, Rotation<double,3> > solver;
solver.setup(*grid,
&mixedAssembler,
deformationFEBasis,
orientationFEBasis,
x,
deformationDirichletDofs,
orientationDirichletDofs, tolerance,
maxSolverSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
}
else
{
//Create BitVector matching the tangential space
const int dimRotationTangent = Rotation<double,3>::TangentVector::dimension;
using VectorForBit = MultiTypeBlockVector<std::vector<FieldVector<double,3> >, std::vector<FieldVector<double,dimRotationTangent> > >;
using BitVector = Solvers::DefaultBitVector_t<VectorForBit>;
BitVector dirichletDofs;
dirichletDofs[_0].resize(compositeBasis.size({0}));
dirichletDofs[_1].resize(compositeBasis.size({1}));
for (size_t i = 0; i < compositeBasis.size({0}); i++) {
for (size_t j = 0; j < 3; j++)
dirichletDofs[_0][i][j] = deformationDirichletDofs[i][j];
}
for (size_t i = 0; i < compositeBasis.size({1}); i++) {
for (int j = 0; j < dimRotationTangent; j++)
dirichletDofs[_1][i][j] = orientationDirichletDofs[i][j];
}
GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,3>, OrientationFEBasis, Rotation<double,3>, BitVector> solver;
solver.setup(*grid,
&mixedAssembler,
x,
dirichletDofs,
tolerance,
maxSolverSteps,
initialRegularization,
instrumented);
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
}
#else
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a ProductManifold.
//Therefore, x and the dirichletDofs are converted to a ProductManifold structure, as well as the Hessian and Gradient that are returned by the assembler
std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
for (std::size_t i = 0; i < compositeBasis.size({0}); i++) {
for (int j = 0; j < 3; j ++) { // Displacement part
xTargetSpace[i][_0].globalCoordinates()[j] = x[_0][i][j];
dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
}
xTargetSpace[i][_1] = x[_1][i]; // Rotation part
for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
}
using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace>;
GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
xTargetSpace,
dirichletDofsTargetSpace,
tolerance,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
} else {
RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
xTargetSpace,
dirichletDofsTargetSpace,
tolerance,
maxSolverSteps,
initialRegularization,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
}
for (std::size_t i = 0; i < xTargetSpace.size(); i++) {
x[_0][i] = xTargetSpace[i][_0];
x[_1][i] = xTargetSpace[i][_1];
} else { //dim != dimworld
using StiffnessType = LocalGeodesicFEADOLCStiffness<CompositeBasis, TargetSpace>;
std::shared_ptr<StiffnessType> localGFEStiffness;
#if HAVE_DUNE_CURVEDGEOMETRY && WORLD_DIM == 3 && GRID_DIM == 2
auto localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<CompositeBasis, 3, adouble, decltype(creator)> >(materialParameters,
&creator,
&neumannBoundary,
neumannFunction,
volumeLoad);
localGFEStiffness = std::make_shared<StiffnessType>(localCosseratEnergy, adolcScalarMode);
MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEStiffness);

Lisa Julia Nebel
committed
#if MIXED_SPACE
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
MixedRiemannianTrustRegionSolver<GridType,
CompositeBasis,
DeformationFEBasis, RealTuple<double,3>,
OrientationFEBasis, Rotation<double,3> > solver;
solver.setup(*grid,
&mixedAssembler,
deformationFEBasis,
orientationFEBasis,
x,
deformationDirichletDofs,
orientationDirichletDofs, tolerance,
maxSolverSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
#else
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a ProductManifold.
//Therefore, x and the dirichletDofs are converted to a ProductManifold structure, as well as the Hessian and Gradient that are returned by the assembler
std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
for (std::size_t i = 0; i < compositeBasis.size({0}); i++) {
for (int j = 0; j < 3; j ++) { // Displacement part
xTargetSpace[i][_0].globalCoordinates()[j] = x[_0][i][j];
dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
}
xTargetSpace[i][_1] = x[_1][i]; // Rotation part
for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
}
using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace>;
GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
xTargetSpace,
dirichletDofsTargetSpace,
tolerance,

Lisa Julia Nebel
committed
maxSolverSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
} else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
xTargetSpace,
dirichletDofsTargetSpace,
tolerance,
maxSolverSteps,
initialRegularization,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
}

Lisa Julia Nebel
committed
for (std::size_t i = 0; i < xTargetSpace.size(); i++) {
x[_0][i] = xTargetSpace[i][_0];
x[_1][i] = xTargetSpace[i][_1];
}
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
if (dim == 2 && dimworld == 2) {
CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
orientationFEBasis,x[_1],
resultPath + "cosserat_homotopy_" + iAsAscii.str());
} else if (dim == 2 && dimworld == 3) {
#if MIXED_SPACE
CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x[_0], resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
#else
CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
} else if (dim == 3 && dimworld == 3) {
BlockVector<FieldVector<double,dimworld> > displacement(x[_0].size());
for (size_t i=0; i<x[_0].size(); i++) {
for (int j = 0; j < 3; j ++)
displacement[i][j] = x[_0][i][j] - identityDeformation[i][j];
}
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(deformationPowerBasis, displacement);
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write(resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
// //////////////////////////////
// Output result
// //////////////////////////////
#if !MIXED_SPACE
// Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
// This data may be used by other applications measuring the discretization error
BlockVector<TargetSpace::CoordinateType> xEmbedded(compositeBasis.size({0}));
for (size_t i=0; i<compositeBasis.size({0}); i++)
for (size_t j=0; j<TargetSpace::TangentVector::dimension; j++)
xEmbedded[i][j] = j < 3 ? x[_0][i][j] : x[_1][i].globalCoordinates()[j-3];
std::ofstream outFile("cosserat-continuum-result-" + std::to_string(numLevels) + ".data", std::ios_base::binary);
MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
if (parameterSet.get<bool>("computeAverageNeumannDeflection", false) and parameterSet.hasKey("neumannValues")) {
// finally: compute the average deformation of the Neumann boundary
// That is what we need for the locking tests
FieldVector<double,3> averageDef(0);
for (size_t i=0; i<x[_0].size(); i++)
if (neumannNodes[i][0])
averageDef += x[_0][i].globalCoordinates();
averageDef /= neumannNodes.count();
if (mpiHelper.rank()==0 and parameterSet.hasKey("neumannValues"))
std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << " "
<< ", average deflection: " << averageDef << std::endl;
}
catch (Exception& e)
std::cout << e.what() << std::endl;