Skip to content
Snippets Groups Projects
Commit 4caf2d4d authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

Reimplementing Basis -- Needed for closure

parent 4f09edba
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 240 additions and 752 deletions
FunctionSpace/BasisScalar.h
+ 3
− 3
View file @ 4caf2d4d
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
class BasisScalar: public Basis{ class BasisScalar: public Basis{
protected: protected:
std::vector<Polynomial>* basis; std::vector<const Polynomial*>* basis;
public: public:
//! Deletes this BasisScalar //! Deletes this BasisScalar
...@@ -29,7 +29,7 @@ class BasisScalar: public Basis{ ...@@ -29,7 +29,7 @@ class BasisScalar: public Basis{
//! @return Returns the set of @em Polynomial%s //! @return Returns the set of @em Polynomial%s
//! defining this (scalar) Basis //! defining this (scalar) Basis
const std::vector<Polynomial>& getFunctions(void) const; const std::vector<const Polynomial*>& getFunctions(void) const;
protected: protected:
//! @internal //! @internal
...@@ -44,7 +44,7 @@ class BasisScalar: public Basis{ ...@@ -44,7 +44,7 @@ class BasisScalar: public Basis{
////////////////////// //////////////////////
inline inline
const std::vector<Polynomial>& BasisScalar:: const std::vector<const Polynomial*>& BasisScalar::
getFunctions(void) const{ getFunctions(void) const{
return *basis; return *basis;
......
FunctionSpace/BasisTest.cpp
+ 3
− 1
View file @ 4caf2d4d
...@@ -32,7 +32,9 @@ int basisTest(int argc, char** argv){ ...@@ -32,7 +32,9 @@ int basisTest(int argc, char** argv){
writer.setDomain(goe.getAll()); writer.setDomain(goe.getAll());
// Plot Basis // // Plot Basis //
TriNodeBasis b(3); HexEdgeBasis b(3);
cout << "Size: " << b.getSize() << endl;
PlotBasis plot(b, goe, writer); PlotBasis plot(b, goe, writer);
plot.plot("basis"); plot.plot("basis");
......
FunctionSpace/BasisVector.cpp
+ 0
− 1
View file @ 4caf2d4d
...@@ -6,4 +6,3 @@ BasisVector::BasisVector(void){ ...@@ -6,4 +6,3 @@ BasisVector::BasisVector(void){
BasisVector::~BasisVector(void){ BasisVector::~BasisVector(void){
} }
FunctionSpace/BasisVector.h
+ 3
− 3
View file @ 4caf2d4d
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
class BasisVector: public Basis{ class BasisVector: public Basis{
protected: protected:
std::vector<std::vector<Polynomial> >* basis; std::vector<const std::vector<Polynomial>*>* basis;
public: public:
//! Deletes this BasisVector //! Deletes this BasisVector
...@@ -29,7 +29,7 @@ class BasisVector: public Basis{ ...@@ -29,7 +29,7 @@ class BasisVector: public Basis{
//! @return Returns the set of @em Polynomial%s //! @return Returns the set of @em Polynomial%s
//! defining this (vectorial) Basis //! defining this (vectorial) Basis
const std::vector<std::vector<Polynomial> >& getFunctions(void) const; const std::vector<const std::vector<Polynomial>*>& getFunctions(void) const;
protected: protected:
//! @internal //! @internal
...@@ -44,7 +44,7 @@ class BasisVector: public Basis{ ...@@ -44,7 +44,7 @@ class BasisVector: public Basis{
////////////////////// //////////////////////
inline inline
const std::vector<std::vector<Polynomial> >& BasisVector:: const std::vector<const std::vector<Polynomial>*>& BasisVector::
getFunctions(void) const{ getFunctions(void) const{
return *basis; return *basis;
......
FunctionSpace/FunctionSpaceEdge.cpp
+ 3
− 3
View file @ 4caf2d4d
...@@ -25,7 +25,7 @@ interpolate(const MElement& element, ...@@ -25,7 +25,7 @@ interpolate(const MElement& element,
const_cast<MElement&>(element); const_cast<MElement&>(element);
// Get Basis Functions // // Get Basis Functions //
const vector<vector<Polynomial> >& fun = getBasis(element).getFunctions(); const vector<const vector<Polynomial>*>& fun = getBasis(element).getFunctions();
const unsigned int nFun = fun.size(); const unsigned int nFun = fun.size();
// Get Reference coordinate // // Get Reference coordinate //
...@@ -47,7 +47,7 @@ interpolate(const MElement& element, ...@@ -47,7 +47,7 @@ interpolate(const MElement& element,
for(unsigned int i = 0; i < nFun; i++){ for(unsigned int i = 0; i < nFun; i++){
fullVector<double> vi = 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); invJac);
vi.scale(coef[i]); vi.scale(coef[i]);
......
FunctionSpace/FunctionSpaceNode.cpp
+ 3
− 3
View file @ 4caf2d4d
...@@ -25,7 +25,7 @@ interpolate(const MElement& element, ...@@ -25,7 +25,7 @@ interpolate(const MElement& element,
const_cast<MElement&>(element); const_cast<MElement&>(element);
// Get Basis Functions // // Get Basis Functions //
const vector<Polynomial>& fun = getBasis(element).getFunctions(); const vector<const Polynomial*>& fun = getBasis(element).getFunctions();
const unsigned int nFun = fun.size(); const unsigned int nFun = fun.size();
// Get Reference coordinate // // Get Reference coordinate //
...@@ -39,7 +39,7 @@ interpolate(const MElement& element, ...@@ -39,7 +39,7 @@ interpolate(const MElement& element,
for(unsigned int i = 0; i < nFun; i++) for(unsigned int i = 0; i < nFun; i++)
val += 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 Interpolated Value //
return val; return val;
......
FunctionSpace/HexEdgeBasis.cpp
+ 64
− 129
View file @ 4caf2d4d
...@@ -148,11 +148,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -148,11 +148,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
liftingSub[i] = lifting[edge1[i]] - lifting[edge2[i]]; liftingSub[i] = lifting[edge1[i]] - lifting[edge2[i]];
// Basis // // Basis (temporary --- *no* const) //
basis = new std::vector<std::vector<Polynomial> >(size); std::vector<std::vector<Polynomial>*> basis(size);
for(int i = 0; i < size; i++)
(*basis)[i].resize(3);
// Edge Based (Nedelec) // // Edge Based (Nedelec) //
...@@ -160,16 +157,16 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -160,16 +157,16 @@ HexEdgeBasis::HexEdgeBasis(const int order){
Polynomial oneHalf(0.5, 0, 0, 0); Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 12; e++){ for(int e = 0; e < 12; e++){
(*basis)[i] = basis[i] =
(liftingSub[e]).gradient(); new std::vector<Polynomial>((liftingSub[e]).gradient());
(*basis)[i][0].mul(lagrangeSum[e]); basis[i]->at(0).mul(lagrangeSum[e]);
(*basis)[i][1].mul(lagrangeSum[e]); basis[i]->at(1).mul(lagrangeSum[e]);
(*basis)[i][2].mul(lagrangeSum[e]); basis[i]->at(2).mul(lagrangeSum[e]);
(*basis)[i][0].mul(oneHalf); basis[i]->at(0).mul(oneHalf);
(*basis)[i][1].mul(oneHalf); basis[i]->at(1).mul(oneHalf);
(*basis)[i][2].mul(oneHalf); basis[i]->at(2).mul(oneHalf);
i++; i++;
} }
...@@ -178,8 +175,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -178,8 +175,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Edge Based (High Order) // // Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){ for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 12; e++){ for(int e = 0; e < 12; e++){
(*basis)[i] = basis[i] =
(intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient(); new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
i++; i++;
} }
...@@ -237,9 +234,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -237,9 +234,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ for(int l2 = 1; l2 < orderPlus; l2++){
for(int f = 0; f < 6; f++){ for(int f = 0; f < 6; f++){
(*basis)[i] = (iLegendreXi[l1][f] * basis[i] = new std::vector<Polynomial>((iLegendreXi[l1][f] *
iLegendreEta[l2][f] * iLegendreEta[l2][f] *
lambda[f]).gradient(); lambda[f]).gradient());
i++; i++;
} }
...@@ -251,6 +248,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -251,6 +248,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ for(int l2 = 1; l2 < orderPlus; l2++){
for(int f = 0; f < 6; f++){ for(int f = 0; f < 6; f++){
basis[i] = new std::vector<Polynomial>(3);
Polynomial tmp1 = Polynomial tmp1 =
legendreXi[l1][f] * legendreXi[l1][f] *
iLegendreEta[l2][f]; iLegendreEta[l2][f];
...@@ -269,9 +268,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -269,9 +268,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){
gr2[1].mul(tmp2); gr2[1].mul(tmp2);
gr2[2].mul(tmp2); gr2[2].mul(tmp2);
(*basis)[i][0] = (gr1[0] - gr2[0]) * lambda[f]; basis[i]->at(0) = (gr1[0] - gr2[0]) * lambda[f];
(*basis)[i][1] = (gr1[1] - gr2[1]) * lambda[f]; basis[i]->at(1) = (gr1[1] - gr2[1]) * lambda[f];
(*basis)[i][2] = (gr1[2] - gr2[2]) * lambda[f]; basis[i]->at(2) = (gr1[2] - gr2[2]) * lambda[f];
i++; i++;
} }
...@@ -284,11 +283,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -284,11 +283,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int f = 0; f < 6; f++){ for(int f = 0; f < 6; f++){
Polynomial tmp = iLegendreEta[l][f] * lambda[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]->at(0).mul(tmp);
(*basis)[i][1].mul(tmp); basis[i]->at(1).mul(tmp);
(*basis)[i][2].mul(tmp); basis[i]->at(2).mul(tmp);
i++; i++;
} }
...@@ -300,11 +299,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -300,11 +299,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int f = 0; f < 6; f++){ for(int f = 0; f < 6; f++){
Polynomial tmp = iLegendreXi[l][f] * lambda[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]->at(0).mul(tmp);
(*basis)[i][1].mul(tmp); basis[i]->at(1).mul(tmp);
(*basis)[i][2].mul(tmp); basis[i]->at(2).mul(tmp);
i++; i++;
} }
...@@ -335,9 +334,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -335,9 +334,9 @@ HexEdgeBasis::HexEdgeBasis(const int order){
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ for(int l2 = 1; l2 < orderPlus; l2++){
for(int l3 = 1; l3 < orderPlus; l3++){ for(int l3 = 1; l3 < orderPlus; l3++){
(*basis)[i] = (iLegendreX[l1] * basis[i] = new std::vector<Polynomial>((iLegendreX[l1] *
iLegendreY[l2] * iLegendreY[l2] *
iLegendreZ[l3]).gradient(); iLegendreZ[l3]).gradient());
i++; i++;
cellNumber++; cellNumber++;
...@@ -348,11 +347,13 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -348,11 +347,13 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 2 -- First Part) // // Cell Based (Type 2 -- First Part) //
for(int j = 0; j < cellNumber; j++){ for(int j = 0; j < cellNumber; j++){
basis[i] = new std::vector<Polynomial>(3);
int off = j + cellStart; int off = j + cellStart;
(*basis)[i][0] = (*basis)[off][0]; basis[i]->at(0) = basis[off]->at(0);
(*basis)[i][1] = (*basis)[off][1] * Polynomial(-1, 0, 0, 0); basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
(*basis)[i][2] = (*basis)[off][2]; basis[i]->at(2) = basis[off]->at(2);
i++; i++;
} }
...@@ -360,11 +361,13 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -360,11 +361,13 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 2 -- Second Part) // // Cell Based (Type 2 -- Second Part) //
for(int j = 0; j < cellNumber; j++){ for(int j = 0; j < cellNumber; j++){
basis[i] = new std::vector<Polynomial>(3);
int off = j + cellStart; int off = j + cellStart;
(*basis)[i][0] = (*basis)[off][0]; basis[i]->at(0) = basis[off]->at(0);
(*basis)[i][1] = (*basis)[off][1] * Polynomial(-1, 0, 0, 0); basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
(*basis)[i][2] = (*basis)[off][2] * Polynomial(-1, 0, 0, 0); basis[i]->at(2) = basis[off]->at(2) * Polynomial(-1, 0, 0, 0);
i++; i++;
} }
...@@ -373,9 +376,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -373,9 +376,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 3 -- First Part) // // Cell Based (Type 3 -- First Part) //
for(int l2 = 1; l2 < orderPlus; l2++){ for(int l2 = 1; l2 < orderPlus; l2++){
for(int l3 = 1; l3 < orderPlus; l3++){ for(int l3 = 1; l3 < orderPlus; l3++){
(*basis)[i][0] = iLegendreY[l2] * iLegendreZ[l3]; basis[i] = new std::vector<Polynomial>(3);
(*basis)[i][1] = zero;
(*basis)[i][2] = zero; basis[i]->at(0) = iLegendreY[l2] * iLegendreZ[l3];
basis[i]->at(1) = zero;
basis[i]->at(2) = zero;
i++; i++;
} }
...@@ -385,9 +390,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -385,9 +390,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 3 -- Second Part) // // Cell Based (Type 3 -- Second Part) //
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l3 = 1; l3 < orderPlus; l3++){ for(int l3 = 1; l3 < orderPlus; l3++){
(*basis)[i][0] = zero; basis[i] = new std::vector<Polynomial>(3);
(*basis)[i][1] = iLegendreX[l1] * iLegendreZ[l3];
(*basis)[i][2] = zero; basis[i]->at(0) = zero;
basis[i]->at(1) = iLegendreX[l1] * iLegendreZ[l3];
basis[i]->at(2) = zero;
i++; i++;
} }
...@@ -397,9 +404,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -397,9 +404,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
// Cell Based (Type 3 -- Thrid Part) // // Cell Based (Type 3 -- Thrid Part) //
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ for(int l2 = 1; l2 < orderPlus; l2++){
(*basis)[i][0] = zero; basis[i] = new std::vector<Polynomial>(3);
(*basis)[i][1] = zero;
(*basis)[i][2] = iLegendreX[l1] * iLegendreY[l2]; basis[i]->at(0) = zero;
basis[i]->at(1) = zero;
basis[i]->at(2) = iLegendreX[l1] * iLegendreY[l2];
i++; i++;
} }
...@@ -438,90 +447,16 @@ HexEdgeBasis::HexEdgeBasis(const int order){ ...@@ -438,90 +447,16 @@ HexEdgeBasis::HexEdgeBasis(const int order){
delete[] iLegendreX; delete[] iLegendreX;
delete[] iLegendreY; delete[] iLegendreY;
delete[] iLegendreZ; delete[] iLegendreZ;
}
HexEdgeBasis::~HexEdgeBasis(void){
delete basis;
}
/* // Set Basis //
#include <cstdio> this->basis = new std::vector<const std::vector<Polynomial>*>
int main(void){ (basis.begin(), basis.end());
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");
} }
printf("\n"); HexEdgeBasis::~HexEdgeBasis(void){
printf("rx = p(i, 1);\n"); for(int i = 0; i < size; i++)
printf("ry = p(i, 2);\n"); delete (*basis)[i];
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;
} }
*/
FunctionSpace/HexNodeBasis.cpp
+ 39
− 151
View file @ 4caf2d4d
...@@ -69,49 +69,49 @@ HexNodeBasis::HexNodeBasis(const int order){ ...@@ -69,49 +69,49 @@ HexNodeBasis::HexNodeBasis(const int order){
// Basis // // Basis //
basis = new std::vector<Polynomial>(size); basis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) // // Vertex Based (Lagrange) //
(*basis)[0] = (*basis)[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) * 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, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)); (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[1] = (*basis)[1] =
(Polynomial(1, 1, 0, 0)) * 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, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)); (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[2] = (*basis)[2] =
(Polynomial(1, 1, 0, 0)) * new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) * (Polynomial(1, 0, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)); (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[3] = (*basis)[3] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) * new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) * (Polynomial(1, 0, 1, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)); (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
(*basis)[4] = (*basis)[4] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) * 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, 1, 0)) *
Polynomial(1, 0, 0, 1); Polynomial(1, 0, 0, 1));
(*basis)[5] = (*basis)[5] =
(Polynomial(1, 1, 0, 0)) * 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, 1, 0)) *
Polynomial(1, 0, 0, 1); Polynomial(1, 0, 0, 1));
(*basis)[6] = (*basis)[6] =
(Polynomial(1, 1, 0, 0)) * new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) * (Polynomial(1, 0, 1, 0)) *
Polynomial(1, 0, 0, 1); Polynomial(1, 0, 0, 1));
(*basis)[7] = (*basis)[7] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) * new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)) * (Polynomial(1, 0, 1, 0)) *
Polynomial(1, 0, 0, 1); Polynomial(1, 0, 0, 1));
// Edge Based // // Edge Based //
...@@ -124,9 +124,9 @@ HexNodeBasis::HexNodeBasis(const int order){ ...@@ -124,9 +124,9 @@ HexNodeBasis::HexNodeBasis(const int order){
for(int l = 1; l < order; l++){ for(int l = 1; l < order; l++){
for(int e = 0; e < 12; e++){ for(int e = 0; e < 12; e++){
(*basis)[i] = (*basis)[i] = new Polynomial(
legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) * legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) *
((*basis)[edge1[e]] + (*basis)[edge2[e]]); (*(*basis)[edge1[e]] + *(*basis)[edge2[e]]));
i++; i++;
} }
...@@ -151,20 +151,20 @@ HexNodeBasis::HexNodeBasis(const int order){ ...@@ -151,20 +151,20 @@ HexNodeBasis::HexNodeBasis(const int order){
// 'Lambda' Functions // 'Lambda' Functions
for(int f = 0; f < 6; f++) for(int f = 0; f < 6; f++)
lambda[f] = lambda[f] =
(*basis)[face1[f]] + *(*basis)[face1[f]] +
(*basis)[face2[f]] + *(*basis)[face2[f]] +
(*basis)[face3[f]] + *(*basis)[face3[f]] +
(*basis)[face4[f]]; *(*basis)[face4[f]];
// Face Based // // Face Based //
for(int l1 = 1; l1 < order; l1++){ for(int l1 = 1; l1 < order; l1++){
for(int l2 = 1; l2 < order; l2++){ for(int l2 = 1; l2 < order; l2++){
for(int f = 0; f < 6; f++){ for(int f = 0; f < 6; f++){
(*basis)[i] = (*basis)[i] = new Polynomial(
legendre[l1].compose(xi[f]) * legendre[l1].compose(xi[f]) *
legendre[l2].compose(eta[f]) * legendre[l2].compose(eta[f]) *
lambda[f]; lambda[f]);
i++; i++;
} }
...@@ -185,9 +185,9 @@ HexNodeBasis::HexNodeBasis(const int order){ ...@@ -185,9 +185,9 @@ HexNodeBasis::HexNodeBasis(const int order){
for(int l2 = 1; l2 < order; l2++){ for(int l2 = 1; l2 < order; l2++){
for(int l3 = 1; l3 < order; l3++){ for(int l3 = 1; l3 < order; l3++){
(*basis)[i] = (*basis)[i] =
legendre[l1].compose(px) * new Polynomial(legendre[l1].compose(px) *
legendre[l2].compose(py) * legendre[l2].compose(py) *
legendre[l3].compose(pz); legendre[l3].compose(pz));
i++; i++;
} }
...@@ -205,120 +205,8 @@ HexNodeBasis::HexNodeBasis(const int order){ ...@@ -205,120 +205,8 @@ HexNodeBasis::HexNodeBasis(const int order){
} }
HexNodeBasis::~HexNodeBasis(void){ HexNodeBasis::~HexNodeBasis(void){
delete basis; for(int i = 0; i < size; i++)
} delete (*basis)[i];
/* delete basis;
#include <cstdio>
int main(void){
const int P = 8;
const double d = 0.05;
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;
} }
*/
FunctionSpace/QuadEdgeBasis.cpp
+ 35
− 104
View file @ 4caf2d4d
...@@ -78,11 +78,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -78,11 +78,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
liftingSub[i] = lifting[j] - lifting[i]; liftingSub[i] = lifting[j] - lifting[i];
// Basis // // Basis (temporary --- *no* const) //
basis = new std::vector<std::vector<Polynomial> >(size); std::vector<std::vector<Polynomial>*> basis(size);
for(int i = 0; i < size; i++)
(*basis)[i].resize(3);
// Edge Based (Nedelec) // // Edge Based (Nedelec) //
...@@ -90,16 +87,15 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -90,16 +87,15 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
Polynomial oneHalf(0.5, 0, 0, 0); Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 4; e++){ for(int e = 0; e < 4; e++){
(*basis)[i] = basis[i] = new std::vector<Polynomial>(liftingSub[e].gradient());
(liftingSub[e]).gradient();
(*basis)[i][0].mul(lagrangeSum[e]); basis[i]->at(0).mul(lagrangeSum[e]);
(*basis)[i][1].mul(lagrangeSum[e]); basis[i]->at(1).mul(lagrangeSum[e]);
(*basis)[i][2].mul(lagrangeSum[e]); basis[i]->at(2).mul(lagrangeSum[e]);
(*basis)[i][0].mul(oneHalf); basis[i]->at(0).mul(oneHalf);
(*basis)[i][1].mul(oneHalf); basis[i]->at(1).mul(oneHalf);
(*basis)[i][2].mul(oneHalf); basis[i]->at(2).mul(oneHalf);
i++; i++;
} }
...@@ -107,8 +103,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -107,8 +103,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Edge Based (High Order) // // Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){ for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 4; e++){ for(int e = 0; e < 4; e++){
(*basis)[i] = basis[i] =
(intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient(); new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
i++; i++;
} }
...@@ -133,7 +129,7 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -133,7 +129,7 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 1) // // Cell Based (Type 1) //
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ 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++; i++;
} }
...@@ -142,9 +138,11 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -142,9 +138,11 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 2) // // Cell Based (Type 2) //
for(int l1 = 1; l1 < orderPlus; l1++){ for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ for(int l2 = 1; l2 < orderPlus; l2++){
(*basis)[i][0] = legendreX[l1] * iLegendreY[l2]; basis[i] = new std::vector<Polynomial>(3);
(*basis)[i][1] = iLegendreX[l1] * legendreY[l2] * -1;
(*basis)[i][2] = zero; basis[i]->at(0) = legendreX[l1] * iLegendreY[l2];
basis[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
basis[i]->at(2) = zero;
i++; i++;
} }
...@@ -152,13 +150,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -152,13 +150,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 3) // // Cell Based (Type 3) //
for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){ for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){
(*basis)[i][0] = iLegendreY[l]; basis[i] = new std::vector<Polynomial>(3);
(*basis)[i][1] = zero; basis[iPlus] = new std::vector<Polynomial>(3);
(*basis)[i][2] = zero;
(*basis)[iPlus][0] = zero; basis[i]->at(0) = iLegendreY[l];
(*basis)[iPlus][1] = iLegendreX[l]; basis[i]->at(1) = zero;
(*basis)[iPlus][2] = zero; basis[i]->at(2) = zero;
basis[iPlus]->at(0) = zero;
basis[iPlus]->at(1) = iLegendreX[l];
basis[iPlus]->at(2) = zero;
i++; i++;
} }
...@@ -177,86 +178,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){ ...@@ -177,86 +178,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
delete[] lifting; delete[] lifting;
delete[] liftingSub; delete[] liftingSub;
}
QuadEdgeBasis::~QuadEdgeBasis(void){
delete basis;
}
/* // Set Basis //
#include <cstdio> this->basis = new std::vector<const std::vector<Polynomial>*>
int main(void){ (basis.begin(), basis.end());
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"); QuadEdgeBasis::~QuadEdgeBasis(void){
printf("rx = p(i, 1);\n"); for(int i = 0; i < size; i++)
printf("ry = p(i, 2);\n"); delete (*basis)[i];
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");
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");
return 0; delete basis;
} }
*/
FunctionSpace/QuadNodeBasis.cpp
+ 16
− 80
View file @ 4caf2d4d
...@@ -42,24 +42,24 @@ QuadNodeBasis::QuadNodeBasis(const int order){ ...@@ -42,24 +42,24 @@ QuadNodeBasis::QuadNodeBasis(const int order){
// Basis // // Basis //
basis = new std::vector<Polynomial>(size); basis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) // // Vertex Based (Lagrange) //
(*basis)[0] = (*basis)[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) * 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, 1, 0)));
(*basis)[1] = (*basis)[1] =
(Polynomial(1, 1, 0, 0)) * 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, 1, 0)));
(*basis)[2] = (*basis)[2] =
(Polynomial(1, 1, 0, 0)) * new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)); (Polynomial(1, 0, 1, 0)));
(*basis)[3] = (*basis)[3] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) * new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)); (Polynomial(1, 0, 1, 0)));
// Edge Based // // Edge Based //
int i = 4; int i = 4;
...@@ -67,7 +67,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){ ...@@ -67,7 +67,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
for(int l = 1; l < order; l++){ for(int l = 1; l < order; l++){
for(int e1 = 0, e2 = 1; e1 < 4; e1++, e2 = (e2 + 1) % 4){ for(int e1 = 0, e2 = 1; e1 < 4; e1++, e2 = (e2 + 1) % 4){
(*basis)[i] = (*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++; i++;
} }
...@@ -82,7 +83,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){ ...@@ -82,7 +83,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
for(int l1 = 1; l1 < order; l1++){ for(int l1 = 1; l1 < order; l1++){
for(int l2 = 1; l2 < order; l2++){ 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++; i++;
} }
...@@ -94,74 +96,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){ ...@@ -94,74 +96,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
} }
QuadNodeBasis::~QuadNodeBasis(void){ QuadNodeBasis::~QuadNodeBasis(void){
delete basis; for(int i = 0; i < size; i++)
} delete (*basis)[i];
/*
#include <cstdio>
int main(void){
const int P = 8;
const double d = 0.05;
QuadNodeBasis 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"); delete basis;
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;
} }
*/
FunctionSpace/TriEdgeBasis.cpp
+ 37
− 108
View file @ 4caf2d4d
...@@ -52,11 +52,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -52,11 +52,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
lagrangeSub[i] = lagrange[i] - lagrange[j]; lagrangeSub[i] = lagrange[i] - lagrange[j];
// Basis // // Basis (temporary --- *no* const) //
basis = new std::vector<std::vector<Polynomial> >(size); std::vector<std::vector<Polynomial>*> basis(size);
for(int i = 0; i < size; i++)
(*basis)[i].resize(3);
// Edge Based (Nedelec) // // Edge Based (Nedelec) //
...@@ -68,15 +65,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -68,15 +65,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){
tmp[1].mul(lagrange[i]); tmp[1].mul(lagrange[i]);
tmp[2].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]->at(0).mul(lagrange[j]);
(*basis)[i][1].mul(lagrange[j]); basis[i]->at(1).mul(lagrange[j]);
(*basis)[i][2].mul(lagrange[j]); basis[i]->at(2).mul(lagrange[j]);
(*basis)[i][0].sub(tmp[0]); basis[i]->at(0).sub(tmp[0]);
(*basis)[i][1].sub(tmp[1]); basis[i]->at(1).sub(tmp[1]);
(*basis)[i][2].sub(tmp[2]); basis[i]->at(2).sub(tmp[2]);
i++; i++;
} }
...@@ -84,8 +81,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -84,8 +81,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Edge Based (High Order) // // Edge Based (High Order) //
for(int l = 1; l < orderPlus; l++){ for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 3; e++){ for(int e = 0; e < 3; e++){
(*basis)[i] = basis[i] =
(intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient(); new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient());
i++; i++;
} }
...@@ -108,15 +105,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -108,15 +105,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){
tmp[1].mul(u[l1]); tmp[1].mul(u[l1]);
tmp[2].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]->at(0).mul(v[l2]);
(*basis)[i][1].mul(v[l2]); basis[i]->at(1).mul(v[l2]);
(*basis)[i][2].mul(v[l2]); basis[i]->at(2).mul(v[l2]);
(*basis)[i][0].add(tmp[0]); basis[i]->at(0).add(tmp[0]);
(*basis)[i][1].add(tmp[1]); basis[i]->at(1).add(tmp[1]);
(*basis)[i][2].add(tmp[2]); basis[i]->at(2).add(tmp[2]);
i++; i++;
} }
...@@ -130,15 +127,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -130,15 +127,15 @@ TriEdgeBasis::TriEdgeBasis(const int order){
tmp[1].mul(u[l1]); tmp[1].mul(u[l1]);
tmp[2].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]->at(0).mul(v[l2]);
(*basis)[i][1].mul(v[l2]); basis[i]->at(1).mul(v[l2]);
(*basis)[i][2].mul(v[l2]); basis[i]->at(2).mul(v[l2]);
(*basis)[i][0].sub(tmp[0]); basis[i]->at(0).sub(tmp[0]);
(*basis)[i][1].sub(tmp[1]); basis[i]->at(1).sub(tmp[1]);
(*basis)[i][2].sub(tmp[2]); basis[i]->at(2).sub(tmp[2]);
i++; i++;
} }
...@@ -146,11 +143,11 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -146,11 +143,11 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Cell Based (Type 3) // // Cell Based (Type 3) //
for(int l = 0; l < orderMinus; l++){ 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]->at(0).mul(v[l]);
(*basis)[i][1].mul(v[l]); basis[i]->at(1).mul(v[l]);
(*basis)[i][2].mul(v[l]); basis[i]->at(2).mul(v[l]);
i++; i++;
} }
...@@ -163,84 +160,16 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -163,84 +160,16 @@ TriEdgeBasis::TriEdgeBasis(const int order){
delete[] lagrangeSum; delete[] lagrangeSum;
delete[] u; delete[] u;
delete[] v; delete[] v;
}
TriEdgeBasis::~TriEdgeBasis(void){
delete basis;
}
/* // Set Basis //
#include <cstdio> this->basis = new std::vector<const std::vector<Polynomial>*>
int main(void){ (basis.begin(), basis.end());
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"); TriEdgeBasis::~TriEdgeBasis(void){
printf("rx = p(i, 1);\n"); for(int i = 0; i < size; i++)
printf("ry = p(i, 2);\n"); delete (*basis)[i];
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;
} }
*/
FunctionSpace/TriNedelecBasis.cpp
+ 18
− 86
View file @ 4caf2d4d
...@@ -26,11 +26,10 @@ TriNedelecBasis::TriNedelecBasis(void){ ...@@ -26,11 +26,10 @@ TriNedelecBasis::TriNedelecBasis(void){
lagrange[2] = lagrange[2] =
Polynomial(1, 0, 1, 0); Polynomial(1, 0, 1, 0);
// Basis //
basis = new std::vector<std::vector<Polynomial> >(size);
for(int i = 0; i < size; i++) // Basis (temporary --- *no* const) //
(*basis)[i].resize(3); std::vector<std::vector<Polynomial>*> basis(size);
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){ for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
std::vector<Polynomial> tmp = lagrange[j].gradient(); std::vector<Polynomial> tmp = lagrange[j].gradient();
...@@ -38,96 +37,29 @@ TriNedelecBasis::TriNedelecBasis(void){ ...@@ -38,96 +37,29 @@ TriNedelecBasis::TriNedelecBasis(void){
tmp[1].mul(lagrange[i]); tmp[1].mul(lagrange[i]);
tmp[2].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]->at(0).mul(lagrange[j]);
(*basis)[i][1].mul(lagrange[j]); basis[i]->at(1).mul(lagrange[j]);
(*basis)[i][2].mul(lagrange[j]); basis[i]->at(2).mul(lagrange[j]);
(*basis)[i][0].sub(tmp[0]); basis[i]->at(0).sub(tmp[0]);
(*basis)[i][1].sub(tmp[1]); basis[i]->at(1).sub(tmp[1]);
(*basis)[i][2].sub(tmp[2]); basis[i]->at(2).sub(tmp[2]);
} }
// Free Temporary Sapce // // Free Temporary Sapce //
delete[] lagrange; delete[] lagrange;
}
TriNedelecBasis::~TriNedelecBasis(void){
delete basis;
}
/* // Set Basis //
#include <cstdio> this->basis = new std::vector<const std::vector<Polynomial>*>
int main(void){ (basis.begin(), basis.end());
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"); TriNedelecBasis::~TriNedelecBasis(void){
printf("rx = p(i, 1);\n"); for(int i = 0; i < size; i++)
printf("ry = p(i, 2);\n"); delete (*basis)[i];
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;
} }
*/
FunctionSpace/TriNodeBasis.cpp
+ 16
− 80
View file @ 4caf2d4d
...@@ -29,26 +29,28 @@ TriNodeBasis::TriNodeBasis(const int order){ ...@@ -29,26 +29,28 @@ TriNodeBasis::TriNodeBasis(const int order){
// Basis // // Basis //
basis = new std::vector<Polynomial>(size); basis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) // // Vertex Based (Lagrange) //
(*basis)[0] = (*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] = (*basis)[1] =
Polynomial(1, 1, 0, 0); new Polynomial(Polynomial(1, 1, 0, 0));
(*basis)[2] = (*basis)[2] =
Polynomial(1, 0, 1, 0); new Polynomial(Polynomial(1, 0, 1, 0));
// Lagrange Sum // // Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3) 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 // // Lagrange Sub //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3) 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 // // Edge Based //
...@@ -56,22 +58,22 @@ TriNodeBasis::TriNodeBasis(const int order){ ...@@ -56,22 +58,22 @@ TriNodeBasis::TriNodeBasis(const int order){
for(int l = 1; l < order; l++){ for(int l = 1; l < order; l++){
for(int e = 0; e < 3; e++){ for(int e = 0; e < 3; e++){
(*basis)[i] = (*basis)[i] = new Polynomial(
intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]); intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]));
i++; i++;
} }
} }
// Cell Based // // 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; const int orderMinusTwo = order - 2;
for(int l1 = 1; l1 < orderMinus; l1++){ for(int l1 = 1; l1 < orderMinus; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){ for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
(*basis)[i] = (*basis)[i] = new Polynomial(
intLegendre[l1].compose(lagrangeSub[0], lagrangeSum[0]) * intLegendre[l1].compose(lagrangeSub[0], lagrangeSum[0]) *
legendre[l2].compose(p) * (*basis)[2]; legendre[l2].compose(p) * *(*basis)[2]);
i++; i++;
} }
...@@ -85,74 +87,8 @@ TriNodeBasis::TriNodeBasis(const int order){ ...@@ -85,74 +87,8 @@ TriNodeBasis::TriNodeBasis(const int order){
} }
TriNodeBasis::~TriNodeBasis(void){ TriNodeBasis::~TriNodeBasis(void){
delete basis; for(int i = 0; i < size; i++)
} delete (*basis)[i];
/*
#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"); delete basis;
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");
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;
} }
*/
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment