From 9033bad71a61e31997ef14432d47392ea271c740 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Thu, 22 Nov 2012 17:25:29 +0000
Subject: [PATCH] Now Lagrange Basis (for Tri) completly in Basis Framework --
Still need a sexy way to choose between Zaglmayr and Lagrange -- Need to put
LagrangeBasis::project elsewhere
---
FunctionSpace/BasisGenerator.cpp | 69 ++++++++++++++++++------
FunctionSpace/BasisGenerator.h | 82 +++++++++++++++++++++++++----
FunctionSpace/CMakeLists.txt | 5 +-
FunctionSpace/FunctionSpace.cpp | 2 +-
FunctionSpace/LagrangeBasis.cpp | 39 ++++++++++++++
FunctionSpace/LagrangeBasis.h | 10 ++++
FunctionSpace/LagrangeGenerator.cpp | 22 --------
FunctionSpace/LagrangeGenerator.h | 52 ------------------
FunctionSpace/TriLagrangeBasis.cpp | 81 ++++++++++++++++++++++++----
9 files changed, 246 insertions(+), 116 deletions(-)
delete mode 100644 FunctionSpace/LagrangeGenerator.cpp
delete mode 100644 FunctionSpace/LagrangeGenerator.h
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 299201fe79..cfb3efbeb4 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -12,6 +12,7 @@
#include "TriNodeBasis.h"
#include "TriEdgeBasis.h"
#include "TriNedelecBasis.h"
+#include "TriLagrangeBasis.h"
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
@@ -28,21 +29,57 @@ BasisGenerator::~BasisGenerator(void){
Basis* BasisGenerator::generate(int elementType,
int basisType,
- int order){
+ int order,
+ std::string family){
+
+ if(!family.compare(std::string("zaglmayr")))
+ return generateZaglmayr(elementType, basisType, order);
+
+ else if(!family.compare(std::string("lagrange")))
+ return generateLagrange(elementType, basisType, order);
+
+ else
+ throw Exception("Unknwown Basis Family: %s", family.c_str());
+}
+
+Basis* BasisGenerator::generateZaglmayr(int elementType,
+ int basisType,
+ int order){
+
+ switch(elementType){
+ case TYPE_LIN: return linZaglmayrGen(basisType, order);
+ case TYPE_TRI: return triZaglmayrGen(basisType, order);
+ case TYPE_QUA: return quaZaglmayrGen(basisType, order);
+ case TYPE_TET: return tetZaglmayrGen(basisType, order);
+ case TYPE_HEX: return hexZaglmayrGen(basisType, order);
+
+ default: throw Exception("Unknown Element Type (%d) for Basis Generation",
+ elementType);
+ }
+}
+
+Basis* BasisGenerator::generateLagrange(int elementType,
+ int basisType,
+ int order){
+ if(basisType != 0)
+ throw
+ Exception("Cannot Have a %d-Form Lagrange Basis (0-Form only)",
+ basisType);
+
switch(elementType){
- case TYPE_LIN: return linGen(basisType, order);
- case TYPE_TRI: return triGen(basisType, order);
- case TYPE_QUA: return quaGen(basisType, order);
- case TYPE_TET: return tetGen(basisType, order);
- case TYPE_HEX: return hexGen(basisType, order);
+ case TYPE_LIN: throw Exception("Lagrange Basis on Lines not Implemented");
+ case TYPE_TRI: return new TriLagrangeBasis(order);
+ case TYPE_QUA: throw Exception("Lagrange Basis on Quads not Implemented");
+ case TYPE_TET: throw Exception("Lagrange Basis on Tets not Implemented");
+ case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
elementType);
}
}
-Basis* BasisGenerator::linGen(int basisType,
- int order){
+Basis* BasisGenerator::linZaglmayrGen(int basisType,
+ int order){
switch(basisType){
case 0: return new LineNodeBasis(order);
case 1:
@@ -56,8 +93,8 @@ Basis* BasisGenerator::linGen(int basisType,
}
}
-Basis* BasisGenerator::triGen(int basisType,
- int order){
+Basis* BasisGenerator::triZaglmayrGen(int basisType,
+ int order){
switch(basisType){
case 0: return new TriNodeBasis(order);
case 1:
@@ -71,8 +108,8 @@ Basis* BasisGenerator::triGen(int basisType,
}
}
-Basis* BasisGenerator::quaGen(int basisType,
- int order){
+Basis* BasisGenerator::quaZaglmayrGen(int basisType,
+ int order){
switch(basisType){
case 0: return new QuadNodeBasis(order);
case 1: return new QuadEdgeBasis(order);
@@ -83,8 +120,8 @@ Basis* BasisGenerator::quaGen(int basisType,
}
}
-Basis* BasisGenerator::tetGen(int basisType,
- int order){
+Basis* BasisGenerator::tetZaglmayrGen(int basisType,
+ int order){
switch(basisType){
case 0: return new TetNodeBasis(order);
case 1: return new TetEdgeBasis(order);
@@ -95,8 +132,8 @@ Basis* BasisGenerator::tetGen(int basisType,
}
}
-Basis* BasisGenerator::hexGen(int basisType,
- int order){
+Basis* BasisGenerator::hexZaglmayrGen(int basisType,
+ int order){
switch(basisType){
case 0: return new HexNodeBasis(order);
case 1: return new HexEdgeBasis(order);
diff --git a/FunctionSpace/BasisGenerator.h b/FunctionSpace/BasisGenerator.h
index 1f071ea3bc..6d3351c319 100644
--- a/FunctionSpace/BasisGenerator.h
+++ b/FunctionSpace/BasisGenerator.h
@@ -1,6 +1,7 @@
#ifndef _BASISGENERATOR_H_
#define _BASISGENERATOR_H_
+#include <string>
#include "Basis.h"
/**
@@ -20,15 +21,29 @@ class BasisGenerator{
BasisGenerator(void);
~BasisGenerator(void);
+ static Basis* generate(int elementType,
+ int basisType,
+ int order,
+ std::string family);
+
static Basis* generate(int elementType,
int basisType,
int order);
- static Basis* linGen(int basisType, int order);
- static Basis* triGen(int basisType, int order);
- static Basis* quaGen(int basisType, int order);
- static Basis* tetGen(int basisType, int order);
- static Basis* hexGen(int basisType, int order);
+ static Basis* linZaglmayrGen(int basisType, int order);
+ static Basis* triZaglmayrGen(int basisType, int order);
+ static Basis* quaZaglmayrGen(int basisType, int order);
+ static Basis* tetZaglmayrGen(int basisType, int order);
+ static Basis* hexZaglmayrGen(int basisType, int order);
+
+ private:
+ static Basis* generateZaglmayr(int elementType,
+ int basisType,
+ int order);
+
+ static Basis* generateLagrange(int elementType,
+ int basisType,
+ int order);
};
@@ -50,8 +65,10 @@ class BasisGenerator{
on which the requested Basis will be created
@param basisType The Basis type
@param order The order or the requested Basis
+ @param family A string
- This method will @em instanciate the requested Basis
+ This method will @em instanciate the requested Basis,
+ of the requested family
@return Returns a @em pointer to a newly
@em instantiated Basis
@@ -68,9 +85,28 @@ class BasisGenerator{
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
+
+ @note Families are:
+ @li @c zaglmayr for
+ <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
+ Basis Functions
+ @li @c lagrange for Lagrange's Basis Functions
+ **
+
+ @fn BasisGenerator::generate
+ @param elementType The type of the element,
+ on which the requested Basis will be created
+ @param basisType The Basis type
+ @param order The order or the requested Basis
+
+ Same as
+ BasisGenerator::generate(@c elementType, @c basisType, @c order, @c zaglmayr)
+
+ @return Returns a @em pointer to a newly
+ @em instantiated Basis
**
- @fn BasisGenerator::linGen
+ @fn BasisGenerator::linZaglmayrGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -85,9 +121,11 @@ class BasisGenerator{
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
+
+ @note The Basis family will be @c zaglmayr
**
- @fn BasisGenerator::triGen
+ @fn BasisGenerator::triZaglmayrGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -102,9 +140,11 @@ class BasisGenerator{
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
+
+ @note The Basis family will be @c zaglmayr
**
- @fn BasisGenerator::quaGen
+ @fn BasisGenerator::quaZaglmayrGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -119,9 +159,11 @@ class BasisGenerator{
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
+
+ @note The Basis family will be @c zaglmayr
**
- @fn BasisGenerator::tetGen
+ @fn BasisGenerator::tetZaglmayrGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -136,9 +178,11 @@ class BasisGenerator{
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
+
+ @note The Basis family will be @c zaglmayr
**
- @fn BasisGenerator::hexGen
+ @fn BasisGenerator::hexZaglmayrGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -153,7 +197,23 @@ class BasisGenerator{
@li @c 1 for 1-Form
@li @c 2 for 2-Form
@li @c 3 for 3-Form
+
+ @note The Basis family will be @c zaglmayr
**
*/
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline Basis* BasisGenerator::generate(int elementType,
+ int basisType,
+ int order){
+
+ return BasisGenerator::generate(elementType,
+ basisType,
+ order,
+ "zaglmayr");
+}
+
#endif
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index d02168cc8d..a6c7d63d7c 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SRC
Legendre.cpp
Basis.cpp
+ LagrangeBasis.cpp
BasisScalar.cpp
BasisVector.cpp
BasisGenerator.cpp
@@ -16,9 +17,6 @@ set(SRC
CurlBasis.cpp
DivBasis.cpp
- LagrangeGenerator.cpp
- LagrangeBasis.cpp
- TriLagrangeBasis.cpp
LineNodeBasis.cpp
LineEdgeBasis.cpp
@@ -30,6 +28,7 @@ set(SRC
TriNodeBasis.cpp
TriEdgeBasis.cpp
TriNedelecBasis.cpp
+ TriLagrangeBasis.cpp
HexNodeBasis.cpp
HexEdgeBasis.cpp
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 9c973bc2af..93cea8990b 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -60,7 +60,7 @@ void FunctionSpace::build(const GroupOfElement& goe,
type = basisType;
basis = BasisGenerator::generate(elementType,
basisType,
- order);
+ order, "lagrange");
// Number of *Per* Entity functions //
fPerVertex = basis->getNVertexBased() / nVertex;
diff --git a/FunctionSpace/LagrangeBasis.cpp b/FunctionSpace/LagrangeBasis.cpp
index 9402df5f4c..4336a6ca16 100644
--- a/FunctionSpace/LagrangeBasis.cpp
+++ b/FunctionSpace/LagrangeBasis.cpp
@@ -58,3 +58,42 @@ project(const MElement& element,
// Return ;
return newCoef;
}
+
+unsigned int** LagrangeBasis::triEdgeOrder(int order){
+ // Check //
+ if(order <= 1)
+ return NULL;
+
+ // Number of Edge Based Basis //
+ const unsigned int orderMinus = order - 1;
+ const unsigned int size = 3 * orderMinus;
+
+ // Alloc //
+ unsigned int i;
+ unsigned int** edgeOrder = new unsigned int*[2];
+ edgeOrder[0] = new unsigned int[size];
+ edgeOrder[1] = new unsigned int[size];
+
+ // Direct Numeration //
+ i = 0;
+
+ for(unsigned int j = 0; j < orderMinus; j++){
+ for(unsigned int k = 0; k < size; k += orderMinus){
+ edgeOrder[0][i] = k + j;
+ i++;
+ }
+ }
+
+ // Invert Numeration //
+ i = 0;
+
+ for(unsigned int j = 0; j < orderMinus; j++){
+ for(unsigned int k = 0; k < size; k += orderMinus){
+ edgeOrder[1][i] = k + orderMinus - 1 - j;
+ i++;
+ }
+ }
+
+ // Retrun //
+ return edgeOrder;
+}
diff --git a/FunctionSpace/LagrangeBasis.h b/FunctionSpace/LagrangeBasis.h
index 0d990a576f..940040348b 100644
--- a/FunctionSpace/LagrangeBasis.h
+++ b/FunctionSpace/LagrangeBasis.h
@@ -67,6 +67,16 @@ class LagrangeBasis: public BasisScalar{
//! Returns a new LagrangeBasis
//!
LagrangeBasis(void);
+
+ //! @param order The ordre of the requested LagrangeBasis
+ //! @return Returns a 2D table
+ //! (of size 2 x (3 * (@c order - 1)))
+ //! with the Edge Triangle Inverted Closure for Lagrange
+ //! @note The first dimension of the returned table
+ //! is the direct closure, and the second dimension
+ //! the invert closure
+ //! @warning The returned table @em must be @em deleted
+ static unsigned int** triEdgeOrder(int order);
};
diff --git a/FunctionSpace/LagrangeGenerator.cpp b/FunctionSpace/LagrangeGenerator.cpp
deleted file mode 100644
index 99c8f0c960..0000000000
--- a/FunctionSpace/LagrangeGenerator.cpp
+++ /dev/null
@@ -1,22 +0,0 @@
-#include "Exception.h"
-#include "TriLagrangeBasis.h"
-#include "LagrangeGenerator.h"
-
-LagrangeGenerator::LagrangeGenerator(void){
-}
-
-LagrangeGenerator::~LagrangeGenerator(void){
-}
-
-LagrangeBasis* LagrangeGenerator::generate(unsigned int elementType,
- unsigned int lagrangeOrder){
- switch(elementType){
- case TYPE_TRI: return new TriLagrangeBasis(lagrangeOrder);
- case TYPE_QUA: throw Exception("Lagrange Basis on Quads not Implemented");
- case TYPE_TET: throw Exception("Lagrange Basis on Tets not Implemented");
- case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented");
-
- default: throw Exception("Unknown Element Type (%d) for Lagrange Basis",
- elementType);
- }
-}
diff --git a/FunctionSpace/LagrangeGenerator.h b/FunctionSpace/LagrangeGenerator.h
deleted file mode 100644
index 233d3d76da..0000000000
--- a/FunctionSpace/LagrangeGenerator.h
+++ /dev/null
@@ -1,52 +0,0 @@
-#ifndef _LAGRANGEGENERATOR_H_
-#define _LAGRANGEGENERATOR_H_
-
-#include "LagrangeBasis.h"
-
-/**
- @class LagrangeGenerator
- @brief Generates a Lagrange Basis
-
- Generates a Lagrange Basis
- */
-
-class LagrangeGenerator{
- public:
- LagrangeGenerator(void);
- ~LagrangeGenerator(void);
-
- static LagrangeBasis* generate(unsigned int elementType,
- unsigned int lagrangeOrder);
-};
-
-/**
- @fn LagrangeGenerator::LagrangeGenerator
- Instantiates a new LagrangeGenerator
-
- @note
- A LagrangeGenerator got @em only @em class @em methods,
- so it is not required to instanciate it.
- **
-
- @fn LagrangeGenerator::~LagrangeGenerator
- Deletes this LagrangeGenerator
- **
-
- @fn LagrangeGenerator::generate
- @param elementType The type of the element,
- on which the requested LagrangeBasis will be created
- @param lagrangeOrder The order or the requested LagrangeBasis
-
- This method will @em instanciate the requested LagrangeBasis
-
- @return Returns a @em pointer to a newly
- @em instantiated LagrangeBasis
-
- @note Element types are:
- @li @c TYPE_TRI for Triangles
- @li @c TYPE_QUA for Quadrangles
- @li @c TYPE_TET for Tetrahedrons
- @li @c TYPE_HEX for Hexahedrons
-*/
-
-#endif
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 99f8a6d7ca..093647adc9 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -66,27 +66,34 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
type = 0;
dim = 2;
- nVertex = l->coefficients.size1();
- nEdge = 0;
+ nVertex = 3;
+ nEdge = 3 * (order - 1);
nFace = 0;
- nCell = 0;
+ nCell = (order - 1) * (order - 2) / 2;
- nEdgeClosure = 0;
+ nEdgeClosure = 2;
nFaceClosure = 0;
- size = nVertex;
-
+ size = nVertex + nEdge + nFace + nCell;
+
+
+ // Alloc Some Stuff //
+ const int nMonomial = l->monomials.size1();
+ unsigned int** edgeOrder;
+ Polynomial* pEdgeClosureLess = new Polynomial[nEdge];
- // Basis (pure nodal) //
+
+ // Basis //
node = new vector<Polynomial*>(nVertex);
- edge = new vector<vector<Polynomial*>*>(0);
+ edge = new vector<vector<Polynomial*>*>(2);
face = new vector<vector<Polynomial*>*>(0);
- cell = new vector<Polynomial*>(0);
+ cell = new vector<Polynomial*>(nCell);
+ (*edge)[0] = new vector<Polynomial*>(nEdge);
+ (*edge)[1] = new vector<Polynomial*>(nEdge);
- // Instanciate Polynomials //
- const int nMonomial = l->monomials.size1();
+ // Vertex Based //
for(int i = 0; i < nVertex; i++){
Polynomial p = Polynomial(0, 0, 0, 0);
@@ -98,6 +105,58 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
(*node)[i] = new Polynomial(p);
}
+
+
+ // Edge Based //
+ // Without Closure
+ for(int i = 0; i < nEdge; i++){
+ int ci = i + nVertex;
+ pEdgeClosureLess[i] = Polynomial(0, 0, 0, 0);
+
+ for(int j = 0; j < nMonomial; j++)
+ pEdgeClosureLess[i] =
+ pEdgeClosureLess[i] +
+ Polynomial(l->coefficients(ci, j), // Coef
+ l->monomials(j, 0), // powerX
+ l->monomials(j, 1), // powerY
+ 0); // powerZ
+ }
+
+ // With Closure
+ edgeOrder = triEdgeOrder(order); // Closure Ordering
+
+ for(int i = 0; i < nEdge; i++){
+ (*(*edge)[0])[i] =
+ new Polynomial(pEdgeClosureLess[edgeOrder[0][i]]);
+
+ (*(*edge)[1])[i] =
+ new Polynomial(pEdgeClosureLess[edgeOrder[1][i]]);
+ }
+
+
+ // Cell Based //
+ for(int i = 0; i < nCell; i++){
+ int ci = i + nVertex + nEdge;
+ Polynomial p = Polynomial(0, 0, 0, 0);
+
+ for(int j = 0; j < nMonomial; j++)
+ p = p + Polynomial(l->coefficients(ci, j), // Coef
+ l->monomials(j, 0), // powerX
+ l->monomials(j, 1), // powerY
+ 0); // powerZ
+
+ (*cell)[i] = new Polynomial(p);
+ }
+
+
+ // Delete Temporary Space //
+ delete[] pEdgeClosureLess;
+
+ if(edgeOrder){
+ delete[] edgeOrder[0];
+ delete[] edgeOrder[1];
+ delete[] edgeOrder;
+ }
}
TriLagrangeBasis::~TriLagrangeBasis(void){
--
GitLab