Update Julia/Python/Rust/Fortran bindings

Author: sebastiangrimberg

interface/ceed-fortran.c

@@ -449,6 +449,40 @@ CEED_EXTERN void fCeedBasisCreateH1(int *ceed, int *topo, int *num_comp, int *nn
   }
 }
 
+#define fCeedBasisCreateHdiv FORTRAN_NAME(ceedbasiscreatehdiv, CEEDBASISCREATEHDIV)
+CEED_EXTERN void fCeedBasisCreateHdiv(int *ceed, int *topo, int *num_comp, int *nnodes, int *nqpts, const CeedScalar *interp, const CeedScalar *div,
+                                      const CeedScalar *qref, const CeedScalar *qweight, int *basis, int *err) {
+  if (CeedBasis_count == CeedBasis_count_max) {
+    CeedBasis_count_max += CeedBasis_count_max / 2 + 1;
+    CeedRealloc(CeedBasis_count_max, &CeedBasis_dict);
+  }
+
+  *err = CeedBasisCreateHdiv(Ceed_dict[*ceed], (CeedElemTopology)*topo, *num_comp, *nnodes, *nqpts, interp, div, qref, qweight,
+                             &CeedBasis_dict[CeedBasis_count]);
+
+  if (*err == 0) {
+    *basis = CeedBasis_count++;
+    CeedBasis_n++;
+  }
+}
+
+#define fCeedBasisCreateHcurl FORTRAN_NAME(ceedbasiscreatehcurl, CEEDBASISCREATEHCURL)
+CEED_EXTERN void fCeedBasisCreateHcurl(int *ceed, int *topo, int *num_comp, int *nnodes, int *nqpts, const CeedScalar *interp, const CeedScalar *curl,
+                                       const CeedScalar *qref, const CeedScalar *qweight, int *basis, int *err) {
+  if (CeedBasis_count == CeedBasis_count_max) {
+    CeedBasis_count_max += CeedBasis_count_max / 2 + 1;
+    CeedRealloc(CeedBasis_count_max, &CeedBasis_dict);
+  }
+
+  *err = CeedBasisCreateHcurl(Ceed_dict[*ceed], (CeedElemTopology)*topo, *num_comp, *nnodes, *nqpts, interp, curl, qref, qweight,
+                              &CeedBasis_dict[CeedBasis_count]);
+
+  if (*err == 0) {
+    *basis = CeedBasis_count++;
+    CeedBasis_n++;
+  }
+}
+
 #define fCeedBasisView FORTRAN_NAME(ceedbasisview, CEEDBASISVIEW)
 CEED_EXTERN void fCeedBasisView(int *basis, int *err) { *err = CeedBasisView(CeedBasis_dict[*basis], stdout); }
 

julia/LibCEED.jl/src/Basis.jl

@@ -17,6 +17,8 @@ created using one of:
 - [`create_tensor_h1_lagrange_basis`](@ref)
 - [`create_tensor_h1_basis`](@ref)
 - [`create_h1_basis`](@ref)
+- [`create_hdiv_basis`](@ref)
+- [`create_hcurl_basis`](@ref)
 """
 mutable struct Basis <: AbstractBasis
     ref::RefValue{C.CeedBasis}
@@ -112,7 +114,7 @@ end
 @doc raw"""
     create_h1_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, grad, qref, qweight)
 
-Create a non tensor-product basis for H^1 discretizations
+Create a non tensor-product basis for $H^1$ discretizations
 
 # Arguments:
 - `ceed`:    A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
