From 06421cbf650d03ca52debf3ce3774a1461576990 Mon Sep 17 00:00:00 2001
From: Jeremy L Thompson <jeremy@jeremylt.org>
Date: Mon, 13 Dec 2021 14:42:41 -0700
Subject: [PATCH] basis - fix scaling for mesh

---
 src/Basis/Base.jl    | 31 ++++++++-----------------------
 src/Operator/Base.jl |  5 +----
 2 files changed, 9 insertions(+), 27 deletions(-)

diff --git a/src/Basis/Base.jl b/src/Basis/Base.jl
index 26caef623..77203135b 100644
--- a/src/Basis/Base.jl
+++ b/src/Basis/Base.jl
@@ -1230,31 +1230,25 @@ function getdXdxgradient(basis::TensorBasis, mesh::Mesh)
 
     # coordinate transformation
     gradient = basis.gradient
-    dxdX = (basis.gradient1d*basis.nodes1d)[1]
 
     # adjust for mesh
     if dimension == 1
         # 1D
-        return gradient / (dxdX * mesh.dx)
+        return gradient / mesh.dx
     elseif dimension == 2
         # 2D
-        scalex = 1 / (dxdX * mesh.dx)
-        scaley = 1 / (dxdX * mesh.dy)
         numberquadraturepoints = basis.numberquadraturepoints
         return [
-            gradient[1:numberquadraturepoints, :] * scalex
-            gradient[numberquadraturepoints+1:end, :] * scaley
+            gradient[1:numberquadraturepoints, :] / mesh.dx
+            gradient[numberquadraturepoints+1:end, :] / mesh.dy
         ]
     elseif dimension == 3
         # 3D
-        scalex = 1 / (dxdX * mesh.dx)
-        scaley = 1 / (dxdX * mesh.dy)
-        scalez = 1 / (dxdX * mesh.dz)
         numberquadraturepoints = basis.numberquadraturepoints
         return [
-            gradient[1:numberquadraturepoints, :] * scalex
-            gradient[numberquadraturepoints+1:2*numberquadraturepoints, :] * scaley
-            gradient[2*numberquadraturepoints+1:end, :] * scalez
+            gradient[1:numberquadraturepoints, :] / mesh.dx
+            gradient[numberquadraturepoints+1:2*numberquadraturepoints, :] / mesh.dy
+            gradient[2*numberquadraturepoints+1:end, :] / mesh.dz
         ]
     end
 end
@@ -1288,7 +1282,7 @@ basis = TensorH1LagrangeBasis(4, 3, 1, 1);
 weights = LFAToolkit.getdxdXquadratureweights(basis, mesh);
 
 # verify
-@assert basis.quadratureweights / 2 ≈ weights
+@assert basis.quadratureweights * mesh.volume / basis.volume ≈ weights
 
 # output
 
@@ -1303,17 +1297,8 @@ function getdxdXquadratureweights(basis::TensorBasis, mesh::Mesh)
         error("mesh dimension must match basis dimension") # COV_EXCL_LINE
     end
 
-    # coordinate transformation
-    dxdX = (basis.gradient1d*basis.nodes1d)[1]
-
     # adjust for mesh
-    if dimension == 1
-        return dxdX * basis.quadratureweights / mesh.dx
-    elseif dimension == 2
-        return dxdX^2 * basis.quadratureweights / (mesh.dx * mesh.dy)
-    elseif dimension == 3
-        return dxdX^3 * basis.quadratureweights / (mesh.dx * mesh.dy * mesh.dz)
-    end
+    return basis.quadratureweights * mesh.volume / basis.volume
 end
 
 function getdxdXquadratureweights(basis::NonTensorBasis, mesh::Mesh)
diff --git a/src/Operator/Base.jl b/src/Operator/Base.jl
index dfc690e5e..7bb7680f1 100644
--- a/src/Operator/Base.jl
+++ b/src/Operator/Base.jl
@@ -357,14 +357,11 @@ function getelementmatrix(operator::Operator)
         end
 
         # quadrature weight input index
-        weightscale = 1.0
         if weightinputindex != 0
             quadratureweights = getdxdXquadratureweights(
                 operator.inputs[weightinputindex].basis,
                 operator.mesh,
             )
-            weightscale =
-                operator.mesh.volume / operator.inputs[weightinputindex].basis.volume
         end
 
         # outputs
@@ -423,7 +420,7 @@ function getelementmatrix(operator::Operator)
             for q = 1:numberquadraturepoints
                 # set quadrature weight
                 if weightinputindex != 0
-                    weakforminputs[weightinputindex][1] = quadratureweights[q] * weightscale
+                    weakforminputs[weightinputindex][1] = quadratureweights[q]
                 end
 
                 # fill sparse matrix