diff --git a/dune/gfe/surfacecosseratstressassembler.hh b/dune/gfe/surfacecosseratstressassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b8246b8fa47a2826bbc63a915bd21924a4fdb15c
--- /dev/null
+++ b/dune/gfe/surfacecosseratstressassembler.hh
@@ -0,0 +1,308 @@
+#ifndef DUNE_GFE_SURFACECOSSERATSTRESSASSEMBLER_HH
+#define DUNE_GFE_SURFACECOSSERATSTRESSASSEMBLER_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/gfe/linearalgebra.hh>
+#include <dune/gfe/rotation.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+
+#include <dune/matrix-vector/transpose.hh>
+
+namespace Dune::GFE {
+ /** \brief An assembler that can calculate the norms of specific stress tensors for each element for an output by film-on-substrate
+
+ \tparam BasisOrderD Basis used for the displacement
+ \tparam BasisOrderR Basis used for the rotation
+ \tparam TargetSpaceD Target space for the Displacement
+ \tparam TargetSpaceR Target space for the Rotation
+ */
+ template <class BasisOrderD, class BasisOrderR, class TargetSpaceD, class TargetSpaceR>
+ class SurfaceCosseratStressAssembler
+ {
+ public:
+ const static int dim = TargetSpaceD::dimension;
+ using GridView = typename BasisOrderD::GridView;
+ using VectorD = std::vector<TargetSpaceD>;
+ using VectorR = std::vector<TargetSpaceR>;
+
+ BasisOrderD basisOrderD_;
+ BasisOrderR basisOrderR_;
+
+ SurfaceCosseratStressAssembler(const BasisOrderD basisOrderD,
+ const BasisOrderR basisOrderR)
+ : basisOrderD_(basisOrderD),
+ basisOrderR_(basisOrderR)
+ {}
+
+ /** \brief Calculate the norm of the 1st-Piola-Kirchhoff-Stress-Tensor and the Cauchy-Stress-Tensor for each element
+
+ The 1st-Piola-Kirchhoff-Stress-Tensor is the derivative of the energy density with respect to the deformation gradient
+ - Calculate the deformation gradient for each element using the basis functions and
+ their gradients; then add them up using the localConfiguration
+ - Evaluate the deformation gradient at each quadrature point using the respective quadrature rule with the given order
+ - Evaluate the density function and tape the evaluation - then use ADOLC to evaluate the derivative (∂/∂F) W(F)
+ - The derivative is then a dim x dim matrix
+ - Then calculate the final stressTensor of the element by averagin over the quadrature points using the quadrature
+ weights and the reference element volume
+ \param x Coefficient vector for the displacement
+ \param elasticDensity Energy density function
+ \param int Order of the quadrature rule
+ \param stressSubstrate1stPiolaKirchhoffTensor Vector containing the the 1st-Piola-Kirchhoff-Stress-Tensor for each element
+ \param stressSubstrateCauchyTensor Vector containing the Cauchy-Stress-Tensor for each element
+ */
+ template <class Density>
+ void assembleSubstrateStress(
+ const VectorD x,
+ const Density* elasticDensity,
+ const int quadOrder,
+ std::vector<FieldMatrix<double,dim,dim>>& stressSubstrate1stPiolaKirchhoffTensor,
+ std::vector<FieldMatrix<double,dim,dim>>& stressSubstrateCauchyTensor)
+ {
+
+ std::cout << "Calculating the Frobenius norm of the 1st-Piola-Kirchhoff-Stress-Tensor ( (∂/∂F) W(F) )" << std::endl
+ << "and the Frobenius norm of the Cauchy-Stress-Tensor (1/det(F) * (∂/∂F) W(F) * F^T) of the substrate..." << std::endl;
+
+ auto xFlat = Functions::istlVectorBackend(x);
+ static constexpr auto partitionType = Partitions::interiorBorder;
+ MultipleCodimMultipleGeomTypeMapper<GridView> elementMapper(basisOrderD_.gridView(),mcmgElementLayout());
+ stressSubstrate1stPiolaKirchhoffTensor.resize(elementMapper.size());
+ stressSubstrateCauchyTensor.resize(elementMapper.size());
+
+ for (const auto& element : elements(basisOrderD_.gridView(), partitionType))
+ {
+ auto localViewOrderD = basisOrderD_.localView();
+ localViewOrderD.bind(element);
+ size_t nDofsOrderD = localViewOrderD.tree().size();
+
+ // Extract values at this element
+ std::vector<double> localConfiguration(nDofsOrderD);
+ for (int i=0; i<nDofsOrderD; i++)
+ localConfiguration[i] = xFlat[localViewOrderD.index(i)]; //localViewOrderD.index(i) is a multi-index
+
+ //Store the reference gradient and the gradients for this element
+ const auto& lFEOrderD = localViewOrderD.tree().child(0).finiteElement();
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients(lFEOrderD.size());
+ std::vector<FieldMatrix<double,1,dim> > gradients(lFEOrderD.size());
+
+ auto evaluateAtPoint = [&](FieldVector<double,3> pointGlobal, FieldVector<double,3> pointLocal) -> std::vector<FieldMatrix<double,dim,dim>>{
+ std::vector<FieldMatrix<double,dim,dim>> stressTensors(2);
+
+ const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(pointLocal);
+
+ // Get gradients of shape functions
+ lFEOrderD.localBasis().evaluateJacobian(pointLocal, referenceGradients);
+
+ // Compute gradients of Base functions
+ for (size_t i=0; i<gradients.size(); ++i)
+ gradients[i] = referenceGradients[i] * transpose(jacobianInverseTransposed);
+ //jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+ // Deformation gradient in vector form
+ size_t nDoubles = dim*dim;
+ std::vector<double> deformationGradientFlat(nDoubles);
+ for (size_t i=0; i<nDoubles; i++)
+ deformationGradientFlat[i] = 0;
+ for (size_t i=0; i<gradients.size(); i++)
+ for (size_t j=0; j<dim; j++)
+ for (size_t k=0; k<dim; k++)
+ deformationGradientFlat[dim*j + k] += localConfiguration[ localViewOrderD.tree().child(j).localIndex(i)]*gradients[i][0][k];
+
+ double pureDensity = 0;
+
+ trace_on(0);
+
+ FieldMatrix<adouble,dim,dim> deformationGradient(0);
+ for (size_t j=0; j<dim; j++)
+ for (size_t k=0; k<dim; k++)
+ deformationGradient[j][k] <<= deformationGradientFlat[dim*j + k];
+
+ // Tape the actual calculation
+ adouble density = 0;
+ try {
+ density = (*elasticDensity)(pointGlobal, deformationGradient);
+ } catch (Exception &e) {
+ trace_off(0);
+ throw e;
+ }
+ density >>= pureDensity;
+
+ trace_off(0);
+
+ // Compute the actual gradient
+ std::vector<double> localStressFlat(nDoubles);
+ gradient(0,nDoubles,deformationGradientFlat.data(),localStressFlat.data());
+
+ FieldMatrix<double,dim,dim> localStress(0);
+ for (size_t j=0; j<dim; j++)
+ for (size_t k=0; k<dim; k++)
+ localStress[j][k] = localStressFlat[dim*j + k];
+
+ stressTensors[0] = localStress; // 1st-Piola-Kirchhoff-Stress-Tensor
+
+ FieldMatrix<double,dim,dim> deformationGradientTransposed(0);
+ for (size_t j=0; j<dim; j++)
+ for (size_t k=0; k<dim; k++)
+ deformationGradient[j][k] >>= deformationGradientTransposed[k][j];
+
+ localStress /= deformationGradientTransposed.determinant();
+ localStress = localStress * deformationGradientTransposed;
+ stressTensors[1] = localStress; // Cauchy-Stress-Tensor
+
+ return stressTensors;
+ };
+
+ //Call evaluateAtPoint for all points in the quadrature rule
+ const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), quadOrder);
+ stressSubstrate1stPiolaKirchhoffTensor[elementMapper.index(element)] = 0;
+ stressSubstrateCauchyTensor[elementMapper.index(element)] = 0;
+
+ for (size_t pt=0; pt<quad.size(); pt++) {
+ auto pointLocal = quad[pt].position();
+ auto pointGlobal = element.geometry().global(pointLocal);
+ auto stressTensors = evaluateAtPoint(pointGlobal, pointLocal);
+ stressSubstrate1stPiolaKirchhoffTensor[elementMapper.index(element)] += stressTensors[0] * quad[pt].weight()/referenceElement(element).volume();
+ stressSubstrateCauchyTensor[elementMapper.index(element)] += stressTensors[1] * quad[pt].weight()/referenceElement(element).volume();
+ }
+ }
+ }
+ /** \brief Calculate the norm of the Biot-Type-Stress-Tensor of the shell
+
+ The formula for Biot-Type-Stress-Tensor of the Cosserat shell given by (4.11) in
+ Ghiba, Bîrsan, Lewintan, Neff, March 2020: "The isotropic Cosserat shell model including terms up to $O(h^5)$. Part I: Derivation in matrix notation"
+ \param rot Coefficient vector for the rotation
+ \param x Coefficient vector for the displacement
+ \param xInitial Coefficient vector for the stress-free configuration of the shell, used to calculate nablaTheta
+ \param lameF Function assinging the lamé parameters to a given point
+ \param mu_c Cosserat couple modulus
+ \param shellBoundary BoundaryPatch containing the elements that actually belong to the shell
+ \param order Order of the quadrature rule
+ \param stressShellBiotTensor Vector containing the Biot-Stress-Tensor for each element
+ */
+ void assembleShellStress(
+ const VectorR rot,
+ const VectorD x,
+ const VectorD xInitial,
+ const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF,
+ const double mu_c,
+ const BoundaryPatch<GridView> shellBoundary,
+ const int quadOrder,
+ std::vector<FieldMatrix<double,dim,dim>>& stressShellBiotTensor)
+ {
+ std::cout << "Calculating the Frobenius norm of the Biot-Type-Stress-Tensor of the shell..." << std::endl;
+ auto xFlat = Functions::istlVectorBackend(x);
+ auto xInitialFlat = Functions::istlVectorBackend(xInitial);
+
+ MultipleCodimMultipleGeomTypeMapper<GridView> elementMapper(basisOrderD_.gridView(),mcmgElementLayout());
+ stressShellBiotTensor.resize(elementMapper.size());
+
+ static constexpr auto partitionType = Partitions::interiorBorder;
+ for (const auto& element : elements(basisOrderD_.gridView(), partitionType))
+ {
+ stressShellBiotTensor[elementMapper.index(element)] = 0;
+
+ int intersectionCounter = 0;
+ for (auto&& it : intersections(shellBoundary.gridView(), element)) {
+ FieldMatrix<double,dim,dim> stressTensorThisIntersection(0);
+
+ //Continue if the element does not intersect with the shell boundary
+ if (not shellBoundary.contains(it))
+ continue;
+
+ //LocalView for the basisOrderD_
+ auto localViewOrderD = basisOrderD_.localView();
+ localViewOrderD.bind(element);
+ size_t nDofsOrderD = localViewOrderD.tree().size();
+
+ // Extract local configuration at this element
+ std::vector<double> localConfiguration(nDofsOrderD);
+ std::vector<double> localConfigurationInitial(nDofsOrderD);
+ for (int i=0; i<nDofsOrderD; i++) {
+ localConfiguration[i] = xFlat[localViewOrderD.index(i)];
+ localConfigurationInitial[i] = xInitialFlat[localViewOrderD.index(i)];
+ }
+
+ const auto& lFEOrderD = localViewOrderD.tree().child(0).finiteElement();
+ //Store the reference gradient and the gradients for *this element*
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients(lFEOrderD.size());
+ std::vector<FieldMatrix<double,1,dim> > gradients(lFEOrderD.size());
+
+ //LocalView for the basisOrderR_
+ auto localViewOrderR = basisOrderR_.localView();
+ localViewOrderR.bind(element);
+ const auto& lFEOrderR = localViewOrderR.tree().child(0).finiteElement();
+ VectorR localConfigurationRot(lFEOrderR.size());
+ for (int i=0; i<localConfigurationRot.size(); i++)
+ localConfigurationRot[i] = rot[localViewOrderR.index(i)[0]];//localViewOrderR.index(i) is a multiindex, its first entry is the actual index
+ typedef LocalGeodesicFEFunction<dim, double, decltype(lFEOrderR), TargetSpaceR> LocalGFEFunctionType;
+ LocalGFEFunctionType localGeodesicFEFunction(lFEOrderR,localConfigurationRot);
+
+ auto evaluateAtPoint = [&](FieldVector<double,3> pointGlobal, FieldVector<double,3> pointLocal3d) -> FieldMatrix<double,dim,dim>{
+ Dune::FieldMatrix<double,dim,dim> nablaTheta;
+
+ const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(pointLocal3d);
+
+ // Get gradients of shape functions
+ lFEOrderD.localBasis().evaluateJacobian(pointLocal3d, referenceGradients);
+
+ // Compute gradients of Base functions at this element
+ for (size_t i=0; i<gradients.size(); i++)
+ gradients[i] = referenceGradients[i] * transpose(jacobianInverseTransposed);
+
+ //jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+ // Deformation gradient - call this U_es_minus_Id already
+ FieldMatrix<double,dim,dim> U_es_minus_Id(0);
+ for (size_t i=0; i<gradients.size(); i++)
+ for (size_t j=0; j<dim; j++){
+ U_es_minus_Id[j].axpy(localConfiguration[ localViewOrderD.tree().child(j).localIndex(i)],gradients[i][0]);
+ nablaTheta[j].axpy(localConfigurationInitial[ localViewOrderD.tree().child(j).localIndex(i)],gradients[i][0]);
+ }
+
+ TargetSpaceR value = localGeodesicFEFunction.evaluate(pointLocal3d);
+ FieldMatrix<double,dim,dim> rotationMatrix(0);
+ FieldMatrix<double,dim,dim> rotationMatrixTransposed(0);
+ value.matrix(rotationMatrix);
+
+ MatrixVector::transpose(rotationMatrix, rotationMatrixTransposed);
+
+ U_es_minus_Id.leftmultiply(rotationMatrixTransposed); // Attention: The rotation here is already is Q_e, we don't have to multiply with Q_0!!!
+
+ nablaTheta.invert();
+ U_es_minus_Id.rightmultiply(nablaTheta);
+
+ for (size_t j = 0; j < dim; j++)
+ U_es_minus_Id[j][j] -= 1.0;
+
+ auto lameConstants = lameF(pointGlobal);
+ double mu = lameConstants[0];
+ double lambda = lameConstants[1];
+
+ FieldMatrix<double,dim,dim> localStressShell(0);
+ localStressShell = 2*mu*GFE::sym(U_es_minus_Id) + 2*mu_c*GFE::skew(U_es_minus_Id);
+ for (size_t j = 0; j < dim; j++)
+ localStressShell[j][j] += lambda*GFE::trace(GFE::sym(U_es_minus_Id));
+ return localStressShell;
+ };
+
+ //Call evaluateAtPoint for all points in the quadrature rule
+ const auto& quad = Dune::QuadratureRules<double, dim-1>::rule(it.type(), quadOrder); //Quad rule on the boundary
+
+ for (size_t pt=0; pt<quad.size(); pt++) {
+ auto pointLocal2d = quad[pt].position();
+ auto pointLocal3d = it.geometryInInside().global(pointLocal2d);
+ auto pointGlobal = it.geometry().global(pointLocal2d);
+
+ auto stressTensor = evaluateAtPoint(pointGlobal, pointLocal3d);
+ stressTensorThisIntersection += stressTensor * quad[pt].weight()/referenceElement(it.inside()).volume();
+ }
+ stressShellBiotTensor[elementMapper.index(element)] += stressTensorThisIntersection;
+ intersectionCounter++;
+ }
+ if (intersectionCounter >= 1)
+ stressShellBiotTensor[elementMapper.index(element)] /= intersectionCounter;
+ }
+ }
+ };
+}
+#endif
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ada923056911e2df35c4ef47120577a3c080cd47..e3e958a6d9fccd3a32402361e945f88671491889 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -5,6 +5,9 @@ set(programs compute-disc-error
rod3d)
# rodobstacle)
+if (${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.8)
+ set(programs film-on-substrate-stress-plot ${programs})
+endif()
if (dune-parmg_FOUND AND dune-curvedgeometry_FOUND AND ${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.8)
set(programs film-on-substrate ${programs})
endif()
diff --git a/src/film-on-substrate-stress-plot.cc b/src/film-on-substrate-stress-plot.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a5a7c37add55d933f4a10eaf903be108e2fb27a7
--- /dev/null
+++ b/src/film-on-substrate-stress-plot.cc
@@ -0,0 +1,335 @@
+#include <iostream>
+#include <fstream>
+
+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
+#include <adolc/taping.h>
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/version.hh>
+
+#include <dune/elasticity/materials/exphenckydensity.hh>
+#include <dune/elasticity/materials/henckydensity.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+#include <dune/elasticity/materials/neohookedensity.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/dunepython.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dune/gfe/filereader.hh>
+#include <dune/gfe/rotation.hh>
+#include <dune/gfe/surfacecosseratstressassembler.hh>
+
+// grid dimension
+#ifndef WORLD_DIM
+# define WORLD_DIM 3
+#endif
+const int dim = WORLD_DIM;
+
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+using namespace Dune;
+using ValueType = adouble;
+
+int main (int argc, char *argv[]) try
+{
+ 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 << "import os"
+ << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
+ << std::endl;
+
+ // parse data file
+ ParameterTree parameterSet;
+
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+ /////////////////////////////////////////////////////////////
+ // CREATE THE GRID
+ /////////////////////////////////////////////////////////////
+ typedef UGGrid<dim> GridType;
+
+ std::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");
+
+ std::array<unsigned int,dim> 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 = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+ }
+
+ grid->setRefinementType(GridType::RefinementType::COPY);
+
+ // Surface Shell Boundary
+ std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("surfaceShellVerticesPredicate", "0") + std::string(")");
+ auto pythonSurfaceShellVertices = Python::make_function<bool>(Python::evaluate(lambda));
+
+ int numLevels = parameterSet.get<int>("numLevels");
+
+ while (numLevels > 0) {
+ for (auto&& e : elements(grid->leafGridView())){
+ bool isSurfaceShell = false;
+ for (int i = 0; i < e.geometry().corners(); i++) {
+ isSurfaceShell = isSurfaceShell || pythonSurfaceShellVertices(e.geometry().corner(i));
+ }
+ grid->mark(isSurfaceShell ? 1 : 0,e);
+ }
+
+ grid->adapt();
+
+ numLevels--;
+ }
+
+ grid->loadBalance();
+
+ if (grid->leafGridView().comm().size() > 1)
+ DUNE_THROW(Exception,
+ std::string("To create a stress plot, please use only one process, now there are ") + std::to_string(grid->leafGridView().comm().size()) + std::string(" procsses."));
+
+ typedef GridType::LeafGridView GridView;
+ GridView gridView = grid->leafGridView();
+
+ /////////////////////////////////////////////////////////////
+ // SURFACE SHELL BOUNDARY
+ /////////////////////////////////////////////////////////////
+
+ const GridView::IndexSet& indexSet = gridView.indexSet();
+ BitSetVector<1> surfaceShellVertices(gridView.size(dim), false);
+ for (auto&& v : vertices(gridView))
+ {
+ bool isSurfaceShell = pythonSurfaceShellVertices(v.geometry().corner(0));
+ surfaceShellVertices[indexSet.index(v)] = isSurfaceShell;
+ }
+ BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
+
+ /////////////////////////////////////////////////////////////
+ // FUNCTION SPACE
+ /////////////////////////////////////////////////////////////
+
+ using namespace Functions::BasisFactory;
+ auto basisOrderD = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<displacementOrder>()
+ ));
+
+ auto basisOrderR = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<rotationOrder>()
+ ));
+
+ /////////////////////////////////////////////////////////////
+ // INITIAL DATA
+ /////////////////////////////////////////////////////////////
+
+ // Read grid deformation information from the file specified in the parameter set via pathToOutput, deformationOutput, rotationOutput and initial deformation
+ const std::string pathToOutput = parameterSet.get("pathToOutput", "");
+
+ std::cout << "Reading in deformation file (" << "order is " << displacementOrder << "): " << pathToOutput + parameterSet.get<std::string>("deformationOutput") << std::endl;
+ auto deformationMap = Dune::GFE::transformFileToMap<dim>(pathToOutput + parameterSet.get<std::string>("deformationOutput"));
+ std::cout << "... done: The basis has " << basisOrderD.size() << " elements and the defomation file has " << deformationMap.size() << " entries." << std::endl;
+
+ const auto dimRotation = Rotation<double,dim>::embeddedDim;
+ std::unordered_map<std::string, FieldVector<double,dimRotation>> rotationMap;
+ if (parameterSet.hasKey("rotationOutput")) {
+ std::cout << "Reading in rotation file (" << "order is " << rotationOrder << "): " << pathToOutput + parameterSet.get<std::string>("rotationOutput") << std::endl;
+ rotationMap = Dune::GFE::transformFileToMap<dimRotation>(pathToOutput + parameterSet.get<std::string>("rotationOutput"));
+ }
+ const bool startFromFile = parameterSet.get<bool>("startFromFile");
+ std::unordered_map<std::string, FieldVector<double,dim>> initialDeformationMap;
+
+ auto gridDeformationLambda = std::string("lambda x: (") + parameterSet.get<std::string>("gridDeformation") + std::string(")");
+ auto gridDeformation = Python::make_function<FieldVector<double,dim> >(Python::evaluate(gridDeformationLambda));
+
+ if (startFromFile) {
+ std::cout << "Reading in file to the stress-free configuration of the shell (" << "order is " << displacementOrder << "): " << parameterSet.get("pathToGridDeformationFile", "") + parameterSet.get<std::string>("gridDeformationFile") << std::endl;
+ initialDeformationMap = Dune::GFE::transformFileToMap<dim>(parameterSet.get("pathToGridDeformationFile", "") + parameterSet.get<std::string>("gridDeformationFile"));
+ }
+
+ using DisplacementVector = std::vector<FieldVector<double,dim> >;
+ DisplacementVector x;
+ x.resize(basisOrderD.size());
+ DisplacementVector xInitial;
+ xInitial.resize(basisOrderD.size());
+ DisplacementVector displacement;
+ displacement.resize(basisOrderD.size());
+
+ Functions::interpolate(basisOrderD, x, [](FieldVector<double,dim> x){ return x; });
+ Functions::interpolate(basisOrderD, xInitial, [](FieldVector<double,dim> x){ return x; });
+
+ for (int i = 0; i < basisOrderD.size(); i++) {
+ std::stringstream stream;
+ stream << x[i];
+ //Look up the displacement for this vertex in the deformationMap
+ displacement[i] = deformationMap.at(stream.str());
+ x[i] += deformationMap.at(stream.str());
+ //In case an a stress-free file was provided: look up the displacement for this vertex in the initialDeformationMap
+ if (startFromFile) {
+ xInitial[i] += initialDeformationMap.at(stream.str());
+ } else {
+ xInitial[i] = gridDeformation(xInitial[i]);
+ }
+ }
+
+ using RotationVector = std::vector<Rotation<double,dim>>;
+ RotationVector rot;
+ rot.resize(basisOrderR.size());
+ DisplacementVector xOrderR;
+ xOrderR.resize(basisOrderR.size());
+ Functions::interpolate(basisOrderR, xOrderR, [](FieldVector<double,dim> x){ return x; });
+
+ using DirectorVector = std::vector<Dune::FieldVector<double,dim>>;
+ std::array<DirectorVector,3> rot_director;
+ rot_director[0].resize(basisOrderR.size());
+ rot_director[1].resize(basisOrderR.size());
+ rot_director[2].resize(basisOrderR.size());
+
+ for (int i = 0; i < basisOrderR.size(); i++) {
+ std::stringstream stream;
+ stream << xOrderR[i];
+ Rotation<double,dim> rotation(rotationMap.at(stream.str()));
+ FieldMatrix<double,dim,dim> rotationMatrix(0);
+ rotation.matrix(rotationMatrix);
+ rot[i].set(rotationMatrix);
+ for (int j = 0; j < 3; j++)
+ rot_director[j][i] = rot[i].director(j);
+ }
+
+ /////////////////////////////////////////////////////////////
+ // STRESS ASSEMBLER
+ /////////////////////////////////////////////////////////////
+ int quadOrder = parameterSet.hasKey("quadOrder") ? parameterSet.get<int>("quadOrder") : 4;
+
+ auto stressAssembler = GFE::SurfaceCosseratStressAssembler<decltype(basisOrderD),decltype(basisOrderR), FieldVector<double,dim>, Rotation<double,dim>>
+ (basisOrderD, basisOrderR);
+
+
+ std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
+ std::shared_ptr<Elasticity::LocalDensity<dim,ValueType>> elasticDensity;
+
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+
+ if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
+ elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "neohooke")
+ elasticDensity = std::make_shared<Elasticity::NeoHookeDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "hencky")
+ elasticDensity = std::make_shared<Elasticity::HenckyDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "exphencky")
+ elasticDensity = std::make_shared<Elasticity::ExpHenckyDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
+ elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(materialParameters);
+
+ if(!elasticDensity)
+ DUNE_THROW(Exception, "Error: Selected energy not available!");
+
+ Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
+ Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
+ Python::Reference pythonObject = surfaceShellCallable();
+ auto fLame = Python::make_function<FieldVector<double, 2> >(pythonObject.get("lame"));
+
+ std::vector<FieldMatrix<double,dim,dim>> stressSubstrate1stPiolaKirchhoffTensor;
+ std::vector<FieldMatrix<double,dim,dim>> stressSubstrateCauchyTensor;
+ std::cout << "Assemble stress for the substrate.." << std::endl;
+ stressAssembler.assembleSubstrateStress<Elasticity::LocalDensity<dim,ValueType>>(x, elasticDensity.get(), quadOrder, stressSubstrate1stPiolaKirchhoffTensor, stressSubstrateCauchyTensor);
+
+ std::vector<FieldMatrix<double,dim,dim>> stressShellBiotTensor;
+ std::cout << "Assemble stress for the shell.." << std::endl;
+ stressAssembler.assembleShellStress(rot, x, xInitial, fLame,/*mu_c*/0, surfaceShellBoundary, quadOrder, stressShellBiotTensor);
+
+ std::vector<double> stressSubstrate1stPiolaKirchhoff(stressSubstrate1stPiolaKirchhoffTensor.size());
+ std::vector<double> stressSubstrateCauchy(stressSubstrate1stPiolaKirchhoffTensor.size());
+ std::vector<double> stressSubstrateVonMises(stressSubstrate1stPiolaKirchhoffTensor.size());
+ std::vector<double> stressShellBiot(stressSubstrate1stPiolaKirchhoffTensor.size());
+
+ for (size_t i = 0; i < stressSubstrate1stPiolaKirchhoffTensor.size(); i++) {
+ stressSubstrate1stPiolaKirchhoff[i] = stressSubstrate1stPiolaKirchhoffTensor[i].frobenius_norm();
+ stressSubstrateCauchy[i] = stressSubstrateCauchyTensor[i].frobenius_norm();
+
+ double vonMises = 0; //von-Mises-Stress
+ for(size_t j=0; j<dim; j++) {
+ int jplus1 = (j+1) % dim;
+ vonMises += 0.5*(stressSubstrate1stPiolaKirchhoffTensor[i][j][j] - stressSubstrate1stPiolaKirchhoffTensor[i][jplus1][jplus1])*(stressSubstrate1stPiolaKirchhoffTensor[i][j][j] - stressSubstrate1stPiolaKirchhoffTensor[i][jplus1][jplus1]);
+ vonMises += 3*stressSubstrate1stPiolaKirchhoffTensor[i][j][jplus1]*stressSubstrate1stPiolaKirchhoffTensor[i][j][jplus1];
+ }
+ stressSubstrateVonMises[i] = std::sqrt(vonMises);
+
+ stressShellBiot[i] = stressShellBiotTensor[i].frobenius_norm();
+ }
+
+
+ auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basisOrderD, displacement);
+
+ auto director0Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basisOrderR, rot_director[0]);
+ auto director1Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basisOrderR, rot_director[1]);
+ auto director2Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basisOrderR, rot_director[2]);
+
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, refinementLevels(displacementOrder-1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ vtkWriter.addVertexData(director0Function, VTK::FieldInfo("director0", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(director1Function, VTK::FieldInfo("director1", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(director2Function, VTK::FieldInfo("director2", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.write("stress_plot_" + parameterSet.get<std::string>("energy"));
+
+ VTKWriter<GridView> vtkWriterElement(gridView);
+ vtkWriterElement.addCellData(stressShellBiot, "stress-shell-element");
+ vtkWriterElement.addCellData(stressShellBiot, "stress-shell-biot");
+ vtkWriterElement.addCellData(stressSubstrate1stPiolaKirchhoff, "stress-substrate-1st-piola-kirchhoff");
+ vtkWriterElement.addCellData(stressSubstrateCauchy, "stress-substrate-cauchy");
+ vtkWriterElement.addCellData(stressSubstrateVonMises, "stress-substrate-von-mises");
+ vtkWriterElement.addVertexData(director0Function, VTK::FieldInfo("director0", VTK::FieldInfo::Type::vector, dim));
+ vtkWriterElement.addVertexData(director1Function, VTK::FieldInfo("director1", VTK::FieldInfo::Type::vector, dim));
+ vtkWriterElement.addVertexData(director2Function, VTK::FieldInfo("director2", VTK::FieldInfo::Type::vector, dim));
+ vtkWriterElement.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ vtkWriterElement.write("stress_plot_" + parameterSet.get<std::string>("energy") + "_element");
+
+ VTKWriter<GridView> vtkWriterElementOnly(gridView);
+ vtkWriterElementOnly.addCellData(stressShellBiot, "stress-shell-biot");
+ vtkWriterElementOnly.addCellData(stressSubstrate1stPiolaKirchhoff, "stress-substrate-1st-piola-kirchhoff");
+ vtkWriterElementOnly.addCellData(stressSubstrateCauchy, "stress-substrate-cauchy");
+ vtkWriterElementOnly.addCellData(stressSubstrateVonMises, "stress-substrate-von-mises");
+ vtkWriterElementOnly.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ vtkWriterElementOnly.write("stress_plot_" + parameterSet.get<std::string>("energy") + "_element_only");
+
+} catch (Exception& e) {
+ std::cout << e.what() << std::endl;
+}