From a88a2028824810a3ee478d3860e77feda03a1c5c Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Thu, 18 Jul 2013 07:50:22 +0000
Subject: [PATCH] Adding Quad and Tet EdgeBasis with nex ReferenceSpace --
Seems OK
---
FunctionSpace/QuadEdgeBasis.cpp | 100 ++++++++++++++---------------
FunctionSpace/QuadEdgeBasis.h | 2 +-
FunctionSpace/QuadNedelecBasis.cpp | 30 ++++-----
FunctionSpace/TetEdgeBasis.cpp | 92 +++++++++++++-------------
FunctionSpace/TetEdgeBasis.h | 2 +-
FunctionSpace/TetNedelecBasis.cpp | 34 +++++-----
6 files changed, 122 insertions(+), 138 deletions(-)
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 5fb906d080..c4dc7d4929 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -4,17 +4,16 @@
using namespace std;
-QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
- /*
+QuadEdgeBasis::QuadEdgeBasis(size_t order){
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
- const vector<const vector<const vector<unsigned int>*>*>&
- faceV = refSpace->getAllFace();
+ const vector<vector<vector<size_t> > >&
+ faceIdx = refSpace->getFaceNodeIndex();
// Set Basis Type //
this->order = order;
@@ -29,7 +28,7 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
nFunction = nVertex + nEdge + nFace + nCell;
// Legendre Polynomial //
- const unsigned int orderPlus = order + 1;
+ const size_t orderPlus = order + 1;
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
@@ -71,23 +70,23 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = 0;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = 0;
- for(unsigned int e = 0; e < 4; e++){
- for(unsigned int l = 0; l < orderPlus; l++){
+ for(size_t e = 0; e < 4; e++){
+ for(size_t l = 0; l < orderPlus; l++){
// Nedelec
if(l == 0){
- Polynomial lambda = (lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]]) * 0.5;
+ Polynomial lambda = (lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]]) * 0.5;
basis[s][i] =
- new vector<Polynomial>((lifting[(*(*edgeV[s])[e])[1]] -
- lifting[(*(*edgeV[s])[e])[0]]).gradient());
+ new vector<Polynomial>((lifting[edgeIdx[s][e][1]] -
+ lifting[edgeIdx[s][e][0]]).gradient());
basis[s][i]->at(0).mul(lambda);
basis[s][i]->at(1).mul(lambda);
@@ -98,11 +97,11 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
else{
basis[s][i] =
new vector<Polynomial>
- ((intLegendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
- lifting[(*(*edgeV[s])[e])[0]])
+ ((intLegendre[l].compose(lifting[edgeIdx[s][e][1]] -
+ lifting[edgeIdx[s][e][0]])
*
- (lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]])).gradient());
+ (lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]])).gradient());
}
i++;
@@ -112,37 +111,37 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// Face Based //
- // NB: We use (*(*faceV[s])[f])[]
+ // NB: We use (*(*faceIdx[s])[f])[]
// where f = 0, because triangles
// have only ONE face: the face '0'
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nEdge;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nEdge;
// Type 1
- for(unsigned int l1 = 1; l1 < orderPlus; l1++){
- for(unsigned int l2 = 1; l2 < orderPlus; l2++){
+ for(size_t l1 = 1; l1 < orderPlus; l1++){
+ for(size_t l2 = 1; l2 < orderPlus; l2++){
basis[s][i] =
new vector<Polynomial>
- ((intLegendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[1]])
+ ((intLegendre[l1].compose(lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][1]])
*
- intLegendre[l2].compose(lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[3]])).gradient());
+ intLegendre[l2].compose(lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][3]])).gradient());
i++;
}
}
// Type 2
- for(unsigned int l1 = 1; l1 < orderPlus; l1++){
- for(unsigned int l2 = 1; l2 < orderPlus; l2++){
- Polynomial pOne = lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[1]];
+ for(size_t l1 = 1; l1 < orderPlus; l1++){
+ for(size_t l2 = 1; l2 < orderPlus; l2++){
+ Polynomial pOne = lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][1]];
- Polynomial pTwo = lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[3]];
+ Polynomial pTwo = lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][3]];
Polynomial lOne = legendre[l1].compose(pOne) *
intLegendre[l2].compose(pTwo);
@@ -172,15 +171,15 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
}
// Type 3.1
- for(unsigned int l1 = 1; l1 < orderPlus; l1++){
+ for(size_t l1 = 1; l1 < orderPlus; l1++){
Polynomial tmp =
- intLegendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[3]]);
+ intLegendre[l1].compose(lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][3]]);
basis[s][i] =
- new vector<Polynomial>((lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[1]]).gradient());
+ new vector<Polynomial>((lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][1]]).gradient());
basis[s][i]->at(0).mul(tmp);
basis[s][i]->at(1).mul(tmp);
@@ -190,15 +189,15 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
}
// Type 3.2
- for(unsigned int l1 = 1; l1 < orderPlus; l1++){
+ for(size_t l1 = 1; l1 < orderPlus; l1++){
Polynomial tmp =
- intLegendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[1]]);
+ intLegendre[l1].compose(lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][1]]);
basis[s][i] =
- new vector<Polynomial>((lifting[(*(*faceV[s])[0])[0]] -
- lifting[(*(*faceV[s])[0])[3]]).gradient());
+ new vector<Polynomial>((lifting[faceIdx[s][0][0]] -
+ lifting[faceIdx[s][0][3]]).gradient());
basis[s][i]->at(0).mul(tmp);
basis[s][i]->at(1).mul(tmp);
@@ -221,8 +220,8 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int i = 0; i < nFunction; i++){
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t i = 0; i < nFunction; i++){
vector<Polynomial>* old;
vector<Polynomial> nxt(3);
@@ -239,22 +238,19 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
- */
}
QuadEdgeBasis::~QuadEdgeBasis(void){
- /*
// ReferenceSpace //
delete refSpace;
// Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
+ for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
- */
}
diff --git a/FunctionSpace/QuadEdgeBasis.h b/FunctionSpace/QuadEdgeBasis.h
index 45b5060176..95d6edb66d 100644
--- a/FunctionSpace/QuadEdgeBasis.h
+++ b/FunctionSpace/QuadEdgeBasis.h
@@ -24,7 +24,7 @@ class QuadEdgeBasis: public BasisHierarchical1From{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Quads of the given order
- QuadEdgeBasis(unsigned int order);
+ QuadEdgeBasis(size_t order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/QuadNedelecBasis.cpp b/FunctionSpace/QuadNedelecBasis.cpp
index 42eda0550c..a91c2e4de1 100644
--- a/FunctionSpace/QuadNedelecBasis.cpp
+++ b/FunctionSpace/QuadNedelecBasis.cpp
@@ -5,13 +5,12 @@
using namespace std;
QuadNedelecBasis::QuadNedelecBasis(void){
- /*
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
// Set Basis Type //
order = 0;
@@ -59,19 +58,19 @@ QuadNedelecBasis::QuadNedelecBasis(void){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int e = 0; e < 4; e++){
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t e = 0; e < 4; e++){
- Polynomial lambda = (lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]]) * 0.5;
+ Polynomial lambda = (lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]]) * 0.5;
basis[s][e] =
- new vector<Polynomial>((lifting[(*(*edgeV[s])[e])[1]] -
- lifting[(*(*edgeV[s])[e])[0]]).gradient());
+ new vector<Polynomial>((lifting[edgeIdx[s][e][1]] -
+ lifting[edgeIdx[s][e][0]]).gradient());
basis[s][e]->at(0).mul(lambda);
basis[s][e]->at(1).mul(lambda);
@@ -92,8 +91,8 @@ QuadNedelecBasis::QuadNedelecBasis(void){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0));
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int i = 0; i < nFunction; i++){
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t i = 0; i < nFunction; i++){
vector<Polynomial>* old;
vector<Polynomial> nxt(3);
@@ -106,22 +105,19 @@ QuadNedelecBasis::QuadNedelecBasis(void){
delete old;
}
}
- */
}
QuadNedelecBasis::~QuadNedelecBasis(void){
- /*
// ReferenceSpace //
delete refSpace;
// Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
+ for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
- */
}
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 92bb3e6649..5d3fd21a79 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -4,17 +4,16 @@
using namespace std;
-TetEdgeBasis::TetEdgeBasis(unsigned int order){
- /*
+TetEdgeBasis::TetEdgeBasis(size_t order){
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
- const vector<const vector<const vector<unsigned int>*>*>&
- faceV = refSpace->getAllFace();
+ const vector<vector<vector<size_t> > >&
+ faceIdx = refSpace->getFaceNodeIndex();
// Set Basis Type //
this->order = order;
@@ -61,27 +60,27 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = 0;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = 0;
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();
+ vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[edgeIdx[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]]);
+ tmp1[0].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[1].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[2].mul(lagrange[edgeIdx[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].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[1].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[2].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[0].sub(tmp1[0]);
tmp2[1].sub(tmp1[1]);
@@ -94,11 +93,11 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
else{
basis[s][i] =
new vector<Polynomial>
- ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
- lagrange[(*(*edgeV[s])[e])[1]]
+ ((intLegendre[l].compose(lagrange[edgeIdx[s][e][0]] -
+ lagrange[edgeIdx[s][e][1]]
,
- lagrange[(*(*edgeV[s])[e])[0]] +
- lagrange[(*(*edgeV[s])[e])[1]])).gradient());
+ lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]])).gradient());
}
i++;
}
@@ -108,27 +107,27 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Face Based //
// TO CHECK: Are Triangles face matching tets ?
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nEdge;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nEdge;
for(int f = 0; f < 4; f++){
- for(unsigned int l1 = 1; l1 < order; l1++){
+ for(size_t 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]];
+ lagrange[faceIdx[s][f][0]] +
+ lagrange[faceIdx[s][f][1]] +
+ lagrange[faceIdx[s][f][2]];
Polynomial u =
- intLegendre[l1].compose(lagrange[(*(*faceV[s])[f])[0]] -
- lagrange[(*(*faceV[s])[f])[1]]
+ intLegendre[l1].compose(lagrange[faceIdx[s][f][0]] -
+ lagrange[faceIdx[s][f][1]]
,
- lagrange[(*(*faceV[s])[f])[0]] +
- lagrange[(*(*faceV[s])[f])[1]]);
+ lagrange[faceIdx[s][f][0]] +
+ lagrange[faceIdx[s][f][1]]);
Polynomial v =
- lagrange[(*(*faceV[s])[f])[2]] *
- sclLegendre[l2].compose(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum);
+ lagrange[faceIdx[s][f][2]] *
+ sclLegendre[l2].compose(lagrange[faceIdx[s][f][2]] * 2 - sum, sum);
// Preliminary Type 2
vector<Polynomial> gradU = u.gradient();
@@ -150,18 +149,18 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
subGradUV[2].sub(uGradV[2]);
// Preliminary Type 3
- vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[f])[0]].gradient();
- vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[f])[1]].gradient();
+ vector<Polynomial> gradL1 = lagrange[faceIdx[s][f][0]].gradient();
+ vector<Polynomial> gradL2 = lagrange[faceIdx[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]]);
+ l2GradL1[0].mul(lagrange[faceIdx[s][f][1]]);
+ l2GradL1[1].mul(lagrange[faceIdx[s][f][1]]);
+ l2GradL1[2].mul(lagrange[faceIdx[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]]);
+ l1GradL2[0].mul(lagrange[faceIdx[s][f][0]]);
+ l1GradL2[1].mul(lagrange[faceIdx[s][f][0]]);
+ l1GradL2[2].mul(lagrange[faceIdx[s][f][0]]);
vector<Polynomial> subGradL1L2V(l2GradL1);
subGradL1L2V[0].sub(l1GradL2[0]);
@@ -200,8 +199,8 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
// Cell Based //
const Polynomial one(1, 0, 0, 0);
- for(unsigned int s = 0; s < nRefSpace; s++){
- unsigned int i = nEdge + nFace;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nEdge + nFace;
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
@@ -322,22 +321,19 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
delete[] legendre;
delete[] sclLegendre;
delete[] intLegendre;
- */
}
TetEdgeBasis::~TetEdgeBasis(void){
- /*
// ReferenceSpace //
delete refSpace;
// Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
+ for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
- */
}
diff --git a/FunctionSpace/TetEdgeBasis.h b/FunctionSpace/TetEdgeBasis.h
index 61e4f982d2..d8e7856fd4 100644
--- a/FunctionSpace/TetEdgeBasis.h
+++ b/FunctionSpace/TetEdgeBasis.h
@@ -20,7 +20,7 @@ class TetEdgeBasis: public BasisHierarchical1From{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Tetrahedra of the given order
- TetEdgeBasis(unsigned int order);
+ TetEdgeBasis(size_t order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TetNedelecBasis.cpp b/FunctionSpace/TetNedelecBasis.cpp
index c803758152..2b06274cb9 100644
--- a/FunctionSpace/TetNedelecBasis.cpp
+++ b/FunctionSpace/TetNedelecBasis.cpp
@@ -4,13 +4,12 @@
using namespace std;
TetNedelecBasis::TetNedelecBasis(void){
- /*
// Reference Space //
refSpace = new TetReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
- const vector<const vector<const vector<unsigned int>*>*>&
- edgeV = refSpace->getAllEdge();
+ const vector<vector<vector<size_t> > >&
+ edgeIdx = refSpace->getEdgeNodeIndex();
// Set Basis Type //
this->order = 0;
@@ -43,22 +42,22 @@ TetNedelecBasis::TetNedelecBasis(void){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
- for(unsigned int s = 0; s < nRefSpace; s++)
+ for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
- for(unsigned int s = 0; s < nRefSpace; s++){
- for(unsigned int e = 0; e < 6; e++){
- vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
- vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t e = 0; e < 6; e++){
+ vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[edgeIdx[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]]);
+ tmp1[0].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[1].mul(lagrange[edgeIdx[s][e][0]]);
+ tmp1[2].mul(lagrange[edgeIdx[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].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[1].mul(lagrange[edgeIdx[s][e][1]]);
+ tmp2[2].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[0].sub(tmp1[0]);
tmp2[1].sub(tmp1[1]);
@@ -67,22 +66,19 @@ TetNedelecBasis::TetNedelecBasis(void){
basis[s][e] = new vector<Polynomial>(tmp2);
}
}
- */
}
TetNedelecBasis::~TetNedelecBasis(void){
- /*
// ReferenceSpace //
delete refSpace;
// Basis //
- for(unsigned int i = 0; i < nRefSpace; i++){
- for(unsigned int j = 0; j < nFunction; j++)
+ for(size_t i = 0; i < nRefSpace; i++){
+ for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
- */
}
--
GitLab