Class to compute and write von mises stress of given displacement or series

dune/elasticity/CMakeLists.txt

@@ -2,3 +2,4 @@ add_subdirectory(assemblers)
 add_subdirectory(common)
 add_subdirectory(estimators)
 add_subdirectory(materials)
+add_subdirectory(utilities)

dune/elasticity/utilities/CMakeLists.txt

+install(FILES
+    stresswriter.hh
+    DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/utilities)

dune/elasticity/utilities/stresswriter.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+#ifndef DUNE_ELASTICITY_UTILITIES_STRESS_WRITER_HH
+#define DUNE_ELASTICITY_UTILITIES_STRESS_WRITER_HH
+
+#include <string>
+#include <vector>
+#include <iomanip>
+#include <algorithm>
+
+#include <dune/common/parametertree.hh>
+#include <dune/grid/io/file/amirameshwriter.hh>
+#include <dune/grid/io/file/amirameshreader.hh>
+
+#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
+#include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p0basis.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
+
+/** \brief Methods to compute von mises stresses and write them to the file system.
+ *
+ * The parameter set the user provides, needs to contain the following flags:
+ *  - "material"  The material type of the body, either "linear" for linear elasticity or
+ *                "geomExactLinear" for the geometrically exact St. Venant-Kirchhof material are supported yet.
+ *  - the material parameter for your material, like "E" and "nu"
+ *  - "format"    The format the stress should be written in. By now only "amira" is supported.
+ *  - "nodalStress" Determines whether the stress is computed as element data or vertex data.
+ *
+ * For writing the stress one additionally needs to specify
+ *  - "resultFile" The path + filename where the stress should be saved
+ *
+ * For computing and writing the stresses of a time series, previously written to the file system you need
+ *
+ * - "pathWithPrefix" The path where the stresses should be written to including the naming for the displacement fields
+ * - "numberTimeSteps" The number of time steps
+ * - "stressName" A name for the stresses to be written
+ *
+ */
+template <class DisplacementField, class Grid>
+struct StressWriter
+{
+  using field_type = typename Dune::FieldTraits<DisplacementField>::field_type;
+  using StressBasis = P0Basis<typename Grid::LeafGridView, field_type>;
+  using LocalStressAssembler = VonMisesStressAssembler<Grid, typename StressBasis::LocalFiniteElement>;
+  using StressFunction = typename LocalStressAssembler::StressFunction;
+
+  using P1Basis = P1NodalBasis<typename Grid::LeafGridView, field_type>;
+  using StressVector = typename LocalStressAssembler::LocalVector;
+
+  static auto getStress(const Dune::ParameterTree &config, const DisplacementField& displace, const Grid& grid);
+
+  static void writeStressToFile(const Dune::ParameterTree& config, const DisplacementField& displace, const Grid& grid);
+
+  static void writeTimeSeriesToFile(const Dune::ParameterTree& config, const Grid& grid);
+
+private:
+  static StressFunction getStressFunction(const Dune::ParameterTree& config) {
+
+    std::string material = config.get<std::string>("material");
+    if (material == "linear")
+      return LocalStressAssembler::linearElasticStressFunction(config.get<field_type>("E"), config.get<field_type>("nu"));
+    else if (material == "geomExactLinear")
+      return GeomExactStVenantMaterial<P1Basis>::stressFunction(config.get<field_type>("E"), config.get<field_type>("nu"));
+    else
+      DUNE_THROW(Dune::NotImplemented, "Stress function for " + material + " not implemented yet");
+  }
+};
+
+template <class DisplacementField, class Grid>
+auto StressWriter<DisplacementField, Grid>::getStress(const Dune::ParameterTree& config, const DisplacementField& displace, const Grid& grid)
+{
+  P1Basis basis(grid.leafGridView());
+  auto displacementFunction = std::make_shared<BasisGridFunction<P1Basis, DisplacementField> >(basis, displace);
+
+  LocalStressAssembler vonMisesAssembler(displacementFunction, getStressFunction(config));
+
+  StressBasis stressBasis(grid.leafGridView());
+  FunctionalAssembler<StressBasis> globalAssembler(stressBasis);
+
+  StressVector stress;
+  globalAssembler.assemble(vonMisesAssembler, stress);
+
+  if (config.get<bool>("nodalStress")) {
+    const auto& indexSet = grid.leafIndexSet();
+
+    StressVector nodalStress(basis.size());
+    nodalStress = 0;
+    std::vector<size_t> neighboringElementsPerVertex(nodalStress.size(), 0);
+
+    for (const auto& e : elements(grid.leafGridView()))
+      for (size_t j = 0; j < e.geometry().corners(); j++) {
+        auto vertexIndex = indexSet.subIndex(e, j, Grid::dimension);
+        neighboringElementsPerVertex[vertexIndex]++;
+        nodalStress[vertexIndex] += stress[indexSet.index(e)];
+      }
+
+    for (size_t i = 0; i < nodalStress.size(); i++)
+      nodalStress[i] /= neighboringElementsPerVertex[i];
+
+    return nodalStress;
+  } else
+    return stress;
+}
+
+template <class DisplacementField, class Grid>
+void StressWriter<DisplacementField, Grid>::writeStressToFile(const Dune::ParameterTree& config, const DisplacementField& displace,
+                                                        const Grid& grid)
+{
+  auto stress = getStress(config, displace, grid);
+
+  std::string format = config.get<std::string>("format");
+  std::transform(std::begin(format), std::end(format), std::begin(format), ::tolower);
+
+  if (format == "amira") {
+#if HAVE_AMIRAMESH
+      Dune::LeafAmiraMeshWriter<Grid>::writeBlockVector(grid, stress, config.get<std::string>("resultFile"));
+#endif
+  } else
+      DUNE_THROW(Dune::NotImplemented, "Format " + format + " is not yet implemented");
+}
+
+template <class DisplacementField, class Grid>
+void StressWriter<DisplacementField, Grid>::writeTimeSeriesToFile(const Dune::ParameterTree& config, const Grid& grid)
+{
+
+  std::string pathWithPrefix = config.get<std::string>("pathWithPrefix");
+  std::string stressName = config.get<std::string>("stressName");
+
+  size_t numberTimeSteps = config.get<size_t>("numberTimeSteps");
+  for (size_t i = 0; i < numberTimeSteps; i++) {
+    std::stringstream step;
+    step << std::setw(2) << std::setfill('0') << i;
+
+#if HAVE_AMIRAMESH
+    DisplacementField displace;
+    Dune::AmiraMeshReader<Grid>::readFunction(displace, pathWithPrefix + step.str());
+    auto stress = getStress(config, displace, grid);
+    Dune::LeafAmiraMeshWriter<Grid>::writeBlockVector(grid, stress, pathWithPrefix + stressName + step.str());
+#endif
+  }
+}
+
+#endif