patchinfo-writer.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/fufem/grid/hierarchic-approximation.hh>
#include <dune/fufem/hdf5/singletonwriter.hh>

#include "patchinfo-writer.hh"

template <class LocalVector, class GridView>
GridEvaluator<LocalVector, GridView>::GridEvaluator(const CuboidGeometry& cuboidGeometry,
    ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
  double const bufferWidth = 0.05 * cuboidGeometry.lengthScale();
  auto const xminmax = std::minmax_element(
      weakPatch.vertices.begin(), weakPatch.vertices.end(),
      [](LocalVector const &a, LocalVector const &b) { return a[0] < b[0]; });
  double const xmin = (*xminmax.first)[0] - bufferWidth;
  double const xspan = (*xminmax.second)[0] + bufferWidth - xmin;

  double const zmin = -cuboidGeometry.depth() / 2.0;
  double const zspan = cuboidGeometry.depth();

  double spatialResolution = 0.01 * cuboidGeometry.lengthScale();
  size_t const xsteps = std::round(xspan / spatialResolution);
  size_t const zsteps = std::round(zspan / spatialResolution);

  xCoordinates.resize(xsteps + 1);
  for (size_t xi = 0; xi <= xsteps; xi++)
    xCoordinates[xi] = xmin + xi * xspan / xsteps;
  zCoordinates.resize(zsteps + 1);
  for (size_t zi = 0; zi <= zsteps; zi++)
    zCoordinates[zi] = zmin + zi * zspan / zsteps;

  HierarchicApproximation<typename GridView::Grid, GridView> const
      hApproximation(gridView.grid(), gridView, 1e-6 * cuboidGeometry.lengthScale());

  LocalVector global(0);
  localInfo.resize(xsteps + 1);
  for (size_t xi = 0; xi < xCoordinates.size(); ++xi) {
    localInfo[xi].resize(zsteps + 1);
    for (size_t zi = 0; zi < zCoordinates.size(); ++zi) {
      global[0] = xCoordinates[xi];
      global[2] = zCoordinates[zi];
      auto const element = hApproximation.findEntity(global);
      localInfo[xi][zi] =
          std::make_pair(element, element.geometry().local(global));
    }
  }
}

template <class LocalVector, class GridView>
template <class Function>
Dune::Matrix<typename Function::RangeType>
GridEvaluator<LocalVector, GridView>::evaluate(Function const &f) const {
  Dune::Matrix<typename Function::RangeType> ret(xCoordinates.size(),
                                                 zCoordinates.size());
  for (size_t xi = 0; xi < localInfo.size(); ++xi) {
    auto const &localInfoX = localInfo[xi];
    for (size_t zi = 0; zi < localInfoX.size(); ++zi) {
      auto const &localInfoXZ = localInfoX[zi];
      f.evaluateLocal(localInfoXZ.first, localInfoXZ.second, ret[xi][zi]);
    }
  }
  return ret;
}

template <class ProgramState, class VertexBasis, class GridView>
PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
    HDF5::Grouplike &file, VertexBasis const &vertexBasis,
    Patch const &frictionalBoundary,
    ConvexPolyhedron<LocalVector> const &weakPatch,
    size_t const id)
    : id_(id),
      group_(file, "weakPatchGrid"+std::to_string(id_)),
      vertexBasis_(vertexBasis),
      gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
      weakPatchGridVelocityWriter_(
          group_, "velocity", gridEvaluator_.xCoordinates.size(),
          gridEvaluator_.zCoordinates.size(), Vector::block_type::dimension) {
  HDF5::SingletonWriter<1> weakPatchGridXCoordinateWriter(
      group_, "xCoordinates", gridEvaluator_.xCoordinates.size());
  setEntry(weakPatchGridXCoordinateWriter, gridEvaluator_.xCoordinates);

  HDF5::SingletonWriter<1> weakPatchGridZCoordinateWriter(
      group_, "zCoordinates", gridEvaluator_.zCoordinates.size());
  setEntry(weakPatchGridZCoordinateWriter, gridEvaluator_.zCoordinates);
}

template <class ProgramState, class VertexBasis, class GridView>
void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
    ProgramState const &programState) {
    // TODO
  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]);
  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
}

#include "patchinfo-writer_tmpl.cc"