Interpolation on Edge -- Projection OK -- Yeahaaaaa \o/ <o> \o/ -- need some cleaning

FunctionSpace/BasisGenerator.cpp

@@ -35,11 +35,11 @@ Basis* BasisGenerator::generate(int elementType,
 Basis* BasisGenerator::TriGen(int basisType, 
 			      int order){
   switch(basisType){
-  case  0: 
+  case  0: return new TriNodeBasis(order);
+  case  1: 
     if (order == 0) return new TriNedelecBasis();
-    else            return new TriNodeBasis(order);
+    else            return new TriEdgeBasis(order);
 
-  case  1: return new TriEdgeBasis(order);
   case  2: throw Exception("2-form not implemented on Triangles");
   case  3: throw Exception("3-form not implemented on Triangles");
 

FunctionSpace/CMakeLists.txt

@@ -22,6 +22,10 @@ set(SRC
   HexEdgeBasis.cpp
 
   FunctionSpace.cpp
+  FunctionSpaceScalar.cpp
+  FunctionSpaceVector.cpp
+  FunctionSpaceNode.cpp
+  FunctionSpaceEdge.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) 

FunctionSpace/FunctionSpace.cpp

@@ -3,8 +3,15 @@
 
 using namespace std;
 
-FunctionSpace::FunctionSpace(const GroupOfElement& goe,
-			     int basisType, int order){
+FunctionSpace::FunctionSpace(void){
+}
+
+FunctionSpace::~FunctionSpace(void){
+  delete basis;
+}
+
+void FunctionSpace::build(const GroupOfElement& goe,
+			  int basisType, int order){
   // Save GroupOfElement & Mesh //
   this->goe  = &goe;
   this->mesh = &(goe.getMesh());
@@ -41,10 +48,6 @@ FunctionSpace::FunctionSpace(const GroupOfElement& goe,
   fPerCell = basis->getNCellBased(); // We always got 1 cell 
 }
 
-FunctionSpace::~FunctionSpace(void){
-  delete basis;
-}
-
 vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{ 
   // Const_Cast //
   MElement& element = const_cast<MElement&>(elem);
@@ -94,9 +97,53 @@ vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{
     }
   }
   
+  // Add Edge Based Dof //
+  for(int i = 0; i < nEdge; i++){
+    for(int j = 0; j < nFEdge; j++){
+      myDof[it].setDof(mesh->getGlobalId(edge[i]), j);
+      it++;
+    }
+  }
+  /*
+  // Add Face Based Dof //
+  for(int i = 0; i < nFace; i++){
+    for(int j = 0; j < nFFace; j++){
+      myDof[it].setDof(mesh->getGlobalId(face[i]), j);
+      it++;
+    }
+  }
+  */
+  // Add Cell Based Dof //
+  for(int j = 0; j < nFCell; j++){
+    myDof[it].setDof(mesh->getGlobalId(element), j);
+    it++;
+  }
+  
+  
   return myDof;
 }
 
+int FunctionSpace::getElementType(const Dof& dof) const{
+  const unsigned int type = dof.getType();
+
+  if(type < fPerVertex)    // Vertex Based
+    return 0; 
+
+  else if(type < fPerEdge) // Edge Based
+    return 1;
+
+  else if(type < fPerFace) // Face Based
+    return 2;
+
+  else                     // Cell Based
+    return 3; 
+}
+
+int FunctionSpace::getElementGlobalId(const Dof& dof) const{
+  return dof.getEntity();
+}
+
+
 /*
 #include "Polynomial.h"
 #include "BasisScalar.h"

FunctionSpace/FunctionSpace.h

@@ -30,37 +30,37 @@ class FunctionSpace{
   const GroupOfElement* goe;
   const Basis*          basis;
 
-  int fPerVertex;
-  int fPerEdge;
-  int fPerFace;
-  int fPerCell;
-  int type;
+  unsigned int fPerVertex;
+  unsigned int fPerEdge;
+  unsigned int fPerFace;
+  unsigned int fPerCell;
+  unsigned int type;
 
  public:
-  FunctionSpace(const GroupOfElement& goe,
-		int basisType, int order);
-
   virtual ~FunctionSpace(void);
 
   const GroupOfElement& getSupport(void) const;
   const Basis&          getBasis(const MElement& element) const;
   int                   getType(void) const;
 
-  int getNFunctionPerVertex(const MElement& element) const;
-  int getNFunctionPerEdge(const MElement& element) const;
-  int getNFunctionPerFace(const MElement& element) const;
-  int getNFunctionPerCell(const MElement& element) const;
+  unsigned int getNFunctionPerVertex(const MElement& element) const;
+  unsigned int getNFunctionPerEdge(const MElement& element) const;
+  unsigned int getNFunctionPerFace(const MElement& element) const;
+  unsigned int getNFunctionPerCell(const MElement& element) const;
 
   std::vector<Dof> getKeys(const MElement& element) const;
   std::vector<Dof> getKeys(const MVertex& element) const;
   std::vector<Dof> getKeys(const MEdge& element) const;
   std::vector<Dof> getKeys(const MFace& element) const;
 
-  /*
-  void interpolateAtNodes(const MElement& element, 
-			  const std::vector<double>& coef,
-			  std::vector<double>& nodeValue) const;
-  */
+  int getElementType(const Dof& dof) const;
+  int getElementGlobalId(const Dof& dof) const;
+
+ protected:
+  FunctionSpace();
+
+  void build(const GroupOfElement& goe,
+	     int basisType, int order);
 };
 
 //////////////////////
@@ -79,19 +79,19 @@ inline int FunctionSpace::getType(void) const{
   return type;
 }
 
-inline int FunctionSpace::getNFunctionPerVertex(const MElement& element) const{
+inline unsigned int FunctionSpace::getNFunctionPerVertex(const MElement& element) const{
   return fPerVertex;
 }
 
-inline int FunctionSpace::getNFunctionPerEdge(const MElement& element) const{
+inline unsigned int FunctionSpace::getNFunctionPerEdge(const MElement& element) const{
   return fPerEdge;
 }
 
-inline int FunctionSpace::getNFunctionPerFace(const MElement& element) const{
+inline unsigned int FunctionSpace::getNFunctionPerFace(const MElement& element) const{
   return fPerFace;
 }
 
-inline int FunctionSpace::getNFunctionPerCell(const MElement& element) const{
+inline unsigned int FunctionSpace::getNFunctionPerCell(const MElement& element) const{
   return fPerCell;
 }
 

FunctionSpace/FunctionSpaceEdge.cpp

+#include <vector>
+
+#include "BasisVector.h"
+#include "Mapper.h"
+#include "FunctionSpaceEdge.h"
+
+using namespace std;
+
+FunctionSpaceEdge::FunctionSpaceEdge(const GroupOfElement& goe,
+				     int order){
+  // Build 1Form Basis //
+  build(goe, 1, order); 
+}
+    
+FunctionSpaceEdge::~FunctionSpaceEdge(void){
+}
+
+#include <cstdio>
+
+fullVector<double> FunctionSpaceEdge::
+interpolate(const MElement& element, 
+	    const std::vector<double>& coef,
+	    const fullVector<double>& xyz) const{
+ 
+  MElement& eelement = 
+    const_cast<MElement&>(element);
+  
+  const vector<vector<Polynomial> >& fun = 
+    static_cast<const BasisVector&>(*basis).getBasis();
+  
+  const unsigned int nFun = fun.size();
+
+  fullVector<double> val(3);
+
+  fullMatrix<double>  invJac(3, 3);        
+  eelement.getJacobian(xyz(0), xyz(1), xyz(2), invJac);
+  invJac.invertInPlace();
+
+  fullVector<double> origin(3);
+  origin(0) = eelement.getVertex(0)->x();
+  origin(1) = eelement.getVertex(0)->y();
+  origin(2) = eelement.getVertex(0)->z();
+
+  fullVector<double> uvw = Mapper::invMap(xyz, origin, invJac);
+  
+  printf("%g, %g, %g\n", uvw(0), uvw(1), uvw(2));
+  
+  //val(0) = 1; val(1) = 1; val(2) = 0;
+  
+  
+  for(unsigned int i = 0; i < nFun; i++){
+    fullVector<double> vi = 
+      Mapper::grad(Polynomial::at(fun[i], uvw(0), uvw(1), uvw(2)),
+		   invJac);
+    
+    vi.scale(coef[i]);
+    val.axpy(vi, 1);
+  }
+
+  return val;
+}

FunctionSpace/FunctionSpaceEdge.h

+#ifndef _FUNCTIONSPACEEDGE_H_
+#define _FUNCTIONSPACEEDGE_H_
+
+#include "FunctionSpaceVector.h"
+
+class FunctionSpaceEdge : public FunctionSpaceVector{
+ public:
+  FunctionSpaceEdge(const GroupOfElement& goe, int order);
+    
+  virtual ~FunctionSpaceEdge(void);
+
+  virtual fullVector<double> 
+    interpolate(const MElement& element, 
+		const std::vector<double>& coef,
+		const fullVector<double>& xyz) const;
+};
+
+#endif

FunctionSpace/FunctionSpaceNode.cpp

+#include "FunctionSpaceNode.h"
+
+FunctionSpaceNode::FunctionSpaceNode(const GroupOfElement& goe,
+				     int order){
+  // Build 0Form Basis //
+  build(goe, 0, order); 
+}
+    
+FunctionSpaceNode::~FunctionSpaceNode(void){
+}
+
+double FunctionSpaceNode::
+interpolate(const MElement& element, 
+	    const std::vector<double>& coef,
+	    const fullVector<double>& xyz) const{
+
+  return 42;
+}

FunctionSpace/FunctionSpaceNode.h

+#ifndef _FUNCTIONSPACENODE_H_
+#define _FUNCTIONSPACENODE_H_
+
+#include "FunctionSpaceScalar.h"
+
+class FunctionSpaceNode : public FunctionSpaceScalar{
+ public:
+  FunctionSpaceNode(const GroupOfElement& goe, int order);
+    
+  virtual ~FunctionSpaceNode(void);
+
+  virtual double 
+    interpolate(const MElement& element, 
+		const std::vector<double>& coef,
+		const fullVector<double>& xyz) const;
+};
+
+#endif

FunctionSpace/FunctionSpaceScalar.cpp

+#include "FunctionSpaceScalar.h"
+
+FunctionSpaceScalar::~FunctionSpaceScalar(void){
+}

FunctionSpace/FunctionSpaceScalar.h

+#ifndef _FUNCTIONSPACESCALAR_H_
+#define _FUNCTIONSPACESCALAR_H_
+
+#include "fullMatrix.h"
+#include "FunctionSpace.h"
+
+class FunctionSpaceScalar : public FunctionSpace{
+ public:
+  virtual ~FunctionSpaceScalar(void);
+
+  virtual double 
+    interpolate(const MElement& element, 
+		const std::vector<double>& coef,
+		const fullVector<double>& xyz) const = 0;
+};
+
+#endif

FunctionSpace/FunctionSpaceVector.cpp

+#include "FunctionSpaceVector.h"
+
+FunctionSpaceVector::~FunctionSpaceVector(void){
+}

FunctionSpace/FunctionSpaceVector.h

+#ifndef _FUNCTIONSPACEVECTOR_H_
+#define _FUNCTIONSPACEVECTOR_H_
+
+#include "fullMatrix.h"
+#include "FunctionSpace.h"
+
+class FunctionSpaceVector : public FunctionSpace{
+ public:
+  virtual ~FunctionSpaceVector(void);
+
+  virtual fullVector<double> 
+    interpolate(const MElement& element, 
+		const std::vector<double>& coef,
+		const fullVector<double>& xyz) const = 0;
+};
+
+#endif