Address PR comments for high-level APIs

Add Rust sibling functions and tests, add Python tests, fix Julia versioning for CI.

.github/workflows/julia-test-with-style.yml

@@ -33,9 +33,9 @@ jobs:
           make -j2
           LIBCEED_LIB=$(find $PWD/lib -name "libceed.*")
           pushd julia/LibCEED.jl
-          echo >> test/Project.toml
-          echo "[preferences.libCEED_jll]" >> test/Project.toml
-          echo "libceed_path = \"$LIBCEED_LIB\"" >> test/Project.toml
+          echo >> Project.toml
+          echo "[preferences.libCEED_jll]" >> Project.toml
+          echo "libceed_path = \"$LIBCEED_LIB\"" >> Project.toml
           [[ "$GITHUB_REF" =~ ^refs/(heads/release|tags/).* ]] || julia --project -e 'import Pkg; Pkg.test("LibCEED"; coverage=true, test_args=["--run-dev-tests"])'
-          git checkout test/Project.toml && julia --project -e 'import Pkg; Pkg.test("LibCEED")'
+          git checkout Project.toml && julia --project -e 'import Pkg; Pkg.test("LibCEED")'
           julia --project=.style/ -e 'import Pkg; Pkg.instantiate()' && julia --project=.style/ .style/ceed_style.jl && git diff --exit-code src test examples

include/ceed/ceed.h

@@ -137,7 +137,7 @@ CEED_EXTERN int CeedResetErrorMessage(Ceed, const char **err_msg);
 #define CEED_VERSION_MAJOR 0
 #define CEED_VERSION_MINOR 11
 #define CEED_VERSION_PATCH 0
-#define CEED_VERSION_RELEASE true
+#define CEED_VERSION_RELEASE false
 
 /// Compile-time check that the the current library version is at least as recent as the specified version.
 /// This macro is typically used in

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)

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 )
 

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",

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

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

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
+
+# -------------------------------------------------------------------------------