From d2d4cb79a8ca72cd0ab73d78d6f3f11ee474c73d Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Thu, 16 Jan 2014 16:38:20 +0000
Subject: [PATCH] Interpolation of Derivative in FunctionSpace
---
FunctionSpace/Basis.h | 4 +++
FunctionSpace/BasisHierarchical0Form.cpp | 18 ++++++++++++
FunctionSpace/BasisHierarchical0Form.h | 4 +++
FunctionSpace/BasisHierarchical1Form.cpp | 6 ++++
FunctionSpace/BasisHierarchical1Form.h | 4 +++
FunctionSpace/BasisLagrange.cpp | 6 ++++
FunctionSpace/BasisLagrange.h | 4 +++
FunctionSpace/FunctionSpaceScalar.cpp | 36 ++++++++++++++++++++++++
FunctionSpace/FunctionSpaceScalar.h | 23 +++++++++++++++
9 files changed, 105 insertions(+)
diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 828044e0b6..49de2e09c2 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 d4137e18fc..351001b9e2 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 ecdc6627a7..c1feaba252 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 8260f3a9b2..783bc19f30 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 48e32d892d..7cd9bc42fe 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 199615fae9..d5e9c2ff9c 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 b3a3bd168b..e3b39b2820 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 c4e96462cf..f04f030ebf 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 33bb161870..95dd1a4b91 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
--
GitLab