diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 82f64d50411061475ed51c89906df8e74639d4a9..06c52de106d1b728dca5a8e243af245ad18ddbc8 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -14,6 +14,7 @@
#include "QuadNodeBasis.h"
#include "QuadEdgeBasis.h"
+#include "QuadLagrangeBasis.h"
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
@@ -71,7 +72,7 @@ BasisLocal* BasisGenerator::generateLagrange(unsigned int elementType,
switch(elementType){
case TYPE_LIN: return new LineLagrangeBasis(order);
case TYPE_TRI: return new TriLagrangeBasis(order);
- case TYPE_QUA: throw Exception("Lagrange Basis on Quads not Implemented");
+ case TYPE_QUA: return new QuadLagrangeBasis(order);
case TYPE_TET: return new TetLagrangeBasis(order);
case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index d93e55d8510b80c6a33084990b8771598648e727..c73a7d795f8bf65c4f26c6683cf899550f24452b 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -33,6 +33,7 @@ set(SRC
QuadNodeBasis.cpp
# QuadEdgeBasis.cpp
+ QuadLagrangeBasis.cpp
TetNodeBasis.cpp
TetEdgeBasis.cpp
diff --git a/FunctionSpace/LineLagrangeBasis.cpp b/FunctionSpace/LineLagrangeBasis.cpp
index ec2eb76cfcfcc9a152ecd5be8bca6f92a0943b5b..26c0bb64bad857f2cccab6e73c382525fbc660f4 100644
--- a/FunctionSpace/LineLagrangeBasis.cpp
+++ b/FunctionSpace/LineLagrangeBasis.cpp
@@ -1,5 +1,6 @@
#include "Exception.h"
#include "LineLagrangeBasis.h"
+#include "pointsGenerators.h"
LineLagrangeBasis::LineLagrangeBasis(unsigned int order){
// If order 0 (Nedelec): use order 1
@@ -22,7 +23,8 @@ LineLagrangeBasis::LineLagrangeBasis(unsigned int order){
lBasis = new polynomialBasis(getTag(order));
// Init Lagrange Point //
- lPoint = new fullMatrix<double>(linePoint(order));
+ lPoint = new fullMatrix<double>
+ (gmshGeneratePointsLine(order));
}
LineLagrangeBasis::~LineLagrangeBasis(void){
@@ -41,23 +43,3 @@ unsigned int LineLagrangeBasis::getTag(unsigned int order){
("Can't instanciate an order %d Lagrangian Basis for a Line",
order);
}
-
-fullMatrix<double> LineLagrangeBasis::
-linePoint(unsigned int order){
- fullMatrix<double> line(order + 1, 3);
- line(0, 0) = 0;
- line(0, 1) = 0;
- line(0, 2) = 0;
-
- if(order > 0){
- line(0, 0) = -1.;
- line(1, 0) = 1.;
-
- double dd = 2. / order;
-
- for(unsigned int i = 2; i < order + 1; i++)
- line(i, 0) = -1. + dd * (i - 1);
- }
-
- return line;
-}
diff --git a/FunctionSpace/LineLagrangeBasis.h b/FunctionSpace/LineLagrangeBasis.h
index a1b7d77faaca5a1aa218ee4e436c990aad4b19e8..af48c318d8f729a36f226c598009437547601a91 100644
--- a/FunctionSpace/LineLagrangeBasis.h
+++ b/FunctionSpace/LineLagrangeBasis.h
@@ -31,11 +31,6 @@ class LineLagrangeBasis: public BasisLagrange{
//! @return Returns the @em tag of a @em Line of
//! the given order
static unsigned int getTag(unsigned int order);
-
- //! @param order A natural number
- //! @return Returns Lagrangian Points on a Line
- //! for the given Order
- static fullMatrix<double> linePoint(unsigned int order);
};
#endif
diff --git a/FunctionSpace/QuadLagrangeBasis.cpp b/FunctionSpace/QuadLagrangeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b8f78f2f7629f48d4bb6e5ff4db424560ad7c324
--- /dev/null
+++ b/FunctionSpace/QuadLagrangeBasis.cpp
@@ -0,0 +1,45 @@
+#include "Exception.h"
+#include "QuadLagrangeBasis.h"
+#include "pointsGenerators.h"
+
+QuadLagrangeBasis::QuadLagrangeBasis(unsigned int order){
+ // If order 0 (Nedelec): use order 1
+ if(order == 0)
+ order = 1;
+
+ // Set Basis Type //
+ this->order = order;
+
+ type = 0;
+ dim = 2;
+
+ nVertex = 4;
+ nEdge = 4 * (order - 1);
+ nFace = (order - 1) * (order - 1);
+ nCell = 0;
+ nFunction = nVertex + nEdge + nFace + nCell;
+
+ // Init polynomialBasis //
+ lBasis = new polynomialBasis(getTag(order));
+
+ // Init Lagrange Point //
+ lPoint = new fullMatrix<double>
+ (gmshGeneratePointsQuadrangle(order, false));
+}
+
+QuadLagrangeBasis::~QuadLagrangeBasis(void){
+ delete lBasis;
+ delete lPoint;
+}
+
+unsigned int QuadLagrangeBasis::getTag(unsigned int order){
+ unsigned int tag = nodalBasis::getTag(TYPE_QUA, order, false);
+
+ if(tag)
+ return tag;
+
+ else
+ throw Exception
+ ("Can't instanciate an order %d Lagrangian Basis for a Quadrangle",
+ order);
+}
diff --git a/FunctionSpace/QuadLagrangeBasis.h b/FunctionSpace/QuadLagrangeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1efbd50d3a0ebbb8296cbd45deaddd480ccd0cf
--- /dev/null
+++ b/FunctionSpace/QuadLagrangeBasis.h
@@ -0,0 +1,36 @@
+#ifndef _QUADLAGRANGEBASIS_H_
+#define _QUADLAGRANGEBASIS_H_
+
+#include "BasisLagrange.h"
+
+/**
+ @class QuadLagrangeBasis
+ @brief Lagrange Basis for Quadrangles
+
+ This class can instantiate a @em Lagrange @em Basis
+ for a Quadrangle and for a given Order.@n
+
+ It uses
+ <a href="http://geuz.org/gmsh/">gmsh</a> Basis.
+ */
+
+class QuadLagrangeBasis: public BasisLagrange{
+ public:
+ //! @param order A natural number
+ //!
+ //! Returns a new QuadLagrangeBasis
+ //! of the given Order
+ QuadLagrangeBasis(unsigned int order);
+
+ //! Deletes this Basis
+ //!
+ virtual ~QuadLagrangeBasis(void);
+
+ private:
+ //! @param order A natural number
+ //! @return Returns the @em tag of a @em Quadrangle of
+ //! the given order
+ static unsigned int getTag(unsigned int order);
+};
+
+#endif
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index a2e89c8ade7857fe0102b8c9be918f85d08b14ca..0a89ba922ba471a3731877f1d1eee98a6d2c6b27 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -290,11 +290,14 @@ inOrder(unsigned int permutation,
unsigned int c,
unsigned int d){
+ unsigned int opposite;
+
unsigned int v;
unsigned int k = 0;
vector<unsigned int>* inorder =
new vector<unsigned int>(4);
+ // Take nodes in order //
for(unsigned int i = 0; i < nVertex; i++){
v = perm[permutation][i];
@@ -304,6 +307,26 @@ inOrder(unsigned int permutation,
}
}
+ // We need node at inoder[2] to be //
+ // opposite to node at inorder[0] //
+ // --> if not make a permutation such //
+ // that node opposite at inorder[2] //
+ // is at inorder[0] //
+
+ opposite = ((*inorder)[0] + 2) % 4;
+ if((*inorder)[2] != opposite){
+ // Find opposite node index
+ unsigned int tmp;
+ unsigned int idx = 1;
+ while((*inorder)[idx] != opposite)
+ idx++;
+
+ // Swap inorder[2] and inorder[idx]
+ tmp = (*inorder)[2];
+ (*inorder)[2] = opposite;
+ (*inorder)[idx] = tmp;
+ }
+
return inorder;
}
@@ -324,24 +347,6 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
// (vertex[i].second->getNum) //
std::sort(vertex.begin(), vertex.end(), sortPredicate);
-
- /*****************************************************************
- * What does that do ? *
- * Tree lookup needs vertices in the range [0 .. nVertex[. *
- * *
- * So we need to convert the node ID from 'element.getVertex()' *
- * to a number between [0 .. nVertex[. *
- * *
- * The convertion is such that the smallest node ID *
- * gets the converted ID '0'. *
- * Then the second smallest the ID '1'. *
- * And so on... *
- * *
- * The sorting of the vector 'vertex' with respect *
- * to the element.getVertex() IDs does that job *
- *****************************************************************/
-
-
// Tree Lookup //
try{
return treeLookup(&pTreeRoot, vertex);
diff --git a/FunctionSpace/TetLagrangeBasis.cpp b/FunctionSpace/TetLagrangeBasis.cpp
index 2c2c382b792fe363b2de559e5495c68d582fe834..6d4fcbc321b7081d8596a5584166813f9e39cd20 100644
--- a/FunctionSpace/TetLagrangeBasis.cpp
+++ b/FunctionSpace/TetLagrangeBasis.cpp
@@ -1,5 +1,6 @@
#include "Exception.h"
#include "TetLagrangeBasis.h"
+#include "pointsGenerators.h"
TetLagrangeBasis::TetLagrangeBasis(unsigned int order){
// If order 0 (Nedelec): use order 1
@@ -22,7 +23,8 @@ TetLagrangeBasis::TetLagrangeBasis(unsigned int order){
lBasis = new polynomialBasis(getTag(order));
// Init Lagrange Point //
- lPoint = new fullMatrix<double>(tetPoint(order));
+ lPoint = new fullMatrix<double>
+ (gmshGeneratePointsTetrahedron(order, false));
}
TetLagrangeBasis::~TetLagrangeBasis(void){
@@ -41,336 +43,3 @@ unsigned int TetLagrangeBasis::getTag(unsigned int order){
("Can't instanciate an order %d Lagrangian Basis for a Tetrahedron",
order);
}
-
-fullMatrix<double> TetLagrangeBasis::
-tetPoint(unsigned int order){
- unsigned int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
- fullMatrix<double> point(nbPoints, 3);
-
- double overOrder = (order == 0 ? 1. : 1. / order);
-
- point(0, 0) = 0.;
- point(0, 1) = 0.;
- point(0, 2) = 0.;
-
- if(order > 0){
- point(1, 0) = order;
- point(1, 1) = 0;
- point(1, 2) = 0;
-
- point(2, 0) = 0.;
- point(2, 1) = order;
- point(2, 2) = 0.;
-
- point(3, 0) = 0.;
- point(3, 1) = 0.;
- point(3, 2) = order;
-
- if(order > 1){
- for(unsigned int k = 0; k < (order - 1); k++){
- point(4 + k, 0) = k + 1;
- point(4 + order - 1 + k, 0) = order - 1 - k;
- point(4 + 2 * (order - 1) + k, 0) = 0.;
- point(4 + 3 * (order - 1) + k, 0) = 0.;
- point(4 + 4 * (order - 1) + k, 0) = 0.;
- point(4 + 5 * (order - 1) + k, 0) = k+1;
-
- point(4 + k, 1) = 0.;
- point(4 + order - 1 + k, 1) = k + 1;
- point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
- point(4 + 3 * (order - 1) + k, 1) = 0.;
- point(4 + 4 * (order - 1) + k, 1) = k+1;
- point(4 + 5 * (order - 1) + k, 1) = 0.;
-
- point(4 + k, 2) = 0.;
- point(4 + order - 1 + k, 2) = 0.;
- point(4 + 2 * (order - 1) + k, 2) = 0.;
- point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
- }
-
- if(order > 2){
- unsigned int ns = 4 + 6 * (order - 1);
- unsigned int nbdofface = nbdoftriangle(order - 3);
-
- double *u = new double[nbdofface];
- double *v = new double[nbdofface];
- double *w = new double[nbdofface];
-
- nodepositionface0(order - 3, u, v, w);
-
- // u-v plane
-
- for(unsigned int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3);
- }
-
- ns = ns + nbdofface;
-
- // u-w plane
-
- nodepositionface1(order - 3, u, v, w);
-
- for(unsigned int i=0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) ;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- // v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface2(order - 3, u, v, w);
-
- for(unsigned int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3);
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- // u-v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface3(order - 3, u, v, w);
-
- for(unsigned int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if(order > 3){
- fullMatrix<double> interior = tetPoint(order - 4);
-
- for(int k = 0; k < interior.size1(); k++){
- point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
- point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
- point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
- }
- }
- }
- }
- }
-
- point.scale(overOrder);
- return point;
-}
-
-unsigned int TetLagrangeBasis::nbdoftriangle(unsigned int order){
- return (order + 1) * (order + 2) / 2;
-}
-
-void TetLagrangeBasis::
-nodepositionface0(unsigned int order, double *u, double *v, double *w){
-
- unsigned int ndofT = nbdoftriangle(order);
-
- if(order == 0){
- u[0] = 0.;
- v[0] = 0.;
- w[0] = 0.;
- return;
- }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= order; v[2]= 0.; w[2] = 0.;
-
- // edges
- for(unsigned int k = 0; k < (order - 1); k++){
- u[3 + k] = 0.;
- v[3 + k] = k + 1;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = k + 1;
- v[3 + order - 1 + k] = order - 1 - k ;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = order - 1 - k;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = 0.;
- }
-
- if(order > 2){
- unsigned int nbdoftemp = nbdoftriangle(order - 3);
-
- nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3* (order - 1)]);
-
- for(unsigned int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
- }
- }
-
- for(unsigned int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-void TetLagrangeBasis::nodepositionface1
-(unsigned int order, double *u, double *v, double *w){
-
- unsigned int ndofT = nbdoftriangle(order);
-
- if(order == 0){
- u[0] = 0.;
- v[0] = 0.;
- w[0] = 0.;
- return;
- }
-
- u[0] = 0.; v[0]= 0.; w[0] = 0.;
- u[1] = order; v[1]= 0.; w[1] = 0.;
- u[2] = 0.; v[2]= 0.; w[2] = order;
-
- // edges
- for(unsigned int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = 0.;
- w[3 + order - 1+ k ] = k + 1;
-
- u[3 + 2 * (order - 1) + k] = 0. ;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
- }
-
- if(order > 2){
- unsigned int nbdoftemp = nbdoftriangle(order - 3);
-
- nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
- &w[3 + 3 * (order - 1)]);
-
- for(unsigned int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
-
- for (unsigned int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-void TetLagrangeBasis::
-nodepositionface2(unsigned int order, double *u, double *v, double *w){
-
- unsigned int ndofT = nbdoftriangle(order);
-
- if(order == 0){
- u[0] = 0.;
- v[0] = 0.;
- w[0] = 0.;
- return;
- }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= 0.; w[1] = order;
- u[2]= 0.; v[2]= order; w[2] = 0.;
-
- // edges
- for(unsigned int k = 0; k < (order - 1); k++){
- u[3 + k] = 0.;
- v[3 + k] = 0.;
- w[3 + k] = k + 1;
-
- u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = order - 1 - k;
-
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
- w[3 + 2 * (order - 1) + k] = 0.;
- }
-
- if(order > 2){
- unsigned int nbdoftemp = nbdoftriangle(order - 3);
-
- nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3 * (order - 1)]);
-
- for(unsigned int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
-
- for(unsigned int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-void TetLagrangeBasis::
-nodepositionface3(unsigned int order, double *u, double *v, double *w){
-
- unsigned int ndofT = nbdoftriangle(order);
-
- if(order == 0){
- u[0] = 0.;
- v[0] = 0.;
- w[0] = 0.;
- return;
- }
-
- u[0]= 0.; v[0]= 0.; w[0] = order;
- u[1]= order; v[1]= 0.; w[1] = 0.;
- u[2]= 0.; v[2]= order; w[2] = 0.;
-
- // edges
- for(unsigned int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = order - 1 - k;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
- w[3 + 2 * (order - 1) + k] = k + 1;
- }
-
- if(order > 2){
- unsigned int nbdoftemp = nbdoftriangle(order - 3);
-
- nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3 * (order - 1)]);
-
- for(unsigned int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
-
- for (unsigned int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
diff --git a/FunctionSpace/TetLagrangeBasis.h b/FunctionSpace/TetLagrangeBasis.h
index 27f7e5a6ecf181025c677e6d020f6d3a9bad7b5a..ef7d8db129286fbf533e4bc173d586dcdd3382bc 100644
--- a/FunctionSpace/TetLagrangeBasis.h
+++ b/FunctionSpace/TetLagrangeBasis.h
@@ -31,30 +31,6 @@ class TetLagrangeBasis: public BasisLagrange{
//! @return Returns the @em tag of a @em Tetrahedron of
//! the given order
static unsigned int getTag(unsigned int order);
-
- //! @param order A natural number
- //! @return Returns Lagrangian Points on a Tetrahedron
- //! for the given Order
- static fullMatrix<double> tetPoint(unsigned int order);
-
- //! Unknown gmsh function
- static unsigned int nbdoftriangle(unsigned int order);
-
- //! Unknown gmsh function
- static void nodepositionface0(unsigned int order,
- double *u, double *v, double *w);
-
- //! Unknown gmsh function
- static void nodepositionface1(unsigned int order,
- double *u, double *v, double *w);
-
- //! Unknown gmsh function
- static void nodepositionface2(unsigned int order,
- double *u, double *v, double *w);
-
- //! Unknown gmsh function
- static void nodepositionface3(unsigned int order,
- double *u, double *v, double *w);
};
#endif
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 0812ff9b9ddff60283426ed2147043e480b57bf6..3da698d557160f36b788150c4c4355dda44378ed 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,5 +1,6 @@
#include "Exception.h"
#include "TriLagrangeBasis.h"
+#include "pointsGenerators.h"
TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
// If order 0 (Nedelec): use order 1
@@ -22,7 +23,8 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
lBasis = new polynomialBasis(getTag(order));
// Init Lagrange Point //
- lPoint = new fullMatrix<double>(triPoint(order));
+ lPoint = new fullMatrix<double>
+ (gmshGeneratePointsTriangle(order, false));
}
TriLagrangeBasis::~TriLagrangeBasis(void){
@@ -41,72 +43,3 @@ unsigned int TriLagrangeBasis::getTag(unsigned int order){
("Can't instanciate an order %d Lagrangian Basis for a Triangle",
order);
}
-
-fullMatrix<double> TriLagrangeBasis::
-triPoint(unsigned int order){
- unsigned int nbPoints = (order + 1) * (order + 2) / 2;
- fullMatrix<double> point(nbPoints, 3);
-
- point(0, 0) = 0;
- point(0, 1) = 0;
- point(0, 2) = 0;
-
- if(order > 0){
- double dd = 1. / order;
-
- point(1, 0) = 1;
- point(1, 1) = 0;
- point(1, 2) = 0;
-
- point(2, 0) = 0;
- point(2, 1) = 1;
- point(2, 2) = 0;
-
- unsigned int index = 3;
-
- if(order > 1){
- double ksi = 0;
- double eta = 0;
-
- for(unsigned int i = 0; i < order - 1; i++, index++){
- ksi += dd;
-
- point(index, 0) = ksi;
- point(index, 1) = eta;
- point(index, 2) = 0;
- }
-
- ksi = 1.;
-
- for(unsigned int i = 0; i < order - 1; i++, index++){
- ksi -= dd;
- eta += dd;
-
- point(index, 0) = ksi;
- point(index, 1) = eta;
- point(index, 2) = 0;
- }
-
- eta = 1.;
- ksi = 0.;
-
- for(unsigned int i = 0; i < order - 1; i++, index++){
- eta -= dd;
-
- point(index, 0) = ksi;
- point(index, 1) = eta;
- point(index, 2) = 0;
- }
-
- if(order > 2){
- fullMatrix<double> inner = triPoint(order - 3);
-
- inner.scale(1. - 3. * dd);
- inner.add(dd);
- point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
-
- return point;
-}
diff --git a/FunctionSpace/TriLagrangeBasis.h b/FunctionSpace/TriLagrangeBasis.h
index 272a21f32d84e27dea05dbe852389e507260aa40..a49ab78c9f842ef98e4323a25928baff68e2f440 100644
--- a/FunctionSpace/TriLagrangeBasis.h
+++ b/FunctionSpace/TriLagrangeBasis.h
@@ -31,11 +31,6 @@ class TriLagrangeBasis: public BasisLagrange{
//! @return Returns the @em tag of a @em Triangle of
//! the given order
static unsigned int getTag(unsigned int order);
-
- //! @param order A natural number
- //! @return Returns Lagrangian Points on a Triangle
- //! for the given Order
- static fullMatrix<double> triPoint(unsigned int order);
};
#endif