BasisLagrange::projection -- BasisGenerator: zaglmayr family becomes...

BasisLagrange::projection -- BasisGenerator: zaglmayr family becomes hierarchical family + Lagrange works again

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);

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

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;
+}

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

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;

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;
 }
-*/