diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index e5c8c2e276dde519c706e274a85f0d3edbb6980a..5a1795328a0a1daf458aa6138b8fd76489c8337d 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -24,7 +24,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-from six.moves import range
+from six.moves import range, intern
import numpy as np
import pymbolic.primitives
@@ -73,6 +73,9 @@ Geometry data
.. autofunction:: mv_nodes
.. autofunction:: forward_metric_derivative
.. autofunction:: inverse_metric_derivative
+.. autofunction:: forward_metric_derivative_mat
+.. autofunction:: inverse_metric_derivative_mat
+.. autofunction:: pseudoscalar
.. autofunction:: area_element
.. autofunction:: mv_normal
.. autofunction:: normal
@@ -286,30 +289,96 @@ def forward_metric_derivative(xyz_axis, rst_axis, where=None,
scope=cse_scope.DISCRETIZATION)
+def forward_metric_derivative_vector(ambient_dim, rst_axis, where=None,
+ quadrature_tag=None):
+ return make_obj_array([
+ forward_metric_derivative(
+ i, rst_axis, where=where, quadrature_tag=quadrature_tag)
+ for i in range(ambient_dim)])
+
+
+def forward_metric_derivative_mv(ambient_dim, rst_axis, where=None,
+ quadrature_tag=None):
+ return MultiVector(
+ forward_metric_derivative_vector(
+ ambient_dim, rst_axis, where=where, quadrature_tag=quadrature_tag))
+
+
def parametrization_derivative(ambient_dim, dim=None, where=None,
quadrature_tag=None):
if dim is None:
dim = ambient_dim
- par_grad = np.zeros((ambient_dim, dim), np.object)
- for i in range(ambient_dim):
- for j in range(dim):
- par_grad[i, j] = forward_metric_derivative(
- i, j, where=where, quadrature_tag=quadrature_tag)
-
from pytools import product
- return product(MultiVector(vec) for vec in par_grad.T)
+ return product(
+ forward_metric_derivative_mv(
+ ambient_dim, rst_axis, where, quadrature_tag)
+ for rst_axis in range(dim))
def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
- quadrature_tag=None):
+ where=None, quadrature_tag=None):
if dim is None:
dim = ambient_dim
if dim != ambient_dim:
raise ValueError("not clear what inverse_metric_derivative means if "
"the derivative matrix is not square")
- raise NotImplementedError()
+
+ par_vecs = [
+ forward_metric_derivative_mv(
+ ambient_dim, rst, where, quadrature_tag)
+ for rst in range(dim)]
+
+ # Yay Cramer's rule! (o rly?)
+ from functools import reduce, partial
+ from operator import xor as outerprod_op
+ outerprod = partial(reduce, outerprod_op)
+
+ def outprod_with_unit(i, at):
+ unit_vec = np.zeros(dim)
+ unit_vec[i] = 1
+
+ vecs = par_vecs[:]
+ vecs[at] = MultiVector(unit_vec)
+
+ return outerprod(vecs)
+
+ volume_pseudoscalar_inv = cse(outerprod(
+ forward_metric_derivative_mv(
+ ambient_dim, rst_axis, where, quadrature_tag)
+ for rst_axis in range(dim)).inv())
+
+ return (outprod_with_unit(xyz_axis, rst_axis)
+ * volume_pseudoscalar_inv
+ ).as_scalar()
+
+
+def forward_metric_derivative_mat(ambient_dim, dim=None,
+ where=None, quadrature_tag=None):
+ if dim is None:
+ dim = ambient_dim
+
+ result = np.zeros((ambient_dim, dim), dtype=np.object)
+ for j in range(dim):
+ result[:, j] = forward_metric_derivative_vector(ambient_dim, j,
+ where=where, quadrature_tag=quadrature_tag)
+ return result
+
+
+def inverse_metric_derivative_mat(ambient_dim, dim=None,
+ where=None, quadrature_tag=None):
+ if dim is None:
+ dim = ambient_dim
+
+ result = np.zeros((ambient_dim, dim), dtype=np.object)
+ for i in range(dim):
+ for j in range(ambient_dim):
+ result[i, j] = inverse_metric_derivative(
+ i, j, ambient_dim, dim,
+ where=where, quadrature_tag=quadrature_tag)
+
+ return result
def pseudoscalar(ambient_dim, dim=None, where=None, quadrature_tag=None):
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 060766a89a2c18d3fa332c0adcac4ba631d1dc83..83f7c78f3c6df6aff416a2c9fa157ad5de6a9446 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -37,9 +37,52 @@ from pyopencl.tools import ( # noqa
import logging
logger = logging.getLogger(__name__)
+from grudge import sym, bind, Discretization
-def test_nothing(ctx_getter):
- pass
+
+@pytest.mark.parametrize("dims", [2, 3])
+def test_inverse_metric(ctx_getter, dims):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(a=(-0.5,)*dims, b=(0.5,)*dims,
+ n=(6,)*dims, order=4)
+
+ def m(x):
+ result = np.empty_like(x)
+ result[0] = (
+ 1.5*x[0] + np.cos(x[0])
+ + 0.1*np.sin(10*x[1]))
+ result[1] = (
+ 0.05*np.cos(10*x[0])
+ + 1.3*x[1] + np.sin(x[1]))
+ if len(x) == 3:
+ result[2] = x[2]
+ return result
+
+ from meshmode.mesh.processing import map_mesh
+ mesh = map_mesh(mesh, m)
+
+ discr = Discretization(cl_ctx, mesh, order=4)
+
+ sym_op = (
+ sym.forward_metric_derivative_mat(mesh.dim)
+ .dot(
+ sym.inverse_metric_derivative_mat(mesh.dim)
+ )
+ .reshape(-1))
+
+ op = bind(discr, sym_op)
+ mat = op(queue).reshape(mesh.dim, mesh.dim)
+
+ for i in range(mesh.dim):
+ for j in range(mesh.dim):
+ tgt = 1 if i == j else 0
+
+ err = np.max(np.abs((mat[i, j] - tgt).get(queue=queue)))
+ print(i, j, err)
+ assert err < 1e-12, (i, j, err)
# You can test individual routines by typing