build file and writer objects only once, report weightedNormalStress possible now

dune/tectonic/io/hdf5-bodywriter.hh

@@ -52,6 +52,10 @@ class HDF5BodyWriter {
         #endif*/
     }
 
+    void reportWeightedNormalStress(const ProgramState& programState) {
+        frictionBoundaryWriter_->writeWeightedNormalStress(programState.timeStep, programState.weightedNormalStress[id_], programState.weights[id_]);
+    }
+
     auto id() const {
         return id_;
     }

dune/tectonic/io/hdf5-writer.hh

@@ -16,6 +16,7 @@ class HDF5Writer {
     using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
     using VertexBases = std::vector<const VertexBasis*>;
     using FrictionPatches = std::vector<const typename HDF5BodyWriter::Patch*>;
+    using ScalarVector = typename ProgramState::ScalarVector;
     //using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
 
     //friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
@@ -26,11 +27,11 @@ class HDF5Writer {
              const FrictionPatches& frictionPatches)
              //const WeakPatches& weakPatches)
       : file_(file),
+        frictionPatches_(frictionPatches),
         iterationWriter_(file_),
-        timeWriter_(file_),
-        frictionPatches_(frictionPatches) {
+        timeWriter_(file_) {
 
-        for (size_t i=0; i<frictionPatches_.size(); i++) {
+        for (size_t i=0; i<vertexCoordinates.size(); i++) {
             if (frictionPatches_[i]->numVertices() > 0)
                 bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i])); //, weakPatches[i]));
         }
@@ -41,7 +42,7 @@ class HDF5Writer {
 
         timeWriter_.write(programState);
 
-        //friction.updateAlpha(programState.alpha);
+        friction.updateAlpha(programState.alpha);
 
         // extract relative velocities
         using Vector = typename ProgramState::Vector;
@@ -54,17 +55,23 @@ class HDF5Writer {
         using ScalarVector = typename ProgramState::ScalarVector;
         const auto frictionCoeff = friction.coefficientOfFriction(mortarV);
 
-        double norm = 0;
+        /*double norm = 0;
         const auto& bodyNodes = *frictionPatches_[0]->getVertices();
         for (size_t i=bodyNodes.size(); i<frictionCoeff.size(); i++) {
             norm += frictionCoeff[i].two_norm();
-        }
-        std::cout << std::setprecision(10) << "friction coefficients norm: " << norm << std::endl;
+        } */
+        //std::cout << std::setprecision(10) << "friction coefficients norm: " << norm << std::endl;
 
         std::vector<ScalarVector> splitCoeff;
         split(frictionCoeff, splitCoeff);
 
-        //print(v_rel, "v_rel: ");
+        /*if (std::isnan(norm) or std::isinf(norm)) {
+            print(programState.alpha, "alpha");
+            print(frictionCoeff, "frictionCoeff");
+            print(splitCoeff, "splitCoeff");
+            print(v_rel, "v_rel: ");
+            DUNE_THROW(Dune::Exception, "invalid state");
+        }*/
 
         for (size_t i=0; i<bodyWriters_.size(); i++) {
             auto bodyID = bodyWriters_[i]->id();
@@ -76,13 +83,21 @@ class HDF5Writer {
         iterationWriter_.write(programState.timeStep, iterationCount);
     }
 
+
+    void reportWeightedNormalStress(const ProgramState& programState) {
+        for (size_t i=0; i<bodyWriters_.size(); i++) {
+            auto bodyID = bodyWriters_[i]->id();
+            bodyWriters_[i]->reportWeightedNormalStress(programState);
+        }
+    }
+
     template <class VectorType>
     void split(const VectorType& v, std::vector<VectorType>& splitV) const {
-
-        const size_t nBodies = frictionPatches_.size();
+        const auto nBodies = frictionPatches_.size();
 
         size_t globalIdx = 0;
         splitV.resize(nBodies);
+
         for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
             const auto& bodyNodes = *frictionPatches_[bodyIdx]->getVertices();
 
@@ -100,11 +115,11 @@ class HDF5Writer {
 
 private:
   HDF5::File& file_;
+  const FrictionPatches& frictionPatches_;
 
   IterationWriter iterationWriter_;
   TimeWriter<ProgramState> timeWriter_;
 
-  const FrictionPatches& frictionPatches_;
   std::vector<std::unique_ptr<HDF5BodyWriter>> bodyWriters_;
 };
 #endif

dune/tectonic/io/hdf5/frictionalboundary-writer.hh

@@ -25,10 +25,15 @@ template <class GridView> class FrictionalBoundaryWriter {
       frictionalBoundaryStateWriter_(group_, "state",
                                      frictionalBoundary.numVertices()),
       frictionalBoundaryCoefficientWriter_(group_, "coefficient",
-                                           frictionalBoundary.numVertices()) {
+                                           frictionalBoundary.numVertices()),
+      frictionalBoundaryWeightsWriter_(group_, "weights",
+                                           frictionalBoundary.numVertices()),
+      frictionalBoundaryWeightedNormalStressWriter_(group_, "weightedNormalStress",
+                                           frictionalBoundary.numVertices()){
 
   auto const frictionalBoundaryCoordinates =
       restrictToSurface(vertexCoordinates, frictionalBoundary);
+
   HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
       group_, "coordinates", frictionalBoundaryCoordinates.size(),
       Vector::block_type::dimension);
@@ -51,6 +56,16 @@ template <class GridView> class FrictionalBoundaryWriter {
       addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
     }
 
+  template <class ScalarVector>
+  void writeWeightedNormalStress(const size_t timeStep, const ScalarVector& weightedNormalStress, const ScalarVector& weights) {
+
+      auto const frictionalBoundaryWeights = restrictToSurface(weights, frictionalBoundary_);
+      addEntry(frictionalBoundaryWeightsWriter_, timeStep, frictionalBoundaryWeights);
+
+      auto const frictionalBoundaryWeightedNormalStress = restrictToSurface(weightedNormalStress, frictionalBoundary_);
+      addEntry(frictionalBoundaryWeightedNormalStressWriter_, timeStep, frictionalBoundaryWeightedNormalStress);
+  }
+
 private:
   HDF5::Group& group_;
 
@@ -60,6 +75,8 @@ template <class GridView> class FrictionalBoundaryWriter {
   HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
   HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
   HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryWeightsWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryWeightedNormalStressWriter_;
 };
 
 #endif

dune/tectonic/io/io-handler.hh

@@ -2,6 +2,7 @@
 #define SRC_IO_HANDLER_HH
 
 #include <filesystem>
+#include <memory>
 
 #include <dune/common/parametertree.hh>
 
@@ -24,7 +25,7 @@
 #include "hdf5-writer.hh"
 #include "vtk.hh"
 
-template <class Assembler, class ContactNetwork>
+template <class Assembler, class ContactNetwork, class ProgramState>
 class IOHandler {
 private:
     using Vector = typename Assembler::Vector;
@@ -32,6 +33,7 @@ class IOHandler {
 
     using MyVertexBasis = typename Assembler::VertexBasis;
     using MyCellBasis = typename Assembler::CellBasis;
+    using HDF5Writer = HDF5Writer<ProgramState, MyVertexBasis, typename Assembler::GV>;
 
     const size_t bodyCount_;
     std::vector<size_t> nVertices_;
@@ -39,8 +41,12 @@ class IOHandler {
     std::vector<Vector> vertexCoordinates_;
     std::vector<const MyVertexBasis* > vertexBases_;
     std::vector<const MyCellBasis* > cellBases_;
+    std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches_;
     std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches_;
 
+    std::unique_ptr<HDF5::File> dataFile_;
+    std::unique_ptr<HDF5Writer> dataWriter_;
+
     bool writeVTK_;
     bool writeData_;
 
@@ -88,9 +94,15 @@ class IOHandler {
             if (!std::filesystem::is_directory("iterates/"))
                 std::filesystem::create_directory("iterates");
         }
+
+
+        contactNetwork.frictionPatches(frictionPatches_);
+        dataFile_ = std::make_unique<HDF5::File>("output.h5");
+        dataWriter_ = std::make_unique<HDF5Writer>(*dataFile_, vertexCoordinates_, vertexBases_,
+                              frictionPatches_);
     }
 
-    template <class ProgramState, class GlobalFriction>
+    template <class GlobalFriction>
     void write(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial = false) {
       if (writeData_) {
         writeData(programState, contactNetwork, friction, iterationCount, initial);
@@ -112,7 +124,6 @@ class IOHandler {
       }
     }
 
-    template <class ProgramState>
     bool read(ProgramState& programState) {
         using ADS = RestartIO<ProgramState>;
 
@@ -128,7 +139,6 @@ class IOHandler {
     }
 
 private:
-    template <class ProgramState>
     void writeVTK(const ProgramState& programState, const ContactNetwork& contactNetwork) {
         using ScalarVector = typename Assembler::ScalarVector;
         std::vector<ScalarVector> stress(bodyCount_);
@@ -146,7 +156,6 @@ class IOHandler {
                         programState.alpha, stress);
     }
 
-    template <class ProgramState>
     void writeRestarts(const ProgramState& programState) {
         if (programState.timeStep % restartSpacing_ == 0) {
           auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READWRITE);
@@ -157,21 +166,16 @@ class IOHandler {
         }
     }
 
-    template <class ProgramState, class GlobalFriction>
+    template <class GlobalFriction>
     void writeData(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial) {
-      std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
-      contactNetwork.frictionPatches(frictionPatches);
-
-      auto dataFile = std::make_unique<HDF5::File>("output.h5");
-      auto dataWriter = std::make_unique<
-                            HDF5Writer<ProgramState, MyVertexBasis, typename Assembler::GV>>(
-                            *dataFile, vertexCoordinates_, vertexBases_,
-                            frictionPatches);
-
-      dataWriter->reportSolution(programState, contactNetwork, friction);
-      if (!initial)
-        dataWriter->reportIterations(programState, iterationCount);
-      dataFile->flush();
+      dataWriter_->reportSolution(programState, contactNetwork, friction);
+      if (!initial) {
+        dataWriter_->reportIterations(programState, iterationCount);
+      } else {
+        dataWriter_->reportWeightedNormalStress(programState);
+      }
+
+      dataFile_->flush();
     }
 
 };