ezsolver.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <net/igabase.hpp>
#include <solver/igasolver.hpp>
#include <splines/bspline.hpp>
#include <utils/matrix.hpp>
 
namespace iganet {
 
template <FunctionSpaceType GeometryMap, FunctionSpaceType Variable>
class EZSolverBase : public iganet::IgASolver<std::tuple<GeometryMap>, std::tuple<Variable>>,
                     public iganet::IgANetCustomizable<std::tuple<GeometryMap>, std::tuple<Variable>> {
  
protected:
  using Base = iganet::IgASolver<std::tuple<GeometryMap>, std::tuple<Variable>>;
  
  Base::template collPts_t<0> collPts_;
 
  using Customizable = iganet::IgANetCustomizable<std::tuple<GeometryMap>, std::tuple<Variable>>;
 
  Customizable::template input_interior_knot_indices_t<0> G_knot_indices_;
 
  Customizable::template input_boundary_knot_indices_t<0> G_knot_indices_boundary_;
 
  Customizable::template output_interior_knot_indices_t<0> var_knot_indices_;
 
  Customizable::template output_boundary_knot_indices_t<0> var_knot_indices_boundary_;
 
public:
  template <std::size_t GeometryMapNumCoeffs, std::size_t VariableNumCoeffs>
  EZSolverBase(const std::array<int64_t, GeometryMapNumCoeffs> &geometryMapNumCoeffs,
               const std::array<int64_t, VariableNumCoeffs> &variableNumCoeffs)
    : Base(std::make_tuple(geometryMapNumCoeffs),
           std::make_tuple(variableNumCoeffs)) {}
  
  auto const &collPts() const { return collPts_; }
  
  
  auto const &G() const { return Base::template input<0>(); }
  
  auto &G() { return Base::template input<0>(); }
  
  auto const &u() const { return Base::template output<0>(); }
  
  auto &u() { return Base::template output<0>(); }
  
  void init() override {
    collPts_ =
      Base::template collPts<0>(iganet::collPts::greville);
    G_knot_indices_ =
      G().template find_knot_indices<iganet::functionspace::interior>(collPts_.first);
    G_knot_indices_boundary_ =
      G().template find_knot_indices<iganet::functionspace::boundary>(collPts_.second);    
    var_knot_indices_ =
      u().template find_knot_indices<iganet::functionspace::interior>(collPts_.first);
    var_knot_indices_boundary_ =
      u().template find_knot_indices<iganet::functionspace::boundary>(collPts_.second);    
  }
};
 
template <FunctionSpaceType GeometryMap, FunctionSpaceType Variable>  
class EZSolver : public EZSolverBase<GeometryMap, Variable> {
private:
  using Base = EZSolverBase<GeometryMap, Variable>;
  
  std::function<
    std::array<torch::Tensor, Variable::template geoDim<0>()>(
                                                              const std::array<torch::Tensor, Variable::template parDim<0>()>&
                                                              )
    > rhs_;
 
public:
  EZSolver(const GeometryMap& geometryMap, const Variable& variable,
           const std::function<
           std::array<torch::Tensor, Variable::template geoDim<0>()>(
                                                                     const std::array<torch::Tensor, Variable::template parDim<0>()>&
                                                                     )
           >& rhs)
    : EZSolverBase<GeometryMap, Variable>(geometryMap.template space<0>().ncoeffs(), variable.template space<0>().ncoeffs()),
      rhs_(rhs) {}
  
  void assembleLhs() override {
 
    auto S_xx = this->u().template eval_basfunc<functionspace::interior,
                                                (deriv::dx^2)>(Base::collPts_.first,
                                                               Base::var_knot_indices_);
    auto S_yy = this->u().template eval_basfunc<functionspace::interior,
                                                (deriv::dy^2)>(Base::collPts_.first,
                                                               Base::var_knot_indices_);
    
    auto M = this->u().template eval_basfunc(Base::collPts_.first,Base::var_knot_indices_);
 
    //M = torch::zeros({16, 100});
        
    auto mask = (Base::collPts_.first[0] == 0.0) | (Base::collPts_.first[0] == 1.0) | (Base::collPts_.first[1] == 0.0) | (Base::collPts_.first[1] == 1.0);
    
    Base::lhs_ = iganet::utils::to_sparseCsrTensor(Base::var_knot_indices_,
                                                   this->u().template space<0>().degrees(),
                                                   this->u().template space<0>().ncoeffs(),
                                                   torch::where(mask, M, S_xx+S_yy).t(),
                                                   { this->u().template space<0>().ncumcoeffs(),
                                                      this->u().template space<0>().ncumcoeffs() });          
  }
  
  void assembleRhs() override {
    Base::rhs_ = rhs_(Base::collPts().first)[0];     
  }
};
  
template <FunctionSpaceType GeometryMap, FunctionSpaceType Variable>  
class EZInterpolation : public EZSolverBase<GeometryMap, Variable> {
private:
  using Base = EZSolverBase<GeometryMap, Variable>;
  
  std::function<
    std::array<torch::Tensor, Variable::template geoDim<0>()>(
                                                              const std::array<torch::Tensor, Variable::template parDim<0>()>&
                                                              )
    > rhs_;
  
public:
  EZInterpolation(const GeometryMap& geometryMap, const Variable& variable,
                  const std::function<
                  std::array<torch::Tensor, Variable::template geoDim<0>()>(
                                                                            const std::array<torch::Tensor, Variable::template parDim<0>()>&
                                                                            )
                  >& rhs)
    : EZSolverBase<GeometryMap, Variable>(geometryMap.template space<0>().ncoeffs(), variable.template space<0>().ncoeffs()), rhs_(rhs) {}
  
  void assembleLhs() override {    
    Base::lhs_ = iganet::utils::to_sparseCsrTensor(Base::var_knot_indices_,
                                                   this->u().template space<0>().degrees(),
                                                   this->u().template space<0>().ncoeffs(),
                                                   this->u().template eval_basfunc(Base::collPts_.first, Base::var_knot_indices_).t(),
                                                   { this->u().template space<0>().ncumcoeffs(),
                                                      this->u().template space<0>().ncumcoeffs() });
  }
  
  void assembleRhs() override {
    Base::rhs_ = rhs_(Base::collPts().first)[0];     
  }
};
  
template <FunctionSpaceType GeometryMap, FunctionSpaceType Variable>
auto ezinterp(const GeometryMap& geometryMap,
              const Variable& variable,
              const std::function<std::array<torch::Tensor, Variable::template geoDim<0>()>(const std::array<torch::Tensor, Variable::template parDim<0>()> &)>
              mapping) {
 
  EZInterpolation interp(geometryMap, variable, mapping);
  interp.init();
  interp.assemble();
  return interp.solve().clone();
}
 
template <FunctionSpaceType GeometryMap, FunctionSpaceType Variable>
auto ezpoisson(const GeometryMap& geometryMap,
               const Variable& variable,
               const std::function<std::array<torch::Tensor, Variable::template geoDim<0>()>(const std::array<torch::Tensor, Variable::template parDim<0>()> &)>
               rhs) {
 
  EZSolver solver(geometryMap, variable, rhs);
  solver.init();
  solver.assemble();
  return solver.solve().clone();
}  
  
} // namespace iganet