Skip to content
Snippets Groups Projects
Commit 8f5ef171 authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

Port to fullVector done

parent 3d1046c2
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 134 additions and 341 deletions
FunctionSpace/BasisVector.h
+ 4
− 4
View file @ 8f5ef171
#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;
}
......
FunctionSpace/CMakeLists.txt
+ 0
− 1
View file @ 8f5ef171
......@@ -4,7 +4,6 @@
# bugs and problems to <gmsh@geuz.org>.
set(SRC
VectorPolynomial.cpp
Polynomial.cpp
Legendre.cpp
......
FunctionSpace/Polynomial.cpp
+ 22
− 8
View file @ 8f5ef171
......@@ -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;
......
FunctionSpace/Polynomial.h
+ 17
− 4
View file @ 8f5ef171
......@@ -3,7 +3,8 @@
#include <string>
#include <stack>
#include "Vector.h"
#include <vector>
#include "fullMatrix.h"
/**
@class Polynomial
......@@ -36,7 +37,7 @@ class Polynomial{
~Polynomial(void);
void derivative(const int dim);
Vector<Polynomial> gradient(void) const;
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
......
FunctionSpace/QuadEdgeBasis.cpp
+ 24
− 15
View file @ 8f5ef171
......@@ -74,7 +74,11 @@ 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;
......@@ -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");
}
......
FunctionSpace/TriEdgeBasis.cpp
+ 48
− 20
View file @ 8f5ef171
......@@ -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].add(tmp);
basis[i][0].mul(v[l2]);
basis[i][1].mul(v[l2]);
basis[i][2].mul(v[l2]);
basis[i][0].add(tmp[0]);
basis[i][1].add(tmp[1]);
basis[i][2].add(tmp[2]);
i++;
}
......@@ -102,13 +120,20 @@ 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");
}
......
FunctionSpace/TriNedelecBasis.cpp
+ 19
− 9
View file @ 8f5ef171
......@@ -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].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]);
}
// 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");
}
......
FunctionSpace/Vector.h deleted 100644 → 0
+ 0
− 227
View file @ 3d1046c2
#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
FunctionSpace/VectorPolynomial.cpp deleted 100644 → 0
+ 0
− 53
View file @ 3d1046c2
#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();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment