Compact C ++ library for programming finite-difference methods in operator style. Part 1. Semantics

 //() 2016-2019  , SharipovTR@gmail.com // As Is #pragma once #include <memory> #include <list> #include "NumericAssert.h" #include "Arrays.h" //#define USE_STD_MOVE template<class T, class DerivCell, class ScalarCell, class VectorCell> class MeshCell { public: int dim_; MeshCell() : dim_(0), x_(0), y_(0), left_(0), right_(0), up_(0), down_(0), op_(OpEqual) { } MeshCell(const DerivCell &rhs) : dim_(rhs.dim_), x_(rhs.x_), y_(rhs.y_), left_(rhs.left_), right_(rhs.right_), up_(rhs.up_), down_(rhs.down_), op_(rhs.op_) {}//todo:op? #if defined(USE_STD_MOVE) && (_MSC_VER > 1500) MeshCell(DerivCell &&rhs) : dim_(std::move(rhs.dim_)), x_(std::move(rhs.x_)), y_(std::move(rhs.y_)), left_(std::move(rhs.left_)), right_(std::move(rhs.right_)), up_(std::move(rhs.up_)), down_(std::move(rhs.down_)), op_(std::move(rhs.op_)) {}//todo:op? #endif virtual ~MeshCell() {} virtual ScalarCell &comp(int i) = 0; virtual const ScalarCell &comp(int i) const = 0; virtual DerivCell *This() = 0; void set_x(double x) { x_ = x; if (dim_ > 1) for(int i = 0; i < dim_; i++) { comp(i).set_x(x); } } double x() const { return x_; } void set_y(double y) { y_ = y; if (dim_ > 1) for(int i = 0; i < dim_; i++) { comp(i).set_y(y); } } double y() const { return y_; } void set_left(DerivCell *left) { left_ = left; if (!left) return; left->right_ = This(); if (dim_ > 1) for(int i = 0; i < dim_; i++) { comp(i).set_left(&left->comp(i)); } } DerivCell *left() { return left_; } void set_right(DerivCell *right) { right_ = right; if (!right) return; right->left_ = This(); if (dim_ > 1) for(int i = 0; i < dim_; i++) { comp(i).set_right(&right->comp(i)); } } DerivCell *right() { return right_; } void set_up(DerivCell *up) { up_ = up; if (!up) return; up->down_ = This(); if (dim_ > 1) for(int i = 0; i < dim_; i++) { comp(i).set_up(&up->comp(i)); } } DerivCell *up() { return up_; } void set_down(DerivCell *down) { down_ = down; if (!down) return; down->up_ = This(); if (dim_ > 1) for(int i = 0; i < dim_; i++) { comp(i).set_down(&down->comp(i)); } } DerivCell *down() { return down_; } DerivCell operator - () { return DerivCell(*This(), -1); } virtual DerivCell &Instance(int i) = 0; virtual ScalarCell &InstanceAsScalar(int i) = 0; virtual VectorCell &InstanceAsVector(int i) = 0; DerivCell &dx_left(bool recur = true); DerivCell &dx_right(bool recur = true); DerivCell &dy_right(bool recur = true); DerivCell &dy_left(bool recur = true); DerivCell &laplacian(bool recur = true); VectorCell &grad_left(bool recur = true); VectorCell &grad_right(bool recur = true); ScalarCell &div_left(bool recur = true); ScalarCell &div_right(bool recur = true); DerivCell &operator = (const DerivCell &rhs) { if (!me().left_ && !me().right_ && !me().up_ && !me().down_) { me().dim_ = rhs.dim_; me().x_ = rhs.x_; me().y_ = rhs.y_; me().left_ = rhs.left_; me().right_ = rhs.right_; me().up_ = rhs.up_; me().down_ = rhs.down_; me().op_ = rhs.op_; } else { for (int i = 0; i < dim_; i++) { me().comp(i).val() = rhs.comp(i).val(); } } return *This(); } #if defined(USE_STD_MOVE) && (_MSC_VER > 1500) DerivCell &operator = (DerivCell &&rhs) { if (me().left_ || me().right_ || me().up_ || me().down_) return *This(); dim_ = std::move(rhs.dim_); x_ = std::move(rhs.x_); y_ = std::move(rhs.y_); left_ = std::move(rhs.left_); right_ = std::move(rhs.right_); up_ = std::move(rhs.up_); down_ = std::move(rhs.down_); op_ = std::move(rhs.op_); return *This(); } #endif DerivCell &operator = (const T &val) { for(int i = 0; i < dim_; i++) { me().comp(i).val() = val; } return *This(); } #define DECLARE_CELL_OP(OPERAND)\ DerivCell &operator OPERAND (const DerivCell &rhs) {\ for(int i = 0; i < dim_; i++) {\ me().comp(i).val() OPERAND rhs.comp(i).val();\ }\ return *This();\ }\ DerivCell &operator OPERAND (const T &val) {\ for(int i = 0; i < dim_; i++) {\ me().comp(i).val() OPERAND val;\ }\ return *This();\ } DECLARE_CELL_OP( += ) DECLARE_CELL_OP( -= ) DECLARE_CELL_OP( *= ) DECLARE_CELL_OP( /= ) #undef DECLARE_CELL_OP DerivCell & PlusEq() { return SetOp(OpPlus, 0); } DerivCell &MinusEq() { return SetOp(OpMinus, 0); } DerivCell &MultEq() { return SetOp(OpMult, 1); } DerivCell &DivideEq() { return SetOp(OpDivide, 1); } DerivCell &Eval() { if (OpEqual == op_) { return *This(); } switch(op_) { case OpPlus: op_ = OpEqual; *This() += ft(); break; case OpMinus: op_ = OpEqual; *This() -= ft(); break; case OpMult: op_ = OpEqual; *This() *= ft(); break; case OpDivide: op_ = OpEqual; *This() /= ft(); break; default: DASSERT(0 && "Unknown arithmetic operation"); } return *This(); } protected: enum Operation { OpEqual, OpPlus, OpMinus, OpMult, OpDivide } op_; double x_; double y_; DerivCell *left_; DerivCell *right_; DerivCell *up_; DerivCell *down_; std::auto_ptr<ScalarCell> scal_[7];//todo:  shared_ptr std::auto_ptr<VectorCell> vec_[7];//todo:  shared_ptr DerivCell &me() { return (OpEqual == op_) ? *This() : ft(); } DerivCell &ft() { static DerivCell res(*This());//todo return res; //return Instance(6); } DerivCell &SetOp(Operation op, int ft_val) { Eval(); op_ = op; ft() = ft_val; return *This(); } }; template<class T, class DerivCell, class ScalarCell, class VectorCell> DerivCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::dx_left(bool recur) { DerivCell &res = Instance(0); for(int i = 0; i < dim_; i++) { DASSERT(left_ != 0); DASSERT(left_->x_ != x_); res.comp(i).val() = (comp(i).val() - left_->comp(i).val()) / (x_ - left_->x_); } if (recur) { if (left_ && left_->left_) res.set_left(&left_->dx_left(false)); if (right_ && right_->left_) res.set_right(&right_->dx_left(false)); if (up_ && up_->left_) res.set_up(&up_->dx_left(false)); if (down_ && down_->left_) res.set_down(&down_->dx_left(false)); } return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> DerivCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::dx_right(bool recur) { DerivCell &res = Instance(1); for(int i = 0; i < dim_; i++) { DASSERT(right_ != 0); DASSERT(right_->x_ != x_); res.comp(i).val() = (right_->comp(i).val() - comp(i).val()) / (right_->x_ - x_); } if (recur) { if (left_ && left_->right_) res.set_left(&left_->dx_right(false)); if (right_ && right_->right_) res.set_right(&right_->dx_right(false)); if (up_ && up_->right_) res.set_up(&up_->dx_right(false)); if (down_ && down_->right_) res.set_down(&down_->dx_right(false)); } return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> DerivCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::dy_left(bool recur) { DerivCell &res = Instance(2); for(int i = 0; i < dim_; i++) { DASSERT(down_ != 0); DASSERT(down_->y_ != y_); res.comp(i).val() = (comp(i).val() - down_->comp(i).val()) / (y_ - down_->y_); } if (recur) { if (left_ && left_->down_) res.set_left(&left_->dy_left(false)); if (right_ && right_->down_) res.set_right(&right_->dy_left(false)); if (up_ && up_->down_) res.set_up(&up_->dy_left(false)); if (down_ && down_->down_) res.set_down(&down_->dy_left(false)); } return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> DerivCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::dy_right(bool recur) { DerivCell &res = Instance(3); for(int i = 0; i < dim_; i++) { DASSERT(up_ != 0); DASSERT(up_->y_ != y_); res.comp(i).val() = (up_->comp(i).val() - comp(i).val()) / (up_->y_ - y_); } if (recur) { if (left_ && left_->up_) res.set_left(&left_->dy_right(false)); if (right_ && right_->up_) res.set_right(&right_->dy_right(false)); if (up_ && up_->up_) res.set_up(&up_->dy_right(false)); if (down_ && down_->up_) res.set_down(&down_->dy_right(false)); } return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> DerivCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::laplacian(bool recur) { DerivCell &res = Instance(4); res = dx_left(recur).dx_right(false) + dy_left(recur).dy_right(false);//todo return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> ScalarCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::div_left(bool recur) { ScalarCell &res = InstanceAsScalar(5); if (dim_ >= 1) res = comp(0).dx_left(recur); if (dim_ >= 2) res += comp(1).dy_left(recur); if (recur) { if (left_ && left_->left_ && left_->down_) res.set_left(&left_->div_left(false)); if (right_ && right_->left_ && right_->down_) res.set_right(&right_->div_left(false)); if (up_ && up_->left_ && up_->down_) res.set_up(&up_->div_left(false)); if (down_ && down_->left_ && down_->down_) res.set_down(&down_->div_left(false)); } return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> ScalarCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::div_right(bool recur) { ScalarCell &res = InstanceAsScalar(6); if (dim_ >= 1) res = comp(0).dx_right(recur); if (dim_ >= 2) res += comp(1).dy_right(recur); return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> VectorCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::grad_left(bool recur) { VectorCell &res = InstanceAsVector(0); res.comp(0) = dx_left(recur); res.comp(1) = dy_left(recur); return res; } template<class T, class DerivCell, class ScalarCell, class VectorCell> VectorCell & MeshCell<T, DerivCell, ScalarCell, VectorCell>::grad_right(bool recur) { VectorCell &res = InstanceAsVector(1); res.comp(0) = dx_right(recur); res.comp(1) = dy_right(recur); return res; } template<class T> class VectorCell; template<class T> class ScalarCell : public MeshCell<T, ScalarCell<T>, ScalarCell<T>, VectorCell<T> > { public: explicit ScalarCell(T v = 0) { dim_ = 1; val() = v; } ScalarCell(const ScalarCell &rhs, double mult = 1) : MeshCell(rhs), val_(mult * rhs.val()) { dim_ = 1; } #if defined(USE_STD_MOVE) && (_MSC_VER > 1500) ScalarCell(ScalarCell &&rhs) : MeshCell(std::move(rhs)), val_(std::move(rhs.val_)) {} #endif virtual ~ScalarCell() {} ScalarCell *This() { return this; } ScalarCell<T> &comp(int i) { DASSERT_LE(i, 1); return *this; } const ScalarCell<T> &comp(int i) const { DASSERT_LE(i, 1); return *this; } const T &val() const { return val_; } T &val() { return val_; } const T &val(int i) const { return comp(i).val(); } T &val(int i) { return comp(i).val(); } operator const T &() const { return val(); } operator T &() { return val(); } ScalarCell &Instance(int i) { std::auto_ptr<ScalarCell> &res = scal_[i]; if (!res.get()) res.reset(new ScalarCell(*this)); return *res; } ScalarCell &InstanceAsScalar(int i) { return Instance(i); } VectorCell<T> &InstanceAsVector(int i) { std::auto_ptr<VectorCell<T> > &res = vec_[i]; if (!res.get()) res.reset(new VectorCell<T>); res->set_x(x_); res->set_y(y_); return *res; } ScalarCell &operator = (const T &val) { return MeshCell::operator=(val); } ScalarCell &operator = (const ScalarCell &rhs) { me().val_ = rhs.val_; return MeshCell::operator=(rhs); } #if defined(USE_STD_MOVE) && (_MSC_VER > 1500) ScalarCell &operator = (ScalarCell &&rhs) { val_ = std::move(rhs.val_); return MeshCell::operator=(std::move(rhs)); } #endif public: T val_; }; //================================================================================== template<class T> class VectorCell : public MeshCell<T, VectorCell<T>, ScalarCell<T>, VectorCell<T> > { public: std::vector<ScalarCell<T> *> comp_;//todo->auto_ptr VectorCell() { Resize(2); for(int i = 0; i < dim_; i++) { comp_[i] = new ScalarCell<T>; } } explicit VectorCell(const T &val) { Resize(2); for(int i = 0; i < dim_; i++) { comp_[i] = new ScalarCell<T>(val); } } VectorCell(const VectorCell &rhs, double mult = 1) : MeshCell(rhs) { DASSERT_EQ(dim_, 2); Resize(dim_); for(int i = 0; i < dim_; i++) { comp_[i] = new ScalarCell<T>(rhs.comp(i), mult); } } #if defined(USE_STD_MOVE) && (_MSC_VER > 1500) VectorCell(VectorCell &&rhs) : MeshCell(std::move(rhs)), comp_(std::move(rhs.comp_)) {} #endif virtual ~VectorCell() { Resize(0); } void Resize(int dim) { for(int i = 0, n = comp_.size(); i < n; i++) { delete comp_[i]; } dim_ = dim; if (dim > 0) { comp_.resize(dim_); } } VectorCell *This() { return this; } ScalarCell<T> &comp(int i) { return *comp_[i]; } const ScalarCell<T> &comp(int i) const { return *comp_[i]; } VectorCell &Instance(int i) { std::auto_ptr<VectorCell> &res = vec_[i]; if (!res.get()) res.reset(new VectorCell(*this)); return *res; } ScalarCell<T> &InstanceAsScalar(int i) { std::auto_ptr<ScalarCell<T> > &res = scal_[i]; if (!res.get()) res.reset(new ScalarCell<T>); res->set_x(x_); res->set_y(y_); return *res; } VectorCell &InstanceAsVector(int i) { return Instance(i); } ScalarCell<T> &operator()(int at) { return comp(at); } const ScalarCell<T> &operator()(int at) const { return comp(at); } VectorCell &operator = (const T &val) { return MeshCell::operator=(val); } VectorCell &operator = (const VectorCell &rhs) { for (int i = 0; i < dim_; i++) { me().comp(i) = rhs.comp(i); } return MeshCell::operator=(rhs); } #if defined(USE_STD_MOVE) && (_MSC_VER > 1500) VectorCell &operator = (VectorCell &&rhs) { me().comp_ = std::move(rhs.comp_); return MeshCell::operator=(std::move(rhs)); } #endif const T &max() const { T ret = comp(0).val(); int i_max = 0; for(int i = 1; i < dim_; i++) { if (ret < comp(i).val()) { ret = comp(i).val(); i_max = i; } } return comp(i_max).val(); } const T &min() const { T ret = comp(0).val(); int i_min = 0; for(int i = 1; i < dim_; i++) { if (comp(i).val() < ret) { ret = comp(i).val(); i_min = i; } } return comp(i_min).val(); } }; #define DECLARE_CELL_OP(OPERAND)\ template<class T>\ ScalarCell<T> operator OPERAND (const ScalarCell<T> &c1, const ScalarCell<T> &c2)\ {\ ScalarCell<T> res(c1);\ res OPERAND= c2;\ return res;\ }\ template<class T>\ VectorCell<T> operator OPERAND (const VectorCell<T> &c1, const VectorCell<T> &c2)\ {\ VectorCell<T> res(c1);\ res OPERAND= c2;\ return res;\ }\ template<class T>\ VectorCell<T> operator OPERAND (const VectorCell<T> &c, double val)\ {\ VectorCell<T> res(c);\ res OPERAND= val;\ return res;\ }\ template<class T>\ VectorCell<T> operator OPERAND (double val, const VectorCell<T> &c)\ {\ VectorCell<T> res(c);\ res OPERAND= val;\ return res;\ } DECLARE_CELL_OP(+) DECLARE_CELL_OP(-) DECLARE_CELL_OP(*) DECLARE_CELL_OP( / ) #undef DECLARE_CELL_OP //---------------------------------------------------------------------------------- template<class T> bool operator < (const VectorCell<T> &lhs, const VectorCell<T> &rhs) { for (int i = 0; i < lhs.dim_; i++) if (lhs(i) >= rhs(i)) return 0; return 1; } //---------------------------------------------------------------------------------- template<class T> bool operator <= (const VectorCell<T> &lhs, const VectorCell<T> &rhs) { for (int i = 0; i < lhs.dim_; i++) if (lhs(i) > rhs(i)) return 0; return 1; } //================================================================================== struct Interval2d { public: int lb_x, ub_x; int lb_y, ub_y; Interval2d() : lb_x(0), ub_x(0), lb_y(0), ub_y(0) {} Interval2d(int _lb_x, int _ub_x, int _lb_y, int _ub_y) { Init(_lb_x, _ub_x, _lb_y, _ub_y); } void Init(int _lb_x, int _ub_x, int _lb_y, int _ub_y) { lb_x = _lb_x; ub_x = _ub_x; lb_y = _lb_y; ub_y = _ub_y; } int nx() const { return ub_x - lb_x + 1; } int ny() const { return ub_y - lb_y + 1; } }; //================================================================================== #define FOREACH_INTERVAL(I,i,j)\ for(int i = I.lb_x, j; i <= I.ub_x; i++)\ for(j = I.lb_y; j <= I.ub_y; j++) template<class MeshCell> class MeshCellStack { public: static MeshCellStack &Instance() { static MeshCellStack instance; return instance; } std::list<MeshCell> cells; }; template<class MeshCell> class MeshOperator { public: Interval2d I; MeshOperator(double val = 0) : val_(val) {} virtual bool IsStencil() const { return false; } virtual int CountOfDependentVars() const { return 1; } virtual MeshCell *GetRef(int i, int j) const { return 0; } virtual MeshCell &Eval(int i, int j) const { MeshCell *r = GetRef(i, j); if (r) { return *r; } static MeshCell res;//todo:    OpenMP Eval(&res, i, j); return res; } virtual void Eval(MeshCell *res, int i, int j) const { *res = val_; } MeshCell &max() const { int i_max = I.lb_x; int j_max = I.lb_y; //todo:  val:    MeshCell ret = Eval(i_max, j_max); FOREACH_INTERVAL(I, i, j) { MeshCell &it = Eval(i, j); if (ret < it) { i_max = i; j_max = j; ret = it; } } return Eval(i_max, j_max); } MeshCell &min() const { int i_min = I.lb_x; int j_min = I.lb_y; //todo:  val:    MeshCell ret = Eval(i_min, j_min); FOREACH_INTERVAL(I, i, j) { MeshCell &it = Eval(i, j); if (it < ret) { i_min = i; j_min = j; ret = it; } } return Eval(i_min, j_min); } private: double val_; }; template<class MeshCellIn, class MeshCellOut> class MeshOperatorBind : public MeshOperator<MeshCellOut> { public: typedef MeshCellOut &(*GlobalCellFunc_ref)(MeshCellIn &, bool recur); typedef MeshCellOut (*GlobalCellFunc_copy)(MeshCellIn &, bool recur); MeshOperatorBind(const MeshOperator<MeshCellIn> &rhs, GlobalCellFunc_copy func) : rhs_(rhs), func_copy_(func), func_ref_(0), mult_(1) { I = rhs.I; } MeshOperatorBind(const MeshOperator<MeshCellIn> &rhs, GlobalCellFunc_ref func) : rhs_(rhs), func_copy_(0), func_ref_(func), mult_(1) { I = rhs.I; } bool IsStencil() const { return true; } MeshCellOut *GetRef(int i, int j) const { if (mult_ != 1.) return 0; MeshCellIn *arg = rhs_.GetRef(i, j); return (arg && func_ref_) ? &func_ref_(*arg, !rhs_.IsStencil()) : 0; } virtual void Eval(MeshCellOut *res, int i, int j) const { MeshCellIn *arg = rhs_.GetRef(i, j); if (arg) { *res = func_copy_ ? func_copy_(*arg, !rhs_.IsStencil()) : func_ref_(*arg, !rhs_.IsStencil()); } else { MeshCellIn &arg = rhs_.Eval(i, j); *res = func_copy_ ? func_copy_(arg, !rhs_.IsStencil()) : func_ref_(arg, !rhs_.IsStencil()); } *res *= mult_; } MeshOperatorBind &operator - () { mult_ = -1; return *this; } private: //MeshProxy     Eval const MeshOperator<MeshCellIn> &rhs_; GlobalCellFunc_copy func_copy_; GlobalCellFunc_ref func_ref_; double mult_; }; //================================================================================== template<class MeshDeriv, class MeshFunc, class MeshFuncProxy, class MeshCellIn, class MeshCellOut> class MeshProxy : public MeshOperator<MeshCellOut> { public: MeshFunc &mesh; MeshProxy(MeshFunc &_mesh, const Interval2d &_I, double mult = 1) : mesh(_mesh), mult_(mult) { I = _I; } virtual ~MeshProxy() {} virtual MeshDeriv *This() = 0; int CountOfDependentVars() const { return 1; } MeshCellOut *GetRef(int i, int j) const { return (1. == mult_) ? &mesh(i, j) : 0; } virtual void Eval(MeshCellOut *res, int i, int j) const { *res = mesh(i, j); if (mult_ != 1.) *res *= mult_; } bool IsEqual(const MeshProxy &rhs) const { if (I.lb_x != rhs.I.lb_x) return 0; if (I.lb_y != rhs.I.lb_y) return 0; FOREACH_INTERVAL(I, i, j) { if (mesh_(i, j) != rhs.mesh_(i, j) || mult_ != rhs.mult_) return 0; } return 1; } MeshDeriv &operator - () { mult_ = -1; return *This(); } MeshDeriv &operator =(const MeshOperator<MeshCellOut> &rhs) { FOREACH_INTERVAL(I, i, j) { MeshCellOut &c = mesh(i, j); if (rhs.CountOfDependentVars() <= 1) rhs.Eval(&c, i, j); else c = rhs.Eval(i, j); } MeshCellStack<MeshCellOut>::Instance().cells.clear(); return *This(); } MeshDeriv &operator +=(const MeshOperator<MeshCellOut> &rhs) { FOREACH_INTERVAL(I, i, j) { MeshCellOut &c = mesh(i, j); c += rhs.Eval(i, j); } return *This(); } MeshDeriv &operator -=(const MeshOperator<MeshCellOut> &rhs) { FOREACH_INTERVAL(I, i, j) { MeshCellOut &c = mesh(i, j); c -= rhs.Eval(i, j); } return *This(); } MeshDeriv &operator *=(const MeshOperator<MeshCellOut> &rhs) { FOREACH_INTERVAL(I, i, j) { MeshCellOut &c = mesh(i, j); c *= rhs.Eval(i, j); } return *This(); } MeshDeriv &operator /=(const MeshOperator<MeshCellOut> &rhs) { FOREACH_INTERVAL(I, i, j) { MeshCellOut &c = mesh(i, j); c /= rhs.Eval(i, j); } return *This(); } #define DECLARE_MESH_OP(OPERAND)\ MeshDeriv &operator OPERAND(double val) {\ FOREACH_INTERVAL(I, i, j) {\ mesh(i, j) OPERAND val;\ }\ return *This();\ } DECLARE_MESH_OP( = ) DECLARE_MESH_OP( += ) DECLARE_MESH_OP( -= ) DECLARE_MESH_OP( *= ) DECLARE_MESH_OP( /= ) #undef DECLARE_MESH_OP void ToVector(dblArray1d *v, int lb = 0) {//todo: template T const int nx = I.nx(); const int n = lb + nx * I.ny(); if (v->size() != n) v->resize(n); int at = lb; FOREACH_INTERVAL(I, i, j) { (*v)[at++] = mesh(i, j); } } void FromVector(const dblArray1d &v, int lb = 0) {//todo: template T if (I.nx() * I.ny() != v.size() - lb) ;//THROW("Size of mesh must and (v.size() - lb) must be equals")//todo int at = lb; FOREACH_INTERVAL(I, i, j) { mesh(i, j) = v[at++]; } } private: double mult_; }; //================================================================================== template<class MeshDeriv, class MeshCell, class MeshProxy, class MeshScalarOperator, class MeshVectorOperator> class MeshFunc : public array2d<MeshCell> { public: double a_; double b_; double hx_; double hy_; MeshFunc(double a = 0, double b = 0, int nx = 0, int ny = 0) { if (nx > 0 && ny > 0) Resize(a, b, nx, ny); } virtual ~MeshFunc() {} void Resize(double a, double b, int nx, int ny); double a() const { return a_; } double b() const { return b_; } int nx() const { return nx_ - 2; } int ny() const { return ny_ - 2; } double hx() const { return hx_; } double hy() const { return hy_; } double x(int i, int j) { return at(i, j).x(); } double y(int i, int j) { return at(i, j).y(); } Interval2d I_all() const { return Interval2d(0, nx_ - 1, 0, ny_ - 1); } Interval2d I_int() const { return Interval2d(1, nx_ - 2, 1, ny_ - 2); } virtual MeshDeriv *This() = 0; operator MeshProxy() { return without_bnd(); } MeshProxy without_bnd() { return MeshProxy(*This(), I_int()); } MeshProxy all() { return MeshProxy(*This(), I_all()); } MeshDeriv &operator =(MeshDeriv &rhs) { all() = rhs.all(); return *This(); } MeshDeriv &operator +=(MeshDeriv &rhs) { all() += rhs.all(); return *This(); } MeshDeriv &operator -=(MeshDeriv &rhs) { all() -= rhs.all(); return *This(); } MeshDeriv &operator *=(MeshDeriv &rhs) { all() *= rhs.all(); return *This(); } MeshDeriv &operator /=(MeshDeriv &rhs) { all() /= rhs.all(); return *This(); } MeshProxy dx_left(const Interval2d &I) { return ::dx_left(MeshProxy(*this, I)); } MeshProxy dx_right(const Interval2d &I) { return ::dx_right(MeshProxy(*this, I)); } MeshProxy laplacian(const Interval2d &I) { return ::laplacian(MeshProxy(*this, I)); } MeshScalarOperator &div_left(const Interval2d &I) { return ::div_left(MeshProxy(*this, I)); } MeshScalarOperator &div_right(const Interval2d &I) { return ::div_right(MeshProxy(*this, I)); } MeshVectorOperator &grad_left(const Interval2d &I) { return ::grad_left(MeshProxy(*this, I)); } MeshVectorOperator &grad_right(const Interval2d &I) { return ::grad_right(MeshProxy(*this, I)); } }; template<class MeshDeriv, class MeshCell, class MeshProxy, class MeshScalarOperator, class MeshVectorOperator> void MeshFunc<MeshDeriv, MeshCell, MeshProxy, MeshScalarOperator, MeshVectorOperator>::Resize(double a, double b, int nx, int ny) { a_ = a; b_ = b; array2d<MeshCell>::Resize(nx + 2, ny + 2); hx_ = (b_ - a_) / (nx + 1); hy_ = (b_ - a_) / (ny + 1); for(int i = 0, j; i <= nx + 1; i++) { for(j = 0; j <= ny + 1; j++) { MeshCell &cell = at(i, j); cell.set_x(a_ + i * hx_); cell.set_y(a_ + j * hy_); if (i > 0) cell.set_left(&at(i - 1, j)); else cell.set_left(0);//,    Resize   ,   if (i <= nx) cell.set_right(&at(i + 1, j)); else cell.set_right(0); if (j > 0) cell.set_down(&at(i, j - 1)); else cell.set_down(0); if (j <= ny) cell.set_up(&at(i, j + 1)); else cell.set_up(0); } } } //================================================================================== template<class T> class ScalarMeshFunc; template<class T> class ScalarMeshProxy : public MeshProxy<ScalarMeshProxy<T>, ScalarMeshFunc<T>, ScalarMeshProxy<T>, ScalarCell<T>, ScalarCell<T> > { public: ScalarMeshProxy(ScalarMeshFunc<T> &mesh, const Interval2d &_I) : MeshProxy(mesh, _I) {} ScalarMeshProxy *This() { return this; } ScalarMeshProxy &operator=(const ScalarMeshProxy &rhs) { DASSERT_GE(I.lb_x, rhs.I.lb_x); DASSERT_GE(I.lb_y, rhs.I.lb_y); DASSERT_LE(I.ub_x, rhs.I.ub_x); DASSERT_LE(I.ub_y, rhs.I.ub_y); return MeshProxy::operator =(rhs); } ScalarMeshProxy &operator=(const MeshOperator<ScalarCell<T> > &rhs) { return MeshProxy::operator =(rhs); } ScalarMeshProxy &operator=(const T &val) { return MeshProxy::operator =(val); } }; template<class T> class Scalar2VectorMeshProxy : public MeshProxy<Scalar2VectorMeshProxy<T>, ScalarMeshFunc<T>, ScalarMeshProxy<T>, ScalarCell<T>, VectorCell<T> > { public: Scalar2VectorMeshProxy *This() { return this; } Scalar2VectorMeshProxy &operator=(const Scalar2VectorMeshProxy &rhs) { return MeshProxy::operator =(rhs); } }; //================================================================================== template<class T> class VectorMeshFunc; template<class T> class ScalarMeshFunc : public MeshFunc<ScalarMeshFunc<T>, ScalarCell<T>, ScalarMeshProxy<T>, ScalarMeshFunc<T>, VectorMeshFunc<T> > { public: ScalarMeshFunc(double a = 0, double b = 0, int nx = 0, int ny = 0) : MeshFunc(a, b, nx, ny) {} ScalarMeshFunc *This() { return this; } ScalarCell<T> &operator()(int i, int j) { return at(i, j); } const ScalarCell<T> &operator()(int i, int j) const { return at(i, j); } ScalarMeshProxy<T> operator()(const Interval2d &I) { return ScalarMeshProxy<T>(*this, I); } }; template<class T> std::ostream &operator << (std::ostream &os, const MeshOperator<T> &mesh) { const Interval2d &I = mesh.I; StreamTable st; st.SetCols(I.ub_y - I.lb_y + 1, 10); FOREACH_INTERVAL(I, i, j) { st << mesh.Eval(i, j); } //os << st.os(); return os; } //================================================================================== template<class T> class VectorMeshProxy : public MeshProxy<VectorMeshProxy<T>, VectorMeshFunc<T>, VectorMeshProxy<T>, VectorCell<T>, VectorCell<T> > { public: VectorMeshProxy(VectorMeshFunc<T> &mesh, const Interval2d &_I) : MeshProxy(mesh, _I) {} VectorMeshProxy *This() { return this; } VectorMeshProxy &operator=(const VectorMeshProxy &rhs) { DASSERT_GE(I.lb_x, rhs.I.lb_x); DASSERT_GE(I.lb_y, rhs.I.lb_y); DASSERT_LE(I.ub_x, rhs.I.ub_x); DASSERT_LE(I.ub_y, rhs.I.ub_y); return MeshProxy::operator =(rhs); } VectorMeshProxy &operator=(const MeshOperator<VectorCell<T> > &rhs) { return MeshProxy::operator =(rhs); } VectorMeshProxy &operator=(const T &val) { return MeshProxy::operator =(val); } }; template<class T> class Vector2ScalarMeshProxy : public MeshProxy<Vector2ScalarMeshProxy<T>, VectorMeshFunc<T>, VectorMeshProxy<T>, VectorCell<T>, ScalarCell<T> > { public: Vector2ScalarMeshProxy *This() { return this; } Vector2ScalarMeshProxy &operator=(const Vector2ScalarMeshProxy &rhs) { return MeshProxy::operator =(rhs); } }; /*template<class T> class Vector2VectorMeshProxy : public MeshProxy<Vector2VectorMeshProxy<T>, VectorMeshFunc<T>, VectorMeshProxy<T>, VectorCell<T>, VectorCell<VectorCell<T> > > { public: Vector2VectorMeshProxy *This() { return this; } Vector2VectorMeshProxy &operator=(const Vector2VectorMeshProxy &rhs) { return MeshProxy::operator =(rhs); } };*/ template<class T> class VectorMeshFunc : public MeshFunc<VectorMeshFunc<T>, VectorCell<T>, VectorMeshProxy<T>, ScalarMeshFunc<T>, VectorMeshFunc<T> > { public: VectorMeshFunc(double a = 0, double b = 0, int nx = 0, int ny = 0) : MeshFunc(a, b, nx, ny) {} VectorMeshFunc *This() { return this; } VectorCell<T> &operator()(int i, int j) { return at(i, j); } const VectorCell<T> &operator()(int i, int j) const { return at(i, j); } VectorMeshProxy<T> operator()(const Interval2d &I) { return VectorMeshProxy<T>(*this, I); } }; //================================================================================== #define DECLARE_MESH_OP(CLASS,OPERAND,OPFUNC)\ template<class MeshCell>\ class CLASS : public MeshOperator<MeshCell> {\ public:\ MeshCell &res_;\ CLASS(const MeshOperator &op1, const MeshOperator &op2, MeshCell &res) : op1_(op1), op2_(op2), res_(res) {\ I = op1.I;/*todo*/\ }\ \ bool IsStencil() const {\ return op1_.IsStencil() || op2_.IsStencil();\ }\ int CountOfDependentVars() const {\ return op1_.CountOfDependentVars() + op2_.CountOfDependentVars() + 1;\ }\ virtual void Eval(MeshCell *res, int i, int j) const {\ *res = op1_.Eval(i, j);\ *res OPERAND= op2_.Eval(i, j);\ }\ virtual MeshCell &Eval(int i, int j) const {\ res_ = op1_.Eval(i, j);\ res_ OPERAND= op2_.Eval(i, j);\ return res_;\ }\ private:\ const MeshOperator &op1_;\ const MeshOperator &op2_;\ };\ inline CLASS<ScalarCell<double> > operator OPERAND (const MeshOperator<ScalarCell<double> > &op1, const MeshOperator<ScalarCell<double> > &op2)\ {\ std::list<ScalarCell<double> > &stack = MeshCellStack<ScalarCell<double> >::Instance().cells;\ stack.push_back(ScalarCell<double>());\ return CLASS<ScalarCell<double> >(op1, op2, stack.back());\ }\ inline CLASS<VectorCell<double> > operator OPERAND (const MeshOperator<VectorCell<double> > &op1, const MeshOperator<VectorCell<double> > &op2)\ {\ std::list<VectorCell<double> > &stack = MeshCellStack<VectorCell<double> >::Instance().cells;\ stack.push_back(VectorCell<double>());\ return CLASS<VectorCell<double> >(op1, op2, stack.back());\ } DECLARE_MESH_OP(OperatorPlus, +, PlusEq()) DECLARE_MESH_OP(OperatorMinus, -, MinusEq()) DECLARE_MESH_OP(OperatorMult, *, MultEq()) DECLARE_MESH_OP(OperatorDivide, / , DivideEq()) #undef DECLARE_MESH_OP //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> fabs(ScalarCell<T> &u, bool recur) { ScalarCell<T> res(u); res.val() = ::fabs(res.val()); return res; } template<class T> VectorCell<T> fabs(VectorCell<T> &u, bool recur) { VectorCell<T> res(u); for(int i = 0; i < res.dim_; i++) { res.comp(i).val() = fabs(res.comp(i).val()); } return res; } #define SCALAR_OPERATOR(T) MeshOperatorBind<ScalarCell<T>, ScalarCell<T> > #define VECTOR_OPERATOR(T) MeshOperatorBind<VectorCell<T>, VectorCell<T> > #define SCALAR_2_VECTOR_OPERATOR(T) MeshOperatorBind<ScalarCell<T>, VectorCell<T> > #define VECTOR_2_SCALAR_OPERATOR(T) MeshOperatorBind<VectorCell<T>, ScalarCell<T> > template<class T> SCALAR_OPERATOR(T) fabs(MeshOperator<ScalarCell<T> > &f) { return SCALAR_OPERATOR(T)(f, &::fabs); } template<class T> VECTOR_OPERATOR(T) fabs(MeshOperator<VectorCell<T> > &f) { return VECTOR_OPERATOR(T)(f, &::fabs); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &dx_left(ScalarCell<T> &u, bool recur = true) { return u.dx_left(recur); } template<class T> VectorCell<T> &dx_left(VectorCell<T> &u, bool recur = true) { return u.dx_left(recur); } template<class T> SCALAR_OPERATOR(T) dx_left(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_OPERATOR(T)(f, &::dx_left); } template<class T> VECTOR_OPERATOR(T) dx_left(const MeshOperator<VectorCell<T> > &f) { return VECTOR_OPERATOR(T)(f, &::dx_left); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &dx_right(ScalarCell<T> &u, bool recur = true) { return u.dx_right(recur); } template<class T> VectorCell<T> &dx_right(VectorCell<T> &u, bool recur = true) { return u.dx_right(recur); } template<class T> SCALAR_OPERATOR(T) dx_right(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_OPERATOR(T)(f, &::dx_right); } template<class T> VECTOR_OPERATOR(T) dx_right(const MeshOperator<VectorCell<T> > &f) { return VECTOR_OPERATOR(T)(f, &::dx_right); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &dy_left(ScalarCell<T> &u, bool recur = true) { return u.dy_left(recur); } template<class T> VectorCell<T> &dy_left(VectorCell<T> &u, bool recur = true) { return u.dy_left(recur); } template<class T> SCALAR_OPERATOR(T) dy_left(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_OPERATOR(T)(f, &::dy_left); } template<class T> VECTOR_OPERATOR(T) dy_left(const MeshOperator<VectorCell<T> > &f) { return VECTOR_OPERATOR(T)(f, &::dy_left); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &dy_right(ScalarCell<T> &u, bool recur = true) { return u.dy_right(recur); } template<class T> VectorCell<T> &dy_right(VectorCell<T> &u, bool recur = true) { return u.dy_right(recur); } template<class T> SCALAR_OPERATOR(T) dy_right(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_OPERATOR(T)(f, &::dy_right); } template<class T> VECTOR_OPERATOR(T) dy_right(const MeshOperator<VectorCell<T> > &f) { return VECTOR_OPERATOR(T)(f, &::dy_right); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &laplacian(ScalarCell<T> &u, bool recur = true) { return u.laplacian(recur); } template<class T> VectorCell<T> &laplacian(VectorCell<T> &u, bool recur = true) { return u.laplacian(recur); } template<class T> SCALAR_OPERATOR(T) laplacian(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_OPERATOR(T)(f, &::laplacian); } template<class T> VECTOR_OPERATOR(T) laplacian(const MeshOperator<VectorCell<T> > &f) { return VECTOR_OPERATOR(T)(f, &::laplacian); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &div_left(VectorCell<T> &v, bool recur = true) { return v.div_left(recur); } template<class T> VECTOR_2_SCALAR_OPERATOR(T) div_left(const MeshOperator<VectorCell<T> > &f) { return VECTOR_2_SCALAR_OPERATOR(T)(f, &::div_left); } //---------------------------------------------------------------------------------- template<class T> ScalarCell<T> &div_right(VectorCell<T> &v, bool recur = true) { return v.div_right(recur); } template<class T> VECTOR_2_SCALAR_OPERATOR(T) div_right(const MeshOperator<VectorCell<T> > &f) { return VECTOR_2_SCALAR_OPERATOR(T)(f, &::div_right); } //---------------------------------------------------------------------------------- template<class T> VectorCell<T> &grad_left(ScalarCell<T> &u, bool recur = true) { return u.grad_left(recur); } template<class T> SCALAR_2_VECTOR_OPERATOR(T) grad_left(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_2_VECTOR_OPERATOR(T)(f, &::grad_left); } /*template<class T> Vector2VectorMeshProxy<T> grad_left(VECTOR_OPERATOR(T) &f) { return Vector2VectorMeshProxy<T>(f, &::grad_left); }*/ //---------------------------------------------------------------------------------- template<class T> VectorCell<T> &grad_right(ScalarCell<T> &u, bool recur = true) { return u.grad_right(recur); } template<class T> SCALAR_2_VECTOR_OPERATOR(T) grad_right(const MeshOperator<ScalarCell<T> > &f) { return SCALAR_2_VECTOR_OPERATOR(T)(f, &::grad_right); } /* template<class T> Vector2VectorMeshProxy<T> grad_right(VECTOR_OPERATOR(T) &f) { return Vector2VectorMeshProxy<T>(f, &::grad_right); }*/ //================================================================================== template<class T> const T &max(const VectorCell<T> &c) { return c.max(); } template<class T> const T &min(const VectorCell<T> &c) { return c.min(); } template<class T> T &max(const MeshOperator<ScalarCell<T> > &m) { return m.max(); } template<class T> T &min(const MeshOperator<ScalarCell<T> > &m) { return m.min(); } template<class T> VectorCell<T> &max(const MeshOperator<VectorCell<T> > &m) { return m.max(); } template<class T> VectorCell<T> &min(const MeshOperator<VectorCell<T> > &m) { return m.min(); } 

 //() 2016-2019  , SharipovTR@gmail.com // As Is #ifndef __ARRAYS_H #define __ARRAYS_H #pragma once #include "NumericAssert.h" #include <vector> /** * Dynamic 0-based 1D array (storage is std::vector) */ template<class T> class array1d : public std::vector<T> { public: typedef std::vector<T> baseClass; typedef std::size_t index; array1d(index n = 0) : baseClass(n) {} const T &operator()(index i) const { return baseClass::operator[](i); } T &operator()(index i) { return baseClass::operator[](i); } void Resize(index /*lb*/, index ub) { baseClass::resize(ub + 1); } }; #define RS2D(i,j,ny) ((i) * (ny) + (j)) /** * Dynamic 0-based row-oriented 2D array (storage is only one std::vector) */ template<class T> class array2d : public std::vector<T> { public: typedef std::vector<T> baseClass; typedef std::size_t index; array2d(index nx = 0, index ny = 0) { Resize(nx, ny); } array2d(int n, const T v[]) { Resize(n, n); std::copy(v, v + n * n, baseClass::begin()); } array2d(const array2d<T> &rhs) { Resize(rhs.nx_, rhs.ny_); std::copy(rhs.begin(), rhs.end(), begin()); } void Resize(index nx, index ny) { if (nx > 0 && ny > 0) baseClass::resize(nx * ny); nx_ = nx; ny_ = ny; } void resize(int nx, int ny, bool reserved) { Resize(nx, ny); } index Size() const { DASSERT_EQ(nx_, ny_); return nx_; } void Clear() { baseClass::clear(); nx_ = 0; ny_ = 0; } const T &at(index i, index j) const { DASSERT(CHR(i, 0, nx_)); DASSERT(CHR(j, 0, ny_)); return baseClass::operator[](RS2D(i, j, ny_)); } virtual T &at(index i, index j) { return baseClass::operator[](RS2D(i, j, ny_)); } const T &operator()(index i, index j) const { return at(i, j); } T &operator()(index i, index j) { return at(i, j); } void FillRow(index irow, T val) { DASSERT(CHR(irow, 0, nx_)); typename baseClass::iterator beg = baseClass::begin() + RS2D(irow, 0, ny_); std::fill(beg, beg + ny_, val); } std::vector<T> &operator =(const std::vector<T> &rhs) { std::copy(rhs.begin(), rhs.end(), baseClass::begin()); return *this; } public: index nx_; index ny_; }; 

 //() 2016-2019  , SharipovTR@gmail.com // As Is #pragma once #include "assert.h" #include <sstream> #include <iostream> #define TO_STR(x) #x #define ASSERT_OP(opname,op,val1,val2)\ if (!(val1 op val2)) do {\ std::stringstream ss;\ ss << TO_STR(ASSERT##opname) << " failed:" << "\n";\ ss << #val1 << ": " << val1 << "\n";\ ss << #val2 << ": " << val2 << "\n";\ ss << TO_STR(val1 - val2) << ": " << val1 - (val2);\ std::string s(ss.str());\ std::cout << s.c_str();\ } while(0) #define ASSERT_EQ(val1, val2) ASSERT_OP(_EQ, ==, val1, val2) #define ASSERT_NE(val1, val2) ASSERT_OP(_NE, !=, val1, val2) #define ASSERT_LE(val1, val2) ASSERT_OP(_LE, <=, val1, val2) #define ASSERT_LT(val1, val2) ASSERT_OP(_LT, < , val1, val2) #define ASSERT_GE(val1, val2) ASSERT_OP(_GE, >=, val1, val2) #define ASSERT_GT(val1, val2) ASSERT_OP(_GT, > , val1, val2) #define ASSERT_NEAR(val1,val2,margin)\ do {\ ASSERT_LE((val1), (val2) + margin);\ ASSERT_GE((val1), (val2) - margin);\ } while(0) #define ASSERT_DBL_EQ(val1,val2)\ ASSERT_NEAR(val1, val2, 1.e-12); #ifdef _DEBUG #define DASSERT assert #define DASSERT_EQ(val1, val2) ASSERT_EQ(val1, val2) #define DASSERT_NE(val1, val2) ASSERT_NE(val1, val2) #define DASSERT_LE(val1, val2) ASSERT_LE(val1, val2) #define DASSERT_LT(val1, val2) ASSERT_LT(val1, val2) #define DASSERT_GE(val1, val2) ASSERT_GE(val1, val2) #define DASSERT_GT(val1, val2) ASSERT_GT(val1, val2) #define DASSERT_NEAR(val1,val2,margin) ASSERT_NEAR(val1,val2,margin) #define DASSERT_NEAR_SMALL(val1,val2) ASSERT_NEAR_SMALL(val1,val2) #define DASSERT_DBL_EQ(val1,val2) ASSERT_DBL_EQ(val1,val2) #else #define DASSERT #define DASSERT_EQ(val1, val2) #define DASSERT_NE(val1, val2) #define DASSERT_LE(val1, val2) #define DASSERT_LT(val1, val2) #define DASSERT_GE(val1, val2) #define DASSERT_GT(val1, val2) #define DASSERT_NEAR(val1,val2,margin) #define DASSERT_NEAR_SMALL(val1,val2) #define DASSERT_DBL_EQ(val1,val2) #endif