From 8f5ef171747468038da0239d44b25865fe689f1d Mon Sep 17 00:00:00 2001 From: Nicolas Marsic <nicolas.marsic@gmail.com> Date: Mon, 4 Jun 2012 09:17:42 +0000 Subject: [PATCH] Port to fullVector done --- FunctionSpace/BasisVector.h | 8 +- FunctionSpace/CMakeLists.txt | 1 - FunctionSpace/Polynomial.cpp | 30 +++- FunctionSpace/Polynomial.h | 21 ++- FunctionSpace/QuadEdgeBasis.cpp | 39 +++-- FunctionSpace/TriEdgeBasis.cpp | 68 ++++++--- FunctionSpace/TriNedelecBasis.cpp | 28 ++-- FunctionSpace/Vector.h | 227 ----------------------------- FunctionSpace/VectorPolynomial.cpp | 53 ------- 9 files changed, 134 insertions(+), 341 deletions(-) delete mode 100644 FunctionSpace/Vector.h delete mode 100644 FunctionSpace/VectorPolynomial.cpp diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h index a932f8e938..293adb2a44 100644 --- a/FunctionSpace/BasisVector.h +++ b/FunctionSpace/BasisVector.h @@ -1,9 +1,9 @@ #ifndef _BASISVECTOR_H_ #define _BASISVECTOR_H_ +#include <vector> #include "Basis.h" #include "Polynomial.h" -#include "Vector.h" /** @class BasisVector @@ -16,12 +16,12 @@ class BasisVector: public Basis{ protected: - Vector<Polynomial>* basis; + std::vector<Polynomial>* basis; public: virtual ~BasisVector(void); - Vector<Polynomial>* getBasis(void) const; + std::vector<Polynomial>* getBasis(void) const; protected: BasisVector(void); @@ -41,7 +41,7 @@ class BasisVector: public Basis{ // Inline Functions // ////////////////////// -inline Vector<Polynomial>* BasisVector::getBasis(void) const{ +inline std::vector<Polynomial>* BasisVector::getBasis(void) const{ return basis; } diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt index 7ac05285d6..fa18efc69d 100644 --- a/FunctionSpace/CMakeLists.txt +++ b/FunctionSpace/CMakeLists.txt @@ -4,7 +4,6 @@ # bugs and problems to <gmsh@geuz.org>. set(SRC - VectorPolynomial.cpp Polynomial.cpp Legendre.cpp diff --git a/FunctionSpace/Polynomial.cpp b/FunctionSpace/Polynomial.cpp index e277d748a7..7d2bf67b2c 100644 --- a/FunctionSpace/Polynomial.cpp +++ b/FunctionSpace/Polynomial.cpp @@ -84,18 +84,18 @@ void Polynomial::derivative(const int dim){ return; } -Vector<Polynomial> Polynomial::gradient(void) const{ - Vector<Polynomial> grad(3); +vector<Polynomial> Polynomial::gradient(void) const{ + vector<Polynomial> grad(3); // Copy Polynomial // - grad(0) = *this; - grad(1) = *this; - grad(2) = *this; + grad[0] = *this; + grad[1] = *this; + grad[2] = *this; // Derivative with respect to each direction // - grad(0).derivative(0); - grad(1).derivative(1); - grad(2).derivative(2); + grad[0].derivative(0); + grad[1].derivative(1); + grad[2].derivative(2); return grad; } @@ -113,6 +113,20 @@ double Polynomial::at return val; } +fullVector<double> Polynomial::at(const vector<Polynomial>& P, + const double x, + const double y, + const double z){ + fullVector<double> val(3); + + val(0) = P[0].at(x, y, z); + val(1) = P[1].at(x, y, z); + val(2) = P[2].at(x, y, z); + + return val; +} + + Polynomial Polynomial::operator+(const Polynomial& other) const{ Polynomial newP; diff --git a/FunctionSpace/Polynomial.h b/FunctionSpace/Polynomial.h index 02d4968b16..a7e73a2a47 100644 --- a/FunctionSpace/Polynomial.h +++ b/FunctionSpace/Polynomial.h @@ -3,7 +3,8 @@ #include <string> #include <stack> -#include "Vector.h" +#include <vector> +#include "fullMatrix.h" /** @class Polynomial @@ -35,8 +36,8 @@ class Polynomial{ Polynomial(void); ~Polynomial(void); - void derivative(const int dim); - Vector<Polynomial> gradient(void) const; + void derivative(const int dim); + std::vector<Polynomial> gradient(void) const; double operator() (const double x, const double y, const double z) const; @@ -44,6 +45,10 @@ class Polynomial{ double at (const double x, const double y, const double z) const; + static fullVector<double> at(const std::vector<Polynomial>& P, + const double x, + const double y, + const double z); Polynomial operator+(const Polynomial& other) const; Polynomial operator-(const Polynomial& other) const; @@ -159,13 +164,21 @@ class Polynomial{ @return Returns the @em evaluation of this Polynomial at (@c x, @c y, @c z) - @fn Polynomial::at + @fn double Polynomial::at(const double, const double, const double) const @param x A value @param y A value @param z A value @return Returns the @em evaluation of this Polynomial at (@c x, @c y, @c z) + @fn fullVector<double> Polynomial::at(const std::vector<Polynomial>&, const double, const double y, const double z) + @param P A vector of Polynomial%s + @param x A value + @param y A value + @param z A value + @return Returns a fullVector with the @em evaluation of + the given vector of Polynomial%s at (@c x, @c y, @c z) + @fn Polynomial Polynomial::operator+(const Polynomial&) const @param other An other Polynomial @return Returns a @em new Polynomial, which is the diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp index dd072beb30..9559139687 100644 --- a/FunctionSpace/QuadEdgeBasis.cpp +++ b/FunctionSpace/QuadEdgeBasis.cpp @@ -74,8 +74,12 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ // Basis // - basis = new Vector<Polynomial>[size]; - + basis = new std::vector<Polynomial>[size]; + + for(int i = 0; i < size; i++) + basis[i].resize(3); + + // Edge Based (Nedelec) // int i = 0; Polynomial oneHalf(0.5, 0, 0, 0); @@ -84,8 +88,13 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ basis[i] = (liftingSub[e]).gradient(); - basis[i].mul(lagrangeSum[e]); - basis[i].mul(oneHalf); + basis[i][0].mul(lagrangeSum[e]); + basis[i][1].mul(lagrangeSum[e]); + basis[i][2].mul(lagrangeSum[e]); + + basis[i][0].mul(oneHalf); + basis[i][1].mul(oneHalf); + basis[i][2].mul(oneHalf); i++; } @@ -128,9 +137,9 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ // Cell Based (Type 2) // for(int l1 = 1; l1 < orderPlus; l1++){ for(int l2 = 1; l2 < orderPlus; l2++){ - basis[i](0) = legendreX[l1] * iLegendreY[l2]; - basis[i](1) = iLegendreX[l1] * legendreY[l2] * -1; - basis[i](2) = zero; + basis[i][0] = legendreX[l1] * iLegendreY[l2]; + basis[i][1] = iLegendreX[l1] * legendreY[l2] * -1; + basis[i][2] = zero; i++; } @@ -138,13 +147,13 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ // Cell Based (Type 3) // for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){ - basis[i](0) = iLegendreY[l]; - basis[i](1) = zero; - basis[i](2) = zero; + basis[i][0] = iLegendreY[l]; + basis[i][1] = zero; + basis[i][2] = zero; - basis[iPlus](0) = zero; - basis[iPlus](1) = iLegendreX[l]; - basis[iPlus](2) = zero; + basis[iPlus][0] = zero; + basis[iPlus][1] = iLegendreX[l]; + basis[iPlus][2] = zero; i++; } @@ -178,7 +187,7 @@ int main(void){ QuadEdgeBasis b(P); - const Vector<Polynomial>* basis = b.getBasis(); + const std::vector<Polynomial>* basis = b.getBasis(); printf("\n"); printf("clear all;\n"); @@ -199,7 +208,7 @@ int main(void){ for(int i = 0; i < b.getSize(); i++){ for(int j = 0; j < 2; j++) - printf("p(%d, %d) = %s;\n", i + 1, j + 1, basis[i](j).toString().c_str()); + printf("p(%d, %d) = %s;\n", i + 1, j + 1, basis[i][j].toString().c_str()); //printf("p(%d) = %s", i, basis[i].toString().c_str()); printf("\n"); } diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp index e2004cc53a..c97102565f 100644 --- a/FunctionSpace/TriEdgeBasis.cpp +++ b/FunctionSpace/TriEdgeBasis.cpp @@ -48,19 +48,30 @@ TriEdgeBasis::TriEdgeBasis(const int order){ // Basis // - basis = new Vector<Polynomial>[size]; + basis = new std::vector<Polynomial>[size]; + + for(int i = 0; i < size; i++) + basis[i].resize(3); + // Edge Based (Nedelec) // int i = 0; for(int j = 1; i < 3; j = (j + 1) % 3){ - Vector<Polynomial> tmp = lagrange[j].gradient(); - tmp.mul(lagrange[i]); + std::vector<Polynomial> tmp = lagrange[j].gradient(); + tmp[0].mul(lagrange[i]); + tmp[1].mul(lagrange[i]); + tmp[2].mul(lagrange[i]); basis[i] = lagrange[i].gradient(); - basis[i].mul(lagrange[j]); - - basis[i].sub(tmp); + + basis[i][0].mul(lagrange[j]); + basis[i][1].mul(lagrange[j]); + basis[i][2].mul(lagrange[j]); + + basis[i][0].sub(tmp[0]); + basis[i][1].sub(tmp[1]); + basis[i][2].sub(tmp[2]); i++; } @@ -87,13 +98,20 @@ TriEdgeBasis::TriEdgeBasis(const int order){ // Cell Based (Type 1) // for(int l1 = 1; l1 < order; l1++){ for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){ - Vector<Polynomial> tmp = v[l2].gradient(); - tmp.mul(u[l1]); + std::vector<Polynomial> tmp = v[l2].gradient(); + tmp[0].mul(u[l1]); + tmp[1].mul(u[l1]); + tmp[2].mul(u[l1]); basis[i] = u[l1].gradient(); - basis[i].mul(v[l2]); + + basis[i][0].mul(v[l2]); + basis[i][1].mul(v[l2]); + basis[i][2].mul(v[l2]); - basis[i].add(tmp); + basis[i][0].add(tmp[0]); + basis[i][1].add(tmp[1]); + basis[i][2].add(tmp[2]); i++; } @@ -102,14 +120,21 @@ TriEdgeBasis::TriEdgeBasis(const int order){ // Cell Based (Type 2) // for(int l1 = 1; l1 < order; l1++){ for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){ - Vector<Polynomial> tmp = v[l2].gradient(); - tmp.mul(u[l1]); + std::vector<Polynomial> tmp = v[l2].gradient(); + tmp[0].mul(u[l1]); + tmp[1].mul(u[l1]); + tmp[2].mul(u[l1]); basis[i] = u[l1].gradient(); - basis[i].mul(v[l2]); - basis[i].sub(tmp); - + basis[i][0].mul(v[l2]); + basis[i][1].mul(v[l2]); + basis[i][2].mul(v[l2]); + + basis[i][0].sub(tmp[0]); + basis[i][1].sub(tmp[1]); + basis[i][2].sub(tmp[2]); + i++; } } @@ -117,7 +142,10 @@ TriEdgeBasis::TriEdgeBasis(const int order){ // Cell Based (Type 3) // for(int l = 0; l < orderMinus; l++){ basis[i] = basis[0]; - basis[i].mul(v[l]); + + basis[i][0].mul(v[l]); + basis[i][1].mul(v[l]); + basis[i][2].mul(v[l]); i++; } @@ -139,13 +167,13 @@ TriEdgeBasis::~TriEdgeBasis(void){ /* #include <cstdio> int main(void){ - const int P = 1; + const int P = 6; const double d = 0.05; const char x[2] = {'X', 'Y'}; TriEdgeBasis b(P); - const Vector<Polynomial>* basis = b.getBasis(); + const std::vector<Polynomial>* basis = b.getBasis(); printf("\n"); printf("clear all;\n"); @@ -165,8 +193,8 @@ int main(void){ for(int i = 0; i < b.getSize(); i++){ //printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str()); - printf("p(%d, 1) = %s;\n", i + 1, basis[i](0).toString().c_str()); - printf("p(%d, 2) = %s;\n", i + 1, basis[i](1).toString().c_str()); + printf("p(%d, 1) = %s;\n", i + 1, basis[i][0].toString().c_str()); + printf("p(%d, 2) = %s;\n", i + 1, basis[i][1].toString().c_str()); printf("\n"); } diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp index fd2805125c..d320b29f6a 100644 --- a/FunctionSpace/TriNedelecBasis.cpp +++ b/FunctionSpace/TriNedelecBasis.cpp @@ -21,16 +21,26 @@ TriNedelecBasis::TriNedelecBasis(void){ Polynomial(1, 0, 1, 0); // Basis // - basis = new Vector<Polynomial>[size]; + basis = new std::vector<Polynomial>[size]; + + for(int i = 0; i < size; i++) + basis[i].resize(3); for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){ - Vector<Polynomial> tmp = lagrange[j].gradient(); - tmp.mul(lagrange[i]); + std::vector<Polynomial> tmp = lagrange[j].gradient(); + tmp[0].mul(lagrange[i]); + tmp[1].mul(lagrange[i]); + tmp[2].mul(lagrange[i]); basis[i] = lagrange[i].gradient(); - basis[i].mul(lagrange[j]); + + basis[i][0].mul(lagrange[j]); + basis[i][1].mul(lagrange[j]); + basis[i][2].mul(lagrange[j]); - basis[i].sub(tmp); + basis[i][0].sub(tmp[0]); + basis[i][1].sub(tmp[1]); + basis[i][2].sub(tmp[2]); } // Free Temporary Sapce // @@ -48,9 +58,9 @@ int main(void){ const double d = 0.05; const char x[2] = {'X', 'Y'}; - TriNedelec b; + TriNedelecBasis b; - const Vector<Polynomial>* basis = b.getBasis(); + const std::vector<Polynomial>* basis = b.getBasis(); printf("\n"); printf("clear all;\n"); @@ -70,8 +80,8 @@ int main(void){ for(int i = 0; i < b.getSize(); i++){ //printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str()); - printf("p(%d, 1) = %s;\n", i + 1, basis[i](0).toString().c_str()); - printf("p(%d, 2) = %s;\n", i + 1, basis[i](1).toString().c_str()); + printf("p(%d, 1) = %s;\n", i + 1, basis[i][0].toString().c_str()); + printf("p(%d, 2) = %s;\n", i + 1, basis[i][1].toString().c_str()); printf("\n"); } diff --git a/FunctionSpace/Vector.h b/FunctionSpace/Vector.h deleted file mode 100644 index 361a16c263..0000000000 --- a/FunctionSpace/Vector.h +++ /dev/null @@ -1,227 +0,0 @@ -#ifndef _VECTOR_H_ -#define _VECTOR_H_ - -#include <string> - -extern "C"{ -#include <cblas.h> -} - -#include "fullMatrix.h" -#include "Exception.h" - -/** - @class Vector - @brief Handles vectors - - This class represents a vector of type < T > - and of dimention @c n. - - @todo - Use Inheritance instead of template@n - |@n - |@n - ---> We have methods that changes between types - (e.g: at(...), 3D automatic allocation, scalar mult) - */ - -class Solver; - -template<class T> -class Vector{ - private: - int N; - T* v; - - friend class Solver; - friend class Matrix; - - public: - Vector(const int a); - Vector(void); - ~Vector(void); - - int dim(void) const; - - T& operator()(const int i); - T operator()(const int i) const; - void operator=(const Vector<T>& other); - - T get(const int i) const; - void set(const int i, const T a); - - fullVector<double> at(const double x, const double y, const double z) const; - - Vector<T> operator+(const Vector<T>& other); - Vector<T> operator-(const Vector<T>& other); - Vector<T> operator*(const T& other); - //Vector<T> operator*(const double alpha); - - void add(const Vector<T>& b); - void sub(const Vector<T>& b); - void mul(const T& other); - //void mul(const double alpha); - double dot(const Vector<T>& v) const; - - void allToZero(void); - - std::string toString(void) const; -}; - - -/** - @fn Vector::Vector(int) - @param a Size of the futur Vector - @return Returns a new Vector of size @c a - - @fn Vector::Vector(void) - @return Returns a new Vector of size @em @c 3 - - @fn Vector::~Vector - @return Deletes this Vector - - @fn Vector::dim - @return Returns the @em dimention of this vector - - @fn T& Vector::operator()(const int) - @param i An index of this Vector - @return Returns a @em reference to - the element at position @c i - - @fn T Vector::operator()(const int) const - @param i An index of this Vector - @return Returns the @em value of - the element at position @c i - - @fn Vector::get - @param i An index in this Vector - @return Returns the @em value of the element - at position @c i - - @fn Vector::set - @param i An index in this Vector - @param a A value - @return Sets the given value at position @c i - in this Vector - - @fn Vector::at - @param x The @em first dimesion of a @c 3D Polynomial - @param y The @em second dimesion of a @c 3D Polynomial - @param z The @em third dimesion of a @c 3D Polynomial - @return Retuns the @c 3D Vector resulting of the @em evaluation - of this @em @c 3D Vector @em of @em Polynomial%s - - @warning - This works only for @em @c 3D Vector @em of @em Polynomial%s - - @fn Vector<T> Vector::operator+(const Vector<T>&) - @param other An other Vector - @return Returns the Vector containing - the @em sum of this Vector with the other - - @fn Vector<T> Vector::operator-(const Vector<T>&) - @param other An other Vector - @return Returns the Vector containing - the @em substraction of this Vector with the other - - @fn Vector<T> Vector::operator*(const T&) - @param other A value - @return Returns the Vector containing - the @em product (@em element @em by @em element) - of this Vector with the given value - - @fn Vector::add - @param b An other Vector - @return Adds the given Vector with this one - @note - The result is stored in this Vector - - @fn Vector::sub - @param b An other Vector - @return Substracts the given Vector with this one - @note - The result is stored in this Vector - - @fn Vector::mul - @param other A value - @return Multiplies all the elements of this Vector - with the given value - @note - The result is stored in this Vector - - @fn Vector::dot - @param v An Vector - @return Returns the @em dot @em product of this - Vector with the given one - - @fn Vector::allToZero - @return Sets all the elements of this Vector - to the @em zero @em element - - @fn Vector::toString - @return Returns a string representing this Vector -*/ - -////////////////////////////////////////////////////////////////////// -// Templated Implementations // -////////////////////////////////////////////////////////////////////// - -template<class T> -Vector<T>::Vector(const int a){ - if(!a) - throw Exception("Vector must by of dimension bigger than 0"); - - N = a; - v = new T[N]; -} - -template<class T> -Vector<T>::Vector(void){ - N = 3; - v = new T[N]; -} - -template<class T> -Vector<T>::~Vector(void){ - delete[] v; -} - -////////////////////////////////////////////////////////////////////// -// Inline Templated Implementations // -////////////////////////////////////////////////////////////////////// - -template<class T> -inline int Vector<T>::dim(void) const{ - return N; -} - -template<class T> -inline T& Vector<T>::operator()(const int i){ - return v[i]; -} - -template<class T> -inline T Vector<T>::operator()(const int i) const{ - return v[i]; -} - -template<class T> -void Vector<T>::operator=(const Vector<T>& other){ - if(N != other.N) - throw Exception("Vectors must be of the same dimension"); - - for(int i = 0; i < N; i++) - v[i] = other.v[i]; -} - -template<class T> -inline T Vector<T>::get(const int i) const{ - return v[i]; -} - -template<class T> -inline void Vector<T>::set(const int i, const T a){ - v[i] = a; -} - -#endif diff --git a/FunctionSpace/VectorPolynomial.cpp b/FunctionSpace/VectorPolynomial.cpp deleted file mode 100644 index 9aa3d249c5..0000000000 --- a/FunctionSpace/VectorPolynomial.cpp +++ /dev/null @@ -1,53 +0,0 @@ -#include <sstream> -#include "Vector.h" -#include "Polynomial.h" - -using namespace std; - -template<> -fullVector<double> Vector<Polynomial>::at(const double x, - const double y, - const double z) const{ - fullVector<double> val(N); - - for(int i = 0; i < N; i++) - val(i) = v[i].at(x, y, z); - - return val; -} - -template<> -void Vector<Polynomial>::add(const Vector<Polynomial>& other){ - for(int i = 0; i < N; i++) - v[i].add(other.v[i]); -} - -template<> -void Vector<Polynomial>::sub(const Vector<Polynomial>& other){ - for(int i = 0; i < N; i++) - v[i].sub(other.v[i]); -} - -template<> -void Vector<Polynomial>::mul(const Polynomial& other){ - for(int i = 0; i < N; i++) - v[i].mul(other); -} -/* -template<> -void Vector<Polynomial>::mul(const double alpha){ - for(int i = 0; i < N; i++) - v[i].mul(alpha); -} -*/ -template<> -string Vector<Polynomial>::toString(void) const{ - stringstream s; - - s << endl; - - for(int i = 0; i < N; i++) - s << "[" << v[i].toString() << "]" << endl; - - return s.str(); -} -- GitLab