From 3e38e65de64d79a185bbadc709ec376c6b13d594 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Mon, 7 Apr 2014 13:30:35 +0000
Subject: [PATCH] Bug in ReferenceSpace::getJacobian -- Has to take sign of
determinant into account when flipping elements ! +
FunctionSpaceVector::interpolateDerivative + Mapper::hDiv -- do not use fabs
in mapping
---
FunctionSpace/BasisHierarchical0Form.cpp | 11 +++++
FunctionSpace/BasisHierarchical1Form.cpp | 25 +++++++++-
FunctionSpace/FunctionSpaceScalar.h | 1 +
FunctionSpace/FunctionSpaceVector.cpp | 34 ++++++++++++++
FunctionSpace/FunctionSpaceVector.h | 60 ++++++++++++++++++++++++
FunctionSpace/ReferenceSpace.cpp | 22 ++++++---
6 files changed, 145 insertions(+), 8 deletions(-)
diff --git a/FunctionSpace/BasisHierarchical0Form.cpp b/FunctionSpace/BasisHierarchical0Form.cpp
index df97eebade..0405beb5fe 100644
--- a/FunctionSpace/BasisHierarchical0Form.cpp
+++ b/FunctionSpace/BasisHierarchical0Form.cpp
@@ -82,10 +82,21 @@ getFunctions(fullMatrix<double>& retValues,
void BasisHierarchical0Form::getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
+
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
// Get Grad //
if(!hasGrad)
getGrad();
+ // Delete if older //
+ if(preEvaluatedGrad){
+ for(size_t i = 0; i < nOrientation; i++)
+ delete preEvaluatedGradFunction[i];
+
+ delete[] preEvaluatedGradFunction;
+ }
+
// Define Orientation //
const size_t orientation = ReferenceSpaceManager::getOrientation(element);
diff --git a/FunctionSpace/BasisHierarchical1Form.cpp b/FunctionSpace/BasisHierarchical1Form.cpp
index 06b51302b2..bd53fa7904 100644
--- a/FunctionSpace/BasisHierarchical1Form.cpp
+++ b/FunctionSpace/BasisHierarchical1Form.cpp
@@ -92,7 +92,30 @@ getFunctions(fullMatrix<double>& retValues,
void BasisHierarchical1Form::getDerivative(fullMatrix<double>& retValues,
const MElement& element,
double u, double v, double w) const{
- throw Exception("Not Implemented");
+
+ const size_t nOrientation = ReferenceSpaceManager::getNOrientation(getType());
+
+ // Build Curl //
+ if(!hasCurl)
+ getCurl();
+
+ // Delete if older //
+ if(preEvaluatedCurl){
+ for(size_t i = 0; i < nOrientation; i++)
+ delete preEvaluatedCurlFunction[i];
+
+ delete[] preEvaluatedCurlFunction;
+ }
+
+ // Define Orientation //
+ const size_t orientation = ReferenceSpaceManager::getOrientation(element);
+
+ // Fill Matrix //
+ for(size_t i = 0; i < nFunction; i++){
+ retValues(i, 0) = curl[orientation][i]->at(0).at(u, v, w);
+ retValues(i, 1) = curl[orientation][i]->at(1).at(u, v, w);
+ retValues(i, 2) = curl[orientation][i]->at(2).at(u, v, w);
+ }
}
void BasisHierarchical1Form::
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 42428d8f0c..0cae40d3ba 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -98,6 +98,7 @@ class FunctionSpaceScalar : public FunctionSpace{
If the given coordinate are not in the given
element @em Bad @em Things may happend
+ **
@fn FunctionSpaceScalar::interpolateDerivative
@param element The MElement to interpolate on
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index 1fb78e7d39..3303886db1 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -57,3 +57,37 @@ interpolateInABC(const MElement& element,
Mapper::hCurl(val, 0, 0, invJac, map);
return map;
}
+
+fullVector<double> FunctionSpaceVector::
+interpolateDerivativeInABC(const MElement& element,
+ const std::vector<double>& coef,
+ double abc[3]) const{
+ // Get Jacobian //
+ fullMatrix<double> jac(3, 3);
+ double det =
+ ReferenceSpaceManager::getJacobian(element, abc[0], abc[1], abc[2], jac);
+
+ // Get Basis Functions //
+ const Basis& basis = getBasis(element);
+ const size_t nFun = basis.getNFunction();
+ fullMatrix<double> fun(nFun, 3);
+
+ basis.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::hDiv(val, 0, 0, jac, det, map);
+ return map;
+}
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 7bdaff6a73..5f2a027760 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -33,11 +33,25 @@ class FunctionSpaceVector : 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;
+
+ fullVector<double>
+ interpolateDerivativeInRefSpace(const MElement& element,
+ const std::vector<double>& coef,
+ const fullVector<double>& uvw) const;
private:
fullVector<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;
};
@@ -91,6 +105,26 @@ class FunctionSpaceVector : public FunctionSpace{
If the given coordinate are not in the given
element @em Bad @em Things may happend
+ **
+
+ @fn FunctionSpaceVector::interpolateDerivative
+ @param element The MElement to interpolate on
+ @param coef The coefficients of the interpolation
+ @param xyz The coordinate (of a point inside the given element)
+ of the interpolation in the @em physical space
+
+ Same as FunctionSpaceVector::interpolate(element, coef, xyz),
+ but this method iterpolates the derivative.
+ **
+
+ @fn FunctionSpaceVector::interpolateDerivativeInRefSpace
+ @param element The MElement to interpolate on
+ @param coef The coefficients of the interpolation
+ @param uvw The coordinate (of a point inside the given element)
+ of the interpolation in the @em reference space
+
+ Same as FunctionSpaceVector::interpolateInRefSpace(element, coef, uvw),
+ but this method iterpolates the derivative.
*/
@@ -124,4 +158,30 @@ interpolateInRefSpace(const MElement& element,
return interpolateInABC(element, coef, abc);
}
+inline fullVector<double> FunctionSpaceVector::
+interpolateDerivative(const MElement& element,
+ const std::vector<double>& coef,
+ const fullVector<double>& xyz) const{
+
+ // Get ABC Space coordinate //
+ double abc[3];
+ ReferenceSpaceManager::mapFromXYZtoABC(element, xyz(0), xyz(1), xyz(2), abc);
+
+ // Interpolate in ABC //
+ return interpolateDerivativeInABC(element, coef, abc);
+}
+
+inline fullVector<double> FunctionSpaceVector::
+interpolateDerivativeInRefSpace(const MElement& element,
+ const std::vector<double>& coef,
+ const fullVector<double>& uvw) const{
+
+ // Get ABC Space coordinate //
+ double abc[3];
+ ReferenceSpaceManager::mapFromUVWtoABC(element, uvw(0), uvw(1), uvw(2), abc);
+
+ // Interpolate in ABC //
+ return interpolateDerivativeInABC(element, coef, abc);
+}
+
#endif
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index 63f1912097..bdd10dc17e 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -244,8 +244,8 @@ isCyclicPermutation(const vector<size_t>& pTest,
// Test if we have the same connectivity
bool isSameConnectivity = isSameEdge(pTest, pRef);
- if(isCyclic && isSameConnectivity){
- //if(isSameConnectivity){
+ //if(isCyclic && isSameConnectivity){
+ if(isSameConnectivity){
tri.first = true;
tri.second = getRefIndexPermutation(pRef, pTest);
tri.third = getReverseIndexPermutation(pRef, pTest);
@@ -779,17 +779,25 @@ double ReferenceSpace::getJacobian(const MElement& element,
double uvw[3];
mapFromABCtoUVW(element, a, b, c, uvw);
- // UVW to XYZ Jacobian + Determinant //
- // NB: Volume is not modified when we go from ABC to UVW //
- // --> Determinant unchanged //
+ // UVW to XYZ Jacobian //
fullMatrix<double> jacUVWtoXYZ(3, 3);
- double det = element.getJacobian(uvw[0], uvw[1], uvw[2], jacUVWtoXYZ);
+ element.getJacobian(uvw[0], uvw[1], uvw[2], jacUVWtoXYZ);
// Product of the two Jacobians & Return //
// Do a naive gemm, so that we do not encounter nested threads
// (limitation of OpenBLAS)
jac.gemm_naive(jacABCtoUVW, jacUVWtoXYZ, 1, 0);
- return det;
+
+ // New Jacobian determinant (same as jacUVWtoXYZ with maybe a SIGN CHANGE) //
+ // ---> Has to be recomputed!
+ return
+ jac(0, 0) * jac(1, 1) * jac(2, 2) +
+ jac(0, 1) * jac(1, 2) * jac(2, 0) +
+ jac(0, 2) * jac(1, 0) * jac(2, 1) -
+
+ jac(0, 2) * jac(1, 1) * jac(2, 0) -
+ jac(0, 1) * jac(1, 0) * jac(2, 2) -
+ jac(0, 0) * jac(1, 2) * jac(2, 1);
}
void ReferenceSpace::regularize(size_t dim, fullMatrix<double>& jac){
--
GitLab