Address PR comments for high-level APIs

sebastiangrimberg · sebastiangrimberg · commit 0b742f503851 · 2023-04-07T11:50:04.000-07:00
Add Rust sibling functions and tests, add Python tests, fix Julia versioning for CI.
diff --git a/julia/LibCEED.jl/src/LibCEED.jl b/julia/LibCEED.jl/src/LibCEED.jl
@@ -150,7 +150,7 @@ include("Request.jl")
 include("Operator.jl")
 include("Misc.jl")
 
-const minimum_libceed_version = v"0.10.0"
+const minimum_libceed_version = v"0.12.0"
 
 function __init__()
     if !ceedversion_ge(minimum_libceed_version)
diff --git a/julia/LibCEED.jl/src/generated/libceed_bindings.jl b/julia/LibCEED.jl/src/generated/libceed_bindings.jl
@@ -1374,7 +1374,7 @@ const CEED_VERSION_MINOR = 11
 
 const CEED_VERSION_PATCH = 0
 
-const CEED_VERSION_RELEASE = true
+const CEED_VERSION_RELEASE = false
 
 # Skipping MacroDefinition: CEED_INTERN extern CEED_VISIBILITY ( hidden )
 
diff --git a/python/__init__.py b/python/__init__.py
@@ -7,7 +7,7 @@
 
 from .ceed import Ceed
 from .ceed_vector import Vector
-from .ceed_basis import Basis, BasisTensorH1, BasisTensorH1Lagrange, BasisH1
+from .ceed_basis import Basis, BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl
 from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction
 from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction
 from .ceed_operator import Operator, CompositeOperator
@@ -18,7 +18,7 @@
 # ------------------------------------------------------------------------------
 __all__ = ["Ceed",
            "Vector",
-           "Basis", "BasisTensorH1", "BasisTensorH1Lagrange", "BasisH1",
+           "Basis", "BasisTensorH1", "BasisTensorH1Lagrange", "BasisH1", "BasisHdiv", "BasisHcurl",
            "ElemRestriction", "StridedElemRestriction", "BlockedElemRestriction", "BlockedStridedelemRestriction",
            "QFunction", "QFunctionByName", "IdentityQFunction",
            "Operator", "CompositeOperator",
diff --git a/python/ceed.py b/python/ceed.py
@@ -13,7 +13,7 @@
 import tempfile
 from abc import ABC
 from .ceed_vector import Vector
-from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1
+from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl
 from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction
 from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction
 from .ceed_qfunctioncontext import QFunctionContext
@@ -356,7 +356,7 @@ def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight):
              *interp: Numpy array holding the row-major (nqpts * nnodes) matrix
                        expressing the values of nodal basis functions at
                        quadrature points
-             *grad: Numpy array holding the row-major (nqpts * dim * nnodes)
+             *grad: Numpy array holding the row-major (dim * nqpts * nnodes)
                      matrix expressing the derivatives of nodal basis functions
                      at quadrature points
              *qref: Numpy array of length (nqpts * dim) holding the locations of
@@ -370,6 +370,59 @@ def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight):
         return BasisH1(self, topo, ncomp, nnodes, nqpts,
                        interp, grad, qref, qweight)
 