@@ -166,6 +168,121 @@ function create_h1_basis(
     Basis(ref)
 end
 
+@doc raw"""
+    create_hdiv_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, div, qref, qweight)
+
+Create a non tensor-product basis for H(div) discretizations
+
+# Arguments:
+- `ceed`:    A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
+- `topo`:    [`Topology`](@ref) of 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`:  Matrix of size `(dim, nqpts, nnodes)` expressing the values of basis functions
+             at quadrature points.
+- `div`:     Array of size `(nqpts, nnodes)` expressing divergence of basis functions at
+             quadrature points.
+- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
+             reference element $[-1, 1]$.
+- `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
+             element.
+"""
+function create_hdiv_basis(
+    c::Ceed,
+    topo::Topology,
+    ncomp,
+    nnodes,
+    nqpts,
+    interp::AbstractArray{CeedScalar},
+    div::AbstractArray{CeedScalar},
+    qref::AbstractArray{CeedScalar},
+    qweight::AbstractArray{CeedScalar},
+)
+    dim = getdimension(topo)
+    @assert size(interp) == (dim, nqpts, nnodes)
+    @assert size(div) == (nqpts, nnodes)
+    @assert length(qref) == nqpts
+    @assert length(qweight) == nqpts
+
+    # Convert from Julia matrices and tensors (column-major) to row-major format
+    interp_rowmajor = permutedims(interp, [3, 2, 1])
+    div_rowmajor = collect(div')
+
+    ref = Ref{C.CeedBasis}()
+    C.CeedBasisCreateHdiv(
+        c[],
+        topo,
+        ncomp,
+        nnodes,
+        nqpts,
+        interp_rowmajor,
+        div_rowmajor,
+        qref,
+        qweight,
+        ref,
+    )
+    Basis(ref)
+end
+
+@doc raw"""
+    create_hdiv_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, curl, qref, qweight)
+
+Create a non tensor-product basis for H(div) discretizations
+
+# Arguments:
+- `ceed`:    A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
+- `topo`:    [`Topology`](@ref) of 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`:  Matrix of size `(dim, nqpts, nnodes)` expressing the values of basis functions
+             at quadrature points.
+- `curl`:    Matrix of size `(curlcomp, nqpts, nnodes)`, `curlcomp = 1 if dim < 3 else dim`)
+             matrix expressing curl of basis functions at quadrature points.
+- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
+             reference element $[-1, 1]$.
+- `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
+             element.
+"""
+function create_hdiv_basis(
+    c::Ceed,
+    topo::Topology,
+    ncomp,
+    nnodes,
+    nqpts,
+    interp::AbstractArray{CeedScalar},
+    curl::AbstractArray{CeedScalar},
+    qref::AbstractArray{CeedScalar},
+    qweight::AbstractArray{CeedScalar},
+)
+    dim = getdimension(topo)
+    curlcomp = dim < 3 ? 1 : dim
+    @assert size(interp) == (dim, nqpts, nnodes)
+    @assert size(curl) == (curlcomp, nqpts, nnodes)
+    @assert length(qref) == nqpts
+    @assert length(qweight) == nqpts
+
+    # Convert from Julia matrices and tensors (column-major) to row-major format
+    interp_rowmajor = permutedims(interp, [3, 2, 1])
+    curl_rowmajor = permutedims(curl, [3, 2, 1])
+
+    ref = Ref{C.CeedBasis}()
+    C.CeedBasisCreateHcurl(
+        c[],
+        topo,
+        ncomp,
+        nnodes,
+        nqpts,
+        interp_rowmajor,
+        curl_rowmajor,
+        qref,
+        qweight,
+        ref,
+    )
+    Basis(ref)
+end
+
 """
     apply!(b::Basis, nelem, tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector)
 
@@ -353,7 +470,13 @@ function getinterp(b::Basis)
     C.CeedBasisGetInterp(b[], ref)
     q = getnumqpts(b)
     p = getnumnodes(b)
-    collect(unsafe_wrap(Array, ref[], (p, q))')
+    qcomp = Ref{CeedInt}()
+    C.CeedBasisGetNumQuadratureComponents(b[], C.CEED_EVAL_INTERP, qcomp)
+    if qcomp == 1
+        collect(unsafe_wrap(Array, ref[], (p, q))')
+    else
+        permutedims(unsafe_wrap(Array, ref[], (p, q, qcomp)), [3, 2, 1])
+    end
 end
 
 """
@@ -399,3 +522,34 @@ function getgrad1d(b::Basis)
     p = getnumnodes1d(b)
     collect(unsafe_wrap(Array, ref[], (p, q))')
 end
+
+"""
+    getdiv(b::Basis)
+
+Get the divergence matrix of the given [`Basis`](@ref). Returns a tensor of size
+`(getnumqpts(b), getnumnodes(b))`.
+"""
+function getdiv(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetDiv(b[], ref)
+    q = getnumqpts(b)
+    p = getnumnodes(b)
+    collect(unsafe_wrap(Array, ref[], (p, q))')
+end
+
+"""
+    getcurl(b::Basis)
+
+Get the curl matrix of the given [`Basis`](@ref). Returns a tensor of size
+`(curlcomp, getnumqpts(b), getnumnodes(b))`, `curlcomp = 1 if getdimension(b) < 3 else
+getdimension(b)`.
+"""
+function getcurl(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetCurl(b[], ref)
+    q = getnumqpts(b)
+    p = getnumnodes(b)
+    qcomp = Ref{CeedInt}()
+    C.CeedBasisGetNumQuadratureComponents(b[], C.CEED_EVAL_CURL, qcomp)
+    permutedims(unsafe_wrap(Array, ref[], (p, q, qcomp)), [3, 2, 1])
+end

julia/LibCEED.jl/src/generated/libceed_bindings.jl

@@ -152,6 +152,10 @@ function CeedVectorReferenceCopy(vec, vec_copy)
     ccall((:CeedVectorReferenceCopy, libceed), Cint, (CeedVector, Ptr{CeedVector}), vec, vec_copy)
 end
 
+function CeedVectorCopy(vec, vec_copy)
+    ccall((:CeedVectorCopy, libceed), Cint, (CeedVector, CeedVector), vec, vec_copy)
+end
+
 function CeedVectorSetArray(vec, mem_type, copy_mode, array)
     ccall((:CeedVectorSetArray, libceed), Cint, (CeedVector, CeedMemType, CeedCopyMode, Ptr{CeedScalar}), vec, mem_type, copy_mode, array)
 end
@@ -200,6 +204,10 @@ function CeedVectorAXPY(y, alpha, x)
     ccall((:CeedVectorAXPY, libceed), Cint, (CeedVector, CeedScalar, CeedVector), y, alpha, x)
 end
 
+function CeedVectorAXPBY(y, alpha, beta, x)
+    ccall((:CeedVectorAXPBY, libceed), Cint, (CeedVector, CeedScalar, CeedScalar, CeedVector), y, alpha, beta, x)
+end
+
 function CeedVectorPointwiseMult(w, x, y)
     ccall((:CeedVectorPointwiseMult, libceed), Cint, (CeedVector, CeedVector, CeedVector), w, x, y)
 end
@@ -208,6 +216,10 @@ function CeedVectorReciprocal(vec)
     ccall((:CeedVectorReciprocal, libceed), Cint, (CeedVector,), vec)
 end
 
+function CeedVectorViewRange(vec, start, stop, step, fp_fmt, stream)
+    ccall((:CeedVectorViewRange, libceed), Cint, (CeedVector, CeedSize, CeedSize, CeedInt, Ptr{Cchar}, Ptr{Libc.FILE}), vec, start, stop, step, fp_fmt, stream)
+end
+
 function CeedVectorView(vec, fp_fmt, stream)
     ccall((:CeedVectorView, libceed), Cint, (CeedVector, Ptr{Cchar}, Ptr{Libc.FILE}), vec, fp_fmt, stream)
 end
@@ -353,6 +365,10 @@ function CeedBasisCreateHdiv(ceed, topo, num_comp, num_nodes, nqpts, interp, div
     ccall((:CeedBasisCreateHdiv, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, num_comp, num_nodes, nqpts, interp, div, q_ref, q_weights, basis)
 end
 
+function CeedBasisCreateHcurl(ceed, topo, num_comp, num_nodes, nqpts, interp, curl, q_ref, q_weights, basis)
+    ccall((:CeedBasisCreateHcurl, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, num_comp, num_nodes, nqpts, interp, curl, q_ref, q_weights, basis)
+end
+
 function CeedBasisCreateProjection(basis_from, basis_to, basis_project)
     ccall((:CeedBasisCreateProjection, libceed), Cint, (CeedBasis, CeedBasis, Ptr{CeedBasis}), basis_from, basis_to, basis_project)
 end
@@ -381,10 +397,6 @@ function CeedBasisGetTopology(basis, topo)
     ccall((:CeedBasisGetTopology, libceed), Cint, (CeedBasis, Ptr{CeedElemTopology}), basis, topo)
 end
 
-function CeedBasisGetNumQuadratureComponents(basis, Q_comp)
-    ccall((:CeedBasisGetNumQuadratureComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q_comp)
-end
-
 function CeedBasisGetNumComponents(basis, num_comp)
     ccall((:CeedBasisGetNumComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, num_comp)
 end
@@ -433,6 +445,10 @@ function CeedBasisGetDiv(basis, div)
     ccall((:CeedBasisGetDiv, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, div)
 end
 
+function CeedBasisGetCurl(basis, curl)
+    ccall((:CeedBasisGetCurl, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, curl)
+end
+
 function CeedBasisDestroy(basis)
     ccall((:CeedBasisDestroy, libceed), Cint, (Ptr{CeedBasis},), basis)
 end
@@ -716,32 +732,36 @@ function CeedOperatorGetFlopsEstimate(op, flops)
     ccall((:CeedOperatorGetFlopsEstimate, libceed), Cint, (CeedOperator, Ptr{CeedSize}), op, flops)
 end
 
-function CeedOperatorContextGetFieldLabel(op, field_name, field_label)
-    ccall((:CeedOperatorContextGetFieldLabel, libceed), Cint, (CeedOperator, Ptr{Cchar}, Ptr{CeedContextFieldLabel}), op, field_name, field_label)
+function CeedOperatorGetContext(op, ctx)
+    ccall((:CeedOperatorGetContext, libceed), Cint, (CeedOperator, Ptr{CeedQFunctionContext}), op, ctx)
 end
 
-function CeedOperatorContextSetDouble(op, field_label, values)
-    ccall((:CeedOperatorContextSetDouble, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cdouble}), op, field_label, values)
+function CeedOperatorGetContextFieldLabel(op, field_name, field_label)
+    ccall((:CeedOperatorGetContextFieldLabel, libceed), Cint, (CeedOperator, Ptr{Cchar}, Ptr{CeedContextFieldLabel}), op, field_name, field_label)
 end
 
-function CeedOperatorContextGetDoubleRead(op, field_label, num_values, values)
-    ccall((:CeedOperatorContextGetDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cdouble}}), op, field_label, num_values, values)
+function CeedOperatorSetContextDouble(op, field_label, values)
+    ccall((:CeedOperatorSetContextDouble, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cdouble}), op, field_label, values)
 end
 
-function CeedOperatorContextRestoreDoubleRead(op, field_label, values)
-    ccall((:CeedOperatorContextRestoreDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cdouble}}), op, field_label, values)
+function CeedOperatorGetContextDoubleRead(op, field_label, num_values, values)
+    ccall((:CeedOperatorGetContextDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cdouble}}), op, field_label, num_values, values)
 end
 
-function CeedOperatorContextSetInt32(op, field_label, values)
-    ccall((:CeedOperatorContextSetInt32, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cint}), op, field_label, values)
+function CeedOperatorRestoreContextDoubleRead(op, field_label, values)
+    ccall((:CeedOperatorRestoreContextDoubleRead, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cdouble}}), op, field_label, values)
 end
 
-function CeedOperatorContextGetInt32Read(op, field_label, num_values, values)
-    ccall((:CeedOperatorContextGetInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cint}}), op, field_label, num_values, values)
+function CeedOperatorSetContextInt32(op, field_label, values)
+    ccall((:CeedOperatorSetContextInt32, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Cint}), op, field_label, values)
 end
 
