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

QuadEdgeBasis: OK -- QuadBasis Convergence: OK \o/ -- Need a sexy way to...

QuadEdgeBasis: OK -- QuadBasis Convergence: OK \o/ -- Need a sexy way to handle the fact that QuadBasis are BILinear: need more integration points
parent 2f4b3de6
No related branches found
No related tags found
No related merge requests found
...@@ -14,6 +14,7 @@ ...@@ -14,6 +14,7 @@
#include "QuadNodeBasis.h" #include "QuadNodeBasis.h"
#include "QuadEdgeBasis.h" #include "QuadEdgeBasis.h"
#include "QuadNedelecBasis.h"
#include "QuadLagrangeBasis.h" #include "QuadLagrangeBasis.h"
#include "TetNodeBasis.h" #include "TetNodeBasis.h"
...@@ -116,8 +117,8 @@ BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType, ...@@ -116,8 +117,8 @@ BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
switch(basisType){ switch(basisType){
case 0: return new QuadNodeBasis(order); case 0: return new QuadNodeBasis(order);
case 1: case 1:
if (order == 0) throw Exception("Nedelec not implemented on Quads"); if (order == 0) return new QuadNedelecBasis();
else throw Exception("1-form not implemented on Quads"); else return new QuadEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Quads"); case 2: throw Exception("2-form not implemented on Quads");
case 3: throw Exception("3-form not implemented on Quads"); case 3: throw Exception("3-form not implemented on Quads");
......
...@@ -32,7 +32,8 @@ set(SRC ...@@ -32,7 +32,8 @@ set(SRC
TriLagrangeBasis.cpp TriLagrangeBasis.cpp
QuadNodeBasis.cpp QuadNodeBasis.cpp
# QuadEdgeBasis.cpp QuadEdgeBasis.cpp
QuadNedelecBasis.cpp
QuadLagrangeBasis.cpp QuadLagrangeBasis.cpp
TetNodeBasis.cpp TetNodeBasis.cpp
......
#include "QuadEdgeBasis.h" #include "QuadEdgeBasis.h"
#include "QuadReferenceSpace.h"
#include "Legendre.h" #include "Legendre.h"
using namespace std; using namespace std;
QuadEdgeBasis::QuadEdgeBasis(int order){ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNPermutation();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
const vector<const vector<const vector<unsigned int>*>*>&
faceV = refSpace->getAllFace();
// Set Basis Type // // Set Basis Type //
this->order = order; this->order = order;
...@@ -14,192 +25,187 @@ QuadEdgeBasis::QuadEdgeBasis(int order){ ...@@ -14,192 +25,187 @@ QuadEdgeBasis::QuadEdgeBasis(int order){
nEdge = 4 * (order + 1); nEdge = 4 * (order + 1);
nFace = 2 * (order + 1) * order; nFace = 2 * (order + 1) * order;
nCell = 0; nCell = 0;
nFunction = nVertex + nEdge + nFace + nCell;
nEdgeClosure = 2; // Legendre Polynomial //
nFaceClosure = 0; const unsigned int orderPlus = order + 1;
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
Polynomial* legendre = new Polynomial[orderPlus]; Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus]; Polynomial* intLegendre = new Polynomial[orderPlus];
Polynomial* iLegendreX = new Polynomial[orderPlus];
Polynomial* iLegendreY = new Polynomial[orderPlus];
Polynomial* legendreX = new Polynomial[orderPlus];
Polynomial* legendreY = new Polynomial[orderPlus];
Polynomial lagrange[4];
Polynomial lifting[4];
Polynomial lagrangeSum[4];
Polynomial liftingSub[2][4];
// Legendre Polynomial //
Legendre::integrated(intLegendre, orderPlus); Legendre::integrated(intLegendre, orderPlus);
Legendre::legendre(legendre, order); Legendre::legendre(legendre, order);
// Vertices definig Edges & Permutations // // Lagrange & Lifting //
const int edgeV[2][4][2] = const Polynomial lagrange[4] =
{ {
{ {0, 1}, {1, 2}, {2, 3}, {3, 0} }, Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
{ {1, 0}, {2, 1}, {3, 2}, {0, 3} } (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)))
};
const Polynomial lifting[4] =
{
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0)))
}; };
// Basis // // Basis //
node = new vector<vector<Polynomial>*>(nVertex); basis = new vector<Polynomial>**[nRefSpace];
edge = new vector<vector<vector<Polynomial>*>*>(2);
face = new vector<vector<vector<Polynomial>*>*>(0); for(unsigned int s = 0; s < nRefSpace; s++)
cell = new vector<vector<Polynomial>*>(nCell); basis[s] = new vector<Polynomial>*[nFunction];
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge); // Edge Based //
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge); for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = 0;
// Lagrange //
lagrange[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
lagrange[1] =
(Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
lagrange[2] =
(Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0));
lagrange[3] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0));
// Lagrange Sum //
for(int e = 0; e < 4; e++)
lagrangeSum[e] =
lagrange[edgeV[0][e][0]] +
lagrange[edgeV[0][e][1]];
// Lifting //
lifting[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
lifting[1] =
(Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
lifting[2] =
(Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0));
lifting[3] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0));
// Lifting Sub //
for(int e = 0; e < 4; e++){
liftingSub[0][e] =
lifting[edgeV[0][e][0]] -
lifting[edgeV[0][e][1]];
liftingSub[1][e] =
lifting[edgeV[1][e][0]] -
lifting[edgeV[1][e][1]];
}
for(unsigned int e = 0; e < 4; e++){
for(unsigned int l = 0; l < orderPlus; l++){
// Nedelec
if(l == 0){
Polynomial lambda = (lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]]) * 0.5;
// Edge Based (Nedelec) // basis[s][i] =
Polynomial oneHalf(0.5, 0, 0, 0); new vector<Polynomial>((lifting[(*(*edgeV[s])[e])[1]] -
lifting[(*(*edgeV[s])[e])[0]]).gradient());
for(int c = 0; c < 2; c++){ basis[s][i]->at(0).mul(lambda);
for(int e = 0; e < 4; e++){ basis[s][i]->at(1).mul(lambda);
(*(*edge)[c])[e] = basis[s][i]->at(2).mul(lambda);
new vector<Polynomial>(liftingSub[c][e].gradient()); }
(*(*edge)[c])[e]->at(0).mul(lagrangeSum[e]); // High Order
(*(*edge)[c])[e]->at(1).mul(lagrangeSum[e]); else{
(*(*edge)[c])[e]->at(2).mul(lagrangeSum[e]); basis[s][i] =
new vector<Polynomial>
((intLegendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
lifting[(*(*edgeV[s])[e])[0]])
*
(lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]])).gradient());
}
(*(*edge)[c])[e]->at(0).mul(oneHalf); i++;
(*(*edge)[c])[e]->at(1).mul(oneHalf); }
(*(*edge)[c])[e]->at(2).mul(oneHalf);
} }
} }
// Face Based //
// Edge Based (High Order) // // NB: We use (*(*faceV[s])[f])[]
for(int c = 0; c < 2; c++){ // where f = 0, because triangles
unsigned int i = 0; // have only ONE face: the face '0'
for(int l = 1; l < orderPlus; l++){ for(unsigned int s = 0; s < nRefSpace; s++){
for(int e = 0; e < 4; e++){ unsigned int i = nEdge;
(*(*edge)[c])[i + 4] =
new vector<Polynomial>((intLegendre[l].compose(liftingSub[c][e]) * // Type 1
lagrangeSum[e]).gradient()); for(unsigned int l1 = 1; l1 < orderPlus; l1++){
for(unsigned int 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[l2].compose(lifting[(*(*faceV[s])[0])[0]] -
lifting[(*(*faceV[s])[0])[3]])).gradient());
i++; 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]];
// Cell Based (Preliminary) // Polynomial pTwo = lifting[(*(*faceV[s])[0])[0]] -
Polynomial px = Polynomial(2, 1, 0, 0); lifting[(*(*faceV[s])[0])[3]];
Polynomial py = Polynomial(2, 0, 1, 0);
Polynomial zero = Polynomial(0, 0, 0, 0);
px = px - Polynomial(1, 0, 0, 0); Polynomial lOne = legendre[l1].compose(pOne) *
py = py - Polynomial(1, 0, 0, 0); intLegendre[l2].compose(pTwo);
unsigned int i = 0; Polynomial lTwo = legendre[l2].compose(pTwo) *
intLegendre[l1].compose(pOne);
for(int l = 0; l < orderPlus; l++){ vector<Polynomial> gradOne = pOne.gradient();
iLegendreX[l] = intLegendre[l].compose(px); vector<Polynomial> gradTwo = pTwo.gradient();
iLegendreY[l] = intLegendre[l].compose(py);
legendreX[l] = legendre[l].compose(px); gradOne[0].mul(lOne);
legendreY[l] = legendre[l].compose(py); gradOne[1].mul(lOne);
} gradOne[2].mul(lOne);
gradTwo[0].mul(lTwo);
gradTwo[1].mul(lTwo);
gradTwo[2].mul(lTwo);
basis[s][i] = new vector<Polynomial>(gradOne);
// Cell Based (Type 1) // basis[s][i]->at(0).sub(gradTwo[0]);
for(int l1 = 1; l1 < orderPlus; l1++){ basis[s][i]->at(1).sub(gradTwo[1]);
for(int l2 = 1; l2 < orderPlus; l2++){ basis[s][i]->at(2).sub(gradTwo[2]);
(*cell)[i] =
new vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
i++; i++;
} }
} }
// Cell Based (Type 2) // // Type 3.1
for(int l1 = 1; l1 < orderPlus; l1++){ for(unsigned int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){ Polynomial tmp =
(*cell)[i] = new vector<Polynomial>(3); intLegendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
lifting[(*(*faceV[s])[0])[3]]);
(*cell)[i]->at(0) = legendreX[l1] * iLegendreY[l2];
(*cell)[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1; basis[s][i] =
(*cell)[i]->at(2) = zero; new vector<Polynomial>((lifting[(*(*faceV[s])[0])[0]] -
lifting[(*(*faceV[s])[0])[1]]).gradient());
basis[s][i]->at(0).mul(tmp);
basis[s][i]->at(1).mul(tmp);
basis[s][i]->at(2).mul(tmp);
i++; i++;
} }
}
// Cell Based (Type 3) // // Type 3.2
for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){ for(unsigned int l1 = 1; l1 < orderPlus; l1++){
(*cell)[i] = new vector<Polynomial>(3); Polynomial tmp =
(*cell)[iPlus] = new vector<Polynomial>(3); intLegendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
lifting[(*(*faceV[s])[0])[1]]);
(*cell)[i]->at(0) = iLegendreY[l]; basis[s][i] =
(*cell)[i]->at(1) = zero; new vector<Polynomial>((lifting[(*(*faceV[s])[0])[0]] -
(*cell)[i]->at(2) = zero; lifting[(*(*faceV[s])[0])[3]]).gradient());
(*cell)[iPlus]->at(0) = zero; basis[s][i]->at(0).mul(tmp);
(*cell)[iPlus]->at(1) = iLegendreX[l]; basis[s][i]->at(1).mul(tmp);
(*cell)[iPlus]->at(2) = zero; basis[s][i]->at(2).mul(tmp);
i++; i++;
} }
}
// Mapping to Gmsh Quad // // Mapping to Gmsh Quad //
// x = (u + 1) / 2 // x = (u + 1) / 2
...@@ -214,56 +220,37 @@ QuadEdgeBasis::QuadEdgeBasis(int order){ ...@@ -214,56 +220,37 @@ QuadEdgeBasis::QuadEdgeBasis(int order){
Polynomial mapY(Polynomial(0.5, 0, 1, 0) + Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
Polynomial(0.5, 0, 0, 0)); Polynomial(0.5, 0, 0, 0));
for(int i = 0; i < nEdgeClosure; i++){ for(unsigned int s = 0; s < nRefSpace; s++){
for(int j = 0; j < nEdge; j++){ for(unsigned int i = 0; i < nFunction; i++){
(*(*edge)[i])[j]->at(0) = (*(*edge)[i])[j]->at(0).compose(mapX, mapY); vector<Polynomial>* old;
(*(*edge)[i])[j]->at(1) = (*(*edge)[i])[j]->at(1).compose(mapX, mapY); vector<Polynomial> nxt(3);
(*(*edge)[i])[j]->at(2) = (*(*edge)[i])[j]->at(2).compose(mapX, mapY);
} old = basis[s][i];
} nxt[0] = (*old)[0].compose(mapX, mapY);
nxt[1] = (*old)[1].compose(mapX, mapY);
nxt[2] = (*old)[2].compose(mapX, mapY);
for(int i = 0; i < nCell; i++){ basis[s][i] = new vector<Polynomial>(nxt);
(*cell)[i]->at(0) = (*cell)[i]->at(0).compose(mapX, mapY); delete old;
(*cell)[i]->at(1) = (*cell)[i]->at(1).compose(mapX, mapY); }
(*cell)[i]->at(2) = (*cell)[i]->at(2).compose(mapX, mapY);
} }
// Free Temporary Sapce // // Free Temporary Sapce //
delete[] legendre; delete[] legendre;
delete[] intLegendre; delete[] intLegendre;
delete[] iLegendreX;
delete[] iLegendreY;
delete[] legendreX;
delete[] legendreY;
} }
QuadEdgeBasis::~QuadEdgeBasis(void){ QuadEdgeBasis::~QuadEdgeBasis(void){
// Vertex Based // // ReferenceSpace //
for(int i = 0; i < nVertex; i++) delete refSpace;
delete (*node)[i];
delete node; // Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
// Edge Based // delete basis[i][j];
for(int c = 0; c < 2; c++){
for(int i = 0; i < nEdge; i++)
delete (*(*edge)[c])[i];
delete (*edge)[c]; delete[] basis[i];
} }
delete edge; delete[] basis;
// Face Based //
delete face;
// Cell Based //
for(int i = 0; i < nCell; i++)
delete (*cell)[i];
delete cell;
} }
...@@ -24,7 +24,7 @@ class QuadEdgeBasis: public BasisHierarchical1From{ ...@@ -24,7 +24,7 @@ class QuadEdgeBasis: public BasisHierarchical1From{
//! @param order The order of the Basis //! @param order The order of the Basis
//! //!
//! Returns a new Edge-Basis for Quads of the given order //! Returns a new Edge-Basis for Quads of the given order
QuadEdgeBasis(int order); QuadEdgeBasis(unsigned int order);
//! Deletes this Basis //! Deletes this Basis
//! //!
......
#include "QuadNedelecBasis.h"
#include "QuadReferenceSpace.h"
#include "Legendre.h"
using namespace std;
QuadNedelecBasis::QuadNedelecBasis(void){
// Reference Space //
refSpace = new QuadReferenceSpace;
nRefSpace = refSpace->getNPermutation();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
// Set Basis Type //
order = 0;
type = 1;
dim = 2;
nVertex = 0;
nEdge = 4;
nFace = 0;
nCell = 0;
nFunction = 4;
// Lagrange & Lifting //
const Polynomial lagrange[4] =
{
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)))
};
const Polynomial lifting[4] =
{
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0))),
Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0)))
};
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int 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++){
Polynomial lambda = (lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]]) * 0.5;
basis[s][e] =
new vector<Polynomial>((lifting[(*(*edgeV[s])[e])[1]] -
lifting[(*(*edgeV[s])[e])[0]]).gradient());
basis[s][e]->at(0).mul(lambda);
basis[s][e]->at(1).mul(lambda);
basis[s][e]->at(2).mul(lambda);
}
}
// Mapping to Gmsh Quad //
// x = (u + 1) / 2
// y = (v + 1) / 2
//
// (x, y) = Zaglmayr Ref Quad
// (u, v) = Gmsh Ref Quad
Polynomial mapX(Polynomial(0.5, 1, 0, 0) +
Polynomial(0.5, 0, 0, 0));
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++){
vector<Polynomial>* old;
vector<Polynomial> nxt(3);
old = basis[s][i];
nxt[0] = (*old)[0].compose(mapX, mapY);
nxt[1] = (*old)[1].compose(mapX, mapY);
nxt[2] = (*old)[2].compose(mapX, mapY);
basis[s][i] = new vector<Polynomial>(nxt);
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++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
}
#ifndef _QUADNEDELECBASIS_H_
#define _QUADNEDELECBASIS_H_
#include "BasisHierarchical1From.h"
/**
@class QuadNedelecBasis
@brief Nedelec Basis for Quads
This class can instantiate a Nedelec Basis
for Quadrangles.@n
*/
class QuadNedelecBasis: public BasisHierarchical1From{
public:
//! Returns a new Nedelec Basis for Quadrangles
//!
QuadNedelecBasis(void);
//! Deletes this Basis
//!
virtual ~QuadNedelecBasis(void);
};
#endif
...@@ -118,14 +118,6 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){ ...@@ -118,14 +118,6 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
} }
} }
cout << (*(*faceV[1])[0])[0] << ", "
<< (*(*faceV[1])[0])[1] << endl
<< (*(*faceV[1])[0])[0] << ", "
<< (*(*faceV[1])[0])[3] << endl
<< (*(*faceV[1])[0])[2] << endl;
// Mapping to Gmsh Quad // // Mapping to Gmsh Quad //
// x = (u + 1) / 2 // x = (u + 1) / 2
// y = (v + 1) / 2 // y = (v + 1) / 2
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment