diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 828044e0b6061dd8bd61d32e3ab261825bdb9a28..49de2e09c269bfbf8a61f56050ad5005b3452c7b 100644
--- a/FunctionSpace/Basis.h
+++ b/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;
diff --git a/FunctionSpace/BasisHierarchical0Form.cpp b/FunctionSpace/BasisHierarchical0Form.cpp
index d4137e18fcd2d5ecc0b407f0575cbd68d363c0aa..351001b9e25e6192ba9869077605074d7fee3251 100644
--- a/FunctionSpace/BasisHierarchical0Form.cpp
+++ b/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 //
diff --git a/FunctionSpace/BasisHierarchical0Form.h b/FunctionSpace/BasisHierarchical0Form.h
index ecdc6627a7cb1490dc1ec6cde8d8faec7431c71a..c1feaba2523e99211760725a942ab511e638c79e 100644
--- a/FunctionSpace/BasisHierarchical0Form.h
+++ b/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;
diff --git a/FunctionSpace/BasisHierarchical1Form.cpp b/FunctionSpace/BasisHierarchical1Form.cpp
index 8260f3a9b263bc4265f2b58ca9d892cd612b5cee..783bc19f30cf34ad71600084ecce4b99d5a0be3f 100644
--- a/FunctionSpace/BasisHierarchical1Form.cpp
+++ b/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 //
diff --git a/FunctionSpace/BasisHierarchical1Form.h b/FunctionSpace/BasisHierarchical1Form.h
index 48e32d892d1f381037532e87bb479145e075d828..7cd9bc42fe1f1d194677925fa1cadf82c43238a5 100644
--- a/FunctionSpace/BasisHierarchical1Form.h
+++ b/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;
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index 199615fae979362078c48c2cc54de97763527f08..d5e9c2ff9c4981565f6183ea6ded6f66a394cc8d 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/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)
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index b3a3bd168bc58750734618d2b90e691e48f103c1..e3b39b28202bca6b21afc3c8a64ec09708da8496 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/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;
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index c4e96462cfce57a45c564b89db381c1978787966..f04f030ebf1284a48471819fcf6b4c7ad6edbf9a 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -1,3 +1,4 @@
+#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;
+}
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 33bb161870d91120cef2561ee7a52b276d9fa80d..95dd1a4b91decdc56fc7a1752f0d64b5aa12822b 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/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