From 44554ea01e90fce366fc2a203c44be15754a38d6 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Tue, 22 Sep 2020 14:47:35 -0700 Subject: [PATCH] Julia: LibCEED.jl Julia interface --- AUTHORS | 1 + julia/LibCEED.jl/.gitignore | 4 + julia/LibCEED.jl/Project.toml | 20 + julia/LibCEED.jl/README.md | 49 ++ julia/LibCEED.jl/deps/build.jl | 18 + julia/LibCEED.jl/src/Basis.jl | 391 +++++++++ julia/LibCEED.jl/src/C.jl | 45 + julia/LibCEED.jl/src/Ceed.jl | 113 +++ julia/LibCEED.jl/src/CeedVector.jl | 308 +++++++ julia/LibCEED.jl/src/Context.jl | 56 ++ julia/LibCEED.jl/src/Cuda.jl | 155 ++++ julia/LibCEED.jl/src/ElemRestriction.jl | 288 +++++++ julia/LibCEED.jl/src/Globals.jl | 117 +++ julia/LibCEED.jl/src/LibCEED.jl | 150 ++++ julia/LibCEED.jl/src/Misc.jl | 119 +++ julia/LibCEED.jl/src/Operator.jl | 149 ++++ julia/LibCEED.jl/src/QFunction.jl | 148 ++++ julia/LibCEED.jl/src/Quadrature.jl | 37 + julia/LibCEED.jl/src/Request.jl | 20 + julia/LibCEED.jl/src/UserQFunction.jl | 278 +++++++ .../src/generated/generate_bindings.jl | 28 + julia/LibCEED.jl/src/generated/libceed_api.jl | 773 ++++++++++++++++++ .../src/generated/libceed_common.jl | 99 +++ 23 files changed, 3366 insertions(+) create mode 100644 julia/LibCEED.jl/.gitignore create mode 100644 julia/LibCEED.jl/Project.toml create mode 100644 julia/LibCEED.jl/README.md create mode 100644 julia/LibCEED.jl/deps/build.jl create mode 100644 julia/LibCEED.jl/src/Basis.jl create mode 100644 julia/LibCEED.jl/src/C.jl create mode 100644 julia/LibCEED.jl/src/Ceed.jl create mode 100644 julia/LibCEED.jl/src/CeedVector.jl create mode 100644 julia/LibCEED.jl/src/Context.jl create mode 100644 julia/LibCEED.jl/src/Cuda.jl create mode 100644 julia/LibCEED.jl/src/ElemRestriction.jl create mode 100644 julia/LibCEED.jl/src/Globals.jl create mode 100644 julia/LibCEED.jl/src/LibCEED.jl create mode 100644 julia/LibCEED.jl/src/Misc.jl create mode 100644 julia/LibCEED.jl/src/Operator.jl create mode 100644 julia/LibCEED.jl/src/QFunction.jl create mode 100644 julia/LibCEED.jl/src/Quadrature.jl create mode 100644 julia/LibCEED.jl/src/Request.jl create mode 100644 julia/LibCEED.jl/src/UserQFunction.jl create mode 100644 julia/LibCEED.jl/src/generated/generate_bindings.jl create mode 100644 julia/LibCEED.jl/src/generated/libceed_api.jl create mode 100644 julia/LibCEED.jl/src/generated/libceed_common.jl diff --git a/AUTHORS b/AUTHORS index c8e8e130f3..40532a43d6 100644 --- a/AUTHORS +++ b/AUTHORS @@ -8,6 +8,7 @@ Yohann Dudouit Leila Ghaffari Tzanio Kolev David Medina +Will Pazner Thilina Rathnayake Jeremy L. Thompson Stan Tomov diff --git a/julia/LibCEED.jl/.gitignore b/julia/LibCEED.jl/.gitignore new file mode 100644 index 0000000000..3e3da309a2 --- /dev/null +++ b/julia/LibCEED.jl/.gitignore @@ -0,0 +1,4 @@ +Manifest.toml +deps/build.log +deps/deps.jl +docs/build diff --git a/julia/LibCEED.jl/Project.toml b/julia/LibCEED.jl/Project.toml new file mode 100644 index 0000000000..63caac695d --- /dev/null +++ b/julia/LibCEED.jl/Project.toml @@ -0,0 +1,20 @@ +name = "LibCEED" +uuid = "2cd74e05-b976-4426-91fa-5f1011f8952b" +authors = ["Will Pazner "] +version = "0.1.0" + +[deps] +CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82" +Cassette = "7057c7e9-c182-5462-911a-8362d720325c" +Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6" +libCEED_jll = "762fde13-7596-547b-826d-8223c52d51c1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/julia/LibCEED.jl/README.md b/julia/LibCEED.jl/README.md new file mode 100644 index 0000000000..38c1032407 --- /dev/null +++ b/julia/LibCEED.jl/README.md @@ -0,0 +1,49 @@ +# LibCEED.jl: Julia Interface for [libCEED](https://github.com/CEED/libCEED) + +## Installation + +When the LibCEED.jl package is built, it requires the environment variable +`JULIA_LIBCEED_LIB` to be set to the location of the compiled libCEED shared +library. + +For example, the package can be installed by: +```julia +% JULIA_LIBCEED_LIB=/path/to/libceed.so julia +julia> # press ] to enter package manager + +(@v1.5) pkg> add LibCEED +``` +or, equivalently, +```julia +% julia + +julia> withenv("JULIA_LIBCEED_LIB" => "/path/to/libceed.so") do + Pkg.add("LibCEED") +end +``` + + +## Usage + +This package provides both a low-level and high-level interface for libCEED. + +### Low-Level Interface + +The low-level interface (provided in the `LibCEED.C` module) is in one-to-one +correspondence with the C libCEED iterface, and is automatically generated (with +some minor manual modifications) using the Julia package Clang.jl. The script +used to generate bindings is available in `generate_bindings.jl`. + +With the low-level interface, the user is responsible for freeing all allocated +memory (calling the appropriate `Ceed*Destroy` functions). This interface is +not type-safe, and calling functions with the wrong arguments can cause libCEED +to crash. + +### High-Level Interface + +The high-level interface provides a more idiomatic Julia interface to the +libCEED library. Objects allocated using the high-level interface will +automatically be destroyed by the garbage collector, so the user does not need +to manually manage memory. + +See the documentation for more information. diff --git a/julia/LibCEED.jl/deps/build.jl b/julia/LibCEED.jl/deps/build.jl new file mode 100644 index 0000000000..118e024a67 --- /dev/null +++ b/julia/LibCEED.jl/deps/build.jl @@ -0,0 +1,18 @@ +const libceed_envvar = "JULIA_LIBCEED_LIB" + +if isfile("deps.jl") + rm("deps.jl") +end + +if haskey(ENV, libceed_envvar) + ceedpath = ENV[libceed_envvar] + @info "Using libCEED library specified by $libceed_envvar at $ceedpath" + open("deps.jl", write=true) do f + println(f, "const libceed = \"$(escape_string(ceedpath))\"") + end +else + @info "Using prebuilt libCEED binaries provided by libCEED_jll" + open("deps.jl", write=true) do f + println(f, "using libCEED_jll") + end +end diff --git a/julia/LibCEED.jl/src/Basis.jl b/julia/LibCEED.jl/src/Basis.jl new file mode 100644 index 0000000000..fdb55b79cf --- /dev/null +++ b/julia/LibCEED.jl/src/Basis.jl @@ -0,0 +1,391 @@ +abstract type AbstractBasis end + +""" + BasisCollocated() + +Returns the singleton object corresponding to libCEED's `CEED_BASIS_COLLOCATED`. +""" +struct BasisCollocated <: AbstractBasis end +Base.getindex(::BasisCollocated) = C.CEED_BASIS_COLLOCATED[] + +""" + Basis + +Wraps a `CeedBasis` object, representing a finite element basis. A `Basis` object can be +created using one of: + +- [`create_tensor_h1_lagrange_basis`](@ref) +- [`create_tensor_h1_basis`](@ref) +- [`create_h1_basis`](@ref) +""" +mutable struct Basis <: AbstractBasis + ref::RefValue{C.CeedBasis} + function Basis(ref) + obj = new(ref) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + destroy(x) + end + return obj + end +end +destroy(b::Basis) = C.CeedBasisDestroy(b.ref) # COV_EXCL_LINE +Base.getindex(b::Basis) = b.ref[] +Base.show(io::IO, ::MIME"text/plain", b::Basis) = ceed_show(io, b, C.CeedBasisView) + +@doc raw""" + create_tensor_h1_lagrange_basis(ceed, dim, ncomp, p, q, qmode) + +Create a tensor-product Lagrange basis. + +# Arguments: +- `ceed`: A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created. +- `dim`: Topological dimension of element. +- `ncomp`: Number of field components (1 for scalar fields). +- `p`: Number of Gauss-Lobatto nodes in one dimension. The polynomial degree of the + resulting $Q_k$ element is $k=p-1$. +- `q`: Number of quadrature points in one dimension. +- `qmode`: Distribution of the $q$ quadrature points (affects order of accuracy for the + quadrature). +""" +function create_tensor_h1_lagrange_basis(c::Ceed, dim, ncomp, p, q, quad_mode::QuadMode) + ref = Ref{C.CeedBasis}() + C.CeedBasisCreateTensorH1Lagrange(c[], dim, ncomp, p, q, quad_mode, ref) + Basis(ref) +end + +@doc raw""" + create_tensor_h1_basis(c::Ceed, dim, ncomp, p, q, interp1d, grad1d, qref1d, qweight1d) + +Create a tensor-product basis for $H^1$ discretizations. + +# Arguments: +- `ceed`: A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created. +- `dim`: Topological dimension. +- `ncomp`: Number of field components (1 for scalar fields). +- `p`: Number of nodes in one dimension. +- `q`: Number of quadrature points in one dimension +- `interp1d`: Matrix of size `(q, p)` expressing the values of nodal basis functions at + quadrature points. +- `grad1d`: Matrix of size `(p, q)` expressing derivatives of nodal basis functions at + quadrature points. +- `qref1d`: Array of length `q` holding the locations of quadrature points on the 1D + reference element $[-1, 1]$. +- `qweight1d`: Array of length `q` holding the quadrature weights on the reference element. +""" +function create_tensor_h1_basis( + c::Ceed, + dim, + ncomp, + p, + q, + interp1d, + grad1d, + qref1d, + qweight1d, +) + @assert size(interp1d) == (q, p) + @assert size(grad1d) == (q, p) + @assert length(qref1d) == q + @assert length(qweight1d) == q + + # Convert from Julia matrices (column-major) to row-major format + interp1d_rowmajor = collect(interp1d') + grad1d_rowmajor = collect(grad1d') + + ref = Ref{C.CeedBasis}() + C.CeedBasisCreateTensorH1( + c[], + dim, + ncomp, + p, + q, + interp1d_rowmajor, + grad1d_rowmajor, + qref1d, + qweight1d, + ref, + ) + Basis(ref) +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 + +# 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 `(nqpts, nnodes)` expressing the values of nodal basis functions + at quadrature points. +- `grad`: Array of size `(dim, nqpts, nnodes)` expressing derivatives of nodal 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_h1_basis( + c::Ceed, + topo::Topology, + ncomp, + nnodes, + nqpts, + interp, + grad, + qref, + qweight, +) + dim = getdimension(topo) + @assert size(interp) == (nqpts, nnodes) + @assert size(grad) == (dim, nqpts, nnodes) + @assert length(qref) == nqpts + @assert length(qweight) == nqpts + + # Convert from Julia matrices and tensors (column-major) to row-major format + interp_rowmajor = collect(interp') + grad_rowmajor = permutedims(grad, [3, 2, 1]) + + ref = Ref{C.CeedBasis}() + C.CeedBasisCreateH1( + c[], + topo, + ncomp, + nnodes, + nqpts, + interp_rowmajor, + grad_rowmajor, + qref, + qweight, + ref, + ) + Basis(ref) +end + +""" + apply!(b::Basis, nelem, tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector) + +Apply basis evaluation from nodes to quadrature points or vice versa, storing the result in +the [`CeedVector`](@ref) `v`. + +`nelem` specifies the number of elements to apply the basis evaluation to; the backend will +specify the ordering in CeedElemRestrictionCreateBlocked() + +Set `tmode` to `CEED_NOTRANSPOSE` to evaluate from nodes to quadrature or to +`CEED_TRANSPOSE` to apply the transpose, mapping from quadrature points to nodes. + +Set the [`EvalMode`](@ref) `emode` to: +- `CEED_EVAL_NONE` to use values directly, +- `CEED_EVAL_INTERP` to use interpolated values, +- `CEED_EVAL_GRAD` to use gradients, +- `CEED_EVAL_WEIGHT` to use quadrature weights. +""" +function apply!( + b::Basis, + nelem, + tmode::TransposeMode, + emode::EvalMode, + u::AbstractCeedVector, + v::AbstractCeedVector, +) + C.CeedBasisApply(b[], nelem, tmode, emode, u[], v[]) +end + +""" + apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP) + +Performs the same function as the above-defined [`apply!`](@ref apply!(b::Basis, nelem, +tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector)), but +automatically convert from Julia arrays to [`CeedVector`](@ref) for convenience. + +The result will be returned in a newly allocated array of the correct size. +""" +function apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP) + ceed_ref = Ref{C.Ceed}() + ccall((:CeedBasisGetCeed, C.libceed), Cint, (C.CeedBasis, Ptr{C.Ceed}), b[], ceed_ref) + c = Ceed(ceed_ref) + + u_vec = CeedVector(c, u) + + len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : getnumqpts(b) + if emode == EVAL_GRAD + len_v *= getdimension(b) + end + + v_vec = CeedVector(c, len_v) + + apply!(b, nelem, tmode, emode, u_vec, v_vec) + Vector(v_vec) +end + +""" + getdimension(b::Basis) + +Return the spatial dimension of the given [`Basis`](@ref). +""" +function getdimension(b::Basis) + dim = Ref{CeedInt}() + C.CeedBasisGetDimension(b[], dim) + dim[] +end + +""" + getdimension(t::Topology) + +Return the spatial dimension of the given [`Topology`](@ref). +""" +function getdimension(t::Topology) + return Int(t) >> 16 +end + +""" + gettopology(b::Basis) + +Return the [`Topology`](@ref) of the given [`Basis`](@ref). +""" +function gettopology(b::Basis) + topo = Ref{Topology}() + C.CeedBasisGetTopology(b[], topo) + topo[] +end + +""" + getnumcomponents(b::Basis) + +Return the number of components of the given [`Basis`](@ref). +""" +function getnumcomponents(b::Basis) + ncomp = Ref{CeedInt}() + C.CeedBasisGetNumComponents(b[], ncomp) + ncomp[] +end + +""" + getnumnodes(b::Basis) + +Return the number of nodes of the given [`Basis`](@ref). +""" +function getnumnodes(b::Basis) + nnodes = Ref{CeedInt}() + C.CeedBasisGetNumNodes(b[], nnodes) + nnodes[] +end + +""" + getnumnodes1d(b::Basis) + + Return the number of 1D nodes of the given (tensor-product) [`Basis`](@ref). +""" +function getnumnodes1d(b::Basis) + nnodes1d = Ref{CeedInt}() + C.CeedBasisGetNumNodes1D(b[], nnodes1d) + nnodes1d[] +end + +""" + getnumqpts(b::Basis) + +Return the number of quadrature points of the given [`Basis`](@ref). +""" +function getnumqpts(b::Basis) + nqpts = Ref{CeedInt}() + C.CeedBasisGetNumQuadraturePoints(b[], nqpts) + nqpts[] +end + +""" + getnumqpts1d(b::Basis) + +Return the number of 1D quadrature points of the given (tensor-product) [`Basis`](@ref). +""" +function getnumqpts1d(b::Basis) + nqpts1d = Ref{CeedInt}() + C.CeedBasisGetNumQuadraturePoints1D(b[], nqpts1d) + nqpts1d[] +end + +""" + getqref(b::Basis) + +Get the reference coordinates of quadrature points (in `dim` dimensions) of the given +[`Basis`](@ref). +""" +function getqref(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetQRef(b[], ref) + copy(unsafe_wrap(Array, ref[], getnumqpts(b))) +end + +""" + getqref(b::Basis) + +Get the quadrature weights of quadrature points (in `dim` dimensions) of the given +[`Basis`](@ref). +""" +function getqweights(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetQWeights(b[], ref) + copy(unsafe_wrap(Array, ref[], getnumqpts(b))) +end + +""" + getinterp(b::Basis) + +Get the interpolation matrix of the given [`Basis`](@ref). Returns a matrix of size +`(getnumqpts(b), getnumnodes(b))`. +""" +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))') +end + +""" + getinterp1d(b::Basis) + +Get the 1D interpolation matrix of the given [`Basis`](@ref). `b` must be a tensor-product +basis, otherwise this function will fail. Returns a matrix of size `(getnumqpts1d(b), +getnumnodes1d(b))`. +""" +function getinterp1d(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetInterp1D(b[], ref) + q = getnumqpts1d(b) + p = getnumnodes1d(b) + collect(unsafe_wrap(Array, ref[], (p, q))') +end + +""" + getgad(b::Basis) + +Get the gradient matrix of the given [`Basis`](@ref). Returns a tensor of size +`(getdimension(b), getnumqpts(b), getnumnodes(b))`. +""" +function getgrad(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetGrad(b[], ref) + dim = getdimension(b) + q = getnumqpts(b) + p = getnumnodes(b) + permutedims(unsafe_wrap(Array, ref[], (p, q, dim)), [3, 2, 1]) +end + +""" + getgrad1d(b::Basis) + +Get the 1D derivative matrix of the given [`Basis`](@ref). Returns a matrix of size +`(getnumqpts(b), getnumnodes(b))`. +""" +function getgrad1d(b::Basis) + ref = Ref{Ptr{CeedScalar}}() + C.CeedBasisGetGrad1D(b[], ref) + q = getnumqpts1d(b) + p = getnumnodes1d(b) + collect(unsafe_wrap(Array, ref[], (p, q))') +end diff --git a/julia/LibCEED.jl/src/C.jl b/julia/LibCEED.jl/src/C.jl new file mode 100644 index 0000000000..40e2114d03 --- /dev/null +++ b/julia/LibCEED.jl/src/C.jl @@ -0,0 +1,45 @@ +# Low-level C API for libCEED + +module C + +using CEnum, Libdl + +# Get the path to the libCEED dynamic library, configured during the build step +# of the LibCEED.jl package. +const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl") +if !isfile(depsfile) + error("LibCEED.jl not properly installed. Please run Pkg.build(\"LibCEED\")") +end +include(depsfile) + +include(joinpath(@__DIR__, "generated", "libceed_common.jl")) +include(joinpath(@__DIR__, "generated", "libceed_api.jl")) + +const CEED_STRIDES_BACKEND = Ref{Ptr{CeedInt}}() +const CEED_BASIS_COLLOCATED = Ref{CeedBasis}() +const CEED_VECTOR_ACTIVE = Ref{CeedVector}() +const CEED_VECTOR_NONE = Ref{CeedVector}() +const CEED_ELEMRESTRICTION_NONE = Ref{CeedElemRestriction}() +const CEED_QFUNCTION_NONE = Ref{CeedQFunction}() +const CEED_REQUEST_IMMEDIATE = Ref{CeedRequest}() +const CEED_REQUEST_ORDERED = Ref{CeedRequest}() + +function __init__() + global libceed_handle = dlopen(libceed) + # some global variables + CEED_STRIDES_BACKEND[] = cglobal((:CEED_STRIDES_BACKEND, libceed)) + CEED_BASIS_COLLOCATED[] = + unsafe_load(cglobal((:CEED_BASIS_COLLOCATED, libceed), CeedBasis)) + CEED_VECTOR_ACTIVE[] = unsafe_load(cglobal((:CEED_VECTOR_ACTIVE, libceed), CeedVector)) + CEED_VECTOR_NONE[] = unsafe_load(cglobal((:CEED_VECTOR_NONE, libceed), CeedVector)) + CEED_ELEMRESTRICTION_NONE[] = + unsafe_load(cglobal((:CEED_ELEMRESTRICTION_NONE, libceed), CeedElemRestriction)) + CEED_QFUNCTION_NONE[] = + unsafe_load(cglobal((:CEED_QFUNCTION_NONE, libceed), CeedQFunction)) + CEED_REQUEST_IMMEDIATE[] = + unsafe_load(cglobal((:CEED_REQUEST_IMMEDIATE, libceed), Ptr{CeedRequest})) + CEED_REQUEST_ORDERED[] = + unsafe_load(cglobal((:CEED_REQUEST_ORDERED, libceed), Ptr{CeedRequest})) +end + +end # module diff --git a/julia/LibCEED.jl/src/Ceed.jl b/julia/LibCEED.jl/src/Ceed.jl new file mode 100644 index 0000000000..d919e75812 --- /dev/null +++ b/julia/LibCEED.jl/src/Ceed.jl @@ -0,0 +1,113 @@ +struct CeedError <: Exception + fname::String + lineno::Int + func::String + ecode::Int + message::String +end + +# COV_EXCL_START +function Base.showerror(io::IO, e::CeedError) + println(io, "libCEED error code ", e.ecode, " in ", e.func) + println(io, e.fname, ':', e.lineno, '\n') + println(io, e.message) +end +# COV_EXCL_STOP + +function handle_ceed_error( + ceed::C.Ceed, + c_fname::Cstring, + lineno::Cint, + c_func::Cstring, + ecode::Cint, + c_format::Cstring, + args::Ptr{Cvoid}, +) + c_message = ccall( + (:CeedErrorFormat, C.libceed), + Cstring, + (C.Ceed, Cstring, Ptr{Cvoid}), + ceed, + c_format, + args, + ) + fname = unsafe_string(c_fname) + func = unsafe_string(c_func) + message = unsafe_string(c_message) + throw(CeedError(fname, lineno, func, ecode, message)) +end + +mutable struct Ceed + ref::RefValue{C.Ceed} +end + +""" + Ceed(spec="/cpu/self") + +Wraps a libCEED `Ceed` object, created with the given resource specification string. +""" +function Ceed(spec::AbstractString="/cpu/self") + obj = Ceed(Ref{C.Ceed}()) + C.CeedInit(spec, obj.ref) + ehandler = @cfunction( + handle_ceed_error, + Cint, + (C.Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}) + ) + C.CeedSetErrorHandler(obj.ref[], ehandler) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + destroy(x) + end + return obj +end +destroy(c::Ceed) = C.CeedDestroy(c.ref) # COV_EXCL_LINE +Base.getindex(c::Ceed) = c.ref[] + +Base.show(io::IO, ::MIME"text/plain", c::Ceed) = ceed_show(io, c, C.CeedView) + +""" + getresource(c::Ceed) + +Returns the resource string associated with the given [`Ceed`](@ref) object. +""" +function getresource(c::Ceed) + res = Ref{Cstring}() + C.CeedGetResource(c[], res) + unsafe_string(res[]) +end + +""" + isdeterministic(c::Ceed) + +Returns true if backend of the given [`Ceed`](@ref) object is deterministic, and false +otherwise. +""" +function isdeterministic(c::Ceed) + isdet = Ref{Bool}() + C.CeedIsDeterministic(c[], isdet) + isdet[] +end + +""" + get_preferred_memtype(c::Ceed) + +Returns the preferred [`MemType`](@ref) (either `MEM_HOST` or `MEM_DEVICE`) of the given +[`Ceed`](@ref) object. +""" +function get_preferred_memtype(c::Ceed) + mtype = Ref{MemType}() + C.CeedGetPreferredMemType(c[], mtype) + mtype[] +end + +""" + iscuda(c::Ceed) + +Returns true if the given [`Ceed`](@ref) object has resource `"/gpu/cuda/*"` and false +otherwise. +""" +function iscuda(c::Ceed) + res_split = split(getresource(c), "/") + length(res_split) >= 3 && res_split[3] == "cuda" +end diff --git a/julia/LibCEED.jl/src/CeedVector.jl b/julia/LibCEED.jl/src/CeedVector.jl new file mode 100644 index 0000000000..7d4dc0304a --- /dev/null +++ b/julia/LibCEED.jl/src/CeedVector.jl @@ -0,0 +1,308 @@ +import LinearAlgebra: norm + +abstract type AbstractCeedVector end + +struct CeedVectorActive <: AbstractCeedVector end +Base.getindex(::CeedVectorActive) = C.CEED_VECTOR_ACTIVE[] + +struct CeedVectorNone <: AbstractCeedVector end +Base.getindex(::CeedVectorNone) = C.CEED_VECTOR_NONE[] + +mutable struct CeedVector <: AbstractCeedVector + ref::RefValue{C.CeedVector} + arr::Union{Nothing,AbstractArray} + CeedVector(ref::Ref{C.CeedVector}) = new(ref, nothing) +end + +""" + CeedVector(c::Ceed, len::Integer) + +Creates a `CeedVector` of given length. +""" +function CeedVector(c::Ceed, len::Integer) + ref = Ref{C.CeedVector}() + C.CeedVectorCreate(c[], len, ref) + obj = CeedVector(ref) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + destroy(x) + end + return obj +end +destroy(v::CeedVector) = C.CeedVectorDestroy(v.ref) # COV_EXCL_LINE +Base.getindex(v::CeedVector) = v.ref[] + +Base.summary(io::IO, v::CeedVector) = print(io, length(v), "-element CeedVector") +function Base.show(io::IO, ::MIME"text/plain", v::CeedVector) + summary(io, v) + println(io, ":") + witharray_read(v, MEM_HOST) do arr + Base.print_array(io, arr) + end +end +Base.show(io::IO, v::CeedVector) = witharray_read(a -> show(io, a), v, MEM_HOST) + +function Base.length(::Type{T}, v::CeedVector) where {T} + len = Ref{C.CeedInt}() + C.CeedVectorGetLength(v[], len) + return T(len[]) +end + +Base.ndims(::CeedVector) = 1 +Base.ndims(::Type{CeedVector}) = 1 +Base.axes(v::CeedVector) = (Base.OneTo(length(v)),) +Base.size(v::CeedVector) = (length(Int, v),) +Base.length(v::CeedVector) = length(Int, v) + +""" + setvalue!(v::CeedVector, val::CeedScalar) + +Set the [`CeedVector`](@ref) to a constant value. +""" +setvalue!(v::CeedVector, val::CeedScalar) = C.CeedVectorSetValue(v[], val) +""" + setindex!(v::CeedVector, val::CeedScalar) + v[] = val + +Set the [`CeedVector`](@ref) to a constant value, synonymous to [`setvalue!`](@ref). +""" +Base.setindex!(v::CeedVector, val::CeedScalar) = setvalue!(v, val) + +""" + norm(v::CeedVector, ntype::NormType) + +Return the norm of the given [`CeedVector`](@ref). + +The norm type can either be specified as one of `NORM_1`, `NORM_2`, `NORM_MAX`. +""" +function norm(v::CeedVector, ntype::NormType) + nrm = Ref{CeedScalar}() + C.CeedVectorNorm(v[], ntype, nrm) + nrm[] +end + +""" + norm(v::CeedVector, p::Real) + +Return the norm of the given [`CeedVector`](@ref), see [`norm(::CeedVector, +::NormType)`](@ref). + +`p` can have value 1, 2, or Inf, corresponding to `NORM_1`, `NORM_2`, and `NORM_MAX`, +respectively. +""" +function norm(v::CeedVector, p::Real) + if p == 1 + ntype = NORM_1 + elseif p == 2 + ntype = NORM_2 + elseif isinf(p) + ntype = NORM_MAX + else + error("norm(v::CeedVector, p): p must be 1, 2, or Inf") + end + norm(v, ntype) +end + +""" + reciprocal!(v::CeedVector) + +Set `v` to be equal to its elementwise reciprocal. +""" +reciprocal!(v::CeedVector) = C.CeedVectorReciprocal(v[]) + +""" + setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr) + +Set the array used by a [`CeedVector`](@ref), freeing any previously allocated array if +applicable. The backend may copy values to a different [`MemType`](@ref). See also +[`syncarray!`](@ref) and [`takearray!`](@ref). + +!!! warning "Avoid OWN_POINTER CopyMode" + The [`CopyMode`](@ref) `OWN_POINTER` is not suitable for use with arrays that are + allocated by Julia, since those cannot be properly freed from libCEED. +""" +function setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr) + C.CeedVectorSetArray(v[], mtype, cmode, arr) + if cmode == USE_POINTER + v.arr = arr + end +end + +""" + syncarray!(v::CeedVector, mtype::MemType) + +Sync the [`CeedVector`](@ref) to a specified [`MemType`](@ref). This function is used to +force synchronization of arrays set with [`setarray!`](@ref). If the requested memtype is +already synchronized, this function results in a no-op. +""" +syncarray!(v::CeedVector, mtype::MemType) = C.CeedVectorSyncArray(v[], mtype) + +""" + takearray!(v::CeedVector, mtype::MemType) + +Take ownership of the [`CeedVector`](@ref) array and remove the array from the +[`CeedVector`](@ref). The caller is responsible for managing and freeing the array. The +array is returns as a `Ptr{CeedScalar}`. +""" +function takearray!(v::CeedVector, mtype::MemType) + ptr = Ref{Ptr{CeedScalar}}() + C.CeedVectorTakeArray(v[], mtype, ptr) + v.arr = nothing + ptr[] +end + +# Helper function to parse arguments of @witharray and @witharray_read +function witharray_parse(assignment, args) + if !Meta.isexpr(assignment, :(=)) + error("@witharray must have first argument of the form v_arr=v") # COV_EXCL_LINE + end + arr = assignment.args[1] + v = assignment.args[2] + mtype = MEM_HOST + sz = :((length($(esc(v))),)) + body = args[end] + for i = 1:length(args)-1 + a = args[i] + if !Meta.isexpr(a, :(=)) + error("Incorrect call to @witharray or @witharray_read") # COV_EXCL_LINE + end + if a.args[1] == :mtype + mtype = a.args[2] + elseif a.args[1] == :size + sz = esc(a.args[2]) + end + end + arr, v, sz, mtype, body +end + +""" + @witharray(v_arr=v, [size=(dims...)], [mtype=MEM_HOST], body) + +Executes `body`, having extracted the contents of the [`CeedVector`](@ref) `v` as an array +with name `v_arr`. If the [`memory type`](@ref MemType) `mtype` is not provided, `MEM_HOST` +will be used. If the size is not specified, a flat vector will be assumed. + +# Examples +Negate the contents of `CeedVector` `v`: +``` +@witharray v_arr=v v_arr *= -1.0 +``` +""" +macro witharray(assignment, args...) + arr, v, sz, mtype, body = witharray_parse(assignment, args) + quote + arr_ref = Ref{Ptr{C.CeedScalar}}() + C.CeedVectorGetArray($(esc(v))[], $(esc(mtype)), arr_ref) + try + $(esc(arr)) = UnsafeArray(arr_ref[], Int.($sz)) + $(esc(body)) + finally + C.CeedVectorRestoreArray($(esc(v))[], arr_ref) + end + end +end + +""" + @witharray_read(v_arr=v, [size=(dims...)], [mtype=MEM_HOST], body) + +Same as [`@witharray`](@ref), but provides read-only access to the data. +""" +macro witharray_read(assignment, args...) + arr, v, sz, mtype, body = witharray_parse(assignment, args) + quote + arr_ref = Ref{Ptr{C.CeedScalar}}() + C.CeedVectorGetArrayRead($(esc(v))[], $(esc(mtype)), arr_ref) + try + $(esc(arr)) = UnsafeArray(arr_ref[], Int.($sz)) + $(esc(body)) + finally + C.CeedVectorRestoreArrayRead($(esc(v))[], arr_ref) + end + end +end + +""" + setindex!(v::CeedVector, v2::AbstractArray) + v[] = v2 + +Sets the values of [`CeedVector`](@ref) `v` equal to those of `v2` using broadcasting. +""" +Base.setindex!(v::CeedVector, v2::AbstractArray) = @witharray(a = v, a .= v2) + +""" + CeedVector(c::Ceed, v2::AbstractVector; mtype=MEM_HOST, cmode=COPY_VALUES) + +Creates a new [`CeedVector`](@ref) using the contents of the given vector `v2`. By default, +the contents of `v2` will be copied to the new [`CeedVector`](@ref), but this behavior can +be changed by specifying a different `cmode`. +""" +function CeedVector(c::Ceed, v2::AbstractVector; mtype=MEM_HOST, cmode=COPY_VALUES) + v = CeedVector(c, length(v2)) + setarray!(v, mtype, cmode, v2) + v +end + +""" + Vector(v::CeedVector) + +Create a new `Vector` by copying the contents of `v`. +""" +function Base.Vector(v::CeedVector) + v2 = Vector{CeedScalar}(undef, length(v)) + @witharray_read(a = v, v2 .= a) +end + +""" + witharray(f, v::CeedVector, mtype=MEM_HOST) + +Calls `f` with an array containing the data of the `CeedVector` `v`, using [`memory +type`](@ref MemType) `mtype`. + +Because of performance issues involving closures, if `f` is a complex operation, it may be +more efficient to use the macro version `@witharray` (cf. the section on "Performance of +captured variable" in the [Julia +documentation](https://docs.julialang.org/en/v1/manual/performance-tips) and related [GitHub +issue](https://github.com/JuliaLang/julia/issues/15276). + +# Examples + +Return the sum of a vector: +``` +witharray(sum, v) +``` +""" +function witharray(f, v::CeedVector, mtype::MemType=MEM_HOST) + arr_ref = Ref{Ptr{C.CeedScalar}}() + C.CeedVectorGetArray(v[], mtype, arr_ref) + arr = UnsafeArray(arr_ref[], (length(v),)) + res = try + f(arr) + finally + C.CeedVectorRestoreArray(v[], arr_ref) + end + return res +end + +""" + witharray_read(f, v::CeedVector, mtype::MemType=MEM_HOST) + +Same as [`witharray`](@ref), but with read-only access to the data. + +# Examples + +Display the contents of a vector: +``` +witharray_read(display, v) +``` +""" +function witharray_read(f, v::CeedVector, mtype::MemType=MEM_HOST) + arr_ref = Ref{Ptr{C.CeedScalar}}() + C.CeedVectorGetArrayRead(v[], mtype, arr_ref) + arr = UnsafeArray(arr_ref[], (length(v),)) + res = try + f(arr) + finally + C.CeedVectorRestoreArrayRead(v[], arr_ref) + end + return res +end diff --git a/julia/LibCEED.jl/src/Context.jl b/julia/LibCEED.jl/src/Context.jl new file mode 100644 index 0000000000..9e3c27a502 --- /dev/null +++ b/julia/LibCEED.jl/src/Context.jl @@ -0,0 +1,56 @@ +mutable struct Context + ref::RefValue{C.CeedQFunctionContext} + data::Any + function Context(ref::Ref{C.CeedQFunctionContext}) + obj = new(ref) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + C.CeedQFunctionContextDestroy(x.ref) + end + return obj + end +end +Base.getindex(ctx::Context) = ctx.ref[] +Base.show(io::IO, ::MIME"text/plain", c::Context) = + ceed_show(io, c, C.CeedQFunctionContextView) + +""" + Context(ceed::Ceed, data; mtype=MEM_HOST, cmode=USE_POINTER) + +Create a `CeedQFunctionContext` object that allows user Q-functions to access an arbitrary +data object. `data` should be an instance of a mutable struct. If the copy mode `cmode` is +`USE_POINTER`, then the data will be preserved from the GC when assigned to a `QFunction` +object using `set_context!`. + +Copy mode `OWN_POINTER` is not supported by this interface because Julia-allocated objects +cannot be freed from C. +""" +function Context(c::Ceed, data; mtype=MEM_HOST, cmode=USE_POINTER) + ref = Ref{C.CeedQFunctionContext}() + C.CeedQFunctionContextCreate(c[], ref) + ctx = Context(ref) + set_data!(ctx, mtype, cmode, data) + return ctx +end + +function set_data!(ctx::Context, mtype, cmode::CopyMode, data) + # Store a reference to the context data so that it will not be GC'd before + # it is accessed in the user Q-function. + # A reference to the context object is stored in the QFunction object, and + # references to the QFunctions are stored in the Operator. + # This means that when `apply!(op, ...)` is called, the context data is + # ensured to be valid. + if cmode == USE_POINTER + ctx.data = data + elseif cmode == OWN_POINTER + error("set_data!: copy mode OWN_POINTER is not supported") + end + + C.CeedQFunctionContextSetData( + ctx[], + mtype, + cmode, + sizeof(data), + pointer_from_objref(data), + ) +end diff --git a/julia/LibCEED.jl/src/Cuda.jl b/julia/LibCEED.jl/src/Cuda.jl new file mode 100644 index 0000000000..372871adf5 --- /dev/null +++ b/julia/LibCEED.jl/src/Cuda.jl @@ -0,0 +1,155 @@ +# COV_EXCL_START +using .CUDA, Cassette + +cuda_is_loaded = true + +#! format: off +const cudafuns = ( + :cos, :cospi, :sin, :sinpi, :tan, + :acos, :asin, :atan, + :cosh, :sinh, :tanh, + :acosh, :asinh, :atanh, + :log, :log10, :log1p, :log2, + :exp, :exp2, :exp10, :expm1, :ldexp, + :abs, + :sqrt, :cbrt, + :ceil, :floor, +) +#! format: on + +Cassette.@context CeedCudaContext + +@inline function Cassette.overdub(::CeedCudaContext, ::typeof(Core.kwfunc), f) + return Core.kwfunc(f) +end +@inline function Cassette.overdub(::CeedCudaContext, ::typeof(Core.apply_type), args...) + return Core.apply_type(args...) +end +@inline function Cassette.overdub( + ::CeedCudaContext, + ::typeof(StaticArrays.Size), + x::Type{<:AbstractArray{<:Any,N}}, +) where {N} + return StaticArrays.Size(x) +end + +for f in cudafuns + @eval @inline function Cassette.overdub( + ::CeedCudaContext, + ::typeof(Base.$f), + x::Union{Float32,Float64}, + ) + return CUDA.$f(x) + end +end + +function setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr::CuArray) + ptr = Ptr{CeedScalar}(UInt64(pointer(arr))) + C.CeedVectorSetArray(v[], mtype, cmode, ptr) + if cmode == USE_POINTER + v.arr = arr + end +end + +struct FieldsCuda + inputs::NTuple{16,Int} + outputs::NTuple{16,Int} +end + +function generate_kernel(qf_name, kf, dims_in, dims_out) + ninputs = length(dims_in) + noutputs = length(dims_out) + + input_sz = prod.(dims_in) + output_sz = prod.(dims_out) + + f_ins = [Symbol("rqi$i") for i = 1:ninputs] + f_outs = [Symbol("rqo$i") for i = 1:noutputs] + + args = Vector{Union{Symbol,Expr}}(undef, ninputs + noutputs) + def_ins = Vector{Expr}(undef, ninputs) + f_ins_j = Vector{Union{Symbol,Expr}}(undef, ninputs) + for i = 1:ninputs + if length(dims_in[i]) == 0 + def_ins[i] = :(local $(f_ins[i])) + f_ins_j[i] = f_ins[i] + args[i] = f_ins[i] + else + def_ins[i] = + :($(f_ins[i]) = LibCEED.MArray{Tuple{$(dims_in[i]...)},Float64}(undef)) + f_ins_j[i] = :($(f_ins[i])[j]) + args[i] = :(LibCEED.SArray{Tuple{$(dims_in[i]...)},Float64}($(f_ins[i]))) + end + end + for i = 1:noutputs + args[ninputs+i] = f_outs[i] + end + + def_outs = [ + :($(f_outs[i]) = LibCEED.MArray{Tuple{$(dims_out[i]...)},Float64}(undef)) + for i = 1:noutputs + ] + + device_ptr_type = Core.LLVMPtr{CeedScalar,LibCEED.AS.Global} + + read_quads_in = [ + :( + for j = 1:$(input_sz[i]) + $(f_ins_j[i]) = unsafe_load( + reinterpret($device_ptr_type, fields.inputs[$i]), + q + (j - 1)*Q, + a, + ) + end + ) for i = 1:ninputs + ] + + write_quads_out = [ + :( + for j = 1:$(output_sz[i]) + unsafe_store!( + reinterpret($device_ptr_type, fields.outputs[$i]), + $(f_outs[i])[j], + q + (j - 1)*Q, + a, + ) + end + ) for i = 1:noutputs + ] + + qf = gensym(qf_name) + quote + function $qf(ctx_ptr, Q, fields) + gd = LibCEED.gridDim() + bi = LibCEED.blockIdx() + bd = LibCEED.blockDim() + ti = LibCEED.threadIdx() + + inc = bd.x*gd.x + + $(def_ins...) + $(def_outs...) + + # Alignment for data read/write + a = Val($(Base.datatype_alignment(CeedScalar))) + + # Cassette context for replacing intrinsics with CUDA versions + ctx = LibCEED.CeedCudaContext() + + for q = (ti.x+(bi.x-1)*bd.x):inc:Q + $(read_quads_in...) + LibCEED.Cassette.overdub(ctx, $kf, ctx_ptr, $(args...)) + $(write_quads_out...) + end + return + end + end +end + +function mk_cufunction(ceed, def_module, qf_name, kf, dims_in, dims_out) + k_fn = Core.eval(def_module, generate_kernel(qf_name, kf, dims_in, dims_out)) + tt = Tuple{Ptr{Nothing},Int32,FieldsCuda} + host_k = cufunction(k_fn, tt; maxregs=64) + return host_k.fun.handle +end +# COV_EXCL_STOP diff --git a/julia/LibCEED.jl/src/ElemRestriction.jl b/julia/LibCEED.jl/src/ElemRestriction.jl new file mode 100644 index 0000000000..8c4ab94c3f --- /dev/null +++ b/julia/LibCEED.jl/src/ElemRestriction.jl @@ -0,0 +1,288 @@ +abstract type AbstractElemRestriction end + +""" + ElemRestrictionNone() + +Returns the singleton object corresponding to libCEED's `CEED_ELEMRESTRICTION_NONE` +""" +struct ElemRestrictionNone <: AbstractElemRestriction end +Base.getindex(::ElemRestrictionNone) = C.CEED_ELEMRESTRICTION_NONE[] + +""" + ElemRestriction + +Wraps a `CeedElemRestriction` object, representing the restriction from local vectors to +elements. An `ElemRestriction` object can be created using [`create_elem_restriction`](@ref) +or [`create_elem_restriction_strided`](@ref). +""" +mutable struct ElemRestriction <: AbstractElemRestriction + ref::RefValue{C.CeedElemRestriction} + function ElemRestriction(ref) + obj = new(ref) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + destroy(x) + end + return obj + end +end +destroy(r::ElemRestriction) = C.CeedElemRestrictionDestroy(r.ref) # COV_EXCL_LINE +Base.getindex(r::ElemRestriction) = r.ref[] +Base.show(io::IO, ::MIME"text/plain", e::ElemRestriction) = + ceed_show(io, e, C.CeedElemRestrictionView) + +@doc raw""" + create_elem_restriction( + ceed::Ceed, + nelem, + elemsize, + ncomp, + compstride, + lsize, + offsets::AbstractArray{CeedInt}, + mtype::MemType=MEM_HOST, + cmode::CopyMode=COPY_VALUES, + ) + +Create a `CeedElemRestriction`. + +!!! warning "Zero-based indexing" + In the below notation, we are using **0-based indexing**. libCEED expects the offset + indices to be 0-based. + +# Arguments: +- `ceed`: The [`Ceed`](@ref) object +- `nelem`: Number of elements described in the `offsets` array +- `elemsize`: Size (number of "nodes") per element +- `ncomp`: Number of field components per interpolation node (1 for scalar fields) +- `compstride`: Stride between components for the same L-vector "node". Data for node $i$, + component $j$, element $k$ can be found in the L-vector at index `offsets[i + + k*elemsize] + j*compstride`. +- `lsize`: The size of the L-vector. This vector may be larger than the elements and + fields given by this restriction. +- `offsets`: Array of shape `(elemsize, nelem)`. Column $i$ holds the ordered list of the + offsets (into the input [`CeedVector`](@ref)) for the unknowns corresponding + to element $i$, where $0 \leq i < \textit{nelem}$. All offsets must be in + the range $[0, \textit{lsize} - 1]$. +- `mtype`: Memory type of the `offsets` array, see [`MemType`](@ref) +- `cmode`: Copy mode for the `offsets` array, see [`CopyMode`](@ref) +""" +function create_elem_restriction( + c::Ceed, + nelem, + elemsize, + ncomp, + compstride, + lsize, + offsets::AbstractArray{CeedInt}; + mtype::MemType=MEM_HOST, + cmode::CopyMode=COPY_VALUES, +) + ref = Ref{C.CeedElemRestriction}() + C.CeedElemRestrictionCreate( + c[], + nelem, + elemsize, + ncomp, + compstride, + lsize, + mtype, + cmode, + offsets, + ref, + ) + ElemRestriction(ref) +end + +@doc raw""" + create_elem_restriction_strided(ceed::Ceed, nelem, elemsize, ncomp, lsize, strides) + +Create a strided `CeedElemRestriction`. + +!!! warning "Zero-based indexing" + In the below notation, we are using **0-based indexing**. libCEED expects the offset + indices to be 0-based. + +# Arguments: +- `ceed`: The [`Ceed`](@ref) object +- `nelem`: Number of elements described by the restriction +- `elemsize`: Size (number of "nodes") per element +- `ncomp`: Number of field components per interpolation node (1 for scalar fields) +- `lsize`: The size of the L-vector. This vector may be larger than the elements and + fields given by this restriction. +- `strides`: Array for strides between [nodes, components, elements]. Data for node $i$, + component $j$, element $k$ can be found in the L-vector at index `i*strides[0] + + j*strides[1] + k*strides[2]`. [`STRIDES_BACKEND`](@ref) may be used with + vectors created by a Ceed backend. +""" +function create_elem_restriction_strided(c::Ceed, nelem, elemsize, ncomp, lsize, strides) + ref = Ref{C.CeedElemRestriction}() + C.CeedElemRestrictionCreateStrided(c[], nelem, elemsize, ncomp, lsize, strides, ref) + ElemRestriction(ref) +end + +""" + apply!( + r::ElemRestriction, + u::CeedVector, + ru::CeedVector; + tmode=NOTRANSPOSE, + request=RequestImmediate(), + ) + +Use the [`ElemRestriction`](@ref) to convert from L-vector to an E-vector (or apply the +tranpose operation). The input [`CeedVector`](@ref) is `u` and the result stored in `ru`. + +If `tmode` is `TRANSPOSE`, then the result is added to `ru`. If `tmode` is `NOTRANSPOSE`, +then `ru` is overwritten with the result. +""" +function apply!( + r::ElemRestriction, + u::CeedVector, + ru::CeedVector; + tmode=NOTRANSPOSE, + request=RequestImmediate(), +) + C.CeedElemRestrictionApply(r[], tmode, u[], ru[], request[]) +end + +""" + apply(r::ElemRestriction, u::AbstractVector; tmode=NOTRANSPOSE) + +Use the [`ElemRestriction`](@ref) to convert from L-vector to an E-vector (or apply the +tranpose operation). The input is given by `u`, and the result is returned as an array of +type `Vector{CeedScalar}`. +""" +function apply(r::ElemRestriction, u::AbstractVector; tmode=NOTRANSPOSE) + ceed_ref = Ref{C.Ceed}() + ccall( + (:CeedElemRestrictionGetCeed, C.libceed), + Cint, + (C.CeedElemRestriction, Ptr{C.Ceed}), + r[], + ceed_ref, + ) + c = Ceed(ceed_ref) + uv = CeedVector(c, u) + if tmode == NOTRANSPOSE + ruv = create_evector(r) + else + ruv = create_lvector(r) + ruv[] = 0.0 + end + apply!(r, uv, ruv; tmode=tmode) + Vector(ruv) +end + +""" + create_evector(r::ElemRestriction) + +Return a new [`CeedVector`](@ref) E-vector. +""" +function create_evector(r::ElemRestriction) + ref = Ref{C.CeedVector}() + C.CeedElemRestrictionCreateVector(r[], C_NULL, ref) + CeedVector(ref) +end + +""" + create_lvector(r::ElemRestriction) + +Return a new [`CeedVector`](@ref) L-vector. +""" +function create_lvector(r::ElemRestriction) + ref = Ref{C.CeedVector}() + C.CeedElemRestrictionCreateVector(r[], ref, C_NULL) + CeedVector(ref) +end + +""" + create_vectors(r::ElemRestriction) + +Return an (L-vector, E-vector) pair. +""" +function create_vectors(r::ElemRestriction) + l_ref = Ref{C.CeedVector}() + e_ref = Ref{C.CeedVector}() + C.CeedElemRestrictionCreateVector(r[], l_ref, e_ref) + CeedVector(l_ref), CeedVector(e_ref) +end + +""" + getcompstride(r::ElemRestriction) + +Get the L-vector component stride. +""" +function getcompstride(r::ElemRestriction) + lsize = Ref{CeedInt}() + C.CeedElemRestrictionGetCompStride(r[], lsize) + lsize[] +end + +""" + getnumelements(r::ElemRestriction) + +Get the total number of elements in the range of an [`ElemRestriction`](@ref). +""" +function getnumelements(r::ElemRestriction) + result = Ref{CeedInt}() + C.CeedElemRestrictionGetNumElements(r[], result) + result[] +end + +""" + getelementsize(r::ElemRestriction) + +Get the size of elements in the given [`ElemRestriction`](@ref). +""" +function getelementsize(r::ElemRestriction) + result = Ref{CeedInt}() + C.CeedElemRestrictionGetElementSize(r[], result) + result[] +end + +""" + getlvectorsize(r::ElemRestriction) + +Get the size of an L-vector for the given [`ElemRestriction`](@ref). +""" +function getlvectorsize(r::ElemRestriction) + result = Ref{CeedInt}() + C.CeedElemRestrictionGetLVectorSize(r[], result) + result[] +end + +""" + getnumcomponents(r::ElemRestriction) + +Get the number of components in the elements of an [`ElemRestriction`](@ref). +""" +function getnumcomponents(r::ElemRestriction) + result = Ref{CeedInt}() + C.CeedElemRestrictionGetNumComponents(r[], result) + result[] +end + +""" + getmultiplicity!(r::ElemRestriction, v::AbstractCeedVector) + +Get the multiplicity of nodes in an [`ElemRestriction`](@ref). The [`CeedVector`](@ref) `v` +should be an L-vector (i.e. `length(v) == getlvectorsize(r)`, see [`create_lvector`](@ref)). +""" +function getmultiplicity!(r::ElemRestriction, v::AbstractCeedVector) + @assert length(v) == getlvectorsize(r) + C.CeedElemRestrictionGetMultiplicity(r[], v[]) +end + +""" + getmultiplicity(r::ElemRestriction) + +Convenience function to get the multiplicity of nodes in the [`ElemRestriction`](@ref), +where the result is returned in a newly allocated Julia `Vector{CeedScalar}` (see also +[`getmultiplicity!`](@ref)). +""" +function getmultiplicity(r::ElemRestriction) + v = create_lvector(r) + getmultiplicity!(r, v) + Vector(v) +end diff --git a/julia/LibCEED.jl/src/Globals.jl b/julia/LibCEED.jl/src/Globals.jl new file mode 100644 index 0000000000..41eae89c93 --- /dev/null +++ b/julia/LibCEED.jl/src/Globals.jl @@ -0,0 +1,117 @@ +""" + CeedScalar + +Scalar (floating point) type. Equivalent to `Float64`. +""" +const CeedScalar = C.CeedScalar +""" + CeedInt + +Integer type, used for indexing. Equivalent to `Int32`. +""" +const CeedInt = C.CeedInt + +""" + QuadMode + +One of `GAUSS` or `GAUSS_LOBATTO`. +""" +const QuadMode = C.CeedQuadMode +const GAUSS = C.CEED_GAUSS +const GAUSS_LOBATTO = C.CEED_GAUSS_LOBATTO + +""" + MemType + +One of `MEM_HOST` or `MEM_DEVICE`. +""" +const MemType = C.CeedMemType +const MEM_HOST = C.CEED_MEM_HOST +const MEM_DEVICE = C.CEED_MEM_DEVICE + +""" + CopyMode + +One of `COPY_VALUES`, `USE_POINTER` or `OWN_POINTER`. + +`OWN_POINTER` is not typically supported for objects created in Julia, because those must be +destroyed by the garbage collector, and cannot be freed from C. +""" +const CopyMode = C.CeedCopyMode +const COPY_VALUES = C.CEED_COPY_VALUES +const USE_POINTER = C.CEED_USE_POINTER +const OWN_POINTER = C.CEED_OWN_POINTER + +""" + EvalMode + +Evaluation mode used in the specification of input and output fields for Q-functions, e.g. +in [`@interior_qf`](@ref). + +One of: +- `EVAL_NONE` +- `EVAL_INTERP` +- `EVAL_GRAD` +- `EVAL_DIV` +- `EVAL_CURL` +- `EVAL_WEIGHT` + +See the [libCEED +documentation](https://libceed.readthedocs.io/en/latest/api/CeedBasis/?highlight=EVAL_MODE#c.CeedEvalMode) +for further information. +""" +const EvalMode = C.CeedEvalMode +const EVAL_NONE = C.CEED_EVAL_NONE +const EVAL_INTERP = C.CEED_EVAL_INTERP +const EVAL_GRAD = C.CEED_EVAL_GRAD +const EVAL_DIV = C.CEED_EVAL_DIV +const EVAL_CURL = C.CEED_EVAL_CURL +const EVAL_WEIGHT = C.CEED_EVAL_WEIGHT + +""" + NormType + +Denotes type of vector norm. One of `NORM_1`, `NORM_2`, or `NORM_MAX`. +""" +const NormType = C.CeedNormType +const NORM_1 = C.CEED_NORM_1 +const NORM_2 = C.CEED_NORM_2 +const NORM_MAX = C.CEED_NORM_MAX + +""" + TransposeMose + +Denotes whether a linear transformation or its transpose should be applied. Either +`NOTRANSPOSE` or `TRANSPOSE`. +""" +const TransposeMode = C.CeedTransposeMode +const NOTRANSPOSE = C.CEED_NOTRANSPOSE +const TRANSPOSE = C.CEED_TRANSPOSE + +""" + Topology + +Type of basis shape to create non-tensor H1 element basis. One of `LINE`, `TRIANGLE`, +`QUAD`, `TET`, `PYRAMID`, `PRISM`, or `HEX`. + +The dimension can be extracted with bitshift: + + dim = Int(topology) >> 16 +""" +const Topology = C.CeedElemTopology +const LINE = C.CEED_LINE +const TRIANGLE = C.CEED_TRIANGLE +const QUAD = C.CEED_QUAD +const TET = C.CEED_TET +const PYRAMID = C.CEED_PYRAMID +const PRISM = C.CEED_PRISM +const HEX = C.CEED_HEX + +function set_globals() + @doc """ + STRIDES_BACKEND + + Indicate that the stride is determined by the backend. + """ + global STRIDES_BACKEND = C.CEED_STRIDES_BACKEND[] +end diff --git a/julia/LibCEED.jl/src/LibCEED.jl b/julia/LibCEED.jl/src/LibCEED.jl new file mode 100644 index 0000000000..6202d7c0e8 --- /dev/null +++ b/julia/LibCEED.jl/src/LibCEED.jl @@ -0,0 +1,150 @@ +module LibCEED + +using StaticArrays, UnsafeArrays, Requires +using Base: RefValue + +# import low-level C interface +include("C.jl") +import .C + +export @interior_qf, + @witharray, + @witharray_read, + Abscissa, + AbscissaAndWeights, + Basis, + BasisCollocated, + COPY_VALUES, + Ceed, + CeedDim, + CeedInt, + CeedScalar, + CeedVector, + CeedVectorActive, + CeedVectorNone, + Context, + CopyMode, + EVAL_CURL, + EVAL_DIV, + EVAL_GRAD, + EVAL_INTERP, + EVAL_NONE, + EVAL_WEIGHT, + ElemRestriction, + ElemRestrictionNone, + EvalMode, + GAUSS, + GAUSS_LOBATTO, + HEX, + LINE, + MEM_DEVICE, + MEM_HOST, + MemType, + NORM_1, + NORM_2, + NORM_MAX, + NOTRANSPOSE, + NormType, + OWN_POINTER, + Operator, + PRISM, + PYRAMID, + QFunction, + QFunctionNone, + QUAD, + QuadMode, + RequestImmediate, + RequestOrdered, + STRIDES_BACKEND, + TET, + TRANSPOSE, + TRIANGLE, + Topology, + TransposeMode, + USE_POINTER, + UserQFunction, + add_input!, + add_output!, + apply!, + apply_add!, + apply, + assemble, + assemble_add_diagonal!, + assemble_diagonal!, + create_composite_operator, + create_elem_restriction, + create_elem_restriction_strided, + create_evector, + create_h1_basis, + create_identity_qfunction, + create_interior_qfunction, + create_lvector, + create_tensor_h1_basis, + create_tensor_h1_lagrange_basis, + create_vectors, + det, + extract_array, + extract_context, + gauss_quadrature, + get_preferred_memtype, + getcompstride, + getnumelements, + getelementsize, + getlvectorsize, + getmultiplicity!, + getmultiplicity, + getdimension, + getgrad, + getgrad1d, + getinterp, + getinterp1d, + getnumcomponents, + getnumnodes, + getnumnodes1d, + getnumqpts, + getnumqpts1d, + getqref, + getqweights, + getresource, + gettopology, + getvoigt!, + getvoigt, + iscuda, + isdeterministic, + lobatto_quadrature, + norm, + reciprocal!, + set_context!, + set_cufunction!, + set_data!, + set_field!, + setarray!, + setvalue!, + setvoigt!, + setvoigt, + syncarray!, + takearray!, + witharray, + witharray_read + +include("Globals.jl") +include("Ceed.jl") +include("CeedVector.jl") +include("Basis.jl") +include("ElemRestriction.jl") +include("Quadrature.jl") +include("Context.jl") +include("UserQFunction.jl") +include("QFunction.jl") +include("Request.jl") +include("Operator.jl") +include("Misc.jl") + +cuda_is_loaded = false + +function __init__() + @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("Cuda.jl") + set_globals() +end + +end # module diff --git a/julia/LibCEED.jl/src/Misc.jl b/julia/LibCEED.jl/src/Misc.jl new file mode 100644 index 0000000000..1ace1ad1b2 --- /dev/null +++ b/julia/LibCEED.jl/src/Misc.jl @@ -0,0 +1,119 @@ +import LinearAlgebra: det + +""" + CeedDim(dim) + +The singleton object of type `CeedDim{dim}`, used for dispatch to linear algebra operations +specialized for small matrices (1, 2, or 3 dimensions). +""" +struct CeedDim{dim} end +@inline CeedDim(dim) = CeedDim{Int(dim)}() + +""" + det(J, ::CeedDim{dim}) + +Specialized determinant calculations for matrices of size 1, 2, or 3. +""" +@inline det(J, ::CeedDim{1}) = @inbounds J[1] +@inline det(J, ::CeedDim{2}) = @inbounds J[1]*J[4] - J[3]*J[2] +#! format: off +@inline det(J, ::CeedDim{3}) = @inbounds ( + J[1]*(J[5]*J[9] - J[6]*J[8]) - + J[2]*(J[4]*J[9] - J[6]*J[7]) + + J[3]*(J[4]*J[8] - J[5]*J[7]) +) +#! format: on + +""" + setvoigt(J::StaticArray{Tuple{D,D},T,2}) + setvoigt(J, ::CeedDim{dim}) + +Given a symmetric matrix `J`, return a `SVector` that encodes `J` using the [Voigt +convention](https://en.wikipedia.org/wiki/Voigt_notation). + +The size of the symmetric matrix `J` must be known statically, either specified using +[`CeedDim`](@ref) or `StaticArray`. +""" +@inline setvoigt(J::StaticArray{Tuple{D,D}}) where {D} = setvoigt(J, CeedDim(D)) +@inline setvoigt(J, ::CeedDim{1}) = @inbounds @SVector [J[1]] +@inline setvoigt(J, ::CeedDim{2}) = @inbounds @SVector [J[1], J[4], J[2]] +@inline setvoigt(J, ::CeedDim{3}) = @inbounds @SVector [J[1], J[5], J[9], J[6], J[3], J[2]] + +@inline function setvoigt!(V, J, ::CeedDim{1}) + @inbounds V[1] = J[1] +end + +@inline function setvoigt!(V, J, ::CeedDim{2}) + @inbounds begin + V[1] = J[1] + V[2] = J[4] + V[3] = J[2] + end +end + +@inline function setvoigt!(V, J, ::CeedDim{3}) + @inbounds begin + V[1] = J[1] + V[2] = J[5] + V[3] = J[9] + V[4] = J[6] + V[5] = J[3] + V[6] = J[2] + end +end + +""" + getvoigt(V, ::CeedDim{dim}) + +Given a vector `V` that encodes a symmetric matrix using the [Voigt +convention](https://en.wikipedia.org/wiki/Voigt_notation), return the corresponding +`SMatrix`. +""" +@inline getvoigt(V, ::CeedDim{1}) = @inbounds @SMatrix [V[1]] +@inline getvoigt(V, ::CeedDim{2}) = @inbounds @SMatrix [V[1] V[3]; V[3] V[2]] +@inline getvoigt(V, ::CeedDim{3}) = @inbounds @SMatrix [ + V[1] V[6] V[5] + V[6] V[2] V[4] + V[5] V[4] V[3] +] +@inline getvoigt(V::StaticArray{Tuple{1}}) = getvoigt(V, CeedDim(1)) +@inline getvoigt(V::StaticArray{Tuple{3}}) = getvoigt(V, CeedDim(2)) +@inline getvoigt(V::StaticArray{Tuple{6}}) = getvoigt(V, CeedDim(3)) + +@inline function getvoigt!(J, V, ::CeedDim{1}) + @inbounds J[1, 1] = V[1] +end + +@inline function getvoigt!(J, V, ::CeedDim{2}) + @inbounds begin + #! format: off + J[1,1] = V[1] ; J[1,2] = V[3] + J[2,1] = V[3] ; J[2,2] = V[2] + #! format: on + end +end + +@inline function getvoigt!(J, V, ::CeedDim{3}) + @inbounds begin + #! format: off + J[1,1] = V[1] ; J[1,2] = V[6] ; J[1,3] = V[5] + J[2,1] = V[6] ; J[2,2] = V[2] ; J[2,3] = V[4] + J[3,1] = V[5] ; J[3,2] = V[4] ; J[3,3] = V[3] + #! format: on + end +end + +function tmp_view(obj, view_fn) + str = mktemp() do fname, f + cf = Libc.FILE(f) + er = view_fn(obj, cf.ptr) + ccall(:fflush, Cint, (Ptr{Cvoid},), cf) + seek(f, 0) + read(f, String) + end + chomp(str) +end + +function ceed_show(io::IO, obj, view_fn) + print(io, tmp_view(obj[], view_fn)) +end diff --git a/julia/LibCEED.jl/src/Operator.jl b/julia/LibCEED.jl/src/Operator.jl new file mode 100644 index 0000000000..74fa56297e --- /dev/null +++ b/julia/LibCEED.jl/src/Operator.jl @@ -0,0 +1,149 @@ +mutable struct Operator + ref::RefValue{C.CeedOperator} + qf::AbstractQFunction + dqf::AbstractQFunction + dqfT::AbstractQFunction + sub_ops::Vector{Operator} + function Operator(ref, qf, dqf, dqfT) + obj = new(ref, qf, dqf, dqfT, []) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + destroy(x) + end + return obj + end +end +destroy(op::Operator) = C.CeedOperatorDestroy(op.ref) # COV_EXCL_LINE +Base.getindex(op::Operator) = op.ref[] +Base.show(io::IO, ::MIME"text/plain", op::Operator) = ceed_show(io, op, C.CeedOperatorView) + +""" + Operator(ceed::Ceed; qf, dqf=QFunctionNone(), dqfT=QFunctionNone(), fields) + +Creates a libCEED `CeedOperator` object using the given Q-function `qf`, and optionally its +derivative and derivative transpose. + +An array of fields must be provided, where each element of the array is a tuple containing +the name of the field (as a string or symbol), the corresponding element restriction, basis, +and vector. + +# Examples + +Create the operator that builds the Q-data associated with the mass matrix. +``` +build_oper = Operator( + ceed, + qf=build_qfunc, + fields=[ + (:J, mesh_restr, mesh_basis, CeedVectorActive()), + (:w, ElemRestrictionNone(), mesh_basis, CeedVectorNone()), + (:qdata, sol_restr_i, BasisCollocated(), CeedVectorActive()) + ] +) +``` +""" +function Operator(c::Ceed; qf, dqf=QFunctionNone(), dqfT=QFunctionNone(), fields) + op = Operator(c, qf, dqf, dqfT) + for f ∈ fields + set_field!(op, String(f[1]), f[2], f[3], f[4]) + end + op +end + +function Operator( + c::Ceed, + qf::AbstractQFunction, + dqf::AbstractQFunction, + dqfT::AbstractQFunction, +) + ref = Ref{C.CeedOperator}() + C.CeedOperatorCreate(c[], qf[], dqf[], dqfT[], ref) + Operator(ref, qf, dqf, dqfT) +end + +""" + create_composite_operator(c::Ceed, ops) + +Create an [`Operator`](@ref) whose action represents the sum of the operators in the +collection `ops`. +""" +function create_composite_operator(c::Ceed, ops) + ref = Ref{C.CeedOperator}() + C.CeedCompositeOperatorCreate(c[], ref) + comp_op = Operator(ref, QFunctionNone(), QFunctionNone(), QFunctionNone()) + comp_op.sub_ops = ops + for op ∈ ops + C.CeedCompositeOperatorAddSub(comp_op[], op[]) + end + comp_op +end + +function set_field!( + op::Operator, + fieldname::AbstractString, + r::AbstractElemRestriction, + b::AbstractBasis, + v::AbstractCeedVector, +) + C.CeedOperatorSetField(op[], fieldname, r[], b[], v[]) +end + +""" + apply!(op::Operator, vin, vout; request=RequestImmediate()) + +Apply the action of the operator `op` to the input vector `vin`, and store the result in the +output vector `vout`. + +For non-blocking application, the user can specify a request object. By default, immediate +(synchronous) completion is requested. +""" +function apply!( + op::Operator, + vin::AbstractCeedVector, + vout::AbstractCeedVector; + request=RequestImmediate(), +) + C.CeedOperatorApply(op[], vin[], vout[], request[]) +end + +""" + apply_add!(op::Operator, vin, vout; request=RequestImmediate()) + +Apply the action of the operator `op` to the input vector `vin`, and add the result to the +output vector `vout`. + +For non-blocking application, the user can specify a request object. By default, immediate +(synchronous) completion is requested. +""" +function apply_add!( + op::Operator, + vin::AbstractCeedVector, + vout::AbstractCeedVector; + request=RequestImmediate(), +) + C.CeedOperatorApplyAdd(op[], vin[], vout[], request[]) +end + +""" + assemble_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate()) + +Overwrites a [`CeedVector`](@ref) with the diagonal of a linear [`Operator`](@ref). + +!!! note "Note:" + Currently only [`Operator`](@ref)s with a single field are supported. +""" +function assemble_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate()) + C.CeedOperatorLinearAssembleDiagonal(op[], diag[], request[]) +end + +""" + assemble_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate()) + +Adds the diagonal of a linear [`Operator`](@ref) to the given [`CeedVector`](@ref). + +!!! note "Note:" + Currently only [`Operator`](@ref)s with a single field are supported. +""" +function assemble_add_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate()) + C.CeedOperatorLinearAssembleAddDiagonal(op[], diag[], request[]) +end diff --git a/julia/LibCEED.jl/src/QFunction.jl b/julia/LibCEED.jl/src/QFunction.jl new file mode 100644 index 0000000000..da5c769cc9 --- /dev/null +++ b/julia/LibCEED.jl/src/QFunction.jl @@ -0,0 +1,148 @@ +abstract type AbstractQFunction end + +struct QFunctionNone <: AbstractQFunction end +Base.getindex(::QFunctionNone) = C.CEED_QFUNCTION_NONE[] + +""" + QFunction + +A libCEED `CeedQFunction` object, typically created using the [`@interior_qf`](@ref) macro. + +A `QFunction` can also be created from the "Q-function gallery" using +[`create_interior_qfunction`](@ref). The identity Q-function can be created using +[`create_identity_qfunction`](@ref). +""" +mutable struct QFunction <: AbstractQFunction + ref::RefValue{C.CeedQFunction} + user_qf::Union{Nothing,UserQFunction} + ctx::Union{Nothing,Context} + function QFunction(ref, user_qf) + obj = new(ref, user_qf, nothing) + finalizer(obj) do x + # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x)) + destroy(x) + end + return obj + end +end +QFunction(ref::Ref{C.CeedQFunction}) = QFunction(ref, nothing) +destroy(qf::QFunction) = C.CeedQFunctionDestroy(qf.ref) # COV_EXCL_LINE +Base.getindex(qf::QFunction) = qf.ref[] +Base.show(io::IO, ::MIME"text/plain", qf::QFunction) = + ceed_show(io, qf, C.CeedQFunctionView) + +function create_interior_qfunction(c::Ceed, f::UserQFunction; vlength=1) + ref = Ref{C.CeedQFunction}() + # Use empty string as source location to indicate to libCEED that there is + # no C source for this Q-function + C.CeedQFunctionCreateInterior(c[], vlength, f.fptr, "", ref) + # COV_EXCL_START + if !isnothing(f.cuf) + C.CeedQFunctionSetCUDAUserFunction(ref[], f.cuf) + elseif iscuda(c) && !cuda_is_loaded + error(string( + "In order to use user Q-functions with a CUDA backend, the CUDA.jl package ", + "must be loaded", + )) + end + # COV_EXCL_STOP + QFunction(ref, f) +end + +""" + create_interior_qfunction(ceed::Ceed, name::AbstractString) + +Create a [`QFunction`](@ref) from the Q-function gallery, using the provided name. + +# Examples + +- Build and apply the 3D mass operator +``` +build_mass_qf = create_interior_qfunction(c, "Mass3DBuild") +apply_mass_qf = create_interior_qfunction(c, "MassApply") +``` +- Build and apply the 3D Poisson operator +``` +build_poi_qf = create_interior_qfunction(c, "Poisson3DBuild") +apply_poi_qf = create_interior_qfunction(c, "Poisson3DApply") +``` +""" +function create_interior_qfunction(c::Ceed, name::AbstractString) + ref = Ref{C.CeedQFunction}() + C.CeedQFunctionCreateInteriorByName(c.ref[], name, ref) + QFunction(ref) +end + +""" + create_identity_qfunction(c::Ceed, size, inmode::EvalMode, outmode::EvalMode) + +Create an identity [`QFunction`](@ref). Inputs are written into outputs in the order given. +This is useful for [`Operators`](@ref Operator) that can be represented with only the action +of a [`ElemRestriction`](@ref) and [`Basis`](@ref), such as restriction and prolongation +operators for p-multigrid. Backends may optimize `CeedOperators` with this Q-function to +avoid the copy of input data to output fields by using the same memory location for both. +""" +function create_identity_qfunction(c::Ceed, size, inmode::EvalMode, outmode::EvalMode) + ref = Ref{C.CeedQFunction}() + C.CeedQFunctionCreateIdentity(c[], size, inmode, outmode, ref) + QFunction(ref) +end + +function add_input!(qf::AbstractQFunction, name::AbstractString, size, emode) + C.CeedQFunctionAddInput(qf[], name, size, emode) +end + +function add_output!(qf::AbstractQFunction, name::AbstractString, size, emode) + C.CeedQFunctionAddOutput(qf[], name, size, emode) +end + +""" + set_context!(qf::QFunction, ctx::Context) + +Associate a [`Context`](@ref) object `ctx` with the given Q-function `qf`. +""" +function set_context!(qf::QFunction, ctx) + # Preserve the context data from the GC by storing a reference + qf.ctx = ctx + C.CeedQFunctionSetContext(qf[], ctx[]) +end + +function get_field_sizes(qf::AbstractQFunction) + ninputs = Ref{CeedInt}() + noutputs = Ref{CeedInt}() + + C.CeedQFunctionGetNumArgs(qf[], ninputs, noutputs) + + inputs = Ref{Ptr{C.CeedQFunctionField}}() + outputs = Ref{Ptr{C.CeedQFunctionField}}() + C.CeedQFunctionGetFields(qf[], inputs, outputs) + + input_sizes = zeros(CeedInt, ninputs[]) + output_sizes = zeros(CeedInt, noutputs[]) + + for i = 1:ninputs[] + field = unsafe_load(inputs[], i) + C.CeedQFunctionFieldGetSize(field, pointer(input_sizes, i)) + end + + for i = 1:noutputs[] + field = unsafe_load(outputs[], i) + C.CeedQFunctionFieldGetSize(field, pointer(output_sizes, i)) + end + + input_sizes, output_sizes +end + +""" + apply!(qf::QFunction, Q, vin, vout) + +Apply the action of a [`QFunction`](@ref) to an array of input vectors, and store the result +in an array of output vectors. +""" +function apply!(qf::QFunction, Q, vin, vout) + vins = map(x -> x[], vin) + vouts = map(x -> x[], vout) + GC.@preserve vin vout begin + C.CeedQFunctionApply(qf[], Q, pointer(vins), pointer(vouts)) + end +end diff --git a/julia/LibCEED.jl/src/Quadrature.jl b/julia/LibCEED.jl/src/Quadrature.jl new file mode 100644 index 0000000000..610fcbbe64 --- /dev/null +++ b/julia/LibCEED.jl/src/Quadrature.jl @@ -0,0 +1,37 @@ +@doc raw""" + gauss_quadrature(q) + +Return the Gauss-Legendre quadrature rule with `q` points (integrates polynomials of degree +$2q-1$ exactly). + +A tuple `(x,w)` is returned. +""" +function gauss_quadrature(q) + x = zeros(CeedScalar, q) + w = zeros(CeedScalar, q) + C.CeedGaussQuadrature(q, x, w) + x, w +end + +struct QuadratureMode{T} end +const Abscissa = QuadratureMode{:Abscissa}() +const AbscissaAndWeights = QuadratureMode{:AbscissaAndWeights}() + +@doc raw""" + lobatto_quadrature(q, mode::Mode=Abscissa) + +Return the Gauss-Lobatto quadrature rule with `q` points (integrates polynomials of degree +$2q-3$ exactly). + +If `mode` is `AbscissaAndWeights`, then both the weights and abscissa are returned as a +tuple `(x,w)`. + +Otherwise, (if `mode` is `Abscissa`), then only the abscissa `x` are returned. +""" +function lobatto_quadrature(q, mode::Mode=Abscissa) where {Mode} + return_weights = (mode == AbscissaAndWeights) + x = zeros(CeedScalar, q) + w = (return_weights) ? zeros(CeedScalar, q) : C_NULL + C.CeedLobattoQuadrature(q, x, w) + return_weights ? (x, w) : x +end diff --git a/julia/LibCEED.jl/src/Request.jl b/julia/LibCEED.jl/src/Request.jl new file mode 100644 index 0000000000..61ea58e17d --- /dev/null +++ b/julia/LibCEED.jl/src/Request.jl @@ -0,0 +1,20 @@ +abstract type AbstractRequest end + +struct RequestImmediate <: AbstractRequest end +Base.getindex(::RequestImmediate) = C.CEED_REQUEST_IMMEDIATE[] + +struct RequestOrdered <: AbstractRequest end +Base.getindex(::RequestOrdered) = C.CEED_REQUEST_ORDERED[] + +#= +# CeedRequest is not fully implemented in libCEED. When it is implemented, the +# following can be used as a starting point for the Julia interface. + +struct Request <: AbstractRequest + ref::RefValue{C.CeedRequest} +end + +Request() = Request(Ref{C.CeedRequest}()) + +Base.wait(req::Request) = C.CeedRequestWait(req[]) +=# diff --git a/julia/LibCEED.jl/src/UserQFunction.jl b/julia/LibCEED.jl/src/UserQFunction.jl new file mode 100644 index 0000000000..17d9603207 --- /dev/null +++ b/julia/LibCEED.jl/src/UserQFunction.jl @@ -0,0 +1,278 @@ +struct UserQFunction{F,K} + f::F + fptr::Ptr{Nothing} + kf::K + cuf::Union{Nothing,Ptr{Nothing}} +end + +@inline function extract_context(ptr, ::Type{T}) where {T} + unsafe_load(Ptr{T}(ptr)) +end + +@inline function extract_array(ptr, idx, dims) + UnsafeArray(Ptr{CeedScalar}(unsafe_load(ptr, idx)), dims) +end + +function generate_user_qfunction( + ceed, + def_module, + qf_name, + constants, + array_names, + ctx, + dims_in, + dims_out, + body, +) + idx = gensym(:i) + Q = gensym(:Q) + ctx_ptr = gensym(:ctx_ptr) + in_ptr = gensym(:in_ptr) + out_ptr = gensym(:out_ptr) + + const_assignments = Vector{Expr}(undef, length(constants)) + for (i, c) ∈ enumerate(constants) + const_assignments[i] = :($(c[1]) = $(c[2])) + end + + narrays = length(array_names) + arrays = Vector{Expr}(undef, narrays) + array_views = Vector{Expr}(undef, narrays) + n_in = length(dims_in) + for (i, arr_name) ∈ enumerate(array_names) + i_inout = (i <= n_in) ? i : i - n_in + dims = (i <= n_in) ? dims_in[i] : dims_out[i-n_in] + ptr = (i <= n_in) ? in_ptr : out_ptr + arr_name_gen = gensym(arr_name) + arrays[i] = :($arr_name_gen = extract_array($ptr, $i_inout, (Int($Q), $(dims...)))) + ndims = length(dims) + slice = Expr(:ref, arr_name_gen, idx, (:(:) for i = 1:ndims)...) + if i <= n_in + if ndims == 0 + array_views[i] = :($arr_name = $slice) + else + S = Tuple{dims...} + array_views[i] = :($arr_name = LibCEED.SArray{$S}(@view $slice)) + end + else + array_views[i] = :($arr_name = @view $slice) + end + end + + if isnothing(ctx) + ctx_assignment = nothing + else + ctx_assignment = :($(ctx.name) = extract_context($ctx_ptr, $(ctx.type))) + end + + qf1 = gensym(qf_name) + f = Core.eval( + def_module, + quote + @inline function $qf1( + $ctx_ptr::Ptr{Cvoid}, + $Q::CeedInt, + $in_ptr::Ptr{Ptr{CeedScalar}}, + $out_ptr::Ptr{Ptr{CeedScalar}}, + ) + $(const_assignments...) + $ctx_assignment + $(arrays...) + @inbounds @simd for $idx = 1:$Q + $(array_views...) + $body + end + CeedInt(0) + end + end, + ) + f_qn = QuoteNode(f) + rt = :CeedInt + at = :(Core.svec(Ptr{Cvoid}, CeedInt, Ptr{Ptr{CeedScalar}}, Ptr{Ptr{CeedScalar}})) + fptr = eval(Expr(:cfunction, Ptr{Cvoid}, f_qn, rt, at, QuoteNode(:ccall))) + + # COV_EXCL_START + if iscuda(ceed) + getresource(ceed) == "/gpu/cuda/gen" && error(string( + "/gpu/cuda/gen is not compatible with user Q-functions defined with ", + "libCEED.jl.\nPlease use a different backend, for example: /gpu/cuda/shared ", + "or /gpu/cuda/ref", + )) + if cuda_is_loaded + !has_cuda() && error("No valid CUDA installation found") + qf2 = gensym(qf_name) + kf = Core.eval( + def_module, + quote + @inline function $qf2($ctx_ptr::Ptr{Cvoid}, $(array_names...)) + $(const_assignments...) + $ctx_assignment + $body + nothing + end + end, + ) + cuf = mk_cufunction(ceed, def_module, qf_name, kf, dims_in, dims_out) + else + error(string( + "User Q-functions with CUDA backends require the CUDA.jl package to be ", + "loaded.\nThe libCEED backend is: $(getresource(ceed))\n", + "Please ensure that the CUDA.jl package is installed and loaded.", + )) + end + else + kf = nothing + cuf = nothing + end + # COV_EXCL_STOP + + UserQFunction(f, fptr, kf, cuf) +end + +function meta_user_qfunction(ceed, def_module, qf, args) + qf_name = Meta.quot(qf) + + ctx = nothing + constants = Expr[] + dims_in = Expr[] + dims_out = Expr[] + names_in = Symbol[] + names_out = Symbol[] + + for a ∈ args[1:end-1] + if Meta.isexpr(a, :(=)) + a1 = Meta.quot(a.args[1]) + a2 = esc(a.args[2]) + push!(constants, :(($a1, $a2))) + elseif Meta.isexpr(a, :tuple) + arr_name = a.args[1] + inout = a.args[2].value + ndim = length(a.args) - 3 + dims = Vector{Expr}(undef, ndim) + for d = 1:ndim + dims[d] = :(Int($(a.args[d+3]))) + end + dims_expr = :(Int[$(esc.(a.args[4:end])...)]) + if inout == :in + push!(dims_in, dims_expr) + push!(names_in, arr_name) + elseif inout == :out + push!(dims_out, dims_expr) + push!(names_out, arr_name) + else + error("Array specification must be either :in or :out. Given $inout.") + end + elseif Meta.isexpr(a, :(::)) + ctx = (name=a.args[1], type=a.args[2]) + else + error("Bad argument to @user_qfunction") + end + end + + body = Meta.quot(args[end]) + + return :(generate_user_qfunction( + $ceed, + $def_module, + $qf_name, + [$(constants...)], + $([names_in; names_out]), + $ctx, + [$(dims_in...)], + [$(dims_out...)], + $body, + )) +end + +""" + @interior_qf name=def + +Creates a user-defined interior (volumetric) Q-function, and assigns it to a variable named +`name`. The definition of the Q-function is given as: +``` +@interior_qf user_qf=( + ceed::CEED, + [const1=val1, const2=val2, ...], + [ctx::ContextType], + (I1, :in, EvalMode, dims...), + (I2, :in, EvalMode, dims...), + (O1, :out, EvalMode, dims...), + body +) +``` +The definitions of form `const=val` are used for definitions which will be compile-time +constants in the Q-function. For example, if `dim` is a variable set to the dimension of the +problem, then `dim=dim` will make `dim` available in the body of the Q-function as a +compile-time constant. + +If the user wants to provide a context struct to the Q-function, that can be achieved by +optionally including `ctx::ContextType`, where `ContextType` is the type of the context +struct, and `ctx` is the name to which is will be bound in the body of the Q-function. + +This is followed by the definition of the input and output arrays, which take the form +`(arr_name, (:in|:out), EvalMode, dims...)`. Each array will be bound to a variable named +`arr_name`. Input arrays should be tagged with :in, and output arrays with :out. An +`EvalMode` should be specified, followed by the dimensions of the array. If the array +consists of scalars (one number per Q-point) then `dims` should be omitted. + +# Examples + +- Q-function to compute the "Q-data" for the mass operator, which is given by the quadrature + weight times the Jacobian determinant. The mesh Jacobian (the gradient of the nodal mesh + points) and the quadrature weights are given as input arrays, and the Q-data is the output + array. `dim` is given as a compile-time constant, and so the array `J` is statically + sized, and therefore `det(J)` will automatically dispatch to an optimized implementation + for the given dimension. +``` +@interior_qf build_qfunc = ( + ceed, dim=dim, + (J, :in, EVAL_GRAD, dim, dim), + (w, :in, EVAL_WEIGHT), + (qdata, :out, EVAL_NONE), + qdata[] = w*det(J) +) +``` +""" +macro interior_qf(args) + if !Meta.isexpr(args, :(=)) + error("@interior_qf must be of form `qf = (body)`") # COV_EXCL_LINE + end + + qf = args.args[1] + user_qf = esc(qf) + args = args.args[2].args + ceed = esc(args[1]) + + # Calculate field sizes + fields_in = Expr[] + fields_out = Expr[] + for a ∈ args + if Meta.isexpr(a, :tuple) + field_name = String(a.args[1]) + inout = a.args[2].value + evalmode = a.args[3] + ndim = length(a.args) - 3 + dims = Vector{Expr}(undef, ndim) + for d = 1:ndim + dims[d] = esc(:(Int($(a.args[d+3])))) + end + sz_expr = :(prod(($(dims...),))) + if inout == :in + push!(fields_in, :(add_input!($user_qf, $field_name, $sz_expr, $evalmode))) + elseif inout == :out + push!( + fields_out, + :(add_output!($user_qf, $field_name, $sz_expr, $evalmode)), + ) + end + end + end + + gen_user_qf = meta_user_qfunction(ceed, __module__, qf, args[2:end]) + + quote + $user_qf = create_interior_qfunction($ceed, $gen_user_qf) + $(fields_in...) + $(fields_out...) + end +end diff --git a/julia/LibCEED.jl/src/generated/generate_bindings.jl b/julia/LibCEED.jl/src/generated/generate_bindings.jl new file mode 100644 index 0000000000..d93ca73a2d --- /dev/null +++ b/julia/LibCEED.jl/src/generated/generate_bindings.jl @@ -0,0 +1,28 @@ +using Clang +using Clang.LibClang.Clang_jll + +""" + generate_ceed_wrapper(ceed_path) + +Generate Julia bindings for the libCEED library, where `ceed_path` is the path +of the libCEED directory. The generated bindings are used (with some manual +modifications) to create the low-level C interface LibCEED.C for the LibCEED.jl +package. +""" +function generate_ceed_wrapper(ceed_path) + headers = ["ceed.h", "ceed-cuda.h", "ceed-backend.h"] + ceed_include = joinpath(ceed_path, "include") + ceed_headers = [joinpath(ceed_include, header) for header in headers] + + wc = init(; + headers=ceed_headers, + output_file=joinpath(@__DIR__, "libceed_api_gen.jl"), + common_file=joinpath(@__DIR__, "libceed_common_gen.jl"), + clang_includes=[ceed_include, CLANG_INCLUDE], + clang_args=map(x -> "-I"*x, find_std_headers()), + header_wrapped=(root, current) -> root == current, + header_library=x -> "libceed", + clang_diagnostics=true, + ) + run(wc, false) +end diff --git a/julia/LibCEED.jl/src/generated/libceed_api.jl b/julia/LibCEED.jl/src/generated/libceed_api.jl new file mode 100644 index 0000000000..a647d5a434 --- /dev/null +++ b/julia/LibCEED.jl/src/generated/libceed_api.jl @@ -0,0 +1,773 @@ +# Julia wrapper for header: ceed.h +# Automatically generated using Clang.jl +#! format: off + +function CeedInit(resource, ceed) + ccall((:CeedInit, libceed), Cint, (Cstring, Ptr{Ceed}), resource, ceed) +end + +function CeedGetResource(ceed, resource) + ccall((:CeedGetResource, libceed), Cint, (Ceed, Ptr{Cstring}), ceed, resource) +end + +function CeedIsDeterministic(ceed, isDeterministic) + ccall((:CeedIsDeterministic, libceed), Cint, (Ceed, Ptr{Bool}), ceed, isDeterministic) +end + +function CeedView(ceed, stream) + ccall((:CeedView, libceed), Cint, (Ceed, Ptr{FILE}), ceed, stream) +end + +function CeedDestroy(ceed) + ccall((:CeedDestroy, libceed), Cint, (Ptr{Ceed},), ceed) +end + +function CeedErrorReturn(arg1, arg2, arg3, arg4, arg5, arg6, arg7) + ccall((:CeedErrorReturn, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7) +end + +function CeedErrorStore(arg1, arg2, arg3, arg4, arg5, arg6, arg7) + ccall((:CeedErrorStore, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7) +end + +function CeedErrorAbort(arg1, arg2, arg3, arg4, arg5, arg6, arg7) + ccall((:CeedErrorAbort, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7) +end + +function CeedErrorExit(arg1, arg2, arg3, arg4, arg5, arg6, arg7) + ccall((:CeedErrorExit, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7) +end + +function CeedSetErrorHandler(ceed, eh) + ccall((:CeedSetErrorHandler, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, eh) +end + +function CeedGetErrorMessage(arg1, errmsg) + ccall((:CeedGetErrorMessage, libceed), Cint, (Ceed, Ptr{Cstring}), arg1, errmsg) +end + +function CeedResetErrorMessage(arg1, errmsg) + ccall((:CeedResetErrorMessage, libceed), Cint, (Ceed, Ptr{Cstring}), arg1, errmsg) +end + +function CeedGetPreferredMemType(ceed, type) + ccall((:CeedGetPreferredMemType, libceed), Cint, (Ceed, Ptr{CeedMemType}), ceed, type) +end + +function CeedVectorCreate(ceed, len, vec) + ccall((:CeedVectorCreate, libceed), Cint, (Ceed, CeedInt, Ptr{CeedVector}), ceed, len, vec) +end + +function CeedVectorSetArray(vec, mtype, cmode, array) + ccall((:CeedVectorSetArray, libceed), Cint, (CeedVector, CeedMemType, CeedCopyMode, Ptr{CeedScalar}), vec, mtype, cmode, array) +end + +function CeedVectorSetValue(vec, value) + ccall((:CeedVectorSetValue, libceed), Cint, (CeedVector, CeedScalar), vec, value) +end + +function CeedVectorSyncArray(vec, mtype) + ccall((:CeedVectorSyncArray, libceed), Cint, (CeedVector, CeedMemType), vec, mtype) +end + +function CeedVectorTakeArray(vec, mtype, array) + ccall((:CeedVectorTakeArray, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array) +end + +function CeedVectorGetArray(vec, mtype, array) + ccall((:CeedVectorGetArray, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array) +end + +function CeedVectorGetArrayRead(vec, mtype, array) + ccall((:CeedVectorGetArrayRead, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array) +end + +function CeedVectorRestoreArray(vec, array) + ccall((:CeedVectorRestoreArray, libceed), Cint, (CeedVector, Ptr{Ptr{CeedScalar}}), vec, array) +end + +function CeedVectorRestoreArrayRead(vec, array) + ccall((:CeedVectorRestoreArrayRead, libceed), Cint, (CeedVector, Ptr{Ptr{CeedScalar}}), vec, array) +end + +function CeedVectorNorm(vec, type, norm) + ccall((:CeedVectorNorm, libceed), Cint, (CeedVector, CeedNormType, Ptr{CeedScalar}), vec, type, norm) +end + +function CeedVectorReciprocal(vec) + ccall((:CeedVectorReciprocal, libceed), Cint, (CeedVector,), vec) +end + +function CeedVectorView(vec, fpfmt, stream) + ccall((:CeedVectorView, libceed), Cint, (CeedVector, Cstring, Ptr{FILE}), vec, fpfmt, stream) +end + +function CeedVectorGetLength(vec, length) + ccall((:CeedVectorGetLength, libceed), Cint, (CeedVector, Ptr{CeedInt}), vec, length) +end + +function CeedVectorDestroy(vec) + ccall((:CeedVectorDestroy, libceed), Cint, (Ptr{CeedVector},), vec) +end + +function CeedRequestWait(req) + ccall((:CeedRequestWait, libceed), Cint, (Ptr{CeedRequest},), req) +end + +function CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr) + ccall((:CeedElemRestrictionCreate, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr) +end + +function CeedElemRestrictionCreateStrided(ceed, nelem, elemsize, ncomp, lsize, strides, rstr) + ccall((:CeedElemRestrictionCreateStrided, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, ncomp, lsize, strides, rstr) +end + +function CeedElemRestrictionCreateBlocked(ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr) + ccall((:CeedElemRestrictionCreateBlocked, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr) +end + +function CeedElemRestrictionCreateBlockedStrided(ceed, nelem, elemsize, blksize, ncomp, lsize, strides, rstr) + ccall((:CeedElemRestrictionCreateBlockedStrided, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, blksize, ncomp, lsize, strides, rstr) +end + +function CeedElemRestrictionCreateVector(rstr, lvec, evec) + ccall((:CeedElemRestrictionCreateVector, libceed), Cint, (CeedElemRestriction, Ptr{CeedVector}, Ptr{CeedVector}), rstr, lvec, evec) +end + +function CeedElemRestrictionApply(rstr, tmode, u, ru, request) + ccall((:CeedElemRestrictionApply, libceed), Cint, (CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, Ptr{CeedRequest}), rstr, tmode, u, ru, request) +end + +function CeedElemRestrictionApplyBlock(rstr, block, tmode, u, ru, request) + ccall((:CeedElemRestrictionApplyBlock, libceed), Cint, (CeedElemRestriction, CeedInt, CeedTransposeMode, CeedVector, CeedVector, Ptr{CeedRequest}), rstr, block, tmode, u, ru, request) +end + +function CeedElemRestrictionGetCompStride(rstr, compstride) + ccall((:CeedElemRestrictionGetCompStride, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, compstride) +end + +function CeedElemRestrictionGetNumElements(rstr, numelem) + ccall((:CeedElemRestrictionGetNumElements, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numelem) +end + +function CeedElemRestrictionGetElementSize(rstr, elemsize) + ccall((:CeedElemRestrictionGetElementSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, elemsize) +end + +function CeedElemRestrictionGetLVectorSize(rstr, lsize) + ccall((:CeedElemRestrictionGetLVectorSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, lsize) +end + +function CeedElemRestrictionGetNumComponents(rstr, numcomp) + ccall((:CeedElemRestrictionGetNumComponents, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numcomp) +end + +function CeedElemRestrictionGetNumBlocks(rstr, numblk) + ccall((:CeedElemRestrictionGetNumBlocks, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numblk) +end + +function CeedElemRestrictionGetBlockSize(rstr, blksize) + ccall((:CeedElemRestrictionGetBlockSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, blksize) +end + +function CeedElemRestrictionGetMultiplicity(rstr, mult) + ccall((:CeedElemRestrictionGetMultiplicity, libceed), Cint, (CeedElemRestriction, CeedVector), rstr, mult) +end + +function CeedElemRestrictionView(rstr, stream) + ccall((:CeedElemRestrictionView, libceed), Cint, (CeedElemRestriction, Ptr{FILE}), rstr, stream) +end + +function CeedElemRestrictionDestroy(rstr) + ccall((:CeedElemRestrictionDestroy, libceed), Cint, (Ptr{CeedElemRestriction},), rstr) +end + +function CeedBasisCreateTensorH1Lagrange(ceed, dim, ncomp, P, Q, qmode, basis) + ccall((:CeedBasisCreateTensorH1Lagrange, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedQuadMode, Ptr{CeedBasis}), ceed, dim, ncomp, P, Q, qmode, basis) +end + +function CeedBasisCreateTensorH1(ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, basis) + ccall((:CeedBasisCreateTensorH1, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, basis) +end + +function CeedBasisCreateH1(ceed, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, basis) + ccall((:CeedBasisCreateH1, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, basis) +end + +function CeedBasisView(basis, stream) + ccall((:CeedBasisView, libceed), Cint, (CeedBasis, Ptr{FILE}), basis, stream) +end + +function CeedBasisApply(basis, nelem, tmode, emode, u, v) + ccall((:CeedBasisApply, libceed), Cint, (CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector), basis, nelem, tmode, emode, u, v) +end + +function CeedBasisGetDimension(basis, dim) + ccall((:CeedBasisGetDimension, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, dim) +end + +function CeedBasisGetTopology(basis, topo) + ccall((:CeedBasisGetTopology, libceed), Cint, (CeedBasis, Ptr{CeedElemTopology}), basis, topo) +end + +function CeedBasisGetNumComponents(basis, numcomp) + ccall((:CeedBasisGetNumComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, numcomp) +end + +function CeedBasisGetNumNodes(basis, P) + ccall((:CeedBasisGetNumNodes, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P) +end + +function CeedBasisGetNumNodes1D(basis, P1d) + ccall((:CeedBasisGetNumNodes1D, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P1d) +end + +function CeedBasisGetNumQuadraturePoints(basis, Q) + ccall((:CeedBasisGetNumQuadraturePoints, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q) +end + +function CeedBasisGetNumQuadraturePoints1D(basis, Q1d) + ccall((:CeedBasisGetNumQuadraturePoints1D, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q1d) +end + +function CeedBasisGetQRef(basis, qref) + ccall((:CeedBasisGetQRef, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, qref) +end + +function CeedBasisGetQWeights(basis, qweight) + ccall((:CeedBasisGetQWeights, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, qweight) +end + +function CeedBasisGetInterp(basis, interp) + ccall((:CeedBasisGetInterp, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, interp) +end + +function CeedBasisGetInterp1D(basis, interp1d) + ccall((:CeedBasisGetInterp1D, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, interp1d) +end + +function CeedBasisGetGrad(basis, grad) + ccall((:CeedBasisGetGrad, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, grad) +end + +function CeedBasisGetGrad1D(basis, grad1d) + ccall((:CeedBasisGetGrad1D, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, grad1d) +end + +function CeedBasisDestroy(basis) + ccall((:CeedBasisDestroy, libceed), Cint, (Ptr{CeedBasis},), basis) +end + +function CeedGaussQuadrature(Q, qref1d, qweight1d) + ccall((:CeedGaussQuadrature, libceed), Cint, (CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), Q, qref1d, qweight1d) +end + +function CeedLobattoQuadrature(Q, qref1d, qweight1d) + ccall((:CeedLobattoQuadrature, libceed), Cint, (CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), Q, qref1d, qweight1d) +end + +function CeedQRFactorization(ceed, mat, tau, m, n) + ccall((:CeedQRFactorization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt), ceed, mat, tau, m, n) +end + +function CeedSymmetricSchurDecomposition(ceed, mat, lambda, n) + ccall((:CeedSymmetricSchurDecomposition, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, mat, lambda, n) +end + +function CeedSimultaneousDiagonalization(ceed, matA, matB, x, lambda, n) + ccall((:CeedSimultaneousDiagonalization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, matA, matB, x, lambda, n) +end + +function CeedQFunctionCreateInterior(ceed, vlength, f, source, qf) + ccall((:CeedQFunctionCreateInterior, libceed), Cint, (Ceed, CeedInt, CeedQFunctionUser, Cstring, Ptr{CeedQFunction}), ceed, vlength, f, source, qf) +end + +function CeedQFunctionCreateInteriorByName(ceed, name, qf) + ccall((:CeedQFunctionCreateInteriorByName, libceed), Cint, (Ceed, Cstring, Ptr{CeedQFunction}), ceed, name, qf) +end + +function CeedQFunctionCreateIdentity(ceed, size, inmode, outmode, qf) + ccall((:CeedQFunctionCreateIdentity, libceed), Cint, (Ceed, CeedInt, CeedEvalMode, CeedEvalMode, Ptr{CeedQFunction}), ceed, size, inmode, outmode, qf) +end + +function CeedQFunctionAddInput(qf, fieldname, size, emode) + ccall((:CeedQFunctionAddInput, libceed), Cint, (CeedQFunction, Cstring, CeedInt, CeedEvalMode), qf, fieldname, size, emode) +end + +function CeedQFunctionAddOutput(qf, fieldname, size, emode) + ccall((:CeedQFunctionAddOutput, libceed), Cint, (CeedQFunction, Cstring, CeedInt, CeedEvalMode), qf, fieldname, size, emode) +end + +function CeedQFunctionSetContext(qf, ctx) + ccall((:CeedQFunctionSetContext, libceed), Cint, (CeedQFunction, CeedQFunctionContext), qf, ctx) +end + +function CeedQFunctionView(qf, stream) + ccall((:CeedQFunctionView, libceed), Cint, (CeedQFunction, Ptr{FILE}), qf, stream) +end + +function CeedQFunctionApply(qf, Q, u, v) + ccall((:CeedQFunctionApply, libceed), Cint, (CeedQFunction, CeedInt, Ptr{CeedVector}, Ptr{CeedVector}), qf, Q, u, v) +end + +function CeedQFunctionDestroy(qf) + ccall((:CeedQFunctionDestroy, libceed), Cint, (Ptr{CeedQFunction},), qf) +end + +function CeedQFunctionContextCreate(ceed, ctx) + ccall((:CeedQFunctionContextCreate, libceed), Cint, (Ceed, Ptr{CeedQFunctionContext}), ceed, ctx) +end + +function CeedQFunctionContextSetData(ctx, mtype, cmode, size, data) + ccall((:CeedQFunctionContextSetData, libceed), Cint, (CeedQFunctionContext, CeedMemType, CeedCopyMode, Csize_t, Ptr{Cvoid}), ctx, mtype, cmode, size, data) +end + +function CeedQFunctionContextGetData(ctx, mtype, data) + ccall((:CeedQFunctionContextGetData, libceed), Cint, (CeedQFunctionContext, CeedMemType, Ptr{Cvoid}), ctx, mtype, data) +end + +function CeedQFunctionContextRestoreData(ctx, data) + ccall((:CeedQFunctionContextRestoreData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data) +end + +function CeedQFunctionContextView(ctx, stream) + ccall((:CeedQFunctionContextView, libceed), Cint, (CeedQFunctionContext, Ptr{FILE}), ctx, stream) +end + +function CeedQFunctionContextDestroy(ctx) + ccall((:CeedQFunctionContextDestroy, libceed), Cint, (Ptr{CeedQFunctionContext},), ctx) +end + +function CeedOperatorCreate(ceed, qf, dqf, dqfT, op) + ccall((:CeedOperatorCreate, libceed), Cint, (Ceed, CeedQFunction, CeedQFunction, CeedQFunction, Ptr{CeedOperator}), ceed, qf, dqf, dqfT, op) +end + +function CeedCompositeOperatorCreate(ceed, op) + ccall((:CeedCompositeOperatorCreate, libceed), Cint, (Ceed, Ptr{CeedOperator}), ceed, op) +end + +function CeedOperatorSetField(op, fieldname, r, b, v) + ccall((:CeedOperatorSetField, libceed), Cint, (CeedOperator, Cstring, CeedElemRestriction, CeedBasis, CeedVector), op, fieldname, r, b, v) +end + +function CeedCompositeOperatorAddSub(compositeop, subop) + ccall((:CeedCompositeOperatorAddSub, libceed), Cint, (CeedOperator, CeedOperator), compositeop, subop) +end + +function CeedOperatorLinearAssembleQFunction(op, assembled, rstr, request) + ccall((:CeedOperatorLinearAssembleQFunction, libceed), Cint, (CeedOperator, Ptr{CeedVector}, Ptr{CeedElemRestriction}, Ptr{CeedRequest}), op, assembled, rstr, request) +end + +function CeedOperatorLinearAssembleDiagonal(op, assembled, request) + ccall((:CeedOperatorLinearAssembleDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request) +end + +function CeedOperatorLinearAssembleAddDiagonal(op, assembled, request) + ccall((:CeedOperatorLinearAssembleAddDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request) +end + +function CeedOperatorLinearAssemblePointBlockDiagonal(op, assembled, request) + ccall((:CeedOperatorLinearAssemblePointBlockDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request) +end + +function CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request) + ccall((:CeedOperatorLinearAssembleAddPointBlockDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request) +end + +function CeedOperatorMultigridLevelCreate(opFine, PMultFine, rstrCoarse, basisCoarse, opCoarse, opProlong, opRestrict) + ccall((:CeedOperatorMultigridLevelCreate, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, opCoarse, opProlong, opRestrict) +end + +function CeedOperatorMultigridLevelCreateTensorH1(opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict) + ccall((:CeedOperatorMultigridLevelCreateTensorH1, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedScalar}, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict) +end + +function CeedOperatorMultigridLevelCreateH1(opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict) + ccall((:CeedOperatorMultigridLevelCreateH1, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedScalar}, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict) +end + +function CeedOperatorCreateFDMElementInverse(op, fdminv, request) + ccall((:CeedOperatorCreateFDMElementInverse, libceed), Cint, (CeedOperator, Ptr{CeedOperator}, Ptr{CeedRequest}), op, fdminv, request) +end + +function CeedOperatorView(op, stream) + ccall((:CeedOperatorView, libceed), Cint, (CeedOperator, Ptr{FILE}), op, stream) +end + +function CeedOperatorApply(op, in, out, request) + ccall((:CeedOperatorApply, libceed), Cint, (CeedOperator, CeedVector, CeedVector, Ptr{CeedRequest}), op, in, out, request) +end + +function CeedOperatorApplyAdd(op, in, out, request) + ccall((:CeedOperatorApplyAdd, libceed), Cint, (CeedOperator, CeedVector, CeedVector, Ptr{CeedRequest}), op, in, out, request) +end + +function CeedOperatorDestroy(op) + ccall((:CeedOperatorDestroy, libceed), Cint, (Ptr{CeedOperator},), op) +end + +function CeedIntPow(base, power) + ccall((:CeedIntPow, libceed), CeedInt, (CeedInt, CeedInt), base, power) +end + +function CeedIntMin(a, b) + ccall((:CeedIntMin, libceed), CeedInt, (CeedInt, CeedInt), a, b) +end + +function CeedIntMax(a, b) + ccall((:CeedIntMax, libceed), CeedInt, (CeedInt, CeedInt), a, b) +end +# Julia wrapper for header: ceed-cuda.h +# Automatically generated using Clang.jl + + +function CeedQFunctionSetCUDAUserFunction(qf, f) + ccall((:CeedQFunctionSetCUDAUserFunction, libceed), Cint, (CeedQFunction, Cint), qf, f) +end +# Julia wrapper for header: ceed-backend.h +# Automatically generated using Clang.jl + + +function CeedMallocArray(n, unit, p) + ccall((:CeedMallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p) +end + +function CeedCallocArray(n, unit, p) + ccall((:CeedCallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p) +end + +function CeedReallocArray(n, unit, p) + ccall((:CeedReallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p) +end + +function CeedFree(p) + ccall((:CeedFree, libceed), Cint, (Ptr{Cvoid},), p) +end + +function CeedRegister(prefix, init, priority) + ccall((:CeedRegister, libceed), Cint, (Cstring, Ptr{Cvoid}, UInt32), prefix, init, priority) +end + +function CeedIsDebug(ceed, isDebug) + ccall((:CeedIsDebug, libceed), Cint, (Ceed, Ptr{Bool}), ceed, isDebug) +end + +function CeedGetParent(ceed, parent) + ccall((:CeedGetParent, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, parent) +end + +function CeedGetDelegate(ceed, delegate) + ccall((:CeedGetDelegate, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, delegate) +end + +function CeedSetDelegate(ceed, delegate) + ccall((:CeedSetDelegate, libceed), Cint, (Ceed, Ceed), ceed, delegate) +end + +function CeedGetObjectDelegate(ceed, delegate, objname) + ccall((:CeedGetObjectDelegate, libceed), Cint, (Ceed, Ptr{Ceed}, Cstring), ceed, delegate, objname) +end + +function CeedSetObjectDelegate(ceed, delegate, objname) + ccall((:CeedSetObjectDelegate, libceed), Cint, (Ceed, Ceed, Cstring), ceed, delegate, objname) +end + +function CeedGetOperatorFallbackResource(ceed, resource) + ccall((:CeedGetOperatorFallbackResource, libceed), Cint, (Ceed, Ptr{Cstring}), ceed, resource) +end + +function CeedSetOperatorFallbackResource(ceed, resource) + ccall((:CeedSetOperatorFallbackResource, libceed), Cint, (Ceed, Cstring), ceed, resource) +end + +function CeedGetOperatorFallbackParentCeed(ceed, parent) + ccall((:CeedGetOperatorFallbackParentCeed, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, parent) +end + +function CeedSetDeterministic(ceed, isDeterministic) + ccall((:CeedSetDeterministic, libceed), Cint, (Ceed, Bool), ceed, isDeterministic) +end + +function CeedSetBackendFunction(ceed, type, object, fname, f) + ccall((:CeedSetBackendFunction, libceed), Cint, (Ceed, Cstring, Ptr{Cvoid}, Cstring, Ptr{Cvoid}), ceed, type, object, fname, f) +end + +function CeedGetData(ceed, data) + ccall((:CeedGetData, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, data) +end + +function CeedSetData(ceed, data) + ccall((:CeedSetData, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, data) +end + +function CeedVectorGetCeed(vec, ceed) + ccall((:CeedVectorGetCeed, libceed), Cint, (CeedVector, Ptr{Ceed}), vec, ceed) +end + +function CeedVectorGetState(vec, state) + ccall((:CeedVectorGetState, libceed), Cint, (CeedVector, Ptr{UInt64}), vec, state) +end + +function CeedVectorAddReference(vec) + ccall((:CeedVectorAddReference, libceed), Cint, (CeedVector,), vec) +end + +function CeedVectorGetData(vec, data) + ccall((:CeedVectorGetData, libceed), Cint, (CeedVector, Ptr{Cvoid}), vec, data) +end + +function CeedVectorSetData(vec, data) + ccall((:CeedVectorSetData, libceed), Cint, (CeedVector, Ptr{Cvoid}), vec, data) +end + +function CeedElemRestrictionGetCeed(rstr, ceed) + ccall((:CeedElemRestrictionGetCeed, libceed), Cint, (CeedElemRestriction, Ptr{Ceed}), rstr, ceed) +end + +function CeedElemRestrictionGetStrides(rstr, strides) + ccall((:CeedElemRestrictionGetStrides, libceed), Cint, (CeedElemRestriction, Ptr{NTuple{3, CeedInt}}), rstr, strides) +end + +function CeedElemRestrictionGetOffsets(rstr, mtype, offsets) + ccall((:CeedElemRestrictionGetOffsets, libceed), Cint, (CeedElemRestriction, CeedMemType, Ptr{Ptr{CeedInt}}), rstr, mtype, offsets) +end + +function CeedElemRestrictionRestoreOffsets(rstr, offsets) + ccall((:CeedElemRestrictionRestoreOffsets, libceed), Cint, (CeedElemRestriction, Ptr{Ptr{CeedInt}}), rstr, offsets) +end + +function CeedElemRestrictionIsStrided(rstr, isstrided) + ccall((:CeedElemRestrictionIsStrided, libceed), Cint, (CeedElemRestriction, Ptr{Bool}), rstr, isstrided) +end + +function CeedElemRestrictionHasBackendStrides(rstr, hasbackendstrides) + ccall((:CeedElemRestrictionHasBackendStrides, libceed), Cint, (CeedElemRestriction, Ptr{Bool}), rstr, hasbackendstrides) +end + +function CeedElemRestrictionGetELayout(rstr, layout) + ccall((:CeedElemRestrictionGetELayout, libceed), Cint, (CeedElemRestriction, Ptr{NTuple{3, CeedInt}}), rstr, layout) +end + +function CeedElemRestrictionSetELayout(rstr, layout) + ccall((:CeedElemRestrictionSetELayout, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, layout) +end + +function CeedElemRestrictionGetData(rstr, data) + ccall((:CeedElemRestrictionGetData, libceed), Cint, (CeedElemRestriction, Ptr{Cvoid}), rstr, data) +end + +function CeedElemRestrictionSetData(rstr, data) + ccall((:CeedElemRestrictionSetData, libceed), Cint, (CeedElemRestriction, Ptr{Cvoid}), rstr, data) +end + +function CeedBasisGetCollocatedGrad(basis, colograd1d) + ccall((:CeedBasisGetCollocatedGrad, libceed), Cint, (CeedBasis, Ptr{CeedScalar}), basis, colograd1d) +end + +function CeedHouseholderApplyQ(A, Q, tau, tmode, m, n, k, row, col) + ccall((:CeedHouseholderApplyQ, libceed), Cint, (Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt), A, Q, tau, tmode, m, n, k, row, col) +end + +function CeedBasisGetCeed(basis, ceed) + ccall((:CeedBasisGetCeed, libceed), Cint, (CeedBasis, Ptr{Ceed}), basis, ceed) +end + +function CeedBasisIsTensor(basis, istensor) + ccall((:CeedBasisIsTensor, libceed), Cint, (CeedBasis, Ptr{Bool}), basis, istensor) +end + +function CeedBasisGetData(basis, data) + ccall((:CeedBasisGetData, libceed), Cint, (CeedBasis, Ptr{Cvoid}), basis, data) +end + +function CeedBasisSetData(basis, data) + ccall((:CeedBasisSetData, libceed), Cint, (CeedBasis, Ptr{Cvoid}), basis, data) +end + +function CeedBasisGetTopologyDimension(topo, dim) + ccall((:CeedBasisGetTopologyDimension, libceed), Cint, (CeedElemTopology, Ptr{CeedInt}), topo, dim) +end + +function CeedBasisGetTensorContract(basis, contract) + ccall((:CeedBasisGetTensorContract, libceed), Cint, (CeedBasis, Ptr{CeedTensorContract}), basis, contract) +end + +function CeedBasisSetTensorContract(basis, contract) + ccall((:CeedBasisSetTensorContract, libceed), Cint, (CeedBasis, Ptr{CeedTensorContract}), basis, contract) +end + +function CeedTensorContractCreate(ceed, basis, contract) + ccall((:CeedTensorContractCreate, libceed), Cint, (Ceed, CeedBasis, Ptr{CeedTensorContract}), ceed, basis, contract) +end + +function CeedTensorContractApply(contract, A, B, C, J, t, tmode, 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, tmode, Add, u, v) +end + +function CeedTensorContractGetCeed(contract, ceed) + ccall((:CeedTensorContractGetCeed, libceed), Cint, (CeedTensorContract, Ptr{Ceed}), contract, ceed) +end + +function CeedTensorContractGetData(contract, data) + ccall((:CeedTensorContractGetData, libceed), Cint, (CeedTensorContract, Ptr{Cvoid}), contract, data) +end + +function CeedTensorContractSetData(contract, data) + ccall((:CeedTensorContractSetData, libceed), Cint, (CeedTensorContract, Ptr{Cvoid}), contract, data) +end + +function CeedTensorContractDestroy(contract) + ccall((:CeedTensorContractDestroy, libceed), Cint, (Ptr{CeedTensorContract},), contract) +end + +function CeedQFunctionRegister(arg1, arg2, arg3, arg4, init) + ccall((:CeedQFunctionRegister, libceed), Cint, (Cstring, Cstring, CeedInt, CeedQFunctionUser, Ptr{Cvoid}), arg1, arg2, arg3, arg4, init) +end + +function CeedQFunctionSetFortranStatus(qf, status) + ccall((:CeedQFunctionSetFortranStatus, libceed), Cint, (CeedQFunction, Bool), qf, status) +end + +function CeedQFunctionGetCeed(qf, ceed) + ccall((:CeedQFunctionGetCeed, libceed), Cint, (CeedQFunction, Ptr{Ceed}), qf, ceed) +end + +function CeedQFunctionGetVectorLength(qf, vlength) + ccall((:CeedQFunctionGetVectorLength, libceed), Cint, (CeedQFunction, Ptr{CeedInt}), qf, vlength) +end + +function CeedQFunctionGetNumArgs(qf, numinputfields, numoutputfields) + ccall((:CeedQFunctionGetNumArgs, libceed), Cint, (CeedQFunction, Ptr{CeedInt}, Ptr{CeedInt}), qf, numinputfields, numoutputfields) +end + +function CeedQFunctionGetSourcePath(qf, source) + ccall((:CeedQFunctionGetSourcePath, libceed), Cint, (CeedQFunction, Ptr{Cstring}), qf, source) +end + +function CeedQFunctionGetUserFunction(qf, f) + ccall((:CeedQFunctionGetUserFunction, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionUser}), qf, f) +end + +function CeedQFunctionGetContext(qf, ctx) + ccall((:CeedQFunctionGetContext, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionContext}), qf, ctx) +end + +function CeedQFunctionGetInnerContext(qf, ctx) + ccall((:CeedQFunctionGetInnerContext, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionContext}), qf, ctx) +end + +function CeedQFunctionIsIdentity(qf, isidentity) + ccall((:CeedQFunctionIsIdentity, libceed), Cint, (CeedQFunction, Ptr{Bool}), qf, isidentity) +end + +function CeedQFunctionGetData(qf, data) + ccall((:CeedQFunctionGetData, libceed), Cint, (CeedQFunction, Ptr{Cvoid}), qf, data) +end + +function CeedQFunctionSetData(qf, data) + ccall((:CeedQFunctionSetData, libceed), Cint, (CeedQFunction, Ptr{Cvoid}), qf, data) +end + +function CeedQFunctionGetFields(qf, inputfields, outputfields) + ccall((:CeedQFunctionGetFields, libceed), Cint, (CeedQFunction, Ptr{Ptr{CeedQFunctionField}}, Ptr{Ptr{CeedQFunctionField}}), qf, inputfields, outputfields) +end + +function CeedQFunctionFieldGetName(qffield, fieldname) + ccall((:CeedQFunctionFieldGetName, libceed), Cint, (CeedQFunctionField, Ptr{Cstring}), qffield, fieldname) +end + +function CeedQFunctionFieldGetSize(qffield, size) + ccall((:CeedQFunctionFieldGetSize, libceed), Cint, (CeedQFunctionField, Ptr{CeedInt}), qffield, size) +end + +function CeedQFunctionFieldGetEvalMode(qffield, emode) + ccall((:CeedQFunctionFieldGetEvalMode, libceed), Cint, (CeedQFunctionField, Ptr{CeedEvalMode}), qffield, emode) +end + +function CeedQFunctionContextGetCeed(cxt, ceed) + ccall((:CeedQFunctionContextGetCeed, libceed), Cint, (CeedQFunctionContext, Ptr{Ceed}), cxt, ceed) +end + +function CeedQFunctionContextGetState(ctx, state) + ccall((:CeedQFunctionContextGetState, libceed), Cint, (CeedQFunctionContext, Ptr{UInt64}), ctx, state) +end + +function CeedQFunctionContextGetContextSize(ctx, ctxsize) + ccall((:CeedQFunctionContextGetContextSize, libceed), Cint, (CeedQFunctionContext, Ptr{Csize_t}), ctx, ctxsize) +end + +function CeedQFunctionContextGetBackendData(ctx, data) + ccall((:CeedQFunctionContextGetBackendData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data) +end + +function CeedQFunctionContextSetBackendData(ctx, data) + ccall((:CeedQFunctionContextSetBackendData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data) +end + +function CeedOperatorGetCeed(op, ceed) + ccall((:CeedOperatorGetCeed, libceed), Cint, (CeedOperator, Ptr{Ceed}), op, ceed) +end + +function CeedOperatorGetNumElements(op, numelem) + ccall((:CeedOperatorGetNumElements, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numelem) +end + +function CeedOperatorGetNumQuadraturePoints(op, numqpts) + ccall((:CeedOperatorGetNumQuadraturePoints, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numqpts) +end + +function CeedOperatorGetNumArgs(op, numargs) + ccall((:CeedOperatorGetNumArgs, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numargs) +end + +function CeedOperatorIsSetupDone(op, issetupdone) + ccall((:CeedOperatorIsSetupDone, libceed), Cint, (CeedOperator, Ptr{Bool}), op, issetupdone) +end + +function CeedOperatorGetQFunction(op, qf) + ccall((:CeedOperatorGetQFunction, libceed), Cint, (CeedOperator, Ptr{CeedQFunction}), op, qf) +end + +function CeedOperatorIsComposite(op, iscomposite) + ccall((:CeedOperatorIsComposite, libceed), Cint, (CeedOperator, Ptr{Bool}), op, iscomposite) +end + +function CeedOperatorGetNumSub(op, numsub) + ccall((:CeedOperatorGetNumSub, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numsub) +end + +function CeedOperatorGetSubList(op, suboperators) + ccall((:CeedOperatorGetSubList, libceed), Cint, (CeedOperator, Ptr{Ptr{CeedOperator}}), op, suboperators) +end + +function CeedOperatorGetData(op, data) + ccall((:CeedOperatorGetData, libceed), Cint, (CeedOperator, Ptr{Cvoid}), op, data) +end + +function CeedOperatorSetData(op, data) + ccall((:CeedOperatorSetData, libceed), Cint, (CeedOperator, Ptr{Cvoid}), op, data) +end + +function CeedOperatorSetSetupDone(op) + ccall((:CeedOperatorSetSetupDone, libceed), Cint, (CeedOperator,), op) +end + +function CeedOperatorGetFields(op, inputfields, outputfields) + ccall((:CeedOperatorGetFields, libceed), Cint, (CeedOperator, Ptr{Ptr{CeedOperatorField}}, Ptr{Ptr{CeedOperatorField}}), op, inputfields, outputfields) +end + +function CeedOperatorFieldGetElemRestriction(opfield, rstr) + ccall((:CeedOperatorFieldGetElemRestriction, libceed), Cint, (CeedOperatorField, Ptr{CeedElemRestriction}), opfield, rstr) +end + +function CeedOperatorFieldGetBasis(opfield, basis) + ccall((:CeedOperatorFieldGetBasis, libceed), Cint, (CeedOperatorField, Ptr{CeedBasis}), opfield, basis) +end + +function CeedOperatorFieldGetVector(opfield, vec) + ccall((:CeedOperatorFieldGetVector, libceed), Cint, (CeedOperatorField, Ptr{CeedVector}), opfield, vec) +end + +function CeedMatrixMultiply(ceed, matA, matB, matC, m, n, kk) + ccall((:CeedMatrixMultiply, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt, CeedInt), ceed, matA, matB, matC, m, n, kk) +end diff --git a/julia/LibCEED.jl/src/generated/libceed_common.jl b/julia/LibCEED.jl/src/generated/libceed_common.jl new file mode 100644 index 0000000000..aad84aa497 --- /dev/null +++ b/julia/LibCEED.jl/src/generated/libceed_common.jl @@ -0,0 +1,99 @@ +# Automatically generated using Clang.jl +#! format: off + +const FILE = Cvoid + +# Skipping MacroDefinition: CEED_QFUNCTION ( name ) static const char name ## _loc [ ] = __FILE__ ":" # name ; static int name + +# Skipping MacroDefinition: CeedError ( ceed , ecode , ... ) ( CeedErrorImpl ( ( ceed ) , __FILE__ , __LINE__ , __func__ , ( ecode ) , __VA_ARGS__ ) ? : ( ecode ) ) + +const CeedInt = Int32 +const CeedScalar = Cdouble +const Ceed_private = Cvoid +const Ceed = Ptr{Ceed_private} +const CeedRequest_private = Cvoid +const CeedRequest = Ptr{CeedRequest_private} +const CeedVector_private = Cvoid +const CeedVector = Ptr{CeedVector_private} +const CeedElemRestriction_private = Cvoid +const CeedElemRestriction = Ptr{CeedElemRestriction_private} +const CeedBasis_private = Cvoid +const CeedBasis = Ptr{CeedBasis_private} +const CeedQFunction_private = Cvoid +const CeedQFunction = Ptr{CeedQFunction_private} +const CeedQFunctionContext_private = Cvoid +const CeedQFunctionContext = Ptr{CeedQFunctionContext_private} +const CeedOperator_private = Cvoid +const CeedOperator = Ptr{CeedOperator_private} + +@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_LINE = 65536 + CEED_TRIANGLE = 131073 + CEED_QUAD = 131074 + CEED_TET = 196611 + CEED_PYRAMID = 196612 + CEED_PRISM = 196613 + CEED_HEX = 196614 +end + + +const CeedQFunctionUser = Ptr{Cvoid} + +# Skipping MacroDefinition: CEED_INTERN CEED_EXTERN __attribute__ ( ( visibility ( "hidden" ) ) ) + +const CEED_MAX_RESOURCE_LEN = 1024 +const CEED_ALIGN = 64 +const CEED_COMPOSITE_MAX = 16 +const CEED_EPSILON = 1.0e-16 +const CEED_DEBUG_COLOR = 0 + +# Skipping MacroDefinition: CeedDebug1 ( ceed , format , ... ) CeedDebugImpl ( ceed , format , ## __VA_ARGS__ ) +# Skipping MacroDefinition: CeedDebug256 ( ceed , color , ... ) CeedDebugImpl256 ( ceed , color , ## __VA_ARGS__ ) +# Skipping MacroDefinition: CeedDebug ( ... ) CeedDebug256 ( ceed , ( unsigned char ) CEED_DEBUG_COLOR , ## __VA_ARGS__ ) +# Skipping MacroDefinition: CeedChk ( ierr ) do { if ( ierr ) return ierr ; } while ( 0 ) +# Skipping MacroDefinition: CeedMalloc ( n , p ) CeedMallocArray ( ( n ) , sizeof ( * * ( p ) ) , p ) +# Skipping MacroDefinition: CeedCalloc ( n , p ) CeedCallocArray ( ( n ) , sizeof ( * * ( p ) ) , p ) +# Skipping MacroDefinition: CeedRealloc ( n , p ) CeedReallocArray ( ( n ) , sizeof ( * * ( p ) ) , p ) + +const CeedTensorContract_private = Cvoid +const CeedTensorContract = Ptr{CeedTensorContract_private} +const CeedQFunctionField_private = Cvoid +const CeedQFunctionField = Ptr{CeedQFunctionField_private} +const CeedOperatorField_private = Cvoid +const CeedOperatorField = Ptr{CeedOperatorField_private}