diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index cb784333e1bc6611c31f209c448beb8e9b4123ec..abdf1d57dcf42733387a061d244c9a816176b1ab 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -402,14 +402,13 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ, bool idea
 
 // Calculate (signed, possibly scaled) Jacobian for one element, with normal vectors to straight element
 // for regularization. Evaluation points depend on the given matrices for shape function gradients.
-template<bool scaling>
-inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
-                                              const fullMatrix<double> &gSMatX,
-                                              const fullMatrix<double> &gSMatY,
-                                              const fullMatrix<double> &gSMatZ,
-                                              const fullMatrix<double> &nodesXYZ,
-                                              bool idealNorm,
-                                              fullVector<double> &jacobian) const
+void JacobianBasis::getJacobianGeneral(int nJacNodes,
+                                       const fullMatrix<double> &gSMatX,
+                                       const fullMatrix<double> &gSMatY,
+                                       const fullMatrix<double> &gSMatZ,
+                                       const fullMatrix<double> &nodesXYZ,
+                                       bool idealNorm, bool scaling,
+                                       fullVector<double> &jacobian) const
 {
   switch (_dim) {
 
@@ -493,45 +492,18 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
 
 }
 
-// Calculate signed Jacobian for one element, with normal vectors to straight element for
-// regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getSignedJacobianGeneral(int nJacNodes,
-                                             const fullMatrix<double> &gSMatX,
-                                             const fullMatrix<double> &gSMatY,
-                                             const fullMatrix<double> &gSMatZ,
-                                             const fullMatrix<double> &nodesXYZ,
-                                             bool idealNorm,
-                                             fullVector<double> &jacobian) const
-{
-  getJacobianGeneral<false>(nJacNodes, gSMatX,  gSMatY, gSMatZ, nodesXYZ, idealNorm, jacobian);
-}
-
-// Calculate (signed) scaled Jacobian for one element, with normal vectors to straight element
-// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getScaledJacobianGeneral(int nJacNodes,
-                                             const fullMatrix<double> &gSMatX,
-                                             const fullMatrix<double> &gSMatY,
-                                             const fullMatrix<double> &gSMatZ,
-                                             const fullMatrix<double> &nodesXYZ,
-                                             bool idealNorm,
-                                             fullVector<double> &jacobian) const
-{
-  getJacobianGeneral<true>(nJacNodes, gSMatX,  gSMatY, gSMatZ, nodesXYZ, idealNorm, jacobian);
-}
-
 // Calculate (signed, possibly scaled) Jacobian for several elements, with normal vectors to straight
 // elements for regularization. Evaluation points depend on the given matrices for shape function gradients.
 // TODO: Optimize and test 1D & 2D
-template<bool scaling>
-inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
-                                              const fullMatrix<double> &gSMatX,
-                                              const fullMatrix<double> &gSMatY,
-                                              const fullMatrix<double> &gSMatZ,
-                                              const fullMatrix<double> &nodesX,
-                                              const fullMatrix<double> &nodesY,
-                                              const fullMatrix<double> &nodesZ,
-                                              bool idealNorm,
-                                              fullMatrix<double> &jacobian) const
+void JacobianBasis::getJacobianGeneral(int nJacNodes,
+                                       const fullMatrix<double> &gSMatX,
+                                       const fullMatrix<double> &gSMatY,
+                                       const fullMatrix<double> &gSMatZ,
+                                       const fullMatrix<double> &nodesX,
+                                       const fullMatrix<double> &nodesY,
+                                       const fullMatrix<double> &nodesZ,
+                                       bool idealNorm, bool scaling,
+                                       fullMatrix<double> &jacobian) const
 {
   switch (_dim) {
 
@@ -636,38 +608,6 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
 
 }
 
-// Calculate signed Jacobian for several elements, with normal vectors to straight elements for
-// regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getSignedJacobianGeneral(int nJacNodes,
-                                             const fullMatrix<double> &gSMatX,
-                                             const fullMatrix<double> &gSMatY,
-                                             const fullMatrix<double> &gSMatZ,
-                                             const fullMatrix<double> &nodesX,
-                                             const fullMatrix<double> &nodesY,
-                                             const fullMatrix<double> &nodesZ,
-                                             bool idealNorm,
-                                             fullMatrix<double> &jacobian) const
-{
-  getJacobianGeneral<false>(nJacNodes, gSMatX,  gSMatY, gSMatZ,
-                            nodesX, nodesY, nodesZ, idealNorm, jacobian);
-}
-
-// Calculate (signed) scaled Jacobian for several elements, with normal vectors to straight elements
-// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getScaledJacobianGeneral(int nJacNodes,
-                                             const fullMatrix<double> &gSMatX,
-                                             const fullMatrix<double> &gSMatY,
-                                             const fullMatrix<double> &gSMatZ,
-                                             const fullMatrix<double> &nodesX,
-                                             const fullMatrix<double> &nodesY,
-                                             const fullMatrix<double> &nodesZ,
-                                             bool idealNorm,
-                                             fullMatrix<double> &jacobian) const
-{
-  getJacobianGeneral<true>(nJacNodes, gSMatX,  gSMatY, gSMatZ,
-                           nodesX, nodesY, nodesZ, idealNorm, jacobian);
-}
-
 // Calculate (signed) Jacobian and its gradients for one element, with normal vectors to straight element
 // for regularization. Evaluation points depend on the given matrices for shape function gradients.
 inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index bbcca1b05957632b700c70599ae8237c8c05aaf0..d0251d7c8a9bae2cc1caff8b1cdd13b3ad497b58 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -112,79 +112,71 @@ public:
                                     _gradBasis->gradShapeIdealMatZ,
                                     nodesXYZ, normals, JDJ);
   }
