multipatch.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <iganet/splines/boundary.hpp>
 
#include <string_view>
 
namespace iganet {
 
template <typename Patch> class PatchInterface {
public:
  PatchInterface(std::shared_ptr<Patch> firstPatch, enum side firstSide,
                 std::shared_ptr<Patch> secondPatch, enum side secondSide)
      : patches_{std::move(firstPatch), std::move(secondPatch)},
        sides_{firstSide, secondSide} {
    if (!patches_[0] || !patches_[1])
      throw std::invalid_argument("An interface requires two valid patches");
    if (sides_[0] == none || sides_[1] == none)
      throw std::invalid_argument("An interface requires two valid sides");
  }
 
  Patch &patch(std::size_t endpoint) {
    assert(endpoint < patches_.size());
    return *patches_[endpoint];
  }
 
  const Patch &patch(std::size_t endpoint) const {
    assert(endpoint < patches_.size());
    return *patches_[endpoint];
  }
 
  const std::shared_ptr<Patch> &patchPtr(std::size_t endpoint) const {
    assert(endpoint < patches_.size());
    return patches_[endpoint];
  }
 
  enum side side(std::size_t endpoint) const {
    assert(endpoint < sides_.size());
    return sides_[endpoint];
  }
 
  Patch &firstPatch() { return patch(0); }
  const Patch &firstPatch() const { return patch(0); }
  Patch &secondPatch() { return patch(1); }
  const Patch &secondPatch() const { return patch(1); }
  enum side firstSide() const { return side(0); }
  enum side secondSide() const { return side(1); }
 
private:
  std::array<std::shared_ptr<Patch>, 2> patches_;
  std::array<enum side, 2> sides_;
};
 
template <typename Patch> class MultiPatch {
 
public:
  using interface_type = PatchInterface<Patch>;
 
  MultiPatch() = default;
 
  MultiPatch(const MultiPatch &other)
      : patches_(other.patches_), interfaces_(other.interfaces_) {}
 
  MultiPatch(MultiPatch &&other) noexcept {
    patches_.swap(other.patches_);
    interfaces_.swap(other.interfaces_);
  }
 
public:
  auto begin() { return patches_.begin(); }
 
  auto begin() const { return patches_.begin(); }
  auto cbegin() const noexcept { return patches_.cbegin(); }
 
  auto end() { return patches_.end(); }
 
  auto end() const { return patches_.end(); }
  auto cend() const noexcept { return patches_.cend(); }
 
  auto rbegin() { return patches_.rbegin(); }
 
  auto rbegin() const { return patches_.rbegin(); }
  auto crbegin() const noexcept { return patches_.crbegin(); }
 
  auto rend() { return patches_.rend(); }
 
  auto rend() const { return patches_.rend(); }
  auto crend() const noexcept { return patches_.crend(); }
 
public:
  [[nodiscard]] std::size_t npatches() const { return patches_.size(); }
 
  [[nodiscard]] std::size_t ninterfaces() const { return interfaces_.size(); }
 
  [[nodiscard]] std::size_t nboundaries() const { return patches_.size(); }
 
public:
  std::size_t addPatch(std::shared_ptr<Patch> patch) {
    std::size_t index = patches_.size();
    patches_.push_back(patch);
    return index;
  }
 
  std::size_t addPatch(std::unique_ptr<Patch> patch) {
    std::size_t index = patches_.size();
    patches_.push_back(std::shared_ptr<Patch>(std::move(patch)));
    return index;
  }
 
  std::size_t addInterface(std::size_t firstPatch, enum side firstSide,
                           std::size_t secondPatch, enum side secondSide) {
    assert(firstPatch < patches_.size());
    assert(secondPatch < patches_.size());
    return addInterface(patches_[firstPatch], firstSide, patches_[secondPatch],
                        secondSide);
  }
 
  std::size_t addInterface(std::shared_ptr<Patch> firstPatch,
                           enum side firstSide,
                           std::shared_ptr<Patch> secondPatch,
                           enum side secondSide) {
    if (std::find(patches_.begin(), patches_.end(), firstPatch) ==
            patches_.end() ||
        std::find(patches_.begin(), patches_.end(), secondPatch) ==
            patches_.end())
      throw std::invalid_argument(
          "Interface patches must belong to the MultiPatch");
 
    const std::size_t index = interfaces_.size();
    interfaces_.emplace_back(std::move(firstPatch), firstSide,
                             std::move(secondPatch), secondSide);
    return index;
  }
 
  std::size_t addInterface(interface_type patchInterface) {
    return addInterface(patchInterface.patchPtr(0), patchInterface.side(0),
                        patchInterface.patchPtr(1), patchInterface.side(1));
  }
 
  void removeInterface(std::size_t index) {
    assert(index < interfaces_.size());
    interfaces_.erase(interfaces_.begin() + index);
  }
 
  void clear() {
    interfaces_.clear();
    patches_.clear();
  }
 
  Patch &patch(std::size_t index) {
    assert(index < patches_.size());
    return *patches_[index];
  }
 
  const Patch &patch(std::size_t index) const {
    assert(index < patches_.size());
    return *patches_[index];
  }
 
  std::vector<std::shared_ptr<Patch>> &patches() { return patches_; }
  const std::vector<std::shared_ptr<Patch>> &patches() const {
    return patches_;
  }
 
  interface_type &interface(std::size_t index) {
    assert(index < interfaces_.size());
    return interfaces_[index];
  }
 
  const interface_type &interface(std::size_t index) const {
    assert(index < interfaces_.size());
    return interfaces_[index];
  }
 
  std::vector<interface_type> &interfaces() { return interfaces_; }
  const std::vector<interface_type> &interfaces() const { return interfaces_; }
 
  std::size_t findPatchIndex(const Patch &patch) const {
    return findPatchIndex(&patch);
  }
 
  std::size_t findPatchIndex(const Patch *patch) const {
    auto it = std::find_if(
        patches_.begin(), patches_.end(),
        [patch](const auto &candidate) { return candidate.get() == patch; });
    if (it == patches_.end())
      throw std::runtime_error("Did not find the patch index");
 
    return it - patches_.begin();
  }
 
  std::size_t findInterfaceIndex(const interface_type &patchInterface) const {
    return findInterfaceIndex(&patchInterface);
  }
 
  std::size_t findInterfaceIndex(const interface_type *patchInterface) const {
    auto it = std::find_if(interfaces_.begin(), interfaces_.end(),
                           [patchInterface](const auto &candidate) {
                             return &candidate == patchInterface;
                           });
    if (it == interfaces_.end())
      throw std::runtime_error("Did not find the patch interface index");
 
    return it - interfaces_.begin();
  }
 
  template <typename P = Patch>
    requires requires { typename P::value_type; }
  void load(const std::string &filename, const std::string &key = "multipatch",
            Options<typename P::value_type> options =
                Options<typename P::value_type>{}) {
    torch::serialize::InputArchive archive;
    archive.load_from(filename);
    read(archive, key, options);
  }
 
  template <typename P = Patch>
    requires requires { typename P::value_type; }
  torch::serialize::InputArchive &read(torch::serialize::InputArchive &archive,
                                       const std::string &key = "multipatch",
                                       Options<typename P::value_type> options =
                                           Options<typename P::value_type>{}) {
    torch::Tensor data;
    archive.read(key + ".json", data);
    data = data.to(torch::kCPU, torch::kUInt8).contiguous();
    const auto *begin =
        reinterpret_cast<const char *>(data.data_ptr<uint8_t>());
    from_json(nlohmann::json::parse(std::string(begin, begin + data.numel())),
              options);
    return archive;
  }
 
  void save(const std::string &filename,
            const std::string &key = "multipatch") const {
    torch::serialize::OutputArchive archive;
    write(archive, key).save_to(filename);
  }
 
  torch::serialize::OutputArchive &
  write(torch::serialize::OutputArchive &archive,
        const std::string &key = "multipatch") const {
    const std::string serialized = to_json().dump();
    const auto data =
        torch::from_blob(const_cast<char *>(serialized.data()),
                         {static_cast<int64_t>(serialized.size())},
                         torch::TensorOptions{}.dtype(torch::kUInt8))
            .clone();
    archive.write(key + ".json", data);
    return archive;
  }
 
  template <typename P = Patch>
    requires requires { typename P::value_type; }
  bool
  isclose(const MultiPatch &other,
          typename P::value_type rtol = typename P::value_type{1e-5},
          typename P::value_type atol = typename P::value_type{1e-8}) const {
    return jsonIsClose(to_json(), other.to_json(), rtol, atol);
  }
 
  bool operator==(const MultiPatch &other) const {
    return to_json() == other.to_json();
  }
 
  bool operator!=(const MultiPatch &other) const { return !(*this == other); }
 
  void pretty_print(std::ostream &os) const noexcept {
    os << "MultiPatch(\nparDim = " << Patch::parDim()
       << ", npatches = " << npatches() << ", ninterfaces = " << ninterfaces()
       << "\n";
    for (std::size_t patchIndex = 0; patchIndex < patches_.size();
         ++patchIndex) {
      const auto json = patches_[patchIndex]->to_json();
      os << "patch[" << patchIndex << "] = {geoDim = " << json["geoDim"]
         << ", degrees = " << json["degrees"] << "}\n";
    }
    for (std::size_t interfaceIndex = 0; interfaceIndex < interfaces_.size();
         ++interfaceIndex) {
      const auto &patchInterface = interfaces_[interfaceIndex];
      os << "interface[" << interfaceIndex << "] = {patch "
         << findPatchIndex(patchInterface.patchPtr(0).get()) << ", side "
         << static_cast<short_t>(patchInterface.side(0)) << " <-> patch "
         << findPatchIndex(patchInterface.patchPtr(1).get()) << ", side "
         << static_cast<short_t>(patchInterface.side(1)) << "}\n";
    }
    os << ")";
  }
 
  [[nodiscard]] nlohmann::json to_json() const {
    nlohmann::json json;
    json["parDim"] = Patch::parDim();
    json["patches"] = nlohmann::json::array();
    json["interfaces"] = nlohmann::json::array();
    json["boundaries"] = nlohmann::json::array();
 
    for (const auto &patch : patches_)
      json["patches"].push_back(patch->to_json());
 
    for (const auto &patchInterface : interfaces_) {
      nlohmann::json interfaceJson;
      interfaceJson["patches"] = {
          findPatchIndex(patchInterface.patchPtr(0).get()),
          findPatchIndex(patchInterface.patchPtr(1).get())};
      interfaceJson["sides"] = {static_cast<short_t>(patchInterface.side(0)),
                                static_cast<short_t>(patchInterface.side(1))};
 
      const short_t firstAxis =
          (static_cast<short_t>(patchInterface.side(0)) - 1) / 2;
      const short_t secondAxis =
          (static_cast<short_t>(patchInterface.side(1)) - 1) / 2;
      interfaceJson["direction"] = nlohmann::json::array();
      interfaceJson["orientation"] = nlohmann::json::array();
      for (short_t axis = 0; axis < Patch::parDim(); ++axis) {
        short_t mappedAxis = axis;
        if (axis == firstAxis)
          mappedAxis = secondAxis;
        else if (axis == secondAxis)
          mappedAxis = firstAxis;
        interfaceJson["direction"].push_back(mappedAxis);
        interfaceJson["orientation"].push_back(1);
      }
      json["interfaces"].push_back(std::move(interfaceJson));
    }
 
    for (std::size_t patchIndex = 0; patchIndex < patches_.size();
         ++patchIndex) {
      for (short_t patchSide = 1; patchSide <= 2 * Patch::parDim();
           ++patchSide) {
        const bool isInterface =
            std::any_of(interfaces_.begin(), interfaces_.end(),
                        [&](const auto &patchInterface) {
                          return (patchInterface.patchPtr(0).get() ==
                                      patches_[patchIndex].get() &&
                                  patchInterface.side(0) == patchSide) ||
                                 (patchInterface.patchPtr(1).get() ==
                                      patches_[patchIndex].get() &&
                                  patchInterface.side(1) == patchSide);
                        });
        if (!isInterface)
          json["boundaries"].push_back(
              {{"patch", patchIndex}, {"side", patchSide}});
      }
    }
    return json;
  }
 
  template <typename P = Patch>
    requires requires { typename P::value_type; }
  MultiPatch &from_json(const nlohmann::json &json,
                        Options<typename P::value_type> options =
                            Options<typename P::value_type>{}) {
    if (json.at("parDim").get<short_t>() != Patch::parDim())
      throw std::runtime_error(
          "MultiPatch JSON provides an incompatible parametric dimension");
 
    const auto &patchJson = json.at("patches");
    if (!patchJson.is_array() ||
        (!patches_.empty() && patchJson.size() != patches_.size()))
      throw std::runtime_error(
          "MultiPatch JSON patch count does not match the patch container");
 
    auto parsedPatches = patches_;
    const bool createPatches = parsedPatches.empty();
    if (createPatches) {
      parsedPatches.reserve(patchJson.size());
      for (const auto &item : patchJson) {
        try {
          parsedPatches.push_back(
              createUniformBSpline<typename Patch::value_type, Patch::geoDim(),
                                   Patch::parDim()>(item, options));
        } catch (const std::runtime_error &) {
          parsedPatches.push_back(
              createNonUniformBSpline<typename Patch::value_type,
                                      Patch::geoDim(), Patch::parDim()>(
                  item, options));
        }
      }
    }
 
    const auto &interfaceJson = json.at("interfaces");
    if (!interfaceJson.is_array())
      throw std::runtime_error("MultiPatch JSON interfaces must be an array");
 
    std::vector<interface_type> parsedInterfaces;
    parsedInterfaces.reserve(interfaceJson.size());
    for (const auto &item : interfaceJson) {
      const auto patchIndices = item.at("patches").get<std::array<size_t, 2>>();
      const auto sides = item.at("sides").get<std::array<short_t, 2>>();
      if (patchIndices[0] >= parsedPatches.size() ||
          patchIndices[1] >= parsedPatches.size() || sides[0] <= none ||
          sides[0] > 2 * Patch::parDim() || sides[1] <= none ||
          sides[1] > 2 * Patch::parDim())
        throw std::runtime_error("MultiPatch JSON has an invalid interface");
 
      if (!item.at("direction").is_array() ||
          item.at("direction").size() != Patch::parDim() ||
          !item.at("orientation").is_array() ||
          item.at("orientation").size() != Patch::parDim())
        throw std::runtime_error(
            "MultiPatch JSON has invalid interface orientation data");
 
      parsedInterfaces.emplace_back(
          parsedPatches[patchIndices[0]], static_cast<enum side>(sides[0]),
          parsedPatches[patchIndices[1]], static_cast<enum side>(sides[1]));
    }
 
    if (!createPatches)
      for (std::size_t patchIndex = 0; patchIndex < parsedPatches.size();
           ++patchIndex)
        parsedPatches[patchIndex]->from_json(patchJson[patchIndex]);
 
    patches_ = std::move(parsedPatches);
    interfaces_ = std::move(parsedInterfaces);
    return *this;
  }
 
  [[nodiscard]] pugi::xml_document
  to_xml(int id = 0, const std::string &label = "", int index = -1) const {
    pugi::xml_document doc;
    pugi::xml_node root = doc.append_child("xml");
    to_xml(root, id, label, index);
    return doc;
  }
 
  pugi::xml_node &to_xml(pugi::xml_node &root, int id = 0,
                         const std::string &label = "", int index = -1) const {
    for (std::size_t patchIndex = 0; patchIndex < patches_.size(); ++patchIndex)
      patches_[patchIndex]->to_xml(root, static_cast<int>(patchIndex));
 
    pugi::xml_node multiPatch = root.append_child("MultiPatch");
    multiPatch.append_attribute("parDim") = Patch::parDim();
    multiPatch.append_attribute("id") = id;
    if (!label.empty())
      multiPatch.append_attribute("label") = label.c_str();
    if (index >= 0)
      multiPatch.append_attribute("index") = index;
 
    pugi::xml_node patches = multiPatch.append_child("patches");
    patches.append_attribute("type") = "id_range";
    if (!patches_.empty()) {
      const std::string range =
          "0 " + std::to_string(static_cast<int64_t>(patches_.size()) - 1);
      patches.append_child(pugi::node_pcdata).set_value(range.c_str());
    }
 
    std::stringstream interfaceData;
    for (const auto &patchInterface : interfaces_) {
      interfaceData << "\n      ";
      const auto firstPatch = findPatchIndex(patchInterface.patchPtr(0).get());
      const auto secondPatch = findPatchIndex(patchInterface.patchPtr(1).get());
      const short_t firstAxis =
          (static_cast<short_t>(patchInterface.side(0)) - 1) / 2;
      const short_t secondAxis =
          (static_cast<short_t>(patchInterface.side(1)) - 1) / 2;
 
      interfaceData << firstPatch << ' '
                    << static_cast<short_t>(patchInterface.side(0)) << ' '
                    << secondPatch << ' '
                    << static_cast<short_t>(patchInterface.side(1));
 
      for (short_t axis = 0; axis < Patch::parDim(); ++axis) {
        short_t mappedAxis = axis;
        if (axis == firstAxis)
          mappedAxis = secondAxis;
        else if (axis == secondAxis)
          mappedAxis = firstAxis;
        interfaceData << ' ' << mappedAxis;
      }
      for (short_t axis = 0; axis < Patch::parDim(); ++axis)
        interfaceData << " 1";
    }
    if (!interfaces_.empty())
      interfaceData << "\n    ";
    multiPatch.append_child("interfaces")
        .append_child(pugi::node_pcdata)
        .set_value(interfaceData.str().c_str());
 
    std::stringstream boundaryData;
    for (std::size_t patchIndex = 0; patchIndex < patches_.size();
         ++patchIndex) {
      for (short_t patchSide = 1; patchSide <= 2 * Patch::parDim();
           ++patchSide) {
        const bool isInterface =
            std::any_of(interfaces_.begin(), interfaces_.end(),
                        [&](const auto &patchInterface) {
                          return (patchInterface.patchPtr(0).get() ==
                                      patches_[patchIndex].get() &&
                                  patchInterface.side(0) == patchSide) ||
                                 (patchInterface.patchPtr(1).get() ==
                                      patches_[patchIndex].get() &&
                                  patchInterface.side(1) == patchSide);
                        });
        if (!isInterface)
          boundaryData << "\n      " << patchIndex << ' ' << patchSide;
      }
    }
    if (!boundaryData.str().empty())
      boundaryData << "\n    ";
    multiPatch.append_child("boundary")
        .append_child(pugi::node_pcdata)
        .set_value(boundaryData.str().c_str());
 
    return root;
  }
 
  template <typename P = Patch>
    requires requires { typename P::value_type; }
  MultiPatch &from_xml(const pugi::xml_document &doc, int id = 0,
                       const std::string &label = "", int index = -1,
                       Options<typename P::value_type> options =
                           Options<typename P::value_type>{}) {
    return from_xml(doc.child("xml"), id, label, index, options);
  }
 
  template <typename P = Patch>
    requires requires { typename P::value_type; }
  MultiPatch &from_xml(const pugi::xml_node &root, int id = 0,
                       const std::string &label = "", int index = -1,
                       Options<typename P::value_type> options =
                           Options<typename P::value_type>{}) {
    pugi::xml_node multiPatch;
    for (const auto &candidate : root.children("MultiPatch")) {
      if ((id >= 0 ? candidate.attribute("id").as_int() == id : true) &&
          (index >= 0 ? candidate.attribute("index").as_int() == index
                      : true) &&
          (!label.empty() ? candidate.attribute("label").value() == label
                          : true)) {
        multiPatch = candidate;
        break;
      }
    }
    if (!multiPatch)
      throw std::runtime_error("Did not find the MultiPatch XML node");
 
    if (multiPatch.attribute("parDim").as_int() != Patch::parDim())
      throw std::runtime_error(
          "MultiPatch XML provides an incompatible parametric dimension");
 
    const pugi::xml_node patchRange = multiPatch.child("patches");
    if (!patchRange ||
        std::string_view{patchRange.attribute("type").value()} != "id_range")
      throw std::runtime_error("MultiPatch XML has no valid patch ID range");
 
    std::stringstream rangeData(patchRange.child_value());
    int64_t firstPatch = 0;
    int64_t lastPatch = -1;
    if (!(rangeData >> firstPatch >> lastPatch) || firstPatch != 0 ||
        (!patches_.empty() &&
         lastPatch + 1 != static_cast<int64_t>(patches_.size())))
      throw std::runtime_error(
          "MultiPatch XML patch range does not match the patch container");
 
    auto parsedPatches = patches_;
    const bool createPatches = parsedPatches.empty();
    if (createPatches) {
      parsedPatches.reserve(static_cast<std::size_t>(lastPatch + 1));
      for (int64_t patchIndex = 0; patchIndex <= lastPatch; ++patchIndex) {
        try {
          parsedPatches.push_back(
              createUniformBSpline<typename Patch::value_type, Patch::geoDim(),
                                   Patch::parDim()>(root, patchIndex, "", -1,
                                                    options));
        } catch (const std::runtime_error &) {
          parsedPatches.push_back(
              createNonUniformBSpline<typename Patch::value_type,
                                      Patch::geoDim(), Patch::parDim()>(
                  root, patchIndex, "", -1, options));
        }
      }
    }
 
    std::vector<interface_type> parsedInterfaces;
    const pugi::xml_node interfaceNode = multiPatch.child("interfaces");
    if (!interfaceNode)
      throw std::runtime_error("MultiPatch XML has no interfaces node");
 
    std::stringstream interfaceData(interfaceNode.child_value());
    std::string line;
    while (std::getline(interfaceData, line)) {
      std::stringstream item(line);
      std::size_t firstPatchIndex;
      std::size_t secondPatchIndex;
      short_t firstSide;
      short_t secondSide;
      if (!(item >> firstPatchIndex >> firstSide >> secondPatchIndex >>
            secondSide))
        continue;
 
      if (firstPatchIndex >= parsedPatches.size() ||
          secondPatchIndex >= parsedPatches.size() || firstSide <= none ||
          firstSide > 2 * Patch::parDim() || secondSide <= none ||
          secondSide > 2 * Patch::parDim())
        throw std::runtime_error("MultiPatch XML has an invalid interface");
 
      short_t topologyEntry;
      for (short_t entry = 0; entry < 2 * Patch::parDim(); ++entry) {
        if (!(item >> topologyEntry))
          throw std::runtime_error(
              "MultiPatch XML has incomplete interface orientation data");
      }
      if (item >> topologyEntry)
        throw std::runtime_error(
            "MultiPatch XML has excess interface orientation data");
 
      parsedInterfaces.emplace_back(
          parsedPatches[firstPatchIndex], static_cast<enum side>(firstSide),
          parsedPatches[secondPatchIndex], static_cast<enum side>(secondSide));
    }
 
    if (!createPatches)
      for (std::size_t patchIndex = 0; patchIndex < parsedPatches.size();
           ++patchIndex)
        parsedPatches[patchIndex]->from_xml(root, static_cast<int>(patchIndex));
 
    patches_ = std::move(parsedPatches);
    interfaces_ = std::move(parsedInterfaces);
    return *this;
  }
 
private:
  template <typename real_t>
  static bool jsonIsClose(const nlohmann::json &first,
                          const nlohmann::json &second, real_t rtol,
                          real_t atol) {
    if (first.type() != second.type())
      return first.is_number() && second.is_number() &&
             std::abs(first.template get<real_t>() -
                      second.template get<real_t>()) <=
                 atol + rtol * std::abs(second.template get<real_t>());
    if (first.is_number())
      return std::abs(first.template get<real_t>() -
                      second.template get<real_t>()) <=
             atol + rtol * std::abs(second.template get<real_t>());
    if (first.is_array()) {
      if (first.size() != second.size())
        return false;
      for (std::size_t i = 0; i < first.size(); ++i)
        if (!jsonIsClose(first[i], second[i], rtol, atol))
          return false;
      return true;
    }
    if (first.is_object()) {
      if (first.size() != second.size())
        return false;
      for (const auto &[key, value] : first.items()) {
        const auto it = second.find(key);
        if (it == second.end() || !jsonIsClose(value, *it, rtol, atol))
          return false;
      }
      return true;
    }
    return first == second;
  }
 
  std::vector<std::shared_ptr<Patch>> patches_;
 
  std::vector<interface_type> interfaces_;
};
 
template <typename Patch>
std::ostream &operator<<(std::ostream &os, const MultiPatch<Patch> &obj) {
  obj.pretty_print(os);
  return os;
}
 
template <typename Patch>
torch::serialize::OutputArchive &
operator<<(torch::serialize::OutputArchive &archive,
           const MultiPatch<Patch> &obj) {
  return obj.write(archive);
}
 
template <typename Patch>
torch::serialize::InputArchive &
operator>>(torch::serialize::InputArchive &archive, MultiPatch<Patch> &obj) {
  return obj.read(archive);
}
 
} // namespace iganet