diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 9d6c4d101d57f8aa842d5a7f2686009b9dc666c1..45a09e9929fa2e21f5abe2511d677218954de0e0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,7 +2,6 @@ set(programs compute-disc-error
cosserat-continuum
gradient-flow
harmonicmaps
- mixed-cosserat-continuum
rod-eoc)
# rodobstacle)
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 7552c374082e5699a714c7e3ae83eec718beb5cc..9029c25248891d0c0e8614d1eca0ad72de6c7911 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -15,6 +15,7 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
+#include <dune/common/tuplevector.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
@@ -26,6 +27,8 @@
#endif
#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
@@ -38,6 +41,7 @@
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
#include <dune/gfe/cosseratenergystiffness.hh>
#include <dune/gfe/nonplanarcosseratshellenergy.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
@@ -47,19 +51,30 @@
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/vertexnormals.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedriemanniantrsolver.hh>
// grid dimension
const int dim = 2;
const int dimworld = 2;
-// Order of the approximation space
-const int order = 2;
+//#define MIXED_SPACE
+
+// Order of the approximation space for the displacement
+const int displacementOrder = 2;
+
+// Order of the approximation space for the microrotations
+const int rotationOrder = 2;
+
+#ifndef MIXED_SPACE
+static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
// Image space of the geodesic fe functions
typedef RigidBodyMotion<double,3> TargetSpace;
// Tangent vector of the image space
const int blocksize = TargetSpace::TangentVector::dimension;
+#endif
using namespace Dune;
@@ -118,7 +133,13 @@ int main (int argc, char *argv[]) try
<< std::endl << "sys.path.append('/home/sander/dune/dune-gfe/problems/')"
<< std::endl;
+ using namespace TypeTree::Indices;
+#ifdef MIXED_SPACE
+ using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
+ std::vector<Rotation<double,3> > >;
+#else
typedef std::vector<TargetSpace> SolutionType;
+#endif
// parse data file
ParameterTree parameterSet;
@@ -194,11 +215,35 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, order> FEBasis;
+ using namespace Dune::Functions::BasisFactory;
+#ifdef MIXED_SPACE
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ lagrange<displacementOrder>(),
+ lagrange<rotationOrder>()
+ )
+ );
+
+ typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
+
+ DeformationFEBasis deformationFEBasis(gridView);
+ OrientationFEBasis orientationFEBasis(gridView);
+
+ // Construct fufem-style function space bases to ease the transition to dune-functions
+ typedef DuneFunctionsBasis<DeformationFEBasis> FufemDeformationFEBasis;
+ FufemDeformationFEBasis fufemDeformationFEBasis(deformationFEBasis);
+
+ typedef DuneFunctionsBasis<OrientationFEBasis> FufemOrientationFEBasis;
+ FufemOrientationFEBasis fufemOrientationFEBasis(orientationFEBasis);
+#else
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, displacementOrder> FEBasis;
FEBasis feBasis(gridView);
typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
FufemFEBasis fufemFeBasis(feBasis);
+#endif
// /////////////////////////////////////////
// Read Dirichlet values
@@ -235,7 +280,28 @@ int main (int argc, char *argv[]) try
if (mpiHelper.rank()==0)
std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+#ifdef MIXED_SPACE
+ BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,fufemDeformationFEBasis,deformationDirichletNodes);
+
+ BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
+ constructBoundaryDofs(neumannBoundary,fufemDeformationFEBasis,neumannNodes);
+
+ BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
+ for (size_t i=0; i<deformationFEBasis.size(); i++)
+ if (deformationDirichletNodes[i][0])
+ for (int j=0; j<3; j++)
+ deformationDirichletDofs[i][j] = true;
+ BitSetVector<1> orientationDirichletNodes(orientationFEBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,fufemOrientationFEBasis,orientationDirichletNodes);
+
+ BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
+ for (size_t i=0; i<orientationFEBasis.size(); i++)
+ if (orientationDirichletNodes[i][0])
+ for (int j=0; j<3; j++)
+ orientationDirichletDofs[i][j] = true;
+#else
BitSetVector<1> dirichletNodes(feBasis.size(), false);
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
@@ -247,11 +313,38 @@ int main (int argc, char *argv[]) try
if (dirichletNodes[i][0])
for (int j=0; j<5; j++)
dirichletDofs[i][j] = true;
+#endif
// //////////////////////////
// Initial iterate
// //////////////////////////
+#ifdef MIXED_SPACE
+ SolutionType x;
+
+ x[_0].resize(deformationFEBasis.size());
+
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+
+ std::vector<FieldVector<double,3> > v;
+ ::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
+
+ for (size_t i=0; i<x[_0].size(); i++)
+ x[_0][i] = v[i];
+
+ x[_1].resize(orientationFEBasis.size());
+#if 0
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+
+ std::vector<FieldVector<double,3> > v;
+ Functions::interpolate(feBasis, v, pythonInitialDeformation);
+
+ for (size_t i=0; i<x.size(); i++)
+ xDisp[i] = v[i];
+#endif
+#else
SolutionType x(feBasis.size());
if (parameterSet.hasKey("startFromFile"))
@@ -285,6 +378,7 @@ int main (int argc, char *argv[]) try
std::vector<FieldVector<double,7> > v;
Dune::Functions::interpolate(feBasis,v,initialFunction);
+ DUNE_THROW(NotImplemented, "Replace scalar basis by power basis!");
for (size_t i=0; i<x.size(); i++)
x[i] = TargetSpace(v[i]);
@@ -299,13 +393,20 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<x.size(); i++)
x[i].r = v[i];
}
+#endif
////////////////////////////////////////////////////////
// Main homotopy loop
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
+#ifdef MIXED_SPACE
+ CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
+ orientationFEBasis,x[_1],
+ resultPath + "mixed-cosserat_homotopy_0");
+#else
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_0");
+#endif
for (int i=0; i<numHomotopySteps; i++) {
@@ -335,6 +436,21 @@ int main (int argc, char *argv[]) try
}
// Assembler using ADOL-C
+#ifdef MIXED_SPACE
+ CosseratEnergyLocalStiffness<decltype(compositeBasis),
+ 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
+ &neumannBoundary,
+ neumannFunction,
+ nullptr);
+
+ MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
+ RealTuple<double,3>,
+ Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
+
+ MixedGFEAssembler<decltype(compositeBasis),
+ RealTuple<double,3>,
+ Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
+#else
using LocalEnergyBase = LocalGeodesicFEStiffness<FEBasis,RigidBodyMotion<adouble,3> >;
std::shared_ptr<LocalEnergyBase> cosseratEnergyADOLCLocalStiffness;
@@ -360,16 +476,33 @@ int main (int argc, char *argv[]) try
TargetSpace> localGFEADOLCStiffness(cosseratEnergyADOLCLocalStiffness.get());
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
+#endif
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
+#ifdef MIXED_SPACE
+ MixedRiemannianTrustRegionSolver<GridType,
+ decltype(compositeBasis),
+ DeformationFEBasis, RealTuple<double,3>,
+ OrientationFEBasis, Rotation<double,3> > solver;
+#else
RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+#endif
solver.setup(*grid,
&assembler,
+#ifdef MIXED_SPACE
+ deformationFEBasis,
+ orientationFEBasis,
+#endif
x,
+#ifdef MIXED_SPACE
+ deformationDirichletDofs,
+ orientationDirichletDofs,
+#else
dirichletDofs,
+#endif
tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
@@ -398,9 +531,21 @@ int main (int argc, char *argv[]) try
PythonFunction<FieldVector<double,dimworld>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
std::vector<FieldVector<double,3> > ddV;
- ::Functions::interpolate(fufemFeBasis, ddV, deformationDirichletValues, dirichletDofs);
-
std::vector<FieldMatrix<double,3,3> > dOV;
+
+#ifdef MIXED_SPACE
+ ::Functions::interpolate(fufemDeformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+ ::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
+
+ for (size_t j=0; j<x[_0].size(); j++)
+ if (deformationDirichletNodes[j][0])
+ x[_0][j] = ddV[j];
+
+ for (size_t j=0; j<x[_1].size(); j++)
+ if (orientationDirichletNodes[j][0])
+ x[_1][j].set(dOV[j]);
+#else
+ ::Functions::interpolate(fufemFeBasis, ddV, deformationDirichletValues, dirichletDofs);
::Functions::interpolate(fufemFeBasis, dOV, orientationDirichletValues, dirichletDofs);
for (size_t j=0; j<x.size(); j++)
@@ -409,6 +554,7 @@ int main (int argc, char *argv[]) try
x[j].r = ddV[j];
x[j].q.set(dOV[j]);
}
+#endif
// /////////////////////////////////////////////////////
// Solve!
@@ -422,7 +568,13 @@ int main (int argc, char *argv[]) try
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
+#ifdef MIXED_SPACE
+ CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
+ orientationFEBasis,x[_1],
+ resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
+#else
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_" + iAsAscii.str());
+#endif
}
@@ -430,6 +582,7 @@ int main (int argc, char *argv[]) try
// Output result
// //////////////////////////////
+#ifndef 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(x.size());
@@ -439,13 +592,20 @@ int main (int argc, char *argv[]) try
std::ofstream outFile("cosserat-continuum-result-" + std::to_string(numLevels) + ".data", std::ios_base::binary);
MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
+#endif
// finally: compute the average deformation of the Neumann boundary
// That is what we need for the locking tests
FieldVector<double,3> averageDef(0);
+#ifdef MIXED_SPACE
+ for (size_t i=0; i<x[_0].size(); i++)
+ if (neumannNodes[i][0])
+ averageDef += x[_0][i].globalCoordinates();
+#else
for (size_t i=0; i<x.size(); i++)
if (neumannNodes[i][0])
averageDef += x[i].r;
+#endif
averageDef /= neumannNodes.count();
if (mpiHelper.rank()==0 and parameterSet.hasKey("neumannValues"))
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
deleted file mode 100644
index e1ca701e53a55eb5723f77114f2d5627733799cd..0000000000000000000000000000000000000000
--- a/src/mixed-cosserat-continuum.cc
+++ /dev/null
@@ -1,410 +0,0 @@
-#include <config.h>
-
-#include <fenv.h>
-#include <array>
-
-// 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>
-
-#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/grid/uggrid.hh>
-#include <dune/grid/onedgrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-
-#include <dune/grid/io/file/gmshreader.hh>
-
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
-#include <dune/functions/functionspacebases/compositebasis.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
-#include <dune/fufem/dunepython.hh>
-
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/norms/energynorm.hh>
-
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
-#include <dune/gfe/cosseratenergystiffness.hh>
-#include <dune/gfe/cosseratvtkwriter.hh>
-#include <dune/gfe/mixedgfeassembler.hh>
-#include <dune/gfe/mixedriemanniantrsolver.hh>
-
-// grid dimension
-const int dim = 2;
-
-using namespace Dune;
-
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
- : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
-{
- NeumannFunction(const FieldVector<double,3> values,
- double homotopyParameter)
- : values_(values),
- homotopyParameter_(homotopyParameter)
- {}
-
- void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const {
- out = 0;
- out.axpy(homotopyParameter_, values_);
- }
-
- FieldVector<double,3> values_;
- double homotopyParameter_;
-};
-
-
-int main (int argc, char *argv[]) try
-{
- // initialize MPI, finalize is done automatically on exit
- Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
-
- // Start Python interpreter
- Python::start();
- Python::Reference main = Python::import("__main__");
- Python::run("import math");
-
- //feenableexcept(FE_INVALID);
- Python::runStream()
- << std::endl << "import sys"
- << std::endl << "sys.path.append('/home/sander/dune/dune-gfe/problems/')"
- << std::endl;
-
- using namespace Dune::TypeTree::Indices;
- typedef Dune::TupleVector<std::vector<RealTuple<double,3> >,
- std::vector<Rotation<double,3> > > SolutionType;
-
- // parse data file
- ParameterTree parameterSet;
- if (argc < 2)
- DUNE_THROW(Exception, "Usage: ./mixed-cosserat-continuum <parameter file>");
-
- ParameterTreeParser::readINITree(argv[1], 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 maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
- 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");
- std::string resultPath = parameterSet.get("resultPath", "");
-
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
- typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
-
- shared_ptr<GridType> grid;
-
- FieldVector<double,dim> lower(0), upper(1);
-
- if (parameterSet.get<bool>("structuredGrid")) {
-
- lower = parameterSet.get<FieldVector<double,dim> >("lower");
- upper = parameterSet.get<FieldVector<double,dim> >("upper");
-
- 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");
- grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
- }
-
- 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::BasisBuilder;
-
- auto compositeBasis = makeBasis(
- gridView,
- composite(
- lagrange<2>(),
- lagrange<1>()
- )
- );
-
- typedef Dune::Functions::LagrangeBasis<GridView,2> DeformationFEBasis;
- typedef Dune::Functions::LagrangeBasis<GridView,1> OrientationFEBasis;
-
- DeformationFEBasis deformationFEBasis(gridView);
- OrientationFEBasis orientationFEBasis(gridView);
-
- // Construct fufem-style function space bases to ease the transition to dune-functions
- typedef DuneFunctionsBasis<DeformationFEBasis> FufemDeformationFEBasis;
- FufemDeformationFEBasis fufemDeformationFEBasis(deformationFEBasis);
-
- typedef DuneFunctionsBasis<OrientationFEBasis> FufemOrientationFEBasis;
- FufemOrientationFEBasis fufemOrientationFEBasis(orientationFEBasis);
-
-
-
- std::cout << "Deformation: " << deformationFEBasis.size() << ", orientation: " << orientationFEBasis.size() << std::endl;
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
-
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
- BitSetVector<1> neumannVertices(gridView.size(dim), false);
-
- GridType::Codim<dim>::LeafIterator vIt = gridView.begin<dim>();
- GridType::Codim<dim>::LeafIterator vEndIt = gridView.end<dim>();
-
- 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(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
-
- // Same for the Neumann boundary
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
-
- for (; vIt!=vEndIt; ++vIt) {
-
- bool isDirichlet;
- pythonDirichletVertices.evaluate(vIt->geometry().corner(0), isDirichlet);
- dirichletVertices[indexSet.index(*vIt)] = isDirichlet;
-
- bool isNeumann;
- pythonNeumannVertices.evaluate(vIt->geometry().corner(0), isNeumann);
- neumannVertices[indexSet.index(*vIt)] = isNeumann;
-
- }
-
- BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
-
- BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
- constructBoundaryDofs(neumannBoundary,fufemDeformationFEBasis,neumannNodes);
-
-
- if (mpiHelper.rank()==0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
-
-
- BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemDeformationFEBasis,deformationDirichletNodes);
-
- BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
- for (size_t i=0; i<deformationFEBasis.size(); i++)
- if (deformationDirichletNodes[i][0])
- for (int j=0; j<3; j++)
- deformationDirichletDofs[i][j] = true;
-
- BitSetVector<1> orientationDirichletNodes(orientationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemOrientationFEBasis,orientationDirichletNodes);
-
- BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
- for (size_t i=0; i<orientationFEBasis.size(); i++)
- if (orientationDirichletNodes[i][0])
- for (int j=0; j<3; j++)
- orientationDirichletDofs[i][j] = true;
-
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
-
- SolutionType x;
-
- x[_0].resize(deformationFEBasis.size());
-
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
-
- std::vector<FieldVector<double,3> > v;
- ::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
-
- for (size_t i=0; i<x[_0].size(); i++)
- x[_0][i] = v[i];
-
- x[_1].resize(orientationFEBasis.size());
-#if 0
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
-
- std::vector<FieldVector<double,3> > v;
- Functions::interpolate(feBasis, v, pythonInitialDeformation);
-
- for (size_t i=0; i<x.size(); i++)
- xDisp[i] = v[i];
-#endif
- ////////////////////////////////////////////////////////
- // Main homotopy loop
- ////////////////////////////////////////////////////////
-
- // Output initial iterate (of homotopy loop)
- CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
- orientationFEBasis,x[_1],
- resultPath + "mixed-cosserat_homotopy_0");
-
- 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");
-
- shared_ptr<NeumannFunction> neumannFunction;
- if (parameterSet.hasKey("neumannValues"))
- neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
- homotopyParameter);
-
-
- if (mpiHelper.rank() == 0) {
- std::cout << "Material parameters:" << std::endl;
- materialParameters.report();
- }
-
- // Assembler using ADOL-C
- CosseratEnergyLocalStiffness<decltype(compositeBasis),
- 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
- &neumannBoundary,
- neumannFunction,
- nullptr);
-
- MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
- RealTuple<double,3>,
- Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
-
- MixedGFEAssembler<decltype(compositeBasis),
- RealTuple<double,3>,
- Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
-
- // /////////////////////////////////////////////////
- // Create a Riemannian trust-region solver
- // /////////////////////////////////////////////////
-
- MixedRiemannianTrustRegionSolver<GridType,
- decltype(compositeBasis),
- DeformationFEBasis, RealTuple<double,3>,
- OrientationFEBasis, Rotation<double,3> > solver;
- solver.setup(*grid,
- &assembler,
- deformationFEBasis,
- orientationFEBasis,
- x,
- deformationDirichletDofs,
- orientationDirichletDofs,
- tolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance,
- instrumented);
-
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
-
- ////////////////////////////////////////////////////////
- // Set Dirichlet values
- ////////////////////////////////////////////////////////
-
- Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
-
- Python::Callable C = dirichletValuesClass.get("DirichletValues");
-
- // Call a constructor.
- Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
-
- // Extract object member functions as Dune functions
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
- PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
-
- std::vector<FieldVector<double,3> > ddV;
- ::Functions::interpolate(fufemDeformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
-
- std::vector<FieldMatrix<double,3,3> > dOV;
- ::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
-
- for (size_t j=0; j<x[_0].size(); j++)
- if (deformationDirichletNodes[j][0])
- x[_0][j] = ddV[j];
-
- for (size_t j=0; j<x[_1].size(); j++)
- if (orientationDirichletNodes[j][0])
- x[_1][j].set(dOV[j]);
-
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
-
- solver.setInitialIterate(x);
- solver.solve();
-
- x = solver.getSol();
-
- // Output result of each homotopy step
- std::stringstream iAsAscii;
- iAsAscii << i+1;
- CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
- orientationFEBasis,x[_1],
- resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
-
- }
-
- // //////////////////////////////
- // Output result
- // //////////////////////////////
-
- // 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)
- {
- std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << " "
- << ", average deflection: " << averageDef << std::endl;
- }
-
- } catch (Exception& e) {
-
- std::cout << e << std::endl;
-
- }