00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #ifndef __VECT_H
00035 #define __VECT_H
00036
00037 #include <valarray>
00038 using std::valarray;
00039
00040 #include <assert.h>
00041 #include <math.h>
00042 #include <stdlib.h>
00043
00044 #include <iostream>
00045 #include <fstream>
00046 using std::ostream;
00047 using std::istream;
00048 using std::endl;
00049
00050 #include <iomanip>
00051 using std::setw;
00052 using std::ios;
00053 using std::setprecision;
00054
00055 #include "OFELI_Config.h"
00056 #include "util.h"
00057 #include "FDF.h"
00058 #include "Mesh.h"
00059 #include "BCVect.h"
00060 #include "Point.h"
00061 #include "ElementVect.h"
00062 #include "NodeVect.h"
00063 #include "SideVect.h"
00064 #include "EdgeVect.h"
00065
00066 namespace OFELI {
00067
00096 enum {
00097 REGULAR_ACCESS = 10,
00098 FIRST_ACCESS = 1,
00099 ACCESS = 0,
00100 LAST_ACCESS = -1
00101 };
00102
00103
00104 template<class T_> class NodeVect;
00105 template<class T_> class ElementVect;
00106 template<class T_> class SideVect;
00107 template<class T_> class BCVect;
00108 template<class T_> class AbsVect;
00109
00110
00111 template<class T_> class Vect : public valarray<T_>
00112 {
00113
00114 private:
00115 size_t _nx, _ny, _nz;
00116 size_t _imin, _imax;
00117
00118 public:
00119
00120 using valarray<T_>::size;
00121 using valarray<T_>::resize;
00122
00124 Vect();
00125
00128 Vect(size_t n);
00129
00135 Vect(size_t nx, size_t ny);
00136
00143 Vect(size_t nx, size_t ny, size_t nz);
00144
00149 Vect(size_t n, T_ *x);
00150
00153 Vect(const Vect<T_>& v);
00154
00165 Vect(const Element *el, const Vect<T_> &v, int opt=0);
00166
00177 Vect(const Side *sd, const Vect<T_> &v, int opt=0);
00178
00185 Vect(const class Mesh &m, const Vect<T_> &v, const BCVect<T_> &bc);
00186
00192 Vect(const NodeVect<T_> &v);
00193
00201 Vect(const NodeVect<T_> &v, size_t nb_dof, size_t first_dof);
00202
00208 Vect(const ElementVect<T_> &v);
00209
00217 Vect(const ElementVect<T_> &v, size_t nb_dof, size_t first_dof);
00218
00224 Vect(const SideVect<T_> &v);
00225
00233 Vect(const SideVect<T_> &v, size_t nb_dof, size_t first_dof);
00234
00236 ~Vect();
00237
00245 void From(const NodeVect<T_> &v, size_t nb_dof=0, size_t first_dof=1);
00246
00254 void From(const ElementVect<T_> &v, size_t nb_dof=0, size_t first_dof=1);
00255
00263 void From(const SideVect<T_> &v, size_t nb_dof=0, size_t first_dof=1);
00264
00267 void setSize(size_t size);
00268
00273 void setSize(size_t nx, size_t ny);
00274
00280 void setSize(size_t nx, size_t ny, size_t nz);
00281
00283 size_t getSize() const;
00284
00286 double getNorm1() const;
00287
00289 double getNorm2() const;
00290
00292 double getWNorm1() const;
00293
00295 double getWNorm2() const;
00296
00298 double getNormMax() const;
00299
00301 size_t getNx() const;
00302
00304 size_t getNy() const;
00305
00307 size_t getNz() const;
00308
00315 void removeBC(const class Mesh &m, const Vect<T_> &v);
00316
00323 void removeBC(const class Mesh &m, NodeVect<T_> *v);
00324
00329 void transferBC(const class Mesh &m, const Vect<T_> &bc);
00330
00336 void insertBC(const class Mesh &m, const Vect<T_> &v, const Vect<T_> &bc);
00337
00343 void insertBC(const class Mesh &m, const Vect<T_> &v);
00344
00345 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00347 void insertBC1(const Mesh &m, const Vect<T_> &v, const Vect<T_> &bc);
00348 #endif
00349
00354 void Assembly(const Element *el, const Vect<T_> &b);
00355
00360 void Assembly(const Element *el, const T_ *b);
00361
00366 void Assembly(const Side *sd, const Vect<T_> &b);
00367
00372 void Assembly(const Side *sd, T_ *b);
00373
00378 Vect<T_> &MultAdd(const Vect<T_> &x, const T_ &a);
00379
00387 T_ &operator()(size_t i);
00388
00396 T_ operator()(size_t i) const;
00397
00403 T_ &operator()(size_t i, size_t j);
00404
00410 T_ operator() (size_t i, size_t j) const;
00411
00418 T_ &operator()(size_t i, size_t j, size_t k);
00419
00426 T_ operator() (size_t i, size_t j, size_t k) const;
00427
00432 Vect<T_> & operator=(const T_ &a);
00433
00439 Vect<T_> & operator=(const NodeVect<T_> &v);
00440
00446 Vect<T_> & operator=(const ElementVect<T_> &v);
00447
00452 Vect<T_> & operator=(const SideVect<T_> &v);
00453
00458 Vect<T_> & operator+=(const Vect<T_> &x);
00459
00464 Vect<T_> & operator+=(const T_ &a);
00465
00468 Vect<T_> & operator-=(const Vect<T_> &x);
00469
00474 Vect<T_> & operator-=(const T_ &a);
00475
00478 Vect<T_> &operator*=(const T_ &a);
00479
00482 Vect<T_> &operator/=(const T_ &a);
00483
00488 void setPrintView(size_t imin, size_t imax);
00489
00494 void getPrintView(size_t &imin, size_t &imax) const;
00495
00496 friend class FDF;
00497 friend class AbsVect<T_>;
00498 };
00499
00500
00502
00504
00505
00506 template<class T_>
00507 Vect<T_>::Vect() : valarray<T_>(), _nx(0), _ny(0), _nz(0)
00508 {
00509 _imin = 1;
00510 _imax = 10000;
00511 #ifdef _OFELI_DEBUG
00512 std::clog << "An instance of class Vect is constructed.\n";
00513 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00514 #endif
00515 }
00516
00517
00518 template<class T_>
00519 Vect<T_>::Vect(size_t n) : valarray<T_>(n), _nx(0), _ny(0), _nz(0)
00520 {
00521 Clear(*this);
00522 _imin = 1;
00523 _imax = size();
00524 #ifdef _OFELI_DEBUG
00525 std::clog << "An instance of class Vect is constructed.\n";
00526 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00527 #endif
00528 }
00529
00530
00531 template<class T_>
00532 Vect<T_>::Vect(size_t nx, size_t ny)
00533 : valarray<T_>(nx*ny), _nx(nx), _ny(ny), _nz(0)
00534 {
00535 Clear(*this);
00536 _imin = 1;
00537 _imax = size();
00538 #ifdef _OFELI_DEBUG
00539 std::clog << "An instance of class Vect is constructed.\n";
00540 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00541 #endif
00542 }
00543
00544
00545 template<class T_>
00546 Vect<T_>::Vect(size_t nx, size_t ny, size_t nz)
00547 : valarray<T_>(nx*ny*nz), _nx(nx), _ny(ny), _nz(nz)
00548 {
00549 Clear(*this);
00550 _imin = 1;
00551 _imax = size();
00552 #ifdef _OFELI_DEBUG
00553 std::clog << "An instance of class Vect is constructed.\n";
00554 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00555 #endif
00556 }
00557
00558
00559 template<class T_>
00560 Vect<T_>::Vect(size_t n, T_ *x) : valarray<T_>(n), _nx(0), _ny(0), _nz(0)
00561 {
00562 for (size_t i=0; i<n; ++i)
00563 (*this)[i] = x[i];
00564 _imin = 1;
00565 _imax = size();
00566 #ifdef _OFELI_DEBUG
00567 std::clog << "An instance of class Vect is constructed.\n";
00568 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00569 #endif
00570 }
00571
00572
00573 template<class T_>
00574 Vect<T_>::Vect(const Vect<T_>& v) : valarray<T_>(v), _nx(v._nx), _ny(v._ny), _nz(v._nz)
00575 {
00576 _imin = v._imin;
00577 _imax = v._imax;
00578 }
00579
00580
00581 template<class T_>
00582 Vect<T_>::Vect(const Element *el, const Vect<T_> &v, int opt)
00583 : valarray<T_>(), _nx(0), _ny(0), _nz(0)
00584 {
00585 size_t i = 0;
00586 if (opt==0) {
00587 resize(el->getNbEq());
00588 for (size_t n=1; n<=el->getNbNodes(); ++n) {
00589 Node *nd = el->getPtrNode(n);
00590 size_t k = nd->getFirstDOF();
00591 for (size_t j=1; j<=nd->getNbDOF(); ++j)
00592 (*this)[i++] = v(k++);
00593 }
00594 }
00595 else {
00596 resize(el->getNbNodes());
00597 for (size_t n=1; n<=el->getNbNodes(); ++n) {
00598 Node *nd = el->getPtrNode(n);
00599 size_t k = nd->getFirstDOF();
00600 for (size_t j=1; j<=nd->getNbDOF(); ++j)
00601 (*this)[i++] = v(k++);
00602 }
00603 }
00604 _imin = 1;
00605 _imax = size();
00606 #ifdef _OFELI_DEBUG
00607 std::clog << "An instance of class Vect is constructed.\n";
00608 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00609 #endif
00610 }
00611
00612
00613 template<class T_>
00614 Vect<T_>::Vect(const Side *sd, const Vect<T_> &v, int opt)
00615 : valarray<T_>(), _nx(0), _ny(0), _nz(0)
00616 {
00617 size_t i = 0;
00618 if (opt==0) {
00619 resize(sd->getNbNodes()*sd->getNbDOF());
00620 for (size_t n=1; n<=sd->getNbNodes(); ++n) {
00621 Node *nd = sd->getPtrNode(n);
00622 size_t k = nd->getFirstDOF();
00623 for (size_t j=1; j<=nd->getNbDOF(); ++j)
00624 (*this)[i++] = v(k++);
00625 }
00626 }
00627 else {
00628 resize(sd->getNbNodes());
00629 for (size_t n=1; n<=sd->getNbNodes(); ++n) {
00630 Node *nd = sd->getPtrNode(n);
00631 size_t k = nd->getFirstDOF();
00632 for (size_t j=1; j<=nd->getNbDOF(); ++j)
00633 (*this)[i++] = v(k++);
00634 }
00635 }
00636 _imin = 1;
00637 _imax = size();
00638 #ifdef _OFELI_DEBUG
00639 std::clog << "An instance of class Vect is constructed.\n";
00640 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00641 #endif
00642 }
00643
00644
00645 template<class T_>
00646 Vect<T_>::Vect(const class Mesh &m, const Vect<T_> &v, const BCVect<T_> &bc)
00647 : valarray<T_>(bc.size()), _nx(0), _ny(0), _nz(0)
00648 {
00649 const Node *nd;
00650 #ifdef _BOUNDS_DEBUG
00651 assert(v.size() == size());
00652 #endif
00653 size_t i=0, n=0;
00654 for (m.topNode(); (nd=m.getNode());) {
00655 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
00656 (*this)[i] = bc[n++];
00657 if (nd->getCode(k) == 0)
00658 (*this)[i] = v[nd->getDOF(k)-1];
00659 i++;
00660 }
00661 }
00662 _imin = 1;
00663 _imax = size();
00664 #ifdef _OFELI_DEBUG
00665 std::clog << "An instance of class Vect is constructed.\n";
00666 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00667 #endif
00668 }
00669
00670
00671 template<class T_>
00672 Vect<T_>::Vect(const NodeVect<T_> &v)
00673 : valarray<T_>(v.getLength()), _nx(0), _ny(0), _nz(0)
00674 {
00675 From(v);
00676 _imin = 1;
00677 _imax = size();
00678 #ifdef _OFELI_DEBUG
00679 std::clog << "An instance of class Vect is constructed.\n";
00680 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00681 #endif
00682 }
00683
00684
00685 template<class T_>
00686 Vect<T_>::Vect(const NodeVect<T_> &v, size_t nb_dof, size_t first_dof)
00687 : valarray<T_>(v.getNbNodes()*nb_dof), _nx(0), _ny(0), _nz(0)
00688 {
00689 From(v,nb_dof,first_dof);
00690 _imin = 1;
00691 _imax = size();
00692 #ifdef _OFELI_DEBUG
00693 std::clog << "An instance of class Vect is constructed.\n";
00694 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00695 #endif
00696 }
00697
00698
00699 template<class T_>
00700 Vect<T_>::Vect(const ElementVect<T_> &v)
00701 : valarray<T_>(v.getLength()), _nx(0), _ny(0), _nz(0)
00702 {
00703 From(v);
00704 _imin = 1;
00705 _imax = size();
00706 #ifdef _OFELI_DEBUG
00707 std::clog << "An instance of class Vect is constructed.\n";
00708 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00709 #endif
00710 }
00711
00712
00713 template<class T_>
00714 Vect<T_>::Vect(const ElementVect<T_> &v, size_t nb_dof, size_t first_dof)
00715 : valarray<T_>(v.getNbElements()*nb_dof), _nx(0), _ny(0), _nz(0)
00716 {
00717 From(v,nb_dof,first_dof);
00718 _imin = 1;
00719 _imax = size();
00720 #ifdef _OFELI_DEBUG
00721 std::clog << "An instance of class Vect is constructed.\n";
00722 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00723 #endif
00724 }
00725
00726
00727 template<class T_>
00728 Vect<T_>::Vect(const SideVect<T_> &v)
00729 : valarray<T_>(v.getLength()), _nx(0), _ny(0), _nz(0)
00730 {
00731 From(v);
00732 _imin = 1;
00733 _imax = size();
00734 #ifdef _OFELI_DEBUG
00735 std::clog << "An instance of class Vect is constructed.\n";
00736 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00737 #endif
00738 }
00739
00740
00741 template<class T_>
00742 Vect<T_>::Vect(const SideVect<T_> &v, size_t nb_dof, size_t first_dof)
00743 : valarray<T_>(v.getNbSides()*nb_dof), _nx(0), _ny(0), _nz(0)
00744 {
00745 From(v,nb_dof,first_dof);
00746 _imin = 1;
00747 _imax = size();
00748 #ifdef _OFELI_DEBUG
00749 std::clog << "An instance of class Vect is constructed.\n";
00750 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00751 #endif
00752 }
00753
00754
00755 template<class T_>
00756 Vect<T_>::~Vect()
00757 {
00758 #ifdef _OFELI_DEBUG
00759 std::clog << "An instance of class Vect is destructed.\n";
00760 std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << std::endl;
00761 #endif
00762 }
00763
00764
00765 template<class T_>
00766 void Vect<T_>::From(const NodeVect<T_> &v, size_t nb_dof, size_t first_dof)
00767 {
00768 if (nb_dof==0)
00769 nb_dof = v.getNbDOF();
00770 resize(nb_dof*v.getNbNodes());
00771 size_t i=0;
00772 for (size_t n=1; n<=v.getNbNodes(); ++n)
00773 for (size_t j=first_dof; j<=nb_dof+first_dof-1; j++)
00774 (*this)[i++] = v(n,j);
00775 }
00776
00777
00778 template<class T_>
00779 void Vect<T_>::From(const ElementVect<T_> &v, size_t nb_dof, size_t first_dof)
00780 {
00781 resize(nb_dof*v.getNbElements());
00782 if (nb_dof==0)
00783 nb_dof = v.getNbDOF();
00784 size_t i = 0;
00785 for (size_t n=1; n<=v.getNbElements(); ++n)
00786 for (size_t j=first_dof; j<=nb_dof+first_dof-1; ++j)
00787 (*this)[i++] = v(n,j);
00788 }
00789
00790
00791 template<class T_>
00792 void Vect<T_>::From(const SideVect<T_> &v, size_t nb_dof, size_t first_dof)
00793 {
00794 resize(nb_dof*v.getNbSides());
00795 if (nb_dof==0)
00796 nb_dof = v.getNbDOF();
00797 size_t i = 0;
00798 for (size_t n=1; n<=v.getNbSides(); ++n)
00799 for (size_t j=first_dof; j<=nb_dof+first_dof-1; j++)
00800 (*this)[i++] = v(n,j);
00801 }
00802
00803
00804 template<class T_>
00805 void Vect<T_>::setSize(size_t size)
00806 {
00807 resize(size);
00808 for (size_t i=0; i<size; ++i)
00809 (*this)[i] = 0;
00810 _imin = 1;
00811 _imax = size;
00812 }
00813
00814
00815 template<class T_>
00816 void Vect<T_>::setSize(size_t nx, size_t ny)
00817 {
00818 _nx = nx;
00819 _ny = ny;
00820 resize(_nx*_ny);
00821 Clear(*this);
00822 }
00823
00824
00825 template<class T_>
00826 void Vect<T_>::setSize(size_t nx, size_t ny, size_t nz)
00827 {
00828 _nx = nx;
00829 _ny = ny;
00830 _nz = nz;
00831 resize(_nx*_ny*_nz);
00832 Clear(*this);
00833 }
00834
00835
00836 template<class T_>
00837 size_t Vect<T_>::getSize() const
00838 {
00839 return size();
00840 }
00841
00842
00843 template<class T_>
00844 double Vect<T_>::getNorm1() const
00845 {
00846 double s = 0.;
00847 for (size_t i=0; i<size(); ++i)
00848 s += Abs((*this)[i]);
00849 return s;
00850 }
00851
00852
00853 template<class T_>
00854 double Vect<T_>::getNorm2() const
00855 {
00856 double s = 0;
00857 for (size_t i=0; i<size(); ++i)
00858 s += Abs2((*this)[i]);
00859 return sqrt(s);
00860 }
00861
00862
00863 template<class T_>
00864 double Vect<T_>::getWNorm1() const
00865 {
00866 return getNorm1()/size();
00867 }
00868
00869
00870 template<class T_>
00871 double Vect<T_>::getWNorm2() const
00872 {
00873 return getNorm2()/sqrt(double(size()));
00874 }
00875
00876
00877 template<class T_>
00878 double Vect<T_>::getNormMax() const
00879 {
00880 double s = 0.;
00881 for (size_t i=0; i<size(); ++i)
00882 s = (Abs((*this)[i]) > s ? Abs((*this)[i]) : s);
00883 return s;
00884 }
00885
00886
00887 template<class T_>
00888 size_t Vect<T_>::getNx() const
00889 {
00890 return _nx;
00891 }
00892
00893
00894 template<class T_>
00895 size_t Vect<T_>::getNy() const
00896 {
00897 return _ny;
00898 }
00899
00900
00901 template<class T_>
00902 size_t Vect<T_>::getNz() const
00903 {
00904 return _nz;
00905 }
00906
00907
00908 template<class T_>
00909 void Vect<T_>::removeBC(const class Mesh &m, const Vect<T_> &v)
00910 {
00911 const Node *nd;
00912 size_t n = 0;
00913 for (m.topNode(); (nd=m.getNode());) {
00914 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
00915 if (nd->getCode(k) == 0)
00916 (*this)[nd->getDOF(k)-1] = v[n];
00917 n++;
00918 }
00919 }
00920 #ifdef _OFELI_DEBUG
00921 std::clog << "An instance of class Vect is constructed.\n";
00922 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00923 #endif
00924 }
00925
00926
00927 template<class T_>
00928 void Vect<T_>::removeBC(const class Mesh &m, NodeVect<T_> *v)
00929 {
00930 resize(m.getNbDOF());
00931 const Node *nd;
00932 size_t n = 0;
00933 for (m.topNode(); (nd=m.getNode());) {
00934 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
00935 if (nd->getCode(k) == 0)
00936 (*this)[n++] = (*v)(nd->getLabel(),k);
00937 }
00938 }
00939 #ifdef _OFELI_DEBUG
00940 std::clog << "An instance of class Vect is constructed.\n";
00941 std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << std::endl;
00942 #endif
00943 }
00944
00945
00946 template<class T_>
00947 void Vect<T_>::transferBC(const class Mesh &m, const Vect<T_> &bc)
00948 {
00949 if (m.NodesAreDOF()) {
00950 const Node *nd;
00951 size_t i=0, n=0;
00952 for (m.topNode(); (nd=m.getNode());)
00953 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
00954 if (nd->getCode(k) > 0)
00955 (*this)[i] = bc[n];
00956 i++; n++;
00957 }
00958 }
00959 else if (m.SidesAreDOF()) {
00960 const Side *sd;
00961 size_t i=0, n=0;
00962 for (m.topSide(); (sd=m.getSide());)
00963 for (size_t k=1; k<=sd->getNbDOF(); ++k) {
00964 if (sd->getCode(k) > 0)
00965 (*this)[i] = bc[n];
00966 i++; n++;
00967 }
00968 }
00969 else
00970 ;
00971 }
00972
00973
00974 template<class T_>
00975 void Vect<T_>::insertBC(const class Mesh &m, const Vect<T_> &v, const Vect<T_> &bc)
00976 {
00977 if (m.NodesAreDOF()) {
00978 const Node *nd;
00979 size_t i=0, n=0;
00980 for (m.topNode(); (nd=m.getNode());)
00981 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
00982 if (nd->getCode(k) == 0)
00983 (*this)[i] = v[nd->getDOF(k)-1];
00984 else
00985 (*this)[i] = bc[n];
00986 i++;
00987 n++;
00988 }
00989 }
00990 else if (m.SidesAreDOF()) {
00991 const Side *sd;
00992 size_t i=0, n=0;
00993 for (m.topSide(); (sd=m.getSide());)
00994 for (size_t k=1; k<=sd->getNbDOF(); ++k) {
00995 if (sd->getCode(k) >= 0)
00996 (*this)[i] = v[sd->getDOF(k)-1];
00997 else
00998 (*this)[i] = bc[n];
00999 i++;
01000 n++;
01001 }
01002 }
01003 else
01004 ;
01005 }
01006
01007
01008 template<class T_>
01009 void Vect<T_>::insertBC(const class Mesh &m, const Vect<T_> &v)
01010 {
01011 if (m.NodesAreDOF()) {
01012 const Node *nd;
01013 size_t i=0;
01014 for (m.topNode(); (nd=m.getNode());) {
01015 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
01016 (*this)[i] = 0;
01017 if (nd->getCode(k) == 0)
01018 (*this)[i] = v[nd->getDOF(k)-1];
01019 i++;
01020 }
01021 }
01022 }
01023 else if (m.SidesAreDOF()) {
01024 const Side *sd;
01025 size_t i=0;
01026 for (m.topSide(); (sd=m.getSide());) {
01027 for (size_t k=1; k<=sd->getNbDOF(); ++k) {
01028 (*this)[i] = 0;
01029 if (sd->getCode(k) == 0)
01030 (*this)[i] = v[sd->getDOF(k)-1];
01031 i++;
01032 }
01033 }
01034 }
01035 else
01036 ;
01037 }
01038
01039
01040 #ifndef DOXYGEN_SHOULD_SKIP_THIS
01041 template<class T_>
01042 void Vect<T_>::insertBC1(const Mesh &m, const Vect<T_> &v, const Vect<T_> &bc)
01043 {
01044 const Node *nd;
01045 size_t i=0, n=0, l=0;
01046 for (m.topNode(); (nd=m.getNode());)
01047 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
01048 (*this)[i] = bc[n++];
01049 if (nd->getCode(k) == 0)
01050 (*this)[i] = v[l++];
01051 i++;
01052 }
01053 }
01054 #endif
01055
01056
01057 template<class T_>
01058 void Vect<T_>::Assembly(const Element *el, const Vect<T_> &b)
01059 {
01060 size_t i = 0;
01061 for (size_t n=1; n<=el->getNbNodes(); ++n) {
01062 Node *nd = el->getPtrNode(n);
01063 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
01064 if (nd->getDOF(k))
01065 (*this)[nd->getDOF(k)-1] += b[i];
01066 i++;
01067 }
01068 }
01069 }
01070
01071
01072 template<class T_>
01073 void Vect<T_>::Assembly(const Element *el, const T_ *b)
01074 {
01075 size_t i=0;
01076 for (size_t n=1; n<=el->getNbNodes(); ++n) {
01077 Node *nd = el->getPtrNode(n);
01078 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
01079 if (nd->getDOF(k))
01080 (*this)(nd->getDOF(k)) += b[i];
01081 i++;
01082 }
01083 }
01084 }
01085
01086
01087 template<class T_>
01088 void Vect<T_>::Assembly(const Side *sd, const Vect<T_> &b)
01089 {
01090 size_t nb_eq = 0;
01091 for (size_t i1=1; i1<=sd->getNbNodes(); ++i1)
01092 for (size_t k1=1; k1<=sd->getPtrNode(i1)->getNbDOF(); ++k1)
01093 nb_eq++;
01094 size_t i = 0;
01095 for (size_t n=1; n<=sd->getNbNodes(); ++n) {
01096 Node *nd = sd->getPtrNode(n);
01097 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
01098 if (nd->getDOF(k))
01099 (*this)[nd->getDOF(k)-1] += b[i];
01100 i++;
01101 }
01102 }
01103 }
01104
01105
01106 template<class T_>
01107 void Vect<T_>::Assembly(const Side *sd, T_ *b)
01108 {
01109 size_t nb_eq = 0;
01110 for (size_t i1=1; i1<=sd->getNbNodes(); ++i1)
01111 for (size_t k1=1; k1<=sd->getPtrNode(i1)->getNbDOF(); ++k1)
01112 nb_eq++;
01113 size_t i = 0;
01114 for (size_t n=1; n<=sd->getNbNodes(); ++n) {
01115 Node *nd = sd->getPtrNode(n);
01116 for (size_t k=1; k<=nd->getNbDOF(); ++k) {
01117 if (nd->getDOF(k))
01118 (*this)[nd->getDOF(k)-1] += b[i];
01119 i++;
01120 }
01121 }
01122 }
01123
01124
01125 template<class T_>
01126 Vect<T_> &Vect<T_>::MultAdd(const Vect<T_> &x, const T_ &a)
01127 {
01128 #ifdef _BOUNDS_DEBUG
01129 assert(x.size() == size());
01130 #endif
01131 for (size_t i=0; i<size(); ++i)
01132 (*this)[i] += a * x[i];
01133 return *this;
01134 }
01135
01136
01137 template<class T_>
01138 T_ &Vect<T_>::operator()(size_t i)
01139 {
01140 #ifdef _BOUNDS_DEBUG
01141 assert(i>0 && i<=size());
01142 #endif
01143 return (*this)[i-1];
01144 }
01145
01146
01147 template<class T_>
01148 T_ Vect<T_>::operator()(size_t i) const
01149 {
01150 return (*this)[i-1];
01151 }
01152
01153
01154 template<class T_>
01155 T_ &Vect<T_>::operator()(size_t i, size_t j)
01156 {
01157 #ifdef _BOUNDS_DEBUG
01158 assert(i>0 && i<=_nx);
01159 assert(j>0 && j<=_ny);
01160 #endif
01161 return (*this)[_nx*(i-1)+j-1];
01162 }
01163
01164
01165 template<class T_>
01166 T_ Vect<T_>::operator()(size_t i, size_t j) const
01167 {
01168 return (*this)[_nx*(i-1)+j-1];
01169 }
01170
01171
01172 template<class T_>
01173 T_ &Vect<T_>::operator()(size_t i, size_t j, size_t k)
01174 {
01175 #ifdef _BOUNDS_DEBUG
01176 assert(i>0 && i<=_nx);
01177 assert(j>0 && j<=_ny);
01178 assert(k>0 && k<=_nz);
01179 #endif
01180 return (*this)[_nx*(i-1)+_ny*(j-1)+k];
01181 }
01182
01183
01184 template<class T_>
01185 T_ Vect<T_>::operator() (size_t i, size_t j, size_t k) const
01186 {
01187 return (*this)[_nx*(i-1)+_ny*(j-1)+k];
01188 }
01189
01190
01191 template<class T_>
01192 Vect<T_> &Vect<T_>::operator=(const T_ &a)
01193 {
01194 for (size_t i=0; i<size(); ++i)
01195 (*this)[i] = a;
01196 return *this;
01197 }
01198
01199
01200 template<class T_>
01201 Vect<T_> &Vect<T_>::operator=(const NodeVect<T_> &v)
01202 {
01203 size_t k = 0;
01204 for (size_t n=1; n<=v.getNbNodes(); ++n)
01205 for (size_t i=1; i<=v.getNbDOF(); ++i)
01206 (*this)[k++] = v(n,i);
01207 return *this;
01208 }
01209
01210
01211 template<class T_>
01212 Vect<T_> &Vect<T_>::operator=(const ElementVect<T_> &v)
01213 {
01214 size_t k = 0;
01215 for (size_t n=1; n<=v.getNbElements(); ++n)
01216 for (size_t i=1; i<=v.getNbDOF(); ++i)
01217 (*this)[k++] = v(n,i);
01218 return *this;
01219 }
01220
01221
01222 template<class T_>
01223 Vect<T_> &Vect<T_>::operator=(const SideVect<T_> &v)
01224 {
01225 size_t k = 0;
01226 for (size_t n=1; n<=v.getNbSides(); ++n)
01227 for (size_t i=1; i<=v.getNbDOF(); ++i)
01228 (*this)[k++] = v(n,i);
01229 return *this;
01230 }
01231
01232
01233 template<class T_>
01234 Vect<T_> &Vect<T_>::operator+=(const Vect<T_> &x)
01235 {
01236 #ifdef _BOUNDS_DEBUG
01237 assert(x.size() == size());
01238 #endif
01239 for (size_t i=0; i<size(); ++i)
01240 (*this)[i] += x[i];
01241 return *this;
01242 }
01243
01244
01245 template<class T_>
01246 Vect<T_> &Vect<T_>::operator+=(const T_ &a)
01247 {
01248 #ifdef _BOUNDS_DEBUG
01249 assert(x.size() == size());
01250 #endif
01251 for (size_t i=0; i<size(); ++i)
01252 (*this)[i] += a;
01253 return *this;
01254 }
01255
01256
01257 template<class T_>
01258 Vect<T_> &Vect<T_>::operator-=(const Vect<T_> &x)
01259 {
01260 #ifdef _BOUNDS_DEBUG
01261 assert(x.size() == size());
01262 #endif
01263 for (size_t i=0; i<size(); ++i)
01264 (*this)[i] -= x[i];
01265 return *this;
01266 }
01267
01268
01269 template<class T_>
01270 Vect<T_> &Vect<T_>::operator-=(const T_ &a)
01271 {
01272 #ifdef _BOUNDS_DEBUG
01273 assert(x.size() == size());
01274 #endif
01275 for (size_t i=0; i<size(); ++i)
01276 (*this)[i] -= a;
01277 return *