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

Closure for Triangles

parent bd4b62ce
Branches
Tags
No related merge requests found
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
class BasisVector: public Basis{ class BasisVector: public Basis{
protected: protected:
std::vector<const std::vector<Polynomial>*>* basis; std::vector<const std::vector<Polynomial>*>* basis;
std::vector<const std::vector<Polynomial>*>* revBasis;
public: public:
//! Deletes this BasisVector //! Deletes this BasisVector
...@@ -31,6 +32,12 @@ class BasisVector: public Basis{ ...@@ -31,6 +32,12 @@ class BasisVector: public Basis{
//! defining this (vectorial) Basis //! defining this (vectorial) Basis
const std::vector<const std::vector<Polynomial>*>& getFunctions(void) const; const std::vector<const std::vector<Polynomial>*>& getFunctions(void) const;
//! @param closure A natural number
//! @return Returns the set of @em Polynomial%s
//! defining this (scalar) Basis, for the given closure
const std::vector<const std::vector<Polynomial>*>&
getFunctions(unsigned int closure) const;
protected: protected:
//! @internal //! @internal
//! Instantiate a new BasisVector //! Instantiate a new BasisVector
...@@ -50,4 +57,15 @@ getFunctions(void) const{ ...@@ -50,4 +57,15 @@ getFunctions(void) const{
return *basis; return *basis;
} }
inline
const std::vector<const std::vector<Polynomial>*>& BasisVector::
getFunctions(unsigned int closure) const{
if(!closure)
return *basis;
else
return *revBasis;
}
#endif #endif
...@@ -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<const vector<Polynomial>*>& fun = getBasis(element).getFunctions(); const vector<const vector<Polynomial>*>& fun = getLocalFunctions(element);
const unsigned int nFun = fun.size(); const unsigned int nFun = fun.size();
// Get Reference coordinate // // Get Reference coordinate //
......
#include "FunctionSpaceVector.h" #include "FunctionSpaceVector.h"
#include "Exception.h"
using namespace std;
FunctionSpaceVector::~FunctionSpaceVector(void){ FunctionSpaceVector::~FunctionSpaceVector(void){
}
const vector<const vector<Polynomial>*> FunctionSpaceVector::
getLocalFunctions(const MElement& element) const{
// Get Closure //
map<const MElement*, vector<bool>*>::iterator it =
edgeClosure->find(&element);
if(it == edgeClosure->end())
throw Exception("Element not found for closure");
vector<bool>* closure = it->second;
const unsigned int size = closure->size();
// Get Basis //
const vector<const vector<Polynomial>*>& basis =
getBasis(element).getFunctions(0);
const vector<const vector<Polynomial>*>& revBasis =
getBasis(element).getFunctions(1);
// Get Functions //
vector<const vector<Polynomial>*> fun(size);
for(unsigned int i = 0; i < size; i++)
if((*closure)[i])
fun[i] = basis[i];
else
fun[i] = revBasis[i];
// Return
return fun;
} }
...@@ -29,6 +29,9 @@ class FunctionSpaceVector : public FunctionSpace{ ...@@ -29,6 +29,9 @@ class FunctionSpaceVector : public FunctionSpace{
const std::vector<double>& coef, const std::vector<double>& coef,
const fullVector<double>& xyz) const = 0; const fullVector<double>& xyz) const = 0;
const std::vector<const std::vector<Polynomial>*>
getLocalFunctions(const MElement& element) const;
const BasisVector& getBasis(const MElement& element) const; const BasisVector& getBasis(const MElement& element) const;
}; };
......
...@@ -21,14 +21,14 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -21,14 +21,14 @@ TriEdgeBasis::TriEdgeBasis(const int order){
Polynomial* legendre = new Polynomial[orderPlus]; Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus]; Polynomial* intLegendre = new Polynomial[orderPlus];
Polynomial* lagrange = new Polynomial[3];
Polynomial* lagrangeSub = new Polynomial[3];
Polynomial* lagrangeSum = new Polynomial[3];
Polynomial* u = new Polynomial[orderPlus]; Polynomial* u = new Polynomial[orderPlus];
Polynomial* v = new Polynomial[orderPlus]; Polynomial* v = new Polynomial[orderPlus];
Polynomial lagrange [3];
Polynomial lagrangeSum [3];
Polynomial lagrangeSub [3];
Polynomial rLagrangeSub [3];
// Classical and Intrated-Scaled Legendre Polynomial // // Classical and Intrated-Scaled Legendre Polynomial //
Legendre::legendre(legendre, order); Legendre::legendre(legendre, order);
Legendre::intScaled(intLegendre, orderPlus); Legendre::intScaled(intLegendre, orderPlus);
...@@ -47,19 +47,23 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -47,19 +47,23 @@ TriEdgeBasis::TriEdgeBasis(const int order){
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] = lagrange[j] + lagrange[i]; lagrangeSum[i] = lagrange[j] + lagrange[i];
// Lagrange Sub // // Lagrange Sub (& revert) //
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] = lagrange[i] - lagrange[j]; lagrangeSub[i] = lagrange[i] - lagrange[j];
rLagrangeSub[i] = lagrange[j] - lagrange[i];
}
// Basis (temporary --- *no* const) // // Basis & Revert (temporary --- *no* const) //
std::vector<std::vector<Polynomial>*> basis(size); std::vector<std::vector<Polynomial>*> basis(size);
std::vector<std::vector<Polynomial>*> revBasis(size);
// Edge Based (Nedelec) // // Edge Based (Nedelec) //
int i = 0; int i = 0;
for(int j = 1; i < 3; j = (j + 1) % 3){ for(int j = 1; i < 3; j = (j + 1) % 3){
// Direct
std::vector<Polynomial> tmp = lagrange[j].gradient(); std::vector<Polynomial> tmp = lagrange[j].gradient();
tmp[0].mul(lagrange[i]); tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]); tmp[1].mul(lagrange[i]);
...@@ -75,6 +79,23 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -75,6 +79,23 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(1).sub(tmp[1]); basis[i]->at(1).sub(tmp[1]);
basis[i]->at(2).sub(tmp[2]); basis[i]->at(2).sub(tmp[2]);
// Revert
std::vector<Polynomial> tmpR = lagrange[i].gradient();
tmpR[0].mul(lagrange[j]);
tmpR[1].mul(lagrange[j]);
tmpR[2].mul(lagrange[j]);
revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient());
revBasis[i]->at(0).mul(lagrange[i]);
revBasis[i]->at(1).mul(lagrange[i]);
revBasis[i]->at(2).mul(lagrange[i]);
revBasis[i]->at(0).sub(tmpR[0]);
revBasis[i]->at(1).sub(tmpR[1]);
revBasis[i]->at(2).sub(tmpR[2]);
// Next
i++; i++;
} }
...@@ -82,8 +103,12 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -82,8 +103,12 @@ TriEdgeBasis::TriEdgeBasis(const int 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] =
new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient()); new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e],
lagrangeSum[e])).gradient());
revBasis[i] =
new std::vector<Polynomial>((intLegendre[l].compose(rLagrangeSub[e],
lagrangeSum[e])).gradient());
i++; i++;
} }
} }
...@@ -115,6 +140,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -115,6 +140,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(1).add(tmp[1]); basis[i]->at(1).add(tmp[1]);
basis[i]->at(2).add(tmp[2]); basis[i]->at(2).add(tmp[2]);
revBasis[i] = basis[i];
i++; i++;
} }
} }
...@@ -137,6 +164,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -137,6 +164,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(1).sub(tmp[1]); basis[i]->at(1).sub(tmp[1]);
basis[i]->at(2).sub(tmp[2]); basis[i]->at(2).sub(tmp[2]);
revBasis[i] = basis[i];
i++; i++;
} }
} }
...@@ -149,15 +178,14 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -149,15 +178,14 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(1).mul(v[l]); basis[i]->at(1).mul(v[l]);
basis[i]->at(2).mul(v[l]); basis[i]->at(2).mul(v[l]);
revBasis[i] = basis[i];
i++; i++;
} }
// Free Temporary Sapce // // Free Temporary Sapce //
delete[] legendre; delete[] legendre;
delete[] intLegendre; delete[] intLegendre;
delete[] lagrange;
delete[] lagrangeSub;
delete[] lagrangeSum;
delete[] u; delete[] u;
delete[] v; delete[] v;
...@@ -165,11 +193,19 @@ TriEdgeBasis::TriEdgeBasis(const int order){ ...@@ -165,11 +193,19 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Set Basis // // Set Basis //
this->basis = new std::vector<const std::vector<Polynomial>*> this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end()); (basis.begin(), basis.end());
this->revBasis = new std::vector<const std::vector<Polynomial>*>
(revBasis.begin(), revBasis.end());
} }
TriEdgeBasis::~TriEdgeBasis(void){ TriEdgeBasis::~TriEdgeBasis(void){
for(int i = 0; i < size; i++) for(int i = 0; i < size; i++){
delete (*basis)[i]; delete (*basis)[i];
if(i >= nVertex && i < nVertex + nEdge)
delete (*revBasis)[i];
}
delete basis; delete basis;
delete revBasis;
} }
...@@ -15,7 +15,7 @@ TriNedelecBasis::TriNedelecBasis(void){ ...@@ -15,7 +15,7 @@ TriNedelecBasis::TriNedelecBasis(void){
size = 3; size = 3;
// Lagrange // // Lagrange //
Polynomial* lagrange = new Polynomial[3]; Polynomial lagrange[3];
lagrange[0] = lagrange[0] =
Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0); Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
...@@ -29,9 +29,11 @@ TriNedelecBasis::TriNedelecBasis(void){ ...@@ -29,9 +29,11 @@ TriNedelecBasis::TriNedelecBasis(void){
// Basis (temporary --- *no* const) // // Basis (temporary --- *no* const) //
std::vector<std::vector<Polynomial>*> basis(size); std::vector<std::vector<Polynomial>*> basis(size);
std::vector<std::vector<Polynomial>*> revBasis(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){
// Direct
std::vector<Polynomial> tmp = lagrange[j].gradient(); std::vector<Polynomial> tmp = lagrange[j].gradient();
tmp[0].mul(lagrange[i]); tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]); tmp[1].mul(lagrange[i]);
...@@ -46,15 +48,31 @@ TriNedelecBasis::TriNedelecBasis(void){ ...@@ -46,15 +48,31 @@ TriNedelecBasis::TriNedelecBasis(void){
basis[i]->at(0).sub(tmp[0]); basis[i]->at(0).sub(tmp[0]);
basis[i]->at(1).sub(tmp[1]); basis[i]->at(1).sub(tmp[1]);
basis[i]->at(2).sub(tmp[2]); basis[i]->at(2).sub(tmp[2]);
}
// Free Temporary Sapce // // Revert
delete[] lagrange; std::vector<Polynomial> tmpR = lagrange[i].gradient();
tmpR[0].mul(lagrange[j]);
tmpR[1].mul(lagrange[j]);
tmpR[2].mul(lagrange[j]);
revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient());
revBasis[i]->at(0).mul(lagrange[i]);
revBasis[i]->at(1).mul(lagrange[i]);
revBasis[i]->at(2).mul(lagrange[i]);
revBasis[i]->at(0).sub(tmpR[0]);
revBasis[i]->at(1).sub(tmpR[1]);
revBasis[i]->at(2).sub(tmpR[2]);
}
// Set Basis // // Set Basis //
this->basis = new std::vector<const std::vector<Polynomial>*> this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end()); (basis.begin(), basis.end());
this->revBasis = new std::vector<const std::vector<Polynomial>*>
(revBasis.begin(), revBasis.end());
} }
TriNedelecBasis::~TriNedelecBasis(void){ TriNedelecBasis::~TriNedelecBasis(void){
...@@ -62,4 +80,9 @@ TriNedelecBasis::~TriNedelecBasis(void){ ...@@ -62,4 +80,9 @@ TriNedelecBasis::~TriNedelecBasis(void){
delete (*basis)[i]; delete (*basis)[i];
delete basis; delete basis;
for(int i = 0; i < size; i++)
delete (*revBasis)[i];
delete revBasis;
} }
...@@ -19,10 +19,9 @@ TriNodeBasis::TriNodeBasis(const int order){ ...@@ -19,10 +19,9 @@ TriNodeBasis::TriNodeBasis(const int order){
Polynomial* legendre = new Polynomial[order]; Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order]; Polynomial* intLegendre = new Polynomial[order];
Polynomial lagrangeSub[3];
Polynomial lagrangeSum[3]; Polynomial lagrangeSum[3];
Polynomial lagrangeSub[3];
Polynomial rLagrangeSub[3]; Polynomial rLagrangeSub[3];
Polynomial rLagrangeSum[3];
// Classical and Intrated-Scaled Legendre Polynomial // // Classical and Intrated-Scaled Legendre Polynomial //
const int orderMinus = order - 1; const int orderMinus = order - 1;
...@@ -51,13 +50,9 @@ TriNodeBasis::TriNodeBasis(const int order){ ...@@ -51,13 +50,9 @@ TriNodeBasis::TriNodeBasis(const int order){
for(int i = 0; i < 3; i++) for(int i = 0; i < 3; i++)
(*revBasis)[i] = (*basis)[i]; (*revBasis)[i] = (*basis)[i];
// Lagrange Sum (& revert) // // Lagrange Sum (& revert) //
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];
rLagrangeSum[i] = *(*basis)[j] + *(*basis)[i];
}
// Lagrange Sub (& revert) // // Lagrange Sub (& revert) //
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){
...@@ -75,7 +70,7 @@ TriNodeBasis::TriNodeBasis(const int order){ ...@@ -75,7 +70,7 @@ TriNodeBasis::TriNodeBasis(const int order){
intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])); intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]));
(*revBasis)[i] = new Polynomial( (*revBasis)[i] = new Polynomial(
intLegendre[l].compose(rLagrangeSub[e], rLagrangeSum[e])); intLegendre[l].compose(rLagrangeSub[e], lagrangeSum[e]));
i++; i++;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment