Компактная С++ библиотека для программирования конечно-разностных методов в операторном стиле. Часть 1. Семантика01.03.2019 23:02
//(с) 2016-2019 Шарипов Тимур, SharipovTR@gmail.com
//Распространяется As Is
#pragma once
#include
#include
#include "NumericAssert.h"
#include "Arrays.h"
//#define USE_STD_MOVE
template
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 scal_[7];//todo: через shared_ptr
std::auto_ptr 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
DerivCell &
MeshCell::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
DerivCell &
MeshCell::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
DerivCell &
MeshCell::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
DerivCell &
MeshCell::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
DerivCell &
MeshCell::laplacian(bool recur) {
DerivCell &res = Instance(4);
res = dx_left(recur).dx_right(false) + dy_left(recur).dy_right(false);//todo
return res;
}
template
ScalarCell &
MeshCell::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
ScalarCell &
MeshCell::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
VectorCell &
MeshCell::grad_left(bool recur) {
VectorCell &res = InstanceAsVector(0);
res.comp(0) = dx_left(recur);
res.comp(1) = dy_left(recur);
return res;
}
template
VectorCell &
MeshCell::grad_right(bool recur) {
VectorCell &res = InstanceAsVector(1);
res.comp(0) = dx_right(recur);
res.comp(1) = dy_right(recur);
return res;
}
template
class VectorCell;
template
class ScalarCell : public MeshCell, ScalarCell, VectorCell > {
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 &comp(int i) {
DASSERT_LE(i, 1);
return *this;
}
const ScalarCell &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 &res = scal_[i];
if (!res.get())
res.reset(new ScalarCell(*this));
return *res;
}
ScalarCell &InstanceAsScalar(int i) {
return Instance(i);
}
VectorCell &InstanceAsVector(int i) {
std::auto_ptr > &res = vec_[i];
if (!res.get())
res.reset(new VectorCell);
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 VectorCell : public MeshCell, ScalarCell, VectorCell > {
public:
std::vector *> comp_;//todo->auto_ptr
VectorCell() {
Resize(2);
for(int i = 0; i < dim_; i++) {
comp_[i] = new ScalarCell;
}
}
explicit VectorCell(const T &val) {
Resize(2);
for(int i = 0; i < dim_; i++) {
comp_[i] = new ScalarCell(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(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 &comp(int i) {
return *comp_[i];
}
const ScalarCell &comp(int i) const {
return *comp_[i];
}
VectorCell &Instance(int i) {
std::auto_ptr &res = vec_[i];
if (!res.get())
res.reset(new VectorCell(*this));
return *res;
}
ScalarCell &InstanceAsScalar(int i) {
std::auto_ptr > &res = scal_[i];
if (!res.get())
res.reset(new ScalarCell);
res->set_x(x_);
res->set_y(y_);
return *res;
}
VectorCell &InstanceAsVector(int i) {
return Instance(i);
}
ScalarCell &operator()(int at) {
return comp(at);
}
const ScalarCell &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\
ScalarCell operator OPERAND (const ScalarCell &c1, const ScalarCell &c2)\
{\
ScalarCell res(c1);\
res OPERAND= c2;\
return res;\
}\
template\
VectorCell operator OPERAND (const VectorCell &c1, const VectorCell &c2)\
{\
VectorCell res(c1);\
res OPERAND= c2;\
return res;\
}\
template\
VectorCell operator OPERAND (const VectorCell &c, double val)\
{\
VectorCell res(c);\
res OPERAND= val;\
return res;\
}\
template\
VectorCell operator OPERAND (double val, const VectorCell &c)\
{\
VectorCell res(c);\
res OPERAND= val;\
return res;\
}
DECLARE_CELL_OP(+)
DECLARE_CELL_OP(-)
DECLARE_CELL_OP(*)
DECLARE_CELL_OP( / )
#undef DECLARE_CELL_OP
//----------------------------------------------------------------------------------
template
bool operator < (const VectorCell &lhs, const VectorCell &rhs) {
for (int i = 0; i < lhs.dim_; i++)
if (lhs(i) >= rhs(i))
return 0;
return 1;
}
//----------------------------------------------------------------------------------
template
bool operator <= (const VectorCell &lhs, const VectorCell &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 MeshCellStack {
public:
static MeshCellStack &Instance() {
static MeshCellStack instance;
return instance;
}
std::list cells;
};
template
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 MeshOperatorBind : public MeshOperator {
public:
typedef MeshCellOut &(*GlobalCellFunc_ref)(MeshCellIn &, bool recur);
typedef MeshCellOut (*GlobalCellFunc_copy)(MeshCellIn &, bool recur);
MeshOperatorBind(const MeshOperator &rhs, GlobalCellFunc_copy func) : rhs_(rhs), func_copy_(func), func_ref_(0), mult_(1) {
I = rhs.I;
}
MeshOperatorBind(const MeshOperator &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 &rhs_;
GlobalCellFunc_copy func_copy_;
GlobalCellFunc_ref func_ref_;
double mult_;
};
//==================================================================================
template
class MeshProxy : public MeshOperator {
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 &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::Instance().cells.clear();
return *This();
}
MeshDeriv &operator +=(const MeshOperator &rhs) {
FOREACH_INTERVAL(I, i, j) {
MeshCellOut &c = mesh(i, j);
c += rhs.Eval(i, j);
}
return *This();
}
MeshDeriv &operator -=(const MeshOperator &rhs) {
FOREACH_INTERVAL(I, i, j) {
MeshCellOut &c = mesh(i, j);
c -= rhs.Eval(i, j);
}
return *This();
}
MeshDeriv &operator *=(const MeshOperator &rhs) {
FOREACH_INTERVAL(I, i, j) {
MeshCellOut &c = mesh(i, j);
c *= rhs.Eval(i, j);
}
return *This();
}
MeshDeriv &operator /=(const MeshOperator &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 MeshFunc : public array2d {
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
void
MeshFunc::Resize(double a, double b, int nx, int ny) {
a_ = a;
b_ = b;
array2d::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 ScalarMeshFunc;
template
class ScalarMeshProxy : public MeshProxy, ScalarMeshFunc, ScalarMeshProxy, ScalarCell, ScalarCell > {
public:
ScalarMeshProxy(ScalarMeshFunc &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 > &rhs) {
return MeshProxy::operator =(rhs);
}
ScalarMeshProxy &operator=(const T &val) {
return MeshProxy::operator =(val);
}
};
template
class Scalar2VectorMeshProxy : public MeshProxy, ScalarMeshFunc, ScalarMeshProxy, ScalarCell, VectorCell > {
public:
Scalar2VectorMeshProxy *This() {
return this;
}
Scalar2VectorMeshProxy &operator=(const Scalar2VectorMeshProxy &rhs) {
return MeshProxy::operator =(rhs);
}
};
//==================================================================================
template
class VectorMeshFunc;
template
class ScalarMeshFunc
: public MeshFunc, ScalarCell, ScalarMeshProxy, ScalarMeshFunc, VectorMeshFunc > {
public:
ScalarMeshFunc(double a = 0, double b = 0, int nx = 0, int ny = 0)
: MeshFunc(a, b, nx, ny) {}
ScalarMeshFunc *This() {
return this;
}
ScalarCell &operator()(int i, int j) {
return at(i, j);
}
const ScalarCell &operator()(int i, int j) const {
return at(i, j);
}
ScalarMeshProxy operator()(const Interval2d &I) {
return ScalarMeshProxy(*this, I);
}
};
template
std::ostream &operator << (std::ostream &os, const MeshOperator &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 VectorMeshProxy : public MeshProxy, VectorMeshFunc, VectorMeshProxy, VectorCell, VectorCell > {
public:
VectorMeshProxy(VectorMeshFunc &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 > &rhs) {
return MeshProxy::operator =(rhs);
}
VectorMeshProxy &operator=(const T &val) {
return MeshProxy::operator =(val);
}
};
template
class Vector2ScalarMeshProxy : public MeshProxy, VectorMeshFunc, VectorMeshProxy, VectorCell, ScalarCell > {
public:
Vector2ScalarMeshProxy *This() {
return this;
}
Vector2ScalarMeshProxy &operator=(const Vector2ScalarMeshProxy &rhs) {
return MeshProxy::operator =(rhs);
}
};
/*template
class Vector2VectorMeshProxy : public MeshProxy, VectorMeshFunc, VectorMeshProxy, VectorCell, VectorCell > > {
public:
Vector2VectorMeshProxy *This() {
return this;
}
Vector2VectorMeshProxy &operator=(const Vector2VectorMeshProxy &rhs) {
return MeshProxy::operator =(rhs);
}
};*/
template
class VectorMeshFunc
: public MeshFunc, VectorCell, VectorMeshProxy, ScalarMeshFunc, VectorMeshFunc > {
public:
VectorMeshFunc(double a = 0, double b = 0, int nx = 0, int ny = 0)
: MeshFunc(a, b, nx, ny) {}
VectorMeshFunc *This() {
return this;
}
VectorCell &operator()(int i, int j) {
return at(i, j);
}
const VectorCell &operator()(int i, int j) const {
return at(i, j);
}
VectorMeshProxy operator()(const Interval2d &I) {
return VectorMeshProxy(*this, I);
}
};
//==================================================================================
#define DECLARE_MESH_OP(CLASS,OPERAND,OPFUNC)\
template\
class CLASS : public MeshOperator {\
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 > operator OPERAND (const MeshOperator > &op1, const MeshOperator > &op2)\
{\
std::list > &stack = MeshCellStack >::Instance().cells;\
stack.push_back(ScalarCell());\
return CLASS >(op1, op2, stack.back());\
}\
inline CLASS > operator OPERAND (const MeshOperator > &op1, const MeshOperator > &op2)\
{\
std::list > &stack = MeshCellStack >::Instance().cells;\
stack.push_back(VectorCell());\
return CLASS >(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
ScalarCell fabs(ScalarCell &u, bool recur) {
ScalarCell res(u);
res.val() = ::fabs(res.val());
return res;
}
template
VectorCell fabs(VectorCell &u, bool recur) {
VectorCell 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 >
#define VECTOR_OPERATOR(T) MeshOperatorBind, VectorCell >
#define SCALAR_2_VECTOR_OPERATOR(T) MeshOperatorBind, VectorCell >
#define VECTOR_2_SCALAR_OPERATOR(T) MeshOperatorBind, ScalarCell >
template
SCALAR_OPERATOR(T) fabs(MeshOperator > &f) {
return SCALAR_OPERATOR(T)(f, &::fabs);
}
template
VECTOR_OPERATOR(T) fabs(MeshOperator > &f) {
return VECTOR_OPERATOR(T)(f, &::fabs);
}
//----------------------------------------------------------------------------------
template
ScalarCell &dx_left(ScalarCell &u, bool recur = true) {
return u.dx_left(recur);
}
template
VectorCell &dx_left(VectorCell &u, bool recur = true) {
return u.dx_left(recur);
}
template
SCALAR_OPERATOR(T) dx_left(const MeshOperator > &f) {
return SCALAR_OPERATOR(T)(f, &::dx_left);
}
template
VECTOR_OPERATOR(T) dx_left(const MeshOperator > &f) {
return VECTOR_OPERATOR(T)(f, &::dx_left);
}
//----------------------------------------------------------------------------------
template
ScalarCell &dx_right(ScalarCell &u, bool recur = true) {
return u.dx_right(recur);
}
template
VectorCell &dx_right(VectorCell &u, bool recur = true) {
return u.dx_right(recur);
}
template
SCALAR_OPERATOR(T) dx_right(const MeshOperator > &f) {
return SCALAR_OPERATOR(T)(f, &::dx_right);
}
template
VECTOR_OPERATOR(T) dx_right(const MeshOperator > &f) {
return VECTOR_OPERATOR(T)(f, &::dx_right);
}
//----------------------------------------------------------------------------------
template
ScalarCell &dy_left(ScalarCell &u, bool recur = true) {
return u.dy_left(recur);
}
template
VectorCell &dy_left(VectorCell &u, bool recur = true) {
return u.dy_left(recur);
}
template
SCALAR_OPERATOR(T) dy_left(const MeshOperator > &f) {
return SCALAR_OPERATOR(T)(f, &::dy_left);
}
template
VECTOR_OPERATOR(T) dy_left(const MeshOperator > &f) {
return VECTOR_OPERATOR(T)(f, &::dy_left);
}
//----------------------------------------------------------------------------------
template
ScalarCell &dy_right(ScalarCell &u, bool recur = true) {
return u.dy_right(recur);
}
template
VectorCell &dy_right(VectorCell &u, bool recur = true) {
return u.dy_right(recur);
}
template
SCALAR_OPERATOR(T) dy_right(const MeshOperator > &f) {
return SCALAR_OPERATOR(T)(f, &::dy_right);
}
template
VECTOR_OPERATOR(T) dy_right(const MeshOperator > &f) {
return VECTOR_OPERATOR(T)(f, &::dy_right);
}
//----------------------------------------------------------------------------------
template
ScalarCell &laplacian(ScalarCell &u, bool recur = true) {
return u.laplacian(recur);
}
template
VectorCell &laplacian(VectorCell &u, bool recur = true) {
return u.laplacian(recur);
}
template
SCALAR_OPERATOR(T) laplacian(const MeshOperator > &f) {
return SCALAR_OPERATOR(T)(f, &::laplacian);
}
template
VECTOR_OPERATOR(T) laplacian(const MeshOperator > &f) {
return VECTOR_OPERATOR(T)(f, &::laplacian);
}
//----------------------------------------------------------------------------------
template
ScalarCell &div_left(VectorCell &v, bool recur = true) {
return v.div_left(recur);
}
template
VECTOR_2_SCALAR_OPERATOR(T) div_left(const MeshOperator > &f) {
return VECTOR_2_SCALAR_OPERATOR(T)(f, &::div_left);
}
//----------------------------------------------------------------------------------
template
ScalarCell &div_right(VectorCell &v, bool recur = true) {
return v.div_right(recur);
}
template
VECTOR_2_SCALAR_OPERATOR(T) div_right(const MeshOperator > &f) {
return VECTOR_2_SCALAR_OPERATOR(T)(f, &::div_right);
}
//----------------------------------------------------------------------------------
template