1# Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors 2# All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3# 4# SPDX-License-Identifier: BSD-2-Clause 5# 6# This file is part of CEED: http://github.com/ceed 7 8# @file 9# Test Ceed ElemRestriction functionality 10 11import os 12import libceed 13import numpy as np 14import check 15 16# ------------------------------------------------------------------------------- 17# Test creation, use, and destruction of an element restriction 18# ------------------------------------------------------------------------------- 19 20 21def test_200(ceed_resource): 22 ceed = libceed.Ceed(ceed_resource) 23 24 num_elem = 3 25 26 x = ceed.Vector(num_elem + 1) 27 a = np.arange(10, 10 + num_elem + 1, dtype=ceed.scalar_type()) 28 x.set_array(a, cmode=libceed.USE_POINTER) 29 30 ind = np.zeros(2 * num_elem, dtype="int32") 31 for i in range(num_elem): 32 ind[2 * i + 0] = i 33 ind[2 * i + 1] = i + 1 34 r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem + 1, ind, 35 cmode=libceed.USE_POINTER) 36 37 y = ceed.Vector(2 * num_elem) 38 y.set_value(0) 39 40 r.apply(x, y) 41 42 with y.array_read() as y_array: 43 for i in range(2 * num_elem): 44 assert 10 + (i + 1) // 2 == y_array[i] 45 46# ------------------------------------------------------------------------------- 47# Test creation, use, and destruction of a strided element restriction 48# ------------------------------------------------------------------------------- 49 50 51def test_201(ceed_resource): 52 ceed = libceed.Ceed(ceed_resource) 53 54 num_elem = 3 55 56 x = ceed.Vector(2 * num_elem) 57 a = np.arange(10, 10 + 2 * num_elem, dtype=ceed.scalar_type()) 58 x.set_array(a, cmode=libceed.USE_POINTER) 59 60 strides = np.array([1, 2, 2], dtype="int32") 61 r = ceed.StridedElemRestriction(num_elem, 2, 1, 2 * num_elem, strides) 62 63 y = ceed.Vector(2 * num_elem) 64 y.set_value(0) 65 66 r.apply(x, y) 67 68 with y.array_read() as y_array: 69 for i in range(2 * num_elem): 70 assert 10 + i == y_array[i] 71 72# ------------------------------------------------------------------------------- 73# Test creation and destruction of a blocked element restriction 74# ------------------------------------------------------------------------------- 75 76 77def test_202(ceed_resource, capsys): 78 ceed = libceed.Ceed(ceed_resource) 79 80 num_elem = 8 81 elem_size = 2 82 num_blk = 2 83 blk_size = 5 84 85 x = ceed.Vector(num_elem + 1) 86 a = np.arange(10, 10 + num_elem + 1, dtype=ceed.scalar_type()) 87 x.set_array(a, cmode=libceed.COPY_VALUES) 88 89 ind = np.zeros(2 * num_elem, dtype="int32") 90 for i in range(num_elem): 91 ind[2 * i + 0] = i 92 ind[2 * i + 1] = i + 1 93 r = ceed.BlockedElemRestriction(num_elem, elem_size, blk_size, 1, 1, 94 num_elem + 1, ind, cmode=libceed.COPY_VALUES) 95 96 y = ceed.Vector(num_blk * blk_size * elem_size) 97 y.set_value(0) 98 99 # NoTranspose 100 r.apply(x, y) 101 layout = r.get_layout() 102 with y.array_read() as y_array: 103 for i in range(elem_size): 104 for j in range(1): 105 for k in range(num_elem): 106 block = int(k / blk_size) 107 elem = k % blk_size 108 indx = (i * blk_size + elem) * \ 109 layout[0] + j * blk_size * layout[1] + \ 110 block * blk_size * layout[2] 111 assert y_array[indx] == a[ind[k * elem_size + i]] 112 113 x.set_value(0) 114 r.T.apply(y, x) 115 with x.array_read() as x_array: 116 for i in range(num_elem + 1): 117 assert x_array[i] == (10 + i) * (2.0 if i > 118 0 and i < num_elem else 1.0) 119 120# ------------------------------------------------------------------------------- 121# Test creation, use, and destruction of a blocked element restriction 122# ------------------------------------------------------------------------------- 123 124 125def test_208(ceed_resource): 126 ceed = libceed.Ceed(ceed_resource) 127 128 num_elem = 8 129 elem_size = 2 130 num_blk = 2 131 blk_size = 5 132 133 x = ceed.Vector(num_elem + 1) 134 a = np.arange(10, 10 + num_elem + 1, dtype=ceed.scalar_type()) 135 x.set_array(a, cmode=libceed.COPY_VALUES) 136 137 ind = np.zeros(2 * num_elem, dtype="int32") 138 for i in range(num_elem): 139 ind[2 * i + 0] = i 140 ind[2 * i + 1] = i + 1 141 r = ceed.BlockedElemRestriction(num_elem, elem_size, blk_size, 1, 1, 142 num_elem + 1, ind, cmode=libceed.COPY_VALUES) 143 144 y = ceed.Vector(blk_size * elem_size) 145 y.set_value(0) 146 147 # NoTranspose 148 r.apply_block(1, x, y) 149 layout = r.get_layout() 150 with y.array_read() as y_array: 151 for i in range(elem_size): 152 for j in range(1): 153 for k in range(blk_size, num_elem): 154 block = int(k / blk_size) 155 elem = k % blk_size 156 indx = (i * blk_size + elem) * layout[0] + j * blk_size * \ 157 layout[1] + block * blk_size * \ 158 layout[2] - blk_size * elem_size 159 assert y_array[indx] == a[ind[k * elem_size + i]] 160 161 x.set_value(0) 162 r.T.apply_block(1, y, x) 163 with x.array_read() as x_array: 164 for i in range(blk_size, num_elem + 1): 165 assert x_array[i] == (10 + i) * (2.0 if i > 166 blk_size and i < num_elem else 1.0) 167 168# ------------------------------------------------------------------------------- 169# Test getting the multiplicity of the indices in an element restriction 170# ------------------------------------------------------------------------------- 171 172 173def test_209(ceed_resource): 174 ceed = libceed.Ceed(ceed_resource) 175 176 num_elem = 3 177 178 ind = np.zeros(4 * num_elem, dtype="int32") 179 for i in range(num_elem): 180 ind[4 * i + 0] = i * 3 + 0 181 ind[4 * i + 1] = i * 3 + 1 182 ind[4 * i + 2] = i * 3 + 2 183 ind[4 * i + 3] = i * 3 + 3 184 r = ceed.ElemRestriction(num_elem, 4, 1, 1, 3 * num_elem + 1, ind, 185 cmode=libceed.USE_POINTER) 186 187 mult = r.get_multiplicity() 188 189 with mult.array_read() as mult_array: 190 for i in range(3 * num_elem + 1): 191 val = 1 + (1 if (i > 0 and i < 3 * num_elem and i % 3 == 0) else 0) 192 assert val == mult_array[i] 193 194# ------------------------------------------------------------------------------- 195# Test creation and view of an element restriction 196# ------------------------------------------------------------------------------- 197 198 199def test_210(ceed_resource, capsys): 200 ceed = libceed.Ceed(ceed_resource) 201 202 num_elem = 3 203 204 ind = np.zeros(2 * num_elem, dtype="int32") 205 for i in range(num_elem): 206 ind[2 * i + 0] = i + 0 207 ind[2 * i + 1] = i + 1 208 r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem + 1, ind, 209 cmode=libceed.USE_POINTER) 210 211 print(r) 212 213 stdout, stderr, ref_stdout = check.output(capsys) 214 assert not stderr 215 assert stdout == ref_stdout 216 217# ------------------------------------------------------------------------------- 218# Test creation and view of a strided element restriction 219# ------------------------------------------------------------------------------- 220 221 222def test_211(ceed_resource, capsys): 223 ceed = libceed.Ceed(ceed_resource) 224 225 num_elem = 3 226 227 strides = np.array([1, 2, 2], dtype="int32") 228 r = ceed.StridedElemRestriction(num_elem, 2, 1, num_elem * 2, strides) 229 230 print(r) 231 232 stdout, stderr, ref_stdout = check.output(capsys) 233 assert not stderr 234 assert stdout == ref_stdout 235 236# ------------------------------------------------------------------------------- 237# Test creation and view of a blocked strided element restriction 238# ------------------------------------------------------------------------------- 239 240 241def test_212(ceed_resource, capsys): 242 ceed = libceed.Ceed(ceed_resource) 243 244 num_elem = 3 245 246 strides = np.array([1, 2, 2], dtype="int32") 247 r = ceed.BlockedStridedElemRestriction( 248 num_elem, 2, 2, 1, num_elem * 2, strides) 249 250 print(r) 251 252 stdout, stderr, ref_stdout = check.output(capsys) 253 assert not stderr 254 assert stdout == ref_stdout 255 256# ------------------------------------------------------------------------------- 257# Test creation, use, and destruction of an oriented element restriction 258# ------------------------------------------------------------------------------- 259 260 261def test_213(ceed_resource): 262 ceed = libceed.Ceed(ceed_resource) 263 264 num_elem = 3 265 266 x = ceed.Vector(num_elem + 1) 267 a = np.arange(10, 10 + num_elem + 1, dtype=ceed.scalar_type()) 268 x.set_array(a, cmode=libceed.USE_POINTER) 269 270 ind = np.zeros(2 * num_elem, dtype="int32") 271 orients = np.zeros(2 * num_elem, dtype="bool") 272 for i in range(num_elem): 273 ind[2 * i + 0] = i 274 ind[2 * i + 1] = i + 1 275 # flip the dofs on element 1, 3, ... 276 orients[2 * i + 0] = (i % 2) > 0 277 orients[2 * i + 1] = (i % 2) > 0 278 r = ceed.OrientedElemRestriction( 279 num_elem, 280 2, 281 1, 282 1, 283 num_elem + 1, 284 ind, 285 orients, 286 cmode=libceed.USE_POINTER) 287 288 y = ceed.Vector(2 * num_elem) 289 y.set_value(0) 290 291 r.apply(x, y) 292 293 with y.array_read() as y_array: 294 for i in range(num_elem): 295 for j in range(2): 296 k = j + 2 * i 297 assert 10 + (k + 1) // 2 == y_array[k] * pow(-1, i % 2) 298 299# ------------------------------------------------------------------------------- 300# Test creation, use, and destruction of a curl-oriented element restriction 301# ------------------------------------------------------------------------------- 302 303 304def test_214(ceed_resource): 305 ceed = libceed.Ceed(ceed_resource) 306 307 num_elem = 3 308 309 x = ceed.Vector(num_elem + 1) 310 a = np.arange(10, 10 + num_elem + 1, dtype=ceed.scalar_type()) 311 x.set_array(a, cmode=libceed.USE_POINTER) 312 313 ind = np.zeros(2 * num_elem, dtype="int32") 314 curl_orients = np.zeros(3 * 2 * num_elem, dtype="int8") 315 for i in range(num_elem): 316 ind[2 * i + 0] = i 317 ind[2 * i + 1] = i + 1 318 curl_orients[3 * 2 * i] = curl_orients[3 * 2 * (i + 1) - 1] = 0 319 if i % 2 > 0: 320 # T = [0 -1] 321 # [-1 0] 322 curl_orients[3 * 2 * i + 1] = 0 323 curl_orients[3 * 2 * i + 2] = -1 324 curl_orients[3 * 2 * i + 3] = -1 325 curl_orients[3 * 2 * i + 4] = 0 326 else: 327 # T = I 328 curl_orients[3 * 2 * i + 1] = 1 329 curl_orients[3 * 2 * i + 2] = 0 330 curl_orients[3 * 2 * i + 3] = 0 331 curl_orients[3 * 2 * i + 4] = 1 332 r = ceed.CurlOrientedElemRestriction( 333 num_elem, 334 2, 335 1, 336 1, 337 num_elem + 1, 338 ind, 339 curl_orients, 340 cmode=libceed.USE_POINTER) 341 342 y = ceed.Vector(2 * num_elem) 343 y.set_value(0) 344 345 r.apply(x, y) 346 347 with y.array_read() as y_array: 348 for i in range(num_elem): 349 for j in range(2): 350 k = j + 2 * i 351 if i % 2 > 0: 352 assert j != 0 or 10 + i + 1 == -y_array[k] 353 assert j != 1 or 10 + i == -y_array[k] 354 else: 355 assert 10 + (k + 1) // 2 == y_array[k] 356 357# ------------------------------------------------------------------------------- 358