diff --git a/dune.module b/dune.module
index 495e6222837193efa1875164c3bf9d45229b7912..6a17dcaf439041470959d85f65e703a9df591a9d 100644
--- a/dune.module
+++ b/dune.module
@@ -8,4 +8,4 @@ Version: svn
Maintainer: oliver.sander@tu-dresden.de
#depending on
Depends: dune-common(>=2.7) dune-grid(>=2.7) dune-uggrid dune-istl dune-localfunctions dune-geometry (>=2.7) dune-functions (>=2.7) dune-solvers dune-fufem dune-elasticity (>=2.7)
-Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgeometry
+Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgrid
diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index c7012e81064f0882911a74944fad7f4aabbb1311..1ccc013a350ed54bfa8685f342a1e22e2b6fbd96 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -109,6 +109,23 @@ class CosseratVTKWriter
}
public:
+ /** \brief Write a configuration given with respect to a scalar function space basis
+ */
+ template <typename Basis>
+ static void write(const Basis& basis,
+ const Dune::TupleVector<std::vector<RealTuple<double,3> >,
+ std::vector<Rotation<double,3> > >& configuration,
+ const std::string& filename)
+ {
+ using namespace Dune::TypeTree::Indices;
+ std::vector<RigidBodyMotion<double,3>> xRBM(basis.size());
+ for (int i = 0; i < basis.size(); i++) {
+ for (int j = 0; j < 3; j ++) // Displacement part
+ xRBM[i].r[j] = configuration[_0][i][j];
+ xRBM[i].q = configuration[_1][i]; // Rotation part
+ }
+ write(basis,xRBM,filename);
+ }
/** \brief Write a configuration given with respect to a scalar function space basis
*/
diff --git a/dune/gfe/geometries/moebiusstrip.hh b/dune/gfe/geometries/moebiusstrip.hh
new file mode 100644
index 0000000000000000000000000000000000000000..793c2283a3dc1cb2cf706a3a6c43327b3a0d9531
--- /dev/null
+++ b/dune/gfe/geometries/moebiusstrip.hh
@@ -0,0 +1,106 @@
+#ifndef DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
+#define DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
+
+#include <cmath>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
+
+namespace Dune {
+
+/// \brief Functor representing a Möbius strip in 3D, where the base circle is the circle in the x-y-plane with radius_ around 0,0,0
+template <class T = double>
+class MoebiusStripProjection
+{
+ static const int dim = 3;
+ using Domain = FieldVector<T,dim>;
+ using Jacobian = FieldMatrix<T,dim,dim>;
+
+ T radius_;
+
+public:
+ MoebiusStripProjection (T radius)
+ : radius_(radius)
+ {}
+
+ /// \brief Project the coordinate to the MS using closest point projection
+ Domain operator() (const Domain& x) const
+ {
+ double nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+ //center for the point x - it lies on the circle around (0,0,0) with radius radius_ in the x-y-plane
+ Domain center{x[0] * radius_/nrm, x[1] * radius_/nrm, 0};
+ double cosu = x[0]/nrm;
+ double sinu = x[1]/nrm;
+ double cosuhalf = std::sqrt((1+cosu)/2);
+ cosuhalf *= (sinu < 0) ? (-1) : 1 ; // u goes from 0 to 2*pi, so cosuhalf is negative if sinu < 0, we multiply that here
+ double sinuhalf = std::sqrt((1-cosu)/2); //u goes from 0 to 2*pi, so sinuhalf is always >= 0
+
+ // now calculate line from center to the new point
+ // the direction is (cosuhalf*cosu,cosuhalf*sinu,sinuhalf), as can be seen from the formula to construct the MS, we normalize this vector
+ Domain centerToPointOnMS{cosuhalf*cosu,cosuhalf*sinu,sinuhalf};
+ centerToPointOnMS /= centerToPointOnMS.two_norm();
+
+ Domain centerToX = center - x;
+ // We need the length, let theta be the angle between the vector centerToX and center to centerToPointOnMS
+ // then cos(theta) = centerToX*centerToPointOnMS/(len(centerToX)*len(centerToPointOnMS))
+ // We want to project the point to MS, s.t. the angle between the MS and the projection is 90deg
+ // Then, length = cos(theta) * len(centerToX) = centerToX*centerToPointOnMS/len(centerToPointOnMS) = centerToX*centerToPointOnMS
+ double length = - centerToX * centerToPointOnMS;
+ centerToPointOnMS *= length;
+ return center + centerToPointOnMS;
+ }
+
+
+ /// \brief derivative of the projection to the MS
+ friend auto derivative (const MoebiusStripProjection& moebiusStrip)
+ {
+ DUNE_THROW(NotImplemented,"The derivative of the projection to the Möbius strip is not implemented yet!");
+ return [radius = moebiusStrip.radius_](const Domain& x)
+ {
+ Jacobian out;
+ return out;
+ };
+ }
+
+ /// \brief Normal Vector of the MS
+ Domain normal (const Domain& x) const
+ {
+ using std::sqrt;
+ Domain nVec = {x[0],x[1],0};
+ double nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+ double cosu = x[0]/nrm;
+ double sinu = x[1]/nrm;
+ double cosuhalf = std::sqrt((1+cosu)/2);
+ cosuhalf *= (sinu < 0) ? (-1) : 1 ; // u goes from 0 to 2*pi, so cosuhalf is negative if sinu < 0, we multiply that here
+ double sinuhalf = std::sqrt((1-cosu)/2);
+ nVec[2] = (-1)*cosuhalf*(cosu*x[0]+sinu*x[1])/sinuhalf;
+ nVec /= nVec.two_norm();
+ return nVec;
+ }
+
+ /// \brief The mean curvature of the MS
+ T mean_curvature (const Domain& /*x*/) const
+ {
+ DUNE_THROW(NotImplemented,"The mean curvature of the Möbius strip is not implemented yet!");
+ return 1;
+ }
+
+ /// \brief The area of the MS
+ T area () const
+ {
+ DUNE_THROW(NotImplemented,"The area of the Möbius strip is not implemented yet!");
+ return 1;
+ }
+};
+
+/// \brief construct a grid function representing the parametrization of a MS
+template <class Grid, class T>
+auto moebiusStripGridFunction (T radius)
+{
+ return analyticGridFunction<Grid>(MoebiusStripProjection<T>{radius});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
diff --git a/dune/gfe/geometries/twistedstrip.hh b/dune/gfe/geometries/twistedstrip.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d7af48436d18dc99f9c8b787d8f94865448a5c6b
--- /dev/null
+++ b/dune/gfe/geometries/twistedstrip.hh
@@ -0,0 +1,104 @@
+#ifndef DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
+#define DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
+
+#include <cmath>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
+
+namespace Dune {
+
+/// \brief Functor representing a twisted strip in 3D, with length length_ and nTwists twists
+template <class T = double>
+class TwistedStripProjection
+{
+ static const int dim = 3;
+ using Domain = FieldVector<T,dim>;
+ using Jacobian = FieldMatrix<T,dim,dim>;
+
+ T length_;
+ double nTwists_;
+
+public:
+ TwistedStripProjection (T length, double nTwists)
+ : length_(length),
+ nTwists_(nTwists)
+ {}
+
+ /// \brief Project the coordinate to the twisted strip using closest point projection
+ Domain operator() (const Domain& x) const
+ {
+ // Angle for the x - position of the point
+ double alpha = std::acos(-1)*nTwists_*x[0]/length_;
+ double cosalpha = std::cos(alpha);
+ double sinalpha = std::sin(alpha);
+
+ // Add the line from middle to the new point - the direction is (0,cosalpha,sinalpha)
+ // We need the distance, calculated by (y^2 + z^2)^0.5
+ double dist = std::sqrt(x[1]*x[1] + x[2]*x[2]);
+ double y = std::abs(cosalpha*dist);
+ if (x[1] < 0 and std::abs(x[1]) > std::abs(x[2])) { // only trust the sign of x[1] if it is larger than x[2]
+ y *= (-1);
+ } else if (std::abs(x[1]) <= std::abs(x[2])) {
+ if (cosalpha * sinalpha >= 0 and x[2] < 0) // if cosalpha*sinalpha > 0, x[1] and x[2] must have the same sign
+ y *= (-1);
+ if (cosalpha * sinalpha <= 0 and x[2] > 0) // if cosalpha*sinalpha < 0, x[1] and x[2] must have opposite signs
+ y *= (-1);
+ }
+ double z = std::abs(sinalpha*dist);
+ if (x[2] < 0 and std::abs(x[2]) > std::abs(x[1])) {
+ z *= (-1);
+ } else if (std::abs(x[2]) <= std::abs(x[1])) {
+ if (cosalpha * sinalpha >= 0 and x[1] < 0)
+ z *= (-1);
+ if (cosalpha * sinalpha <= 0 and x[1] > 0)
+ z *= (-1);
+ }
+ return {x[0], y , z};
+ }
+
+ /// \brief derivative of the projection to the twistedstrip
+ friend auto derivative (const TwistedStripProjection& twistedStrip)
+ {
+ DUNE_THROW(NotImplemented,"The derivative of the projection to the twisted strip is not implemented yet!");
+ return [length = twistedStrip.length_, nTwists = twistedStrip.nTwists_](const Domain& x)
+ {
+ Jacobian out;
+ return out;
+ };
+ }
+
+ /// \brief Normal Vector of the twisted strip
+ Domain normal (const Domain& x) const
+ {
+ // Angle for the x - position of the point
+ double alpha = std::acos(-1)*nTwists_*x[0]/length_;
+ std::cout << "normal at " << x[0] << "is " << -std::sin(alpha) << ", " << std::cos(alpha) << std::endl;
+ return {0,-std::sin(alpha), std::cos(alpha)};
+ }
+
+ /// \brief The mean curvature of the twisted strip
+ T mean_curvature (const Domain& /*x*/) const
+ {
+ DUNE_THROW(NotImplemented,"The mean curvature of the twisted strip is not implemented yet!");
+ return 1;
+ }
+
+ /// \brief The area of the twisted strip
+ T area (T width) const
+ {
+ return length_*width;
+ }
+};
+
+/// \brief construct a grid function representing the parametrization of a twisted strip
+template <class Grid, class T>
+auto twistedStripGridFunction (T length, double nTwists)
+{
+ return analyticGridFunction<Grid>(TwistedStripProjection<T>{length, nTwists});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index 8fd0aceb4de0effe76dc7d854a990bf5cd7396d2..21a1f2a9e390bff5db91d6984bfcce7b03c56343 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -9,6 +9,8 @@
#include <dune/fufem/boundarypatch.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
#include <dune/gfe/localenergy.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/rigidbodymotion.hh>
@@ -21,7 +23,15 @@
#include <dune/localfunctions/lagrange/lfecache.hh>
#endif
-template<class Basis, int dim, class field_type=double>
+/** \brief Assembles the cosserat energy for a single element.
+ *
+ * \tparam Basis Type of the Basis used for assembling
+ * \tparam dim Dimension of the Targetspace, 3
+ * \tparam field_type The coordinate type of the TargetSpace
+ * \tparam StressFreeStateGridFunction Type of the GridFunction representing the Cosserat shell in a stress free state
+ */
+template<class Basis, int dim, class field_type=double, class StressFreeStateGridFunction =
+ Dune::Functions::DiscreteGlobalBasisFunction<Basis,std::vector<Dune::FieldVector<double, Basis::GridView::dimensionworld>> > >
class NonplanarCosseratShellEnergy
: public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >
{
@@ -40,13 +50,16 @@ class NonplanarCosseratShellEnergy
public:
/** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
+ * \param parameters The material parameters
+ * \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
*/
NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
+ const StressFreeStateGridFunction* stressFreeStateGridFunction,
const BoundaryPatch<GridView>* neumannBoundary,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
- : neumannBoundary_(neumannBoundary),
+ : stressFreeStateGridFunction_(stressFreeStateGridFunction),
+ neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction),
volumeLoad_(volumeLoad)
{
@@ -111,6 +124,9 @@ public:
/** \brief Curvature parameters */
double b1_, b2_, b3_;
+ /** \brief The geometry used for assembling */
+ const StressFreeStateGridFunction* stressFreeStateGridFunction_;
+
/** \brief The Neumann boundary */
const BoundaryPatch<GridView>* neumannBoundary_;
@@ -121,15 +137,34 @@ public:
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
};
-template <class Basis, int dim, class field_type>
-typename NonplanarCosseratShellEnergy<Basis,dim,field_type>::RT
-NonplanarCosseratShellEnergy<Basis,dim,field_type>::
+template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
+typename NonplanarCosseratShellEnergy<Basis, dim, field_type, StressFreeStateGridFunction>::RT
+NonplanarCosseratShellEnergy<Basis,dim,field_type, StressFreeStateGridFunction>::
energy(const typename Basis::LocalView& localView,
const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
{
// The element geometry
auto element = localView.element();
+
+#if HAVE_DUNE_CURVEDGEOMETRY
+ // Construct a curved geometry of this element of the Cosserat shell in stress-free state
+ // When using element.geometry(), then the curvatures on the element are zero, when using a curved geometry, they are not
+ // If a parametrization representing the Cosserat shell in a stress-free state is given,
+ // this is used for the curved geometry approximation.
+ // The variable local holds the local coordinates in the reference element
+ // and localGeometry.global maps them to the world coordinates
+ Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> geometry(referenceElement(element),
+ [this,element](const auto& local) {
+ if (not stressFreeStateGridFunction_) {
+ return element.geometry().global(local);
+ }
+ auto localGridFunction = localFunction(*stressFreeStateGridFunction_);
+ localGridFunction.bind(element);
+ return localGridFunction(local);
+ }, 2); /*order*/
+#else
auto geometry = element.geometry();
+#endif
// The set of shape functions on this element
const auto& localFiniteElement = localView.tree().finiteElement();
@@ -215,17 +250,8 @@ energy(const typename Basis::LocalView& localView,
c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
#if HAVE_DUNE_CURVEDGEOMETRY
- // Construct a curved geometry to evaluate the derivative of the normal field on each quadrature point
- // The variable local holds the local coordinates in the reference element
- // and localGeometry.global maps them to the world coordinates
- // we want to take the derivative of the normal field on the element in world coordinates
- Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> curvedGeometry(referenceElement(element),
- [localGeometry=element.geometry()](const auto& local) {
- return localGeometry.global(local);
- }, 1); //order = 1
-
- // Second fundamental form: The derivative of the normal field
- auto normalDerivative = curvedGeometry.normalGradient(quad[pt].position());
+ // Second fundamental form: The derivative of the normal field, on each quadrature point
+ auto normalDerivative = geometry.normalGradient(quad[pt].position());
#else
//In case dune-curvedgeometry is not installed, the normal derivative is set to zero.
Dune::FieldMatrix<double,3,3> normalDerivative(0);
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 8e0ca8bc5dd8478f53a1123463ae8353c2619d1b..bb7dc9c666d9d94e14b4a825ec0c99887dab3826 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -67,7 +67,7 @@ setup(const GridType& grid,
// //////////////////////////////////////////////////////////////////////////////////////
typedef DuneFunctionsBasis<Basis> FufemBasis;
- FufemBasis basis(assembler_->basis_);
+ FufemBasis basis(assembler_->getBasis());
OperatorAssembler<FufemBasis,FufemBasis> operatorAssembler(basis, basis);
LaplaceAssembler<GridType, typename FufemBasis::LocalFiniteElement, typename FufemBasis::LocalFiniteElement> laplaceStiffness;
diff --git a/problems/cosserat-continuum-cantilever.parset b/problems/cosserat-continuum-cantilever.parset
index a697675dde85aed1bb666aa657d391ded7affca2..9440301af2c2e3a244b9799310e33768d16e1212 100644
--- a/problems/cosserat-continuum-cantilever.parset
+++ b/problems/cosserat-continuum-cantilever.parset
@@ -24,7 +24,7 @@ numHomotopySteps = 1
tolerance = 1e-3
# Max number of steps of the trust region solver
-maxTrustRegionSteps = 200
+maxSolverSteps = 200
trustRegionScaling = 1 1 1 0.01 0.01 0.01
diff --git a/problems/cosserat-continuum-wriggers-l-shape.parset b/problems/cosserat-continuum-wriggers-l-shape.parset
index eb108bf0f5cbdf7474cd46f1234f78c263009278..12c715bd9e9618fa4ce069b8409584294f0f5660 100644
--- a/problems/cosserat-continuum-wriggers-l-shape.parset
+++ b/problems/cosserat-continuum-wriggers-l-shape.parset
@@ -20,7 +20,7 @@ numHomotopySteps = 1
tolerance = 1e-8
# Max number of steps of the trust region solver
-maxTrustRegionSteps = 1000
+maxSolverSteps = 1000
trustRegionScaling = 1 1 1 0.01 0.01 0.01
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index 55b999fca1bb4703c13bd3126015711877869661..8ccce0d1c12ba5b5145e813910dade3ba1281cea 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -59,7 +59,7 @@ initialDeformation = "x"
#############################################
# Inner solver, cholmod or multigrid
-solvertype = multigrid
+solvertype = trustRegion
# Number of homotopy steps for the Dirichlet boundary conditions
numHomotopySteps = 1
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index fb68dbcd00eb541d27049b8843c6eafc6612b472..19fcf4135b2debf46457c937552555e858481b33 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -1,3 +1,5 @@
+#define MIXED_SPACE 0
+
#include <config.h>
#include <fenv.h>
@@ -37,7 +39,6 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
-#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
#include <dune/gfe/cosseratenergystiffness.hh>
@@ -45,10 +46,17 @@
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/cosseratvtkreader.hh>
#include <dune/gfe/geodesicfeassembler.hh>
-#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/mixedgfeassembler.hh>
+
+#if MIXED_SPACE
#include <dune/gfe/mixedriemanniantrsolver.hh>
+#else
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
+#include <dune/gfe/riemannianpnsolver.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+#endif
#if HAVE_DUNE_VTK
#include <dune/vtk/vtkreader.hh>
@@ -59,22 +67,17 @@
const int dim = 2;
const int dimworld = 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
+#if !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;
+using TargetSpace = RigidBodyMotion<double, 3>;
#endif
using namespace Dune;
@@ -84,6 +87,15 @@ 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;
+ #else
+ std::cout << "MIXED_SPACE = 0" << std::endl;
+ #endif
+ }
// Start Python interpreter
Python::start();
Python::Reference main = Python::import("__main__");
@@ -97,12 +109,8 @@ int main (int argc, char *argv[]) try
<< 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;
@@ -117,7 +125,8 @@ int main (int argc, char *argv[]) try
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 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");
@@ -145,6 +154,8 @@ int main (int argc, char *argv[]) try
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");
@@ -183,7 +194,7 @@ int main (int argc, char *argv[]) try
GridView gridView = grid->leafGridView();
using namespace Dune::Functions::BasisFactory;
-#ifdef MIXED_SPACE
+
const int dimRotation = Rotation<double,dim>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
@@ -195,6 +206,13 @@ int main (int argc, char *argv[]) try
lagrange<rotationOrder>()
)
));
+ using CompositeBasis = decltype(compositeBasis);
+
+ auto deformationPowerBasis = makeBasis(
+ gridView,
+ power<3>(
+ lagrange<displacementOrder>()
+ ));
typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
@@ -202,10 +220,6 @@ int main (int argc, char *argv[]) try
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
-#else
- typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, displacementOrder> FEBasis;
- FEBasis feBasis(gridView);
-#endif
// /////////////////////////////////////////
// Read Dirichlet values
@@ -237,10 +251,10 @@ int main (int argc, char *argv[]) try
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
- if (mpiHelper.rank()==0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
-
-#ifdef MIXED_SPACE
+ std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces()
+ << " faces and " << dirichletVertices.count() << " nodes.\n";
+ std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces() << " faces.\n";
+
BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
constructBoundaryDofs(dirichletBoundary,deformationFEBasis,deformationDirichletNodes);
@@ -261,43 +275,26 @@ int main (int argc, char *argv[]) try
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);
-
- BitSetVector<1> neumannNodes(feBasis.size(), false);
- constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
-
-
- BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
- for (size_t i=0; i<feBasis.size(); i++)
- 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());
+ 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(deformationFEBasis, v, pythonInitialDeformation);
+ Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
for (size_t i=0; i<x[_0].size(); i++)
x[_0][i] = v[i];
- x[_1].resize(orientationFEBasis.size());
-#else
- SolutionType x(feBasis.size());
-
+ x[_1].resize(compositeBasis.size({1}));
+#if !MIXED_SPACE
if (parameterSet.hasKey("startFromFile"))
{
std::shared_ptr<GridType> initialIterateGrid;
@@ -314,35 +311,32 @@ int main (int argc, char *argv[]) try
initialIterateGrid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + initialIterateGridFilename));
}
- SolutionType initialIterate;
+ 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 = LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+ using LocalInterpolationRule = LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
#else
- using LocalInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::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;
- 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]);
-
- } else {
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
+ //TODO: Interpolate does not work with an GFE:EmbeddedGlobalGFEFunction
+ //Dune::Functions::interpolate(powerBasis,v,initialFunction);
- std::vector<FieldVector<double,3> > v;
- Functions::interpolate(feBasis, v, pythonInitialDeformation);
-
- for (size_t i=0; i<x.size(); i++)
- x[i].r = v[i];
+ for (size_t i=0; i<x.size(); i++) {
+ auto vTargetSpace = TargetSpace(v[i]);
+ x[_0][i] = vTargetSpace.r;
+ x[_1][i] = vTargetSpace.q;
+ }
}
#endif
@@ -351,14 +345,30 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
-#ifdef MIXED_SPACE
+ BlockVector<FieldVector<double,3> > identity(compositeBasis.size({0}));
+ Dune::Functions::interpolate(deformationPowerBasis, identity, [&](FieldVector<double,dimworld> x){ return x;});
+ BlockVector<FieldVector<double,3> > displacement(compositeBasis.size({0}));
+
+ if (dim == dimworld) {
CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
orientationFEBasis,x[_1],
resultPath + "mixed-cosserat_homotopy_0");
+ } else {
+#if MIXED_SPACE
+ for (int i = 0; i < displacement.size(); i++) {
+ for (int j = 0; j < 3; j++)
+ displacement[i][j] = x[_0][i][j];
+ displacement[i] -= identity[i];
+ }
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationPowerBasis, displacement);
+
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dimworld));
+ vtkWriter.write(resultPath + "mixed-cosserat_homotopy_0");
#else
- CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_0");
-#endif
-
+ CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_0");
+#endif
+ }
for (int i=0; i<numHomotopySteps; i++) {
double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
@@ -375,7 +385,7 @@ int main (int argc, char *argv[]) try
if (parameterSet.hasKey("neumannValues"))
neumannValues = parameterSet.get<FieldVector<double,3> >("neumannValues");
- auto neumannFunction = [&]( FieldVector<double,dim> ) {
+ auto neumannFunction = [&]( FieldVector<double,dimworld> ) {
auto nV = neumannValues;
nV *= homotopyParameter;
return nV;
@@ -385,7 +395,7 @@ int main (int argc, char *argv[]) try
if (parameterSet.hasKey("volumeLoad"))
volumeLoadValues = parameterSet.get<FieldVector<double,3> >("volumeLoad");
- auto volumeLoad = [&]( FieldVector<double,dim>) {
+ auto volumeLoad = [&]( FieldVector<double,dimworld>) {
auto vL = volumeLoadValues;
vL *= homotopyParameter;
return vL;
@@ -396,82 +406,6 @@ int main (int argc, char *argv[]) try
materialParameters.report();
}
- // Assembler using ADOL-C
-#ifdef MIXED_SPACE
- CosseratEnergyLocalStiffness<decltype(compositeBasis),
- 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
- &neumannBoundary,
- neumannFunction,
- volumeLoad);
-
- MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
- RealTuple<double,3>,
- Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
-
- MixedGFEAssembler<decltype(compositeBasis),
- RealTuple<double,3>,
- Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
-#else
- std::shared_ptr<GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble,3> > > localCosseratEnergy;
-
- if (dim==dimworld)
- {
- localCosseratEnergy = std::make_shared<CosseratEnergyLocalStiffness<FEBasis,3,adouble> >(materialParameters,
- &neumannBoundary,
- neumannFunction,
- volumeLoad);
- }
- else
- {
- localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
- &neumannBoundary,
- neumannFunction,
- volumeLoad);
- }
-
- LocalGeodesicFEADOLCStiffness<FEBasis,
- TargetSpace> localGFEADOLCStiffness(localCosseratEnergy.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,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance,
- instrumented);
-
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
-
////////////////////////////////////////////////////////
// Set Dirichlet values
////////////////////////////////////////////////////////
@@ -484,85 +418,227 @@ int main (int argc, char *argv[]) try
Python::Reference dirichletValuesPythonObject = C(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"));
-
- std::vector<FieldVector<double,3> > ddV;
- std::vector<FieldMatrix<double,3,3> > dOV;
-
-#ifdef MIXED_SPACE
- Functions::interpolate(deformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
- Functions::interpolate(orientationFEBasis, 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]);
+ 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, deformationDirichletDofs);
+
+ BlockVector<FieldMatrix<double,3,3> > dOV;
+ Dune::Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
+
+ for (int i = 0; i < compositeBasis.size({0}); i++) {
+ if (deformationDirichletDofs[i][0])
+ x[_0][i] = ddV[i];
+ }
+ for (int i = 0; i < compositeBasis.size({1}); i++)
+ if (orientationDirichletDofs[i][0])
+ x[_1][i].set(dOV[i]);
+
+ if (dim==dimworld) {
+ CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> localCosseratEnergy(materialParameters,
+ &neumannBoundary,
+ neumannFunction,
+ volumeLoad);
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,3>,
+ Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,3>,
+ Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
+#if MIXED_SPACE
+ 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> >("trustRegionScaling"));
+
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
#else
- Functions::interpolate(feBasis, ddV, deformationDirichletValues, dirichletDofs);
- Functions::interpolate(feBasis, dOV, orientationDirichletValues, dirichletDofs);
-
- for (size_t j=0; j<x.size(); j++)
- if (dirichletNodes[j][0])
- {
- x[j].r = ddV[j];
- x[j].q.set(dOV[j]);
- }
+ //The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
+ //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
+ //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion 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 (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < 3; j ++) { // Displacement part
+ xTargetSpace[i].r[j] = x[_0][i][j];
+ dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
+ }
+ xTargetSpace[i].q = 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, RealTuple<double, 3>, Rotation<double,3>>;
+ GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+ RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
+ solver.setup(*grid,
+ &assembler,
+ xTargetSpace,
+ dirichletDofsTargetSpace,
+ tolerance,
+ maxSolverSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setInitialIterate(xTargetSpace);
+ solver.solve();
+ xTargetSpace = solver.getSol();
+ for (int i = 0; i < xTargetSpace.size(); i++) {
+ x[_0][i] = xTargetSpace[i].r;
+ x[_1][i] = xTargetSpace[i].q;
+ }
#endif
+ } else {
+#if MIXED_SPACE
+ DUNE_THROW(Exception, "Problems with dim != dimworld do not work for MIXED_SPACE = 1!");
+#else
+ std::shared_ptr<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>> localGFEADOLCStiffness;
+
+ NonplanarCosseratShellEnergy<DeformationFEBasis, 3, adouble> localCosseratEnergyPlanar(materialParameters,
+ nullptr,
+ &neumannBoundary,
+ neumannFunction,
+ volumeLoad);
+ localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>>(&localCosseratEnergyPlanar);
+
+ GeodesicFEAssembler<DeformationFEBasis,TargetSpace> assembler(gridView, *localGFEADOLCStiffness);
+ //The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
+ //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
+ //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion 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 (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < 3; j ++) { // Displacement part
+ xTargetSpace[i].r[j] = x[_0][i][j];
+ dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
+ }
+ xTargetSpace[i].q = x[_1][i]; // Rotation part
+ for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
+ dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
+ }
+ if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
+
+ RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace> solver;
+ solver.setup(*grid,
+ &assembler,
+ xTargetSpace,
+ dirichletDofsTargetSpace,
+ tolerance,
+ maxSolverSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setInitialIterate(xTargetSpace);
+ solver.solve();
+ xTargetSpace = solver.getSol();
+ } else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
+#if DUNE_VERSION_LT(DUNE_COMMON, 2, 8)
+ DUNE_THROW(Exception, "Please install dune-solvers >= 2.8 to use the Proximal Newton Solver with Cholmod!");
+#else
+ RiemannianProximalNewtonSolver<DeformationFEBasis, TargetSpace> solver;
+ solver.setup(*grid,
+ &assembler,
+ xTargetSpace,
+ dirichletDofsTargetSpace,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ solver.setInitialIterate(xTargetSpace);
+ solver.solve();
+ xTargetSpace = solver.getSol();
+#endif
+ }
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
-
- solver.setInitialIterate(x);
- solver.solve();
-
- x = solver.getSol();
-
+ for (int i = 0; i < xTargetSpace.size(); i++) {
+ x[_0][i] = xTargetSpace[i].r;
+ x[_1][i] = xTargetSpace[i].q;
+ }
+#endif
+ }
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
-#ifdef MIXED_SPACE
+ if (dim == dimworld) {
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());
+ } else {
+#if MIXED_SPACE
+ for (int i = 0; i < displacement.size(); i++) {
+ for (int j = 0; j < 3; j++) {
+ displacement[i][j] = x[_0][i][j];
+ }
+ displacement[i] -= identity[i];
+ }
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationPowerBasis, displacement);
+
+ // We need to subsample, because VTK cannot natively display real second-order functions
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+ vtkWriter.write(resultPath + "cosserat_homotopy_" + std::to_string(i+1));
+#else
+ CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_" + std::to_string(i+1));
#endif
-
+ }
}
// //////////////////////////////
// Output result
// //////////////////////////////
-#ifndef MIXED_SPACE
+#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(x.size());
- for (size_t i=0; i<x.size(); i++)
- xEmbedded[i] = x[i].globalCoordinates();
+ 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();
#endif
+ 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);
-#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"))
@@ -570,6 +646,7 @@ int main (int argc, char *argv[]) try
std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << " "
<< ", average deflection: " << averageDef << std::endl;
}
+ }
} catch (Exception& e)
{
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 5f5134fe7d641f802b120830cb58d9892fd5f731..e52201554d6630a4c50dde4f5283c60bd5d9d592 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -554,7 +554,7 @@ int main (int argc, char *argv[]) try
// Create the chosen solver and solve!
///////////////////////////////////////////////////
- if (parameterSet.get<std::string>("solvertype") == "multigrid") {
+ if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion") {
#if MIXED_SPACE
MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>> solver;
solver.setup(*grid,
@@ -603,7 +603,7 @@ int main (int argc, char *argv[]) try
x[_1][i] = xRBM[i].q;
}
#endif
- } else { //parameterSet.get<std::string>("solvertype") == "cholmod"
+ } else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
#if MIXED_SPACE
DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");