Address PR comments for high-level APIs

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

Author: sebastiangrimberg

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

@@ -24,6 +24,57 @@ end
     CEED_ERROR_UNSUPPORTED = -3
 end
 
+@cenum CeedMemType::UInt32 begin
+    CEED_MEM_HOST = 0
+    CEED_MEM_DEVICE = 1
+end
+
+@cenum CeedCopyMode::UInt32 begin
+    CEED_COPY_VALUES = 0
+    CEED_USE_POINTER = 1
+    CEED_OWN_POINTER = 2
+end
+
+@cenum CeedNormType::UInt32 begin
+    CEED_NORM_1 = 0
+    CEED_NORM_2 = 1
+    CEED_NORM_MAX = 2
+end
+
+@cenum CeedTransposeMode::UInt32 begin
+    CEED_NOTRANSPOSE = 0
+    CEED_TRANSPOSE = 1
+end
+
+@cenum CeedEvalMode::UInt32 begin
+    CEED_EVAL_NONE = 0
+    CEED_EVAL_INTERP = 1
+    CEED_EVAL_GRAD = 2
+    CEED_EVAL_DIV = 4
+    CEED_EVAL_CURL = 8
+    CEED_EVAL_WEIGHT = 16
+end
+
+@cenum CeedQuadMode::UInt32 begin
+    CEED_GAUSS = 0
+    CEED_GAUSS_LOBATTO = 1
+end
+
+@cenum CeedElemTopology::UInt32 begin
+    CEED_TOPOLOGY_LINE = 65536
+    CEED_TOPOLOGY_TRIANGLE = 131073
+    CEED_TOPOLOGY_QUAD = 131074
+    CEED_TOPOLOGY_TET = 196611
+    CEED_TOPOLOGY_PYRAMID = 196612
+    CEED_TOPOLOGY_PRISM = 196613
+    CEED_TOPOLOGY_HEX = 196614
+end
+
+@cenum CeedContextFieldType::UInt32 begin
+    CEED_CONTEXT_FIELD_DOUBLE = 1
+    CEED_CONTEXT_FIELD_INT32 = 2
+end
+
 mutable struct Ceed_private end
 
 const Ceed = Ptr{Ceed_private}
@@ -123,27 +174,10 @@ function CeedGetScalarType(scalar_type)
     ccall((:CeedGetScalarType, libceed), Cint, (Ptr{CeedScalarType},), scalar_type)
 end
 
-@cenum CeedMemType::UInt32 begin
-    CEED_MEM_HOST = 0
-    CEED_MEM_DEVICE = 1
-end
-
 function CeedGetPreferredMemType(ceed, type)
     ccall((:CeedGetPreferredMemType, libceed), Cint, (Ceed, Ptr{CeedMemType}), ceed, type)
 end
 
-@cenum CeedCopyMode::UInt32 begin
-    CEED_COPY_VALUES = 0
-    CEED_USE_POINTER = 1
-    CEED_OWN_POINTER = 2
-end
-
-@cenum CeedNormType::UInt32 begin
-    CEED_NORM_1 = 0
-    CEED_NORM_2 = 1
-    CEED_NORM_MAX = 2
-end
-
 function CeedVectorCreate(ceed, len, vec)
     ccall((:CeedVectorCreate, libceed), Cint, (Ceed, CeedSize, Ptr{CeedVector}), ceed, len, vec)
 end
@@ -240,11 +274,6 @@ function CeedRequestWait(req)
     ccall((:CeedRequestWait, libceed), Cint, (Ptr{CeedRequest},), req)
 end
 
-@cenum CeedTransposeMode::UInt32 begin
-    CEED_NOTRANSPOSE = 0
-    CEED_TRANSPOSE = 1
-end
-
 function CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr)
     ccall((:CeedElemRestrictionCreate, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedSize, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr)
 end
@@ -325,30 +354,6 @@ function CeedElemRestrictionDestroy(rstr)
     ccall((:CeedElemRestrictionDestroy, libceed), Cint, (Ptr{CeedElemRestriction},), rstr)
 end
 
-@cenum CeedEvalMode::UInt32 begin
-    CEED_EVAL_NONE = 0
-    CEED_EVAL_INTERP = 1
-    CEED_EVAL_GRAD = 2
-    CEED_EVAL_DIV = 4
-    CEED_EVAL_CURL = 8
-    CEED_EVAL_WEIGHT = 16
-end
-
-@cenum CeedQuadMode::UInt32 begin
-    CEED_GAUSS = 0
-    CEED_GAUSS_LOBATTO = 1
-end
-
-@cenum CeedElemTopology::UInt32 begin
-    CEED_TOPOLOGY_LINE = 65536
-    CEED_TOPOLOGY_TRIANGLE = 131073
-    CEED_TOPOLOGY_QUAD = 131074
-    CEED_TOPOLOGY_TET = 196611
-    CEED_TOPOLOGY_PYRAMID = 196612
-    CEED_TOPOLOGY_PRISM = 196613
-    CEED_TOPOLOGY_HEX = 196614
-end
-
 function CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P, Q, quad_mode, basis)
     ccall((:CeedBasisCreateTensorH1Lagrange, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedQuadMode, Ptr{CeedBasis}), ceed, dim, num_comp, P, Q, quad_mode, basis)
 end
@@ -532,11 +537,6 @@ function CeedQFunctionFieldGetEvalMode(qf_field, eval_mode)
     ccall((:CeedQFunctionFieldGetEvalMode, libceed), Cint, (CeedQFunctionField, Ptr{CeedEvalMode}), qf_field, eval_mode)
 end
 
-@cenum CeedContextFieldType::UInt32 begin
-    CEED_CONTEXT_FIELD_DOUBLE = 1
-    CEED_CONTEXT_FIELD_INT32 = 2
-end
-
 # typedef int ( * CeedQFunctionContextDataDestroyUser ) ( void * data )
 const CeedQFunctionContextDataDestroyUser = Ptr{Cvoid}
 
@@ -1364,7 +1364,9 @@ end
 
 # Skipping MacroDefinition: CEED_EXTERN extern CEED_VISIBILITY ( default )
 
-# Skipping MacroDefinition: CEED_QFUNCTION_HELPER CEED_QFUNCTION_ATTR static inline
+# Skipping MacroDefinition: CEED_QFUNCTION_HELPER_ATTR CEED_QFUNCTION_ATTR __attribute__ ( ( always_inline ) )
+
+# Skipping MacroDefinition: CEED_QFUNCTION_HELPER CEED_QFUNCTION_HELPER_ATTR static inline
 
 const CeedInt_FMT = "d"
 
@@ -1374,7 +1376,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