From 34209e7d1f4a7e2cfd475cad96fd58b71af699e6 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 15 Jan 2013 09:48:41 +0000
Subject: [PATCH] Reorganisation of FunctionSpace -- Basis can be ONLY
evaluated + Derivatives -- New Dof ordering: Iterate on Entity and then on
Function Order
---
FunctionSpace/Basis.h | 18 +-
FunctionSpace/BasisGenerator.cpp | 10 +-
FunctionSpace/BasisGenerator.h | 3 +-
FunctionSpace/BasisHierarchicalScalar.cpp | 217 +++++++++++++++
FunctionSpace/BasisHierarchicalScalar.h | 56 ++++
FunctionSpace/BasisHierarchicalVector.cpp | 325 ++++++++++++++++++++++
FunctionSpace/BasisHierarchicalVector.h | 68 +++++
FunctionSpace/BasisLagrange.cpp | 50 ++++
FunctionSpace/BasisLagrange.h | 33 +++
FunctionSpace/BasisScalar.cpp | 31 ---
FunctionSpace/BasisScalar.h | 62 +----
FunctionSpace/BasisVector.cpp | 39 ---
FunctionSpace/BasisVector.h | 68 ++---
FunctionSpace/CMakeLists.txt | 9 +-
FunctionSpace/CurlBasis.cpp | 45 ---
FunctionSpace/CurlBasis.h | 26 --
FunctionSpace/DivBasis.cpp | 45 ---
FunctionSpace/DivBasis.h | 27 --
FunctionSpace/EvaluatedBasis.cpp | 75 -----
FunctionSpace/EvaluatedBasis.h | 106 -------
FunctionSpace/FunctionSpace.cpp | 26 +-
FunctionSpace/FunctionSpace.h | 44 ---
FunctionSpace/FunctionSpaceEdge.cpp | 47 ++--
FunctionSpace/FunctionSpaceNode.cpp | 30 +-
FunctionSpace/FunctionSpaceScalar.cpp | 51 ----
FunctionSpace/FunctionSpaceScalar.h | 63 ++---
FunctionSpace/FunctionSpaceVector.cpp | 81 ------
FunctionSpace/FunctionSpaceVector.h | 102 +++----
FunctionSpace/GradBasis.cpp | 45 ---
FunctionSpace/GradBasis.h | 27 --
FunctionSpace/HexEdgeBasis.h | 4 +-
FunctionSpace/HexNodeBasis.h | 4 +-
FunctionSpace/LagrangeBasis.cpp | 99 -------
FunctionSpace/LagrangeBasis.h | 97 -------
FunctionSpace/LineEdgeBasis.cpp | 16 +-
FunctionSpace/LineEdgeBasis.h | 4 +-
FunctionSpace/LineNedelecBasis.cpp | 14 +-
FunctionSpace/LineNedelecBasis.h | 4 +-
FunctionSpace/LineNodeBasis.cpp | 16 +-
FunctionSpace/LineNodeBasis.h | 4 +-
FunctionSpace/QuadEdgeBasis.h | 4 +-
FunctionSpace/QuadNodeBasis.h | 4 +-
FunctionSpace/TetEdgeBasis.cpp | 312 ++++++++-------------
FunctionSpace/TetEdgeBasis.h | 4 +-
FunctionSpace/TetNodeBasis.cpp | 34 +--
FunctionSpace/TetNodeBasis.h | 4 +-
FunctionSpace/TriEdgeBasis.cpp | 85 +++---
FunctionSpace/TriEdgeBasis.h | 4 +-
FunctionSpace/TriLagrangeBasis.cpp | 178 ++----------
FunctionSpace/TriLagrangeBasis.h | 14 +-
FunctionSpace/TriNedelecBasis.cpp | 12 +-
FunctionSpace/TriNedelecBasis.h | 4 +-
FunctionSpace/TriNodeBasis.cpp | 24 +-
FunctionSpace/TriNodeBasis.h | 4 +-
54 files changed, 1173 insertions(+), 1605 deletions(-)
create mode 100644 FunctionSpace/BasisHierarchicalScalar.cpp
create mode 100644 FunctionSpace/BasisHierarchicalScalar.h
create mode 100644 FunctionSpace/BasisHierarchicalVector.cpp
create mode 100644 FunctionSpace/BasisHierarchicalVector.h
create mode 100644 FunctionSpace/BasisLagrange.cpp
create mode 100644 FunctionSpace/BasisLagrange.h
delete mode 100644 FunctionSpace/CurlBasis.cpp
delete mode 100644 FunctionSpace/CurlBasis.h
delete mode 100644 FunctionSpace/DivBasis.cpp
delete mode 100644 FunctionSpace/DivBasis.h
delete mode 100644 FunctionSpace/EvaluatedBasis.cpp
delete mode 100644 FunctionSpace/EvaluatedBasis.h
delete mode 100644 FunctionSpace/GradBasis.cpp
delete mode 100644 FunctionSpace/GradBasis.h
delete mode 100644 FunctionSpace/LagrangeBasis.cpp
delete mode 100644 FunctionSpace/LagrangeBasis.h
diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 6a3badc137..27db1ac355 100644
--- a/FunctionSpace/Basis.h
+++ b/FunctionSpace/Basis.h
@@ -1,9 +1,6 @@
#ifndef _BASIS_H_
#define _BASIS_H_
-#include <string>
-#include "ReferenceSpace.h"
-
/**
@interface Basis
@brief Common Interface of all Basis
@@ -19,9 +16,7 @@
class Basis{
protected:
- const ReferenceSpace* refSpace;
- unsigned int nRefSpace;
- bool scalar;
+ bool scalar;
unsigned int order;
unsigned int type;
@@ -39,9 +34,6 @@ class Basis{
//!
virtual ~Basis(void);
- //! @return Returns the ReferenceSpace of this Basis
- const ReferenceSpace& getReferenceSpace(void) const;
-
//! @return Returns:
//! @li @c true, if this is a
//! @em scalar Basis (BasisScalar)
@@ -90,9 +82,6 @@ class Basis{
//! in this Basis
unsigned int getNFunction(void) const;
- //! @return Returns the Basis String
- virtual std::string toString(void) const = 0;
-
protected:
//! @internal
//! Instantiate a new Basis
@@ -105,11 +94,6 @@ class Basis{
// Inline Functions //
//////////////////////
-inline const ReferenceSpace&
-Basis::getReferenceSpace(void) const{
- return *refSpace;
-}
-
inline bool Basis::isScalar(void) const{
return scalar;
}
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index ee7f4f46b6..94acc23326 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -12,7 +12,7 @@
#include "TriNodeBasis.h"
#include "TriEdgeBasis.h"
#include "TriNedelecBasis.h"
-#include "TriLagrangeBasis.h"
+//#include "TriLagrangeBasis.h"
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
@@ -34,10 +34,10 @@ Basis* BasisGenerator::generate(unsigned int elementType,
if(!family.compare(std::string("zaglmayr")))
return generateZaglmayr(elementType, basisType, order);
-
+ /*
else if(!family.compare(std::string("lagrange")))
return generateLagrange(elementType, basisType, order);
-
+ */
else
throw Exception("Unknwown Basis Family: %s", family.c_str());
}
@@ -56,7 +56,7 @@ Basis* BasisGenerator::generateZaglmayr(unsigned int elementType,
elementType);
}
}
-
+/*
Basis* BasisGenerator::generateLagrange(unsigned int elementType,
unsigned int basisType,
unsigned int order){
@@ -76,7 +76,7 @@ Basis* BasisGenerator::generateLagrange(unsigned int elementType,
elementType);
}
}
-
+*/
Basis* BasisGenerator::linZaglmayrGen(unsigned int basisType,
unsigned int order){
switch(basisType){
diff --git a/FunctionSpace/BasisGenerator.h b/FunctionSpace/BasisGenerator.h
index 2a6dc8c919..ee2584d2f8 100644
--- a/FunctionSpace/BasisGenerator.h
+++ b/FunctionSpace/BasisGenerator.h
@@ -40,10 +40,11 @@ class BasisGenerator{
static Basis* generateZaglmayr(unsigned int elementType,
unsigned int basisType,
unsigned int order);
-
+ /*
static Basis* generateLagrange(unsigned int elementType,
unsigned int basisType,
unsigned int order);
+ */
};
diff --git a/FunctionSpace/BasisHierarchicalScalar.cpp b/FunctionSpace/BasisHierarchicalScalar.cpp
new file mode 100644
index 0000000000..ea6946e9f0
--- /dev/null
+++ b/FunctionSpace/BasisHierarchicalScalar.cpp
@@ -0,0 +1,217 @@
+#include <sstream>
+#include "Exception.h"
+#include "BasisHierarchicalScalar.h"
+
+using namespace std;
+
+BasisHierarchicalScalar::BasisHierarchicalScalar(void){
+ // Grad Basis //
+ hasGrad = false;
+ grad = NULL;
+
+ // PreEvaluation //
+ preEvaluated = false;
+ preEvaluatedGrad = false;
+
+ preEvaluatedFunction = NULL;
+ preEvaluatedGradFunction = NULL;
+}
+
+BasisHierarchicalScalar::~BasisHierarchicalScalar(void){
+ // Grad Basis //
+ if(hasGrad){
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete grad[i][j];
+
+ delete[] grad[i];
+ }
+
+ delete[] grad;
+ }
+
+ // PreEvaluation //
+ if(preEvaluated){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedFunction[i];
+
+ delete[] preEvaluatedFunction;
+ }
+
+ if(preEvaluatedGrad){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedGradFunction[i];
+
+ delete[] preEvaluatedGradFunction;
+ }
+}
+
+fullMatrix<double>* BasisHierarchicalScalar::
+getFunctions(const MElement& element,
+ double u, double v, double w) const{
+
+ // Alloc Memory //
+ fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
+
+ // Define Permutation //
+ unsigned int permutation = refSpace->getPermutation(element);
+
+ // Fill Matrix //
+ for(unsigned int i = 0; i < nFunction; i++)
+ (*values)(i, 0) = basis[permutation][i]->at(u, v, w);
+
+ // Return //
+ return values;
+}
+
+fullMatrix<double>* BasisHierarchicalScalar::
+getFunctions(unsigned int permutation,
+ double u, double v, double w) const{
+
+ // Alloc Memory //
+ fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
+
+ // Fill Matrix //
+ for(unsigned int i = 0; i < nFunction; i++)
+ (*values)(i, 0) = basis[permutation][i]->at(u, v, w);
+
+ // Return //
+ return values;
+}
+
+void BasisHierarchicalScalar::
+preEvaluateFunctions(const fullMatrix<double>& point) const{
+ // Delete if older //
+ if(preEvaluated){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedFunction[i];
+
+ delete[] preEvaluatedFunction;
+ }
+
+ // Alloc //
+ const unsigned int nPoint = point.size1();
+ preEvaluatedFunction = new fullMatrix<double>*[nRefSpace];
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ preEvaluatedFunction[i] =
+ new fullMatrix<double>(nFunction, nPoint);
+
+ // Fill Matrix //
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ for(unsigned int j = 0; j < nFunction; j++)
+ for(unsigned int k = 0; k < nPoint; k++)
+ (*preEvaluatedFunction[i])(j, k) =
+ basis[i][j]->at(point(k, 0),
+ point(k, 1),
+ point(k, 2));
+
+ // PreEvaluated //
+ preEvaluated = true;
+}
+
+void BasisHierarchicalScalar::
+preEvaluateGradFunctions(const fullMatrix<double>& point) const{
+ // Build Grad //
+ if(!hasGrad)
+ getGrad();
+
+ // Delete if older //
+ if(preEvaluatedGrad){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedGradFunction[i];
+
+ delete[] preEvaluatedGradFunction;
+ }
+
+ // Alloc //
+ const unsigned int nPoint = point.size1();
+ const unsigned int nPoint3 = nPoint * 3;
+ preEvaluatedGradFunction = new fullMatrix<double>*[nRefSpace];
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ preEvaluatedGradFunction[i] =
+ new fullMatrix<double>(nFunction, nPoint3);
+
+ // Fill Matrix //
+ fullVector<double> tmp(3);
+
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++){
+ for(unsigned int k = 0; k < nPoint; k++){
+ tmp = Polynomial::at(*grad[i][j],
+ point(k, 0),
+ point(k, 1),
+ point(k, 2));
+
+ (*preEvaluatedGradFunction[i])(j, 3 * k) = tmp(0);
+ (*preEvaluatedGradFunction[i])(j, 3 * k + 1) = tmp(1);
+ (*preEvaluatedGradFunction[i])(j, 3 * k + 2) = tmp(2);
+ }
+ }
+ }
+
+ // PreEvaluated //
+ preEvaluatedGrad = true;
+}
+
+const fullMatrix<double>& BasisHierarchicalScalar::
+getPreEvaluatedFunctions(const MElement& element) const{
+ if(!preEvaluated)
+ throw Exception("getPreEvaluatedFunction: function has not been preEvaluated");
+
+ return *preEvaluatedFunction[refSpace->getPermutation(element)];
+}
+
+const fullMatrix<double>& BasisHierarchicalScalar::
+getPreEvaluatedGradFunctions(const MElement& element) const{
+ if(!preEvaluatedGrad)
+ throw Exception("getPreEvaluatedGradFunction: gradient has not been preEvaluated");
+
+ return *preEvaluatedGradFunction[refSpace->getPermutation(element)];
+}
+
+void BasisHierarchicalScalar::getGrad(void) const{
+ // Alloc //
+ grad = new vector<Polynomial>**[nRefSpace];
+
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ grad[s] = new vector<Polynomial>*[nFunction];
+
+ // Grad //
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ for(unsigned int f = 0 ; f < nFunction; f++)
+ grad[s][f] =
+ new vector<Polynomial>(basis[s][f]->gradient());
+
+ // Has Grad //
+ hasGrad = true;
+}
+
+string BasisHierarchicalScalar::toString(void) const{
+ stringstream stream;
+ unsigned int i = 0;
+ const unsigned int refSpace = 0;
+
+ stream << "Vertex Based:" << endl;
+ for(; i < nVertex; i++)
+ stream << "f(" << i + 1 << ") = "
+ << basis[refSpace][i]->toString() << endl;
+
+ stream << "Edge Based:" << endl;
+ for(; i < nVertex + nEdge; i++)
+ stream << "f(" << i + 1 << ") = "
+ << basis[refSpace][i]->toString() << endl;
+
+ stream << "Face Based:" << endl;
+ for(; i < nVertex + nEdge + nFace; i++)
+ stream << "f(" << i + 1 << ") = "
+ << basis[refSpace][i]->toString() << endl;
+
+ stream << "Cell Based:" << endl;
+ for(; i < nVertex + nEdge + nFace + nCell; i++)
+ stream << "f(" << i + 1 << ") = "
+ << basis[refSpace][i]->toString() << endl;
+
+ return stream.str();
+}
diff --git a/FunctionSpace/BasisHierarchicalScalar.h b/FunctionSpace/BasisHierarchicalScalar.h
new file mode 100644
index 0000000000..d8103a3276
--- /dev/null
+++ b/FunctionSpace/BasisHierarchicalScalar.h
@@ -0,0 +1,56 @@
+#ifndef _BASISHIERARCHICALSCALAR_H_
+#define _BASISHIERARCHICALSCALAR_H_
+
+#include <string>
+#include "BasisScalar.h"
+#include "Polynomial.h"
+#include "ReferenceSpace.h"
+
+class BasisHierarchicalScalar: public BasisScalar{
+ protected:
+ // Permutation //
+ ReferenceSpace* refSpace;
+ unsigned int nRefSpace;
+
+ // Basis //
+ Polynomial*** basis;
+
+ // Grad Basis //
+ mutable bool hasGrad;
+ mutable std::vector<Polynomial>*** grad;
+
+ // PreEvaluation //
+ mutable bool preEvaluated;
+ mutable bool preEvaluatedGrad;
+
+ mutable fullMatrix<double>** preEvaluatedFunction;
+ mutable fullMatrix<double>** preEvaluatedGradFunction;
+
+ public:
+ virtual ~BasisHierarchicalScalar(void);
+
+ virtual fullMatrix<double>* getFunctions(const MElement& element,
+ double u, double v, double w) const;
+
+ virtual fullMatrix<double>* getFunctions(unsigned int permutation,
+ double u, double v, double w) const;
+
+ virtual void preEvaluateFunctions (const fullMatrix<double>& point) const;
+ virtual void preEvaluateGradFunctions(const fullMatrix<double>& point) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedFunctions(const MElement& element) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedGradFunctions(const MElement& element) const;
+
+ std::string toString(void) const;
+
+ protected:
+ BasisHierarchicalScalar(void);
+
+ private:
+ void getGrad(void) const;
+};
+
+#endif
diff --git a/FunctionSpace/BasisHierarchicalVector.cpp b/FunctionSpace/BasisHierarchicalVector.cpp
new file mode 100644
index 0000000000..29e865eaf5
--- /dev/null
+++ b/FunctionSpace/BasisHierarchicalVector.cpp
@@ -0,0 +1,325 @@
+#include <sstream>
+#include "Exception.h"
+#include "BasisHierarchicalVector.h"
+
+using namespace std;
+
+BasisHierarchicalVector::BasisHierarchicalVector(void){
+ // Curl Basis //
+ hasCurl = false;
+ curl = NULL;
+
+ // Div Basis //
+ hasDiv = false;
+ div = NULL;
+
+ // PreEvaluation //
+ preEvaluated = false;
+ preEvaluatedCurl = false;
+ preEvaluatedDiv = false;
+
+ preEvaluatedFunction = NULL;
+ preEvaluatedCurlFunction = NULL;
+ preEvaluatedDivFunction = NULL;
+}
+
+BasisHierarchicalVector::~BasisHierarchicalVector(void){
+ // Curl Basis //
+ if(hasCurl){
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete curl[i][j];
+
+ delete[] curl[i];
+ }
+
+ delete[] curl;
+ }
+
+ // Div Basis //
+ if(hasDiv){
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete div[i][j];
+
+ delete[] div[i];
+ }
+
+ delete[] div;
+ }
+
+ // PreEvaluation //
+ if(preEvaluated){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedFunction[i];
+
+ delete[] preEvaluatedFunction;
+ }
+
+ if(preEvaluatedCurl){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedCurlFunction[i];
+
+ delete[] preEvaluatedCurlFunction;
+ }
+}
+
+fullMatrix<double>* BasisHierarchicalVector::
+getFunctions(const MElement& element,
+ double u, double v, double w) const{
+
+ // Alloc Memory //
+ fullMatrix<double>* values = new fullMatrix<double>(nFunction, 3);
+
+ // Define Permutation //
+ unsigned int permutation = refSpace->getPermutation(element);
+
+ // Fill Vector //
+ for(unsigned int i = 0; i < nFunction; i++){
+ fullVector<double> eval =
+ Polynomial::at(*basis[permutation][i], u, v, w);
+
+ (*values)(i, 0) = eval(0);
+ (*values)(i, 1) = eval(1);
+ (*values)(i, 2) = eval(2);
+ }
+
+ // Return //
+ return values;
+}
+
+fullMatrix<double>* BasisHierarchicalVector::
+getFunctions(unsigned int permutation,
+ double u, double v, double w) const{
+
+ // Alloc Memory //
+ fullMatrix<double>* values = new fullMatrix<double>(nFunction, 3);
+
+ // Fill Vector //
+ for(unsigned int i = 0; i < nFunction; i++){
+ fullVector<double> eval =
+ Polynomial::at(*basis[permutation][i], u, v, w);
+
+ (*values)(i, 0) = eval(0);
+ (*values)(i, 1) = eval(1);
+ (*values)(i, 2) = eval(2);
+ }
+
+ // Return //
+ return values;
+}
+
+void BasisHierarchicalVector::
+preEvaluateFunctions(const fullMatrix<double>& point) const{
+ // Delete if older //
+ if(preEvaluated){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedFunction[i];
+
+ delete[] preEvaluatedFunction;
+ }
+
+ // Alloc //
+ const unsigned int nPoint = point.size1();
+ const unsigned int nPoint3 = nPoint * 3;
+ preEvaluatedFunction = new fullMatrix<double>*[nRefSpace];
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ preEvaluatedFunction[i] =
+ new fullMatrix<double>(nFunction, nPoint3);
+
+ // Fill Matrix //
+ fullVector<double> tmp(3);
+
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++){
+ for(unsigned int k = 0; k < nPoint; k++){
+ tmp = Polynomial::at(*basis[i][j],
+ point(k, 0),
+ point(k, 1),
+ point(k, 2));
+
+ (*preEvaluatedFunction[i])(j, 3 * k) = tmp(0);
+ (*preEvaluatedFunction[i])(j, 3 * k + 1) = tmp(1);
+ (*preEvaluatedFunction[i])(j, 3 * k + 2) = tmp(2);
+ }
+ }
+ }
+
+ // PreEvaluated //
+ preEvaluated = true;
+}
+
+void BasisHierarchicalVector::
+preEvaluateCurlFunctions(const fullMatrix<double>& point) const{
+ // Build Curl //
+ if(!hasCurl)
+ getCurl();
+
+ // Delete if older //
+ if(preEvaluatedCurl){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedCurlFunction[i];
+
+ delete[] preEvaluatedCurlFunction;
+ }
+
+ // Alloc //
+ const unsigned int nPoint = point.size1();
+ const unsigned int nPoint3 = nPoint * 3;
+ preEvaluatedCurlFunction = new fullMatrix<double>*[nRefSpace];
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ preEvaluatedCurlFunction[i] =
+ new fullMatrix<double>(nFunction, nPoint3);
+
+ // Fill Matrix //
+ fullVector<double> tmp(3);
+
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++){
+ for(unsigned int k = 0; k < nPoint; k++){
+ tmp = Polynomial::at(*curl[i][j],
+ point(k, 0),
+ point(k, 1),
+ point(k, 2));
+
+ (*preEvaluatedCurlFunction[i])(j, 3 * k) = tmp(0);
+ (*preEvaluatedCurlFunction[i])(j, 3 * k + 1) = tmp(1);
+ (*preEvaluatedCurlFunction[i])(j, 3 * k + 2) = tmp(2);
+ }
+ }
+ }
+
+ // PreEvaluated //
+ preEvaluatedCurl = true;
+}
+
+void BasisHierarchicalVector::
+preEvaluateDivFunctions(const fullMatrix<double>& point) const{
+ // Build Div //
+ if(!hasDiv)
+ getDiv();
+
+ // Delete if older //
+ if(preEvaluatedDiv){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete preEvaluatedDivFunction[i];
+
+ delete[] preEvaluatedDivFunction;
+ }
+
+ // Alloc //
+ const unsigned int nPoint = point.size1();
+ preEvaluatedDivFunction = new fullMatrix<double>*[nRefSpace];
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ preEvaluatedDivFunction[i] =
+ new fullMatrix<double>(nFunction, nPoint);
+
+ // Fill Matrix //
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ for(unsigned int j = 0; j < nFunction; j++)
+ for(unsigned int k = 0; k < nPoint; k++)
+ (*preEvaluatedDivFunction[i])(j, k) =
+ div[i][j]->at(point(k, 0),
+ point(k, 1),
+ point(k, 2));
+
+ // PreEvaluated //
+ preEvaluatedDiv = true;
+}
+
+const fullMatrix<double>& BasisHierarchicalVector::
+getPreEvaluatedFunctions(const MElement& element) const{
+ if(!preEvaluated)
+ throw Exception("getPreEvaluatedFunction: function has not been preEvaluated");
+
+ return *preEvaluatedFunction[refSpace->getPermutation(element)];
+}
+
+const fullMatrix<double>& BasisHierarchicalVector::
+getPreEvaluatedCurlFunctions(const MElement& element) const{
+ if(!preEvaluatedCurl)
+ throw Exception("getPreEvaluatedCurlFunction: curl has not been preEvaluated");
+
+ return *preEvaluatedCurlFunction[refSpace->getPermutation(element)];
+}
+
+const fullMatrix<double>& BasisHierarchicalVector::
+getPreEvaluatedDivFunctions(const MElement& element) const{
+ if(!preEvaluatedDiv)
+ throw Exception("getPreEvaluatedDivFunction: divergence has not been preEvaluated");
+
+ return *preEvaluatedDivFunction[refSpace->getPermutation(element)];
+}
+
+void BasisHierarchicalVector::getCurl(void) const{
+ // Alloc //
+ curl = new vector<Polynomial>**[nRefSpace];
+
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ curl[s] = new vector<Polynomial>*[nFunction];
+
+ // Curl //
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ for(unsigned int f = 0 ; f < nFunction; f++)
+ curl[s][f] =
+ new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
+
+ // Has Curl //
+ hasCurl = true;
+}
+
+void BasisHierarchicalVector::getDiv(void) const{
+ // Alloc //
+ div = new Polynomial**[nRefSpace];
+
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ div[s] = new Polynomial*[nFunction];
+
+ // Div //
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ for(unsigned int f = 0 ; f < nFunction; f++)
+ div[s][f] =
+ new Polynomial(Polynomial::divergence(*basis[s][f]));
+
+ // Has Div //
+ hasDiv = true;
+}
+
+string BasisHierarchicalVector::toString(void) const{
+ stringstream stream;
+ unsigned int i = 0;
+ const unsigned int refSpace = 0;
+
+ stream << "Vertex Based:" << endl;
+ for(; i < nVertex; i++)
+ stream << "f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
+
+ stream << "Edge Based:" << endl;
+ for(; i < nVertex + nEdge; i++)
+ stream << " f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
+
+ stream << "Face Based:" << endl;
+ for(; i < nVertex + nEdge + nFace; i++)
+ stream << " f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
+
+ stream << "Cell Based:" << endl;
+ for(; i < nVertex + nEdge + nFace + nCell; i++)
+ stream << " f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]" << endl
+ << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]" << endl;
+
+ return stream.str();
+}
diff --git a/FunctionSpace/BasisHierarchicalVector.h b/FunctionSpace/BasisHierarchicalVector.h
new file mode 100644
index 0000000000..f951fbc84a
--- /dev/null
+++ b/FunctionSpace/BasisHierarchicalVector.h
@@ -0,0 +1,68 @@
+#ifndef _BASISHIERARCHICALVECTOR_H_
+#define _BASISHIERARCHICALVECTOR_H_
+
+#include <string>
+#include "BasisVector.h"
+#include "Polynomial.h"
+#include "ReferenceSpace.h"
+
+class BasisHierarchicalVector: public BasisVector{
+ protected:
+ // Permutation //
+ ReferenceSpace* refSpace;
+ unsigned int nRefSpace;
+
+ // Basis //
+ std::vector<Polynomial>*** basis;
+
+ // Curl Basis //
+ mutable bool hasCurl;
+ mutable std::vector<Polynomial>*** curl;
+
+ // Div Basis //
+ mutable bool hasDiv;
+ mutable Polynomial*** div;
+
+ // PreEvaluation //
+ mutable bool preEvaluated;
+ mutable bool preEvaluatedCurl;
+ mutable bool preEvaluatedDiv;
+
+ mutable fullMatrix<double>** preEvaluatedFunction;
+ mutable fullMatrix<double>** preEvaluatedCurlFunction;
+ mutable fullMatrix<double>** preEvaluatedDivFunction;
+
+ public:
+ virtual ~BasisHierarchicalVector(void);
+
+ virtual fullMatrix<double>* getFunctions(const MElement& element,
+ double u, double v, double w) const;
+
+ virtual fullMatrix<double>* getFunctions(unsigned int permutation,
+ double u, double v, double w) const;
+
+ virtual void preEvaluateFunctions (const fullMatrix<double>& point) const;
+ virtual void preEvaluateCurlFunctions(const fullMatrix<double>& point) const;
+ virtual void preEvaluateDivFunctions (const fullMatrix<double>& point) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedFunctions(const MElement& element) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedCurlFunctions(const MElement& element) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedDivFunctions (const MElement& element)const;
+
+ std::string toString(void) const;
+
+ protected:
+ BasisHierarchicalVector(void);
+
+ private:
+ void getCurl(void) const;
+ void getDiv(void) const;
+};
+
+#endif
+
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
new file mode 100644
index 0000000000..47f1ae99f2
--- /dev/null
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -0,0 +1,50 @@
+#include "BasisLagrange.h"
+
+BasisLagrange::BasisLagrange(void){
+}
+
+BasisLagrange::~BasisLagrange(void){
+}
+
+fullMatrix<double>* BasisLagrange::
+getFunctions(const MElement& element,
+ double u, double v, double w) const{
+ return new fullMatrix<double>(nFunction, 1);
+}
+
+fullMatrix<double>* BasisLagrange::getFunctions(unsigned int permutation,
+ double u, double v, double w) const{
+
+ // Alloc Memory //
+ fullMatrix<double> tmp(nFunction, 1);
+ fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
+
+ // Fill Matrix //
+ fullMatrix<double> point(1, 3);
+ point(0, 0) = u;
+ point(0, 1) = v;
+ point(0, 2) = w;
+
+ basis->f(point, tmp);
+
+ *values = tmp.transpose(); // Otherwise not coherent with df !!
+
+ // Return //
+ return values;
+}
+
+void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
+}
+
+void BasisLagrange::preEvaluateGradFunctions(const fullMatrix<double>& point) const{
+}
+
+const fullMatrix<double>&
+BasisLagrange::getPreEvaluatedFunctions(const MElement& element) const{
+ return fullMatrix<double>(nFunction, 1);
+}
+
+const fullMatrix<double>&
+BasisLagrange::getPreEvaluatedGradFunctions(const MElement& element) const{
+ return fullMatrix<double>(nFunction, 1);
+}
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
new file mode 100644
index 0000000000..ee9fa43bf4
--- /dev/null
+++ b/FunctionSpace/BasisLagrange.h
@@ -0,0 +1,33 @@
+#ifndef _BASISLAGRANGE_H_
+#define _BASISLAGRANGE_H_
+
+#include "BasisScalar.h"
+#include "polynomialBasis.h"
+
+class BasisLagrange: public BasisScalar{
+ protected:
+ polynomialBasis* basis;
+
+ public:
+ virtual ~BasisLagrange(void);
+
+ virtual fullMatrix<double>* getFunctions(const MElement& element,
+ double u, double v, double w) const;
+
+ virtual fullMatrix<double>* getFunctions(unsigned int permutation,
+ double u, double v, double w) const;
+
+ virtual void preEvaluateFunctions (const fullMatrix<double>& point) const;
+ virtual void preEvaluateGradFunctions(const fullMatrix<double>& point) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedFunctions(const MElement& element) const;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedGradFunctions(const MElement& element) const;
+
+ protected:
+ BasisLagrange(void);
+};
+
+#endif
diff --git a/FunctionSpace/BasisScalar.cpp b/FunctionSpace/BasisScalar.cpp
index 55a09412a6..e902527ff8 100644
--- a/FunctionSpace/BasisScalar.cpp
+++ b/FunctionSpace/BasisScalar.cpp
@@ -1,39 +1,8 @@
-#include <sstream>
#include "BasisScalar.h"
-using namespace std;
-
BasisScalar::BasisScalar(void){
scalar = true;
}
BasisScalar::~BasisScalar(void){
}
-
-string BasisScalar::toString(void) const{
- stringstream stream;
- unsigned int i = 0;
- const unsigned int refSpace = 0;
-
- stream << "Vertex Based:" << endl;
- for(; i < nVertex; i++)
- stream << "f(" << i + 1 << ") = "
- << (*(*basis)[refSpace])[i]->toString() << endl;
-
- stream << "Edge Based:" << endl;
- for(; i < nVertex + nEdge; i++)
- stream << "f(" << i + 1 << ") = "
- << (*(*basis)[refSpace])[i]->toString() << endl;
-
- stream << "Face Based:" << endl;
- for(; i < nVertex + nEdge + nFace; i++)
- stream << "f(" << i + 1 << ") = "
- << (*(*basis)[refSpace])[i]->toString() << endl;
-
- stream << "Cell Based:" << endl;
- for(; i < nVertex + nEdge + nFace + nCell; i++)
- stream << "f(" << i + 1 << ") = "
- << (*(*basis)[refSpace])[i]->toString() << endl;
-
- return stream.str();
-}
diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h
index 01e295baa0..fa9cad9be8 100644
--- a/FunctionSpace/BasisScalar.h
+++ b/FunctionSpace/BasisScalar.h
@@ -1,9 +1,9 @@
#ifndef _BASISSCALAR_H_
#define _BASISSCALAR_H_
-#include <vector>
#include "Basis.h"
-#include "Polynomial.h"
+#include "MElement.h"
+#include "fullMatrix.h"
/**
@interface BasisScalar
@@ -12,42 +12,28 @@
This class is the @em common @em interface for all
@em scalar Basis.@n
-
- @note
- A BasisScalar is an @em interface,
- so it @em can't be instanciated
*/
class BasisScalar: public Basis{
- protected:
- std::vector<std::vector<const Polynomial*>*>* basis;
-
public:
//! Deletes this BasisScalar
//!
virtual ~BasisScalar(void);
- //! @param refSpace A natural number
- //! @param i A natural number
- //! @return Returns the @c i%th @em
- //! Basis Function of the @c refSpace%th ReferenceSpace
- const Polynomial&
- getFunction(unsigned int refSpace, unsigned int i) const;
+ virtual fullMatrix<double>* getFunctions(const MElement& element,
+ double u, double v, double w) const = 0;
- //! @param refSpace A natural number
- //! @return Returns the @em all
- //! Basis Function of the @c refSpace%th ReferenceSpace
- const std::vector<const Polynomial*>&
- getFunction(unsigned int refSpace) const;
+ virtual fullMatrix<double>* getFunctions(unsigned int permutation,
+ double u, double v, double w) const = 0;
- //! @param element An Element
- //! @return Returns the @em all
- //! Basis Function in the @c given element
- //! @em ReferenceSpace
- const std::vector<const Polynomial*>&
- getFunction(const MElement& element) const;
+ virtual void preEvaluateFunctions (const fullMatrix<double>& point) const = 0;
+ virtual void preEvaluateGradFunctions(const fullMatrix<double>& point) const = 0;
- virtual std::string toString(void) const;
+ virtual const fullMatrix<double>&
+ getPreEvaluatedFunctions(const MElement& element) const = 0;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedGradFunctions(const MElement& element) const = 0;
protected:
//! @internal
@@ -57,26 +43,4 @@ class BasisScalar: public Basis{
BasisScalar(void);
};
-//////////////////////
-// Inline Function //
-//////////////////////
-
-inline
-const Polynomial&
-BasisScalar::getFunction(unsigned int refSpace, unsigned int i) const{
- return *(*(*basis)[refSpace])[i];
-}
-
-inline
-const std::vector<const Polynomial*>&
-BasisScalar::getFunction(unsigned int refSpace) const{
- return *(*basis)[refSpace];
-}
-
-inline
-const std::vector<const Polynomial*>&
-BasisScalar::getFunction(const MElement& element) const{
- return *(*basis)[refSpace->getPermutation(element)];
-}
-
#endif
diff --git a/FunctionSpace/BasisVector.cpp b/FunctionSpace/BasisVector.cpp
index c331c191d5..e1d389822e 100644
--- a/FunctionSpace/BasisVector.cpp
+++ b/FunctionSpace/BasisVector.cpp
@@ -1,47 +1,8 @@
-#include <sstream>
#include "BasisVector.h"
-using namespace std;
-
BasisVector::BasisVector(void){
scalar = false;
}
BasisVector::~BasisVector(void){
}
-
-string BasisVector::toString(void) const{
- stringstream stream;
- unsigned int i = 0;
- const unsigned int refSpace = 0;
-
- stream << "Vertex Based:" << endl;
- for(; i < nVertex; i++)
- stream << "f(" << i + 1 << ") = " << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
-
- stream << "Edge Based:" << endl;
- for(; i < nVertex + nEdge; i++)
- stream << " f(" << i + 1 << ") = " << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
-
- stream << "Face Based:" << endl;
- for(; i < nVertex + nEdge + nFace; i++)
- stream << " f(" << i + 1 << ") = " << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
-
- stream << "Cell Based:" << endl;
- for(; i < nVertex + nEdge + nFace + nCell; i++)
- stream << " f(" << i + 1 << ") = " << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[0].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[1].toString() << " ]" << endl
- << "\t[ " << (*(*(*basis)[refSpace])[i])[2].toString() << " ]" << endl;
-
- return stream.str();
-}
diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h
index eefdead46c..cf2615d515 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 "MElement.h"
+#include "fullMatrix.h"
/**
@interface BasisVector
@@ -12,71 +12,39 @@
This class is the @em common @em interface for all
@em vectorial Basis.@n
-
- @note
- A BasisVector is an @em interface,
- so it @em can't be instanciated
*/
class BasisVector: public Basis{
- protected:
- std::vector<std::vector<const std::vector<Polynomial>*>*>* basis;
-
public:
//! Deletes this BasisVector
//!
virtual ~BasisVector(void);
- //! @param refSpace A natural number
- //! @param i A natural number
- //! @return Returns the @c i%th @em
- //! Basis Function of the @c refSpace%th ReferenceSpace
- const std::vector<Polynomial>&
- getFunction(unsigned int refSpace, unsigned int i) const;
+ virtual fullMatrix<double>* getFunctions(const MElement& element,
+ double u, double v, double w) const = 0;
+
+ virtual fullMatrix<double>* getFunctions(unsigned int permutation,
+ double u, double v, double w) const = 0;
- //! @param refSpace A natural number
- //! @return Returns the @em all
- //! Basis Function of the @c refSpace%th ReferenceSpace
- const std::vector<const std::vector<Polynomial>*>&
- getFunction(unsigned int refSpace) const;
+ virtual void preEvaluateFunctions (const fullMatrix<double>& point) const = 0;
+ virtual void preEvaluateCurlFunctions(const fullMatrix<double>& point) const = 0;
+ virtual void preEvaluateDivFunctions (const fullMatrix<double>& point) const = 0;
+
+ virtual const fullMatrix<double>&
+ getPreEvaluatedFunctions(const MElement& element) const = 0;
- //! @param element An Element
- //! @return Returns the @em all
- //! Basis Function in the @c given element
- //! @em ReferenceSpace
- const std::vector<const std::vector<Polynomial>*>&
- getFunction(const MElement& element) const;
+ virtual const fullMatrix<double>&
+ getPreEvaluatedCurlFunctions(const MElement& element) const = 0;
- virtual std::string toString(void) const;
+ virtual const fullMatrix<double>&
+ getPreEvaluatedDivFunctions(const MElement& element) const = 0;
protected:
//! @internal
- //! Instantiate a new BasisVector
+ //! Instantiates a new BasisVector
//!
//! @endinternal
BasisVector(void);
};
-//////////////////////
-// Inline Functions //
-//////////////////////
-
-inline
-const std::vector<Polynomial>&
-BasisVector::getFunction(unsigned int refSpace, unsigned int i) const{
- return *(*(*basis)[refSpace])[i];
-}
-
-inline
-const std::vector<const std::vector<Polynomial>*>&
-BasisVector::getFunction(unsigned int refSpace) const{
- return *(*basis)[refSpace];
-}
-
-inline
-const std::vector<const std::vector<Polynomial>*>&
-BasisVector::getFunction(const MElement& element) const{
- return *(*basis)[refSpace->getPermutation(element)];
-}
-
#endif
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index 932379819d..21a18aa747 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -13,16 +13,13 @@ set(SRC
TetReferenceSpace.cpp
Basis.cpp
- LagrangeBasis.cpp
BasisScalar.cpp
BasisVector.cpp
BasisGenerator.cpp
- GradBasis.cpp
- CurlBasis.cpp
- DivBasis.cpp
-
- EvaluatedBasis.cpp
+ BasisLagrange.cpp
+ BasisHierarchicalScalar.cpp
+ BasisHierarchicalVector.cpp
LineNodeBasis.cpp
LineEdgeBasis.cpp
diff --git a/FunctionSpace/CurlBasis.cpp b/FunctionSpace/CurlBasis.cpp
deleted file mode 100644
index ba2c320185..0000000000
--- a/FunctionSpace/CurlBasis.cpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "CurlBasis.h"
-
-using namespace std;
-
-CurlBasis::CurlBasis(const BasisVector& other){
- // Reference Space //
- refSpace = &other.getReferenceSpace();
- nRefSpace = other.getReferenceSpace().getNPermutation();
-
- // Set Basis Type //
- order = other.getOrder() - 1;
-
- type = other.getType();
- dim = other.getDim();
-
- nVertex = other.getNVertexBased();
- nEdge = other.getNEdgeBased();
- nFace = other.getNFaceBased();
- nCell = other.getNCellBased();
- nFunction = other.getNFunction();
-
- // Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- for(unsigned int i = 0; i < nFunction; i++)
- (*(*basis)[s])[i] =
- new vector<Polynomial>
- (Polynomial::curl(other.getFunction(s, i)));
-}
-
-CurlBasis::~CurlBasis(void){
- // Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
-
- delete (*basis)[i];
- }
-
- delete basis;
-}
diff --git a/FunctionSpace/CurlBasis.h b/FunctionSpace/CurlBasis.h
deleted file mode 100644
index ffedacf3fe..0000000000
--- a/FunctionSpace/CurlBasis.h
+++ /dev/null
@@ -1,26 +0,0 @@
-#ifndef _CURLBASIS_H_
-#define _CURLBASIS_H_
-
-#include "BasisVector.h"
-
-/**
- @class CurlBasis
- @brief The Curl of an other Basis
-
- This class can instantiate a Curl Basis.
- */
-
-class CurlBasis: public BasisVector{
- public:
- //! @param other An Other Basis
- //!
- //! Returns a new Basis, which is the curl of
- //! the given Basis
- CurlBasis(const BasisVector& other);
-
- //! Deletes this Basis
- //!
- virtual ~CurlBasis(void);
-};
-
-#endif
diff --git a/FunctionSpace/DivBasis.cpp b/FunctionSpace/DivBasis.cpp
deleted file mode 100644
index 774d411da5..0000000000
--- a/FunctionSpace/DivBasis.cpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "DivBasis.h"
-
-using namespace std;
-
-DivBasis::DivBasis(const BasisVector& other){
- // Reference Space //
- refSpace = &other.getReferenceSpace();
- nRefSpace = other.getReferenceSpace().getNPermutation();
-
- // Set Basis Type //
- order = other.getOrder() - 1;
-
- type = other.getType();
- dim = other.getDim();
-
- nVertex = other.getNVertexBased();
- nEdge = other.getNEdgeBased();
- nFace = other.getNFaceBased();
- nCell = other.getNCellBased();
- nFunction = other.getNFunction();
-
- // Basis //
- basis = new vector<vector<const Polynomial*>*>(nRefSpace);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const Polynomial*>(nFunction);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- for(unsigned int i = 0; i < nFunction; i++)
- (*(*basis)[s])[i] =
- new Polynomial
- (Polynomial::divergence(other.getFunction(s, i)));
-}
-
-DivBasis::~DivBasis(void){
- // Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
-
- delete (*basis)[i];
- }
-
- delete basis;
-}
diff --git a/FunctionSpace/DivBasis.h b/FunctionSpace/DivBasis.h
deleted file mode 100644
index 2ebf45e9bd..0000000000
--- a/FunctionSpace/DivBasis.h
+++ /dev/null
@@ -1,27 +0,0 @@
-#ifndef _DIVBASIS_H_
-#define _DIVBASIS_H_
-
-#include "BasisScalar.h"
-#include "BasisVector.h"
-
-/**
- @class DivBasis
- @brief The Divergence of an other Basis
-
- This class can instantiate a Divergence Basis.
- */
-
-class DivBasis: public BasisScalar{
- public:
- //! @param other An Other Basis
- //!
- //! Returns a new Basis, which is the divergence of
- //! the given Basis
- DivBasis(const BasisVector& other);
-
- //! Deletes this Basis
- //!
- virtual ~DivBasis(void);
-};
-
-#endif
diff --git a/FunctionSpace/EvaluatedBasis.cpp b/FunctionSpace/EvaluatedBasis.cpp
deleted file mode 100644
index 4bc2debdda..0000000000
--- a/FunctionSpace/EvaluatedBasis.cpp
+++ /dev/null
@@ -1,75 +0,0 @@
-#include "Polynomial.h"
-#include "EvaluatedBasis.h"
-
-using namespace std;
-
-EvaluatedBasis::
-EvaluatedBasis(const BasisScalar& basis,
- const fullMatrix<double>& point){
- // Data //
- refSpace = &basis.getReferenceSpace();
- nRefSpace = refSpace->getNPermutation();
- scalar = true;
- nFunction = basis.getNFunction();
- nPoint = point.size1();
-
- // Alloc //
- eBasis = new vector<fullMatrix<double>*>(nRefSpace);
-
- for(unsigned int i = 0; i < nRefSpace; i++)
- (*eBasis)[i] = new fullMatrix<double>(nFunction, nPoint);
-
- // Evaluate //
- for(unsigned int i = 0; i < nRefSpace; i++)
- for(unsigned int j = 0; j < nFunction; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*eBasis)[i])(j, k) =
- basis.getFunction(i, j).at(point(k, 0),
- point(k, 1),
- point(k, 2));
-}
-
-EvaluatedBasis::
-EvaluatedBasis(const BasisVector& basis,
- const fullMatrix<double>& point){
- // Data //
- refSpace = &basis.getReferenceSpace();
- nRefSpace = refSpace->getNPermutation();
- scalar = false;
- nFunction = basis.getNFunction();
- nPoint = point.size1();
-
- // Alloc //
- eBasis = new vector<fullMatrix<double>*>(nRefSpace);
-
- // WARNING Each Evaluation Got *3* Component !
- const unsigned int nPointThree = nPoint * 3;
-
- for(unsigned int i = 0; i < nRefSpace; i++)
- (*eBasis)[i] = new fullMatrix<double>(nFunction, nPointThree);
-
- // Evaluate //
- fullVector<double> tmp(3);
-
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++){
- for(unsigned int k = 0; k < nPoint; k++){
- tmp = Polynomial::at(basis.getFunction(i, j),
- point(k, 0),
- point(k, 1),
- point(k, 2));
-
- (*(*eBasis)[i])(j, 3 * k) = tmp(0);
- (*(*eBasis)[i])(j, 3 * k + 1) = tmp(1);
- (*(*eBasis)[i])(j, 3 * k + 2) = tmp(2);
- }
- }
- }
-}
-
-EvaluatedBasis::~EvaluatedBasis(void){
- for(unsigned int i = 0; i < nRefSpace; i++)
- delete (*eBasis)[i];
-
- delete eBasis;
-}
diff --git a/FunctionSpace/EvaluatedBasis.h b/FunctionSpace/EvaluatedBasis.h
deleted file mode 100644
index 3ff8a8d2aa..0000000000
--- a/FunctionSpace/EvaluatedBasis.h
+++ /dev/null
@@ -1,106 +0,0 @@
-#ifndef _EVALUATEDBASIS_H_
-#define _EVALUATEDBASIS_H_
-
-#include <vector>
-#include "fullMatrix.h"
-#include "BasisScalar.h"
-#include "BasisVector.h"
-#include "MElement.h"
-#include "ReferenceSpace.h"
-
-/**
- @class EvaluatedBasis
- @brief A Basis that has been Evaluated
-
- This class is A basis that has been Evaluated at
- some given points.@n
- */
-
-class EvaluatedBasis{
- private:
- const ReferenceSpace* refSpace;
- unsigned int nRefSpace;
- bool scalar;
-
- unsigned int nFunction;
- unsigned int nPoint;
-
- std::vector<fullMatrix<double>*>* eBasis;
-
- public:
- //! @param basis the Basis to Evaluate
- //! @param point the Evaluation Points
- //!
- //! Instanciates the requested EvaluatedBasisScalar
- //!
- EvaluatedBasis(const BasisScalar& basis,
- const fullMatrix<double>& point);
-
- //! @param basis the Basis to Evaluate
- //! @param point the Evaluation Points
- //!
- //! Instanciates the requested EvaluatedBasisVector
- //!
- EvaluatedBasis(const BasisVector& basis,
- const fullMatrix<double>& point);
-
- //! Deletes this EvaluatedBasis
- //!
- ~EvaluatedBasis(void);
-
- //! @return Returns:
- //! @li @c true, if the evaluated
- //! values are @em scalar
- //! @li @c false, otherwise
- bool isScalar(void) const;
-
- //! @return Returns the number of
- //! evaluated @em Functions
- unsigned int getNFunction(void) const;
-
- //! @return Returns the number of
- //! evaluation @em Points
- unsigned int getNPoint(void) const;
-
- //! @param refSpace A natural number
- //! @return Returns the evaluation of all Basis Functions
- //! (at every evaluation Points)
- //! for the @c refSpace%th @em ReferenceSpace
- const fullMatrix<double>&
- getEvaluation(unsigned int refSpace) const;
-
- //! @param element A MElement
- //! @return Returns the evaluation of all Basis Functions
- //! (at every evaluation Points)
- //! in the ReferenceSpace of the given element
- const fullMatrix<double>&
- getEvaluation(const MElement& element) const;
-};
-
-//////////////////////
-// Inline Function //
-//////////////////////
-
-inline bool EvaluatedBasis::isScalar(void) const{
- return scalar;
-}
-
-inline unsigned int EvaluatedBasis::getNFunction(void) const{
- return nFunction;
-}
-
-inline unsigned int EvaluatedBasis::getNPoint(void) const{
- return nPoint;
-}
-
-inline const fullMatrix<double>&
-EvaluatedBasis::getEvaluation(unsigned int refSpace) const{
- return *(*eBasis)[refSpace];
-}
-
-inline const fullMatrix<double>&
-EvaluatedBasis::getEvaluation(const MElement& element) const{
- return *(*eBasis)[refSpace->getPermutation(element)];
-}
-
-#endif
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index a68b2a9ad9..713be01704 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -177,32 +177,32 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
int it = 0;
// Add Vertex Based Dof //
- for(int i = 0; i < nFVertex; i++){
- for(int j = 0; j < nVertex; j++){
- myDof[it].setDof(mesh->getGlobalId(*vertex[j]), i);
+ for(int i = 0; i < nVertex; i++){
+ for(int j = 0; j < nFVertex; j++){
+ myDof[it].setDof(mesh->getGlobalId(*vertex[i]), j);
it++;
}
}
// Add Edge Based Dof //
- for(int i = 0; i < nFEdge; i++){
- for(int j = 0; j < nEdge; j++){
- myDof[it].setDof(mesh->getGlobalId(edge[j]), i);
+ for(int i = 0; i < nEdge; i++){
+ for(int j = 0; j < nFEdge; j++){
+ myDof[it].setDof(mesh->getGlobalId(edge[i]), j);
it++;
}
}
// Add Face Based Dof //
- for(int i = 0; i < nFFace; i++){
- for(int j = 0; j < nFace; j++){
- myDof[it].setDof(mesh->getGlobalId(face[j]), i);
+ for(int i = 0; i < nFace; i++){
+ for(int j = 0; j < nFFace; j++){
+ myDof[it].setDof(mesh->getGlobalId(face[i]), j);
it++;
}
}
// Add Cell Based Dof //
- for(int i = 0; i < nFCell; i++){
- myDof[it].setDof(mesh->getGlobalId(element), i);
+ for(int j = 0; j < nFCell; j++){
+ myDof[it].setDof(mesh->getGlobalId(element), j);
it++;
}
@@ -220,7 +220,3 @@ const GroupOfDof& FunctionSpace::getGoDFromElement(const MElement& element) cons
else
return *(it->second);
}
-
-string FunctionSpace::toString(void) const{
- return basis->toString();
-}
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index f73d663cff..b9951ee070 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -3,12 +3,10 @@
#include <map>
#include <vector>
-#include <string>
#include "Basis.h"
#include "BasisScalar.h"
#include "BasisVector.h"
-#include "EvaluatedBasis.h"
#include "Comparators.h"
#include "Dof.h"
@@ -92,8 +90,6 @@ class FunctionSpace{
unsigned int dofNumber(void) const;
unsigned int groupNumber(void) const;
- std::string toString(void) const;
-
protected:
// Init
FunctionSpace(void);
@@ -104,21 +100,6 @@ class FunctionSpace{
// Dof
void buildDof(void);
void insertDof(Dof& d, GroupOfDof* god);
-
- // Local Basis
- static
- const std::vector<const Polynomial*>&
- locBasis(const MElement& element,
- const BasisScalar& basis);
- static
- const std::vector<const std::vector<Polynomial>*>&
- locBasis(const MElement& element,
- const BasisVector& basis);
-
- static
- const fullMatrix<double>&
- locEvalBasis(const MElement& element,
- const EvaluatedBasis& evalBasis);
};
@@ -229,10 +210,6 @@ class FunctionSpace{
@return Returns the number of GroupOfDof%s
given by FunctionSpace::getAllGroups()
**
-
- @fn FunctionSpace::toString
- @return Returns the FunctionSpace string
- **
*/
@@ -288,25 +265,4 @@ inline const std::vector<GroupOfDof*>& FunctionSpace::getAllGroups(void) const{
return *group;
}
-inline const std::vector<const Polynomial*>&
-FunctionSpace::locBasis(const MElement& element,
- const BasisScalar& basis){
-
- return basis.getFunction(element);
-}
-
-inline const std::vector<const std::vector<Polynomial>*>&
-FunctionSpace::locBasis(const MElement& element,
- const BasisVector& basis){
-
- return basis.getFunction(element);
-}
-
-inline const fullMatrix<double>&
-FunctionSpace::locEvalBasis(const MElement& element,
- const EvaluatedBasis& evalBasis){
-
- return evalBasis.getEvaluation(element);
-}
-
#endif
diff --git a/FunctionSpace/FunctionSpaceEdge.cpp b/FunctionSpace/FunctionSpaceEdge.cpp
index dd3bd0dda0..5047d3bd9c 100644
--- a/FunctionSpace/FunctionSpaceEdge.cpp
+++ b/FunctionSpace/FunctionSpaceEdge.cpp
@@ -23,15 +23,11 @@ fullVector<double> FunctionSpaceEdge::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
-
+
// Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
- // Get Basis Functions //
- const vector<const vector<Polynomial>*>& fun = getLocalFunctions(element);
- const unsigned int nFun = fun.size();
-
// Get Reference coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
double uvw[3];
@@ -43,6 +39,13 @@ interpolate(const MElement& element,
eelement.getJacobian(uvw[0], uvw[1], uvw[2], invJac);
invJac.invertInPlace();
+ // Get Basis Functions //
+ fullMatrix<double>* fun = basisVector->getFunctions(element,
+ uvw[0],
+ uvw[1],
+ uvw[2]);
+ const unsigned int nFun = fun->size1();
+
// Interpolate (in Reference Place) //
fullVector<double> val(3);
val(0) = 0;
@@ -50,14 +53,13 @@ interpolate(const MElement& element,
val(2) = 0;
for(unsigned int i = 0; i < nFun; i++){
- fullVector<double> vi =
- Polynomial::at(*fun[i], uvw[0], uvw[1], uvw[2]);
-
- vi.scale(coef[i]);
- val.axpy(vi, 1);
+ val(0) += (*fun)(i, 0) * coef[i];
+ val(1) += (*fun)(i, 1) * coef[i];
+ val(2) += (*fun)(i, 2) * coef[i];
}
// Return Interpolated Value //
+ delete fun;
return Mapper::grad(val, invJac);
}
@@ -65,20 +67,24 @@ fullVector<double> FunctionSpaceEdge::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
-
+
// Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
- // Get Basis Functions //
- const vector<const vector<Polynomial>*>& fun = getLocalFunctions(element);
- const unsigned int nFun = fun.size();
-
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
eelement.getJacobian(uvw(0), uvw(1), uvw(2), invJac);
invJac.invertInPlace();
-
+
+ // Get Basis Functions //
+ fullMatrix<double>* fun = basisVector->getFunctions(element,
+ uvw(0),
+ uvw(1),
+ uvw(2));
+ const unsigned int nFun = fun->size1();
+
+
// Interpolate (in Reference Place) //
fullVector<double> val(3);
val(0) = 0;
@@ -86,13 +92,12 @@ interpolateInRefSpace(const MElement& element,
val(2) = 0;
for(unsigned int i = 0; i < nFun; i++){
- fullVector<double> vi =
- Polynomial::at(*fun[i], uvw(0), uvw(1), uvw(2));
-
- vi.scale(coef[i]);
- val.axpy(vi, 1);
+ val(0) += (*fun)(i, 0) * coef[i];
+ val(1) += (*fun)(i, 1) * coef[i];
+ val(2) += (*fun)(i, 2) * coef[i];
}
// Return Interpolated Value //
+ delete fun;
return Mapper::grad(val, invJac);
}
diff --git a/FunctionSpace/FunctionSpaceNode.cpp b/FunctionSpace/FunctionSpaceNode.cpp
index 0adb1198b6..c2dc4227de 100644
--- a/FunctionSpace/FunctionSpaceNode.cpp
+++ b/FunctionSpace/FunctionSpaceNode.cpp
@@ -23,29 +23,32 @@ double FunctionSpaceNode::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
-
+
// Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
-
- // Get Basis Functions //
- const vector<const Polynomial*> fun = getLocalFunctions(element);
- const unsigned int nFun = fun.size();
-
+
// Get Reference coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
double uvw[3];
eelement.xyz2uvw(phys, uvw);
+ // Get Basis Functions //
+ fullMatrix<double>* fun = basisScalar->getFunctions(element,
+ uvw[0],
+ uvw[1],
+ uvw[2]);
+ const unsigned int nFun = fun->size1();
+
// Interpolate (in Reference Place) //
double val = 0;
for(unsigned int i = 0; i < nFun; i++)
- val +=
- fun[i]->at(uvw[0], uvw[1], uvw[2]) * coef[i];
+ val += (*fun)(i, 0) * coef[i];
// Return Interpolated Value //
+ delete fun;
return val;
}
@@ -55,16 +58,19 @@ interpolateInRefSpace(const MElement& element,
const fullVector<double>& uvw) const{
// Get Basis Functions //
- const vector<const Polynomial*> fun = getLocalFunctions(element);
- const unsigned int nFun = fun.size();
+ fullMatrix<double>* fun = basisScalar->getFunctions(element,
+ uvw(0),
+ uvw(1),
+ uvw(2));
+ const unsigned int nFun = fun->size1();
// Interpolate (in Reference Place) //
double val = 0;
for(unsigned int i = 0; i < nFun; i++)
- val +=
- fun[i]->at(uvw(0), uvw(1), uvw(2)) * coef[i];
+ val += (*fun)(i, 0) * coef[i];
// Return Interpolated Value //
+ delete fun;
return val;
}
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index cf828ab28e..e7b173dc8d 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -1,58 +1,7 @@
#include "FunctionSpaceScalar.h"
-using namespace std;
-
FunctionSpaceScalar::FunctionSpaceScalar(void){
- hasGrad = false;
- gradBasis = NULL;
-
- locPreEvaluated = false;
- gradPreEvaluated = false;
- evalLoc = NULL;
- evalGrad = NULL;
}
FunctionSpaceScalar::~FunctionSpaceScalar(void){
- if(hasGrad)
- delete gradBasis;
-
- if(locPreEvaluated)
- delete evalLoc;
-
- if(gradPreEvaluated)
- delete evalGrad;
-}
-
-void FunctionSpaceScalar::
-preEvaluateLocalFunctions(fullMatrix<double>& points){
- // Delete Old Struct (if any) //
- if(locPreEvaluated)
- delete evalLoc;
-
- // New Struct //
- evalLoc = new EvaluatedBasis(*basisScalar, points);
-
- // PreEvaluated //
- locPreEvaluated = true;
-}
-
-void FunctionSpaceScalar::
-preEvaluateGradLocalFunctions(fullMatrix<double>& points){
- // Got Grad Basis ? //
- // --> mutable data
- // --> Just a 'cache memory'
- if(!hasGrad){
- gradBasis = new GradBasis(*basisScalar);
- hasGrad = true;
- }
-
- // Delete Old Struct (if any) //
- if(gradPreEvaluated)
- delete evalGrad;
-
- // New Struct //
- evalGrad = new EvaluatedBasis(*gradBasis, points);
-
- // PreEvaluated //
- gradPreEvaluated = true;
}
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 979252dbb6..19812163ab 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -1,11 +1,8 @@
#ifndef _FUNCTIONSPACESCALAR_H_
#define _FUNCTIONSPACESCALAR_H_
-#include "fullMatrix.h"
#include "Exception.h"
#include "BasisScalar.h"
-#include "GradBasis.h"
-#include "EvaluatedBasis.h"
#include "FunctionSpace.h"
/**
@@ -29,15 +26,6 @@ class FunctionSpaceScalar : public FunctionSpace{
protected:
const BasisScalar* basisScalar;
- mutable bool hasGrad;
- mutable GradBasis* gradBasis;
-
- bool locPreEvaluated;
- bool gradPreEvaluated;
-
- EvaluatedBasis* evalLoc;
- EvaluatedBasis* evalGrad;
-
public:
virtual ~FunctionSpaceScalar(void);
@@ -51,14 +39,8 @@ class FunctionSpaceScalar : public FunctionSpace{
const std::vector<double>& coef,
const fullVector<double>& uvw) const = 0;
- const std::vector<const Polynomial*>&
- getLocalFunctions(const MElement& element) const;
-
- const std::vector<const std::vector<Polynomial>*>&
- getGradLocalFunctions(const MElement& element) const;
-
- void preEvaluateLocalFunctions(fullMatrix<double>& points);
- void preEvaluateGradLocalFunctions(fullMatrix<double>& points);
+ void preEvaluateLocalFunctions(const fullMatrix<double>& point) const;
+ void preEvaluateGradLocalFunctions(const fullMatrix<double>& point) const;
const fullMatrix<double>&
getEvaluatedLocalFunctions(const MElement& element) const;
@@ -177,39 +159,36 @@ class FunctionSpaceScalar : public FunctionSpace{
// Inline Functions //
//////////////////////
-inline const std::vector<const Polynomial*>&
-FunctionSpaceScalar::getLocalFunctions(const MElement& element) const{
- return locBasis(element, *basisScalar);
+inline void FunctionSpaceScalar::
+preEvaluateLocalFunctions(const fullMatrix<double>& point) const{
+ basisScalar->preEvaluateFunctions(point);
}
-inline const std::vector<const std::vector<Polynomial>*>&
-FunctionSpaceScalar::getGradLocalFunctions(const MElement& element) const{
-
- // Got Grad Basis ? //
- // --> mutable data
- // --> Just a 'cache memory'
- if(!hasGrad){
- gradBasis = new GradBasis(*basisScalar);
- hasGrad = true;
- }
-
- return locBasis(element, *gradBasis);
+inline void FunctionSpaceScalar::
+preEvaluateGradLocalFunctions(const fullMatrix<double>& point) const{
+ basisScalar->preEvaluateGradFunctions(point);
}
inline const fullMatrix<double>&
FunctionSpaceScalar::getEvaluatedLocalFunctions(const MElement& element) const{
- if(!locPreEvaluated)
+ try{
+ return basisScalar->getPreEvaluatedFunctions(element);
+ }
+
+ catch(Exception& any){
throw Exception("Local Basis Functions not PreEvaluated");
-
- return locEvalBasis(element, *evalLoc);
+ }
}
inline const fullMatrix<double>&
FunctionSpaceScalar::getEvaluatedGradLocalFunctions(const MElement& element) const{
- if(!gradPreEvaluated)
- throw Exception("Gradients of Local Basis Functions not PreEvaluated");
-
- return locEvalBasis(element, *evalGrad);
+ try{
+ return basisScalar->getPreEvaluatedGradFunctions(element);
+ }
+
+ catch(Exception& any){
+ throw Exception("Gradient of Local Basis Functions not PreEvaluated");
+ }
}
#endif
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index fb4dbbd7e9..e3686be909 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -1,90 +1,9 @@
#include "FunctionSpaceVector.h"
-using namespace std;
-
FunctionSpaceVector::FunctionSpaceVector(void){
- hasCurl = false;
- curlBasis = NULL;
-
- hasDiv = false;
- divBasis = NULL;
-
- locPreEvaluated = false;
- curlPreEvaluated = false;
- divPreEvaluated = false;
- evalLoc = NULL;
- evalCurl = NULL;
- evalDiv = NULL;
}
FunctionSpaceVector::~FunctionSpaceVector(void){
- if(hasCurl)
- delete curlBasis;
-
- if(hasDiv)
- delete divBasis;
-
- if(locPreEvaluated)
- delete evalLoc;
-
- if(curlPreEvaluated)
- delete evalCurl;
-
- if(divPreEvaluated)
- delete evalDiv;
}
-void FunctionSpaceVector::
-preEvaluateLocalFunctions(fullMatrix<double>& points){
- // Delete Old Struct (if any) //
- if(locPreEvaluated)
- delete evalLoc;
-
- // New Struct //
- evalLoc = new EvaluatedBasis(*basisVector, points);
-
- // PreEvaluated //
- locPreEvaluated = true;
-}
-void FunctionSpaceVector::
-preEvaluateCurlLocalFunctions(fullMatrix<double>& points){
- // Got Curl Basis ? //
- // --> mutable data
- // --> Just a 'cache memory'
- if(!hasCurl){
- curlBasis = new CurlBasis(*basisVector);
- hasCurl = true;
- }
-
- // Delete Old Struct (if any) //
- if(curlPreEvaluated)
- delete evalCurl;
-
- // New Struct //
- evalCurl = new EvaluatedBasis(*curlBasis, points);
-
- // PreEvaluated //
- curlPreEvaluated = true;
-}
-
-void FunctionSpaceVector::
-preEvaluateDivLocalFunctions(fullMatrix<double>& points){
- // Got Div Basis ? //
- // --> mutable data
- // --> Just a 'cache memory'
- if(!hasDiv){
- divBasis = new DivBasis(*basisVector);
- hasDiv = true;
- }
-
- // Delete Old Struct (if any) //
- if(divPreEvaluated)
- delete evalDiv;
-
- // New Struct //
- evalDiv = new EvaluatedBasis(*divBasis, points);
-
- // PreEvaluated //
- divPreEvaluated = true;
-}
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 134ff39b3b..2a029b6e7c 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -1,12 +1,8 @@
#ifndef _FUNCTIONSPACEVECTOR_H_
#define _FUNCTIONSPACEVECTOR_H_
-#include "fullMatrix.h"
#include "Exception.h"
#include "BasisVector.h"
-#include "CurlBasis.h"
-#include "DivBasis.h"
-#include "EvaluatedBasis.h"
#include "FunctionSpace.h"
/**
@@ -30,20 +26,6 @@ class FunctionSpaceVector : public FunctionSpace{
protected:
const BasisVector* basisVector;
- mutable bool hasCurl;
- mutable CurlBasis* curlBasis;
-
- mutable bool hasDiv;
- mutable DivBasis* divBasis;
-
- bool locPreEvaluated;
- bool curlPreEvaluated;
- bool divPreEvaluated;
-
- EvaluatedBasis* evalLoc;
- EvaluatedBasis* evalCurl;
- EvaluatedBasis* evalDiv;
-
public:
virtual ~FunctionSpaceVector(void);
@@ -57,19 +39,10 @@ class FunctionSpaceVector : public FunctionSpace{
const std::vector<double>& coef,
const fullVector<double>& uvw) const = 0;
- const std::vector<const std::vector<Polynomial>*>&
- getLocalFunctions(const MElement& element) const;
-
- const std::vector<const std::vector<Polynomial>*>&
- getCurlLocalFunctions(const MElement& element) const;
-
- const std::vector<const Polynomial*>&
- getDivLocalFunctions(const MElement& element) const;
-
- void preEvaluateLocalFunctions(fullMatrix<double>& points);
- void preEvaluateCurlLocalFunctions(fullMatrix<double>& points);
- void preEvaluateDivLocalFunctions(fullMatrix<double>& points);
-
+ void preEvaluateLocalFunctions(const fullMatrix<double>& point) const;
+ void preEvaluateCurlLocalFunctions(const fullMatrix<double>& point) const;
+ void preEvaluateDivLocalFunctions(const fullMatrix<double>& point) const;
+
const fullMatrix<double>&
getEvaluatedLocalFunctions(const MElement& element) const;
@@ -220,61 +193,52 @@ class FunctionSpaceVector : public FunctionSpace{
// Inline Functions //
//////////////////////
-inline const std::vector<const std::vector<Polynomial>*>&
-FunctionSpaceVector::getLocalFunctions(const MElement& element) const{
- return locBasis(element, *basisVector);
+inline void FunctionSpaceVector::
+preEvaluateLocalFunctions(const fullMatrix<double>& point) const{
+ basisVector->preEvaluateFunctions(point);
}
-inline const std::vector<const std::vector<Polynomial>*>&
-FunctionSpaceVector::getCurlLocalFunctions(const MElement& element) const{
-
- // Got Curl Basis ? //
- // --> mutable data
- // --> Just a 'cache memory'
- if(!hasCurl){
- curlBasis = new CurlBasis(*basisVector);
- hasCurl = true;
- }
-
- return locBasis(element, *curlBasis);
+inline void FunctionSpaceVector::
+preEvaluateCurlLocalFunctions(const fullMatrix<double>& point) const{
+ basisVector->preEvaluateCurlFunctions(point);
}
-inline const std::vector<const Polynomial*>&
-FunctionSpaceVector::getDivLocalFunctions(const MElement& element) const{
-
- // Got Div Basis ? //
- // --> mutable data
- // --> Just a 'cache memory'
- if(!hasDiv){
- divBasis = new DivBasis(*basisVector);
- hasDiv = true;
- }
-
- return locBasis(element, *divBasis);
+inline void FunctionSpaceVector::
+preEvaluateDivLocalFunctions(const fullMatrix<double>& point) const{
+ basisVector->preEvaluateDivFunctions(point);
}
inline const fullMatrix<double>&
FunctionSpaceVector::getEvaluatedLocalFunctions(const MElement& element) const{
- if(!locPreEvaluated)
+ try{
+ return basisVector->getPreEvaluatedFunctions(element);
+ }
+
+ catch(Exception& any){
throw Exception("Local Basis Functions not PreEvaluated");
-
- return locEvalBasis(element, *evalLoc);
+ }
}
inline const fullMatrix<double>&
FunctionSpaceVector::getEvaluatedCurlLocalFunctions(const MElement& element) const{
- if(!curlPreEvaluated)
- throw Exception("Curls of Local Basis Functions not PreEvaluated");
-
- return locEvalBasis(element, *evalCurl);
+ try{
+ return basisVector->getPreEvaluatedCurlFunctions(element);
+ }
+
+ catch(Exception& any){
+ throw Exception("Curl of Local Basis Functions not PreEvaluated");
+ }
}
inline const fullMatrix<double>&
FunctionSpaceVector::getEvaluatedDivLocalFunctions(const MElement& element) const{
- if(!divPreEvaluated)
- throw Exception("Divergences of Local Basis Functions not PreEvaluated");
-
- return locEvalBasis(element, *evalDiv);
+ try{
+ return basisVector->getPreEvaluatedFunctions(element);
+ }
+
+ catch(Exception& any){
+ throw Exception("Divergence of Local Basis Functions not PreEvaluated");
+ }
}
#endif
diff --git a/FunctionSpace/GradBasis.cpp b/FunctionSpace/GradBasis.cpp
deleted file mode 100644
index 7507d11b92..0000000000
--- a/FunctionSpace/GradBasis.cpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "GradBasis.h"
-
-using namespace std;
-
-GradBasis::GradBasis(const BasisScalar& other){
- // Reference Space //
- refSpace = &other.getReferenceSpace();
- nRefSpace = other.getReferenceSpace().getNPermutation();
-
- // Set Basis Type //
- order = other.getOrder() - 1;
-
- type = other.getType();
- dim = other.getDim();
-
- nVertex = other.getNVertexBased();
- nEdge = other.getNEdgeBased();
- nFace = other.getNFaceBased();
- nCell = other.getNCellBased();
- nFunction = other.getNFunction();
-
- // Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- for(unsigned int i = 0; i < nFunction; i++)
- (*(*basis)[s])[i] =
- new vector<Polynomial>
- (other.getFunction(s, i).gradient());
-}
-
-GradBasis::~GradBasis(void){
- // Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
-
- delete (*basis)[i];
- }
-
- delete basis;
-}
diff --git a/FunctionSpace/GradBasis.h b/FunctionSpace/GradBasis.h
deleted file mode 100644
index 34b32974c1..0000000000
--- a/FunctionSpace/GradBasis.h
+++ /dev/null
@@ -1,27 +0,0 @@
-#ifndef _GRADBASIS_H_
-#define _GRADBASIS_H_
-
-#include "BasisVector.h"
-#include "BasisScalar.h"
-
-/**
- @class GradBasis
- @brief The Gradient of an other Basis
-
- This class can instantiate a Gradient Basis.
- */
-
-class GradBasis: public BasisVector{
- public:
- //! @param other An Other Basis
- //!
- //! Returns a new Basis, which is the gradient of
- //! the given Basis
- GradBasis(const BasisScalar& other);
-
- //! Deletes this Basis
- //!
- virtual ~GradBasis(void);
-};
-
-#endif
diff --git a/FunctionSpace/HexEdgeBasis.h b/FunctionSpace/HexEdgeBasis.h
index 8990e7c7cd..0daa7609ac 100644
--- a/FunctionSpace/HexEdgeBasis.h
+++ b/FunctionSpace/HexEdgeBasis.h
@@ -1,7 +1,7 @@
#ifndef _HEXEDGEBASIS_H_
#define _HEXEDGEBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class HexEdgeBasis
@@ -15,7 +15,7 @@
Basis for @em high @em order Polynomial%s generation.@n
*/
-class HexEdgeBasis: public BasisVector{
+class HexEdgeBasis: public BasisHierarchicalVector{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/HexNodeBasis.h b/FunctionSpace/HexNodeBasis.h
index 3aced98883..68d4bc1693 100644
--- a/FunctionSpace/HexNodeBasis.h
+++ b/FunctionSpace/HexNodeBasis.h
@@ -1,7 +1,7 @@
#ifndef _HEXNODEBASIS_H_
#define _HEXNODEBASIS_H_
-#include "BasisScalar.h"
+#include "BasisHierarchicalScalar.h"
/**
@class HexNodeBasis
@@ -15,7 +15,7 @@
Basis for @em high @em order Polynomial%s generation.@n
*/
-class HexNodeBasis: public BasisScalar{
+class HexNodeBasis: public BasisHierarchicalScalar{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/LagrangeBasis.cpp b/FunctionSpace/LagrangeBasis.cpp
deleted file mode 100644
index 4336a6ca16..0000000000
--- a/FunctionSpace/LagrangeBasis.cpp
+++ /dev/null
@@ -1,99 +0,0 @@
-#include "Exception.h"
-#include "LagrangeBasis.h"
-
-using namespace std;
-
-LagrangeBasis::LagrangeBasis(void){
-}
-
-LagrangeBasis::~LagrangeBasis(void){
-}
-
-vector<double> LagrangeBasis::
-project(const MElement& element,
- const vector<double>& coef,
- const FunctionSpaceScalar& fSpace){
-
- // Init New Coefs //
- const unsigned int size = point->size1();
- vector<double> newCoef(size);
-
- // Interpolation at Lagrange Points //
- for(unsigned int i = 0; i < size; i++){
- fullVector<double> uvw(3);
- uvw(0) = (*point)(i, 0);
- uvw(1) = (*point)(i, 1);
- uvw(2) = (*point)(i, 2);
-
- newCoef[i] = fSpace.interpolateInRefSpace(element,
- coef,
- uvw);
- }
-
- // Return ;
- return newCoef;
-}
-
-vector<fullVector<double> > LagrangeBasis::
-project(const MElement& element,
- const vector<double>& coef,
- const FunctionSpaceVector& fSpace){
-
- // Init New Coefs //
- const unsigned int size = point->size1();
- vector<fullVector<double> > newCoef(size);
-
- // Interpolation at Lagrange Points //
- for(unsigned int i = 0; i < size; i++){
- fullVector<double> uvw(3);
- uvw(0) = (*point)(i, 0);
- uvw(1) = (*point)(i, 1);
- uvw(2) = (*point)(i, 2);
-
- newCoef[i] = fSpace.interpolateInRefSpace(element,
- coef,
- uvw);
- }
-
- // Return ;
- return newCoef;
-}
-
-unsigned int** LagrangeBasis::triEdgeOrder(int order){
- // Check //
- if(order <= 1)
- return NULL;
-
- // Number of Edge Based Basis //
- const unsigned int orderMinus = order - 1;
- const unsigned int size = 3 * orderMinus;
-
- // Alloc //
- unsigned int i;
- unsigned int** edgeOrder = new unsigned int*[2];
- edgeOrder[0] = new unsigned int[size];
- edgeOrder[1] = new unsigned int[size];
-
- // Direct Numeration //
- i = 0;
-
- for(unsigned int j = 0; j < orderMinus; j++){
- for(unsigned int k = 0; k < size; k += orderMinus){
- edgeOrder[0][i] = k + j;
- i++;
- }
- }
-
- // Invert Numeration //
- i = 0;
-
- for(unsigned int j = 0; j < orderMinus; j++){
- for(unsigned int k = 0; k < size; k += orderMinus){
- edgeOrder[1][i] = k + orderMinus - 1 - j;
- i++;
- }
- }
-
- // Retrun //
- return edgeOrder;
-}
diff --git a/FunctionSpace/LagrangeBasis.h b/FunctionSpace/LagrangeBasis.h
deleted file mode 100644
index 940040348b..0000000000
--- a/FunctionSpace/LagrangeBasis.h
+++ /dev/null
@@ -1,97 +0,0 @@
-#ifndef _LAGRANGEBASIS_H_
-#define _LAGRANGEBASIS_H_
-
-#include <vector>
-
-#include "polynomialBasis.h"
-#include "BasisScalar.h"
-
-#include "MElement.h"
-#include "FunctionSpaceScalar.h"
-#include "FunctionSpaceVector.h"
-
-/**
- @interface LagrangeBasis
- @brief Interface for Lagrange Basis
-
- This is an interface for Lagrange Basis.@n
-
- These Scalar Basis allow a @em Coefficient Matrix,
- and a Monomial Matrix, to be consulted.@n
-
- A vector from an Other Basis (set of Functions)
- can also be projected into a Lagrange Basis.@n
-
- @todo
- Add a method to erase polynomialBasis in polynomialBasis@n
- Add a method to get lagrange Point in polynomialBasis
- */
-
-class LagrangeBasis: public BasisScalar{
- protected:
- fullMatrix<double>* point;
- const polynomialBasis* l;
-
- public:
- //! Deletes this Basis
- //!
- virtual ~LagrangeBasis(void);
-
- //! @return Returns the Coefficient Matrix
- const fullMatrix<double>& getCoefficient(void) const;
-
- //! @return Returns the Monomial Matrix
- const fullMatrix<double>& getMonomial(void) const;
-
- //! @param element A MElement
- //! @param coef A vector of coefficient associated
- //! to the given Element
- //! @param fSpace The (scalar) Function Space
- //! of the given Coefficients
- //! @return Projects the given Coefficients in this LagrangeBasis
- std::vector<double> project(const MElement& element,
- const std::vector<double>& coef,
- const FunctionSpaceScalar& fSpace);
-
- //! @param element A MElement
- //! @param coef A vector of coefficient associated
- //! to the given Element
- //! @param fSpace The (vectorial) Function Space
- //! of the given Coefficients
- //! @return Projects the given Coefficients in this LagrangeBasis
- std::vector<fullVector<double> >
- project(const MElement& element,
- const std::vector<double>& coef,
- const FunctionSpaceVector& fSpace);
- protected:
- //! Returns a new LagrangeBasis
- //!
- LagrangeBasis(void);
-
- //! @param order The ordre of the requested LagrangeBasis
- //! @return Returns a 2D table
- //! (of size 2 x (3 * (@c order - 1)))
- //! with the Edge Triangle Inverted Closure for Lagrange
- //! @note The first dimension of the returned table
- //! is the direct closure, and the second dimension
- //! the invert closure
- //! @warning The returned table @em must be @em deleted
- static unsigned int** triEdgeOrder(int order);
-};
-
-
-//////////////////////
-// Inline Functions //
-//////////////////////
-
-inline const fullMatrix<double>& LagrangeBasis::
-getCoefficient(void) const{
- return l->coefficients;
-}
-
-inline const fullMatrix<double>& LagrangeBasis::
-getMonomial(void) const{
- return l->monomials;
-}
-
-#endif
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
index d822a81b84..408061a8c1 100644
--- a/FunctionSpace/LineEdgeBasis.cpp
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -47,21 +47,21 @@ LineEdgeBasis::LineEdgeBasis(unsigned int order){
Legendre::integrated(intLegendre, orderPlus);
// Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
+ basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- (*(*basis)[0])[0] = new vector<Polynomial>(first);
- (*(*basis)[1])[0] = new vector<Polynomial>(second);
+ basis[0][0] = new vector<Polynomial>(first);
+ basis[1][0] = new vector<Polynomial>(second);
// Edge Based (High Order) //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = 1;
for(unsigned int l = 1; l < orderPlus; l++){
- (*(*basis)[s])[i] =
+ basis[s][i] =
new vector<Polynomial>((intLegendre[l].compose
(x[(*(*edgeV[s])[0])[0]])).gradient());
@@ -80,10 +80,10 @@ LineEdgeBasis::~LineEdgeBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/LineEdgeBasis.h b/FunctionSpace/LineEdgeBasis.h
index 7986b2a7df..50519874e3 100644
--- a/FunctionSpace/LineEdgeBasis.h
+++ b/FunctionSpace/LineEdgeBasis.h
@@ -1,7 +1,7 @@
#ifndef _LINEEDGEBASIS_H_
#define _LINEEDGEBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class LineEdgeBasis
@@ -19,7 +19,7 @@
It also uses the following mapping: @f$x = \frac{u + 1}{2}@f$.
*/
-class LineEdgeBasis: public BasisVector{
+class LineEdgeBasis: public BasisHierarchicalVector{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
index 3dd3545bb9..7f3d9fd8aa 100644
--- a/FunctionSpace/LineNedelecBasis.cpp
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -33,14 +33,14 @@ LineNedelecBasis::LineNedelecBasis(void){
second[2] = Polynomial( 0 , 0, 0, 0);
// Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
+ basis[s] = new vector<Polynomial>*[nFunction];
// Nedelec //
- (*(*basis)[0])[0] = new vector<Polynomial>(first);
- (*(*basis)[1])[0] = new vector<Polynomial>(second);
+ basis[0][0] = new vector<Polynomial>(first);
+ basis[1][0] = new vector<Polynomial>(second);
}
LineNedelecBasis::~LineNedelecBasis(void){
@@ -50,10 +50,10 @@ LineNedelecBasis::~LineNedelecBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/LineNedelecBasis.h b/FunctionSpace/LineNedelecBasis.h
index 8579b1a470..a23ee315bc 100644
--- a/FunctionSpace/LineNedelecBasis.h
+++ b/FunctionSpace/LineNedelecBasis.h
@@ -1,7 +1,7 @@
#ifndef _LINENEDELECBASIS_H_
#define _LINENEDELECBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class LineNedelecBasis
@@ -11,7 +11,7 @@
for Lines.@n
*/
-class LineNedelecBasis: public BasisVector{
+class LineNedelecBasis: public BasisHierarchicalVector{
public:
//! Returns a new Nedelec Basis for Lines
//!
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index 0c80278d59..6e9c8933f0 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -36,18 +36,18 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
Legendre::integrated(intLegendre, order);
// Basis //
- basis = new vector<vector<const Polynomial*>*>(nRefSpace);
+ basis = new Polynomial**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const Polynomial*>(nFunction);
+ basis[s] = new Polynomial*[nFunction];
// Vertex Based (Lagrange) //
for(unsigned int s = 0; s < nRefSpace; s++){
- (*(*basis)[s])[0] =
+ basis[s][0] =
new Polynomial(Polynomial(0.5, 0, 0, 0) -
Polynomial(0.5, 1, 0, 0));
- (*(*basis)[s])[1] =
+ basis[s][1] =
new Polynomial(Polynomial(0.5, 0, 0, 0) +
Polynomial(0.5, 1, 0, 0));
}
@@ -57,7 +57,7 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
unsigned int i = nVertex;
for(unsigned int l = 1; l < order; l++){
- (*(*basis)[s])[i] =
+ basis[s][i] =
new Polynomial(intLegendre[l].compose(x[(*(*edgeV[s])[0])[0]]));
i++;
@@ -75,10 +75,10 @@ LineNodeBasis::~LineNodeBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/LineNodeBasis.h b/FunctionSpace/LineNodeBasis.h
index 9bc6bbea9f..cc555d92d8 100644
--- a/FunctionSpace/LineNodeBasis.h
+++ b/FunctionSpace/LineNodeBasis.h
@@ -1,7 +1,7 @@
#ifndef _LINENODEBASIS_H_
#define _LINENODEBASIS_H_
-#include "BasisScalar.h"
+#include "BasisHierarchicalScalar.h"
/**
@class LineNodeBasis
@@ -19,7 +19,7 @@
It also uses the following mapping: @f$x = \frac{u + 1}{2}@f$.
*/
-class LineNodeBasis: public BasisScalar{
+class LineNodeBasis: public BasisHierarchicalScalar{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/QuadEdgeBasis.h b/FunctionSpace/QuadEdgeBasis.h
index 40db8df987..b5a624b5eb 100644
--- a/FunctionSpace/QuadEdgeBasis.h
+++ b/FunctionSpace/QuadEdgeBasis.h
@@ -1,7 +1,7 @@
#ifndef _QUADEDGEBASIS_H_
#define _QUADEDGEBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class QuadEdgeBasis
@@ -19,7 +19,7 @@
@li @f$y = \frac{v + 1}{2}@f$
*/
-class QuadEdgeBasis: public BasisVector{
+class QuadEdgeBasis: public BasisHierarchicalVector{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/QuadNodeBasis.h b/FunctionSpace/QuadNodeBasis.h
index 80f4ff0432..16c2eed893 100644
--- a/FunctionSpace/QuadNodeBasis.h
+++ b/FunctionSpace/QuadNodeBasis.h
@@ -1,7 +1,7 @@
#ifndef _QUADNODEBASIS_H_
#define _QUADNODEBASIS_H_
-#include "BasisScalar.h"
+#include "BasisHierarchicalScalar.h"
/**
@class QuadNodeBasis
@@ -19,7 +19,7 @@
@li @f$y = \frac{v + 1}{2}@f$
*/
-class QuadNodeBasis: public BasisScalar{
+class QuadNodeBasis: public BasisHierarchicalScalar{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 6d51adbfd9..e29bcb5599 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -58,196 +58,142 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
+ basis[s] = new vector<Polynomial>*[nFunction];
- // Edge Based (Nedelec) //
+ // Edge Based //
for(unsigned int s = 0; s < nRefSpace; s++){
- for(int e = 0; e < 6; e++){
- vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
- vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
-
- tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
-
- tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
-
- tmp1[0].sub(tmp2[0]);
- tmp1[1].sub(tmp2[1]);
- tmp1[2].sub(tmp2[2]);
+ unsigned int i = 0;
- (*(*basis)[s])[e] = new vector<Polynomial>(tmp1);
- }
- }
-
- // Edge Based (High Order) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = 6;
+ for(int e = 0; e < 6; e++){
+ for(int l = 0; l < orderPlus; l++){
+ // Nedelec
+ if(l == 0){
+ vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+
+ tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+
+ tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+
+ tmp1[0].sub(tmp2[0]);
+ tmp1[1].sub(tmp2[1]);
+ tmp1[2].sub(tmp2[2]);
+
+ basis[s][i] = new vector<Polynomial>(tmp1);
+ }
- for(int l = 1; l < orderPlus; l++){
- for(int e = 0; e < 6; e++){
- (*(*basis)[s])[i] =
- new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
- ,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]])).gradient());
+ // High Order
+ else{
+ basis[s][i] =
+ new vector<Polynomial>
+ ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
+ lagrange[(*(*edgeV[s])[e])[1]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[1]])).gradient());
+ }
i++;
}
}
}
- // Face Based (Preliminary) //
- unsigned int faceOffset = 0;
-
- // Sum
- Polynomial** sum;
- sum = new Polynomial*[nRefSpace];
-
- for(unsigned int s = 0; s < nRefSpace; s++){
- sum[s] = new Polynomial[4];
-
- for(unsigned int f = 0; f < 4; f++)
- sum[s][f] =
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]] +
- lagrange[(*(*faceV[s])[f])[2]];
- }
-
- // Polynomial u
- Polynomial*** u;
- u = new Polynomial**[nRefSpace];
-
- for(unsigned int s = 0; s < nRefSpace; s++){
- u[s] = new Polynomial*[orderPlus];
-
- for(int l = 0; l < orderPlus; l++){
- u[s][l] = new Polynomial[4];
-
- for(unsigned int f = 0; f < 4; f++)
- u[s][l][f] =
- intLegendre[l].compose(lagrange[(*(*faceV[s])[f])[0]] -
- lagrange[(*(*faceV[s])[f])[1]]
- ,
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]]);
- }
- }
-
- // Polynomial v
- Polynomial*** v;
- v = new Polynomial**[nRefSpace];
-
- for(unsigned int s = 0; s < nRefSpace; s++){
- v[s] = new Polynomial*[orderMinus];
-
- for(int l = 0; l < orderMinus; l++){
- v[s][l] = new Polynomial[4];
-
- for(unsigned int f = 0; f < 4; f++)
- v[s][l][f] =
- lagrange[(*(*faceV[s])[f])[2]] *
- sclLegendre[l].compose
- (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum[s][f], sum[s][f]);
- }
- }
-
- // Face Based (Type 1) //
+ // Face Based //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
-
- for(unsigned int l1 = 1; l1 < order; l1++){
- for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
- for(int f = 0; f < 4; f++){
+
+ for(int f = 0; f < 4; f++){
+ for(unsigned int l1 = 1; l1 < order; l1++){
+ for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
+ // Preliminary Type 1
+ Polynomial sum =
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]] +
+ lagrange[(*(*faceV[s])[f])[2]];
- (*(*basis)[s])[i] =
- new vector<Polynomial>((u[s][l1][f] * v[s][l2][f]).gradient());
+ Polynomial u =
+ intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
+ lagrange[(*(*faceV[s])[f])[1]]
+ ,
+ lagrange[(*(*faceV[s])[f])[0]] +
+ lagrange[(*(*faceV[s])[f])[1]]);
+ Polynomial v =
+ lagrange[(*(*faceV[s])[f])[2]] *
+ sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
- i++;
+ // Preliminary Type 2
+ vector<Polynomial> gradU = u.gradient();
+ vector<Polynomial> gradV = v.gradient();
- if(s == 0)
- faceOffset++;
- }
- }
- }
- }
-
- // Face Based (Type 2) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nEdge + faceOffset;
-
- for(unsigned int l1 = 1; l1 < order; l1++){
- for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
- for(int f = 0; f < 4; f++){
- vector<Polynomial> gradU = u[s][l1][f].gradient();
- vector<Polynomial> gradV = v[s][l2][f].gradient();
-
vector<Polynomial> vGradU(gradU);
- vGradU[0].mul(v[s][l2][f]);
- vGradU[1].mul(v[s][l2][f]);
- vGradU[2].mul(v[s][l2][f]);
-
+ vGradU[0].mul(v);
+ vGradU[1].mul(v);
+ vGradU[2].mul(v);
+
vector<Polynomial> uGradV(gradV);
- uGradV[0].mul(u[s][l1][f]);
- uGradV[1].mul(u[s][l1][f]);
- uGradV[2].mul(u[s][l1][f]);
-
+ uGradV[0].mul(u);
+ uGradV[1].mul(u);
+ uGradV[2].mul(u);
+
vector<Polynomial> subGradUV(vGradU);
subGradUV[0].sub(uGradV[0]);
subGradUV[1].sub(uGradV[1]);
subGradUV[2].sub(uGradV[2]);
-
- (*(*basis)[s])[i] =
+
+ // Preliminary Type 3
+ vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
+ vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
+
+ vector<Polynomial> l2GradL1(gradL1);
+ l2GradL1[0].mul(lagrange[(*(*faceV[s])[f])[1]]);
+ l2GradL1[1].mul(lagrange[(*(*faceV[s])[f])[1]]);
+ l2GradL1[2].mul(lagrange[(*(*faceV[s])[f])[1]]);
+
+ vector<Polynomial> l1GradL2(gradL2);
+ l1GradL2[0].mul(lagrange[(*(*faceV[s])[f])[0]]);
+ l1GradL2[1].mul(lagrange[(*(*faceV[s])[f])[0]]);
+ l1GradL2[2].mul(lagrange[(*(*faceV[s])[f])[0]]);
+
+ vector<Polynomial> subGradL1L2V(l2GradL1);
+ subGradL1L2V[0].sub(l1GradL2[0]);
+ subGradL1L2V[1].sub(l1GradL2[1]);
+ subGradL1L2V[2].sub(l1GradL2[2]);
+
+ subGradL1L2V[0].mul(v);
+ subGradL1L2V[1].mul(v);
+ subGradL1L2V[2].mul(v);
+
+
+ // Type 1
+ basis[s][i] =
+ new vector<Polynomial>((u * v).gradient());
+
+ i++;
+
+ // Type 2
+ basis[s][i] =
new vector<Polynomial>(subGradUV);
-
+
i++;
- }
- }
- }
- }
-
- // Face Based (Type 3) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nEdge + faceOffset * 2;
-
- for(int l2 = 0; l2 < orderMinus; l2++){
- for(int f = 0; f < 4; f++){
- vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
- vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
-
- vector<Polynomial> l2GradL1(gradL1);
- l2GradL1[0].mul(lagrange[(*(*faceV[s])[f])[1]]);
- l2GradL1[1].mul(lagrange[(*(*faceV[s])[f])[1]]);
- l2GradL1[2].mul(lagrange[(*(*faceV[s])[f])[1]]);
-
- vector<Polynomial> l1GradL2(gradL2);
- l1GradL2[0].mul(lagrange[(*(*faceV[s])[f])[0]]);
- l1GradL2[1].mul(lagrange[(*(*faceV[s])[f])[0]]);
- l1GradL2[2].mul(lagrange[(*(*faceV[s])[f])[0]]);
-
- vector<Polynomial> subGradL1L2V(l2GradL1);
- subGradL1L2V[0].sub(l1GradL2[0]);
- subGradL1L2V[1].sub(l1GradL2[1]);
- subGradL1L2V[2].sub(l1GradL2[2]);
-
- subGradL1L2V[0].mul(v[s][l2][f]);
- subGradL1L2V[1].mul(v[s][l2][f]);
- subGradL1L2V[2].mul(v[s][l2][f]);
- (*(*basis)[s])[i] =
- new vector<Polynomial>(subGradL1L2V);
-
- i++;
+ // Type 3
+ if(l1 == 1){
+ basis[s][i] =
+ new vector<Polynomial>(subGradL1L2V);
+
+ i++;
+ }
+ }
}
}
- }
-
+ }
+
// Cell Based //
const Polynomial one(1, 0, 0, 0);
@@ -347,20 +293,20 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Type 1
- (*(*basis)[s])[i] = new vector<Polynomial>((u * v * w).gradient());
+ basis[s][i] = new vector<Polynomial>((u * v * w).gradient());
i++;
// Type 2 -- Part 1
- (*(*basis)[s])[i] = new vector<Polynomial>(term1);
+ basis[s][i] = new vector<Polynomial>(term1);
i++;
// Type 2 -- Part 2
- (*(*basis)[s])[i] = new vector<Polynomial>(term2);
+ basis[s][i] = new vector<Polynomial>(term2);
i++;
// Type 3
if(l1 == 1){
- (*(*basis)[s])[i] = new vector<Polynomial>(subGradL1L2VW);
+ basis[s][i] = new vector<Polynomial>(subGradL1L2VW);
i++;
}
}
@@ -373,32 +319,6 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
-
- // Sum
- for(unsigned int s = 0; s < nRefSpace; s++)
- delete[] sum[s];
-
- delete[] sum;
-
- // Polynomial u
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(int l = 0; l < orderPlus; l++)
- delete[] u[s][l];
-
- delete[] u[s];
- }
-
- delete[] u;
-
- // Polynomial v
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(int l = 0; l < orderMinus; l++)
- delete[] v[s][l];
-
- delete[] v[s];
- }
-
- delete[] v;
}
TetEdgeBasis::~TetEdgeBasis(void){
@@ -408,11 +328,11 @@ TetEdgeBasis::~TetEdgeBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/TetEdgeBasis.h b/FunctionSpace/TetEdgeBasis.h
index d5e5444b93..e0f58df669 100644
--- a/FunctionSpace/TetEdgeBasis.h
+++ b/FunctionSpace/TetEdgeBasis.h
@@ -1,7 +1,7 @@
#ifndef _TETEDGEBASIS_H_
#define _TETEDGEBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class TetEdgeBasis
@@ -15,7 +15,7 @@
Basis for @em high @em order Polynomial%s generation.@n
*/
-class TetEdgeBasis: public BasisVector{
+class TetEdgeBasis: public BasisHierarchicalVector{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 118c0bd5c2..4f36f7a8d2 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -58,26 +58,26 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
// Basis //
- basis = new vector<vector<const Polynomial*>*>(nRefSpace);
+ basis = new Polynomial**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const Polynomial*>(nFunction);
+ basis[s] = new Polynomial*[nFunction];
// Vertex Based //
for(unsigned int s = 0; s < nRefSpace; s++){
- (*(*basis)[s])[0] = new Polynomial(lagrange[0]);
- (*(*basis)[s])[1] = new Polynomial(lagrange[1]);
- (*(*basis)[s])[2] = new Polynomial(lagrange[2]);
- (*(*basis)[s])[3] = new Polynomial(lagrange[3]);
+ basis[s][0] = new Polynomial(lagrange[0]);
+ basis[s][1] = new Polynomial(lagrange[1]);
+ basis[s][2] = new Polynomial(lagrange[2]);
+ basis[s][3] = new Polynomial(lagrange[3]);
}
// Edge Based //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex;
- for(unsigned int l = 1; l < order; l++){
- for(int e = 0; e < 6; e++){
- (*(*basis)[s])[i] =
+ for(int e = 0; e < 6; e++){
+ for(unsigned int l = 1; l < order; l++){
+ basis[s][i] =
new Polynomial(intLegendre[l].compose
(lagrange[(*(*edgeV[s])[e])[0]] -
lagrange[(*(*edgeV[s])[e])[1]]
@@ -93,15 +93,15 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex + nEdge;
- for(int l1 = 1; l1 < orderMinus; l1++){
- for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
- for(int f = 0; f < 4; f++){
+ for(int f = 0; f < 4; f++){
+ for(int l1 = 1; l1 < orderMinus; l1++){
+ for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
Polynomial sum =
lagrange[(*(*faceV[s])[f])[0]] +
lagrange[(*(*faceV[s])[f])[1]] +
lagrange[(*(*faceV[s])[f])[2]];
- (*(*basis)[s])[i] =
+ basis[s][i] =
new Polynomial(intLegendre[l1].compose
(lagrange[(*(*faceV[s])[f])[0]] -
lagrange[(*(*faceV[s])[f])[1]]
@@ -135,7 +135,7 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
- (*(*basis)[s])[i] =
+ basis[s][i] =
new Polynomial(intLegendre[l1].compose(sub, add) *
lagrange[2] *
sclLegendre[l2].compose(twoThreeOneMinusFour,
@@ -161,10 +161,10 @@ TetNodeBasis::~TetNodeBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/TetNodeBasis.h b/FunctionSpace/TetNodeBasis.h
index 9ade551271..efa5043dae 100644
--- a/FunctionSpace/TetNodeBasis.h
+++ b/FunctionSpace/TetNodeBasis.h
@@ -1,7 +1,7 @@
#ifndef _TETNODEBASIS_H_
#define _TETNODEBASIS_H_
-#include "BasisScalar.h"
+#include "BasisHierarchicalScalar.h"
/**
@class TetNodeBasis
@@ -15,7 +15,7 @@
Basis for @em high @em order Polynomial%s generation.@n
*/
-class TetNodeBasis: public BasisScalar{
+class TetNodeBasis: public BasisHierarchicalScalar{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 99b56254ae..73e4f83ab3 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -50,47 +50,48 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
};
// Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
+ basis[s] = new vector<Polynomial>*[nFunction];
- // Edge Based (Nedelec) //
+ // Edge Based //
for(unsigned int s = 0; s < nRefSpace; s++){
- for(int e = 0; e < 3; e++){
- vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
- vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
-
- tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
- tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
-
-
- tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
- tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
-
- tmp2[0].sub(tmp1[0]);
- tmp2[1].sub(tmp1[1]);
- tmp2[2].sub(tmp1[2]);
+ unsigned int i = 0;
- (*(*basis)[s])[e] = new vector<Polynomial>(tmp2);
- }
- }
-
- // Edge Based (High Order) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = 3;
-
- for(int l = 1; l < orderPlus; l++){
- for(int e = 0; e < 3; e++){
- (*(*basis)[s])[i] =
- new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
- ,
- lagrange[(*(*edgeV[s])[e])[1]] +
- lagrange[(*(*edgeV[s])[e])[0]])).gradient());
+ for(int e = 0; e < 3; e++){
+ for(int l = 0; l < orderPlus; l++){
+ // Nedelec
+ if(l == 0){
+ vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+
+ tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+
+
+ tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+
+ basis[s][i] = new vector<Polynomial>(tmp2);
+ }
+
+ // High Order
+ else{
+ basis[s][i] =
+ new vector<Polynomial>
+ ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
+ lagrange[(*(*edgeV[s])[e])[1]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[1]] +
+ lagrange[(*(*edgeV[s])[e])[0]])).gradient());
+ }
i++;
}
}
@@ -147,7 +148,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
tmp2[1].add(tmp1[1]);
tmp2[2].add(tmp1[2]);
- (*(*basis)[s])[i] = new vector<Polynomial>(tmp2);
+ basis[s][i] = new vector<Polynomial>(tmp2);
i++;
}
@@ -171,7 +172,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
tmp2[1].sub(tmp1[1]);
tmp2[2].sub(tmp1[2]);
- (*(*basis)[s])[i] = new vector<Polynomial>(tmp2);
+ basis[s][i] = new vector<Polynomial>(tmp2);
i++;
}
@@ -185,7 +186,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
subGradL1L2V[1].mul(v[l]);
subGradL1L2V[2].mul(v[l]);
- (*(*basis)[s])[i] = new vector<Polynomial>(subGradL1L2V);
+ basis[s][i] = new vector<Polynomial>(subGradL1L2V);
i++;
}
@@ -205,10 +206,10 @@ TriEdgeBasis::~TriEdgeBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/TriEdgeBasis.h b/FunctionSpace/TriEdgeBasis.h
index 4f9a1c40b7..2dbcb9bbe7 100644
--- a/FunctionSpace/TriEdgeBasis.h
+++ b/FunctionSpace/TriEdgeBasis.h
@@ -1,7 +1,7 @@
#ifndef _TRIEDGEBASIS_H_
#define _TRIEDGEBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class TriEdgeBasis
@@ -15,7 +15,7 @@
Basis for @em high @em order Polynomial%s generation.@n
*/
-class TriEdgeBasis: public BasisVector{
+class TriEdgeBasis: public BasisHierarchicalVector{
public:
//! @param order The order of the Basis
//!
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index cfd97bb7ff..40192124c8 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,72 +1,7 @@
-#include "BasisFactory.h"
#include "Exception.h"
-#include "TriReferenceSpace.h"
#include "TriLagrangeBasis.h"
-using namespace std;
-
TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
- throw Exception("TriLagrangeBasis: Removed");
- /*
- // Call polynomialBasis procedure //
- unsigned int tag;
-
- switch(order){
- case 1:
- tag = MSH_TRI_3;
- break;
-
- case 2:
- tag = MSH_TRI_6;
- break;
-
- case 3:
- tag = MSH_TRI_10;
- break;
-
- case 4:
- tag = MSH_TRI_15;
- break;
-
- case 5:
- tag = MSH_TRI_21;
- break;
-
- case 6:
- tag = MSH_TRI_28;
- break;
-
- case 7:
- tag = MSH_TRI_36;
- break;
-
- case 8:
- tag = MSH_TRI_45;
- break;
-
- case 9:
- tag = MSH_TRI_55;
- break;
-
- case 10:
- tag = MSH_TRI_66;
- break;
-
- default:
- throw Exception
- ("Can't instanciate an order %d Lagrangian Basis for a Triangle",
- order);
-
- break;
- }
-
- l = (polynomialBasis*)BasisFactory::create(tag);
- point = new fullMatrix<double>(triPoint(order));
-
- // Reference Space //
- refSpace = new TriReferenceSpace;
- nRefSpace = refSpace->getNReferenceSpace();
-
// Set Basis Type //
this->order = order;
@@ -79,101 +14,35 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
nCell = (order - 1) * (order - 2) / 2;
nFunction = nVertex + nEdge + nFace + nCell;
- // Alloc Some Stuff //
- const unsigned int nMonomial = l->monomials.size1();
- unsigned int** edgeOrder;
- Polynomial* pEdgeClosureLess = new Polynomial[nEdge];
-
- // Basis //
- basis = new vector<vector<const Polynomial*>*>(nRefSpace);
-
- for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const Polynomial*>(nFunction);
-
- // Vertex Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int i = 0; i < nVertex; i++){
- Polynomial p = Polynomial(0, 0, 0, 0);
-
- for(unsigned int j = 0; j < nMonomial; j++)
- p = p + Polynomial(l->coefficients(i, j), // Coef
- l->monomials(j, 0), // powerX
- l->monomials(j, 1), // powerY
- 0); // powerZ
-
- (*(*basis)[s])[i] = new Polynomial(p);
- }
- }
-
- // Edge Based //
- // Without Closure
- for(unsigned int i = 0; i < nEdge; i++){
- unsigned int ci = i + nVertex;
- pEdgeClosureLess[i] = Polynomial(0, 0, 0, 0);
-
- for(unsigned int j = 0; j < nMonomial; j++)
- pEdgeClosureLess[i] =
- pEdgeClosureLess[i] +
- Polynomial(l->coefficients(ci, j), // Coef
- l->monomials(j, 0), // powerX
- l->monomials(j, 1), // powerY
- 0); // powerZ
- }
-
- // With Closure
- edgeOrder = triEdgeOrder(order); // Closure Ordering
-
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int i = nVertex; i < nEdge; i++){
- (*(*basis)[s])[i] =
- new Polynomial(pEdgeClosureLess[edgeOrder[s][i]]);
- }
- }
-
- // Cell Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int i = nVertex + nEdge; i < nCell; i++){
- Polynomial p = Polynomial(0, 0, 0, 0);
-
- for(unsigned int j = 0; j < nMonomial; j++)
- p = p + Polynomial(l->coefficients(i, j), // Coef
- l->monomials(j, 0), // powerX
- l->monomials(j, 1), // powerY
- 0); // powerZ
-
- (*(*basis)[s])[i] = new Polynomial(p);
- }
- }
-
- // Delete Temporary Space //
- delete[] pEdgeClosureLess;
-
- if(edgeOrder){
- delete[] edgeOrder[0];
- delete[] edgeOrder[1];
- delete[] edgeOrder;
- }
- */
+ // Init polynomialBasis //
+ basis = new polynomialBasis(getTag(order));
}
TriLagrangeBasis::~TriLagrangeBasis(void){
- // Delete gmsh polynomial Basis //
- // --> no method to do so :-(
- // --> erased ??
- // ReferenceSpace //
- delete refSpace;
-
- // Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
-
- delete (*basis)[i];
- }
-
delete basis;
}
+unsigned int TriLagrangeBasis::getTag(unsigned int order){
+ switch(order){
+ case 1: return MSH_TRI_3;
+ case 2: return MSH_TRI_6;
+ case 3: return MSH_TRI_10;
+ case 4: return MSH_TRI_15;
+ case 5: return MSH_TRI_21;
+ case 6: return MSH_TRI_28;
+ case 7: return MSH_TRI_36;
+ case 8: return MSH_TRI_45;
+ case 9: return MSH_TRI_55;
+ case 10: return MSH_TRI_66;
+
+ default:
+ throw Exception
+ ("Can't instanciate an order %d Lagrangian Basis for a Triangle",
+ order);
+ }
+}
+
+/*
fullMatrix<double> TriLagrangeBasis::
triPoint(unsigned int order){
unsigned int nbPoints = (order + 1) * (order + 2) / 2;
@@ -242,3 +111,4 @@ triPoint(unsigned int order){
return point;
}
+*/
diff --git a/FunctionSpace/TriLagrangeBasis.h b/FunctionSpace/TriLagrangeBasis.h
index 39b1b2f976..a008a69fb2 100644
--- a/FunctionSpace/TriLagrangeBasis.h
+++ b/FunctionSpace/TriLagrangeBasis.h
@@ -1,7 +1,7 @@
#ifndef _TRILAGRANGEBASIS_H_
#define _TRILAGRANGEBASIS_H_
-#include "LagrangeBasis.h"
+#include "BasisLagrange.h"
/**
@class TriLagrangeBasis
@@ -11,13 +11,10 @@
for a Triangle and for a given Order.@n
It uses
- <a href="http://geuz.org/gmsh/">gmsh</a> Basis.@n
-
- @todo
- Add method to erase polynomialBasis in polynomialBasis
+ <a href="http://geuz.org/gmsh/">gmsh</a> Basis.
*/
-class TriLagrangeBasis: public LagrangeBasis{
+class TriLagrangeBasis: public BasisLagrange{
public:
//! @param order A natural number
//!
@@ -30,6 +27,11 @@ class TriLagrangeBasis: public LagrangeBasis{
virtual ~TriLagrangeBasis(void);
private:
+ //! @param order A natural number
+ //! @return Returns the @em tag of a @em Triangle of
+ //! the given order
+ static unsigned int getTag(unsigned int order);
+
//! @param order A natural number
//! @return Returns Lagrangian Points on a Triangle
//! for the given Order
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index d88d0be855..6ce5cd4f4b 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -36,10 +36,10 @@ TriNedelecBasis::TriNedelecBasis(void){
};
// Basis //
- basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
+ basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
for(unsigned int s = 0; s < nRefSpace; s++){
@@ -60,7 +60,7 @@ TriNedelecBasis::TriNedelecBasis(void){
tmp2[1].sub(tmp1[1]);
tmp2[2].sub(tmp1[2]);
- (*(*basis)[s])[e] = new vector<Polynomial>(tmp2);
+ basis[s][e] = new vector<Polynomial>(tmp2);
}
}
}
@@ -72,10 +72,10 @@ TriNedelecBasis::~TriNedelecBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/TriNedelecBasis.h b/FunctionSpace/TriNedelecBasis.h
index f23a77e64d..47a4d480c6 100644
--- a/FunctionSpace/TriNedelecBasis.h
+++ b/FunctionSpace/TriNedelecBasis.h
@@ -1,7 +1,7 @@
#ifndef _TRINEDELECBASIS_H_
#define _TRINEDELECBASIS_H_
-#include "BasisVector.h"
+#include "BasisHierarchicalVector.h"
/**
@class TriNedelecBasis
@@ -11,7 +11,7 @@
for Triangles.@n
*/
-class TriNedelecBasis: public BasisVector{
+class TriNedelecBasis: public BasisHierarchicalVector{
public:
//! Returns a new Nedelec Basis for Triangles
//!
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index d4556c5575..55d769c73d 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -47,25 +47,25 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
};
// Basis //
- basis = new vector<vector<const Polynomial*>*>(nRefSpace);
+ basis = new Polynomial**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
- (*basis)[s] = new vector<const Polynomial*>(nFunction);
+ basis[s] = new Polynomial*[nFunction];
// Vertex Based //
for(unsigned int s = 0; s < nRefSpace; s++){
- (*(*basis)[s])[0] = new Polynomial(lagrange[0]);
- (*(*basis)[s])[1] = new Polynomial(lagrange[1]);
- (*(*basis)[s])[2] = new Polynomial(lagrange[2]);
+ basis[s][0] = new Polynomial(lagrange[0]);
+ basis[s][1] = new Polynomial(lagrange[1]);
+ basis[s][2] = new Polynomial(lagrange[2]);
}
// Edge Based //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex;
- for(unsigned int l = 1; l < order; l++){
- for(int e = 0; e < 3; e++){
- (*(*basis)[s])[i] =
+ for(int e = 0; e < 3; e++){
+ for(unsigned int l = 1; l < order; l++){
+ basis[s][i] =
new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
lagrange[(*(*edgeV[s])[e])[0]]
,
@@ -85,7 +85,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
- (*(*basis)[s])[i] =
+ basis[s][i] =
new Polynomial(intLegendre[l1].compose(lagrange[1] - lagrange[0],
lagrange[1] + lagrange[0])
*
@@ -108,10 +108,10 @@ TriNodeBasis::~TriNodeBasis(void){
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
- delete (*(*basis)[i])[j];
+ delete basis[i][j];
- delete (*basis)[i];
+ delete[] basis[i];
}
- delete basis;
+ delete[] basis;
}
diff --git a/FunctionSpace/TriNodeBasis.h b/FunctionSpace/TriNodeBasis.h
index 215731cb51..129ca4ce2c 100644
--- a/FunctionSpace/TriNodeBasis.h
+++ b/FunctionSpace/TriNodeBasis.h
@@ -1,7 +1,7 @@
#ifndef _TRINODEBASIS_H_
#define _TRINODEBASIS_H_
-#include "BasisScalar.h"
+#include "BasisHierarchicalScalar.h"
/**
@class TriNodeBasis
@@ -15,7 +15,7 @@
Basis for @em high @em order Polynomial%s generation.@n
*/
-class TriNodeBasis: public BasisScalar{
+class TriNodeBasis: public BasisHierarchicalScalar{
public:
//! @param order The order of the Basis
//!
--
GitLab