igabase.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <utils/tuple.hpp>
 
#include <filesystem>
 
namespace iganet {
 
enum class collPts : short_t {
  greville = 0,          
  greville_interior = 1, 
  greville_ref1 = 2,     
  greville_interior_ref1 =
      3,             
  greville_ref2 = 4, 
  greville_interior_ref2 =
      5, 
};
 
namespace detail {
  template <typename T>
  concept HasAsTensor = requires(T a) {
    { a.as_tensor() };
  };
} 
 
template <typename, typename, typename = void>
class IgABase2;
  
template <detail::HasAsTensor... Inputs,
          detail::HasAsTensor... Outputs,
          detail::HasAsTensor... CollPts>
class IgABase2<std::tuple<Inputs...>,
               std::tuple<Outputs...>,
               std::tuple<CollPts...>> {
public:
  using value_type =
    typename std::common_type<typename Inputs::value_type..., typename Outputs::value_type...>::type;
 
  using inputs_type = std::tuple<Inputs...>;
 
  using outputs_type = std::tuple<Outputs...>;
 
  using collPts_type =
    std::tuple<std::pair<typename CollPts::eval_type,
                         typename CollPts::boundary_eval_type>...>;
 
protected:
  inputs_type inputs_;
 
  outputs_type outputs_;
 
  collPts_type collPts_;
 
private:
  template <typename... Objs, std::size_t... NumCoeffs>
  auto construct_tuple_from_arrays(const std::tuple<std::array<int64_t, NumCoeffs>...>& numCoeffs,
                                   enum init init,
                                   iganet::Options<value_type> options) {
    return std::make_tuple(
                           std::apply(
                                      [&](auto&&... args) { return Objs(std::forward<decltype(args)>(args)..., init, options); },
                                      numCoeffs
                                      )...
                           );
  }
 
  template <typename... Objs, typename... NumCoeffsTuples, std::size_t... Is>
  auto construct_tuple_from_arrays_impl(const std::tuple<NumCoeffsTuples...>& numCoeffs,
                                        enum init init,
                                        iganet::Options<value_type> options,
                                        std::index_sequence<Is...>) {
    static_assert(sizeof...(Objs) == sizeof...(NumCoeffsTuples));
    return std::make_tuple(
                           std::apply(
                                      [&](auto&&... args) { return Objs(std::forward<decltype(args)>(args)..., init, options); },
                                      std::get<Is>(numCoeffs)
                                      )...
                           );
  }
    
  template <typename... Objs, typename... NumCoeffsTuples>
  auto construct_tuple_from_tuples(const std::tuple<NumCoeffsTuples...>& numCoeffs,
                                   enum init init,
                                   iganet::Options<value_type> options) {
    return construct_tuple_from_tuples_impl<Objs...>(numCoeffs,
                                                     init,
                                                     options,
                                                     std::index_sequence_for<Objs...>{});
  }
 
public:
  explicit IgABase2(
      iganet::Options<value_type> options = iganet::Options<value_type>{})
    : inputs_(), outputs_(), collPts_() {}
 
  template <std::size_t NumCoeffs>
  IgABase2(const std::array<int64_t, NumCoeffs>& ncoeffs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : IgABase2(std::tuple{ncoeffs}, std::tuple{ncoeffs}, std::tuple{ncoeffs}, init, options) {}
 
  template <std::size_t NumCoeffsInputs, std::size_t NumCoeffsOutputs, std::size_t NumCoeffsCollPts>
  IgABase2(const std::array<int64_t, NumCoeffsInputs>& ncoeffsInputs,
           const std::array<int64_t, NumCoeffsOutputs>& ncoeffsOutputs,
           const std::array<int64_t, NumCoeffsCollPts>& ncoeffsCollPts,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : IgABase2(std::tuple{ncoeffsInputs}, std::tuple{ncoeffsOutputs},
               std::tuple{ncoeffsCollPts}, init, options) {}
 
  template <std::size_t... NumCoeffs>
  IgABase2(const std::tuple<std::array<int64_t, NumCoeffs>...>& ncoeffs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : IgABase2(ncoeffs, ncoeffs, ncoeffs, init, options)
  {}
  
  template <std::size_t... NumCoeffsInputs, std::size_t... NumCoeffsOutputs, std::size_t... NumCoeffsCollPts>
  IgABase2(const std::tuple<std::array<int64_t, NumCoeffsInputs>...>& ncoeffsInputs,
           const std::tuple<std::array<int64_t, NumCoeffsOutputs>...>& ncoeffsOutputs,
           const std::tuple<std::array<int64_t, NumCoeffsCollPts>...>& ncoeffsCollPts,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : inputs_(construct_tuple_from_arrays<Inputs...>(ncoeffsInputs, init, options)),
      outputs_(construct_tuple_from_arrays<Outputs...>(ncoeffsOutputs, init, options)),
      collPts_(construct_tuple_from_arrays<Outputs...>(ncoeffsCollPts, init, options))
  {}
 
  template<typename... CoeffsInputs, typename... CoeffsOutputs, typename... CoeffsCollPts>
  IgABase2(const std::tuple<CoeffsInputs...>& coeffsInputs,
           const std::tuple<CoeffsOutputs...>& coeffsOutputs,
           const std::tuple<CoeffsCollPts...>& coeffsCollPts,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : inputs_(construct_tuple_from_tuples<Inputs...>(coeffsInputs, init, options)),
      outputs_(construct_tuple_from_tuples<Outputs...>(coeffsOutputs, init, options)),
      collPts_(construct_tuple_from_tuples<CollPts...>(coeffsCollPts, init, options))
  {}
 
  inline constexpr const auto &inputs() const {
    return inputs_;
  }
 
  inline constexpr auto &inputs() { 
    return inputs_;
  }
  
  template <std::size_t index>
  inline constexpr const auto &input() const {
    static_assert(index >= 0 && index < sizeof...(Inputs));
    return std::get<index>(inputs_);
  }
 
  template <std::size_t index>
  inline constexpr auto &input() { 
    static_assert(index >= 0 && index < sizeof...(Inputs));
    return std::get<index>(inputs_);
  }
 
  inline constexpr const auto &outputs() const {
    return outputs_;
  }
 
  inline constexpr auto &outputs() { 
    return outputs_;
  }
  
  template <std::size_t index>
  inline constexpr const auto &output() const {
    static_assert(index >= 0 && index < sizeof...(Outputs));
    return std::get<index>(outputs_);
  }
 
  template <std::size_t index>
  inline constexpr auto &output() { 
    static_assert(index >= 0 && index < sizeof...(Outputs));
    return std::get<index>(outputs_);
  }
 
  inline constexpr const auto &collPts() const {
    return collPts_;
  }
 
  inline constexpr auto &collPts() { 
    return collPts_;
  }
  
private:
  template <std::size_t index, std::size_t... Is>
  std::tuple_element_t<index, collPts_type>
  collPts(enum collPts collPtsType, std::index_sequence<Is...>) const {
    std::tuple_element_t<index, collPts_type> collPts;
 
    switch (collPtsType) {
 
    case collPts::greville:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
        std::get<index>(collPts_).template space<Is>().greville(/* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(collPts_).template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_interior:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            std::get<index>(collPts_).template space<Is>().greville(/* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(collPts_).template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_ref1:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            std::get<index>(collPts_).template space<Is>().clone().uniform_refine().greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            std::get<index>(collPts_).template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref1:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            std::get<index>(collPts_).template space<Is>().clone().uniform_refine().greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            std::get<index>(collPts_).template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_ref2:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            std::get<index>(collPts_).template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(collPts_).template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref2:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            std::get<index>(collPts_).template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(collPts_).template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    default:
      throw std::runtime_error("Invalid collocation point specifier");
    }
 
    return collPts;
  }  
 
public:
  template<std::size_t index>
  std::tuple_element_t<index, collPts_type>
  collPts(enum collPts collPts) const {
    if constexpr (std::tuple_element_t<index, inputs_type>::nspaces() == 1)
 
      switch (collPts) {
 
      case collPts::greville:
        return {std::get<index>(collPts_).space().greville(/* interior */ false),
                std::get<index>(collPts_).boundary().greville()};
 
      case collPts::greville_interior:
        return {std::get<index>(collPts_).space().greville(/* interior */ true),
                std::get<index>(collPts_).boundary().greville()};
 
      case collPts::greville_ref1:
        return {
            std::get<index>(collPts_).space().clone().uniform_refine().greville(/* interior */ false),
            std::get<index>(collPts_).boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_interior_ref1:
        return {
            std::get<index>(collPts_).space().clone().uniform_refine().greville(/* interior */ true),
            std::get<index>(collPts_).boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_ref2:
        return {std::get<index>(collPts_).space().clone().uniform_refine(2, -1).greville(
                    /* interior */ false),
                std::get<index>(collPts_).boundary().clone().uniform_refine(2, -1).greville()};
 
      case collPts::greville_interior_ref2:
        return {std::get<index>(collPts_).space().clone().uniform_refine(2, -1).greville(
                    /* interior */ true),
                std::get<index>(collPts_).boundary().clone().uniform_refine(2, -1).greville()};
 
      default:
        throw std::runtime_error("Invalid collocation point specifier");
      }
 
    else
      return collPts(
                     collPts, std::make_index_sequence<std::tuple_element_t<index, inputs_type>::nspaces()>{});
  }
  
};
 
template <detail::HasAsTensor... Inputs,
          detail::HasAsTensor... Outputs>
class IgABase2<std::tuple<Inputs...>,
               std::tuple<Outputs...>, void> {
public:
  using value_type =
    typename std::common_type<typename Inputs::value_type..., typename Outputs::value_type...>::type;
 
  using inputs_type = std::tuple<Inputs...>;
 
  template<std::size_t index>
  using input_t = typename std::tuple_element_t<index, inputs_type>;
 
  using outputs_type = std::tuple<Outputs...>;
 
  template<std::size_t index>
  using output_t = typename std::tuple_element_t<index, outputs_type>;
  
  using collPts_type =
    std::tuple<std::pair<typename Outputs::eval_type,
                         typename Outputs::boundary_eval_type>...>;
 
  template<std::size_t index>
  using collPts_t = typename std::tuple_element_t<index, collPts_type>;
  
protected:
  inputs_type inputs_;
 
  outputs_type outputs_;
 
private:
  template <typename... Objs, std::size_t... NumCoeffs>
  auto construct_tuple_from_arrays(const std::tuple<std::array<int64_t, NumCoeffs>...>& numCoeffs,
                                   enum init init,
                                   iganet::Options<value_type> options) {
    return std::make_tuple(
                           std::apply(
                                      [&](auto&&... args) { return Objs(std::forward<decltype(args)>(args)..., init, options); },
                                      numCoeffs
                                      )...
                           );
  }
 
  template <typename... Objs, typename... NumCoeffsTuples, std::size_t... Is>
  auto construct_tuple_from_arrays_impl(const std::tuple<NumCoeffsTuples...>& numCoeffs,
                                        enum init init,
                                        iganet::Options<value_type> options,
                                        std::index_sequence<Is...>) {
    static_assert(sizeof...(Objs) == sizeof...(NumCoeffsTuples));
    return std::make_tuple(
                           std::apply(
                                      [&](auto&&... args) { return Objs(std::forward<decltype(args)>(args)..., init, options); },
                                      std::get<Is>(numCoeffs)
                                      )...
                           );
  }
    
  template <typename... Objs, typename... NumCoeffsTuples>
  auto construct_tuple_from_tuples(const std::tuple<NumCoeffsTuples...>& numCoeffs,
                                   enum init init,
                                   iganet::Options<value_type> options) {
    return construct_tuple_from_tuples_impl<Objs...>(numCoeffs,
                                                     init,
                                                     options,
                                                     std::index_sequence_for<Objs...>{});
  }
 
public:
  explicit IgABase2(
      iganet::Options<value_type> options = iganet::Options<value_type>{})
    : inputs_(), outputs_() {}
 
  template <std::size_t NumCoeffs>
  IgABase2(const std::array<int64_t, NumCoeffs>& ncoeffs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : IgABase2(std::tuple{ncoeffs}, std::tuple{ncoeffs}, init, options) {}
 
  template <std::size_t NumCoeffsInputs, std::size_t NumCoeffsOutputs>
  IgABase2(const std::array<int64_t, NumCoeffsInputs>& ncoeffsInputs,
           const std::array<int64_t, NumCoeffsOutputs>& ncoeffsOutputs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : IgABase2(std::tuple{ncoeffsInputs}, std::tuple{ncoeffsOutputs},
               init, options) {}
 
  template <std::size_t... NumCoeffs>
  IgABase2(const std::tuple<std::array<int64_t, NumCoeffs>...>& ncoeffs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : IgABase2(ncoeffs, ncoeffs, init, options)
  {}
  
  template <std::size_t... NumCoeffsInputs, std::size_t... NumCoeffsOutputs>
  IgABase2(const std::tuple<std::array<int64_t, NumCoeffsInputs>...>& ncoeffsInputs,
           const std::tuple<std::array<int64_t, NumCoeffsOutputs>...>& ncoeffsOutputs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : inputs_(construct_tuple_from_arrays<Inputs...>(ncoeffsInputs, init, options)),
      outputs_(construct_tuple_from_arrays<Outputs...>(ncoeffsOutputs, init, options))
  {}
 
  template<typename... CoeffsInputs, typename... CoeffsOutputs>
  IgABase2(const std::tuple<CoeffsInputs...>& coeffsInputs,
           const std::tuple<CoeffsOutputs...>& coeffsOutputs,
           enum init init = init::greville,
           iganet::Options<value_type> options = iganet::Options<value_type>{})
    : inputs_(construct_tuple_from_tuples<Inputs...>(coeffsInputs, init, options)),
      outputs_(construct_tuple_from_tuples<Outputs...>(coeffsOutputs, init, options))
  {}
 
  inline constexpr const auto &inputs() const {
    return inputs_;
  }
 
  inline constexpr auto &inputs() { 
    return inputs_;
  }
  
  template <std::size_t index>
  inline constexpr const auto &input() const {
    static_assert(index >= 0 && index < sizeof...(Inputs));
    return std::get<index>(inputs_);
  }
 
  template <std::size_t index>
  inline constexpr auto &input() { 
    static_assert(index >= 0 && index < sizeof...(Inputs));
    return std::get<index>(inputs_);
  }
 
  inline constexpr const auto &outputs() const {
    return outputs_;
  }
 
  inline constexpr auto &outputs() { 
    return outputs_;
  }
  
  template <std::size_t index>
  inline constexpr const auto &output() const {
    static_assert(index >= 0 && index < sizeof...(Outputs));
    return std::get<index>(outputs_);
  }
 
  template <std::size_t index>
  inline constexpr auto &output() { 
    static_assert(index >= 0 && index < sizeof...(Outputs));
    return std::get<index>(outputs_);
  }
 
  inline constexpr const auto &collPts() const {
    return outputs_;
  }
 
  inline constexpr auto &collPts() { 
    return outputs_;
  }
  
private:
  template <std::size_t index, std::size_t... Is>
  std::tuple_element_t<index, collPts_type>
  collPts(enum collPts collPtsType, std::index_sequence<Is...>) const {
    std::tuple_element_t<index, collPts_type> collPts;
 
    switch (collPtsType) {
 
    case collPts::greville:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
        std::get<index>(outputs_).template space<Is>().greville(/* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(outputs_).template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_interior:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            std::get<index>(outputs_).template space<Is>().greville(/* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(outputs_).template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_ref1:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            std::get<index>(outputs_).template space<Is>().clone().uniform_refine().greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            std::get<index>(outputs_).template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref1:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            std::get<index>(outputs_).template space<Is>().clone().uniform_refine().greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            std::get<index>(outputs_).template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_ref2:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            std::get<index>(outputs_).template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(outputs_).template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref2:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            std::get<index>(outputs_).template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = std::get<index>(outputs_).template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    default:
      throw std::runtime_error("Invalid collocation point specifier");
    }
 
    return collPts;
  }  
 
public:
  template<std::size_t index>
  std::tuple_element_t<index, collPts_type>
  collPts(enum collPts collPts) const {
    if constexpr (std::tuple_element_t<index, outputs_type>::nspaces() == 1)
 
      switch (collPts) {
 
      case collPts::greville:
        return {std::get<index>(outputs_).space().greville(/* interior */ false),
                std::get<index>(outputs_).boundary().greville()};
 
      case collPts::greville_interior:
        return {std::get<index>(outputs_).space().greville(/* interior */ true),
                std::get<index>(outputs_).boundary().greville()};
 
      case collPts::greville_ref1:
        return {
            std::get<index>(outputs_).space().clone().uniform_refine().greville(/* interior */ false),
            std::get<index>(outputs_).boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_interior_ref1:
        return {
            std::get<index>(outputs_).space().clone().uniform_refine().greville(/* interior */ true),
            std::get<index>(outputs_).boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_ref2:
        return {std::get<index>(outputs_).space().clone().uniform_refine(2, -1).greville(
                    /* interior */ false),
                std::get<index>(outputs_).boundary().clone().uniform_refine(2, -1).greville()};
 
      case collPts::greville_interior_ref2:
        return {std::get<index>(outputs_).space().clone().uniform_refine(2, -1).greville(
                    /* interior */ true),
                std::get<index>(outputs_).boundary().clone().uniform_refine(2, -1).greville()};
 
      default:
        throw std::runtime_error("Invalid collocation point specifier");
      }
 
    else
      return collPts(
                     collPts, std::make_index_sequence<std::tuple_element_t<index, outputs_type>::nspaces()>{});
  }  
};
 
 
 
  
 
 
 
  
  template <typename GeometryMap, typename Variable>
  requires FunctionSpaceType<GeometryMap> && FunctionSpaceType<Variable>
  class IgABaseNoRefData {
public:
  using value_type =
      typename std::common_type<typename GeometryMap::value_type,
                                typename Variable::value_type>::type;
 
  using geometryMap_type = GeometryMap;
 
  using variable_type = Variable;
 
  using geometryMap_collPts_type =
      std::pair<typename GeometryMap::eval_type,
                typename GeometryMap::boundary_eval_type>;
 
  using variable_collPts_type =
      std::pair<typename Variable::eval_type,
                typename Variable::boundary_eval_type>;
 
  bool static constexpr has_GeometryMap = true;
 
  bool static constexpr has_RefData = false;
 
  bool static constexpr has_Solution = true;
 
protected:
  GeometryMap G_;
 
  Variable u_;
 
private:
  template <std::size_t... GeometryMapNumCoeffs, std::size_t... Is,
            std::size_t... VariableNumCoeffs, std::size_t... Js>
  IgABaseNoRefData(
      std::tuple<std::array<int64_t, GeometryMapNumCoeffs>...>
          geometryMapNumCoeffs,
      std::index_sequence<Is...>,
      std::tuple<std::array<int64_t, VariableNumCoeffs>...> variableNumCoeffs,
      std::index_sequence<Js...>,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : // Construct the different spline objects individually
        G_(std::get<Is>(geometryMapNumCoeffs)..., init::greville, options),
        u_(std::get<Js>(variableNumCoeffs)..., init::random, options) {}
 
public:
  explicit IgABaseNoRefData(
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : G_(), u_() {}
 
  template <std::size_t NumCoeffs>
  IgABaseNoRefData(
      std::array<int64_t, NumCoeffs> numCoeffs,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABaseNoRefData(std::tuple{numCoeffs}, std::tuple{numCoeffs}, options) {}
 
  template <std::size_t... NumCoeffs>
  IgABaseNoRefData(
      std::tuple<std::array<int64_t, NumCoeffs>...> numCoeffs,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABaseNoRefData(numCoeffs, numCoeffs, options) {}
 
  template <std::size_t GeometryMapNumCoeffs, std::size_t VariableNumCoeffs>
  IgABaseNoRefData(
      std::array<int64_t, GeometryMapNumCoeffs> geometryMapNumCoeffs,
      std::array<int64_t, VariableNumCoeffs> variableNumCoeffs,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABaseNoRefData(std::tuple{geometryMapNumCoeffs},
                         std::tuple{variableNumCoeffs}, options) {}
 
  template <std::size_t... GeometryMapNumCoeffs,
            std::size_t... VariableNumCoeffs>
  IgABaseNoRefData(
      std::tuple<std::array<int64_t, GeometryMapNumCoeffs>...>
          geometryMapNumCoeffs,
      std::tuple<std::array<int64_t, VariableNumCoeffs>...> variableNumCoeffs,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABaseNoRefData(
            geometryMapNumCoeffs,
            std::make_index_sequence<sizeof...(GeometryMapNumCoeffs)>{},
            variableNumCoeffs,
            std::make_index_sequence<sizeof...(VariableNumCoeffs)>{}, options) {
  }
 
  inline const GeometryMap &G() const { return G_; }
 
  inline GeometryMap &G() { return G_; }
 
  inline const Variable &u() const { return u_; }
 
  inline Variable &u() { return u_; }
 
private:
  template <std::size_t... Is>
  geometryMap_collPts_type
  geometryMap_collPts(enum collPts collPtsType, std::index_sequence<Is...>) const {
    geometryMap_collPts_type collPts;
 
    switch (collPtsType) {
 
    case collPts::greville:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            G_.template space<Is>().greville(/* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = G_.template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_interior:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            G_.template space<Is>().greville(/* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = G_.template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_ref1:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            G_.template space<Is>().clone().uniform_refine().greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            G_.template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref1:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            G_.template space<Is>().clone().uniform_refine().greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            G_.template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_ref2:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            G_.template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = G_.template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref2:
      // Get Greville abscissae inside the domain
      ((std::get<Is>(collPts.first) =
            G_.template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = G_.template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    default:
      throw std::runtime_error("Invalid collocation point specifier");
    }
 
    return collPts;
  }
 
  template <std::size_t... Is>
  variable_collPts_type variable_collPts(enum collPts collPtsType,
                                         std::index_sequence<Is...>) const {
    variable_collPts_type collPts;
 
    switch (collPtsType) {
 
    case collPts::greville:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            u_.template space<Is>().greville(/* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = u_.template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_interior:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            u_.template space<Is>().greville(/* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = u_.template boundary<Is>().greville()),
       ...);
      break;
 
    case collPts::greville_ref1:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            u_.template space<Is>().clone().uniform_refine().greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            u_.template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref1:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            u_.template space<Is>().clone().uniform_refine().greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) =
            u_.template boundary<Is>().clone().uniform_refine().greville()),
       ...);
      break;
 
    case collPts::greville_ref2:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            u_.template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ false)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = u_.template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    case collPts::greville_interior_ref2:
      // Get Greville abscissae inside the domain and at the boundary
      ((std::get<Is>(collPts.first) =
            u_.template space<Is>().clone().uniform_refine(2, -1).greville(
                /* interior */ true)),
       ...);
 
      // Get Greville abscissae at the domain
      ((std::get<Is>(collPts.second) = u_.template boundary<Is>()
                                            .clone()
                                            .uniform_refine(2, -1)
                                            .greville()),
       ...);
      break;
 
    default:
      throw std::runtime_error("Invalid collocation point specifier");
    }
 
    return collPts;
  }
 
public:
  virtual geometryMap_collPts_type
  geometryMap_collPts(enum collPts collPts) const {
    if constexpr (GeometryMap::nspaces() == 1)
 
      switch (collPts) {
 
      case collPts::greville:
        return {G_.space().greville(/* interior */ false),
                G_.boundary().greville()};
 
      case collPts::greville_interior:
        return {G_.space().greville(/* interior */ true),
                G_.boundary().greville()};
 
      case collPts::greville_ref1:
        return {
            G_.space().clone().uniform_refine().greville(/* interior */ false),
            G_.boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_interior_ref1:
        return {
            G_.space().clone().uniform_refine().greville(/* interior */ true),
            G_.boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_ref2:
        return {G_.space().clone().uniform_refine(2, -1).greville(
                    /* interior */ false),
                G_.boundary().clone().uniform_refine(2, -1).greville()};
 
      case collPts::greville_interior_ref2:
        return {G_.space().clone().uniform_refine(2, -1).greville(
                    /* interior */ true),
                G_.boundary().clone().uniform_refine(2, -1).greville()};
 
      default:
        throw std::runtime_error("Invalid collocation point specifier");
      }
 
    else
      return geometryMap_collPts(
          collPts, std::make_index_sequence<GeometryMap::nspaces()>{});
  }
 
  virtual variable_collPts_type variable_collPts(enum collPts collPts) const {
    if constexpr (Variable::nspaces() == 1)
 
      switch (collPts) {
 
      case collPts::greville:
        return {u_.space().greville(/* interior */ false),
                u_.boundary().greville()};
 
      case collPts::greville_interior:
        return {u_.space().greville(/* interior */ true),
                u_.boundary().greville()};
 
      case collPts::greville_ref1:
        return {
            u_.space().clone().uniform_refine().greville(/* interior */ false),
            u_.boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_interior_ref1:
        return {
            u_.space().clone().uniform_refine().greville(/* interior */ true),
            u_.boundary().clone().uniform_refine().greville()};
 
      case collPts::greville_ref2:
        return {u_.space().clone().uniform_refine(2, -1).greville(
                    /* interior */ false),
                u_.boundary().clone().uniform_refine(2, -1).greville()};
 
      case collPts::greville_interior_ref2:
        return {u_.space().clone().uniform_refine(2, -1).greville(
                    /* interior */ true),
                u_.boundary().clone().uniform_refine(2, -1).greville()};
 
      default:
        throw std::runtime_error("Invalid collocation point specifier");
      }
 
    else
      return variable_collPts(collPts,
                              std::make_index_sequence<Variable::nspaces()>{});
  }
};
 
  template <typename GeometryMap, typename Variable>
    requires FunctionSpaceType<GeometryMap> && FunctionSpaceType<Variable>
class IgABase : public IgABaseNoRefData<GeometryMap, Variable> {
public:
  using Base = IgABaseNoRefData<GeometryMap, Variable>;
 
  using value_type = typename Base::value_type;
 
  using geometryMap_type = GeometryMap;
 
  using variable_type = Variable;
 
  using geometryMap_collPts_type = typename Base::geometryMap_collPts_type;
 
  using variable_collPts_type = typename Base::variable_collPts_type;
 
  bool static constexpr has_GeometryMap = true;
 
  bool static constexpr has_RefData = true;
 
  bool static constexpr has_Solution = true;
 
protected:
  Variable f_;
 
private:
  template <std::size_t... GeometryMapNumCoeffs, std::size_t... Is,
            std::size_t... VariableNumCoeffs, std::size_t... Js>
  IgABase(
      std::tuple<std::array<int64_t, GeometryMapNumCoeffs>...>
          geometryMapNumCoeffs,
      std::index_sequence<Is...>,
      std::tuple<std::array<int64_t, VariableNumCoeffs>...> variableNumCoeffs,
      std::index_sequence<Js...>,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : // Construct the different spline objects individually
        Base(geometryMapNumCoeffs, variableNumCoeffs, options),
        f_(std::get<Js>(variableNumCoeffs)..., init::zeros, options) {}
 
public:
  explicit IgABase(
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : Base::G_(), f_(), Base::u_() {}
 
  template <std::size_t NumCoeffs>
  IgABase(std::array<int64_t, NumCoeffs> numCoeffs,
          iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABase(std::tuple{numCoeffs}, std::tuple{numCoeffs}, options) {}
 
  template <std::size_t... NumCoeffs>
  IgABase(std::tuple<std::array<int64_t, NumCoeffs>...> numCoeffs,
          iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABase(numCoeffs, numCoeffs, options) {}
 
  template <std::size_t GeometryMapNumCoeffs, std::size_t VariableNumCoeffs>
  IgABase(std::array<int64_t, GeometryMapNumCoeffs> geometryMapNumCoeffs,
          std::array<int64_t, VariableNumCoeffs> variableNumCoeffs,
          iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABase(std::tuple{geometryMapNumCoeffs}, std::tuple{variableNumCoeffs},
                options) {}
 
  template <std::size_t... GeometryMapNumCoeffs,
            std::size_t... VariableNumCoeffs>
  IgABase(
      std::tuple<std::array<int64_t, GeometryMapNumCoeffs>...>
          geometryMapNumCoeffs,
      std::tuple<std::array<int64_t, VariableNumCoeffs>...> variableNumCoeffs,
      iganet::Options<value_type> options = iganet::Options<value_type>{})
      : IgABase(geometryMapNumCoeffs,
                std::make_index_sequence<sizeof...(GeometryMapNumCoeffs)>{},
                variableNumCoeffs,
                std::make_index_sequence<sizeof...(VariableNumCoeffs)>{},
                options) {}
 
  inline const Variable &f() const { return f_; }
 
  inline Variable &f() { return f_; }
};
 
class IgADatasetBase {
protected:
  template <typename T>
  inline void read_from_xml(std::string location, T &obj,
                            std::vector<torch::Tensor> &v) {
 
    std::filesystem::path path(location);
 
    if (std::filesystem::exists(path)) {
      if (std::filesystem::is_regular_file(path)) {
        try {
          pugi::xml_document doc;
          doc.load_file(path.c_str());
          v.emplace_back(obj.from_xml(doc).as_tensor());
        } catch (...) {
        }
      } else if (std::filesystem::is_directory(path)) {
        for (const auto &file : std::filesystem::directory_iterator(path)) {
          if (file.is_regular_file() && file.path().extension() == ".xml") {
            try {
              pugi::xml_document doc;
              doc.load_file(file.path().c_str());
              v.emplace_back(obj.from_xml(doc).as_tensor());
            } catch (...) {
            }
          }
        }
      } else
        throw std::runtime_error(
            "The path refers to neither a file nor a directory");
    } else
      throw std::runtime_error("The path does not exist");
  }
};
 
template <bool solution = false> class IgADataset;
 
template <>
class IgADataset<false>
    : public IgADatasetBase,
      public torch::data::Dataset<
          IgADataset<false>,
          torch::data::Example<torch::Tensor, torch::data::example::NoTarget>> {
private:
  std::vector<torch::Tensor> G_;
 
  std::vector<torch::Tensor> f_;
 
public:
  using example_type =
      torch::data::Example<torch::Tensor, torch::data::example::NoTarget>;
 
  template <typename T> void add_geometryMap(T &obj, std::string location) {
    read_from_xml(location, obj, G_);
  }
 
  template <typename T> void add_geometryMap(T &&obj, std::string location) {
    read_from_xml(location, obj, G_);
  }
 
  template <typename T>
  void add_geometryMap(T &obj, const pugi::xml_document &doc, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_geometryMap(T &&obj, const pugi::xml_document &doc, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_geometryMap(T &obj, const pugi::xml_node &root, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T>
  void add_geometryMap(T &&obj, const pugi::xml_node &root, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T> void add_referenceData(T &obj, std::string location) {
    read_from_xml(location, obj, f_);
  }
 
  template <typename T> void add_referenceData(T &&obj, std::string location) {
    read_from_xml(location, obj, f_);
  }
 
  template <typename T>
  void add_referenceData(T &obj, const pugi::xml_document &doc, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_referenceData(T &&obj, const pugi::xml_document &doc, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_referenceData(T &obj, const pugi::xml_node &root, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T>
  void add_referenceData(T &&obj, const pugi::xml_node &root, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T, typename Func>
  void add_referenceData(T &obj, Func func) {
    f_.emplace_back(obj.transform(func).as_tensor());
  }
 
  template <typename T, typename Func>
  void add_referenceData(T &&obj, Func func) {
    f_.emplace_back(obj.transform(func).as_tensor());
  }
 
  inline example_type get(std::size_t index) override {
 
    std::size_t geo_index = index / (f_.empty() ? 1 : f_.size());
    std::size_t ref_index = index - geo_index * f_.size();
 
    if (!G_.empty()) {
      if (!f_.empty())
        return torch::cat({G_.at(geo_index), f_.at(ref_index)});
      else
        return G_.at(geo_index);
    } else {
      if (!f_.empty())
        return f_.at(ref_index);
      else
        throw std::runtime_error("No geometry maps and reference data");
    }
  };
 
  // @brief Return the total size of the data set
  inline torch::optional<std::size_t> size() const override {
    return (G_.empty() ? 1 : G_.size()) * (f_.empty() ? 1 : f_.size());
  }
};
 
template <>
class IgADataset<true>
    : public IgADatasetBase,
      public torch::data::Dataset<IgADataset<true>, torch::data::Example<>> {
private:
  std::vector<torch::Tensor> G_;
 
  std::vector<torch::Tensor> f_;
 
  std::vector<torch::Tensor> u_;
 
public:
  template <typename T> void add_geometryMap(T &obj, std::string location) {
    read_from_xml(location, obj, G_);
  }
 
  template <typename T> void add_geometryMap(T &&obj, std::string location) {
    read_from_xml(location, obj, G_);
  }
 
  template <typename T>
  void add_geometryMap(T &obj, const pugi::xml_document &doc, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_geometryMap(T &&obj, const pugi::xml_document &doc, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_geometryMap(T &obj, const pugi::xml_node &root, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T>
  void add_geometryMap(T &&obj, const pugi::xml_node &root, int id = 0,
                       std::string label = "") {
    G_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T> void add_referenceData(T &obj, std::string location) {
    read_from_xml(location, obj, f_);
  }
 
  template <typename T> void add_referenceData(T &&obj, std::string location) {
    read_from_xml(location, obj, f_);
  }
 
  template <typename T>
  void add_referenceData(T &obj, const pugi::xml_document &doc, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_referenceData(T &&obj, const pugi::xml_document &doc, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_referenceData(T &obj, const pugi::xml_node &root, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T>
  void add_referenceData(T &&obj, const pugi::xml_node &root, int id = 0,
                         std::string label = "") {
    f_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T, typename Func>
  void add_referenceData(T &obj, Func func) {
    f_.emplace_back(obj.transform(func).as_tensor());
  }
 
  template <typename T, typename Func>
  void add_referenceData(T &&obj, Func func) {
    f_.emplace_back(obj.transform(func).as_tensor());
  }
 
  template <typename T> void add_solution(T &obj, std::string location) {
    read_from_xml(location, obj, u_);
  }
 
  template <typename T> void add_solution(T &&obj, std::string location) {
    read_from_xml(location, obj, u_);
  }
 
  template <typename T>
  void add_solution(T &obj, const pugi::xml_document &doc, int id = 0,
                    std::string label = "") {
    u_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_solution(T &&obj, const pugi::xml_document &doc, int id = 0,
                    std::string label = "") {
    u_.emplace_back(obj.from_xml(doc.child("xml"), id, label).as_tensor());
  }
 
  template <typename T>
  void add_solution(T &obj, const pugi::xml_node &root, int id = 0,
                    std::string label = "") {
    u_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  template <typename T>
  void add_solution(T &&obj, const pugi::xml_node &root, int id = 0,
                    std::string label = "") {
    u_.emplace_back(obj.from_xml(root, id, label).as_tensor());
  }
 
  inline torch::data::Example<> get(std::size_t index) override {
 
    std::size_t geo_index = index / (f_.empty() ? 1 : f_.size());
    std::size_t ref_index = index - geo_index * f_.size();
 
    if (!G_.empty()) {
      if (!f_.empty())
        return {torch::cat({G_.at(geo_index), f_.at(ref_index)}), u_.at(index)};
      else
        return {G_.at(geo_index), u_.at(index)};
    } else {
      if (!f_.empty())
        return {f_.at(ref_index), u_.at(index)};
      else
        throw std::runtime_error("No geometry maps and reference data");
    }
  };
 
  // @brief Return the total size of the data set
  inline torch::optional<std::size_t> size() const override {
    return (G_.empty() ? 1 : G_.size()) * (f_.empty() ? 1 : f_.size());
  }
};
 
} // namespace iganet