From e66e9e63f95b0a64c4d9a2dc3aa683edea4bca18 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 13 Nov 2012 14:21:44 +0000
Subject: [PATCH] LineEdgeBasis: Convergence Test -- Seems OK
---
FunctionSpace/BasisGenerator.cpp | 7 +-
FunctionSpace/CMakeLists.txt | 2 +
FunctionSpace/LineEdgeBasis.cpp | 106 +++++++++++++++++++++++++++++
FunctionSpace/LineEdgeBasis.h | 34 +++++++++
FunctionSpace/LineNedelecBasis.cpp | 75 ++++++++++++++++++++
FunctionSpace/LineNedelecBasis.h | 25 +++++++
FunctionSpace/LineNodeBasis.cpp | 18 +++--
7 files changed, 259 insertions(+), 8 deletions(-)
create mode 100644 FunctionSpace/LineEdgeBasis.cpp
create mode 100644 FunctionSpace/LineEdgeBasis.h
create mode 100644 FunctionSpace/LineNedelecBasis.cpp
create mode 100644 FunctionSpace/LineNedelecBasis.h
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 58376a0404..299201fe79 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -3,6 +3,8 @@
#include "Exception.h"
#include "LineNodeBasis.h"
+#include "LineEdgeBasis.h"
+#include "LineNedelecBasis.h"
#include "QuadNodeBasis.h"
#include "QuadEdgeBasis.h"
@@ -43,7 +45,10 @@ Basis* BasisGenerator::linGen(int basisType,
int order){
switch(basisType){
case 0: return new LineNodeBasis(order);
- case 1: throw Exception("1-form not implemented on Lines");
+ case 1:
+ if (order == 0) return new LineNedelecBasis();
+ else return new LineEdgeBasis(order);
+
case 2: throw Exception("2-form not implemented on Lines");
case 3: throw Exception("3-form not implemented on Lines");
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index 04d39e2f1e..d02168cc8d 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -21,6 +21,8 @@ set(SRC
TriLagrangeBasis.cpp
LineNodeBasis.cpp
+ LineEdgeBasis.cpp
+ LineNedelecBasis.cpp
QuadNodeBasis.cpp
QuadEdgeBasis.cpp
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
new file mode 100644
index 0000000000..e881054042
--- /dev/null
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -0,0 +1,106 @@
+#include "LineEdgeBasis.h"
+#include "Legendre.h"
+
+using namespace std;
+
+LineEdgeBasis::LineEdgeBasis(int order){
+ // Set Basis Type //
+ this->order = order;
+
+ type = 1;
+ dim = 1;
+
+ nVertex = 0;
+ nEdge = (order + 1);
+ nFace = 0;
+ nCell = 0;
+
+ nEdgeClosure = 2;
+ nFaceClosure = 0;
+
+ size = nVertex + nEdge + nFace + nCell;
+
+ // Alloc Temporary Space //
+ const int orderPlus = order + 1;
+ Polynomial* intLegendre = new Polynomial[orderPlus];
+
+ const Polynomial plusOneHalf(+0.5, 0, 0, 0);
+ const Polynomial minusOneHalf(-0.5, 0, 0, 0);
+ const Polynomial zero(0, 0, 0, 0);
+
+ const Polynomial x[2] = {
+ Polynomial(+1, 1, 0, 0),
+ Polynomial(-1, 1, 0, 0)
+ };
+
+ // Legendre Polynomial //
+ Legendre::integrated(intLegendre, orderPlus);
+
+ // Permutations //
+ const int permutation[2] = {0, 1};
+
+ // Basis //
+ node = new vector<vector<Polynomial>*>(nVertex);
+ edge = new vector<vector<vector<Polynomial>*>*>(nEdgeClosure);
+ face = new vector<vector<vector<Polynomial>*>*>(nFaceClosure);
+ cell = new vector<vector<Polynomial>*>(nCell);
+
+ (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
+ (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+
+
+ // Edge Based (Nedelec) //
+ (*(*edge)[0])[0] = new vector<Polynomial>(3);
+ (*(*edge)[0])[0]->at(0) = plusOneHalf;
+ (*(*edge)[0])[0]->at(1) = zero;
+ (*(*edge)[0])[0]->at(2) = zero;
+
+ (*(*edge)[1])[0] = new vector<Polynomial>(3);
+ (*(*edge)[1])[0]->at(0) = minusOneHalf;
+ (*(*edge)[1])[0]->at(1) = zero;
+ (*(*edge)[1])[0]->at(2) = zero;
+
+
+ // Edge Based (High Order) //
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
+
+ for(int l = 1; l < orderPlus; l++){
+ (*(*edge)[c])[i + 1] =
+ new vector<Polynomial>((intLegendre[l].compose(x[permutation[c]])).gradient());
+
+ i++;
+ }
+ }
+
+
+ // Free Temporary Space //
+ delete[] intLegendre;
+}
+
+LineEdgeBasis::~LineEdgeBasis(void){
+ // Vertex Based //
+ for(int i = 0; i < nVertex; i++)
+ delete (*node)[i];
+
+ delete node;
+
+
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[i];
+
+ delete (*edge)[c];
+ }
+
+ delete edge;
+
+
+ // Face Based //
+ delete face;
+
+
+ // Cell Based //
+ delete cell;
+}
diff --git a/FunctionSpace/LineEdgeBasis.h b/FunctionSpace/LineEdgeBasis.h
new file mode 100644
index 0000000000..c5e7e794f4
--- /dev/null
+++ b/FunctionSpace/LineEdgeBasis.h
@@ -0,0 +1,34 @@
+#ifndef _LINEEDGEBASIS_H_
+#define _LINEEDGEBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+ @class LineEdgeBasis
+ @brief An Edge Basis for Lines
+
+ This class can instantiate an Edge-Based Basis
+ (high or low order) for Lines.@n
+
+ It uses an @em adaptation of
+ <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
+ Basis for @em high @em order Polynomial%s generation.@n
+
+ This Basis is a restriction of a Quad Basis to @f$y = 0@f$.@n
+
+ It also uses the following mapping: @f$x = \frac{u + 1}{2}@f$.
+*/
+
+class LineEdgeBasis: public BasisVector{
+ public:
+ //! @param order The order of the Basis
+ //!
+ //! Returns a new Edge-Basis for Lines of the given order
+ LineEdgeBasis(int order);
+
+ //! Deletes this Basis
+ //!
+ virtual ~LineEdgeBasis(void);
+};
+
+#endif
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
new file mode 100644
index 0000000000..04b9a79533
--- /dev/null
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -0,0 +1,75 @@
+#include "LineNedelecBasis.h"
+#include "Legendre.h"
+
+using namespace std;
+
+LineNedelecBasis::LineNedelecBasis(void){
+ // Set Basis Type //
+ order = 1;
+
+ type = 1;
+ dim = 1;
+
+ nVertex = 0;
+ nEdge = 1;
+ nFace = 0;
+ nCell = 0;
+
+ nEdgeClosure = 2;
+ nFaceClosure = 0;
+
+ size = 1;
+
+ // Alloc Temporary Space //
+ const Polynomial plusOneHalf(+0.5, 0, 0, 0);
+ const Polynomial minusOneHalf(-0.5, 0, 0, 0);
+ const Polynomial zero(0, 0, 0, 0);
+
+ // Basis //
+ node = new vector<vector<Polynomial>*>(nVertex);
+ edge = new vector<vector<vector<Polynomial>*>*>(nEdgeClosure);
+ face = new vector<vector<vector<Polynomial>*>*>(nFaceClosure);
+ cell = new vector<vector<Polynomial>*>(nCell);
+
+ (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
+ (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+
+
+ // Nedelec //
+ (*(*edge)[0])[0] = new vector<Polynomial>(3);
+ (*(*edge)[0])[0]->at(0) = plusOneHalf;
+ (*(*edge)[0])[0]->at(1) = zero;
+ (*(*edge)[0])[0]->at(2) = zero;
+
+ (*(*edge)[1])[0] = new vector<Polynomial>(3);
+ (*(*edge)[1])[0]->at(0) = minusOneHalf;
+ (*(*edge)[1])[0]->at(1) = zero;
+ (*(*edge)[1])[0]->at(2) = zero;
+}
+
+LineNedelecBasis::~LineNedelecBasis(void){
+ // Vertex Based //
+ for(int i = 0; i < nVertex; i++)
+ delete (*node)[i];
+
+ delete node;
+
+
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[i];
+
+ delete (*edge)[c];
+ }
+
+ delete edge;
+
+
+ // Face Based //
+ delete face;
+
+
+ // Cell Based //
+ delete cell;
+}
diff --git a/FunctionSpace/LineNedelecBasis.h b/FunctionSpace/LineNedelecBasis.h
new file mode 100644
index 0000000000..8579b1a470
--- /dev/null
+++ b/FunctionSpace/LineNedelecBasis.h
@@ -0,0 +1,25 @@
+#ifndef _LINENEDELECBASIS_H_
+#define _LINENEDELECBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+ @class LineNedelecBasis
+ @brief Nedelec Basis for Lines
+
+ This class can instantiate a Nedelec Basis
+ for Lines.@n
+*/
+
+class LineNedelecBasis: public BasisVector{
+ public:
+ //! Returns a new Nedelec Basis for Lines
+ //!
+ LineNedelecBasis(void);
+
+ //! Deletes this Basis
+ //!
+ virtual ~LineNedelecBasis(void);
+};
+
+#endif
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index 762d68b31f..4b94c8c5c2 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -20,10 +20,20 @@ LineNodeBasis::LineNodeBasis(int order){
size = nVertex + nEdge + nFace + nCell;
- // Legendre Polynomial //
+ // Alloc Temporary Space //
Polynomial* intLegendre = new Polynomial[order];
+
+ const Polynomial x[2] = {
+ Polynomial(+1, 1, 0, 0),
+ Polynomial(-1, 1, 0, 0)
+ };
+
+ // Legendre Polynomial //
Legendre::integrated(intLegendre, order);
+ // Permutations //
+ const int permutation[2] = {0, 1};
+
// Basis //
node = new vector<Polynomial*>(nVertex);
edge = new vector<vector<Polynomial*>*>(nEdgeClosure);
@@ -45,12 +55,6 @@ LineNodeBasis::LineNodeBasis(int order){
// Edge Based //
- const int permutation[2] = {0, 1};
- const Polynomial x[2] = {
- Polynomial(+1, 1, 0, 0),
- Polynomial(-1, 1, 0, 0)
- };
-
for(int c = 0; c < 2; c++){
unsigned int i = 0;
--
GitLab