Pustaka C ++ yang ringkas untuk memprogram metode perbedaan terbatas ala operator. Bagian 1. Semantik

Semantik dari pde ++ library yang dikembangkan untuk pemrograman metode hingga-perbedaan dalam gaya operator disajikan. Objek utama perpustakaan adalah fungsi kisi, sel kisi dan operator kisi, hubungan aritmatika di antaranya membawa kode program sedekat mungkin dengan notasi matematikanya. Pde ++ library diwakili oleh hanya beberapa file header, tidak memiliki dependensi eksternal, dan menggunakan konsep komputasi malas.

Sejumlah besar masalah pemodelan matematika direduksi menjadi solusi numerik persamaan diferensial parsial (PDE) dengan metode grid. Dalam teori skema perbedaan (Samarsky A.A.), operator grid yang sesuai membentuk ruang linear di atas fungsi grid yang tidak ada representasi langsung dalam bahasa pemrograman untuk tujuan umum seperti C ++. Akibatnya, praktik merekam hasil penerapan operator kisi ke fungsi kisi menggunakan array atau matriks multidimensi banyak digunakan dalam implementasi perangkat lunak.

Praktek menunjukkan bahwa pendekatan yang disebutkan di atas sangat berguna dalam menguasai keterampilan menerapkan metode numerik, terutama karena visibilitasnya ketika bekerja dengan perkiraan pra-tertulis dari PDE dalam bentuk indeks. Masalah signifikan juga tidak muncul ketika memperluas teknik ini untuk PDE umum, jika dimaksudkan untuk menerapkan skema perbedaan dengan parameter sekali dan menggunakan kembali kode program yang sesuai tanpa modifikasi lebih lanjut.

Dalam kasus umum, program komputasi dapat dimodifikasi ke arah yang berbeda, sehingga teknik yang dijelaskan di atas akan memerlukan penulisan sejumlah besar kode program, dan ini, pada gilirannya, akan meningkatkan kemungkinan kesalahan ketik dan pencatatan yang tidak konsisten dari operator jaringan yang sama dalam modul program yang berbeda. Perlu juga dicatat masalah duplikasi kode program dengan variabilitas dimensi spasial (1D, 2D, 3D) dan metode untuk memperkirakan PDE.

Jadi, alternatifnya adalah pengembangan dan penggunaan perpustakaan perangkat lunak khusus dengan abstraksi domain tingkat tinggi yang membawa kode program lebih dekat ke notasi matematisnya. Dalam pustaka Blitz ++, abstraksi semacam itu adalah perhitungan tensor pada templat perbedaan, diimplementasikan berdasarkan penggunaan teknik metaprogramming templat. Pustaka freePOOMA memperluas konsep Blitz ++ dengan analog perbedaan operator perbedaan divergensi dan gradien dan kemampuan untuk bekerja pada sistem komputasi multiprosesor. Sayangnya, perpustakaan-perpustakaan ini tidak didukung untuk waktu yang lama, memiliki sejumlah keterbatasan (mereka akan dibahas pada bagian selanjutnya) ketika mereka digunakan untuk perkiraan perbedaan-hingga-perbedaan yang cukup klasik dari PDE yang dipertimbangkan dalam tulisan ini.

Pde ++ open-source library, yang dikembangkan oleh penulis, idealnya terinspirasi oleh pustaka freePOOMA dan dirancang untuk merekam dalam bentuk skema beda hingga untuk fungsi skalar dan vektor yang didefinisikan dalam pengaturan 2D (1D dan 3D dalam operasi) pada grid persegi panjang yang seragam.

Peringatan: kode diuji hanya di bawah 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<double> u(a, b, n, n); ScalarMeshFunc<double> r(a, b, n, n); //    for (int i = 0; i <= n + 1; i++) { for (int j = 0; j <= n + 1; j++) { ScalarCell<double> &c = u(i, j); c.val() = sln_u(cx(), cy()); } } //  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<double> 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<double> U(a, b, n, n); //    for (int i = 0; i <= n + 1; i++) { for (int j = 0; j <= n + 1; j++) { VectorCell<double> &c = U(i, j); c.comp(0).val() = sln_u(cx(), cy()); c.comp(1).val() = sln_v(cx(), cy()); } } //   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<double> 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<double> &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<double> f(a, b, n, n); f(I) = div_left(grad_right(u(I))); f(I) = div_right(grad_left(u(I))); VectorMeshFunc<double> F(a, b, n, n); F(I) = grad_left(div_right(U(I))); F(I) = grad_right(div_left(U(I))); } 

File pdepp.h
 //() 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(); } 


Arrays.h
 //() 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_; }; 


NumericAssert.h
 //() 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 

Source: https://habr.com/ru/post/id442316/


All Articles