Bug in ReferenceSpace::getJacobian -- Has to take sign of determinant into...

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

@@ -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);
 

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::

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

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;
+}

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

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){