From 4caf2d4db43f937c52a6280d3c34f28bedfb3192 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Mon, 10 Sep 2012 15:53:22 +0000
Subject: [PATCH] Reimplementing Basis -- Needed for closure
---
FunctionSpace/BasisScalar.h | 6 +-
FunctionSpace/BasisTest.cpp | 4 +-
FunctionSpace/BasisVector.cpp | 1 -
FunctionSpace/BasisVector.h | 6 +-
FunctionSpace/FunctionSpaceEdge.cpp | 6 +-
FunctionSpace/FunctionSpaceNode.cpp | 6 +-
FunctionSpace/HexEdgeBasis.cpp | 193 +++++++++-------------------
FunctionSpace/HexNodeBasis.cpp | 190 ++++++---------------------
FunctionSpace/QuadEdgeBasis.cpp | 139 +++++---------------
FunctionSpace/QuadNodeBasis.cpp | 96 +++-----------
FunctionSpace/TriEdgeBasis.cpp | 145 ++++++---------------
FunctionSpace/TriNedelecBasis.cpp | 104 +++------------
FunctionSpace/TriNodeBasis.cpp | 96 +++-----------
13 files changed, 240 insertions(+), 752 deletions(-)
diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h
index a5bd26afdf..6f156f4dbe 100644
--- a/FunctionSpace/BasisScalar.h
+++ b/FunctionSpace/BasisScalar.h
@@ -20,7 +20,7 @@
class BasisScalar: public Basis{
protected:
- std::vector<Polynomial>* basis;
+ std::vector<const Polynomial*>* basis;
public:
//! Deletes this BasisScalar
@@ -29,7 +29,7 @@ class BasisScalar: public Basis{
//! @return Returns the set of @em Polynomial%s
//! defining this (scalar) Basis
- const std::vector<Polynomial>& getFunctions(void) const;
+ const std::vector<const Polynomial*>& getFunctions(void) const;
protected:
//! @internal
@@ -44,7 +44,7 @@ class BasisScalar: public Basis{
//////////////////////
inline
-const std::vector<Polynomial>& BasisScalar::
+const std::vector<const Polynomial*>& BasisScalar::
getFunctions(void) const{
return *basis;
diff --git a/FunctionSpace/BasisTest.cpp b/FunctionSpace/BasisTest.cpp
index 8f7de4f3a9..906af7f59d 100644
--- a/FunctionSpace/BasisTest.cpp
+++ b/FunctionSpace/BasisTest.cpp
@@ -32,8 +32,10 @@ int basisTest(int argc, char** argv){
writer.setDomain(goe.getAll());
// Plot Basis //
- TriNodeBasis b(3);
+ HexEdgeBasis b(3);
+ cout << "Size: " << b.getSize() << endl;
+
PlotBasis plot(b, goe, writer);
plot.plot("basis");
diff --git a/FunctionSpace/BasisVector.cpp b/FunctionSpace/BasisVector.cpp
index 97c61e7209..e1d389822e 100644
--- a/FunctionSpace/BasisVector.cpp
+++ b/FunctionSpace/BasisVector.cpp
@@ -6,4 +6,3 @@ BasisVector::BasisVector(void){
BasisVector::~BasisVector(void){
}
-
diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h
index 32ef5b9193..85e3491147 100644
--- a/FunctionSpace/BasisVector.h
+++ b/FunctionSpace/BasisVector.h
@@ -20,7 +20,7 @@
class BasisVector: public Basis{
protected:
- std::vector<std::vector<Polynomial> >* basis;
+ std::vector<const std::vector<Polynomial>*>* basis;
public:
//! Deletes this BasisVector
@@ -29,7 +29,7 @@ class BasisVector: public Basis{
//! @return Returns the set of @em Polynomial%s
//! defining this (vectorial) Basis
- const std::vector<std::vector<Polynomial> >& getFunctions(void) const;
+ const std::vector<const std::vector<Polynomial>*>& getFunctions(void) const;
protected:
//! @internal
@@ -44,7 +44,7 @@ class BasisVector: public Basis{
//////////////////////
inline
-const std::vector<std::vector<Polynomial> >& BasisVector::
+const std::vector<const std::vector<Polynomial>*>& BasisVector::
getFunctions(void) const{
return *basis;
diff --git a/FunctionSpace/FunctionSpaceEdge.cpp b/FunctionSpace/FunctionSpaceEdge.cpp
index 8386525cea..6673e9f467 100644
--- a/FunctionSpace/FunctionSpaceEdge.cpp
+++ b/FunctionSpace/FunctionSpaceEdge.cpp
@@ -25,8 +25,8 @@ interpolate(const MElement& element,
const_cast<MElement&>(element);
// Get Basis Functions //
- const vector<vector<Polynomial> >& fun = getBasis(element).getFunctions();
- const unsigned int nFun = fun.size();
+ const vector<const vector<Polynomial>*>& fun = getBasis(element).getFunctions();
+ const unsigned int nFun = fun.size();
// Get Reference coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
@@ -47,7 +47,7 @@ interpolate(const MElement& element,
for(unsigned int i = 0; i < nFun; i++){
fullVector<double> vi =
- Mapper::grad(Polynomial::at(fun[i], uvw[0], uvw[1], uvw[2]),
+ Mapper::grad(Polynomial::at(*fun[i], uvw[0], uvw[1], uvw[2]),
invJac);
vi.scale(coef[i]);
diff --git a/FunctionSpace/FunctionSpaceNode.cpp b/FunctionSpace/FunctionSpaceNode.cpp
index 7db4f3f71d..baaadc1332 100644
--- a/FunctionSpace/FunctionSpaceNode.cpp
+++ b/FunctionSpace/FunctionSpaceNode.cpp
@@ -25,8 +25,8 @@ interpolate(const MElement& element,
const_cast<MElement&>(element);
// Get Basis Functions //
- const vector<Polynomial>& fun = getBasis(element).getFunctions();
- const unsigned int nFun = fun.size();
+ const vector<const Polynomial*>& fun = getBasis(element).getFunctions();
+ const unsigned int nFun = fun.size();
// Get Reference coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
@@ -39,7 +39,7 @@ interpolate(const MElement& element,
for(unsigned int i = 0; i < nFun; i++)
val +=
- fun[i].at(uvw[0], uvw[1], uvw[2]) * coef[i];
+ fun[i]->at(uvw[0], uvw[1], uvw[2]) * coef[i];
// Return Interpolated Value //
return val;
diff --git a/FunctionSpace/HexEdgeBasis.cpp b/FunctionSpace/HexEdgeBasis.cpp
index c9a57d09c9..22ca9fddde 100644
--- a/FunctionSpace/HexEdgeBasis.cpp
+++ b/FunctionSpace/HexEdgeBasis.cpp
@@ -148,11 +148,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
liftingSub[i] = lifting[edge1[i]] - lifting[edge2[i]];
- // Basis //
- basis = new std::vector<std::vector<Polynomial> >(size);
-
- for(int i = 0; i < size; i++)
- (*basis)[i].resize(3);
+ // Basis (temporary --- *no* const) //
+ std::vector<std::vector<Polynomial>*> basis(size);
// Edge Based (Nedelec) //
@@ -160,16 +157,16 @@ HexEdgeBasis::HexEdgeBasis(const int order){
Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 12; e++){
- (*basis)[i] =
- (liftingSub[e]).gradient();
+ basis[i] =
+ new std::vector<Polynomial>((liftingSub[e]).gradient());
- (*basis)[i][0].mul(lagrangeSum[e]);
- (*basis)[i][1].mul(lagrangeSum[e]);
- (*basis)[i][2].mul(lagrangeSum[e]);
+ basis[i]->at(0).mul(lagrangeSum[e]);
+ basis[i]->at(1).mul(lagrangeSum[e]);
+ basis[i]->at(2).mul(lagrangeSum[e]);
- (*basis)[i][0].mul(oneHalf);
- (*basis)[i][1].mul(oneHalf);
- (*basis)[i][2].mul(oneHalf);
+ basis[i]->at(0).mul(oneHalf);
+ basis[i]->at(1).mul(oneHalf);
+ basis[i]->at(2).mul(oneHalf);
i++;
}
@@ -178,8 +175,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 12; e++){
- (*basis)[i] =
- (intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient();
+ basis[i] =
+ new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
i++;
}
@@ -237,9 +234,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
for(int f = 0; f < 6; f++){
- (*basis)[i] = (iLegendreXi[l1][f] *
- iLegendreEta[l2][f] *
- lambda[f]).gradient();
+ basis[i] = new std::vector<Polynomial>((iLegendreXi[l1][f] *
+ iLegendreEta[l2][f] *
+ lambda[f]).gradient());
i++;
}
@@ -251,6 +248,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
for(int f = 0; f < 6; f++){
+ basis[i] = new std::vector<Polynomial>(3);
+
Polynomial tmp1 =
legendreXi[l1][f] *
iLegendreEta[l2][f];
@@ -269,9 +268,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){
gr2[1].mul(tmp2);
gr2[2].mul(tmp2);
- (*basis)[i][0] = (gr1[0] - gr2[0]) * lambda[f];
- (*basis)[i][1] = (gr1[1] - gr2[1]) * lambda[f];
- (*basis)[i][2] = (gr1[2] - gr2[2]) * lambda[f];
+ basis[i]->at(0) = (gr1[0] - gr2[0]) * lambda[f];
+ basis[i]->at(1) = (gr1[1] - gr2[1]) * lambda[f];
+ basis[i]->at(2) = (gr1[2] - gr2[2]) * lambda[f];
i++;
}
@@ -284,11 +283,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int f = 0; f < 6; f++){
Polynomial tmp = iLegendreEta[l][f] * lambda[f];
- (*basis)[i] = grXi[f];
+ basis[i] = new std::vector<Polynomial>(grXi[f]);
- (*basis)[i][0].mul(tmp);
- (*basis)[i][1].mul(tmp);
- (*basis)[i][2].mul(tmp);
+ basis[i]->at(0).mul(tmp);
+ basis[i]->at(1).mul(tmp);
+ basis[i]->at(2).mul(tmp);
i++;
}
@@ -300,11 +299,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int f = 0; f < 6; f++){
Polynomial tmp = iLegendreXi[l][f] * lambda[f];
- (*basis)[i] = grEta[f];
+ basis[i] = new std::vector<Polynomial>(grEta[f]);
- (*basis)[i][0].mul(tmp);
- (*basis)[i][1].mul(tmp);
- (*basis)[i][2].mul(tmp);
+ basis[i]->at(0).mul(tmp);
+ basis[i]->at(1).mul(tmp);
+ basis[i]->at(2).mul(tmp);
i++;
}
@@ -335,9 +334,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
for(int l3 = 1; l3 < orderPlus; l3++){
- (*basis)[i] = (iLegendreX[l1] *
- iLegendreY[l2] *
- iLegendreZ[l3]).gradient();
+ basis[i] = new std::vector<Polynomial>((iLegendreX[l1] *
+ iLegendreY[l2] *
+ iLegendreZ[l3]).gradient());
i++;
cellNumber++;
@@ -348,11 +347,13 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 2 -- First Part) //
for(int j = 0; j < cellNumber; j++){
+ basis[i] = new std::vector<Polynomial>(3);
+
int off = j + cellStart;
- (*basis)[i][0] = (*basis)[off][0];
- (*basis)[i][1] = (*basis)[off][1] * Polynomial(-1, 0, 0, 0);
- (*basis)[i][2] = (*basis)[off][2];
+ basis[i]->at(0) = basis[off]->at(0);
+ basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
+ basis[i]->at(2) = basis[off]->at(2);
i++;
}
@@ -360,11 +361,13 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 2 -- Second Part) //
for(int j = 0; j < cellNumber; j++){
+ basis[i] = new std::vector<Polynomial>(3);
+
int off = j + cellStart;
- (*basis)[i][0] = (*basis)[off][0];
- (*basis)[i][1] = (*basis)[off][1] * Polynomial(-1, 0, 0, 0);
- (*basis)[i][2] = (*basis)[off][2] * Polynomial(-1, 0, 0, 0);
+ basis[i]->at(0) = basis[off]->at(0);
+ basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
+ basis[i]->at(2) = basis[off]->at(2) * Polynomial(-1, 0, 0, 0);
i++;
}
@@ -373,9 +376,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 3 -- First Part) //
for(int l2 = 1; l2 < orderPlus; l2++){
for(int l3 = 1; l3 < orderPlus; l3++){
- (*basis)[i][0] = iLegendreY[l2] * iLegendreZ[l3];
- (*basis)[i][1] = zero;
- (*basis)[i][2] = zero;
+ basis[i] = new std::vector<Polynomial>(3);
+
+ basis[i]->at(0) = iLegendreY[l2] * iLegendreZ[l3];
+ basis[i]->at(1) = zero;
+ basis[i]->at(2) = zero;
i++;
}
@@ -385,9 +390,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 3 -- Second Part) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l3 = 1; l3 < orderPlus; l3++){
- (*basis)[i][0] = zero;
- (*basis)[i][1] = iLegendreX[l1] * iLegendreZ[l3];
- (*basis)[i][2] = zero;
+ basis[i] = new std::vector<Polynomial>(3);
+
+ basis[i]->at(0) = zero;
+ basis[i]->at(1) = iLegendreX[l1] * iLegendreZ[l3];
+ basis[i]->at(2) = zero;
i++;
}
@@ -397,9 +404,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 3 -- Thrid Part) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
- (*basis)[i][0] = zero;
- (*basis)[i][1] = zero;
- (*basis)[i][2] = iLegendreX[l1] * iLegendreY[l2];
+ basis[i] = new std::vector<Polynomial>(3);
+
+ basis[i]->at(0) = zero;
+ basis[i]->at(1) = zero;
+ basis[i]->at(2) = iLegendreX[l1] * iLegendreY[l2];
i++;
}
@@ -438,90 +447,16 @@ HexEdgeBasis::HexEdgeBasis(const int order){
delete[] iLegendreX;
delete[] iLegendreY;
delete[] iLegendreZ;
-}
-HexEdgeBasis::~HexEdgeBasis(void){
- delete basis;
+
+ // Set Basis //
+ this->basis = new std::vector<const std::vector<Polynomial>*>
+ (basis.begin(), basis.end());
}
-/*
-#include <cstdio>
-int main(void){
- const int P = 1;
- const double d = 0.1;
- const char x[3] = {'X', 'Y', 'Z'};
-
- HexEdgeBasis b(P);
-
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<std::vector<Polynomial> >& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("close all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function [rx ry rz] = p(i, x, y, z)\n");
- printf("p = zeros(%d, 3);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++){
- for(int j = 0; j < 3; j++)
- printf("p(%d, %d) = %s;\n", i + 1, j + 1, basis[i][j].toString().c_str());
- //printf("p(%d) = %s", i, basis[i].toString().c_str());
- printf("\n");
- }
+HexEdgeBasis::~HexEdgeBasis(void){
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
- printf("\n");
- printf("rx = p(i, 1);\n");
- printf("ry = p(i, 2);\n");
- printf("rz = p(i, 3);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %f;\nx = [0:d:1];\ny = x;\nz = x;\n\nlx = length(x);\nly = length(y);\nlz = length(z);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- for(int j = 0; j < 3; j++)
- printf("p%d%c = zeros(lx, ly, lz);\n", i + 1, x[j]);
-
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly\n");
- printf("for k = 1:lz\n");
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("[p%dX(j, i, k), p%dY(j, i, k), p%dZ(j, i, k)] = p(%d, x(i), y(j), z(k));\n", i + 1, i + 1, i + 1, i + 1);
-
- printf("\n");
- printf("end\n");
- printf("end\n");
- printf("end\n");
-
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
-
- printf("\n");
- for(int i = b.getSize() - 1; i >= 0; i--)
- printf("figure;\nquiver3(x, y, z, p%dX, p%dY, p%dZ);\n", i + 1, i + 1, i + 1);
-
- printf("\n");
-
- return 0;
+ delete basis;
}
-*/
-
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 89f3331cec..fef2ee1cf4 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -69,49 +69,49 @@ HexNodeBasis::HexNodeBasis(const int order){
// Basis //
- basis = new std::vector<Polynomial>(size);
+ basis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) //
(*basis)[0] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[1] =
- (Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ new Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[2] =
- (Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ new Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[3] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+ new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[4] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1);
+ new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ Polynomial(1, 0, 0, 1));
(*basis)[5] =
- (Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1);
+ new Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ Polynomial(1, 0, 0, 1));
(*basis)[6] =
- (Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1);
+ new Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ Polynomial(1, 0, 0, 1));
(*basis)[7] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0)) *
- Polynomial(1, 0, 0, 1);
+ new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ Polynomial(1, 0, 0, 1));
// Edge Based //
@@ -124,9 +124,9 @@ HexNodeBasis::HexNodeBasis(const int order){
for(int l = 1; l < order; l++){
for(int e = 0; e < 12; e++){
- (*basis)[i] =
+ (*basis)[i] = new Polynomial(
legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) *
- ((*basis)[edge1[e]] + (*basis)[edge2[e]]);
+ (*(*basis)[edge1[e]] + *(*basis)[edge2[e]]));
i++;
}
@@ -151,20 +151,20 @@ HexNodeBasis::HexNodeBasis(const int order){
// 'Lambda' Functions
for(int f = 0; f < 6; f++)
lambda[f] =
- (*basis)[face1[f]] +
- (*basis)[face2[f]] +
- (*basis)[face3[f]] +
- (*basis)[face4[f]];
+ *(*basis)[face1[f]] +
+ *(*basis)[face2[f]] +
+ *(*basis)[face3[f]] +
+ *(*basis)[face4[f]];
// Face Based //
for(int l1 = 1; l1 < order; l1++){
for(int l2 = 1; l2 < order; l2++){
for(int f = 0; f < 6; f++){
- (*basis)[i] =
+ (*basis)[i] = new Polynomial(
legendre[l1].compose(xi[f]) *
legendre[l2].compose(eta[f]) *
- lambda[f];
+ lambda[f]);
i++;
}
@@ -185,9 +185,9 @@ HexNodeBasis::HexNodeBasis(const int order){
for(int l2 = 1; l2 < order; l2++){
for(int l3 = 1; l3 < order; l3++){
(*basis)[i] =
- legendre[l1].compose(px) *
- legendre[l2].compose(py) *
- legendre[l3].compose(pz);
+ new Polynomial(legendre[l1].compose(px) *
+ legendre[l2].compose(py) *
+ legendre[l3].compose(pz));
i++;
}
@@ -205,120 +205,8 @@ HexNodeBasis::HexNodeBasis(const int order){
}
HexNodeBasis::~HexNodeBasis(void){
- delete basis;
-}
-
-/*
-#include <cstdio>
-int main(void){
-
- const int P = 8;
- const double d = 0.05;
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
- HexNodeBasis b(P);
-
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<Polynomial>& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function r = p(i, x, y, z)\n");
- printf("p = zeros(%d, 1);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str());
-
- printf("\n");
- printf("r = p(i, 1);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %f;\nx = [0:d:1];\ny = x;\nz = x;\n\nlx = length(x);\nly = length(y);\nlz = length(z);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p%d = zeros(lx, ly, lz);\n", i + 1);
-
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly\n");
- printf("for k = 1:lz\n");
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p%d(j, i, k) = p(%d, x(i), y(j), z(k));\n", i + 1, i + 1);
-
- printf("\nend\n");
- printf("end\n");
- printf("end\n");
-
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
-
- printf("\n");
- for(int i = b.getSize(); i > 0; i--){
- printf("figure;\n");
-
- printf("subplot(3, 2, 1);\n");
- printf("contourf(x, y, squeeze(p%d(:, :, 1)));\n", i);
- printf("colorbar;\n");
- printf("title('z = 0');\n");
- printf("ylabel('x');\n");
- printf("xlabel('y');\n\n");
-
- printf("subplot(3, 2, 2);\n");
- printf("contourf(x, y, squeeze(p%d(:, :, end)));\n", i);
- printf("colorbar;\n");
- printf("title('z = 1');\n");
- printf("ylabel('x');\n");
- printf("xlabel('y');\n\n");
-
- printf("subplot(3, 2, 3);\n");
- printf("contourf(x, z, squeeze(p%d(:, 1, :)));\n", i);
- printf("colorbar;\n");
- printf("title('y = 0');\n");
- printf("ylabel('x');\n");
- printf("xlabel('z');\n\n");
-
- printf("subplot(3, 2, 4);\n");
- printf("contourf(x, z, squeeze(p%d(:, end, :)));\n", i);
- printf("colorbar;\n");
- printf("title('y = 1');\n");
- printf("ylabel('x');\n");
- printf("xlabel('z');\n\n");
-
- printf("subplot(3, 2, 5);\n");
- printf("contourf(y, z, squeeze(p%d(1, :, :)));\n", i);
- printf("colorbar;\n");
- printf("title('x = 0');\n");
- printf("ylabel('y');\n");
- printf("xlabel('z');\n\n");
-
- printf("subplot(3, 2, 6);\n");
- printf("contourf(y, z, squeeze(p%d(end, :, :)));\n", i);
- printf("colorbar;\n");
- printf("title('x = 1');\n");
- printf("ylabel('y');\n");
- printf("xlabel('z');\n\n\n");
-}
-
- printf("\n");
-
- return 0;
+ delete basis;
}
-*/
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 91d530f5e1..9d53d02e5f 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -78,11 +78,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
liftingSub[i] = lifting[j] - lifting[i];
- // Basis //
- basis = new std::vector<std::vector<Polynomial> >(size);
-
- for(int i = 0; i < size; i++)
- (*basis)[i].resize(3);
+ // Basis (temporary --- *no* const) //
+ std::vector<std::vector<Polynomial>*> basis(size);
// Edge Based (Nedelec) //
@@ -90,16 +87,15 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 4; e++){
- (*basis)[i] =
- (liftingSub[e]).gradient();
+ basis[i] = new std::vector<Polynomial>(liftingSub[e].gradient());
- (*basis)[i][0].mul(lagrangeSum[e]);
- (*basis)[i][1].mul(lagrangeSum[e]);
- (*basis)[i][2].mul(lagrangeSum[e]);
-
- (*basis)[i][0].mul(oneHalf);
- (*basis)[i][1].mul(oneHalf);
- (*basis)[i][2].mul(oneHalf);
+ basis[i]->at(0).mul(lagrangeSum[e]);
+ basis[i]->at(1).mul(lagrangeSum[e]);
+ basis[i]->at(2).mul(lagrangeSum[e]);
+
+ basis[i]->at(0).mul(oneHalf);
+ basis[i]->at(1).mul(oneHalf);
+ basis[i]->at(2).mul(oneHalf);
i++;
}
@@ -107,8 +103,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 4; e++){
- (*basis)[i] =
- (intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient();
+ basis[i] =
+ new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
i++;
}
@@ -133,7 +129,7 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 1) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
- (*basis)[i] = (iLegendreX[l1] * iLegendreY[l2]).gradient();
+ basis[i] = new std::vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
i++;
}
@@ -142,9 +138,11 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 2) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
- (*basis)[i][0] = legendreX[l1] * iLegendreY[l2];
- (*basis)[i][1] = iLegendreX[l1] * legendreY[l2] * -1;
- (*basis)[i][2] = zero;
+ basis[i] = new std::vector<Polynomial>(3);
+
+ basis[i]->at(0) = legendreX[l1] * iLegendreY[l2];
+ basis[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
+ basis[i]->at(2) = zero;
i++;
}
@@ -152,13 +150,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 3) //
for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){
- (*basis)[i][0] = iLegendreY[l];
- (*basis)[i][1] = zero;
- (*basis)[i][2] = zero;
+ basis[i] = new std::vector<Polynomial>(3);
+ basis[iPlus] = new std::vector<Polynomial>(3);
+
+ basis[i]->at(0) = iLegendreY[l];
+ basis[i]->at(1) = zero;
+ basis[i]->at(2) = zero;
- (*basis)[iPlus][0] = zero;
- (*basis)[iPlus][1] = iLegendreX[l];
- (*basis)[iPlus][2] = zero;
+ basis[iPlus]->at(0) = zero;
+ basis[iPlus]->at(1) = iLegendreX[l];
+ basis[iPlus]->at(2) = zero;
i++;
}
@@ -177,86 +178,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
delete[] lifting;
delete[] liftingSub;
-}
-QuadEdgeBasis::~QuadEdgeBasis(void){
- delete basis;
-}
-
-/*
-#include <cstdio>
-int main(void){
- const int P = 8;
- const double d = 0.05;
- const char x[2] = {'X', 'Y'};
- QuadEdgeBasis b(P);
-
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<std::vector<Polynomial> >& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("close all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function [rx ry] = p(i, x, y)\n");
- printf("p = zeros(%d, 2);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++){
- for(int j = 0; j < 2; j++)
- printf("p(%d, %d) = %s;\n", i + 1, j + 1, basis[i][j].toString().c_str());
- //printf("p(%d) = %s", i, basis[i].toString().c_str());
- printf("\n");
- }
-
- printf("\n");
- printf("rx = p(i, 1);\n");
- printf("ry = p(i, 2);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- for(int j = 0; j < 2; j++)
- printf("p%d%c = zeros(lx, ly);\n", i + 1, x[j]);
-
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly\n");
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("[p%dX(j, i), p%dY(j, i)] = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
-
- printf("\n");
- printf("end\n");
- printf("end\n");
-
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
+ // Set Basis //
+ this->basis = new std::vector<const std::vector<Polynomial>*>
+ (basis.begin(), basis.end());
+}
- printf("\n");
- for(int i = 0; i < b.getSize(); i++)
- printf("figure;\nquiver(x, y, p%dX, p%dY);\n", i + 1, i + 1);
-
- printf("\n");
+QuadEdgeBasis::~QuadEdgeBasis(void){
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
- return 0;
+ delete basis;
}
-*/
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 281bc6624f..861ff3eae0 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -42,24 +42,24 @@ QuadNodeBasis::QuadNodeBasis(const int order){
// Basis //
- basis = new std::vector<Polynomial>(size);
+ basis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) //
(*basis)[0] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+ new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
(*basis)[1] =
- (Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+ new Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
(*basis)[2] =
- (Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0));
+ new Polynomial((Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)));
(*basis)[3] =
- (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
- (Polynomial(1, 0, 1, 0));
+ new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)));
// Edge Based //
int i = 4;
@@ -67,7 +67,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
for(int l = 1; l < order; l++){
for(int e1 = 0, e2 = 1; e1 < 4; e1++, e2 = (e2 + 1) % 4){
(*basis)[i] =
- legendre[l].compose(lifting[e2] - lifting[e1]) * ((*basis)[e1] + (*basis)[e2]);
+ new Polynomial(legendre[l].compose(lifting[e2] - lifting[e1]) *
+ (*(*basis)[e1] + *(*basis)[e2]));
i++;
}
@@ -82,7 +83,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
for(int l1 = 1; l1 < order; l1++){
for(int l2 = 1; l2 < order; l2++){
- (*basis)[i] = legendre[l1].compose(px) * legendre[l2].compose(py);
+ (*basis)[i] =
+ new Polynomial(legendre[l1].compose(px) * legendre[l2].compose(py));
i++;
}
@@ -94,74 +96,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
}
QuadNodeBasis::~QuadNodeBasis(void){
- delete basis;
-}
-
-/*
-#include <cstdio>
-int main(void){
- const int P = 8;
- const double d = 0.05;
-
- QuadNodeBasis b(P);
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<Polynomial>& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function r = p(i, x, y)\n");
- printf("p = zeros(%d, 1);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str());
-
- printf("\n");
- printf("r = p(i, 1);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p%d = zeros(lx, ly);\n", i + 1);
-
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly\n");
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p%d(j, i) = p(%d, x(i), y(j));\n", i + 1, i + 1);
-
- printf("end\n");
- printf("end\n");
-
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
-
- printf("\n");
- for(int i = b.getSize(); i > 0; i--)
- printf("figure;\ncontourf(x, y, p%d);\ncolorbar;\n", i);
-
- printf("\n");
-
- return 0;
+ delete basis;
}
-*/
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index aa1492858f..8b70bc1fd2 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -52,11 +52,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
lagrangeSub[i] = lagrange[i] - lagrange[j];
- // Basis //
- basis = new std::vector<std::vector<Polynomial> >(size);
-
- for(int i = 0; i < size; i++)
- (*basis)[i].resize(3);
+ // Basis (temporary --- *no* const) //
+ std::vector<std::vector<Polynomial>*> basis(size);
// Edge Based (Nedelec) //
@@ -68,15 +65,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){
tmp[1].mul(lagrange[i]);
tmp[2].mul(lagrange[i]);
- (*basis)[i] = lagrange[i].gradient();
+ basis[i] = new std::vector<Polynomial>(lagrange[i].gradient());
- (*basis)[i][0].mul(lagrange[j]);
- (*basis)[i][1].mul(lagrange[j]);
- (*basis)[i][2].mul(lagrange[j]);
+ basis[i]->at(0).mul(lagrange[j]);
+ basis[i]->at(1).mul(lagrange[j]);
+ basis[i]->at(2).mul(lagrange[j]);
- (*basis)[i][0].sub(tmp[0]);
- (*basis)[i][1].sub(tmp[1]);
- (*basis)[i][2].sub(tmp[2]);
+ basis[i]->at(0).sub(tmp[0]);
+ basis[i]->at(1).sub(tmp[1]);
+ basis[i]->at(2).sub(tmp[2]);
i++;
}
@@ -84,8 +81,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 3; e++){
- (*basis)[i] =
- (intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient();
+ basis[i] =
+ new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient());
i++;
}
@@ -108,15 +105,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){
tmp[1].mul(u[l1]);
tmp[2].mul(u[l1]);
- (*basis)[i] = u[l1].gradient();
+ basis[i] = new std::vector<Polynomial>(u[l1].gradient());
- (*basis)[i][0].mul(v[l2]);
- (*basis)[i][1].mul(v[l2]);
- (*basis)[i][2].mul(v[l2]);
+ basis[i]->at(0).mul(v[l2]);
+ basis[i]->at(1).mul(v[l2]);
+ basis[i]->at(2).mul(v[l2]);
- (*basis)[i][0].add(tmp[0]);
- (*basis)[i][1].add(tmp[1]);
- (*basis)[i][2].add(tmp[2]);
+ basis[i]->at(0).add(tmp[0]);
+ basis[i]->at(1).add(tmp[1]);
+ basis[i]->at(2).add(tmp[2]);
i++;
}
@@ -130,15 +127,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){
tmp[1].mul(u[l1]);
tmp[2].mul(u[l1]);
- (*basis)[i] = u[l1].gradient();
+ basis[i] = new std::vector<Polynomial>(u[l1].gradient());
- (*basis)[i][0].mul(v[l2]);
- (*basis)[i][1].mul(v[l2]);
- (*basis)[i][2].mul(v[l2]);
+ basis[i]->at(0).mul(v[l2]);
+ basis[i]->at(1).mul(v[l2]);
+ basis[i]->at(2).mul(v[l2]);
- (*basis)[i][0].sub(tmp[0]);
- (*basis)[i][1].sub(tmp[1]);
- (*basis)[i][2].sub(tmp[2]);
+ basis[i]->at(0).sub(tmp[0]);
+ basis[i]->at(1).sub(tmp[1]);
+ basis[i]->at(2).sub(tmp[2]);
i++;
}
@@ -146,11 +143,11 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Cell Based (Type 3) //
for(int l = 0; l < orderMinus; l++){
- (*basis)[i] = (*basis)[0];
+ basis[i] = new std::vector<Polynomial>(*basis[0]);
- (*basis)[i][0].mul(v[l]);
- (*basis)[i][1].mul(v[l]);
- (*basis)[i][2].mul(v[l]);
+ basis[i]->at(0).mul(v[l]);
+ basis[i]->at(1).mul(v[l]);
+ basis[i]->at(2).mul(v[l]);
i++;
}
@@ -163,84 +160,16 @@ TriEdgeBasis::TriEdgeBasis(const int order){
delete[] lagrangeSum;
delete[] u;
delete[] v;
-}
-
-TriEdgeBasis::~TriEdgeBasis(void){
- delete basis;
-}
-
-/*
-#include <cstdio>
-int main(void){
- const int P = 8;
- const double d = 0.05;
- const char x[2] = {'X', 'Y'};
-
- TriEdgeBasis b(P);
-
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<std::vector<Polynomial> >& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function [rx ry] = p(i, x, y)\n");
- printf("p = zeros(%d, 2);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++){
- //printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str());
- printf("p(%d, 1) = %s;\n", i + 1, basis[i][0].toString().c_str());
- printf("p(%d, 2) = %s;\n", i + 1, basis[i][1].toString().c_str());
- printf("\n");
- }
-
- printf("\n");
- printf("rx = p(i, 1);\n");
- printf("ry = p(i, 2);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- for(int j = 0; j < 2; j++)
- printf("p%d%c = zeros(lx, ly);\n", i + 1, x[j]);
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly - i + 1\n");
- printf("\n");
- for(int i = 0; i < b.getSize(); i++)
- printf("[p%dX(j, i), p%dY(j, i)] = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
-
- printf("end\n");
- printf("end\n");
+ // Set Basis //
+ this->basis = new std::vector<const std::vector<Polynomial>*>
+ (basis.begin(), basis.end());
+}
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
+TriEdgeBasis::~TriEdgeBasis(void){
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
- printf("\n");
- for(int i = b.getSize() - 1; i >= 0; i--)
- printf("figure;\nquiver(x, y, p%dX, p%dY);\n", i + 1, i + 1);
-
- printf("\n");
-
- return 0;
+ delete basis;
}
-*/
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index 0e0beb05de..7a81191200 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -26,11 +26,10 @@ TriNedelecBasis::TriNedelecBasis(void){
lagrange[2] =
Polynomial(1, 0, 1, 0);
- // Basis //
- basis = new std::vector<std::vector<Polynomial> >(size);
- for(int i = 0; i < size; i++)
- (*basis)[i].resize(3);
+ // Basis (temporary --- *no* const) //
+ std::vector<std::vector<Polynomial>*> basis(size);
+
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
std::vector<Polynomial> tmp = lagrange[j].gradient();
@@ -38,96 +37,29 @@ TriNedelecBasis::TriNedelecBasis(void){
tmp[1].mul(lagrange[i]);
tmp[2].mul(lagrange[i]);
- (*basis)[i] = lagrange[i].gradient();
+ basis[i] = new std::vector<Polynomial>(lagrange[i].gradient());
- (*basis)[i][0].mul(lagrange[j]);
- (*basis)[i][1].mul(lagrange[j]);
- (*basis)[i][2].mul(lagrange[j]);
+ basis[i]->at(0).mul(lagrange[j]);
+ basis[i]->at(1).mul(lagrange[j]);
+ basis[i]->at(2).mul(lagrange[j]);
- (*basis)[i][0].sub(tmp[0]);
- (*basis)[i][1].sub(tmp[1]);
- (*basis)[i][2].sub(tmp[2]);
+ basis[i]->at(0).sub(tmp[0]);
+ basis[i]->at(1).sub(tmp[1]);
+ basis[i]->at(2).sub(tmp[2]);
}
// Free Temporary Sapce //
delete[] lagrange;
+
+
+ // Set Basis //
+ this->basis = new std::vector<const std::vector<Polynomial>*>
+ (basis.begin(), basis.end());
}
TriNedelecBasis::~TriNedelecBasis(void){
- delete basis;
-}
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
-/*
-#include <cstdio>
-int main(void){
- const double d = 0.05;
- const char x[2] = {'X', 'Y'};
-
- TriNedelecBasis b;
-
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<std::vector<Polynomial> >& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function [rx ry] = p(i, x, y)\n");
- printf("p = zeros(%d, 2);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++){
- //printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str());
- printf("p(%d, 1) = %s;\n", i + 1, basis[i][0].toString().c_str());
- printf("p(%d, 2) = %s;\n", i + 1, basis[i][1].toString().c_str());
- printf("\n");
- }
-
- printf("\n");
- printf("rx = p(i, 1);\n");
- printf("ry = p(i, 2);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- for(int j = 0; j < 2; j++)
- printf("p%d%c = zeros(lx, ly);\n", i + 1, x[j]);
-
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly - i + 1\n");
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("[p%dX(j, i), p%dY(j, i)] = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
-
- printf("end\n");
- printf("end\n");
-
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
-
- printf("\n");
- for(int i = b.getSize() - 1; i >= 0; i--)
- printf("figure;\nquiver(x, y, p%dX, p%dY);\n", i + 1, i + 1);
-
- printf("\n");
-
- return 0;
+ delete basis;
}
-*/
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index d034614cbc..362bb46421 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -29,26 +29,28 @@ TriNodeBasis::TriNodeBasis(const int order){
// Basis //
- basis = new std::vector<Polynomial>(size);
+ basis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) //
(*basis)[0] =
- Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
+ new Polynomial(Polynomial(1, 0, 0, 0) -
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0));
(*basis)[1] =
- Polynomial(1, 1, 0, 0);
+ new Polynomial(Polynomial(1, 1, 0, 0));
(*basis)[2] =
- Polynomial(1, 0, 1, 0);
+ new Polynomial(Polynomial(1, 0, 1, 0));
// Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = (*basis)[i] + (*basis)[j];
+ lagrangeSum[i] = *(*basis)[i] + *(*basis)[j];
// Lagrange Sub //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSub[i] = (*basis)[j] - (*basis)[i];
+ lagrangeSub[i] = *(*basis)[j] - *(*basis)[i];
// Edge Based //
@@ -56,22 +58,22 @@ TriNodeBasis::TriNodeBasis(const int order){
for(int l = 1; l < order; l++){
for(int e = 0; e < 3; e++){
- (*basis)[i] =
- intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]);
+ (*basis)[i] = new Polynomial(
+ intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]));
i++;
}
}
// Cell Based //
- Polynomial p = (*basis)[2] * 2 - Polynomial(1, 0, 0, 0);
+ Polynomial p = *(*basis)[2] * 2 - Polynomial(1, 0, 0, 0);
const int orderMinusTwo = order - 2;
for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
- (*basis)[i] =
+ (*basis)[i] = new Polynomial(
intLegendre[l1].compose(lagrangeSub[0], lagrangeSum[0]) *
- legendre[l2].compose(p) * (*basis)[2];
+ legendre[l2].compose(p) * *(*basis)[2]);
i++;
}
@@ -85,74 +87,8 @@ TriNodeBasis::TriNodeBasis(const int order){
}
TriNodeBasis::~TriNodeBasis(void){
- delete basis;
-}
-
-/*
-#include <cstdio>
-int main(void){
- const int P = 8;
- const double d = 0.01;
-
- TriNodeBasis b(P);
-
- printf("%d = %d + %d + %d + %d = %d\n",
- b.getSize(),
- b.getNVertex(), b.getNEdge(), b.getNFace(), b.getNCell(),
- b.getNVertex() + b.getNEdge() + b.getNFace() + b.getNCell());
-
- const std::vector<Polynomial>& basis = b.getBasis();
-
- printf("\n");
- printf("clear all;\n");
- printf("\n");
-
- printf("\n");
- printf("Order = %d\n", b.getOrder());
- printf("Type = %d\n", b.getType());
- printf("Size = %d\n", b.getSize());
- printf("NodeNumber = %d\n", b.getNodeNbr());
- printf("Dimension = %d\n", b.getDim());
- printf("\n");
-
- printf("function r = p(i, x, y)\n");
- printf("p = zeros(%d, 1);\n", b.getSize());
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str());
-
- printf("\n");
- printf("r = p(i, 1);\n");
- printf("end\n");
- printf("\n");
-
- printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p%d = zeros(lx, ly);\n", i + 1);
-
- printf("\n");
- printf("for i = 1:lx\n");
- printf("for j = 1:ly - i + 1\n");
- printf("\n");
-
- for(int i = 0; i < b.getSize(); i++)
- printf("p%d(j, i) = p(%d, x(i), y(j));\n", i + 1, i + 1);
-
- printf("end\n");
- printf("end\n");
+ for(int i = 0; i < size; i++)
+ delete (*basis)[i];
- printf("\n");
- printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
- printf("\n");
-
- printf("\n");
- for(int i = 0; i < b.getSize(); i++)
- printf("figure;\ncontourf(x, y, p%d);\ncolorbar;\n", i + 1);
-
- printf("\n");
-
- return 0;
+ delete basis;
}
-*/
--
GitLab