From 4a222ceacb8a468cef32aa945b921088179f13ee Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 5 Apr 2023 17:41:46 -0700 Subject: [PATCH] Address PR comments for high-level APIs Add Rust sibling functions and tests, add Python tests, fix Julia versioning for CI. --- .github/workflows/julia-test-with-style.yml | 8 +- include/ceed/ceed.h | 2 +- julia/LibCEED.jl/src/LibCEED.jl | 2 +- .../src/generated/libceed_bindings.jl | 2 +- julia/LibCEED.jl/test/Project.toml | 11 - python/__init__.py | 4 +- python/ceed.py | 57 ++++- python/tests/buildmats.py | 77 +++++++ python/tests/test-3-basis.py | 83 ++++++++ rust/libceed/src/lib.rs | 194 +++++++++++++++++- 10 files changed, 417 insertions(+), 23 deletions(-) delete mode 100644 julia/LibCEED.jl/test/Project.toml diff --git a/.github/workflows/julia-test-with-style.yml b/.github/workflows/julia-test-with-style.yml index d576a900ce..5ef7c22d9c 100644 --- a/.github/workflows/julia-test-with-style.yml +++ b/.github/workflows/julia-test-with-style.yml @@ -33,9 +33,9 @@ jobs: make -j2 LIBCEED_LIB=$(find $PWD/lib -name "libceed.*") pushd julia/LibCEED.jl - echo >> test/Project.toml - echo "[preferences.libCEED_jll]" >> test/Project.toml - echo "libceed_path = \"$LIBCEED_LIB\"" >> test/Project.toml + echo >> Project.toml + echo "[preferences.libCEED_jll]" >> Project.toml + echo "libceed_path = \"$LIBCEED_LIB\"" >> Project.toml [[ "$GITHUB_REF" =~ ^refs/(heads/release|tags/).* ]] || julia --project -e 'import Pkg; Pkg.test("LibCEED"; coverage=true, test_args=["--run-dev-tests"])' - git checkout test/Project.toml && julia --project -e 'import Pkg; Pkg.test("LibCEED")' + git checkout Project.toml && julia --project -e 'import Pkg; Pkg.test("LibCEED")' julia --project=.style/ -e 'import Pkg; Pkg.instantiate()' && julia --project=.style/ .style/ceed_style.jl && git diff --exit-code src test examples diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index a5f1ac15b6..dfe5f4b51b 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -137,7 +137,7 @@ CEED_EXTERN int CeedResetErrorMessage(Ceed, const char **err_msg); #define CEED_VERSION_MAJOR 0 #define CEED_VERSION_MINOR 11 #define CEED_VERSION_PATCH 0 -#define CEED_VERSION_RELEASE true +#define CEED_VERSION_RELEASE false /// Compile-time check that the the current library version is at least as recent as the specified version. /// This macro is typically used in diff --git a/julia/LibCEED.jl/src/LibCEED.jl b/julia/LibCEED.jl/src/LibCEED.jl index aaa8d5d88d..243415c6fe 100644 --- a/julia/LibCEED.jl/src/LibCEED.jl +++ b/julia/LibCEED.jl/src/LibCEED.jl @@ -150,7 +150,7 @@ include("Request.jl") include("Operator.jl") include("Misc.jl") -const minimum_libceed_version = v"0.10.0" +const minimum_libceed_version = v"0.12.0" function __init__() if !ceedversion_ge(minimum_libceed_version) diff --git a/julia/LibCEED.jl/src/generated/libceed_bindings.jl b/julia/LibCEED.jl/src/generated/libceed_bindings.jl index 2f65570b36..946a797c42 100644 --- a/julia/LibCEED.jl/src/generated/libceed_bindings.jl +++ b/julia/LibCEED.jl/src/generated/libceed_bindings.jl @@ -1374,7 +1374,7 @@ const CEED_VERSION_MINOR = 11 const CEED_VERSION_PATCH = 0 -const CEED_VERSION_RELEASE = true +const CEED_VERSION_RELEASE = false # Skipping MacroDefinition: CEED_INTERN extern CEED_VISIBILITY ( hidden ) diff --git a/julia/LibCEED.jl/test/Project.toml b/julia/LibCEED.jl/test/Project.toml deleted file mode 100644 index 96e406909d..0000000000 --- a/julia/LibCEED.jl/test/Project.toml +++ /dev/null @@ -1,11 +0,0 @@ -# A bug in Julia 1.6.0's Pkg causes Preferences to be dropped during `Pkg.test()`, so we work around -# it by explicitly creating a `test/Project.toml` which will correctly communicate any preferences -# through to the child Julia process. X-ref: https://github.com/JuliaLang/Pkg.jl/issues/2500 -# See also test/Project.toml in FFTW.jl. - -[deps] -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Preferences = "21216c6a-2e73-6563-6e65-726566657250" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -libCEED_jll = "762fde13-7596-547b-826d-8223c52d51c1" diff --git a/python/__init__.py b/python/__init__.py index 311b1af3b7..40cc969118 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -7,7 +7,7 @@ from .ceed import Ceed from .ceed_vector import Vector -from .ceed_basis import Basis, BasisTensorH1, BasisTensorH1Lagrange, BasisH1 +from .ceed_basis import Basis, BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction from .ceed_operator import Operator, CompositeOperator @@ -18,7 +18,7 @@ # ------------------------------------------------------------------------------ __all__ = ["Ceed", "Vector", - "Basis", "BasisTensorH1", "BasisTensorH1Lagrange", "BasisH1", + "Basis", "BasisTensorH1", "BasisTensorH1Lagrange", "BasisH1", "BasisHdiv", "BasisHcurl", "ElemRestriction", "StridedElemRestriction", "BlockedElemRestriction", "BlockedStridedelemRestriction", "QFunction", "QFunctionByName", "IdentityQFunction", "Operator", "CompositeOperator", diff --git a/python/ceed.py b/python/ceed.py index 6274c66f8a..9af5c34346 100644 --- a/python/ceed.py +++ b/python/ceed.py @@ -13,7 +13,7 @@ import tempfile from abc import ABC from .ceed_vector import Vector -from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1 +from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction from .ceed_qfunctioncontext import QFunctionContext @@ -356,7 +356,7 @@ def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): *interp: Numpy array holding the row-major (nqpts * nnodes) matrix expressing the values of nodal basis functions at quadrature points - *grad: Numpy array holding the row-major (nqpts * dim * nnodes) + *grad: Numpy array holding the row-major (dim * nqpts * nnodes) matrix expressing the derivatives of nodal basis functions at quadrature points *qref: Numpy array of length (nqpts * dim) holding the locations of @@ -370,6 +370,59 @@ def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): return BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight) + def BasisHdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight): + """Ceed Hdiv Basis: finite element non tensor-product basis for H(div) + discretizations. + + Args: + topo: topology of the element, e.g. hypercube, simplex, etc + ncomp: number of field components (1 for scalar fields) + nnodes: total number of nodes + nqpts: total number of quadrature points + *interp: Numpy array holding the row-major (dim * nqpts * nnodes) + matrix expressing the values of basis functions at + quadrature points + *div: Numpy array holding the row-major (nqpts * nnodes) matrix + expressing the divergence of basis functions at + quadrature points + *qref: Numpy array of length (nqpts * dim) holding the locations of + quadrature points on the reference element [-1, 1] + *qweight: Numpy array of length nnodes holding the quadrature + weights on the reference element + + Returns: + basis: Ceed Basis""" + + return BasisHdiv(self, topo, ncomp, nnodes, nqpts, + interp, div, qref, qweight) + + def BasisHcurl(self, topo, ncomp, nnodes, nqpts, + interp, curl, qref, qweight): + """Ceed Hcurl Basis: finite element non tensor-product basis for H(curl) + discretizations. + + Args: + topo: topology of the element, e.g. hypercube, simplex, etc + ncomp: number of field components (1 for scalar fields) + nnodes: total number of nodes + nqpts: total number of quadrature points + *interp: Numpy array holding the row-major (dim * nqpts * nnodes) + matrix expressing the values of basis functions at + quadrature points + *curl: Numpy array holding the row-major (curlcomp * nqpts * nnodes), + curlcomp = 1 if dim < 3 else dim, matrix expressing the curl + of basis functions at quadrature points + *qref: Numpy array of length (nqpts * dim) holding the locations of + quadrature points on the reference element [-1, 1] + *qweight: Numpy array of length nnodes holding the quadrature + weights on the reference element + + Returns: + basis: Ceed Basis""" + + return BasisHcurl(self, topo, ncomp, nnodes, nqpts, + interp, curl, qref, qweight) + # CeedQFunction def QFunction(self, vlength, f, source): """Ceed QFunction: point-wise operation at quadrature points for diff --git a/python/tests/buildmats.py b/python/tests/buildmats.py index e444c7d6e1..a361c1def9 100644 --- a/python/tests/buildmats.py +++ b/python/tests/buildmats.py @@ -47,3 +47,80 @@ def buildmats(qref, qweight, mat_dtype="float64"): grad[(i + Q) * P + 5] = 2. * (1. * (x2 - 1. / 2.) + x2 * 1.) return interp, grad + + +def buildmatshdiv(qref, qweight, mat_dtype="float64"): + P, Q, dim = 4, 4, 2 + interp = np.empty(dim * P * Q, dtype=mat_dtype) + div = np.empty(P * Q, dtype=mat_dtype) + + qref[0] = -1. / np.sqrt(3.) + qref[1] = qref[0] + qref[2] = qref[0] + qref[3] = -qref[0] + qref[4] = -qref[0] + qref[5] = -qref[0] + qref[6] = qref[0] + qref[7] = qref[0] + qweight[0] = 1. + qweight[1] = 1. + qweight[2] = 1. + qweight[3] = 1. + + # Loop over quadrature points + for i in range(Q): + x1 = qref[0 * Q + i] + x2 = qref[1 * Q + i] + # Interp + interp[(i + 0) * P + 0] = 0. + interp[(i + Q) * P + 0] = 1. - x2 + interp[(i + 0) * P + 1] = x1 - 1. + interp[(i + Q) * P + 1] = 0. + interp[(i + 0) * P + 2] = -x1 + interp[(i + Q) * P + 2] = 0. + interp[(i + 0) * P + 3] = 0. + interp[(i + Q) * P + 3] = x2 + # Div + div[i * P + 0] = -1. + div[i * P + 1] = 1. + div[i * P + 2] = -1. + div[i * P + 3] = 1. + + return interp, div + + +def buildmatshcurl(qref, qweight, mat_dtype="float64"): + P, Q, dim = 3, 4, 2 + interp = np.empty(dim * P * Q, dtype=mat_dtype) + curl = np.empty(P * Q, dtype=mat_dtype) + + qref[0] = 0.2 + qref[1] = 0.6 + qref[2] = 1. / 3. + qref[3] = 0.2 + qref[4] = 0.2 + qref[5] = 0.2 + qref[6] = 1. / 3. + qref[7] = 0.6 + qweight[0] = 25. / 96. + qweight[1] = 25. / 96. + qweight[2] = -27. / 96. + qweight[3] = 25. / 96. + + # Loop over quadrature points + for i in range(Q): + x1 = qref[0 * Q + i] + x2 = qref[1 * Q + i] + # Interp + interp[(i + 0) * P + 0] = -x2 + interp[(i + Q) * P + 0] = x1 + interp[(i + 0) * P + 1] = x2 + interp[(i + Q) * P + 1] = 1. - x1 + interp[(i + 0) * P + 2] = 1. - x2 + interp[(i + Q) * P + 2] = x1 + # Curl + curl[i * P + 0] = 2. + curl[i * P + 1] = -2. + curl[i * P + 2] = -2. + + return interp, curl diff --git a/python/tests/test-3-basis.py b/python/tests/test-3-basis.py index 2a08c23ed3..fe77428dde 100644 --- a/python/tests/test-3-basis.py +++ b/python/tests/test-3-basis.py @@ -303,3 +303,86 @@ def test_323(ceed_resource): assert math.fabs(out_array[1 * Q + i] - value) < 10. * TOL # ------------------------------------------------------------------------------- +# Test interpolation with a 2D Quad non-tensor H(div) basis +# ------------------------------------------------------------------------------- + + +def test_330(ceed_resource): + ceed = libceed.Ceed(ceed_resource) + + P, Q, dim = 4, 4, 2 + + in_array = np.ones(P, dtype=ceed.scalar_type()) + qref = np.empty(dim * Q, dtype=ceed.scalar_type()) + qweight = np.empty(Q, dtype=ceed.scalar_type()) + + interp, div = bm.buildmatshdiv(qref, qweight, libceed.scalar_types[ + libceed.lib.CEED_SCALAR_TYPE]) + + b = ceed.BasisHdiv(libceed.QUAD, 1, P, Q, interp, div, qref, qweight) + + # Interpolate function to quadrature points + in_vec = ceed.Vector(P) + in_vec.set_array(in_array, cmode=libceed.USE_POINTER) + out_vec = ceed.Vector(Q * dim) + out_vec.set_value(0) + + b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) + + # Check values at quadrature points + with out_vec.array_read() as out_array: + for i in range(Q): + assert math.fabs(out_array[0 * Q + i] + 1.) < 10. * TOL + assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL + +# ------------------------------------------------------------------------------- +# Test interpolation with a 2D Simplex non-tensor H(curl) basis +# ------------------------------------------------------------------------------- + + +def test_340(ceed_resource): + ceed = libceed.Ceed(ceed_resource) + + P, Q, dim = 3, 4, 2 + + in_array = np.empty(P, dtype=ceed.scalar_type()) + qref = np.empty(dim * Q, dtype=ceed.scalar_type()) + qweight = np.empty(Q, dtype=ceed.scalar_type()) + + interp, curl = bm.buildmatshcurl(qref, qweight, libceed.scalar_types[ + libceed.lib.CEED_SCALAR_TYPE]) + + b = ceed.BasisHcurl(libceed.TRIANGLE, 1, P, Q, interp, curl, qref, qweight) + + # Interpolate function to quadrature points + in_array[0] = 1. + in_array[1] = 2. + in_array[2] = 1. + + in_vec = ceed.Vector(P) + in_vec.set_array(in_array, cmode=libceed.USE_POINTER) + out_vec = ceed.Vector(Q * dim) + out_vec.set_value(0) + + b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) + + # Check values at quadrature points + with out_vec.array_read() as out_array: + for i in range(Q): + assert math.fabs(out_array[0 * Q + i] - 1.) < 10. * TOL + + # Interpolate function to quadrature points + in_array[0] = -1. + in_array[1] = 1. + in_array[2] = 2. + + out_vec.set_value(0) + + b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) + + # Check values at quadrature points + with out_vec.array_read() as out_array: + for i in range(Q): + assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL + +# ------------------------------------------------------------------------------- diff --git a/rust/libceed/src/lib.rs b/rust/libceed/src/lib.rs index c11bc63da0..a884585169 100755 --- a/rust/libceed/src/lib.rs +++ b/rust/libceed/src/lib.rs @@ -586,7 +586,7 @@ impl Ceed { /// * `nqpts` - Total number of quadrature points /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of /// nodal basis functions at quadrature points - /// * `grad` - Row-major `(nqpts * dim * nnodes)` matrix expressing + /// * `grad` - Row-major `(dim * nqpts * nnodes)` matrix 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]` @@ -707,6 +707,198 @@ impl Ceed { ) } + /// Returns an $H(div)$ Basis + /// + /// # arguments + /// + /// * `topo` - Topology of element, e.g. hypercube, simplex, ect + /// * `ncomp` - Number of field components (1 for scalar fields) + /// * `nnodes` - Total number of nodes + /// * `nqpts` - Total number of quadrature points + /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the + /// values of basis functions at quadrature points + /// * `div` - Row-major `(nqpts * nnodes)` matrix expressing the + /// divergence of basis functions at quadrature points + /// * `qref` - Array of length `nqpts` holding the locations of quadrature + /// points on the reference element `[-1, 1]` + /// * `qweight` - Array of length `nqpts` holding the quadrature weights on + /// the reference element + /// + /// ``` + /// # use libceed::prelude::*; + /// # fn main() -> libceed::Result<()> { + /// # let ceed = libceed::Ceed::default_init(); + /// let interp = [ + /// 0.00000000, + /// -1.57735027, + /// 0.57735027, + /// 0.00000000, + /// 0.00000000, + /// -1.57735027, + /// 0.57735027, + /// 0.00000000, + /// 0.00000000, + /// -1.57735027, + /// 0.57735027, + /// 0.00000000, + /// 0.00000000, + /// -0.42264973, + /// -0.57735027, + /// 0.00000000, + /// 0.42264973, + /// 0.00000000, + /// 0.00000000, + /// 0.57735027, + /// 0.42264973, + /// 0.00000000, + /// 0.00000000, + /// 0.57735027, + /// 1.57735027, + /// 0.00000000, + /// 0.00000000, + /// -0.57735027, + /// 1.57735027, + /// 0.00000000, + /// 0.00000000, + /// -0.57735027, + /// ]; + /// let div = [ + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// -1.00000000, + /// 1.00000000, + /// ]; + /// let qref = [ + /// 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, + /// 0.57735026, + /// ]; + /// let qweight = [1.00000000, 1.00000000, 1.00000000, 1.00000000]; + /// let b = ceed.basis_Hdiv(ElemTopology::Quad, 1, 4, 4, &interp, &div, &qref, &qweight)?; + /// # Ok(()) + /// # } + /// ``` + pub fn basis_Hdiv<'a>( + &self, + topo: ElemTopology, + ncomp: usize, + nnodes: usize, + nqpts: usize, + interp: &[crate::Scalar], + div: &[crate::Scalar], + qref: &[crate::Scalar], + qweight: &[crate::Scalar], + ) -> Result> { + Basis::create_Hdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight) + } + + /// Returns an $H(curl)$ Basis + /// + /// # arguments + /// + /// * `topo` - Topology of element, e.g. hypercube, simplex, ect + /// * `ncomp` - Number of field components (1 for scalar fields) + /// * `nnodes` - Total number of nodes + /// * `nqpts` - Total number of quadrature points + /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the + /// values of basis functions at quadrature points + /// * `curl` - Row-major `(curl_comp * nqpts * nnodes)`, `curl_comp = 1 if + /// dim < 3 else dim` matrix expressing the curl of basis + /// functions at quadrature points + /// * `qref` - Array of length `nqpts` holding the locations of quadrature + /// points on the reference element `[-1, 1]` + /// * `qweight` - Array of length `nqpts` holding the quadrature weights on + /// the reference element + /// + /// ``` + /// # use libceed::prelude::*; + /// # fn main() -> libceed::Result<()> { + /// # let ceed = libceed::Ceed::default_init(); + /// let interp = [ + /// -0.20000000, + /// 0.20000000, + /// 0.80000000, + /// -0.20000000, + /// 0.20000000, + /// 0.80000000, + /// -0.33333333, + /// 0.33333333, + /// 0.66666667, + /// -0.60000000, + /// 0.60000000, + /// 0.40000000, + /// 0.20000000, + /// 0.80000000, + /// 0.20000000, + /// 0.60000000, + /// 0.40000000, + /// 0.60000000, + /// 0.33333333, + /// 0.66666667, + /// 0.33333333, + /// 0.20000000, + /// 0.80000000, + /// 0.20000000, + /// ]; + /// let curl = [ + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// 2.00000000, + /// -2.00000000, + /// -2.00000000, + /// ]; + /// let qref = [ + /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, + /// 0.60000000, + /// ]; + /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; + /// let b = ceed.basis_Hcurl( + /// ElemTopology::Triangle, + /// 1, + /// 3, + /// 4, + /// &interp, + /// &curl, + /// &qref, + /// &qweight, + /// )?; + /// # Ok(()) + /// # } + /// ``` + pub fn basis_Hcurl<'a>( + &self, + topo: ElemTopology, + ncomp: usize, + nnodes: usize, + nqpts: usize, + interp: &[crate::Scalar], + curl: &[crate::Scalar], + qref: &[crate::Scalar], + qweight: &[crate::Scalar], + ) -> Result> { + Basis::create_Hcurl( + self, topo, ncomp, nnodes, nqpts, interp, curl, qref, qweight, + ) + } + /// Returns a QFunction for evaluating interior (volumetric) terms /// of a weak form corresponding to the $L^2$ inner product ///