diff --git a/python/build_ceed_cffi.py b/python/build_ceed_cffi.py index 6ffcfe0137..e4c6b99122 100644 --- a/python/build_ceed_cffi.py +++ b/python/build_ceed_cffi.py @@ -36,6 +36,7 @@ header = '\n'.join(lines) header = header.split("static inline CeedInt CeedIntPow", 1)[0] header += '\nextern int CeedVectorGetState(CeedVector, uint64_t*);' + header += '\nextern int CeedElemRestrictionGetELayout(CeedElemRestriction, CeedInt *layout);' # Note: cffi cannot handle vargs header = re.sub("va_list", "const char *", header) ffibuilder.cdef(header) diff --git a/python/ceed_elemrestriction.py b/python/ceed_elemrestriction.py index d1ae83fbe4..a08f461394 100644 --- a/python/ceed_elemrestriction.py +++ b/python/ceed_elemrestriction.py @@ -131,6 +131,29 @@ def get_multiplicity(self): # Return return mult + # Get ElemRestrition Layout + def get_layout(self): + """Get the element vector layout of an ElemRestriction. + + Returns: + layout: Vector containing layout array, stored as [nodes, components, elements]. + The data for node i, component j, element k in the element + vector is given by i*layout[0] + j*layout[1] + k*layout[2].""" + + # Create output array + layout = np.zeros(3, dtype="int32") + array_pointer = ffi.cast( + "CeedInt *", + layout.__array_interface__['data'][0]) + + # libCEED call + err_code = lib.CeedElemRestrictionGetELayout( + self._pointer[0], array_pointer) + self._ceed._check_error(err_code) + + # Return + return layout + # ------------------------------------------------------------------------------ diff --git a/python/tests/output/test_202.out b/python/tests/output/test_202.out deleted file mode 100644 index 7e9ed6855d..0000000000 --- a/python/tests/output/test_202.out +++ /dev/null @@ -1,33 +0,0 @@ -CeedVector length 20 - 10.000000 - 11.000000 - 12.000000 - 13.000000 - 14.000000 - 11.000000 - 12.000000 - 13.000000 - 14.000000 - 15.000000 - 15.000000 - 16.000000 - 17.000000 - 17.000000 - 17.000000 - 16.000000 - 17.000000 - 18.000000 - 18.000000 - 18.000000 - -CeedVector length 9 - 10.000000 - 22.000000 - 24.000000 - 26.000000 - 28.000000 - 30.000000 - 32.000000 - 34.000000 - 18.000000 - diff --git a/python/tests/output/test_208.out b/python/tests/output/test_208.out deleted file mode 100644 index e40943e09b..0000000000 --- a/python/tests/output/test_208.out +++ /dev/null @@ -1,23 +0,0 @@ -CeedVector length 10 - 15.000000 - 16.000000 - 17.000000 - 17.000000 - 17.000000 - 16.000000 - 17.000000 - 18.000000 - 18.000000 - 18.000000 - -CeedVector length 9 - 0.000000 - 0.000000 - 0.000000 - 0.000000 - 0.000000 - 15.000000 - 32.000000 - 34.000000 - 18.000000 - diff --git a/python/tests/output/test_301.out b/python/tests/output/test_301.out deleted file mode 100644 index 74a87f6bc9..0000000000 --- a/python/tests/output/test_301.out +++ /dev/null @@ -1,15 +0,0 @@ - -2.00000000 - -3.00000000 - -2.00000000 - 0.33333333 - -5.00000000 - 2.00000000 - 0.33333333 - 0.40000000 - -4.00000000 - 0.33333333 - -0.20000000 - -0.50000000 - 1.50000000 - 1.66666667 - 1.60000000 diff --git a/python/tests/output/test_304.out b/python/tests/output/test_304.out deleted file mode 100644 index 6361fddd1e..0000000000 --- a/python/tests/output/test_304.out +++ /dev/null @@ -1,22 +0,0 @@ -Q: - 0.69153918 - -0.70710678 - 0.14755868 - 0.00004347 - -0.14721060 - 0.00047835 - 0.69240821 - -0.70632831 - 0.14790715 - 0.00047882 - -0.69066919 - -0.70788369 - -0.69153892 - -0.70710646 - -0.14755804 - -0.00100064 -lambda: - 0.13487964 - 0.23325338 - 0.86513545 - 1.16666509 diff --git a/python/tests/test-2-elemrestriction.py b/python/tests/test-2-elemrestriction.py index 95de43d5bf..d3ccbe1dc6 100644 --- a/python/tests/test-2-elemrestriction.py +++ b/python/tests/test-2-elemrestriction.py @@ -30,26 +30,26 @@ def test_200(ceed_resource): ceed = libceed.Ceed(ceed_resource) - ne = 3 + num_elem = 3 - x = ceed.Vector(ne + 1) - a = np.arange(10, 10 + ne + 1, dtype="float64") + x = ceed.Vector(num_elem + 1) + a = np.arange(10, 10 + num_elem + 1, dtype="float64") x.set_array(a, cmode=libceed.USE_POINTER) - ind = np.zeros(2 * ne, dtype="int32") - for i in range(ne): + ind = np.zeros(2 * num_elem, dtype="int32") + for i in range(num_elem): ind[2 * i + 0] = i ind[2 * i + 1] = i + 1 - r = ceed.ElemRestriction(ne, 2, 1, 1, ne + 1, ind, + r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem + 1, ind, cmode=libceed.USE_POINTER) - y = ceed.Vector(2 * ne) + y = ceed.Vector(2 * num_elem) y.set_value(0) r.apply(x, y) with y.array_read() as y_array: - for i in range(2 * ne): + for i in range(2 * num_elem): assert 10 + (i + 1) // 2 == y_array[i] # ------------------------------------------------------------------------------- @@ -60,22 +60,22 @@ def test_200(ceed_resource): def test_201(ceed_resource): ceed = libceed.Ceed(ceed_resource) - ne = 3 + num_elem = 3 - x = ceed.Vector(2 * ne) - a = np.arange(10, 10 + 2 * ne, dtype="float64") + x = ceed.Vector(2 * num_elem) + a = np.arange(10, 10 + 2 * num_elem, dtype="float64") x.set_array(a, cmode=libceed.USE_POINTER) strides = np.array([1, 2, 2], dtype="int32") - r = ceed.StridedElemRestriction(ne, 2, 1, 2 * ne, strides) + r = ceed.StridedElemRestriction(num_elem, 2, 1, 2 * num_elem, strides) - y = ceed.Vector(2 * ne) + y = ceed.Vector(2 * num_elem) y.set_value(0) r.apply(x, y) with y.array_read() as y_array: - for i in range(2 * ne): + for i in range(2 * num_elem): assert 10 + i == y_array[i] # ------------------------------------------------------------------------------- @@ -86,71 +86,93 @@ def test_201(ceed_resource): def test_202(ceed_resource, capsys): ceed = libceed.Ceed(ceed_resource) - ne = 8 - blksize = 5 + num_elem = 8 + elem_size = 2 + num_blk = 2 + blk_size = 5 - x = ceed.Vector(ne + 1) - a = np.arange(10, 10 + ne + 1, dtype="float64") - x.set_array(a, cmode=libceed.USE_POINTER) + x = ceed.Vector(num_elem + 1) + a = np.arange(10, 10 + num_elem + 1, dtype="float64") + x.set_array(a, cmode=libceed.COPY_VALUES) - ind = np.zeros(2 * ne, dtype="int32") - for i in range(ne): + ind = np.zeros(2 * num_elem, dtype="int32") + for i in range(num_elem): ind[2 * i + 0] = i ind[2 * i + 1] = i + 1 - r = ceed.BlockedElemRestriction(ne, 2, blksize, 1, 1, ne + 1, ind, - cmode=libceed.USE_POINTER) + r = ceed.BlockedElemRestriction(num_elem, elem_size, blk_size, 1, 1, + num_elem + 1, ind, cmode=libceed.COPY_VALUES) - y = ceed.Vector(2 * blksize * 2) + y = ceed.Vector(num_blk * blk_size * elem_size) y.set_value(0) + # NoTranspose r.apply(x, y) - - print(y) + layout = r.get_layout() + with y.array_read() as y_array: + for i in range(elem_size): + for j in range(1): + for k in range(num_elem): + block = int(k / blk_size) + elem = k % blk_size + indx = (i * blk_size + elem) * \ + layout[0] + j * blk_size * layout[1] + \ + block * blk_size * layout[2] + assert y_array[indx] == a[ind[k * elem_size + i]] x.set_value(0) r.T.apply(y, x) - print(x) - - stdout, stderr, ref_stdout = check.output(capsys) - assert not stderr - assert stdout == ref_stdout + with x.array_read() as x_array: + for i in range(num_elem + 1): + assert x_array[i] == (10 + i) * (2.0 if i > + 0 and i < num_elem else 1.0) # ------------------------------------------------------------------------------- # Test creation, use, and destruction of a blocked element restriction # ------------------------------------------------------------------------------- -def test_208(ceed_resource, capsys): +def test_208(ceed_resource): ceed = libceed.Ceed(ceed_resource) - ne = 8 - blksize = 5 + num_elem = 8 + elem_size = 2 + num_blk = 2 + blk_size = 5 - x = ceed.Vector(ne + 1) - a = np.arange(10, 10 + ne + 1, dtype="float64") - x.set_array(a, cmode=libceed.USE_POINTER) + x = ceed.Vector(num_elem + 1) + a = np.arange(10, 10 + num_elem + 1, dtype="float64") + x.set_array(a, cmode=libceed.COPY_VALUES) - ind = np.zeros(2 * ne, dtype="int32") - for i in range(ne): + ind = np.zeros(2 * num_elem, dtype="int32") + for i in range(num_elem): ind[2 * i + 0] = i ind[2 * i + 1] = i + 1 - r = ceed.BlockedElemRestriction(ne, 2, blksize, 1, 1, ne + 1, ind, - cmode=libceed.USE_POINTER) + r = ceed.BlockedElemRestriction(num_elem, elem_size, blk_size, 1, 1, + num_elem + 1, ind, cmode=libceed.COPY_VALUES) - y = ceed.Vector(blksize * 2) + y = ceed.Vector(blk_size * elem_size) y.set_value(0) + # NoTranspose r.apply_block(1, x, y) - - print(y) + layout = r.get_layout() + with y.array_read() as y_array: + for i in range(elem_size): + for j in range(1): + for k in range(blk_size, num_elem): + block = int(k / blk_size) + elem = k % blk_size + indx = (i * blk_size + elem) * layout[0] + j * blk_size * \ + layout[1] + block * blk_size * \ + layout[2] - blk_size * elem_size + assert y_array[indx] == a[ind[k * elem_size + i]] x.set_value(0) r.T.apply_block(1, y, x) - print(x) - - stdout, stderr, ref_stdout = check.output(capsys) - assert not stderr - assert stdout == ref_stdout + with x.array_read() as x_array: + for i in range(blk_size, num_elem + 1): + assert x_array[i] == (10 + i) * (2.0 if i > + blk_size and i < num_elem else 1.0) # ------------------------------------------------------------------------------- # Test getting the multiplicity of the indices in an element restriction @@ -160,22 +182,22 @@ def test_208(ceed_resource, capsys): def test_209(ceed_resource): ceed = libceed.Ceed(ceed_resource) - ne = 3 + num_elem = 3 - ind = np.zeros(4 * ne, dtype="int32") - for i in range(ne): + ind = np.zeros(4 * num_elem, dtype="int32") + for i in range(num_elem): ind[4 * i + 0] = i * 3 + 0 ind[4 * i + 1] = i * 3 + 1 ind[4 * i + 2] = i * 3 + 2 ind[4 * i + 3] = i * 3 + 3 - r = ceed.ElemRestriction(ne, 4, 1, 1, 3 * ne + 1, ind, + r = ceed.ElemRestriction(num_elem, 4, 1, 1, 3 * num_elem + 1, ind, cmode=libceed.USE_POINTER) mult = r.get_multiplicity() with mult.array_read() as mult_array: - for i in range(3 * ne + 1): - val = 1 + (1 if (i > 0 and i < 3 * ne and i % 3 == 0) else 0) + for i in range(3 * num_elem + 1): + val = 1 + (1 if (i > 0 and i < 3 * num_elem and i % 3 == 0) else 0) assert val == mult_array[i] # ------------------------------------------------------------------------------- @@ -186,13 +208,13 @@ def test_209(ceed_resource): def test_210(ceed_resource, capsys): ceed = libceed.Ceed(ceed_resource) - ne = 3 + num_elem = 3 - ind = np.zeros(2 * ne, dtype="int32") - for i in range(ne): + ind = np.zeros(2 * num_elem, dtype="int32") + for i in range(num_elem): ind[2 * i + 0] = i + 0 ind[2 * i + 1] = i + 1 - r = ceed.ElemRestriction(ne, 2, 1, 1, ne + 1, ind, + r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem + 1, ind, cmode=libceed.USE_POINTER) print(r) @@ -209,10 +231,10 @@ def test_210(ceed_resource, capsys): def test_211(ceed_resource, capsys): ceed = libceed.Ceed(ceed_resource) - ne = 3 + num_elem = 3 strides = np.array([1, 2, 2], dtype="int32") - r = ceed.StridedElemRestriction(ne, 2, 1, ne + 1, strides) + r = ceed.StridedElemRestriction(num_elem, 2, 1, num_elem + 1, strides) print(r) @@ -228,10 +250,11 @@ def test_211(ceed_resource, capsys): def test_212(ceed_resource, capsys): ceed = libceed.Ceed(ceed_resource) - ne = 3 + num_elem = 3 strides = np.array([1, 2, 2], dtype="int32") - r = ceed.BlockedStridedElemRestriction(ne, 2, 2, 1, ne + 1, strides) + r = ceed.BlockedStridedElemRestriction( + num_elem, 2, 2, 1, num_elem + 1, strides) print(r) diff --git a/python/tests/test-3-basis.py b/python/tests/test-3-basis.py index a9366f74ce..0b83aeefe0 100644 --- a/python/tests/test-3-basis.py +++ b/python/tests/test-3-basis.py @@ -71,34 +71,32 @@ def test_300(ceed_resource, capsys): # ------------------------------------------------------------------------------- -def test_301(ceed_resource, capsys): +def test_301(ceed_resource): ceed = libceed.Ceed(ceed_resource) + m = 4 + n = 3 + a = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype="float64") qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype="float64") tau = np.empty(3, dtype="float64") - qr, tau = ceed.qr_factorization(qr, tau, 4, 3) + qr, tau = ceed.qr_factorization(qr, tau, m, n) + np_qr, np_tau = np.linalg.qr(a.reshape(m, n), mode="raw") - for i in range(len(qr)): - if qr[i] <= TOL and qr[i] >= -TOL: - qr[i] = 0 - print("%12.8f" % qr[i]) + for i in range(n): + assert tau[i] == np_tau[i] - for i in range(len(tau)): - if tau[i] <= TOL and tau[i] >= -TOL: - tau[i] = 0 - print("%12.8f" % tau[i]) - - stdout, stderr, ref_stdout = check.output(capsys) - assert not stderr - assert stdout == ref_stdout + qr = qr.reshape(m, n) + for i in range(m): + for j in range(n): + assert round(qr[i, j] - np_qr[j, i], 10) == 0 # ------------------------------------------------------------------------------- # Test Symmetric Schur Decomposition # ------------------------------------------------------------------------------- -def test_304(ceed_resource, capsys): +def test_304(ceed_resource): ceed = libceed.Ceed(ceed_resource) A = np.array([0.2, 0.0745355993, -0.0745355993, 0.0333333333, @@ -120,7 +118,7 @@ def test_304(ceed_resource, capsys): # ------------------------------------------------------------------------------- -def test_305(ceed_resource, capsys): +def test_305(ceed_resource): ceed = libceed.Ceed(ceed_resource) M = np.array([0.2, 0.0745355993, -0.0745355993, 0.0333333333,