From 8d6b49c3985e8c851a09505aaae3fdca21d7061d Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 16 Oct 2012 15:50:10 +0000
Subject: [PATCH] Preliminary for 3D closure -- broken basis (except
TriNodeBasis)
---
FunctionSpace/BasisScalar.h | 37 ++++++++-
FunctionSpace/FunctionSpaceScalar.cpp | 35 +++++++--
FunctionSpace/TriNodeBasis.cpp | 108 +++++++++++++++-----------
3 files changed, 125 insertions(+), 55 deletions(-)
diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h
index 2b53803771..16b95fdd14 100644
--- a/FunctionSpace/BasisScalar.h
+++ b/FunctionSpace/BasisScalar.h
@@ -20,6 +20,11 @@
class BasisScalar: public Basis{
protected:
+ std::vector <const Polynomial*>* node;
+ std::vector<std::vector<const Polynomial*>*>* edge;
+ std::vector<std::vector<const Polynomial*>*>* face;
+ std::vector <const Polynomial*>* cell;
+
std::vector<const Polynomial*>* basis;
std::vector<const Polynomial*>* revBasis;
@@ -28,15 +33,16 @@ class BasisScalar: public Basis{
//!
virtual ~BasisScalar(void);
- //! @return Returns the set of @em Polynomial%s
- //! defining this (scalar) Basis
- const std::vector<const Polynomial*>& getFunctions(void) const;
-
//! @param closure A natural number
//! @return Returns the set of @em Polynomial%s
//! defining this (scalar) Basis, for the given closure
const std::vector<const Polynomial*>& getFunctions(unsigned int closure) const;
+ const std::vector <const Polynomial*>& getNodeFunctions(void) const;
+ const std::vector<std::vector<const Polynomial*>*>& getEdgeFunctions(void) const;
+ const std::vector<std::vector<const Polynomial*>*>& getFaceFunctions(void) const;
+ const std::vector <const Polynomial*>& getCellFunctions(void) const;
+
virtual std::string toString(void) const;
protected:
@@ -61,5 +67,28 @@ getFunctions(unsigned int closure) const{
return *revBasis;
}
+inline
+const std::vector<const Polynomial*>&
+BasisScalar::getNodeFunctions(void) const{
+ return *node;
+}
+
+inline
+const std::vector<std::vector<const Polynomial*>*>&
+BasisScalar::getEdgeFunctions(void) const{
+ return *edge;
+}
+
+inline
+const std::vector<std::vector<const Polynomial*>*>&
+BasisScalar::getFaceFunctions(void) const{
+ return *face;
+}
+
+inline
+const std::vector<const Polynomial*>&
+BasisScalar::getCellFunctions(void) const{
+ return *cell;
+}
#endif
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index a664302211..e5b7887f7c 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -17,21 +17,42 @@ getLocalFunctions(const MElement& element) const{
throw Exception("Element not found for closure");
vector<bool>* closure = it->second;
- const unsigned int size = closure->size();
// Get Basis //
- const vector<const Polynomial*>& basis = getBasis(element).getFunctions(0);
- const vector<const Polynomial*>& revBasis = getBasis(element).getFunctions(1);
+ const vector <const Polynomial*>& node = getBasis(element).getNodeFunctions();
+ const vector<vector<const Polynomial*>*>& edge = getBasis(element).getEdgeFunctions();
+ const vector <const Polynomial*>& cell = getBasis(element).getCellFunctions();
// Get Functions //
- vector<const Polynomial*> fun(size);
+ unsigned int i = 0;
+ const unsigned int nNode = node.size();
+ const unsigned int nEdge = edge[0]->size();
+ const unsigned int nCell = cell.size();
+
+ vector<const Polynomial*> fun(getBasis(element).getSize());
- for(unsigned int i = 0; i < size; i++)
+ // Vertex Based //
+ for(unsigned int j = 0; j < nNode; j++){
+ fun[i] = node[j];
+ i++;
+ }
+
+ // Edge Based //
+ for(unsigned int j = 0; j < nEdge; j++){
if((*closure)[i])
- fun[i] = basis[i];
+ fun[i] = (*edge[0])[j];
else
- fun[i] = revBasis[i];
+ fun[i] = (*edge[1])[j];
+
+ i++;
+ }
+
+ // Vertex Based //
+ for(unsigned int j = 0; j < nCell; j++){
+ fun[i] = cell[j];
+ i++;
+ }
// Return
return fun;
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 2b0ef40f39..13260e0f7e 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -1,6 +1,8 @@
#include "TriNodeBasis.h"
#include "Legendre.h"
+using namespace std;
+
TriNodeBasis::TriNodeBasis(const int order){
// Set Basis Type //
this->order = order;
@@ -20,8 +22,7 @@ TriNodeBasis::TriNodeBasis(const int order){
Polynomial* intLegendre = new Polynomial[order];
Polynomial lagrangeSum[3];
- Polynomial lagrangeSub[3];
- Polynomial rLagrangeSub[3];
+ Polynomial lagrangeSub[2][3];
// Classical and Intrated-Scaled Legendre Polynomial //
const int orderMinus = order - 1;
@@ -30,83 +31,102 @@ TriNodeBasis::TriNodeBasis(const int order){
Legendre::intScaled(intLegendre, order);
- // Basis (& revert) //
- basis = new std::vector<const Polynomial*>(size);
- revBasis = new std::vector<const Polynomial*>(size);
+ // Basis //
+ node = new vector<const Polynomial*>(nVertex);
+ edge = new vector<vector<const Polynomial*>*>(2);
+ face = new vector<vector<const Polynomial*>*>(0);
+ cell = new vector<const Polynomial*>(nCell);
+
+ (*edge)[0] = new vector<const Polynomial*>(nEdge);
+ (*edge)[1] = new vector<const Polynomial*>(nEdge);
+
- // Vertex Based (Lagrange) //
- (*basis)[0] =
+ // Vertex Based (Lagrange) //
+ (*node)[0] =
new Polynomial(Polynomial(1, 0, 0, 0) -
Polynomial(1, 1, 0, 0) -
Polynomial(1, 0, 1, 0));
- (*basis)[1] =
+ (*node)[1] =
new Polynomial(Polynomial(1, 1, 0, 0));
- (*basis)[2] =
+ (*node)[2] =
new Polynomial(Polynomial(1, 0, 1, 0));
+
- // Vertex Based (revert) //
- for(int i = 0; i < 3; i++)
- (*revBasis)[i] = (*basis)[i];
-
- // Lagrange Sum (& revert) //
+ // Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = *(*basis)[i] + *(*basis)[j];
+ lagrangeSum[i] = *(*node)[i] + *(*node)[j];
- // Lagrange Sub (& revert) //
+ // Lagrange Sub //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
- lagrangeSub[i] = *(*basis)[j] - *(*basis)[i];
- rLagrangeSub[i] = *(*basis)[i] - *(*basis)[j];
+ lagrangeSub[0][i] = *(*node)[j] - *(*node)[i];
+ lagrangeSub[1][i] = *(*node)[i] - *(*node)[j];
}
-
+
// Edge Based //
- int i = 3;
-
- for(int l = 1; l < order; l++){
- for(int e = 0; e < 3; e++){
- (*basis)[i] = new Polynomial(
- intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]));
-
- (*revBasis)[i] = new Polynomial(
- intLegendre[l].compose(rLagrangeSub[e], lagrangeSum[e]));
-
- i++;
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
+
+ for(int l = 1; l < order; l++){
+ for(int e = 0; e < 3; e++){
+ (*(*edge)[c])[i] =
+ new Polynomial(intLegendre[l].compose(lagrangeSub[c][e], lagrangeSum[e]));
+
+ i++;
+ }
}
}
// Cell Based //
- Polynomial p = *(*basis)[2] * 2 - Polynomial(1, 0, 0, 0);
+ Polynomial p = *(*node)[2] * 2 - Polynomial(1, 0, 0, 0);
const int orderMinusTwo = order - 2;
+
+ unsigned int i = 0;
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
- (*basis)[i] = new Polynomial(
- intLegendre[l1].compose(lagrangeSub[0], lagrangeSum[0]) *
- legendre[l2].compose(p) * *(*basis)[2]);
-
- (*revBasis)[i] = (*basis)[i];
+ (*cell)[i] =
+ new Polynomial(intLegendre[l1].compose(lagrangeSub[0][0], lagrangeSum[0]) *
+ legendre[l2].compose(p) * *(*node)[2]);
i++;
}
}
-
-
+
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
}
TriNodeBasis::~TriNodeBasis(void){
- for(int i = 0; i < size; i++){
- delete (*basis)[i];
+ // 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];
- if(i >= nVertex && i < nVertex + nEdge)
- delete (*revBasis)[i];
+ delete (*edge)[c];
}
+
+ delete edge;
+
+
+ // Face Based //
+ delete face;
+
+
+ // Cell Based //
+ for(int i = 0; i < nCell; i++)
+ delete (*cell)[i];
- delete basis;
- delete revBasis;
+ delete cell;
}
--
GitLab