-//  inline void getSignedIdealJacAndGradientsFast(const fullMatrix<double> &nodesXYZ,
-//                                                const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
-//    getSignedIdealJacAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
-//                                         gradShapeMatZFast,nodesXYZ,normals,JDJ);
-//  }
   void getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
                                 const fullMatrix<double> &nodesXYZStraight,
                                 fullVector<double> &lambdaJ,
                                 fullMatrix<double> &gradLambdaJ) const;
   inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
                                 fullVector<double> &jacobian) const {
-    getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
-                             _gradBasis->gradShapeMatY,
-                             _gradBasis->gradShapeMatZ,
-                             nodesXYZ, false, jacobian);
+    getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+                       _gradBasis->gradShapeMatY,
+                       _gradBasis->gradShapeMatZ,
+                       nodesXYZ, false, false, jacobian);
   }
   inline void getSignedJacobian(const fullMatrix<double> &nodesX,
                                 const fullMatrix<double> &nodesY,
                                 const fullMatrix<double> &nodesZ,
                                 fullMatrix<double> &jacobian) const {
-    getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
-                             _gradBasis->gradShapeMatY,
-                             _gradBasis->gradShapeMatZ,
-                             nodesX, nodesY, nodesZ, false, jacobian);
+    getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+                       _gradBasis->gradShapeMatY,
+                       _gradBasis->gradShapeMatZ,
+                       nodesX, nodesY, nodesZ, false, false, jacobian);
   }
   inline void getSignedIdealJacobian(const fullMatrix<double> &nodesXYZ,
                                      fullVector<double> &jacobian) const {
-    getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
-                             _gradBasis->gradShapeIdealMatY,
-                             _gradBasis->gradShapeIdealMatZ,
-                             nodesXYZ, true, jacobian);
+    getJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
+                       _gradBasis->gradShapeIdealMatY,
+                       _gradBasis->gradShapeIdealMatZ,
+                       nodesXYZ, true, false, jacobian);
   }
   inline void getSignedIdealJacobian(const fullMatrix<double> &nodesX,
                                      const fullMatrix<double> &nodesY,
                                      const fullMatrix<double> &nodesZ,
                                      fullMatrix<double> &jacobian) const {
-    getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
-                             _gradBasis->gradShapeIdealMatY,
-                             _gradBasis->gradShapeIdealMatZ,
-                             nodesX, nodesY, nodesZ, true, jacobian);
+    getJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
+                       _gradBasis->gradShapeIdealMatY,
+                       _gradBasis->gradShapeIdealMatZ,
+                       nodesX, nodesY, nodesZ, true, false, jacobian);
   }
   inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ,
                                 fullVector<double> &jacobian) const {
-    getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
-                             _gradBasis->gradShapeMatY,
-                             _gradBasis->gradShapeMatZ,
-                             nodesXYZ, false, jacobian);
+    getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+                       _gradBasis->gradShapeMatY,
+                       _gradBasis->gradShapeMatZ,
+                       nodesXYZ, false, true, jacobian);
   }
   inline void getScaledJacobian(const fullMatrix<double> &nodesX,
                                 const fullMatrix<double> &nodesY,
                                 const fullMatrix<double> &nodesZ,
                                 fullMatrix<double> &jacobian) const {
-    getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
-                             _gradBasis->gradShapeMatY,
-                             _gradBasis->gradShapeMatZ,
-                             nodesX, nodesY, nodesZ, false, jacobian);
+    getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+                       _gradBasis->gradShapeMatY,
+                       _gradBasis->gradShapeMatZ,
+                       nodesX, nodesY, nodesZ, false, true, jacobian);
   }
   inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ,
                                     fullVector<double> &jacobian) const {
-    getSignedJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
-                             gradShapeMatYFast, gradShapeMatZFast,
-                             nodesXYZ, false, jacobian);
+    getJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
+                       gradShapeMatYFast, gradShapeMatZFast,
+                       nodesXYZ, false, false, jacobian);
   }
