hdf5-writer.hh

#ifndef SRC_HDF5_WRITER_HH
#define SRC_HDF5_WRITER_HH

#include <vector>

#include "hdf5/iteration-writer.hh"
#include "hdf5/time-writer.hh"
#include "hdf5-bodywriter.hh"

#include <dune/fufem/hdf5/file.hh>

template <class ProgramState, class VertexBasis, class GridView>
class HDF5Writer {
public:
    using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
    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>;

    HDF5Writer(HDF5::File& file,
             const VertexCoordinates& vertexCoordinates,
             const VertexBases& vertexBases,
             const FrictionPatches& frictionPatches)
             //const WeakPatches& weakPatches)
      : file_(file),
        frictionPatches_(frictionPatches),
        iterationWriter_(file_),
        timeWriter_(file_) {

        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]));
        }
    }

    template <class ContactNetwork, class Friction>
    void reportSolution(const ProgramState& programState, const ContactNetwork& contactNetwork, const Friction& friction) {

        timeWriter_.write(programState);

        friction.updateAlpha(programState.alpha);

        // extract relative velocities
        using Vector = typename ProgramState::Vector;
        Vector mortarV;
        contactNetwork.nBodyAssembler().nodalToTransformed(programState.v, mortarV);

        std::vector<Vector> v_rel;
        split(mortarV, v_rel);

        using ScalarVector = typename ProgramState::ScalarVector;
        const auto frictionCoeff = friction.coefficientOfFriction(mortarV);

        /*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::vector<ScalarVector> splitCoeff;
        split(frictionCoeff, splitCoeff);

        /*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();
            bodyWriters_[i]->reportSolution(programState, v_rel[bodyID], splitCoeff[bodyID]);
        }
    }

    void reportIterations(const ProgramState& programState, const IterationRegister& iterationCount) {
        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 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();

            auto& splitV_ = splitV[bodyIdx];
            splitV_.resize(bodyNodes.size());

            for (size_t i=0; i<splitV_.size(); i++) {
                if (toBool(bodyNodes[i])) {
                    splitV_[i] = v[globalIdx];
                }
                globalIdx++;
            }
        }
    }

private:
  HDF5::File& file_;
  const FrictionPatches& frictionPatches_;

  IterationWriter iterationWriter_;
  TimeWriter<ProgramState> timeWriter_;

  std::vector<std::unique_ptr<HDF5BodyWriter>> bodyWriters_;
};
#endif