Interpolation -- Node

FunctionSpace/BasisTest.cpp

@@ -25,7 +25,8 @@ int ain(int argc, char** argv){
   Mesh msh(argv[1]);
 
   // Writer for .msh
-  WriterMsh writer(msh.getGroup(2).getAll());
+  WriterMsh writer;
+  writer.setDomain(msh.getGroup(2).getAll());
 
   // Plot Basis //
   HexNodeBasis b(1);

FunctionSpace/CMakeLists.txt

@@ -22,6 +22,7 @@ set(SRC
   HexEdgeBasis.cpp
 
   FunctionSpace.cpp
+  FunctionSpaceScalar.cpp
   FunctionSpaceNode.cpp
 )
 

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):
-FunctionSpace(goe, 0, order){
+FunctionSpaceScalar(goe, 0, order){
   // Just Call Super Constructor with basisType = 0
 }
 
@@ -14,9 +16,10 @@ FunctionSpace(goe, 0, order){
 FunctionSpaceNode::~FunctionSpaceNode(void){
 }
 
-vector<double> FunctionSpaceNode::
+void FunctionSpaceNode::
 interpolateAtNodes(const MElement& elem, 
-		   const vector<double>& coef) const{
+		   const vector<double>& coef,
+		   std::vector<double>& nodeValue) const{
   // Check
   unsigned int wS = coef.size();
   unsigned int bS = basis->getSize();
@@ -42,20 +45,33 @@ interpolateAtNodes(const MElement& elem,
   const vector<Polynomial>& fun = 
     static_cast<const BasisScalar*>(basis)->getBasis();
 
-  // Init Vector
-  vector<double> value(N);
+  // 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++){
-    const double x = node[n]->x();
-    const double y = node[n]->y();
-    const double z = node[n]->z();
+    // Map from physical to reference space 
+    xyz(0) = node[n]->x();
+    xyz(1) = node[n]->y();
+    xyz(2) = node[n]->z();
 
-    value[n] = 0;
+    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++)
-      value[n] += fun[i].at(x, y, z) * coef[i];
+      nodeValue[id] += fun[i].at(uvw(0), uvw(1), uvw(2)) * coef[i];
   }
-
-  return value;
 }

FunctionSpace/FunctionSpaceNode.h

 #ifndef _FUNCTIONSPACENODE_H_
 #define _FUNCTIONSPACENODE_H_
 
-#include "FunctionSpace.h"
+#include "FunctionSpaceScalar.h"
 
-class FunctionSpaceNode: public FunctionSpace{
+class FunctionSpaceNode: public FunctionSpaceScalar{
+ public:
   FunctionSpaceNode(const GroupOfElement& goe, int order);
-
+  
   virtual ~FunctionSpaceNode(void);
-
-  std::vector<double> interpolateAtNodes(const MElement& element, 
-					 const std::vector<double>& coef) const;
+  
+  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