-//  inline void getSignedIdealJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
-//    getSignedIdealJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
-//  }
   inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ,
                                     fullVector<double> &jacobian) const {
-    getScaledJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
-                             gradShapeMatYFast, gradShapeMatZFast,
-                             nodesXYZ, false, jacobian);
+    getJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
+                       gradShapeMatYFast, gradShapeMatZFast,
+                       nodesXYZ, false, true, jacobian);
   }
-  //
+
   inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
     getBezier()->matrixLag2Bez.mult(jac,bez);
   }
@@ -194,26 +186,24 @@ public:
   inline void primJac2Jac(const fullVector<double> &primJac, fullVector<double> &jac) const {
     matrixPrimJac2Jac.mult(primJac,jac);
   }
-  //
+
   void interpolate(const fullVector<double> &jacobian,
                    const fullMatrix<double> &uvw,
                    fullMatrix<double> &result, bool areBezier = false) const;
 
-  //
   static int jacobianOrder(int tag);
   static int jacobianOrder(int parentType, int order);
   static FuncSpaceData jacobianMatrixSpace(int type, int order);
 
 
  private :
-  template<bool scaling>
   void getJacobianGeneral(int nJacNodes,
                           const fullMatrix<double> &gSMatX,
                           const fullMatrix<double> &gSMatY,
                           const fullMatrix<double> &gSMatZ,
                           const fullMatrix<double> &nodesXYZ,
-                          bool idealNorm, fullVector<double> &jacobian) const;
-  template<bool scaling>
+                          bool idealNorm, bool scaling,
+                          fullVector<double> &jacobian) const;
   void getJacobianGeneral(int nJacNodes,
                           const fullMatrix<double> &gSMatX,
                           const fullMatrix<double> &gSMatY,
@@ -221,35 +211,8 @@ public:
                           const fullMatrix<double> &nodesX,
                           const fullMatrix<double> &nodesY,
                           const fullMatrix<double> &nodesZ,
-                          bool idealNorm, fullMatrix<double> &jacobian) const;
-  void getSignedJacobianGeneral(int nJacNodes,
-                                const fullMatrix<double> &gSMatX,
-                                const fullMatrix<double> &gSMatY,
-                                const fullMatrix<double> &gSMatZ,
-                                const fullMatrix<double> &nodesXYZ,
-                                bool idealNorm, fullVector<double> &jacobian) const;
-  void getSignedJacobianGeneral(int nJacNodes,
-                                const fullMatrix<double> &gSMatX,
-                                const fullMatrix<double> &gSMatY,
-                                const fullMatrix<double> &gSMatZ,
-                                const fullMatrix<double> &nodesX,
-                                const fullMatrix<double> &nodesY,
-                                const fullMatrix<double> &nodesZ,
-                                bool idealNorm, fullMatrix<double> &jacobian) const;
-  void getScaledJacobianGeneral(int nJacNodes,
-                                const fullMatrix<double> &gSMatX,
-                                const fullMatrix<double> &gSMatY,
-                                const fullMatrix<double> &gSMatZ,
-                                const fullMatrix<double> &nodesXYZ,
-                                bool idealNorm, fullVector<double> &jacobian) const;
-  void getScaledJacobianGeneral(int nJacNodes,
-                                const fullMatrix<double> &gSMatX,
-                                const fullMatrix<double> &gSMatY,
-                                const fullMatrix<double> &gSMatZ,
-                                const fullMatrix<double> &nodesX,
-                                const fullMatrix<double> &nodesY,
-                                const fullMatrix<double> &nodesZ,
-                                bool idealNorm, fullMatrix<double> &jacobian) const;
+                          bool idealNorm, bool scaling,
+                          fullMatrix<double> &jacobian) const;
 
   void getSignedJacAndGradientsGeneral(int nJacNodes,
                                        const fullMatrix<double> &gSMatX,