-function CeedOperatorContextRestoreInt32Read(op, field_label, values)
-    ccall((:CeedOperatorContextRestoreInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cint}}), op, field_label, values)
+function CeedOperatorGetContextInt32Read(op, field_label, num_values, values)
+    ccall((:CeedOperatorGetContextInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Csize_t}, Ptr{Ptr{Cint}}), op, field_label, num_values, values)
+end
+
+function CeedOperatorRestoreContextInt32Read(op, field_label, values)
+    ccall((:CeedOperatorRestoreContextInt32Read, libceed), Cint, (CeedOperator, CeedContextFieldLabel, Ptr{Ptr{Cint}}), op, field_label, values)
 end
 
 function CeedOperatorApply(op, in, out, request)
@@ -987,6 +1007,7 @@ end
 @cenum CeedFESpace::UInt32 begin
     CEED_FE_SPACE_H1 = 1
     CEED_FE_SPACE_HDIV = 2
+    CEED_FE_SPACE_HCURL = 3
 end
 
 function CeedBasisGetCollocatedGrad(basis, colo_grad_1d)
@@ -1009,10 +1030,18 @@ function CeedBasisReference(basis)
     ccall((:CeedBasisReference, libceed), Cint, (CeedBasis,), basis)
 end
 
+function CeedBasisGetNumQuadratureComponents(basis, eval_mode, q_comp)
+    ccall((:CeedBasisGetNumQuadratureComponents, libceed), Cint, (CeedBasis, CeedEvalMode, Ptr{CeedInt}), basis, eval_mode, q_comp)
+end
+
 function CeedBasisGetFlopsEstimate(basis, t_mode, eval_mode, flops)
     ccall((:CeedBasisGetFlopsEstimate, libceed), Cint, (CeedBasis, CeedTransposeMode, CeedEvalMode, Ptr{CeedSize}), basis, t_mode, eval_mode, flops)
 end
 
+function CeedBasisGetFESpace(basis, fe_space)
+    ccall((:CeedBasisGetFESpace, libceed), Cint, (CeedBasis, Ptr{CeedFESpace}), basis, fe_space)
+end
+
 function CeedBasisGetTopologyDimension(topo, dim)
     ccall((:CeedBasisGetTopologyDimension, libceed), Cint, (CeedElemTopology, Ptr{CeedInt}), topo, dim)
 end
@@ -1033,6 +1062,10 @@ function CeedTensorContractApply(contract, A, B, C, J, t, t_mode, Add, u, v)
     ccall((:CeedTensorContractApply, libceed), Cint, (CeedTensorContract, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), contract, A, B, C, J, t, t_mode, Add, u, v)
 end
 
+function CeedTensorContractStridedApply(contract, A, B, C, D, J, t, t_mode, add, u, v)
+    ccall((:CeedTensorContractStridedApply, libceed), Cint, (CeedTensorContract, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), contract, A, B, C, D, J, t, t_mode, add, u, v)
+end
+
 function CeedTensorContractGetCeed(contract, ceed)
     ccall((:CeedTensorContractGetCeed, libceed), Cint, (CeedTensorContract, Ptr{Ceed}), contract, ceed)
 end