Interpolation of Derivative in FunctionSpace

FunctionSpace/Basis.h

@@ -93,6 +93,10 @@ class Basis{
                             size_t orientation,
                             double u, double v, double w) const = 0;
 
+  virtual void getDerivative(fullMatrix<double>& retValues,
+                             const MElement& element,
+                             double u, double v, double w) const = 0;
+
   // Precompute Functions //
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const = 0;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const = 0;

FunctionSpace/BasisHierarchical0Form.cpp

@@ -72,6 +72,24 @@ getFunctions(fullMatrix<double>& retValues,
     retValues(i, 0) = basis[orientation][i]->at(u, v, w);
 }
 
+void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
+                                           const MElement& element,
+                                           double u, double v, double w) const{
+  // Get Grad //
+  if(!hasGrad)
+    getGrad();
+
+  // Define Orientation //
+  const size_t orientation = refSpace->getReferenceSpace(element);
+
+  // Fill Matrix //
+  for(size_t i = 0; i < nFunction; i++){
+    retValues(i, 0) = grad[orientation][i]->at(0).at(u, v, w);
+    retValues(i, 1) = grad[orientation][i]->at(1).at(u, v, w);
+    retValues(i, 2) = grad[orientation][i]->at(2).at(u, v, w);
+  }
+}
+
 void BasisHierarchical0Form::
 preEvaluateFunctions(const fullMatrix<double>& point) const{
   // Delete if older //

FunctionSpace/BasisHierarchical0Form.h

@@ -38,6 +38,10 @@ class BasisHierarchical0Form: public BasisLocal{
                             size_t orientation,
                             double u, double v, double w) const;
 
+  virtual void getDerivative(fullMatrix<double>& retValues,
+                             const MElement& element,
+                             double u, double v, double w) const;
+
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
 

FunctionSpace/BasisHierarchical1Form.cpp

@@ -84,6 +84,12 @@ getFunctions(fullMatrix<double>& retValues,
   }
 }
 
+void BasisHierarchical1Form::getDerivative(fullMatrix<double>& retValues,
+                                           const MElement& element,
+                                           double u, double v, double w) const{
+  throw Exception("Not Implemented");
+}
+
 void BasisHierarchical1Form::
 preEvaluateFunctions(const fullMatrix<double>& point) const{
   // Delete if older //

FunctionSpace/BasisHierarchical1Form.h

@@ -38,6 +38,10 @@ class BasisHierarchical1Form: public BasisLocal{
                             size_t orientation,
                             double u, double v, double w) const;
 
+  virtual void getDerivative(fullMatrix<double>& retValues,
+                             const MElement& element,
+                             double u, double v, double w) const;
+
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
 

FunctionSpace/BasisLagrange.cpp

@@ -62,6 +62,12 @@ getFunctions(fullMatrix<double>& retValues,
   permutation(orientation, retValues);
 }
 
+void BasisLagrange::getDerivative(fullMatrix<double>& retValues,
+                                  const MElement& element,
+                                  double u, double v, double w) const{
+  throw Exception("Not Implemented");
+}
+
 void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
   // Delete if older //
   if(preEvaluated)

FunctionSpace/BasisLagrange.h

@@ -43,6 +43,10 @@ class BasisLagrange: public BasisLocal{
                             size_t orientation,
                             double u, double v, double w) const;
 
+  virtual void getDerivative(fullMatrix<double>& retValues,
+                             const MElement& element,
+                             double u, double v, double w) const;
+
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
 

FunctionSpace/FunctionSpaceScalar.cpp

+#include "Mapper.h"
 #include "FunctionSpaceScalar.h"
 
 FunctionSpaceScalar::FunctionSpaceScalar(GroupOfElement& goe,
@@ -30,3 +31,38 @@ interpolateInABC(const MElement& element,
   // Return Interpolated Value //
   return val;
 }
+
+fullVector<double> FunctionSpaceScalar::
+interpolateDerivativeInABC(const MElement& element,
+                           const std::vector<double>& coef,
+                           double abc[3]) const{
+  // Get Jacobian //
+  fullMatrix<double> invJac(3, 3);
+  (*basis)[0]->getReferenceSpace().getJacobian(element,
+                                               abc[0], abc[1], abc[2],
+                                               invJac);
+  invJac.invertInPlace();
+
+  // Get Basis Functions //
+  const size_t nFun = (*basis)[0]->getNFunction();
+  fullMatrix<double> fun(nFun, 3);
+
+  (*basis)[0]->getDerivative(fun, element, abc[0], abc[1], abc[2]);
+
+  // Interpolate (in Reference Place) //
+  fullMatrix<double> val(1, 3);
+  val(0, 0) = 0;
+  val(0, 1) = 0;
+  val(0, 2) = 0;
+
+  for(size_t i = 0; i < nFun; i++){
+    val(0, 0) += fun(i, 0) * coef[i];
+    val(0, 1) += fun(i, 1) * coef[i];
+    val(0, 2) += fun(i, 2) * coef[i];
+  }
+
+  // Return Interpolated Value //
+  fullVector<double> map(3);
+  Mapper::hCurl(val, 0, 0, invJac, map);
+  return map;
+}

FunctionSpace/FunctionSpaceScalar.h

@@ -31,10 +31,19 @@ class FunctionSpaceScalar : public FunctionSpace{
                           const std::vector<double>& coef,
                           const fullVector<double>& uvw) const;
 
+  fullVector<double>
+    interpolateDerivative(const MElement& element,
+                          const std::vector<double>& coef,
+                          const fullVector<double>& xyz) const;
+
  private:
   double interpolateInABC(const MElement& element,
                           const std::vector<double>& coef,
                           double abc[3]) const;
+  fullVector<double>
+    interpolateDerivativeInABC(const MElement& element,
+                               const std::vector<double>& coef,
+                               double abc[3]) const;
 };
 
 
@@ -122,4 +131,18 @@ interpolateInRefSpace(const MElement& element,
   return interpolateInABC(element, coef, abc);
 }
 
+inline fullVector<double> FunctionSpaceScalar::
+interpolateDerivative(const MElement& element,
+                      const std::vector<double>& coef,
+                      const fullVector<double>& xyz) const{
+
+  // Get ABC Space coordinate //
+  double abc[3];
+  (*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
+                                                   xyz(0), xyz(1), xyz(2),
+                                                   abc);
+  // Interpolate in ABC //
+  return interpolateDerivativeInABC(element, coef, abc);
+}
+
 #endif