A compact C ++ library for programming operator-style finite difference methods. Part 1. Semantics
The semantics of the developed pde ++ library for programming finite-difference methods in the operator style is presented. The main objects of the library are a grid function, a grid cell and grid operators, the arithmetic relations between which bring the program code as close as possible to its mathematical notation. The pde ++ library is represented by just a few header files, has no external dependencies, and uses the concept of lazy computing.
A large number of mathematical modeling problems is reduced to the numerical solution of partial differential equations (PDEs) by grid methods. In the theory of difference schemes (Samarsky A.A.), the corresponding grid operators form a linear space over grid functions for which there is no direct representation in general-purpose programming languages such as C ++. As a result, the practice of recording the result of applying grid operators to grid functions using multidimensional arrays or matrices is widely used in software implementation.
Practice shows that the approach noted above is very useful in mastering the skills of implementing numerical methods, primarily due to its visibility when working with pre-written approximations of PDEs in index form. Significant problems also do not arise when extending this technique to generalized PDEs, if it is intended to implement a difference scheme with parameters once and reuse the appropriate program code without further modifications.
In the general case, a computing program can be modified in different directions, so the technique described above will require writing a significant amount of program code, and this, in turn, will increase the likelihood of typos and inconsistent recording of the same grid operators in different program modules. It is also worth noting the problem of duplication of program code with spatial dimensions variability (1D, 2D, 3D) and methods for approximating PDEs.
Thus, an alternative is the development and use of specialized software libraries with high-level domain abstractions that bring the program code closer to its mathematical notation. In the Blitz ++ library, such an abstraction is tensor computations on difference templates, implemented on the basis of the use of the template metaprogramming technique. FreePOOMA libraryexpands the Blitz ++ concept with difference analogues of differential divergence and gradient operators and the ability to work on multiprocessor computing systems. Unfortunately, these libraries have not been supported for a long time, possessing a number of limitations (they will be discussed in the next part) when they are used for fairly classical finite-difference approximations of the PDEs considered in this paper.
The open-source library pde ++, developed by the author, is ideally inspired by the freePOOMA library and is designed to record in the form of finite-difference schemes for scalar and vector grid functions defined in a 2D setting (1D and 3D in operation) on uniform rectangular grids.
Warning: the code was tested only under Windows.
A large number of mathematical modeling problems is reduced to the numerical solution of partial differential equations (PDEs) by grid methods. In the theory of difference schemes (Samarsky A.A.), the corresponding grid operators form a linear space over grid functions for which there is no direct representation in general-purpose programming languages such as C ++. As a result, the practice of recording the result of applying grid operators to grid functions using multidimensional arrays or matrices is widely used in software implementation.
Practice shows that the approach noted above is very useful in mastering the skills of implementing numerical methods, primarily due to its visibility when working with pre-written approximations of PDEs in index form. Significant problems also do not arise when extending this technique to generalized PDEs, if it is intended to implement a difference scheme with parameters once and reuse the appropriate program code without further modifications.
In the general case, a computing program can be modified in different directions, so the technique described above will require writing a significant amount of program code, and this, in turn, will increase the likelihood of typos and inconsistent recording of the same grid operators in different program modules. It is also worth noting the problem of duplication of program code with spatial dimensions variability (1D, 2D, 3D) and methods for approximating PDEs.
Thus, an alternative is the development and use of specialized software libraries with high-level domain abstractions that bring the program code closer to its mathematical notation. In the Blitz ++ library, such an abstraction is tensor computations on difference templates, implemented on the basis of the use of the template metaprogramming technique. FreePOOMA libraryexpands the Blitz ++ concept with difference analogues of differential divergence and gradient operators and the ability to work on multiprocessor computing systems. Unfortunately, these libraries have not been supported for a long time, possessing a number of limitations (they will be discussed in the next part) when they are used for fairly classical finite-difference approximations of the PDEs considered in this paper.
The open-source library pde ++, developed by the author, is ideally inspired by the freePOOMA library and is designed to record in the form of finite-difference schemes for scalar and vector grid functions defined in a 2D setting (1D and 3D in operation) on uniform rectangular grids.
Warning: the code was tested only under Windows.
#include "pdepp.h"
double sln_u(double x, double y)
{
return x * x + y * y;
}
//скалярные сеточные-функции
void test_pdepp_1()
{
//2d-область [a, b] x [a, b]
double a = 0;
double b = 1;
//число внутренних узлов равномерной сетки
int n = 10;
double h = (b - a) / (n + 1);
//инициализация сеточных функции (автоматически добавляется по одной граничной ячейке)
ScalarMeshFunc u(a, b, n, n);
ScalarMeshFunc r(a, b, n, n);
//поиндексный доступ к ячейкам
for (int i = 0; i <= n + 1; i++) {
for (int j = 0; j <= n + 1; j++) {
ScalarCell &c = u(i, j);
c.val() = sln_u(c.x(), c.y());
}
}
//сеточные операторы
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
ASSERT_DBL_EQ(dx_left(u(i, j)), (u(i, j) - u(i - 1, j)) / h);
ASSERT_DBL_EQ(dx_right(u(i, j)), (u(i + 1, j) - u(i, j)) / h);
ASSERT_DBL_EQ(dy_left(u(i, j)), (u(i, j) - u(i, j - 1)) / h);
ASSERT_DBL_EQ(dy_right(u(i, j)), (u(i, j + 1) - u(i, j)) / h);
ASSERT_DBL_EQ(laplacian(u(i, j)), dx_left(dx_right(u(i, j))) + dy_left(dy_right(u(i, j))));
ASSERT_DBL_EQ(laplacian(u(i, j)), u(i, j).dx_left().dx_right() + u(i, j).dy_left().dy_right());
ASSERT_DBL_EQ(laplacian(u(i, j)), dx_right(dx_left(u(i, j))) + dy_right(dy_left(u(i, j))));
ASSERT_DBL_EQ(laplacian(u(i, j)), u(i, j).dx_right().dx_left() + u(i, j).dy_right().dy_left());
r(i, j) = laplacian(u(i, j));
ASSERT_DBL_EQ(r(i, j), 4.);
}
}
//итераторы
//все значения сетки
ScalarMeshFunc f(a, b, n, n);
f.all() = 0;
//итератор по внутренним ячейкам в диапазоне индексов
Interval2d I(1, n, 1, n);
f(I) = 4;
r(I) = laplacian(u(I)) - f(I);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
ASSERT_DBL_EQ(r(i, j), 0.);
}
}
}
#include "pdepp.h"
double sln_u(double x, double y)
{
return x * x + y * y;
}
double sln_v(double x, double y)
{
return x * x * x + y * y * y;
}
//векторные сеточные-функции
void test_pdepp_2()
{
//2d-область [a, b] x [a, b]
double a = 0;
double b = 1;
//число внутренних узлов равномерной сетки
int n = 10;
double h = (b - a) / (n + 1);
//инициализация сеточной функции
VectorMeshFunc U(a, b, n, n);
//поиндексный доступ к ячейкам
for (int i = 0; i <= n + 1; i++) {
for (int j = 0; j <= n + 1; j++) {
VectorCell &c = U(i, j);
c.comp(0).val() = sln_u(c.x(), c.y());
c.comp(1).val() = sln_v(c.x(), c.y());
}
}
//сеточный оператор дивергенции
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
ASSERT_DBL_EQ(div_left(U(i, j)), (U(i, j).comp(0) - U(i - 1, j).comp(0)) / h + (U(i, j).comp(1) - U(i, j - 1).comp(1)) / h);
ASSERT_DBL_EQ(div_left(U(i, j)), U(i, j).comp(0).dx_left() + U(i, j).comp(1).dy_left());
ASSERT_DBL_EQ(div_left(U(i, j)), dx_left(U(i, j).comp(0)) + dy_left(U(i, j).comp(1)));
ASSERT_DBL_EQ(div_right(U(i, j)), (U(i + 1, j).comp(0) - U(i, j).comp(0)) / h + (U(i, j + 1).comp(1) - U(i, j).comp(1)) / h);
ASSERT_DBL_EQ(div_right(U(i, j)), U(i, j).comp(0).dx_right() + U(i, j).comp(1).dy_right());
ASSERT_DBL_EQ(div_right(U(i, j)), dx_right(U(i, j).comp(0)) + dy_right(U(i, j).comp(1)));
}
}
//преобразование векторной функции в скалярную через оператор дивергенции
ScalarMeshFunc u(a, b, n, n);
Interval2d I(1, n, 1, n);
u(I) = div_left(U(I));
u(I) = div_right(U(I));
//сеточный оператор градиента
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
VectorCell &c = grad_left(u(i, j));
ASSERT_DBL_EQ(c.comp(0), (u(i, j) - u(i - 1, j)) / h);
ASSERT_DBL_EQ(c.comp(0), u(i, j).dx_left());
ASSERT_DBL_EQ(c.comp(0), dx_left(u(i, j)));
ASSERT_DBL_EQ(c.comp(1), (u(i, j) - u(i, j - 1)) / h);
ASSERT_DBL_EQ(c.comp(1), u(i, j).dy_left());
ASSERT_DBL_EQ(c.comp(1), dy_left(u(i, j)));
c = grad_right(u(i, j));
ASSERT_DBL_EQ(c.comp(0), (u(i + 1, j) - u(i, j)) / h);
ASSERT_DBL_EQ(c.comp(0), u(i, j).dx_right());
ASSERT_DBL_EQ(c.comp(0), dx_right(u(i, j)));
ASSERT_DBL_EQ(c.comp(1), (u(i, j + 1) - u(i, j)) / h);
ASSERT_DBL_EQ(c.comp(1), u(i, j).dy_right());
ASSERT_DBL_EQ(c.comp(1), dy_right(u(i, j)));
}
}
//преобразование скалярной функции в векторную через оператор градиента
U(I) = grad_left(u(I));
U(I) = grad_right(u(I));
ScalarMeshFunc f(a, b, n, n);
f(I) = div_left(grad_right(u(I)));
f(I) = div_right(grad_left(u(I)));
VectorMeshFunc F(a, b, n, n);
F(I) = grad_left(div_right(U(I)));
F(I) = grad_right(div_left(U(I)));
}
Pdepp.h file
//(с) 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