Solution & interpolation -- BAD (but works)

FunctionSpace/CMakeLists.txt

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

FunctionSpace/FunctionSpace.cpp

@@ -129,3 +129,70 @@ FunctionSpace::~FunctionSpace(void){
   delete basis;
 }
 
+/*
+#include "Polynomial.h"
+#include "BasisScalar.h"
+#include "fullMatrix.h"
+#include "Mapper.h"
+#include "Exception.h"
+void FunctionSpace::
+interpolateAtNodes(const MElement& elem, 
+		   const vector<double>& coef,
+		   std::vector<double>& nodeValue) const{
+  // Check
+  unsigned int wS = coef.size();
+  unsigned int bS = basis->getSize();
+
+  if(wS < bS)
+    throw Exception
+      ("Not enough coefs for interpolation:\nBasis: %d -- coefs: %d",
+       bS, wS);
+
+  if(wS > bS)
+    throw Exception
+      ("Too much coefs for interpolation:\nBasis: %d -- coefs: %d",
+       bS, wS);
+
+  // Get Nodes
+  MElement& element = const_cast<MElement&>(elem);
+
+  vector<MVertex*> node;
+  element.getVertices(node);
+  unsigned int N = node.size();
+
+  // Get Functions
+  const vector<Polynomial>& fun = 
+    static_cast<const BasisScalar*>(basis)->getBasis();
+
+  // Init some stuff
+  fullMatrix<double> invJac(3, 3);
+  fullVector<double> xyz(3);
+  fullVector<double> origin(3);
+  
+  origin(0) = node[0]->x();
+  origin(1) = node[0]->y();
+  origin(2) = node[0]->z();
+
+  // Interpolate
+  for(unsigned int n = 0; n < N; n++){
+    // Map from physical to reference space 
+    xyz(0) = node[n]->x();
+    xyz(1) = node[n]->y();
+    xyz(2) = node[n]->z();
+
+    element.getJacobian(xyz(0), xyz(1), xyz(2), invJac);
+    invJac.invertInPlace();
+    
+    const fullVector<double> uvw = 
+      Mapper::invMap(xyz, origin, invJac);
+
+    printf("(%lf\t%lf\t%lf)\n", uvw(0), uvw(1), uvw(2));
+
+    // Interpolate
+    const int id = node[n]->getNum() - 1;
+
+    for(unsigned int i = 0; i < bS; i++)
+      nodeValue[id] += fun[i].at(uvw(0), uvw(1), uvw(2)) * coef[i];
+  }
+}
+*/

FunctionSpace/FunctionSpace.h

@@ -52,6 +52,12 @@ class FunctionSpace{
   int getNFunctionPerEdge(const MElement& element) const;
   int getNFunctionPerFace(const MElement& element) const;
   int getNFunctionPerCell(const MElement& element) const;
+
+  /*
+  void interpolateAtNodes(const MElement& element, 
+			  const std::vector<double>& coef,
+			  std::vector<double>& nodeValue) const;
+  */
 };
 
 //////////////////////

FunctionSpace/FunctionSpaceNode.cpp

-#include "FunctionSpaceNode.h"
-#include "Polynomial.h"
-#include "BasisScalar.h"
-#include "fullMatrix.h"
-#include "Mapper.h"
-#include "Exception.h"
-
-using namespace std;
-
-FunctionSpaceNode::FunctionSpaceNode(const GroupOfElement& goe, int order):
-FunctionSpaceScalar(goe, 0, order){
-  // Just Call Super Constructor with basisType = 0
-}
-
-
-FunctionSpaceNode::~FunctionSpaceNode(void){
-}
-
-void FunctionSpaceNode::
-interpolateAtNodes(const MElement& elem, 
-		   const vector<double>& coef,
-		   std::vector<double>& nodeValue) const{
-  // Check
-  unsigned int wS = coef.size();
-  unsigned int bS = basis->getSize();
-
-  if(wS < bS)
-    throw Exception
-      ("Not enough coefs for interpolation:\nBasis: %d -- coefs: %d",
-       bS, wS);
-
-  if(wS > bS)
-    throw Exception
-      ("Too much coefs for interpolation:\nBasis: %d -- coefs: %d",
-       bS, wS);
-
-  // Get Nodes
-  MElement& element = const_cast<MElement&>(elem);
-
-  vector<MVertex*> node;
-  element.getVertices(node);
-  unsigned int N = node.size();
-
-  // Get Functions
-  const vector<Polynomial>& fun = 
-    static_cast<const BasisScalar*>(basis)->getBasis();
-
-  // Init some stuff
-  fullMatrix<double> invJac(3, 3);
-  fullVector<double> xyz(3);
-  fullVector<double> origin(3);
-  
-  origin(0) = node[0]->x();
-  origin(1) = node[0]->y();
-  origin(2) = node[0]->z();
-
-  // Interpolate
-  for(unsigned int n = 0; n < N; n++){
-    // Map from physical to reference space 
-    xyz(0) = node[n]->x();
-    xyz(1) = node[n]->y();
-    xyz(2) = node[n]->z();
-
-    element.getJacobian(xyz(0), xyz(1), xyz(2), invJac);
-    invJac.invertInPlace();
-    
-    const fullVector<double> uvw = 
-      Mapper::invMap(xyz, origin, invJac);
-
-
-    // Interpolate
-    const int id = node[n]->getNum();
-
-    for(unsigned int i = 0; i < bS; i++)
-      nodeValue[id] += fun[i].at(uvw(0), uvw(1), uvw(2)) * coef[i];
-  }
-}

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 void interpolateAtNodes(const MElement& element, 
-				  const std::vector<double>& coef,
-				  std::vector<double>& nodeValue) const;
-};
-
-#endif

FunctionSpace/FunctionSpaceScalar.cpp

-#include "FunctionSpaceScalar.h"
-
-FunctionSpaceScalar::FunctionSpaceScalar(const GroupOfElement& goe, 
-					  int basisType, int order):
-FunctionSpace(goe, basisType, order){
-}
-
-FunctionSpaceScalar::~FunctionSpaceScalar(void){
-}

FunctionSpace/FunctionSpaceScalar.h

-#ifndef _FUNCTIONSPACESCALAR_H_
-#define _FUNCTIONSPACESCALAR_H_
-
-#include "FunctionSpace.h"
-
-class FunctionSpaceScalar: public FunctionSpace{
- public:
-  FunctionSpaceScalar(const GroupOfElement& goe, 
-		      int basisType, int order);
-
-  virtual ~FunctionSpaceScalar(void);
-
-  virtual void interpolateAtNodes(const MElement& element, 
-				  const std::vector<double>& coef,
-				  std::vector<double>& nodeValue) const = 0;
-};
-
-#endif