Une bibliothèque C ++ compacte pour programmer des méthodes de différences finies de type opérateur. Partie 1. Sémantique

La sémantique de la bibliothèque pde ++ développée pour la programmation de méthodes de différences finies dans le style opérateur est présentée. Les principaux objets de la bibliothèque sont une fonction de grille, une cellule de grille et des opérateurs de grille, les relations arithmétiques entre lesquelles rapprochent le code du programme au plus près de sa notation mathématique. La bibliothèque pde ++ n'est représentée que par quelques fichiers d'en-tête, n'a pas de dépendances externes et utilise le concept de l'informatique paresseuse.

Un grand nombre de problèmes de modélisation mathématique est réduit à la solution numérique des équations aux dérivées partielles (EDP) par des méthodes de grille. Dans la théorie des schémas de différence (Samarsky A.A.), les opérateurs de grille correspondants forment un espace linéaire sur les fonctions de grille pour lesquelles il n'y a pas de représentation directe dans les langages de programmation à usage général tels que C ++. Par conséquent, la pratique consistant à enregistrer le résultat de l'application d'opérateurs de grille à des fonctions de grille à l'aide de matrices ou matrices multidimensionnelles est largement utilisée dans l'implémentation logicielle.

La pratique montre que l'approche notée ci-dessus est très utile pour maîtriser les compétences de mise en œuvre de méthodes numériques, principalement en raison de sa visibilité lors de l'utilisation d'approximations pré-écrites de PDE sous forme d'index. Des problèmes importants ne se posent pas non plus lors de l'extension de cette technique aux PDE généralisés, si elle est destinée à implémenter une fois un schéma de différence avec des paramètres et à réutiliser le code de programme approprié sans autres modifications.

Dans le cas général, un programme informatique peut être modifié dans différentes directions, de sorte que la technique décrite ci-dessus nécessitera l'écriture d'une quantité importante de code de programme, ce qui, à son tour, augmentera la probabilité de fautes de frappe et d'enregistrement incohérent des mêmes opérateurs de réseau dans différents modules de programme. Il convient également de noter le problème de la duplication du code du programme avec la variabilité des dimensions spatiales (1D, 2D, 3D) et les méthodes d'approximation des PDE.

Ainsi, une alternative est le développement et l'utilisation de bibliothèques de logiciels spécialisés avec des abstractions de domaine de haut niveau qui rapprochent le code du programme de sa notation mathématique. Dans la bibliothèque Blitz ++, une telle abstraction est des calculs tensoriels sur des modèles de différence, implémentés sur la base de l'utilisation de la technique de métaprogrammation de modèles. La bibliothèque freePOOMA étend le concept Blitz ++ avec des analogues de différence d'opérateurs de divergence différentielle et de gradient et la possibilité de travailler sur des systèmes informatiques multiprocesseurs. Malheureusement, ces bibliothèques n'ont pas été prises en charge depuis longtemps, possédant un certain nombre de limitations (elles seront discutées dans la partie suivante) lorsqu'elles sont utilisées pour des approximations de différences finies assez classiques des PDE considérés dans cet article.

La bibliothèque open-source pde ++, développée par l'auteur, s'inspire idéalement de la bibliothèque freePOOMA et est conçue pour enregistrer sous la forme de schémas de différences finies pour les fonctions de grille scalaire et vectorielle définies dans un cadre 2D (1D et 3D en fonctionnement) sur des grilles rectangulaires uniformes.

Attention: le code n'a été testé que sous 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))); } 

Fichier 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/fr442316/


All Articles