From dfc425fb7562a3fb5a7c20eff7e140adc400bdbb Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 15 Jan 2013 10:23:28 +0000
Subject: [PATCH] BasisLagrange::projection -- BasisGenerator: zaglmayr family
becomes hierarchical family + Lagrange works again
---
FunctionSpace/BasisGenerator.cpp | 46 +++++++++++-----------
FunctionSpace/BasisGenerator.h | 45 +++++++++++----------
FunctionSpace/BasisLagrange.cpp | 52 +++++++++++++++++++++++-
FunctionSpace/BasisLagrange.h | 63 +++++++++++++++++++++++++++++-
FunctionSpace/FunctionSpace.cpp | 2 +-
FunctionSpace/TriLagrangeBasis.cpp | 10 +++--
6 files changed, 165 insertions(+), 53 deletions(-)
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index 94acc23326..7f0626745e 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -12,7 +12,7 @@
#include "TriNodeBasis.h"
#include "TriEdgeBasis.h"
#include "TriNedelecBasis.h"
-//#include "TriLagrangeBasis.h"
+#include "TriLagrangeBasis.h"
#include "TetNodeBasis.h"
#include "TetEdgeBasis.h"
@@ -32,31 +32,31 @@ Basis* BasisGenerator::generate(unsigned int elementType,
unsigned int order,
std::string family){
- if(!family.compare(std::string("zaglmayr")))
- return generateZaglmayr(elementType, basisType, order);
- /*
+ if(!family.compare(std::string("hierarchical")))
+ return generateHierarchical(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(unsigned int elementType,
+Basis* BasisGenerator::generateHierarchical(unsigned int elementType,
unsigned int basisType,
unsigned 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);
+ case TYPE_LIN: return linHierarchicalGen(basisType, order);
+ case TYPE_TRI: return triHierarchicalGen(basisType, order);
+ case TYPE_QUA: return quaHierarchicalGen(basisType, order);
+ case TYPE_TET: return tetHierarchicalGen(basisType, order);
+ case TYPE_HEX: return hexHierarchicalGen(basisType, order);
default: throw Exception("Unknown Element Type (%d) for Basis Generation",
elementType);
}
}
-/*
+
Basis* BasisGenerator::generateLagrange(unsigned int elementType,
unsigned int basisType,
unsigned int order){
@@ -76,9 +76,9 @@ Basis* BasisGenerator::generateLagrange(unsigned int elementType,
elementType);
}
}
-*/
-Basis* BasisGenerator::linZaglmayrGen(unsigned int basisType,
- unsigned int order){
+
+Basis* BasisGenerator::linHierarchicalGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
case 0: return new LineNodeBasis(order);
case 1:
@@ -92,8 +92,8 @@ Basis* BasisGenerator::linZaglmayrGen(unsigned int basisType,
}
}
-Basis* BasisGenerator::triZaglmayrGen(unsigned int basisType,
- unsigned int order){
+Basis* BasisGenerator::triHierarchicalGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
case 0: return new TriNodeBasis(order);
case 1:
@@ -107,8 +107,8 @@ Basis* BasisGenerator::triZaglmayrGen(unsigned int basisType,
}
}
-Basis* BasisGenerator::quaZaglmayrGen(unsigned int basisType,
- unsigned int order){
+Basis* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
//case 0: return new QuadNodeBasis(order);
//case 1: return new QuadEdgeBasis(order);
@@ -119,8 +119,8 @@ Basis* BasisGenerator::quaZaglmayrGen(unsigned int basisType,
}
}
-Basis* BasisGenerator::tetZaglmayrGen(unsigned int basisType,
- unsigned int order){
+Basis* BasisGenerator::tetHierarchicalGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
case 0: return new TetNodeBasis(order);
case 1: return new TetEdgeBasis(order);
@@ -131,8 +131,8 @@ Basis* BasisGenerator::tetZaglmayrGen(unsigned int basisType,
}
}
-Basis* BasisGenerator::hexZaglmayrGen(unsigned int basisType,
- unsigned int order){
+Basis* BasisGenerator::hexHierarchicalGen(unsigned int basisType,
+ unsigned 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 ee2584d2f8..453d99a526 100644
--- a/FunctionSpace/BasisGenerator.h
+++ b/FunctionSpace/BasisGenerator.h
@@ -30,21 +30,20 @@ class BasisGenerator{
unsigned int basisType,
unsigned int order);
- static Basis* linZaglmayrGen(unsigned int basisType, unsigned int order);
- static Basis* triZaglmayrGen(unsigned int basisType, unsigned int order);
- static Basis* quaZaglmayrGen(unsigned int basisType, unsigned int order);
- static Basis* tetZaglmayrGen(unsigned int basisType, unsigned int order);
- static Basis* hexZaglmayrGen(unsigned int basisType, unsigned int order);
+ static Basis* linHierarchicalGen(unsigned int basisType, unsigned int order);
+ static Basis* triHierarchicalGen(unsigned int basisType, unsigned int order);
+ static Basis* quaHierarchicalGen(unsigned int basisType, unsigned int order);
+ static Basis* tetHierarchicalGen(unsigned int basisType, unsigned int order);
+ static Basis* hexHierarchicalGen(unsigned int basisType, unsigned int order);
private:
- static Basis* generateZaglmayr(unsigned int elementType,
- unsigned int basisType,
- unsigned int order);
- /*
+ static Basis* generateHierarchical(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order);
+
static Basis* generateLagrange(unsigned int elementType,
unsigned int basisType,
unsigned int order);
- */
};
@@ -88,7 +87,7 @@ class BasisGenerator{
@li @c 3 for 3-Form
@note Families are:
- @li @c zaglmayr for
+ @li @c hierarchical for
<a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
Basis Functions
@li @c lagrange for Lagrange's Basis Functions
@@ -101,13 +100,13 @@ class BasisGenerator{
@param order The order or the requested Basis
Same as
- BasisGenerator::generate(@c elementType, @c basisType, @c order, @c zaglmayr)
+ BasisGenerator::generate(@c elementType, @c basisType, @c order, @c hierarchical)
@return Returns a @em pointer to a newly
@em instantiated Basis
**
- @fn BasisGenerator::linZaglmayrGen
+ @fn BasisGenerator::linHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -123,10 +122,10 @@ class BasisGenerator{
@li @c 2 for 2-Form
@li @c 3 for 3-Form
- @note The Basis family will be @c zaglmayr
+ @note The Basis family will be @c hierarchical
**
- @fn BasisGenerator::triZaglmayrGen
+ @fn BasisGenerator::triHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -142,10 +141,10 @@ class BasisGenerator{
@li @c 2 for 2-Form
@li @c 3 for 3-Form
- @note The Basis family will be @c zaglmayr
+ @note The Basis family will be @c hierarchical
**
- @fn BasisGenerator::quaZaglmayrGen
+ @fn BasisGenerator::quaHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -161,10 +160,10 @@ class BasisGenerator{
@li @c 2 for 2-Form
@li @c 3 for 3-Form
- @note The Basis family will be @c zaglmayr
+ @note The Basis family will be @c hierarchical
**
- @fn BasisGenerator::tetZaglmayrGen
+ @fn BasisGenerator::tetHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -180,10 +179,10 @@ class BasisGenerator{
@li @c 2 for 2-Form
@li @c 3 for 3-Form
- @note The Basis family will be @c zaglmayr
+ @note The Basis family will be @c hierarchical
**
- @fn BasisGenerator::hexZaglmayrGen
+ @fn BasisGenerator::hexHierarchicalGen
@param basisType The Basis type
@param order The order or the requested Basis
@@ -199,7 +198,7 @@ class BasisGenerator{
@li @c 2 for 2-Form
@li @c 3 for 3-Form
- @note The Basis family will be @c zaglmayr
+ @note The Basis family will be @c hierarchical
**
*/
@@ -214,7 +213,7 @@ inline Basis* BasisGenerator::generate(unsigned int elementType,
return BasisGenerator::generate(elementType,
basisType,
order,
- "zaglmayr");
+ "hierarchical");
}
#endif
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index 47f1ae99f2..8c7db4dfd2 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -25,7 +25,7 @@ fullMatrix<double>* BasisLagrange::getFunctions(unsigned int permutation,
point(0, 1) = v;
point(0, 2) = w;
- basis->f(point, tmp);
+ lBasis->f(point, tmp);
*values = tmp.transpose(); // Otherwise not coherent with df !!
@@ -48,3 +48,53 @@ const fullMatrix<double>&
BasisLagrange::getPreEvaluatedGradFunctions(const MElement& element) const{
return fullMatrix<double>(nFunction, 1);
}
+
+std::vector<double> BasisLagrange::
+project(const MElement& element,
+ const std::vector<double>& coef,
+ const FunctionSpaceScalar& fSpace){
+
+ // Init New Coefs //
+ const unsigned int size = lPoint->size1();
+ std::vector<double> newCoef(size);
+
+ // Interpolation at Lagrange Points //
+ for(unsigned int i = 0; i < size; i++){
+ fullVector<double> uvw(3);
+ uvw(0) = (*lPoint)(i, 0);
+ uvw(1) = (*lPoint)(i, 1);
+ uvw(2) = (*lPoint)(i, 2);
+
+ newCoef[i] = fSpace.interpolateInRefSpace(element,
+ coef,
+ uvw);
+ }
+
+ // Return ;
+ return newCoef;
+}
+
+std::vector<fullVector<double> > BasisLagrange::
+project(const MElement& element,
+ const std::vector<double>& coef,
+ const FunctionSpaceVector& fSpace){
+
+ // Init New Coefs //
+ const unsigned int size = lPoint->size1();
+ std::vector<fullVector<double> > newCoef(size);
+
+ // Interpolation at Lagrange Points //
+ for(unsigned int i = 0; i < size; i++){
+ fullVector<double> uvw(3);
+ uvw(0) = (*lPoint)(i, 0);
+ uvw(1) = (*lPoint)(i, 1);
+ uvw(2) = (*lPoint)(i, 2);
+
+ newCoef[i] = fSpace.interpolateInRefSpace(element,
+ coef,
+ uvw);
+ }
+
+ // Return ;
+ return newCoef;
+}
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index ee9fa43bf4..7554ade178 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -2,11 +2,31 @@
#define _BASISLAGRANGE_H_
#include "BasisScalar.h"
+#include "FunctionSpaceScalar.h"
+#include "FunctionSpaceVector.h"
+#include "fullMatrix.h"
#include "polynomialBasis.h"
+/**
+ @interface BasisLagrange
+ @brief Interface for Lagrange Basis
+
+ This is an interface for Lagrange Basis.@n
+
+ These Scalar Basis allow a @em Coefficient Matrix,
+ and a Monomial Matrix, to be consulted.@n
+
+ A vector from an Other Basis (set of Functions)
+ can also be projected into a Lagrange Basis.@n
+
+ @todo
+ Add a method to get lagrange Point in polynomialBasis
+*/
+
class BasisLagrange: public BasisScalar{
protected:
- polynomialBasis* basis;
+ polynomialBasis* lBasis; // Lagrange Basis
+ fullMatrix<double>* lPoint; // Lagrange Points
public:
virtual ~BasisLagrange(void);
@@ -26,8 +46,49 @@ class BasisLagrange: public BasisScalar{
virtual const fullMatrix<double>&
getPreEvaluatedGradFunctions(const MElement& element) const;
+ //! @return Returns the Coefficient Matrix
+ const fullMatrix<double>& getCoefficient(void) const;
+
+ //! @return Returns the Monomial Matrix
+ const fullMatrix<double>& getMonomial(void) const;
+
+ //! @param element A MElement
+ //! @param coef A vector of coefficient associated
+ //! to the given Element
+ //! @param fSpace The (scalar) Function Space
+ //! of the given Coefficients
+ //! @return Projects the given Coefficients in this BasisLagrange
+ std::vector<double> project(const MElement& element,
+ const std::vector<double>& coef,
+ const FunctionSpaceScalar& fSpace);
+
+ //! @param element A MElement
+ //! @param coef A vector of coefficient associated
+ //! to the given Element
+ //! @param fSpace The (vectorial) Function Space
+ //! of the given Coefficients
+ //! @return Projects the given Coefficients in this BasisLagrange
+ std::vector<fullVector<double> >
+ project(const MElement& element,
+ const std::vector<double>& coef,
+ const FunctionSpaceVector& fSpace);
+
protected:
BasisLagrange(void);
};
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline const fullMatrix<double>& BasisLagrange::
+getCoefficient(void) const{
+ return lBasis->coefficients;
+}
+
+inline const fullMatrix<double>& BasisLagrange::
+getMonomial(void) const{
+ return lBasis->monomials;
+}
+
#endif
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 713be01704..154b0ac8c6 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, "zaglmayr");
+ order, "hierarchical");
// Number of *Per* Entity functions //
fPerVertex = basis->getNVertexBased() / nVertex;
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 40192124c8..99ff0469da 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -15,11 +15,15 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
nFunction = nVertex + nEdge + nFace + nCell;
// Init polynomialBasis //
- basis = new polynomialBasis(getTag(order));
+ lBasis = new polynomialBasis(getTag(order));
+
+ // Init Lagrange Point //
+ lPoint = new fullMatrix<double>(triPoint(order));
}
TriLagrangeBasis::~TriLagrangeBasis(void){
- delete basis;
+ delete lBasis;
+ delete lPoint;
}
unsigned int TriLagrangeBasis::getTag(unsigned int order){
@@ -42,7 +46,6 @@ unsigned int TriLagrangeBasis::getTag(unsigned int order){
}
}
-/*
fullMatrix<double> TriLagrangeBasis::
triPoint(unsigned int order){
unsigned int nbPoints = (order + 1) * (order + 2) / 2;
@@ -111,4 +114,3 @@ triPoint(unsigned int order){
return point;
}
-*/
--
GitLab