ImplicitOperators_fwd.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_OPERATORS_IMPLICITOPERATORS_FWD_HPP_
#define MEDUSA_BITS_OPERATORS_IMPLICITOPERATORS_FWD_HPP_
 
#include <medusa/Config.hpp>
#include <medusa/bits/approximations/Operators_fwd.hpp>
#include <medusa/bits/utils/assert.hpp>
 
namespace mm {
 
template <class shape_storage_type, class matrix_type, class rhs_type>
class ImplicitOperators {
  public:
    typedef shape_storage_type shape_storage_t;  
    typedef matrix_type matrix_t;  
    typedef rhs_type rhs_t;  
    typedef typename shape_storage_t::vector_t vector_t;  
    typedef typename matrix_t::Scalar scalar_t;
    enum {  dim = shape_storage_t::dim };
 
  private:
    const shape_storage_t* ss;
    matrix_t* M;  
    rhs_t* rhs;  
    int row_offset;  
    int col_offset;  
 
    void addToM(int i, int j, scalar_t v);
    void addToRhs(int i, scalar_t v);
 
    template <typename derived_t> class OpBase;  // Forward declaration needed.
 
    class RowOp {
        friend class ImplicitOperators;
      protected:
        ImplicitOperators& op;  
        int row_;  
        RowOp(ImplicitOperators& op, int row) : op(op), row_(row) {}
        int row() const { return op.row_offset + row_; }
      public:
        RowOp operator+(const RowOp& other) const {
            assert_msg(row() == other.row(),
                       "Cannot add together terms for different matrix rows %d and %d.",
                       row(), other.row());
            return *this;
        }
        template <typename derived_t>
        RowOp operator+(const OpBase<derived_t>& right) { return *this + right.eval(); }
 
        void operator=(scalar_t value) {
            op.addToRhs(row_, value);
        }
    };
 
    template <typename derived_t>
    class OpBase {
      protected:
        ImplicitOperators& op;  
        int node;  
        int row;  
        OpBase(ImplicitOperators& op, int node, int row)
                : op(op), node(node), row(row) {}
      public:
        RowOp eval(scalar_t alpha) const {
            return static_cast<const derived_t &>(*this).eval(alpha);
        }
        RowOp eval() const { return eval(1.0); }
        void operator=(scalar_t x) { eval() = x; }
        RowOp operator*(scalar_t alpha) { return eval(alpha); }
        RowOp operator+(RowOp right) { return eval() + right; }
        RowOp operator-() { return eval(-1.0); }
        template <typename other_derived_t>
        RowOp operator+(const OpBase<other_derived_t>& right) { return eval() + right.eval(); }
        friend RowOp operator*(scalar_t alpha, const OpBase& o) { return o.eval(alpha); }
    };
 
#define USE_BASE(Name) \
  private: \
    using OpBase<Name>::op; \
    using OpBase<Name>::node; \
    using OpBase<Name>::row; \
    using OpBase<Name>::OpBase; \
    friend class ImplicitOperators; \
  public: \
    using OpBase<Name>::eval; \
    using OpBase<Name>::operator=;
 
    class ValueOp : public OpBase<ValueOp> {
        USE_BASE(ValueOp)
        RowOp eval(scalar_t alpha) const {
            op.addToM(row, node, alpha);
            return RowOp(op, row);
        }
    };
    template <typename op_family_t>
    class BasicOp : public OpBase<BasicOp<op_family_t>> {
        USE_BASE(BasicOp<op_family_t>)
      private:
        typename op_family_t::operator_t o;  
      public:
        BasicOp(ImplicitOperators& op, int node, int row, typename op_family_t::operator_t o) :
                OpBase<BasicOp<op_family_t>>(op, node, row), o(o) {}
        RowOp eval(scalar_t alpha) const {
            // There is some room for minor optimization here: the index does not need to be
            // computed for 1-element families. A special Op could be added for that,
            // but is probably not worth the code complexity.
            int idx = op_family_t::index(o);
            for (int i = 0; i < op.ss->supportSize(node); ++i) {
                op.addToM(row, op.ss->support(node, i), alpha *
                    op.ss->template get<op_family_t>(idx, node, i));
            }
            return RowOp(op, row);
        }
    };
    class GradOp : public OpBase<GradOp> {
        USE_BASE(GradOp)
        vector_t vec;  
 
