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/docs/src/Basis.md

@@ -14,11 +14,13 @@ BasisCollocated
 create_tensor_h1_lagrange_basis
 create_tensor_h1_basis
 create_h1_basis
+create_hdiv_basis
+create_hcurl_basis
 apply!(b::Basis, nelem, tmode::TransposeMode, emode::EvalMode, u::LibCEED.AbstractCeedVector, v::LibCEED.AbstractCeedVector)
 apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP)
 getdimension
 gettopology
-getnumcomponents(b::Basis)
+getnumcomponents
 getnumnodes
 getnumnodes1d
 getnumqpts
@@ -29,4 +31,6 @@ getinterp
 getinterp1d
 getgrad
 getgrad1d
+getdiv
+getcurl
 ```

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_hcurl_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, curl, qref, qweight)
+
+Create a non tensor-product basis for H(curl) 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_hcurl_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)
 
@@ -342,18 +459,25 @@ function getqweights(b::Basis)
     copy(unsafe_wrap(Array, ref[], istensor[] ? getnumqpts1d(b) : getnumqpts(b)))
 end
 
-"""
+@doc raw"""
     getinterp(b::Basis)
 
 Get the interpolation matrix of the given [`Basis`](@ref). Returns a matrix of size
-`(getnumqpts(b), getnumnodes(b))`.
+`(getnumqpts(b), getnumnodes(b))` for a given $H^1$ basis or `(getdimension(b),
+getnumqpts(b), getnumnodes(b))` for a given vector $H(div)$ or $H(curl)$ basis.
 """
 function getinterp(b::Basis)
     ref = Ref{Ptr{CeedScalar}}()
     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 +523,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