From 08e637de44ff49435fe58d317cc67c0d82409d5b Mon Sep 17 00:00:00 2001
From: AlexanderMueller <rathet@gmail.com>
Date: Tue, 30 Nov 2021 19:53:34 +0100
Subject: [PATCH] Add Simo-Fox shell model
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
This commit adds the energy implementat for the nonlinear Reissner-Mindlin
shell model known as the Simo and Fox shell model.
Details may be found in
"A consistent finite element formulation of the geometrically
non-linear Reissner-Mindlin shell model [Müller, Bischoff]"
---
dune/gfe/simofoxenergy.hh | 358 ++++++++++++++++
.../simofox-cantilever-dirichlet-values.py | 16 +
problems/simofox-cantilever.parset | 99 +++++
src/CMakeLists.txt | 3 +-
src/simofoxshell.cc | 390 ++++++++++++++++++
5 files changed, 865 insertions(+), 1 deletion(-)
create mode 100644 dune/gfe/simofoxenergy.hh
create mode 100644 problems/simofox-cantilever-dirichlet-values.py
create mode 100644 problems/simofox-cantilever.parset
create mode 100644 src/simofoxshell.cc
diff --git a/dune/gfe/simofoxenergy.hh b/dune/gfe/simofoxenergy.hh
new file mode 100644
index 00000000..853a5fe6
--- /dev/null
+++ b/dune/gfe/simofoxenergy.hh
@@ -0,0 +1,358 @@
+#ifndef DUNE_GFE_SIMOFOX_ENERGY_HH
+#define DUNE_GFE_SIMOFOX_ENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parametertree.hh>
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 8)
+#include <dune/common/transpose.hh>
+#endif
+#include <dune/common/tuplevector.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/gfe/linearalgebra.hh>
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
+#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/unitvector.hh>
+
+namespace Dune::GFE {
+
+ /** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
+ *
+ * Copied from CosseratEnergyLocalStiffness.hh
+ */
+ template <class Basis, std::size_t i>
+ class LocalFiniteElementFactory {
+ public:
+ static auto get(const typename Basis::LocalView &localView, std::integral_constant<std::size_t, i> iType)
+ -> decltype(localView.tree().child(iType, 0).finiteElement()) {
+ return localView.tree().child(iType, 0).finiteElement();
+ }
+ };
+
+ /** \brief Specialize for scalar bases, here we cannot call tree().child() */
+ template <class GridView, int order, std::size_t i>
+ class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView, order>, i> {
+ public:
+ static auto get(const typename Dune::Functions::LagrangeBasis<GridView, order>::LocalView &localView,
+ std::integral_constant<std::size_t, i> iType) -> decltype(localView.tree().finiteElement()) {
+ return localView.tree().finiteElement();
+ }
+ };
+
+
+ /** \brief Implements the energy of a Simo-Fox shell
+ *
+ * Described in:
+ * "A consistent finite element formulation of the geometrically non-linear Reissner-Mindlin shell model [Müller,
+ * Bischoff]"
+ *
+ * \tparam LocalFEFunction Decides which interpolation function is used for the interpolation of the midsurface
+ * and director field.
+ */
+ template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type = double>
+ class SimoFoxEnergyLocalStiffness
+ : public Dune::GFE::LocalEnergy<Basis, RealTuple<field_type, 3>,
+ UnitVector<field_type, 3>>, // inheritance to allow usage with LocalGeodesicFEADOLCStiffness
+ public MixedLocalGeodesicFEStiffness<Basis, RealTuple<field_type, 3>,
+ UnitVector<field_type, 3>> // inheritance to allow usage with MixedGFEAssembler
+ {
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype DT;
+ typedef field_type RT;
+ typedef typename GridView::template Codim<0>::Entity Entity;
+
+ // some other sizes
+ static constexpr int gridDim = GridView::dimension;
+ static constexpr int dimworld = GridView::dimensionworld;
+
+ // The local finite element type used for midsurface position and displacement interpolation
+ using MidSurfaceElement =
+ typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 0>::get)(typename Basis::LocalView, decltype(Dune::Indices::_0))>::type;
+ // The local finite element type used for the director interpolation
+ using DirectorElement =
+ typename std::result_of<decltype (&LocalFiniteElementFactory<Basis, 1>::get)(typename Basis::LocalView, decltype(Dune::Indices::_1))>::type;
+
+ // The local finite element function type to evaluate the midsurface position
+ using LocalMidSurfaceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<field_type, 3>>;
+ // The local finite element function type to evaluate the director
+ using LocalDirectorFunctionType = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<field_type, 3>>;
+
+ // Extra function type for reference quantities since they are unconditionally doubles, i.e. no ADOL-C types
+ using LocalMidSurfaceReferenceFunctionType = LocalFEFunction<gridDim, DT, MidSurfaceElement, RealTuple<double, 3>>;
+ using LocalDirectorReferenceFunctionType = LocalFEFunction<gridDim, DT, DirectorElement, UnitVector<double, 3>>;
+
+ public:
+ /** \brief Constructor with a set of material parameters
+ * \param parameters The material parameters
+ * \param x0 reference configuration
+ */
+ SimoFoxEnergyLocalStiffness(const Dune::ParameterTree ¶meters, const BoundaryPatch<GridView> *neumannBoundary,
+ const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> neumannFunction,
+ const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> volumeLoad,
+ const Dune::TupleVector<std::vector<RealTuple<double, 3>>, std::vector<UnitVector<double, 3>>> &x0)
+ : neumannBoundary_(neumannBoundary),
+ neumannFunction_(neumannFunction),
+ thickness_{parameters.template get<double>("thickness")}, // The sheeqll thickness
+ mu_{parameters.template get<double>("mu")}, // Lame constant 1
+ lambda_{parameters.template get<double>("lambda")}, // Lame constant 2
+ kappa_{parameters.template get<double>("kappa")}, // Shear correction factor
+ midSurfaceRefConfig{x0[Dune::Indices::_0]},
+ directorRefConfig{x0[Dune::Indices::_1]}
+ {
+ // Calculate the over the thickness preintegrated St.Venant Kirchhoff material matrix, Paper equation 10.1 */
+ const double Emodul = mu_ * (3 * lambda_ + 2 * mu_) / (lambda_ + mu_); // Young's modulus
+ const double nu = lambda_ / (2 * (lambda_ + mu_)); // Poisson ratio
+
+ // membrane
+ const double fac1 = thickness_ * Emodul / (1 - nu * nu);
+ CMat_[0][0] = CMat_[1][1] = fac1;
+ CMat_[2][2] = fac1 * (1 - nu) * 0.5;
+ CMat_[1][0] = CMat_[0][1] = fac1 * nu;
+
+ // bending
+ const double fac2 = thickness_ * thickness_ * thickness_ / 12 * Emodul / (1 - nu * nu);
+ CMat_[3][3] = CMat_[4][4] = fac2;
+ CMat_[5][5] = fac2 * (1 - nu) * 0.5;
+ CMat_[3][4] = CMat_[4][3] = fac2 * nu;
+
+ // transverse shear
+ const double fac3 = kappa_ * thickness_ * Emodul * 0.5 / (1 + nu);
+ CMat_[6][6] = CMat_[7][7] = fac3;
+ }
+
+ /** \brief Assemble the energy for a single element */
+ RT energy(const typename Basis::LocalView &localView, const std::vector<RealTuple<field_type, 3>> &localMidSurfaceConfiguration,
+ const std::vector<UnitVector<field_type, 3>> &localDirectorConfiguration) const override;
+
+ private:
+ /** \brief A structure that contains all quantities to calculate the Lagrangian strains */
+ struct KinematicVariables {
+ // current configuration
+ FieldMatrix<field_type, 2, 3> a1anda2; // first and second tangent base vector on the current geometry
+ FieldMatrix<field_type, 2, 3> td1Andtd2; // partial derivative of the director on the current geometry in first and second direction
+ FieldMatrix<field_type, 2, 3> ud1andud2; // partial derivative of the displacement function in first and second direction
+ FieldVector<field_type, 3> t; // unit director on the current geometry
+
+ // reference configuration
+ FieldMatrix<double, 2, 3> A1andA2; // first and second tangent base vector on the reference geometry
+ FieldMatrix<double, 2, 3> t0d1Andt0d2; // partial derivative of the director on the reference geometry in first and second direction
+ FieldVector<double, 3> t0; // unit director on the reference geometry
+ };
+
+ auto getReferenceLocalConfigurations(const typename Basis::LocalView &localView) const;
+
+ /** \brief Calculates all kinematic quantities */
+ template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
+ typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
+ static auto kinematicVariablesFactory(const Element &element, const LocalDirectorFunction &directorFunction,
+ const LocalDirectorReferenceFunction &directorReferenceFunction,
+ const LocalMidSurfaceFunction &midSurfaceFunction,
+ const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
+ const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos);
+
+ /** \brief Calculate the Green-Lagrange strain components */
+ static Dune::FieldVector<RT, 8> calculateGreenLagrangianStrains(const KinematicVariables &kin);
+
+ /** \brief Save all tangent base matrices for all nodes in one place*/
+ Dune::BlockVector<Dune::FieldMatrix<field_type, 2, 3>> directorTangentSpaces;
+
+ /** \brief The Neumann boundary */
+ const BoundaryPatch<GridView> *neumannBoundary_;
+
+ /** \brief The function implementing the Neumann data */
+ const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> neumannFunction_;
+
+ /** \brief The shell thickness */
+ double thickness_;
+
+ /** \brief Lame constants */
+ double mu_, lambda_;
+
+ /** \brief Shear correction factor */
+ double kappa_;
+
+ /** \brief Material tangent matrix */
+ Dune::FieldMatrix<double, 8, 8> CMat_;
+
+ /** \brief The function implementing a volume load */
+ const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> volumeLoad_;
+
+ /** \brief Stores the reference configuration of the midsurface and the director field */
+ const std::vector<RealTuple<double, 3>> &midSurfaceRefConfig;
+ const std::vector<UnitVector<double, 3>> &directorRefConfig;
+ };
+
+ /** \brief Calculate the Green-Lagrange strain components for the thickness-integrated setting
+ *
+ * The components are returned as [a_11, a_22, a_12, b_11, b_22, b_12, gamma_1, gamma_2]
+ * where a is the midsurface metric,
+ * b is similar to the second fundamental form of the surface, but the unit normal is replaced with the shell director
+ * gamma_1 and gamma_2 are the transverse shear in the two parametric directions
+ * Paper Equation 4.10 */
+ template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+ Dune::FieldVector<field_type, 8> SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::calculateGreenLagrangianStrains(
+ const KinematicVariables &kin)
+ {
+ Dune::FieldVector<RT, 8> egl;
+ // membrane
+ egl[0] = kin.A1andA2[0] * kin.ud1andud2[0] + 0.5 * kin.ud1andud2[0].two_norm2();
+ egl[1] = kin.A1andA2[1] * kin.ud1andud2[1] + 0.5 * kin.ud1andud2[1].two_norm2();
+ egl[2] = kin.a1anda2[0] * kin.a1anda2[1];
+
+ // bending
+ egl[3] = kin.a1anda2[0] * kin.td1Andtd2[0] - kin.A1andA2[0] * kin.t0d1Andt0d2[0];
+ egl[4] = kin.a1anda2[1] * kin.td1Andtd2[1] - kin.A1andA2[1] * kin.t0d1Andt0d2[1];
+ egl[5] = kin.a1anda2[0] * kin.td1Andtd2[1] + kin.a1anda2[1] * kin.td1Andtd2[0] - kin.A1andA2[0] * kin.t0d1Andt0d2[1]
+ - kin.A1andA2[1] * kin.t0d1Andt0d2[0];
+
+ // transverse shear
+ egl[6] = kin.a1anda2[0] * kin.t - kin.A1andA2[0] * kin.t0;
+ egl[7] = kin.a1anda2[1] * kin.t - kin.A1andA2[1] * kin.t0;
+
+ return egl;
+ }
+
+ /** \brief Calculates all kinematic quantities that are needed for strain calculation
+ *
+ * Paper Equations 6.4 - 6.8
+ */
+ template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+ template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
+ typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
+ auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::kinematicVariablesFactory(
+ const Element &element, const LocalDirectorFunction &directorFunction, const LocalDirectorReferenceFunction &directorReferenceFunction,
+ const LocalMidSurfaceFunction &midSurfaceFunction, const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
+ const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos)
+ {
+ KinematicVariables kin{};
+
+ const auto jInvT = element.geometry().jacobianInverseTransposed(quadPos);
+
+ kin.t = directorFunction.evaluate(quadPos).globalCoordinates();
+ kin.t0 = directorReferenceFunction.evaluate(quadPos).globalCoordinates();
+ kin.a1anda2 = transpose(midSurfaceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+ kin.A1andA2 = transpose(midSurfaceReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+ kin.ud1andud2 = transpose(midSurfaceDisplacementFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+ kin.t0d1Andt0d2 = transpose(directorReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+ kin.td1Andtd2 = transpose(directorFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+
+ return kin;
+ }
+
+ template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+ auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::getReferenceLocalConfigurations(
+ const typename Basis::LocalView &localView) const {
+ using namespace Dune::TypeTree::Indices;
+ const int nDofs0 = localView.tree().child(_0, 0).finiteElement().size();
+ const int nDofs1 = localView.tree().child(_1, 0).finiteElement().size();
+
+ std::vector<RealTuple<double, 3>> localConfiguration0(nDofs0);
+ std::vector<UnitVector<double, 3>> localConfiguration1(nDofs1);
+
+ for (int i = 0; i < nDofs0 + nDofs1; i++) {
+ int localIndexI = 0;
+ if (i < nDofs0) {
+ auto &node = localView.tree().child(_0, 0);
+ localIndexI = node.localIndex(i);
+ } else {
+ auto &node = localView.tree().child(_1, 0);
+ localIndexI = node.localIndex(i - nDofs0);
+ }
+ auto multiIndex = localView.index(localIndexI);
+
+ // The CompositeBasis number is contained in multiIndex[0]
+ // multiIndex[1] contains the actual index
+ if (multiIndex[0] == 0)
+ localConfiguration0[i] = midSurfaceRefConfig[multiIndex[1]];
+ else if (multiIndex[0] == 1)
+ localConfiguration1[i - nDofs0] = directorRefConfig[multiIndex[1]];
+ }
+ return std::make_tuple(localConfiguration0, localConfiguration1);
+ }
+
+ /** \brief Computes the Simo-Fox shell energy
+ *
+ * The energy is 0.5 * S * E = 0.5 * transpose(E) * Cmat * E, where
+ * S are the stress resultants [membrane forces, bending moments, transverse shear forces]
+ * E are the components of the Green-Lagrangian strains [membrane strains, bending , tranverse shear]
+ * see for details Paper Equation 4.11,4.10 and 10.1
+ */
+ template <class Basis, template <int, typename, typename, typename> typename LocalFEFunction, typename field_type>
+ typename SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::RT
+ SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::energy(
+ const typename Basis::LocalView &localView, const std::vector<RealTuple<field_type, 3>> &localMidSurfaceConfiguration,
+ const std::vector<UnitVector<field_type, 3>> &localDirectorConfiguration) const
+ {
+ auto element = localView.element();
+
+ using namespace Dune::Indices;
+ const auto &midSurfaceElement = LocalFiniteElementFactory<Basis, 0>::get(localView, _0);
+ const auto &directorElement = LocalFiniteElementFactory<Basis, 1>::get(localView, _1);
+
+ const auto [localRefMidSurfaceConfiguration, localRefDirectorConfiguration] = getReferenceLocalConfigurations(localView);
+
+ std::vector<RealTuple<field_type, 3>> displacements;
+ displacements.reserve(localMidSurfaceConfiguration.size());
+ for (size_t i = 0; i < localMidSurfaceConfiguration.size(); ++i)
+ displacements.emplace_back(localMidSurfaceConfiguration[i].globalCoordinates() - localRefMidSurfaceConfiguration[i].globalCoordinates());
+
+ const LocalMidSurfaceFunctionType localMidSurfaceFunction(midSurfaceElement, localMidSurfaceConfiguration);
+ const LocalMidSurfaceFunctionType localMidSurfaceDisplacementFunction(midSurfaceElement, displacements);
+ const LocalMidSurfaceReferenceFunctionType localMidSurfaceReferenceFunction(midSurfaceElement, localRefMidSurfaceConfiguration);
+ const LocalDirectorFunctionType localDirectorFunction(directorElement, localDirectorConfiguration);
+ const LocalDirectorReferenceFunctionType localDirectorReferenceFunction(directorElement, localRefDirectorConfiguration);
+
+ const int quadOrder = (element.type().isSimplex()) ? std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order())
+ : std::max(midSurfaceElement.localBasis().order(), directorElement.localBasis().order()) + 1;
+
+ const auto &quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+ RT energy = 0;
+ for (const auto &curQuad : quad) {
+ const Dune::FieldVector<DT, gridDim> &quadPos = curQuad.position();
+
+ const KinematicVariables kin
+ = kinematicVariablesFactory(element, localDirectorFunction, localDirectorReferenceFunction, localMidSurfaceFunction,
+ localMidSurfaceReferenceFunction, localMidSurfaceDisplacementFunction, quadPos);
+
+ const DT integrationElement = element.geometry().integrationElement(quadPos);
+ const FieldVector<field_type, 8> Egl = calculateGreenLagrangianStrains(kin);
+ energy += 0.5 * Egl * (CMat_ * Egl) * curQuad.weight() * integrationElement;
+ }
+
+ //////////////////////////////////////////////////////////////////////////////
+ // Assemble boundary contributions
+ //////////////////////////////////////////////////////////////////////////////
+
+ if (not neumannFunction_) return energy;
+
+ for (auto &&it : intersections(neumannBoundary_->gridView(), element))
+ {
+ if (not neumannBoundary_ or not neumannBoundary_->contains(it)) continue;
+
+ const auto &quadLine = QuadratureRules<DT, gridDim - 1>::rule(it.type(), quadOrder);
+
+ for (const auto &curQuad : quadLine)
+ {
+ // Local position of the quadrature point
+ const FieldVector<DT, gridDim> &quadPos = it.geometryInInside().global(curQuad.position());
+
+ const DT integrationElement = it.geometry().integrationElement(curQuad.position());
+
+ // The value of the local function
+ RealTuple<field_type, 3> deformationValue = localMidSurfaceFunction.evaluate(quadPos);
+
+ // Value of the Neumann data at the current position
+ auto neumannValue = neumannFunction_(it.geometry().global(curQuad.position()));
+
+ // Only translational dofs are affected by the Neumann force
+ for (size_t i = 0; i < neumannValue.size(); i++)
+ energy -= (neumannValue[i] * deformationValue.globalCoordinates()[i]) * curQuad.weight() * integrationElement;
+ }
+ }
+ return energy;
+ }
+} // namespace Dune::GFE
+#endif // DUNE_GFE_SIMOFOX_ENERGY_HH
diff --git a/problems/simofox-cantilever-dirichlet-values.py b/problems/simofox-cantilever-dirichlet-values.py
new file mode 100644
index 00000000..79695fd4
--- /dev/null
+++ b/problems/simofox-cantilever-dirichlet-values.py
@@ -0,0 +1,16 @@
+import math
+
+class DirichletValues:
+ def __init__(self, homotopyParameter):
+ self.homotopyParameter = homotopyParameter
+
+ def deformation(self, x):
+ # Dirichlet b.c. simply clamp the shell in the reference configuration
+ out = [x[0], x[1], 0]
+
+ return out
+
+ def director(self, x):
+ dir = [0, 0, 1]
+ return dir
+
diff --git a/problems/simofox-cantilever.parset b/problems/simofox-cantilever.parset
new file mode 100644
index 00000000..1eab7e5b
--- /dev/null
+++ b/problems/simofox-cantilever.parset
@@ -0,0 +1,99 @@
+#############################################
+# Grid parameters
+#############################################
+
+structuredGrid = true
+
+# bounding box
+lower = 0 0
+upper = 10 5
+
+elements = 10 1
+
+# Number of grid levels
+numLevels = 1
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 10
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+trustRegionScaling = 1 1 1 0.01 0.01
+
+# Initial trust-region radius
+initialTrustRegionRadius = 10
+
+# Number of multigrid iterations per trust-region step
+numIt = 400
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 1
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-6
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-6
+
+# Measure convergence
+instrumented = 0
+
+############################
+# Material parameters
+############################
+
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+# shell thickness
+thickness = 1
+
+# Lame parameters
+mu = 100000
+lambda = 0
+
+# Shear correction factor
+kappa = 1
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = simofox-cantilever
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+### Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 9.99"
+
+### Neumann values
+neumannValues = 0 0 320
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0]"
+
+#startFromFile = yes
+#initialIterateFilename = initial_iterate.vtu
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 33b44cda..4962aa5d 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,7 +1,8 @@
set(programs compute-disc-error
cosserat-continuum
gradient-flow
- harmonicmaps)
+ harmonicmaps
+ simofoxshell)
# rodobstacle)
if (${DUNE_COMMON_VERSION} VERSION_GREATER_EQUAL 2.8)
set (programs rod3d ${programs})
diff --git a/src/simofoxshell.cc b/src/simofoxshell.cc
new file mode 100644
index 00000000..eda66c50
--- /dev/null
+++ b/src/simofoxshell.cc
@@ -0,0 +1,390 @@
+#include <config.h>
+
+#include <array>
+#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/utility/structuredgridfactory.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/uggrid.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/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/simofoxenergy.hh>
+#include <dune/gfe/cosseratvtkwriter.hh>
+#include <dune/gfe/embeddedglobalgfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+
+#if HAVE_DUNE_VTK
+#include <dune/vtk/vtkreader.hh>
+#endif
+
+template <int dim, class ctype, class LocalFiniteElement, class TS>
+using LocalFEFunction = Dune::GFE::LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TS>;
+//using LocalFEFunction = LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TS>;
+
+// Order of the approximation space for the midsurface position
+const int midsurfaceOrder = 1;
+
+// Order of the approximation space for the director
+const int directorOrder = 1;
+
+using namespace Dune;
+
+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");
+
+ Python::runStream()
+ << std::endl << "import sys"
+ << std::endl << "import os"
+ << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
+ << std::endl;
+
+ using namespace TypeTree::Indices;
+ using SolutionType = TupleVector<std::vector<RealTuple<double,3> >, std::vector<UnitVector<double,3> > >;
+
+ // parse data file
+ ParameterTree parameterSet;
+ if (argc < 2)
+ DUNE_THROW(Exception, "Usage: ./simofoxshell inputfile.dat");
+
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+ // read solver settings
+ const auto numLevels = parameterSet.get<int>("numLevels");
+ const auto totalLoadSteps = parameterSet.get<int>("numHomotopySteps");
+ const auto tolerance = parameterSet.get<double>("tolerance");
+ const auto maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
+ const auto initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+ const auto multigridIterations = parameterSet.get<int>("numIt");
+ const auto nu1 = parameterSet.get<int>("nu1");
+ const auto nu2 = parameterSet.get<int>("nu2");
+ const auto mu = parameterSet.get<int>("mu");
+ const auto baseIterations = parameterSet.get<int>("baseIt");
+ const auto mgTolerance = parameterSet.get<double>("mgTolerance");
+ const auto baseTolerance = parameterSet.get<double>("baseTolerance");
+ const auto instrumented = parameterSet.get<bool>("instrumented");
+ std::string resultPath = parameterSet.get("resultPath", "");
+
+ /////////////////////////////////////////
+ // Create the grid
+ /////////////////////////////////////////
+ using Grid = UGGrid<2>;
+
+ std::shared_ptr<Grid> grid;
+
+ FieldVector<double, 2> lower(0), upper(1);
+
+ std::string structuredGridType = parameterSet["structuredGrid"];
+ if (parameterSet.get<bool>("structuredGrid"))
+ {
+ lower = parameterSet.get<FieldVector<double, 2> >("lower");
+ upper = parameterSet.get<FieldVector<double, 2> >("upper");
+
+ auto elements = parameterSet.get<std::array<unsigned int, 2> >("elements");
+ grid = StructuredGridFactory<Grid>::createCubeGrid(lower, upper, elements);
+
+ } else {
+ auto path = parameterSet.get<std::string>("path");
+ auto 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")
+ grid = std::shared_ptr<Grid>(GmshReader<Grid>::read(path + "/" + gridFile));
+ else if (suffix == ".vtu" or suffix == ".vtp")
+#if HAVE_DUNE_VTK
+ grid = VtkReader<Grid>::createGridFromFile(path + "/" + gridFile);
+#else
+ DUNE_THROW(NotImplemented, "Please install dune-vtk for VTK reading support!");
+#endif
+ }
+
+ grid->globalRefine(numLevels - 1);
+ grid->loadBalance();
+
+ if (mpiHelper.rank() == 0)
+ std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
+
+ using GridView = Grid::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ using namespace Dune::Functions::BasisFactory;
+
+ auto compositeBasis = makeBasis(gridView,
+ composite(
+ power<2>(lagrange<midsurfaceOrder>()),
+ power<2>(lagrange<directorOrder>())));
+
+ typedef Dune::Functions::LagrangeBasis<GridView, midsurfaceOrder> MidsurfaceFEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView, directorOrder> DirectorFEBasis;
+
+ MidsurfaceFEBasis midsurfaceFEBasis(gridView);
+ DirectorFEBasis directorFEBasis(gridView);
+
+ ///////////////////////////////////////////
+ // Read Dirichlet values
+ ///////////////////////////////////////////
+
+ BitSetVector<1> dirichletVertices(gridView.size(2), false);
+ BitSetVector<1> neumannVertices(gridView.size(2), 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<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));
+
+ for (auto &&vertex: vertices(gridView))
+ {
+ bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
+ dirichletVertices[indexSet.index(vertex)] = isDirichlet;
+
+ bool isNeumann = pythonNeumannVertices(vertex.geometry().corner(0));
+ neumannVertices[indexSet.index(vertex)] = isNeumann;
+ }
+
+ BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+ BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+
+ if (mpiHelper.rank() == 0)
+ std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+
+ BitSetVector<1> deformationDirichletNodes(midsurfaceFEBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary, midsurfaceFEBasis, deformationDirichletNodes);
+
+ BitSetVector<1> neumannNodes(midsurfaceFEBasis.size(), false);
+ constructBoundaryDofs(neumannBoundary, directorFEBasis, neumannNodes);
+
+ BitSetVector<3> deformationDirichletDofs(midsurfaceFEBasis.size(), false);
+ for (size_t i = 0; i < midsurfaceFEBasis.size(); i++)
+ if (deformationDirichletNodes[i][0])
+ for (int j = 0; j < 3; j++)
+ deformationDirichletDofs[i][j] = true;
+
+ BitSetVector<1> orientationDirichletNodes(directorFEBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary, directorFEBasis, orientationDirichletNodes);
+
+ BitSetVector<2> orientationDirichletDofs(directorFEBasis.size(), false);
+ for (size_t i = 0; i < directorFEBasis.size(); i++)
+ if (orientationDirichletNodes[i][0])
+ for (int j = 0; j < 2; j++)
+ orientationDirichletDofs[i][j] = true;
+
+ ////////////////////////////
+ // Initial iterate
+ ////////////////////////////
+
+ SolutionType x;
+
+ x[_0].resize(midsurfaceFEBasis.size());
+
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ auto pythonInitialDeformation = Python::make_function<FieldVector<double, 3> >(Python::evaluate(lambda));
+
+ auto deformationPowerBasis = makeBasis(gridView,power<3>(lagrange<midsurfaceOrder>()));
+
+ std::vector<FieldVector<double, 3> > v;
+ Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
+ std::copy(v.begin(), v.end(), x[_0].begin());
+
+ x[_1].resize(directorFEBasis.size());
+
+ // The code currently assumes that the grid resides in the x-y plane
+ // Therefore, all references directors point into z direction
+ // If the reference is not planar one has to calculate the reference nodal directors.
+ // E.g. for bilinear grid elements this can be just the average of the 4 element normals
+ const FieldVector<double, 3> referenceDirectorVector({0, 0, 1}); //only plain reference!
+ std::fill(x[_1].begin(), x[_1].end(), referenceDirectorVector);
+
+ const SolutionType x0 = x;
+
+ ////////////////////////////////////////////////////////
+ // Main homotopy loop
+ ////////////////////////////////////////////////////////
+
+ for (int loadStep = 0; loadStep < totalLoadSteps; loadStep++)
+ {
+ const double loadFactor = (loadStep + 1) * (1.0 / totalLoadSteps);
+ if (mpiHelper.rank() == 0)
+ std::cout << "Homotopy step: " << loadStep << ", parameter: " << loadFactor << 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");
+
+ auto neumannFunction = [&](FieldVector<double, 2>) {
+ auto nV = neumannValues;
+ nV *= loadFactor;
+ return nV;
+ };
+
+ // Output initial iterate (of homotopy loop)
+ auto directorPowerBasis = makeBasis(gridView, power<3>(lagrange<directorOrder>()));
+
+ BlockVector<FieldVector<double, 3> > displacementInitial(compositeBasis.size({0}));
+ //calculates the global coordinates of midsurface displacements into displacementInitial
+ std::transform(x[_0].begin(),x[_0].end(),x0[_0].begin(),displacementInitial.begin(),[](const auto& x_0Entry,const auto& x0_0Entry){
+ return x_0Entry.globalCoordinates()-x0_0Entry.globalCoordinates();
+ });
+
+ //copies the global coordinates of the directors to directorInitial for a VTKWriter usable format
+ BlockVector<FieldVector<double, 3> > directorInitial(compositeBasis.size({1}));
+ std::transform(x[_1].begin(),x[_1].end(),directorInitial.begin(),[](const auto& x_1Entry){
+ return x_1Entry.globalCoordinates();
+ });
+
+ auto displacementFunctionInitial = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(deformationPowerBasis,
+ displacementInitial);
+ auto directorFunctionInitial = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(directorPowerBasis,
+ directorInitial);
+ // We need to subsample, because VTK cannot natively display real second-order functions
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(midsurfaceOrder - 1));
+ vtkWriter.addVertexData(displacementFunctionInitial, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+ vtkWriter.addVertexData(directorFunctionInitial, VTK::FieldInfo("director", VTK::FieldInfo::Type::scalar, 3));
+ vtkWriter.write(resultPath + "simo-fox_homotopy_" + "_" + std::to_string(neumannValues[2]) + "_" + std::to_string(0));
+
+ if (mpiHelper.rank() == 0) {
+ std::cout << "Material parameters:" << std::endl;
+ materialParameters.report();
+ }
+
+ // Assembler using ADOL-C
+ Dune::GFE::SimoFoxEnergyLocalStiffness<decltype(compositeBasis), LocalFEFunction,adouble> simoFoxEnergyADOLCLocalStiffness(materialParameters,
+ &neumannBoundary,
+ neumannFunction,
+ nullptr, x0);
+
+ MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
+ RealTuple<double,3>,
+ UnitVector<double,3> > localGFEADOLCStiffness(&simoFoxEnergyADOLCLocalStiffness);
+
+ MixedGFEAssembler<decltype(compositeBasis),
+ RealTuple<double,3>, UnitVector<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
+
+ // /////////////////////////////////////////////////
+ // Create a Riemannian trust-region solver
+ // /////////////////////////////////////////////////
+
+ MixedRiemannianTrustRegionSolver<Grid,
+ decltype(compositeBasis),
+ MidsurfaceFEBasis, RealTuple<double,3>,
+ DirectorFEBasis, UnitVector<double,3> > solver;
+
+ solver.setup(*grid,
+ &assembler,
+ midsurfaceFEBasis,
+ directorFEBasis,
+ x,
+ deformationDirichletDofs,
+ orientationDirichletDofs,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ solver.setScaling(parameterSet.get<FieldVector<double, 5> >("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(loadFactor);
+
+ // Extract object member functions as Dune functions
+ auto deformationDirichletValues = Python::make_function<FieldVector<double, 3> >(dirichletValuesPythonObject.get("deformation"));
+
+ std::vector<FieldVector<double, 3> > ddV;
+ Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+
+ for (size_t j = 0; j < x[_0].size(); j++)
+ if (deformationDirichletNodes[j][0])
+ x[_0][j] = ddV[j];
+
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ solver.setInitialIterate(x);
+ solver.solve();
+
+ x = solver.getSol();
+
+ // Output result of each load step
+ std::stringstream iAsAscii;
+ iAsAscii << loadStep + 1;
+
+ //calculates the global coordinates of midsurface displacements into displacementInitial
+ BlockVector<FieldVector<double, 3> > displacement(compositeBasis.size({0}));
+ std::transform(x[_0].begin(),x[_0].end(),x0[_0].begin(),displacement.begin(),[](const auto& x_0Entry,const auto& x0_0Entry){
+ return x_0Entry.globalCoordinates()-x0_0Entry.globalCoordinates();
+ });
+
+ //copies the global coordinates of the directors to directorInitial for a VTKWriter usable format
+ BlockVector<FieldVector<double, 3> > director(compositeBasis.size({1}));
+ std::transform(x[_1].begin(),x[_1].end(),director.begin(),[](const auto& x_1Entry){
+ return x_1Entry.globalCoordinates();
+ });
+
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(deformationPowerBasis,
+ displacement);
+ auto directorFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double, 3> >(directorPowerBasis,
+ director);
+ // We need to subsample, because VTK cannot natively display real second-order functions
+ // SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(midsurfaceOrder - 1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
+ vtkWriter.addVertexData(directorFunction, VTK::FieldInfo("director", VTK::FieldInfo::Type::scalar, 3));
+ vtkWriter.write(resultPath + "simo-fox_homotopy_" + "_" + std::to_string(neumannValues[2]) + "_" + std::to_string(loadStep + 1));
+ }
+}
+catch (Exception &e) {
+ std::cout << e.what() << std::endl;
+}
--
GitLab