ImplicitVectorOperators_fwd.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_OPERATORS_IMPLICITVECTOROPERATORS_FWD_HPP_
#define MEDUSA_BITS_OPERATORS_IMPLICITVECTOROPERATORS_FWD_HPP_
 
#include <medusa/Config.hpp>
#include <medusa/bits/utils/assert.hpp>
#include "ImplicitOperators_fwd.hpp"
 
namespace mm {
 
template <class shape_storage_type, class matrix_type, class rhs_type>
class ImplicitVectorOperators {
  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 domain_vector_t;
    typedef typename matrix_t::Scalar scalar_t;
    enum {  dim = shape_storage_t::dim };
    typedef Eigen::Matrix<scalar_t, dim, 1> vector_t;  
 
  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) {
        assert_msg(v == v, "You tried to set M(%d, %d) to NAN. Check your computes shapes.",
                   row_offset+i, col_offset+j);
        M->coeffRef(row_offset+i, col_offset+j) += v;
    }
    void addToRhs(int i, const vector_t& v) {
        assert_msg(!v.array().isNaN().any(),
                   "You tried to set rhs(%d) to NaN. Check you computed shapes.", row_offset+i);
        for (int d = 0; d < dim; ++d) {
            rhs->operator[](row_offset + d * ss->size() + i) += v[d];
        }
    }
 
    template <typename derived_t> class VecOpBase;  // Forward declaration needed.
 
    class RowVecOp {
        friend class ImplicitVectorOperators;
      protected:
        ImplicitVectorOperators& op;  
        int row_;  
        RowVecOp(ImplicitVectorOperators& op, int row) : op(op), row_(row) {}
        int row() const { return op.row_offset + row_; }
      public:
        RowVecOp operator+(const RowVecOp& 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>
        RowVecOp operator+(const VecOpBase<derived_t>& right) { return *this + right.eval(); }
        void operator=(const vector_t& value) {
            op.addToRhs(row_, value);
        }
    };
 
    template <typename derived_t>
    class VecOpBase {
      protected:
        ImplicitVectorOperators& op;  
        int node;  
        int row;  
        VecOpBase(ImplicitVectorOperators& op, int node, int row)
                : op(op), node(node), row(row) {}
      public:
        RowVecOp eval(scalar_t alpha) const {
            return static_cast<const derived_t &>(*this).eval(alpha);
        }
        RowVecOp eval() const { return eval(1.0); }
        void operator=(const vector_t& x) { eval() = x; }
        RowVecOp operator*(scalar_t alpha) { return eval(alpha); }
        RowVecOp operator+(RowVecOp right) { return eval() + right; }
        RowVecOp operator-() { return eval(-1.0); }
        template <typename other_derived_t>
        RowVecOp operator+(const VecOpBase<other_derived_t>& right) { return eval()+right.eval(); }
        friend RowVecOp operator*(scalar_t alpha, const VecOpBase& op) { return op.eval(alpha); }
    };
 
#define USE_VECTOR_BASE(Name) \
  private: \
    using VecOpBase<Name>::op; \
    using VecOpBase<Name>::node; \
    using VecOpBase<Name>::row; \
    using VecOpBase<Name>::VecOpBase; \
    friend class ImplicitVectorOperators; \
  public: \
    using VecOpBase<Name>::eval; \
    using VecOpBase<Name>::operator=;
 
    class ValueOp : public VecOpBase<ValueOp> {
        USE_VECTOR_BASE(ValueOp)
        RowVecOp eval(scalar_t alpha) const {
            for (int d = 0; d < dim; ++d) {
                op.addToM(row + d*op.ss->size(), node + d*op.ss->size(), alpha);
            }
            return RowVecOp(op, row);
        }
    };
 
    template <typename op_family_t>
    class BasicOp : public VecOpBase<BasicOp<op_family_t>> {
      USE_VECTOR_BASE(BasicOp<op_family_t>)
      private:
        typename op_family_t::operator_t o;  
      public:
        BasicOp(ImplicitVectorOperators& op, int node, int row, typename op_family_t::operator_t o)
                : VecOpBase<BasicOp<op_family_t>>(op, node, row), o(o) {}
        RowVecOp eval(scalar_t alpha) const {
            int idx = op_family_t::index(o);
            for (int d = 0; d < dim; ++d) {
                for (int i = 0; i < op.ss->supportSize(node); ++i) {
                    op.addToM(d*op.ss->size() + row, d*op.ss->size() + op.ss->support(node, i),
                      alpha * op.ss->template get<op_family_t>(idx, node, i));
                }
            }
            return RowVecOp(op, row);
        }
    };
 
    class GradOp : public VecOpBase<GradOp> {
        USE_VECTOR_BASE(GradOp)
      private:
        domain_vector_t vec;  
        GradOp(ImplicitVectorOperators& op, int node, int row, const domain_vector_t& vec) :
                VecOpBase<GradOp>(op, node, row), vec(vec) {}
      public:
        RowVecOp eval(scalar_t alpha) const {
            for (int i = 0; i < op.ss->supportSize(node); ++i) {
                scalar_t shape_value = 0;
                for (int var = 0; var < dim; ++var) {
                    shape_value += vec[var] * op.ss->d1(var, node, i);
                }
                for (int d = 0; d < dim; ++d) {
                    op.addToM(d*op.ss->size() + row, d*op.ss->size() + op.ss->support(node, i),
                              alpha*shape_value);
                }
            }
            return RowVecOp(op, row);
        }
    };
 
