BasicPostProc.hpp

Go to the documentation of this file.

/** \file BasicPostProc.hpp
 * \brief Helpers for post-processing basic finite element fields.
 */
 
#ifndef __BASIC_POST_PROC_HPP__
#define __BASIC_POST_PROC_HPP__
 
/**
 * \brief Core broken-mesh post-processor adapted to the fat-prism element.
 *
 * FatPrismElementForcesAndSourcesCore builds its integration points as a
 * Cartesian product of points on a triangle and points through the prism
 * thickness.  The core post-processor stores the same points as one reference
 * prism.  This adapter exposes that reference-prism point set in the form
 * expected by the fat-prism element and keeps the point-to-vertex map in the
 * resulting Cartesian-product order.
 */
struct BasicPostProcFatPrism
    : public PostProcBrokenMeshInMoabBase<
          FatPrismElementForcesAndSourcesCore> {
 
  using Base =
      PostProcBrokenMeshInMoabBase<FatPrismElementForcesAndSourcesCore>;
  using Base::Base;
 
  int getRuleTrianglesOnly(int) override { return -1; }
  int getRuleThroughThickness(int) override { return -1; }
 
  MoFEMErrorCode setGaussPtsTrianglesOnly(int order) override {
    MoFEMFunctionBegin;
 
    CHKERR Base::setGaussPts(order);
 
    constexpr double eps = 100 * std::numeric_limits<double>::epsilon();
    std::vector<std::array<double, 2>> triangle_points;
    std::vector<double> thickness_points;
 
    auto same = [eps](const double a, const double b) {
      return std::abs(a - b) <= eps;
    };
    auto has_triangle_point = [&](const double ksi, const double eta) {
      return std::find_if(
                 triangle_points.begin(), triangle_points.end(),
                 [&](const auto &p) {
                   return same(p[0], ksi) && same(p[1], eta);
                 }) != triangle_points.end();
    };
    auto has_thickness_point = [&](const double zeta) {
      return std::find_if(thickness_points.begin(), thickness_points.end(),
                          [&](const double z) { return same(z, zeta); }) !=
             thickness_points.end();
    };
 
    const auto reference_points = this->gaussPts;
    const auto original_map = this->mapGaussPts;
    for (size_t gg = 0; gg != reference_points.size2(); ++gg) {
      const double ksi = reference_points(0, gg);
      const double eta = reference_points(1, gg);
      const double zeta = reference_points(2, gg);
      if (!has_triangle_point(ksi, eta))
        triangle_points.push_back({ksi, eta});
      if (!has_thickness_point(zeta))
        thickness_points.push_back(zeta);
    }
 
    const size_t expected_points =
        triangle_points.size() * thickness_points.size();
    if (expected_points != reference_points.size2())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Reference prism points do not form a Cartesian product");
 
    this->gaussPtsTrianglesOnly.resize(3, triangle_points.size(), false);
    this->gaussPtsTrianglesOnly.clear();
    for (size_t tt = 0; tt != triangle_points.size(); ++tt) {
      this->gaussPtsTrianglesOnly(0, tt) = triangle_points[tt][0];
      this->gaussPtsTrianglesOnly(1, tt) = triangle_points[tt][1];
    }
 
    this->gaussPtsThroughThickness.resize(2, thickness_points.size(), false);
    this->gaussPtsThroughThickness.clear();
    for (size_t zz = 0; zz != thickness_points.size(); ++zz)
      this->gaussPtsThroughThickness(0, zz) = thickness_points[zz];
 
    this->mapGaussPts.clear();
    this->mapGaussPts.reserve(expected_points);
    for (const auto &p : triangle_points) {
      for (const double zeta : thickness_points) {
        auto gg = size_t{0};
        for (; gg != reference_points.size2(); ++gg) {
          if (same(reference_points(0, gg), p[0]) &&
              same(reference_points(1, gg), p[1]) &&
              same(reference_points(2, gg), zeta))
            break;
        }
        if (gg == reference_points.size2())
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Could not map a fat-prism post-processing point");
        this->mapGaussPts.push_back(original_map[gg]);
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setGaussPtsThroughThickness(int) override {
    MoFEMFunctionBegin;
    if (!this->gaussPtsThroughThickness.size2())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Fat-prism post-processing points are not initialised");
    MoFEMFunctionReturn(0);
  }
};
 
#endif // __BASIC_POST_PROC_HPP__