        GradOp(ImplicitOperators& op, int node, int row, vector_t vec) :
                OpBase<GradOp>(op, node, row), vec(vec) {}
      public:
        RowOp eval(scalar_t alpha) const {
            for (int i = 0; i < op.ss->supportSize(node); ++i) {
                scalar_t value = 0;
                for (int var = 0; var < op.dim; ++var) {
                    value += vec[var] * op.ss->d1(var, node, i);
                }
                op.addToM(row, op.ss->support(node, i), alpha * value);
            }
            return RowOp(op, row);
        }
    };
 
  public:
    ImplicitOperators() : ss(nullptr), M(nullptr), rhs(nullptr), row_offset(0), col_offset(0) {}
    explicit ImplicitOperators(const shape_storage_t& ss) : ss(&ss), M(nullptr), rhs(nullptr),
                                                            row_offset(0), col_offset(0) {}
    ImplicitOperators(const shape_storage_t& ss, matrix_t& M, rhs_t& rhs,
                      int row_offset = 0, int col_offset = 0);
 
    void setProblem(matrix_t& M, rhs_t& rhs, int row_offset = 0, int col_offset = 0) {
        ImplicitOperators::M = &M;
        ImplicitOperators::rhs = &rhs;
        setRowOffset(row_offset);
        setColOffset(col_offset);
    }
    void setRowOffset(int row_offset) {
        assert_msg(0 <= row_offset, "Row offset cannot be negative, got %d.", row_offset);
        ImplicitOperators::row_offset = row_offset;
    }
    void setColOffset(int col_offset) {
        assert_msg(0 <= col_offset, "Col offset cannot be negative, got %d.", col_offset);
        ImplicitOperators::col_offset = col_offset;
    }
 
    bool hasShapes() const { return ss != nullptr; }
    bool hasMatrix() const { return M != nullptr; }
    bool hasRhs() const { return rhs != nullptr; }
 
    ValueOp value(int node, int row) { return ValueOp(*this, node, row); }
    ValueOp value(int node) { return value(node, node); }
 
    BasicOp<Lap<dim>> lap(int node, int row) { return BasicOp<Lap<dim>>(*this, node, row, {}); }
    BasicOp<Lap<dim>> lap(int node) { return lap(node, node); }
 
    GradOp grad(int node, const vector_t& v, int row) { return GradOp(*this, node, row, v); }
    GradOp grad(int node, const vector_t& v) { return grad(node, v, node); }
 
    GradOp neumann(int node, const vector_t& unit_normal, int row) {
        assert_msg(std::abs(unit_normal.norm() - 1.0) < 1e-13,
                   "Normal %s is not of unit length, got %.16f.", unit_normal, unit_normal.norm());
        return grad(node, unit_normal, row);
    }
    GradOp neumann(int node, const vector_t& n) { return neumann(node, n, node); }
 
    BasicOp<Der1s<dim>> der1(int node, int var, int row) {
        return BasicOp<Der1s<dim>>(*this, node, row, {var});
    }
    BasicOp<Der1s<dim>> der1(int node, int var) { return der1(node, var, node); }
 
    BasicOp<Der2s<dim>> der2(int node, int varmin, int varmax, int row) {
        return BasicOp<Der2s<dim>>(*this, node, row, {varmin, varmax});
    }
    BasicOp<Der2s<dim>> der2(int node, int varmin, int varmax) {
        return der2(node, varmin, varmax, node);
    }
 
 
    template <typename op_family_t>
    BasicOp<op_family_t> apply(int node, typename op_family_t::operator_t o, int row) {
        return BasicOp<op_family_t>(*this, node, row, o);
    }
 
    template <typename op_family_t>
    BasicOp<op_family_t> apply(int node, typename op_family_t::operator_t o) {
        return apply<op_family_t>(node, o, node);
    }
 
    template <typename op_family_t>
    BasicOp<op_family_t> apply(int node) { return apply<op_family_t>(node, {}, node); }
 
 
    template <typename S, typename M, typename R>
    friend std::ostream& operator<<(std::ostream& os, const ImplicitOperators<S, M, R>& op);
};
 
}  // namespace mm
 
#endif  // MEDUSA_BITS_OPERATORS_IMPLICITOPERATORS_FWD_HPP_