+    def BasisHdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight):
+        """Ceed Hdiv Basis: finite element non tensor-product basis for H(div)
+             discretizations.
+
+           Args:
+             topo: topology of the element, e.g. hypercube, simplex, etc
+             ncomp: number of field components (1 for scalar fields)
+             nnodes: total number of nodes
+             nqpts: total number of quadrature points
+             *interp: Numpy array holding the row-major (dim * nqpts * nnodes)
+                       matrix expressing the values of basis functions at
+                       quadrature points
+             *div: Numpy array holding the row-major (nqpts * nnodes) matrix
+                    expressing the divergence of basis functions at
+                    quadrature points
+             *qref: Numpy array of length (nqpts * dim) holding the locations of
+                     quadrature points on the reference element [-1, 1]
+             *qweight: Numpy array of length nnodes holding the quadrature
+                        weights on the reference element
+
+           Returns:
+             basis: Ceed Basis"""
+
+        return BasisHdiv(self, topo, ncomp, nnodes, nqpts,
+                         interp, div, qref, qweight)
+
+    def BasisHcurl(self, topo, ncomp, nnodes, nqpts,
+                   interp, curl, qref, qweight):
+        """Ceed Hcurl Basis: finite element non tensor-product basis for H(curl)
+             discretizations.
+
+           Args:
+             topo: topology of the element, e.g. hypercube, simplex, etc
+             ncomp: number of field components (1 for scalar fields)
+             nnodes: total number of nodes
+             nqpts: total number of quadrature points
+             *interp: Numpy array holding the row-major (dim * nqpts * nnodes)
+                       matrix expressing the values of basis functions at
+                       quadrature points
+             *curl: Numpy array holding the row-major (curlcomp * nqpts * nnodes),
+                     curlcomp = 1 if dim < 3 else dim, matrix expressing the curl
+                     of basis functions at quadrature points
+             *qref: Numpy array of length (nqpts * dim) holding the locations of
+                     quadrature points on the reference element [-1, 1]
+             *qweight: Numpy array of length nnodes holding the quadrature
+                        weights on the reference element
+
+           Returns:
+             basis: Ceed Basis"""
+
+        return BasisHcurl(self, topo, ncomp, nnodes, nqpts,
+                          interp, curl, qref, qweight)
+
     # CeedQFunction
     def QFunction(self, vlength, f, source):
         """Ceed QFunction: point-wise operation at quadrature points for
diff --git a/python/tests/buildmats.py b/python/tests/buildmats.py
@@ -47,3 +47,80 @@ def buildmats(qref, qweight, mat_dtype="float64"):
         grad[(i + Q) * P + 5] = 2. * (1. * (x2 - 1. / 2.) + x2 * 1.)
 
     return interp, grad
+
+
+def buildmatshdiv(qref, qweight, mat_dtype="float64"):
+    P, Q, dim = 4, 4, 2
+    interp = np.empty(dim * P * Q, dtype=mat_dtype)
+    div = np.empty(P * Q, dtype=mat_dtype)
+
+    qref[0] = -1. / np.sqrt(3.)
+    qref[1] = qref[0]
+    qref[2] = qref[0]
+    qref[3] = -qref[0]
+    qref[4] = -qref[0]
+    qref[5] = -qref[0]
+    qref[6] = qref[0]
+    qref[7] = qref[0]
+    qweight[0] = 1.
+    qweight[1] = 1.
+    qweight[2] = 1.
+    qweight[3] = 1.
+
+    # Loop over quadrature points
+    for i in range(Q):
+        x1 = qref[0 * Q + i]
+        x2 = qref[1 * Q + i]
+        # Interp
+        interp[(i + 0) * P + 0] = 0.
+        interp[(i + Q) * P + 0] = 1. - x2
+        interp[(i + 0) * P + 1] = x1 - 1.
+        interp[(i + Q) * P + 1] = 0.
+        interp[(i + 0) * P + 2] = -x1
+        interp[(i + Q) * P + 2] = 0.
+        interp[(i + 0) * P + 3] = 0.
+        interp[(i + Q) * P + 3] = x2
+        # Div
+        div[i * P + 0] = -1.
+        div[i * P + 1] = 1.
+        div[i * P + 2] = -1.
+        div[i * P + 3] = 1.
+
+    return interp, div
+
+
+def buildmatshcurl(qref, qweight, mat_dtype="float64"):
+    P, Q, dim = 3, 4, 2
+    interp = np.empty(dim * P * Q, dtype=mat_dtype)
+    curl = np.empty(P * Q, dtype=mat_dtype)
+
+    qref[0] = 0.2
+    qref[1] = 0.6
+    qref[2] = 1. / 3.
+    qref[3] = 0.2
+    qref[4] = 0.2
+    qref[5] = 0.2
+    qref[6] = 1. / 3.
+    qref[7] = 0.6
+    qweight[0] = 25. / 96.
+    qweight[1] = 25. / 96.
+    qweight[2] = -27. / 96.
+    qweight[3] = 25. / 96.
+
+    # Loop over quadrature points
+    for i in range(Q):
+        x1 = qref[0 * Q + i]
+        x2 = qref[1 * Q + i]
+        # Interp
+        interp[(i + 0) * P + 0] = -x2
+        interp[(i + Q) * P + 0] = x1
+        interp[(i + 0) * P + 1] = x2
+        interp[(i + Q) * P + 1] = 1. - x1
+        interp[(i + 0) * P + 2] = 1. - x2
+        interp[(i + Q) * P + 2] = x1
+        # Curl
+        curl[i * P + 0] = 2.
+        curl[i * P + 1] = -2.
+        curl[i * P + 2] = -2.
+
+    return interp, curl
diff --git a/python/tests/test-3-basis.py b/python/tests/test-3-basis.py
@@ -303,3 +303,86 @@ def test_323(ceed_resource):
             assert math.fabs(out_array[1 * Q + i] - value) < 10. * TOL
 
 # -------------------------------------------------------------------------------
+# Test interpolation with a 2D Quad non-tensor H(div) basis
+# -------------------------------------------------------------------------------
+
+
+def test_330(ceed_resource):
+    ceed = libceed.Ceed(ceed_resource)
+
+    P, Q, dim = 4, 4, 2
+
+    in_array = np.ones(P, dtype=ceed.scalar_type())
+    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
+    qweight = np.empty(Q, dtype=ceed.scalar_type())
+
+    interp, div = bm.buildmatshdiv(qref, qweight, libceed.scalar_types[
+        libceed.lib.CEED_SCALAR_TYPE])
+
+    b = ceed.BasisHdiv(libceed.QUAD, 1, P, Q, interp, div, qref, qweight)
+
+    # Interpolate function to quadrature points
+    in_vec = ceed.Vector(P)
+    in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
+    out_vec = ceed.Vector(Q * dim)
+    out_vec.set_value(0)
+
+    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
+
+    # Check values at quadrature points
+    with out_vec.array_read() as out_array:
+        for i in range(Q):
+            assert math.fabs(out_array[0 * Q + i] + 1.) < 10. * TOL
+            assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL
+
+# -------------------------------------------------------------------------------
+# Test interpolation with a 2D Simplex non-tensor H(curl) basis
+# -------------------------------------------------------------------------------
+
+
+def test_340(ceed_resource):
+    ceed = libceed.Ceed(ceed_resource)
+
+    P, Q, dim = 3, 4, 2
+
+    in_array = np.empty(P, dtype=ceed.scalar_type())
+    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
+    qweight = np.empty(Q, dtype=ceed.scalar_type())
+
+    interp, curl = bm.buildmatshcurl(qref, qweight, libceed.scalar_types[
+        libceed.lib.CEED_SCALAR_TYPE])
+
+    b = ceed.BasisHcurl(libceed.TRIANGLE, 1, P, Q, interp, curl, qref, qweight)
+
+    # Interpolate function to quadrature points
+    in_array[0] = 1.
+    in_array[1] = 2.
+    in_array[2] = 1.
+
+    in_vec = ceed.Vector(P)
+    in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
+    out_vec = ceed.Vector(Q * dim)
+    out_vec.set_value(0)
+
+    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
+
+    # Check values at quadrature points
+    with out_vec.array_read() as out_array:
+        for i in range(Q):
+            assert math.fabs(out_array[0 * Q + i] - 1.) < 10. * TOL
+
+    # Interpolate function to quadrature points
+    in_array[0] = -1.
+    in_array[1] = 1.
+    in_array[2] = 2.
+
+    out_vec.set_value(0)
+
+    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
+
+    # Check values at quadrature points
+    with out_vec.array_read() as out_array:
+        for i in range(Q):
+            assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL
+
+# -------------------------------------------------------------------------------
diff --git a/rust/libceed/src/lib.rs b/rust/libceed/src/lib.rs