    class GradDivOp : public VecOpBase<GradDivOp> {
        USE_VECTOR_BASE(GradDivOp)
        RowVecOp eval(scalar_t alpha) const {
            // Implicit equation is for all d: \sum_k \sum_j \chi_djk u_j(s_k) = v_d(s)
            for (int d = 0; d < dim; ++d) {
                for (int k = 0; k < op.ss->supportSize(node); ++k) {
                    for (int j = 0; j < dim; ++j) {  // loop over dimensions
                        int dmin = std::min(d, j);
                        int dmax = std::max(d, j);
                        op.addToM(d*op.ss->size() + row,
                                  j*op.ss->size() + op.ss->support(node, k),
                                  alpha * op.ss->d2(dmin, dmax, node, k));
                    }
                }
            }
            return RowVecOp(op, row);
        }
    };
 
    class TractionOp : public VecOpBase<TractionOp> {
      USE_VECTOR_BASE(TractionOp)
      private:
        scalar_t lam,  
                 mu;  
        domain_vector_t normal;  
 
        TractionOp(ImplicitVectorOperators& op, int node, int row, scalar_t lam, scalar_t mu,
                   const domain_vector_t& normal) :
                VecOpBase<TractionOp>(op, node, row), lam(lam), mu(mu), normal(normal) {}
      public:
        RowVecOp eval(scalar_t alpha) const {
            for (int i = 0; i < dim; ++i) {  // which equation of sigma.n = t
                for (int j = 0; j < dim; ++j) {  // computation of sigma (lam tr(eps) I) part
                    op.eq(i).c(j).der1(node, j, row).eval(alpha*lam*normal(i));
                }
                for (int j = 0; j < dim; ++j) {  // computation of sigma (2 mu eps) part
                    op.eq(i).c(i).der1(node, j, row).eval(alpha*mu*normal(j));
                    op.eq(i).c(j).der1(node, i, row).eval(alpha*mu*normal(j));
                }
            }
            return RowVecOp(op, row);
        }
    };
 
    class Equation {
        ImplicitVectorOperators& op;  
        int eq_ixd;  
        friend class ImplicitVectorOperators;
 
        Equation(ImplicitVectorOperators& op, int idx) : op(op), eq_ixd(idx) {
            assert_msg(0 <= idx && idx < op.dim, "Equation number %d out of range [0, %d).",
                       idx, dim);
        }
 
      public:
        ImplicitOperators<shape_storage_t, matrix_t, rhs_t> c(int comp) {
            assert_msg(0 <= comp && comp < dim,
                       "Dimension %d out of range, expected in range [0, %d).", comp, dim);
            auto sop = ImplicitOperators<shape_storage_t, matrix_t, rhs_t>(*op.ss);
            sop.setProblem(*op.M, *op.rhs, op.row_offset + eq_ixd*op.ss->size(),
                    op.col_offset + comp*op.ss->size());
            return sop;
        }
    };
 
  public:
    ImplicitVectorOperators() : ss(nullptr), M(nullptr), rhs(nullptr), row_offset(0),
                                col_offset(0) {}
    explicit ImplicitVectorOperators(const shape_storage_t& ss) : ss(&ss), M(nullptr), rhs(nullptr),
                                                                  row_offset(0), col_offset(0) {}
    ImplicitVectorOperators(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) {
        ImplicitVectorOperators::M = &M;
        ImplicitVectorOperators::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);
        ImplicitVectorOperators::row_offset = row_offset;
    }
    void setColOffset(int col_offset) {
        assert_msg(0 <= col_offset, "Col offset cannot be negative, got %d.", col_offset);
        ImplicitVectorOperators::col_offset = col_offset;
    }
 
    bool hasShapes() const { return ss != nullptr; }
    bool hasMatrix() const { return M != nullptr; }
    bool hasRhs() const { return rhs != nullptr; }
    Equation eq(int num) {
        assert_msg(0 <= num && num < dim, "Equation must be in range [0, %d), got %d", num, dim);
        return Equation(*this, num);
    }
 
    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); }
 
    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); }
 
    GradOp grad(int node, const domain_vector_t& v, int row) { return GradOp(*this, node, row, v); }
    GradOp grad(int node, const domain_vector_t& v) { return grad(node, v, node); }
 
    GradDivOp graddiv(int node, int row) { return GradDivOp(*this, node, row); }
    GradDivOp graddiv(int node) { return graddiv(node, node); }
 
    GradOp neumann(int node, const domain_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 domain_vector_t& n) { return neumann(node, n, node); }
 
    TractionOp traction(int node, scalar_t lam, scalar_t mu, const domain_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 TractionOp(*this, node, row, lam, mu, unit_normal);
    }
    TractionOp traction(int node, scalar_t lam, scalar_t mu, const domain_vector_t& n) {
        return traction(node, lam, mu, n, node);
    }
 
    template <typename S, typename M, typename R>
    friend std::ostream& operator<<(std::ostream& os, const ImplicitVectorOperators<S, M, R>& op);
};
 
}  // namespace mm
 
#endif  // MEDUSA_BITS_OPERATORS_IMPLICITVECTOROPERATORS_FWD_HPP_