1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3# reserved. See files LICENSE and NOTICE for details. 4# 5# This file is part of CEED, a collection of benchmarks, miniapps, software 6# libraries and APIs for efficient high-order finite element and spectral 7# element discretizations for exascale applications. For more information and 8# source code availability see http://github.com/ceed. 9# 10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11# a collaborative effort of two U.S. Department of Energy organizations (Office 12# of Science and the National Nuclear Security Administration) responsible for 13# the planning and preparation of a capable exascale ecosystem, including 14# software, applications, hardware, advanced system engineering and early 15# testbed platforms, in support of the nation's exascale computing imperative. 16 17# @file 18# Test Ceed Operator functionality 19 20import os 21import libceed 22import numpy as np 23import check 24import buildmats as bm 25 26TOL = np.finfo(float).eps * 256 27 28# ------------------------------------------------------------------------------- 29# Utility 30# ------------------------------------------------------------------------------- 31 32 33def load_qfs_so(): 34 from distutils.sysconfig import get_config_var 35 import ctypes 36 37 file_dir = os.path.dirname(os.path.abspath(__file__)) 38 qfs_so = os.path.join( 39 file_dir, 40 "libceed_qfunctions" + get_config_var("EXT_SUFFIX")) 41 42 # Load library 43 return ctypes.cdll.LoadLibrary(qfs_so) 44 45# ------------------------------------------------------------------------------- 46# Test creation, action, and destruction for mass matrix operator 47# ------------------------------------------------------------------------------- 48 49 50def test_500(ceed_resource): 51 ceed = libceed.Ceed(ceed_resource) 52 53 nelem = 15 54 p = 5 55 q = 8 56 nx = nelem + 1 57 nu = nelem * (p - 1) + 1 58 59 # Vectors 60 x = ceed.Vector(nx) 61 x_array = np.zeros(nx) 62 for i in range(nx): 63 x_array[i] = i / (nx - 1.0) 64 x.set_array(x_array, cmode=libceed.USE_POINTER) 65 66 qdata = ceed.Vector(nelem * q) 67 u = ceed.Vector(nu) 68 v = ceed.Vector(nu) 69 70 # Restrictions 71 indx = np.zeros(nx * 2, dtype="int32") 72 for i in range(nx): 73 indx[2 * i + 0] = i 74 indx[2 * i + 1] = i + 1 75 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 76 cmode=libceed.USE_POINTER) 77 78 indu = np.zeros(nelem * p, dtype="int32") 79 for i in range(nelem): 80 for j in range(p): 81 indu[p * i + j] = i * (p - 1) + j 82 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 83 cmode=libceed.USE_POINTER) 84 strides = np.array([1, q, q], dtype="int32") 85 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 86 87 # Bases 88 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 89 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 90 91 # QFunctions 92 file_dir = os.path.dirname(os.path.abspath(__file__)) 93 qfs = load_qfs_so() 94 95 qf_setup = ceed.QFunction(1, qfs.setup_mass, 96 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 97 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 98 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 99 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 100 101 qf_mass = ceed.QFunction(1, qfs.apply_mass, 102 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 103 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 104 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 105 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 106 107 # Operators 108 op_setup = ceed.Operator(qf_setup) 109 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 110 libceed.VECTOR_NONE) 111 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 112 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 113 libceed.VECTOR_ACTIVE) 114 115 op_mass = ceed.Operator(qf_mass) 116 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 117 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 118 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 119 120 # Setup 121 op_setup.apply(x, qdata) 122 123 # Apply mass matrix 124 u.set_value(0) 125 op_mass.apply(u, v) 126 127 # Check 128 with v.array_read() as v_array: 129 for i in range(q): 130 assert abs(v_array[i]) < TOL 131 132# ------------------------------------------------------------------------------- 133# Test creation, action, and destruction for mass matrix operator 134# ------------------------------------------------------------------------------- 135 136 137def test_501(ceed_resource): 138 ceed = libceed.Ceed(ceed_resource) 139 140 nelem = 15 141 p = 5 142 q = 8 143 nx = nelem + 1 144 nu = nelem * (p - 1) + 1 145 146 # Vectors 147 x = ceed.Vector(nx) 148 x_array = np.zeros(nx) 149 for i in range(nx): 150 x_array[i] = i / (nx - 1.0) 151 x.set_array(x_array, cmode=libceed.USE_POINTER) 152 153 qdata = ceed.Vector(nelem * q) 154 u = ceed.Vector(nu) 155 v = ceed.Vector(nu) 156 157 # Restrictions 158 indx = np.zeros(nx * 2, dtype="int32") 159 for i in range(nx): 160 indx[2 * i + 0] = i 161 indx[2 * i + 1] = i + 1 162 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 163 cmode=libceed.USE_POINTER) 164 165 indu = np.zeros(nelem * p, dtype="int32") 166 for i in range(nelem): 167 for j in range(p): 168 indu[p * i + j] = i * (p - 1) + j 169 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 170 cmode=libceed.USE_POINTER) 171 strides = np.array([1, q, q], dtype="int32") 172 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 173 174 # Bases 175 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 176 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 177 178 # QFunctions 179 file_dir = os.path.dirname(os.path.abspath(__file__)) 180 qfs = load_qfs_so() 181 182 qf_setup = ceed.QFunction(1, qfs.setup_mass, 183 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 184 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 185 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 186 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 187 188 qf_mass = ceed.QFunction(1, qfs.apply_mass, 189 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 190 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 191 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 192 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 193 194 # Operators 195 op_setup = ceed.Operator(qf_setup) 196 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 197 libceed.VECTOR_NONE) 198 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 199 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 200 libceed.VECTOR_ACTIVE) 201 202 op_mass = ceed.Operator(qf_mass) 203 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 204 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 205 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 206 207 # Setup 208 op_setup.apply(x, qdata) 209 210 # Apply mass matrix 211 u.set_value(1.) 212 op_mass.apply(u, v) 213 214 # Check 215 with v.array_read() as v_array: 216 total = 0.0 217 for i in range(nu): 218 total = total + v_array[i] 219 assert abs(total - 1.0) < TOL 220 221# ------------------------------------------------------------------------------- 222# Test creation, action, and destruction for mass matrix operator with multiple 223# components 224# ------------------------------------------------------------------------------- 225 226 227def test_502(ceed_resource): 228 ceed = libceed.Ceed(ceed_resource) 229 230 nelem = 15 231 p = 5 232 q = 8 233 nx = nelem + 1 234 nu = nelem * (p - 1) + 1 235 236 # Vectors 237 x = ceed.Vector(nx) 238 x_array = np.zeros(nx) 239 for i in range(nx): 240 x_array[i] = i / (nx - 1.0) 241 x.set_array(x_array, cmode=libceed.USE_POINTER) 242 243 qdata = ceed.Vector(nelem * q) 244 u = ceed.Vector(2 * nu) 245 v = ceed.Vector(2 * nu) 246 247 # Restrictions 248 indx = np.zeros(nx * 2, dtype="int32") 249 for i in range(nx): 250 indx[2 * i + 0] = i 251 indx[2 * i + 1] = i + 1 252 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 253 cmode=libceed.USE_POINTER) 254 255 indu = np.zeros(nelem * p, dtype="int32") 256 for i in range(nelem): 257 for j in range(p): 258 indu[p * i + j] = 2 * (i * (p - 1) + j) 259 ru = ceed.ElemRestriction(nelem, p, 2, 1, 2 * nu, indu, 260 cmode=libceed.USE_POINTER) 261 strides = np.array([1, q, q], dtype="int32") 262 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 263 264 # Bases 265 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 266 bu = ceed.BasisTensorH1Lagrange(1, 2, p, q, libceed.GAUSS) 267 268 # QFunctions 269 file_dir = os.path.dirname(os.path.abspath(__file__)) 270 qfs = load_qfs_so() 271 272 qf_setup = ceed.QFunction(1, qfs.setup_mass, 273 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 274 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 275 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 276 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 277 278 qf_mass = ceed.QFunction(1, qfs.apply_mass_two, 279 os.path.join(file_dir, "test-qfunctions.h:apply_mass_two")) 280 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 281 qf_mass.add_input("u", 2, libceed.EVAL_INTERP) 282 qf_mass.add_output("v", 2, libceed.EVAL_INTERP) 283 284 # Operators 285 op_setup = ceed.Operator(qf_setup) 286 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 287 libceed.VECTOR_NONE) 288 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 289 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 290 libceed.VECTOR_ACTIVE) 291 292 op_mass = ceed.Operator(qf_mass) 293 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 294 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 295 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 296 297 # Setup 298 op_setup.apply(x, qdata) 299 300 # Apply mass matrix 301 with u.array() as u_array: 302 for i in range(nu): 303 u_array[2 * i] = 1. 304 u_array[2 * i + 1] = 2. 305 op_mass.apply(u, v) 306 307 # Check 308 with v.array_read() as v_array: 309 total_1 = 0.0 310 total_2 = 0.0 311 for i in range(nu): 312 total_1 = total_1 + v_array[2 * i] 313 total_2 = total_2 + v_array[2 * i + 1] 314 assert abs(total_1 - 1.0) < 1E-13 315 assert abs(total_2 - 2.0) < 1E-13 316 317# ------------------------------------------------------------------------------- 318# Test creation, action, and destruction for mass matrix operator with passive 319# inputs and outputs 320# ------------------------------------------------------------------------------- 321 322 323def test_503(ceed_resource): 324 ceed = libceed.Ceed(ceed_resource) 325 326 nelem = 15 327 p = 5 328 q = 8 329 nx = nelem + 1 330 nu = nelem * (p - 1) + 1 331 332 # Vectors 333 x = ceed.Vector(nx) 334 x_array = np.zeros(nx) 335 for i in range(nx): 336 x_array[i] = i / (nx - 1.0) 337 x.set_array(x_array, cmode=libceed.USE_POINTER) 338 339 qdata = ceed.Vector(nelem * q) 340 u = ceed.Vector(nu) 341 v = ceed.Vector(nu) 342 343 # Restrictions 344 indx = np.zeros(nx * 2, dtype="int32") 345 for i in range(nx): 346 indx[2 * i + 0] = i 347 indx[2 * i + 1] = i + 1 348 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 349 cmode=libceed.USE_POINTER) 350 351 indu = np.zeros(nelem * p, dtype="int32") 352 for i in range(nelem): 353 for j in range(p): 354 indu[p * i + j] = i * (p - 1) + j 355 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 356 cmode=libceed.USE_POINTER) 357 strides = np.array([1, q, q], dtype="int32") 358 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 359 360 # Bases 361 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 362 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 363 364 # QFunctions 365 file_dir = os.path.dirname(os.path.abspath(__file__)) 366 qfs = load_qfs_so() 367 368 qf_setup = ceed.QFunction(1, qfs.setup_mass, 369 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 370 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 371 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 372 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 373 374 qf_mass = ceed.QFunction(1, qfs.apply_mass, 375 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 376 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 377 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 378 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 379 380 # Operators 381 op_setup = ceed.Operator(qf_setup) 382 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 383 libceed.VECTOR_NONE) 384 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 385 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 386 libceed.VECTOR_ACTIVE) 387 388 op_mass = ceed.Operator(qf_mass) 389 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 390 op_mass.set_field("u", ru, bu, u) 391 op_mass.set_field("v", ru, bu, v) 392 393 # Setup 394 op_setup.apply(x, qdata) 395 396 # Apply mass matrix 397 u.set_value(1) 398 op_mass.apply(libceed.VECTOR_NONE, libceed.VECTOR_NONE) 399 400 # Check 401 with v.array_read() as v_array: 402 total = 0.0 403 for i in range(nu): 404 total = total + v_array[i] 405 assert abs(total - 1.0) < 1E-13 406 407# ------------------------------------------------------------------------------- 408# Test viewing of mass matrix operator 409# ------------------------------------------------------------------------------- 410 411 412def test_504(ceed_resource, capsys): 413 ceed = libceed.Ceed(ceed_resource) 414 415 nelem = 15 416 p = 5 417 q = 8 418 nx = nelem + 1 419 nu = nelem * (p - 1) + 1 420 421 # Vectors 422 qdata = ceed.Vector(nelem * q) 423 424 # Restrictions 425 indx = np.zeros(nx * 2, dtype="int32") 426 for i in range(nx): 427 indx[2 * i + 0] = i 428 indx[2 * i + 1] = i + 1 429 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 430 cmode=libceed.USE_POINTER) 431 432 indu = np.zeros(nelem * p, dtype="int32") 433 for i in range(nelem): 434 for j in range(p): 435 indu[p * i + j] = i * (p - 1) + j 436 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 437 cmode=libceed.USE_POINTER) 438 strides = np.array([1, q, q], dtype="int32") 439 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 440 441 # Bases 442 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 443 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 444 445 # QFunctions 446 file_dir = os.path.dirname(os.path.abspath(__file__)) 447 qfs = load_qfs_so() 448 449 qf_setup = ceed.QFunction(1, qfs.setup_mass, 450 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 451 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 452 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 453 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 454 455 qf_mass = ceed.QFunction(1, qfs.apply_mass, 456 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 457 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 458 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 459 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 460 461 # Operators 462 op_setup = ceed.Operator(qf_setup) 463 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 464 libceed.VECTOR_NONE) 465 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 466 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 467 libceed.VECTOR_ACTIVE) 468 469 op_mass = ceed.Operator(qf_mass) 470 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 471 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 472 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 473 474 # View 475 print(op_setup) 476 print(op_mass) 477 478 stdout, stderr, ref_stdout = check.output(capsys) 479 assert not stderr 480 assert stdout == ref_stdout 481 482# ------------------------------------------------------------------------------- 483# Test CeedOperatorApplyAdd 484# ------------------------------------------------------------------------------- 485 486 487def test_505(ceed_resource): 488 ceed = libceed.Ceed(ceed_resource) 489 490 nelem = 15 491 p = 5 492 q = 8 493 nx = nelem + 1 494 nu = nelem * (p - 1) + 1 495 496 # Vectors 497 x = ceed.Vector(nx) 498 x_array = np.zeros(nx) 499 for i in range(nx): 500 x_array[i] = i / (nx - 1.0) 501 x.set_array(x_array, cmode=libceed.USE_POINTER) 502 503 qdata = ceed.Vector(nelem * q) 504 u = ceed.Vector(nu) 505 v = ceed.Vector(nu) 506 507 # Restrictions 508 indx = np.zeros(nx * 2, dtype="int32") 509 for i in range(nx): 510 indx[2 * i + 0] = i 511 indx[2 * i + 1] = i + 1 512 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 513 cmode=libceed.USE_POINTER) 514 515 indu = np.zeros(nelem * p, dtype="int32") 516 for i in range(nelem): 517 for j in range(p): 518 indu[p * i + j] = i * (p - 1) + j 519 ru = ceed.ElemRestriction(nelem, p, 1, 1, nu, indu, 520 cmode=libceed.USE_POINTER) 521 strides = np.array([1, q, q], dtype="int32") 522 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 523 524 # Bases 525 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 526 bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS) 527 528 # QFunctions 529 file_dir = os.path.dirname(os.path.abspath(__file__)) 530 qfs = load_qfs_so() 531 532 qf_setup = ceed.QFunction(1, qfs.setup_mass, 533 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 534 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 535 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 536 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 537 538 qf_mass = ceed.QFunction(1, qfs.apply_mass, 539 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 540 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 541 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 542 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 543 544 # Operators 545 op_setup = ceed.Operator(qf_setup) 546 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 547 libceed.VECTOR_NONE) 548 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 549 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 550 libceed.VECTOR_ACTIVE) 551 552 op_mass = ceed.Operator(qf_mass) 553 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 554 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 555 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 556 557 # Setup 558 op_setup.apply(x, qdata) 559 560 # Apply mass matrix with v = 0 561 u.set_value(1.) 562 v.set_value(0.) 563 op_mass.apply_add(u, v) 564 565 # Check 566 with v.array_read() as v_array: 567 total = 0.0 568 for i in range(nu): 569 total = total + v_array[i] 570 assert abs(total - 1.0) < TOL 571 572 # Apply mass matrix with v = 0 573 v.set_value(1.) 574 op_mass.apply_add(u, v) 575 576 # Check 577 with v.array_read() as v_array: 578 total = -nu 579 for i in range(nu): 580 total = total + v_array[i] 581 assert abs(total - 1.0) < 1E-10 582 583# ------------------------------------------------------------------------------- 584# Test creation, action, and destruction for mass matrix operator 585# ------------------------------------------------------------------------------- 586 587 588def test_510(ceed_resource): 589 ceed = libceed.Ceed(ceed_resource) 590 591 nelem = 12 592 dim = 2 593 p = 6 594 q = 4 595 nx, ny = 3, 2 596 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 597 nqpts = nelem * q 598 599 # Vectors 600 x = ceed.Vector(dim * ndofs) 601 x_array = np.zeros(dim * ndofs) 602 for i in range(ndofs): 603 x_array[i] = (1. / (nx * 2)) * (i % (nx * 2 + 1)) 604 x_array[i + ndofs] = (1. / (ny * 2)) * (i / (nx * 2 + 1)) 605 x.set_array(x_array, cmode=libceed.USE_POINTER) 606 607 qdata = ceed.Vector(nqpts) 608 u = ceed.Vector(ndofs) 609 v = ceed.Vector(ndofs) 610 611 # Restrictions 612 indx = np.zeros(nelem * p, dtype="int32") 613 for i in range(nelem // 2): 614 col = i % nx 615 row = i // nx 616 offset = col * 2 + row * (nx * 2 + 1) * 2 617 618 indx[i * 2 * p + 0] = 2 + offset 619 indx[i * 2 * p + 1] = 9 + offset 620 indx[i * 2 * p + 2] = 16 + offset 621 indx[i * 2 * p + 3] = 1 + offset 622 indx[i * 2 * p + 4] = 8 + offset 623 indx[i * 2 * p + 5] = 0 + offset 624 625 indx[i * 2 * p + 6] = 14 + offset 626 indx[i * 2 * p + 7] = 7 + offset 627 indx[i * 2 * p + 8] = 0 + offset 628 indx[i * 2 * p + 9] = 15 + offset 629 indx[i * 2 * p + 10] = 8 + offset 630 indx[i * 2 * p + 11] = 16 + offset 631 632 rx = ceed.ElemRestriction(nelem, p, dim, ndofs, dim * ndofs, indx, 633 cmode=libceed.USE_POINTER) 634 635 ru = ceed.ElemRestriction(nelem, p, 1, 1, ndofs, indx, 636 cmode=libceed.USE_POINTER) 637 strides = np.array([1, q, q], dtype="int32") 638 rui = ceed.StridedElemRestriction(nelem, q, 1, nqpts, strides) 639 640 # Bases 641 qref = np.empty(dim * q, dtype="float64") 642 qweight = np.empty(q, dtype="float64") 643 interp, grad = bm.buildmats(qref, qweight) 644 645 bx = ceed.BasisH1(libceed.TRIANGLE, dim, p, q, interp, grad, qref, qweight) 646 bu = ceed.BasisH1(libceed.TRIANGLE, 1, p, q, interp, grad, qref, qweight) 647 648 # QFunctions 649 file_dir = os.path.dirname(os.path.abspath(__file__)) 650 qfs = load_qfs_so() 651 652 qf_setup = ceed.QFunction(1, qfs.setup_mass_2d, 653 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 654 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 655 qf_setup.add_input("dx", dim * dim, libceed.EVAL_GRAD) 656 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 657 658 qf_mass = ceed.QFunction(1, qfs.apply_mass, 659 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 660 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 661 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 662 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 663 664 # Operators 665 op_setup = ceed.Operator(qf_setup) 666 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 667 libceed.VECTOR_NONE) 668 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 669 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 670 libceed.VECTOR_ACTIVE) 671 672 op_mass = ceed.Operator(qf_mass) 673 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 674 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 675 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 676 677 # Setup 678 op_setup.apply(x, qdata) 679 680 # Apply mass matrix 681 u.set_value(0.) 682 op_mass.apply(u, v) 683 684 # Check 685 with v.array_read() as v_array: 686 for i in range(ndofs): 687 assert abs(v_array[i]) < TOL 688 689# ------------------------------------------------------------------------------- 690# Test creation, action, and destruction for mass matrix operator 691# ------------------------------------------------------------------------------- 692 693 694def test_511(ceed_resource): 695 ceed = libceed.Ceed(ceed_resource) 696 697 nelem = 12 698 dim = 2 699 p = 6 700 q = 4 701 nx, ny = 3, 2 702 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 703 nqpts = nelem * q 704 705 # Vectors 706 x = ceed.Vector(dim * ndofs) 707 x_array = np.zeros(dim * ndofs) 708 for i in range(ndofs): 709 x_array[i] = (1. / (nx * 2)) * (i % (nx * 2 + 1)) 710 x_array[i + ndofs] = (1. / (ny * 2)) * (i / (nx * 2 + 1)) 711 x.set_array(x_array, cmode=libceed.USE_POINTER) 712 713 qdata = ceed.Vector(nqpts) 714 u = ceed.Vector(ndofs) 715 v = ceed.Vector(ndofs) 716 717 # Restrictions 718 indx = np.zeros(nelem * p, dtype="int32") 719 for i in range(nelem // 2): 720 col = i % nx 721 row = i // nx 722 offset = col * 2 + row * (nx * 2 + 1) * 2 723 724 indx[i * 2 * p + 0] = 2 + offset 725 indx[i * 2 * p + 1] = 9 + offset 726 indx[i * 2 * p + 2] = 16 + offset 727 indx[i * 2 * p + 3] = 1 + offset 728 indx[i * 2 * p + 4] = 8 + offset 729 indx[i * 2 * p + 5] = 0 + offset 730 731 indx[i * 2 * p + 6] = 14 + offset 732 indx[i * 2 * p + 7] = 7 + offset 733 indx[i * 2 * p + 8] = 0 + offset 734 indx[i * 2 * p + 9] = 15 + offset 735 indx[i * 2 * p + 10] = 8 + offset 736 indx[i * 2 * p + 11] = 16 + offset 737 738 rx = ceed.ElemRestriction(nelem, p, dim, ndofs, dim * ndofs, indx, 739 cmode=libceed.USE_POINTER) 740 741 ru = ceed.ElemRestriction(nelem, p, 1, 1, ndofs, indx, 742 cmode=libceed.USE_POINTER) 743 strides = np.array([1, q, q], dtype="int32") 744 rui = ceed.StridedElemRestriction(nelem, q, 1, nqpts, strides) 745 746 # Bases 747 qref = np.empty(dim * q, dtype="float64") 748 qweight = np.empty(q, dtype="float64") 749 interp, grad = bm.buildmats(qref, qweight) 750 751 bx = ceed.BasisH1(libceed.TRIANGLE, dim, p, q, interp, grad, qref, qweight) 752 bu = ceed.BasisH1(libceed.TRIANGLE, 1, p, q, interp, grad, qref, qweight) 753 754 # QFunctions 755 file_dir = os.path.dirname(os.path.abspath(__file__)) 756 qfs = load_qfs_so() 757 758 qf_setup = ceed.QFunction(1, qfs.setup_mass_2d, 759 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 760 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 761 qf_setup.add_input("dx", dim * dim, libceed.EVAL_GRAD) 762 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 763 764 qf_mass = ceed.QFunction(1, qfs.apply_mass, 765 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 766 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 767 qf_mass.add_input("u", 1, libceed.EVAL_INTERP) 768 qf_mass.add_output("v", 1, libceed.EVAL_INTERP) 769 770 # Operators 771 op_setup = ceed.Operator(qf_setup) 772 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 773 libceed.VECTOR_NONE) 774 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 775 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 776 libceed.VECTOR_ACTIVE) 777 778 op_mass = ceed.Operator(qf_mass) 779 op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 780 op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE) 781 op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE) 782 783 # Setup 784 op_setup.apply(x, qdata) 785 786 # Apply mass matrix 787 u.set_value(1.) 788 op_mass.apply(u, v) 789 790 # Check 791 with v.array_read() as v_array: 792 total = 0.0 793 for i in range(ndofs): 794 total = total + v_array[i] 795 assert abs(total - 1.0) < 1E-10 796 797# ------------------------------------------------------------------------------- 798# Test creation, action, and destruction for composite mass matrix operator 799# ------------------------------------------------------------------------------- 800 801 802def test_520(ceed_resource): 803 ceed = libceed.Ceed(ceed_resource) 804 805 nelem_tet, p_tet, q_tet = 6, 6, 4 806 nelem_hex, p_hex, q_hex = 6, 3, 4 807 nx, ny = 3, 3 808 dim = 2 809 nx_tet, ny_tet, nx_hex = 3, 1, 3 810 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 811 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 812 813 # Vectors 814 x = ceed.Vector(dim * ndofs) 815 x_array = np.zeros(dim * ndofs) 816 for i in range(ny * 2 + 1): 817 for j in range(nx * 2 + 1): 818 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 819 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 820 x.set_array(x_array, cmode=libceed.USE_POINTER) 821 822 qdata_hex = ceed.Vector(nqpts_hex) 823 qdata_tet = ceed.Vector(nqpts_tet) 824 u = ceed.Vector(ndofs) 825 v = ceed.Vector(ndofs) 826 827 # ------------------------- Tet Elements ------------------------- 828 829 # Restrictions 830 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 831 for i in range(nelem_tet // 2): 832 col = i % nx 833 row = i // nx 834 offset = col * 2 + row * (nx * 2 + 1) * 2 835 836 indx_tet[i * 2 * p_tet + 0] = 2 + offset 837 indx_tet[i * 2 * p_tet + 1] = 9 + offset 838 indx_tet[i * 2 * p_tet + 2] = 16 + offset 839 indx_tet[i * 2 * p_tet + 3] = 1 + offset 840 indx_tet[i * 2 * p_tet + 4] = 8 + offset 841 indx_tet[i * 2 * p_tet + 5] = 0 + offset 842 843 indx_tet[i * 2 * p_tet + 6] = 14 + offset 844 indx_tet[i * 2 * p_tet + 7] = 7 + offset 845 indx_tet[i * 2 * p_tet + 8] = 0 + offset 846 indx_tet[i * 2 * p_tet + 9] = 15 + offset 847 indx_tet[i * 2 * p_tet + 10] = 8 + offset 848 indx_tet[i * 2 * p_tet + 11] = 16 + offset 849 850 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 851 indx_tet, cmode=libceed.USE_POINTER) 852 853 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 854 cmode=libceed.USE_POINTER) 855 strides = np.array([1, q_tet, q_tet], dtype="int32") 856 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 857 strides) 858 859 # Bases 860 qref = np.empty(dim * q_tet, dtype="float64") 861 qweight = np.empty(q_tet, dtype="float64") 862 interp, grad = bm.buildmats(qref, qweight) 863 864 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 865 qweight) 866 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 867 qweight) 868 869 # QFunctions 870 file_dir = os.path.dirname(os.path.abspath(__file__)) 871 qfs = load_qfs_so() 872 873 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 874 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 875 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 876 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 877 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 878 879 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 880 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 881 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 882 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 883 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 884 885 # Operators 886 op_setup_tet = ceed.Operator(qf_setup_tet) 887 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 888 libceed.VECTOR_NONE) 889 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 890 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 891 qdata_tet) 892 893 op_mass_tet = ceed.Operator(qf_mass_tet) 894 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 895 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 896 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 897 898 # ------------------------- Hex Elements ------------------------- 899 900 # Restrictions 901 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 902 for i in range(nelem_hex): 903 col = i % nx_hex 904 row = i // nx_hex 905 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 906 907 for j in range(p_hex): 908 for k in range(p_hex): 909 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 910 k * (nx_hex * 2 + 1) + j 911 912 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 913 dim * ndofs, indx_hex, cmode=libceed.USE_POINTER) 914 915 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 916 indx_hex, cmode=libceed.USE_POINTER) 917 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 918 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 919 nqpts_hex, strides) 920 921 # Bases 922 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 923 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 924 925 # QFunctions 926 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 927 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 928 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 929 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 930 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 931 932 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 933 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 934 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 935 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 936 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 937 938 # Operators 939 op_setup_hex = ceed.Operator(qf_setup_tet) 940 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 941 libceed.VECTOR_NONE) 942 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 943 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 944 qdata_hex) 945 946 op_mass_hex = ceed.Operator(qf_mass_hex) 947 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 948 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 949 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 950 951 # ------------------------- Composite Operators ------------------------- 952 953 # Setup 954 op_setup = ceed.CompositeOperator() 955 op_setup.add_sub(op_setup_tet) 956 op_setup.add_sub(op_setup_hex) 957 op_setup.apply(x, libceed.VECTOR_NONE) 958 959 # Apply mass matrix 960 op_mass = ceed.CompositeOperator() 961 op_mass.add_sub(op_mass_tet) 962 op_mass.add_sub(op_mass_hex) 963 964 u.set_value(0.) 965 op_mass.apply(u, v) 966 967 # Check 968 with v.array_read() as v_array: 969 for i in range(ndofs): 970 assert abs(v_array[i]) < TOL 971 972# ------------------------------------------------------------------------------- 973# Test creation, action, and destruction for composite mass matrix operator 974# ------------------------------------------------------------------------------- 975 976 977def test_521(ceed_resource): 978 ceed = libceed.Ceed(ceed_resource) 979 980 nelem_tet, p_tet, q_tet = 6, 6, 4 981 nelem_hex, p_hex, q_hex = 6, 3, 4 982 nx, ny = 3, 3 983 dim = 2 984 nx_tet, ny_tet, nx_hex = 3, 1, 3 985 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 986 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 987 988 # Vectors 989 x = ceed.Vector(dim * ndofs) 990 x_array = np.zeros(dim * ndofs) 991 for i in range(ny * 2 + 1): 992 for j in range(nx * 2 + 1): 993 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 994 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 995 x.set_array(x_array, cmode=libceed.USE_POINTER) 996 997 qdata_hex = ceed.Vector(nqpts_hex) 998 qdata_tet = ceed.Vector(nqpts_tet) 999 u = ceed.Vector(ndofs) 1000 v = ceed.Vector(ndofs) 1001 1002 # ------------------------- Tet Elements ------------------------- 1003 1004 # Restrictions 1005 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1006 for i in range(nelem_tet // 2): 1007 col = i % nx 1008 row = i // nx 1009 offset = col * 2 + row * (nx * 2 + 1) * 2 1010 1011 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1012 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1013 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1014 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1015 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1016 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1017 1018 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1019 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1020 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1021 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1022 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1023 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1024 1025 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1026 indx_tet, cmode=libceed.USE_POINTER) 1027 1028 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1029 cmode=libceed.USE_POINTER) 1030 strides = np.array([1, q_tet, q_tet], dtype="int32") 1031 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1032 strides) 1033 1034 # Bases 1035 qref = np.empty(dim * q_tet, dtype="float64") 1036 qweight = np.empty(q_tet, dtype="float64") 1037 interp, grad = bm.buildmats(qref, qweight) 1038 1039 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 1040 qweight) 1041 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 1042 qweight) 1043 1044 # QFunctions 1045 file_dir = os.path.dirname(os.path.abspath(__file__)) 1046 qfs = load_qfs_so() 1047 1048 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 1049 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1050 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 1051 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1052 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 1053 1054 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 1055 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1056 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 1057 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 1058 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 1059 1060 # Operators 1061 op_setup_tet = ceed.Operator(qf_setup_tet) 1062 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1063 libceed.VECTOR_NONE) 1064 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1065 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1066 qdata_tet) 1067 1068 op_mass_tet = ceed.Operator(qf_mass_tet) 1069 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1070 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1071 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1072 1073 # ------------------------- Hex Elements ------------------------- 1074 1075 # Restrictions 1076 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1077 for i in range(nelem_hex): 1078 col = i % nx_hex 1079 row = i // nx_hex 1080 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1081 1082 for j in range(p_hex): 1083 for k in range(p_hex): 1084 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1085 k * (nx_hex * 2 + 1) + j 1086 1087 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1088 dim * ndofs, indx_hex, cmode=libceed.USE_POINTER) 1089 1090 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1091 indx_hex, cmode=libceed.USE_POINTER) 1092 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1093 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1094 nqpts_hex, strides) 1095 1096 # Bases 1097 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1098 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1099 1100 # QFunctions 1101 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1102 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1103 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1104 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1105 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1106 1107 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1108 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1109 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1110 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1111 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1112 1113 # Operators 1114 op_setup_hex = ceed.Operator(qf_setup_tet) 1115 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1116 libceed.VECTOR_NONE) 1117 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1118 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1119 qdata_hex) 1120 1121 op_mass_hex = ceed.Operator(qf_mass_hex) 1122 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1123 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1124 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1125 1126 # ------------------------- Composite Operators ------------------------- 1127 1128 # Setup 1129 op_setup = ceed.CompositeOperator() 1130 op_setup.add_sub(op_setup_tet) 1131 op_setup.add_sub(op_setup_hex) 1132 op_setup.apply(x, libceed.VECTOR_NONE) 1133 1134 # Apply mass matrix 1135 op_mass = ceed.CompositeOperator() 1136 op_mass.add_sub(op_mass_tet) 1137 op_mass.add_sub(op_mass_hex) 1138 u.set_value(1.) 1139 op_mass.apply(u, v) 1140 1141 # Check 1142 with v.array_read() as v_array: 1143 total = 0.0 1144 for i in range(ndofs): 1145 total = total + v_array[i] 1146 assert abs(total - 1.0) < 1E-10 1147 1148# ------------------------------------------------------------------------------- 1149# Test viewing of composite mass matrix operator 1150# ------------------------------------------------------------------------------- 1151 1152 1153def test_523(ceed_resource, capsys): 1154 ceed = libceed.Ceed(ceed_resource) 1155 1156 nelem_tet, p_tet, q_tet = 6, 6, 4 1157 nelem_hex, p_hex, q_hex = 6, 3, 4 1158 nx, ny = 3, 3 1159 dim = 2 1160 nx_tet, ny_tet, nx_hex = 3, 1, 3 1161 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1162 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 1163 1164 # Vectors 1165 qdata_hex = ceed.Vector(nqpts_hex) 1166 qdata_tet = ceed.Vector(nqpts_tet) 1167 1168 # ------------------------- Tet Elements ------------------------- 1169 1170 # Restrictions 1171 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1172 for i in range(nelem_tet // 2): 1173 col = i % nx 1174 row = i // nx 1175 offset = col * 2 + row * (nx * 2 + 1) * 2 1176 1177 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1178 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1179 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1180 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1181 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1182 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1183 1184 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1185 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1186 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1187 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1188 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1189 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1190 1191 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1192 indx_tet, cmode=libceed.USE_POINTER) 1193 1194 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1195 cmode=libceed.USE_POINTER) 1196 strides = np.array([1, q_tet, q_tet], dtype="int32") 1197 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1198 strides) 1199 1200 # Bases 1201 qref = np.empty(dim * q_tet, dtype="float64") 1202 qweight = np.empty(q_tet, dtype="float64") 1203 interp, grad = bm.buildmats(qref, qweight) 1204 1205 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 1206 qweight) 1207 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 1208 qweight) 1209 1210 # QFunctions 1211 file_dir = os.path.dirname(os.path.abspath(__file__)) 1212 qfs = load_qfs_so() 1213 1214 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 1215 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1216 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 1217 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1218 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 1219 1220 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 1221 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1222 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 1223 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 1224 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 1225 1226 # Operators 1227 op_setup_tet = ceed.Operator(qf_setup_tet) 1228 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1229 libceed.VECTOR_NONE) 1230 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1231 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1232 qdata_tet) 1233 1234 op_mass_tet = ceed.Operator(qf_mass_tet) 1235 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1236 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1237 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1238 1239 # ------------------------- Hex Elements ------------------------- 1240 1241 # Restrictions 1242 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1243 for i in range(nelem_hex): 1244 col = i % nx_hex 1245 row = i // nx_hex 1246 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1247 1248 for j in range(p_hex): 1249 for k in range(p_hex): 1250 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1251 k * (nx_hex * 2 + 1) + j 1252 1253 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1254 dim * ndofs, indx_hex, 1255 cmode=libceed.USE_POINTER) 1256 1257 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1258 indx_hex, cmode=libceed.USE_POINTER) 1259 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1260 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1261 nqpts_hex, strides) 1262 1263 # Bases 1264 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1265 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1266 1267 # QFunctions 1268 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1269 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1270 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1271 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1272 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1273 1274 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1275 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1276 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1277 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1278 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1279 1280 # Operators 1281 op_setup_hex = ceed.Operator(qf_setup_tet) 1282 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1283 libceed.VECTOR_NONE) 1284 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1285 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1286 qdata_hex) 1287 1288 op_mass_hex = ceed.Operator(qf_mass_hex) 1289 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1290 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1291 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1292 1293 # ------------------------- Composite Operators ------------------------- 1294 1295 # Setup 1296 op_setup = ceed.CompositeOperator() 1297 op_setup.add_sub(op_setup_tet) 1298 op_setup.add_sub(op_setup_hex) 1299 1300 # Apply mass matrix 1301 op_mass = ceed.CompositeOperator() 1302 op_mass.add_sub(op_mass_tet) 1303 op_mass.add_sub(op_mass_hex) 1304 1305 # View 1306 print(op_setup) 1307 print(op_mass) 1308 1309 stdout, stderr, ref_stdout = check.output(capsys) 1310 assert not stderr 1311 assert stdout == ref_stdout 1312 1313# ------------------------------------------------------------------------------- 1314# CeedOperatorApplyAdd for composite operator 1315# ------------------------------------------------------------------------------- 1316 1317 1318def test_524(ceed_resource): 1319 ceed = libceed.Ceed(ceed_resource) 1320 1321 nelem_tet, p_tet, q_tet = 6, 6, 4 1322 nelem_hex, p_hex, q_hex = 6, 3, 4 1323 nx, ny = 3, 3 1324 dim = 2 1325 nx_tet, ny_tet, nx_hex = 3, 1, 3 1326 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1327 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 1328 1329 # Vectors 1330 x = ceed.Vector(dim * ndofs) 1331 x_array = np.zeros(dim * ndofs) 1332 for i in range(ny * 2 + 1): 1333 for j in range(nx * 2 + 1): 1334 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 1335 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 1336 x.set_array(x_array, cmode=libceed.USE_POINTER) 1337 1338 qdata_hex = ceed.Vector(nqpts_hex) 1339 qdata_tet = ceed.Vector(nqpts_tet) 1340 u = ceed.Vector(ndofs) 1341 v = ceed.Vector(ndofs) 1342 1343 # ------------------------- Tet Elements ------------------------- 1344 1345 # Restrictions 1346 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1347 for i in range(nelem_tet // 2): 1348 col = i % nx 1349 row = i // nx 1350 offset = col * 2 + row * (nx * 2 + 1) * 2 1351 1352 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1353 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1354 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1355 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1356 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1357 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1358 1359 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1360 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1361 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1362 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1363 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1364 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1365 1366 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1367 indx_tet, cmode=libceed.USE_POINTER) 1368 1369 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1370 cmode=libceed.USE_POINTER) 1371 strides = np.array([1, q_tet, q_tet], dtype="int32") 1372 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1373 strides) 1374 1375 # Bases 1376 qref = np.empty(dim * q_tet, dtype="float64") 1377 qweight = np.empty(q_tet, dtype="float64") 1378 interp, grad = bm.buildmats(qref, qweight) 1379 1380 bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref, 1381 qweight) 1382 bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref, 1383 qweight) 1384 1385 # QFunctions 1386 file_dir = os.path.dirname(os.path.abspath(__file__)) 1387 qfs = load_qfs_so() 1388 1389 qf_setup_tet = ceed.QFunction(1, qfs.setup_mass_2d, 1390 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1391 qf_setup_tet.add_input("weights", 1, libceed.EVAL_WEIGHT) 1392 qf_setup_tet.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1393 qf_setup_tet.add_output("rho", 1, libceed.EVAL_NONE) 1394 1395 qf_mass_tet = ceed.QFunction(1, qfs.apply_mass, 1396 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1397 qf_mass_tet.add_input("rho", 1, libceed.EVAL_NONE) 1398 qf_mass_tet.add_input("u", 1, libceed.EVAL_INTERP) 1399 qf_mass_tet.add_output("v", 1, libceed.EVAL_INTERP) 1400 1401 # Operators 1402 op_setup_tet = ceed.Operator(qf_setup_tet) 1403 op_setup_tet.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1404 libceed.VECTOR_NONE) 1405 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1406 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1407 qdata_tet) 1408 1409 op_mass_tet = ceed.Operator(qf_mass_tet) 1410 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1411 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1412 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1413 1414 # ------------------------- Hex Elements ------------------------- 1415 1416 # Restrictions 1417 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1418 for i in range(nelem_hex): 1419 col = i % nx_hex 1420 row = i // nx_hex 1421 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1422 1423 for j in range(p_hex): 1424 for k in range(p_hex): 1425 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1426 k * (nx_hex * 2 + 1) + j 1427 1428 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1429 dim * ndofs, indx_hex, 1430 cmode=libceed.USE_POINTER) 1431 1432 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1433 indx_hex, cmode=libceed.USE_POINTER) 1434 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1435 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1436 nqpts_hex, strides) 1437 1438 # Bases 1439 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1440 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1441 1442 # QFunctions 1443 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1444 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1445 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1446 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1447 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1448 1449 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1450 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1451 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1452 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1453 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1454 1455 # Operators 1456 op_setup_hex = ceed.Operator(qf_setup_tet) 1457 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1458 libceed.VECTOR_NONE) 1459 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1460 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1461 qdata_hex) 1462 1463 op_mass_hex = ceed.Operator(qf_mass_hex) 1464 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1465 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1466 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1467 1468 # ------------------------- Composite Operators ------------------------- 1469 1470 # Setup 1471 op_setup = ceed.CompositeOperator() 1472 op_setup.add_sub(op_setup_tet) 1473 op_setup.add_sub(op_setup_hex) 1474 op_setup.apply(x, libceed.VECTOR_NONE) 1475 1476 # Apply mass matrix 1477 op_mass = ceed.CompositeOperator() 1478 op_mass.add_sub(op_mass_tet) 1479 op_mass.add_sub(op_mass_hex) 1480 u.set_value(1.) 1481 op_mass.apply(u, v) 1482 1483 # Check 1484 with v.array_read() as v_array: 1485 total = 0.0 1486 for i in range(ndofs): 1487 total = total + v_array[i] 1488 assert abs(total - 1.0) < 1E-10 1489 1490 # ApplyAdd mass matrix 1491 v.set_value(1.) 1492 op_mass.apply_add(u, v) 1493 1494 # Check 1495 with v.array_read() as v_array: 1496 total = -ndofs 1497 for i in range(ndofs): 1498 total = total + v_array[i] 1499 assert abs(total - 1.0) < 1E-10 1500 1501# ------------------------------------------------------------------------------- 1502# Test assembly of mass matrix operator diagonal 1503# ------------------------------------------------------------------------------- 1504 1505 1506def test_533(ceed_resource): 1507 ceed = libceed.Ceed(ceed_resource) 1508 1509 nelem = 6 1510 p = 3 1511 q = 4 1512 dim = 2 1513 nx = 3 1514 ny = 2 1515 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1516 nqpts = nelem * q * q 1517 1518 # Vectors 1519 x = ceed.Vector(dim * ndofs) 1520 x_array = np.zeros(dim * ndofs) 1521 for i in range(nx * 2 + 1): 1522 for j in range(ny * 2 + 1): 1523 x_array[i + j * (nx * 2 + 1) + 0 * ndofs] = i / (2 * nx) 1524 x_array[i + j * (nx * 2 + 1) + 1 * ndofs] = j / (2 * ny) 1525 x.set_array(x_array, cmode=libceed.USE_POINTER) 1526 1527 qdata = ceed.Vector(nqpts) 1528 u = ceed.Vector(ndofs) 1529 v = ceed.Vector(ndofs) 1530 1531 # Restrictions 1532 indx = np.zeros(nelem * p * p, dtype="int32") 1533 for i in range(nelem): 1534 col = i % nx 1535 row = i // nx 1536 offset = col * (p - 1) + row * (nx * 2 + 1) * (p - 1) 1537 for j in range(p): 1538 for k in range(p): 1539 indx[p * (p * i + k) + j] = offset + k * (nx * 2 + 1) + j 1540 rx = ceed.ElemRestriction(nelem, p * p, dim, ndofs, dim * ndofs, 1541 indx, cmode=libceed.USE_POINTER) 1542 1543 ru = ceed.ElemRestriction(nelem, p * p, 1, 1, ndofs, indx, 1544 cmode=libceed.USE_POINTER) 1545 strides = np.array([1, q * q, q * q], dtype="int32") 1546 rui = ceed.StridedElemRestriction(nelem, q * q, 1, nqpts, strides) 1547 1548 # Bases 1549 bx = ceed.BasisTensorH1Lagrange(dim, dim, p, q, libceed.GAUSS) 1550 bu = ceed.BasisTensorH1Lagrange(dim, 1, p, q, libceed.GAUSS) 1551 1552 # QFunctions 1553 qf_setup = ceed.QFunctionByName("Mass2DBuild") 1554 1555# ------------------------------------------------------------------------------- 1556# Test creation, action, and destruction for mass matrix operator with 1557# multigrid level, tensor basis and interpolation basis generation 1558# ------------------------------------------------------------------------------- 1559 1560 1561def test_550(ceed_resource): 1562 ceed = libceed.Ceed(ceed_resource) 1563 1564 nelem = 15 1565 p_coarse = 3 1566 p_fine = 5 1567 q = 8 1568 ncomp = 2 1569 nx = nelem + 1 1570 nu_coarse = nelem * (p_coarse - 1) + 1 1571 nu_fine = nelem * (p_fine - 1) + 1 1572 1573 # Vectors 1574 x = ceed.Vector(nx) 1575 x_array = np.zeros(nx) 1576 for i in range(nx): 1577 x_array[i] = i / (nx - 1.0) 1578 x.set_array(x_array, cmode=libceed.USE_POINTER) 1579 1580 qdata = ceed.Vector(nelem * q) 1581 u_coarse = ceed.Vector(ncomp * nu_coarse) 1582 v_coarse = ceed.Vector(ncomp * nu_coarse) 1583 u_fine = ceed.Vector(ncomp * nu_fine) 1584 v_fine = ceed.Vector(ncomp * nu_fine) 1585 1586 # Restrictions 1587 indx = np.zeros(nx * 2, dtype="int32") 1588 for i in range(nx): 1589 indx[2 * i + 0] = i 1590 indx[2 * i + 1] = i + 1 1591 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 1592 cmode=libceed.USE_POINTER) 1593 1594 indu_coarse = np.zeros(nelem * p_coarse, dtype="int32") 1595 for i in range(nelem): 1596 for j in range(p_coarse): 1597 indu_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j 1598 ru_coarse = ceed.ElemRestriction(nelem, p_coarse, ncomp, nu_coarse, 1599 ncomp * nu_coarse, indu_coarse, 1600 cmode=libceed.USE_POINTER) 1601 1602 indu_fine = np.zeros(nelem * p_fine, dtype="int32") 1603 for i in range(nelem): 1604 for j in range(p_fine): 1605 indu_fine[p_fine * i + j] = i * (p_fine - 1) + j 1606 ru_fine = ceed.ElemRestriction(nelem, p_fine, ncomp, nu_fine, 1607 ncomp * nu_fine, indu_fine, 1608 cmode=libceed.USE_POINTER) 1609 strides = np.array([1, q, q], dtype="int32") 1610 1611 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 1612 1613 # Bases 1614 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 1615 bu_coarse = ceed.BasisTensorH1Lagrange(1, ncomp, p_coarse, q, libceed.GAUSS) 1616 bu_fine = ceed.BasisTensorH1Lagrange(1, ncomp, p_fine, q, libceed.GAUSS) 1617 1618 # QFunctions 1619 file_dir = os.path.dirname(os.path.abspath(__file__)) 1620 qfs = load_qfs_so() 1621 1622 qf_setup = ceed.QFunction(1, qfs.setup_mass, 1623 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 1624 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 1625 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 1626 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 1627 1628 qf_mass = ceed.QFunction(1, qfs.apply_mass_two, 1629 os.path.join(file_dir, "test-qfunctions.h:apply_mass_two")) 1630 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 1631 qf_mass.add_input("u", ncomp, libceed.EVAL_INTERP) 1632 qf_mass.add_output("v", ncomp, libceed.EVAL_INTERP) 1633 1634 # Operators 1635 op_setup = ceed.Operator(qf_setup) 1636 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 1637 libceed.VECTOR_NONE) 1638 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 1639 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 1640 libceed.VECTOR_ACTIVE) 1641 1642 op_mass_fine = ceed.Operator(qf_mass) 1643 op_mass_fine.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 1644 op_mass_fine.set_field("u", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1645 op_mass_fine.set_field("v", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1646 1647 # Setup 1648 op_setup.apply(x, qdata) 1649 1650 # Create multigrid level 1651 p_mult_fine = ceed.Vector(ncomp * nu_fine) 1652 p_mult_fine.set_value(1.0) 1653 [op_mass_coarse, op_prolong, op_restrict] = op_mass_fine.multigrid_create(p_mult_fine, 1654 ru_coarse, 1655 bu_coarse) 1656 1657 # Apply coarse mass matrix 1658 u_coarse.set_value(1.0) 1659 op_mass_coarse.apply(u_coarse, v_coarse) 1660 1661 # Check 1662 with v_coarse.array_read() as v_array: 1663 total = 0.0 1664 for i in range(nu_coarse * ncomp): 1665 total = total + v_array[i] 1666 assert abs(total - 2.0) < 1E-13 1667 1668 # Prolong coarse u 1669 op_prolong.apply(u_coarse, u_fine) 1670 1671 # Apply mass matrix 1672 op_mass_fine.apply(u_fine, v_fine) 1673 1674 # Check 1675 with v_fine.array_read() as v_array: 1676 total = 0.0 1677 for i in range(nu_fine * ncomp): 1678 total = total + v_array[i] 1679 assert abs(total - 2.0) < 1E-13 1680 1681 # Restrict state to coarse grid 1682 op_restrict.apply(v_fine, v_coarse) 1683 1684 # Check 1685 with v_coarse.array_read() as v_array: 1686 total = 0.0 1687 for i in range(nu_coarse * ncomp): 1688 total = total + v_array[i] 1689 assert abs(total - 2.0) < 1E-13 1690 1691# ------------------------------------------------------------------------------- 1692# Test creation, action, and destruction for mass matrix operator with 1693# multigrid level, tensor basis 1694# ------------------------------------------------------------------------------- 1695 1696 1697def test_552(ceed_resource): 1698 ceed = libceed.Ceed(ceed_resource) 1699 1700 nelem = 15 1701 p_coarse = 3 1702 p_fine = 5 1703 q = 8 1704 ncomp = 2 1705 nx = nelem + 1 1706 nu_coarse = nelem * (p_coarse - 1) + 1 1707 nu_fine = nelem * (p_fine - 1) + 1 1708 1709 # Vectors 1710 x = ceed.Vector(nx) 1711 x_array = np.zeros(nx) 1712 for i in range(nx): 1713 x_array[i] = i / (nx - 1.0) 1714 x.set_array(x_array, cmode=libceed.USE_POINTER) 1715 1716 qdata = ceed.Vector(nelem * q) 1717 u_coarse = ceed.Vector(ncomp * nu_coarse) 1718 v_coarse = ceed.Vector(ncomp * nu_coarse) 1719 u_fine = ceed.Vector(ncomp * nu_fine) 1720 v_fine = ceed.Vector(ncomp * nu_fine) 1721 1722 # Restrictions 1723 indx = np.zeros(nx * 2, dtype="int32") 1724 for i in range(nx): 1725 indx[2 * i + 0] = i 1726 indx[2 * i + 1] = i + 1 1727 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 1728 cmode=libceed.USE_POINTER) 1729 1730 indu_coarse = np.zeros(nelem * p_coarse, dtype="int32") 1731 for i in range(nelem): 1732 for j in range(p_coarse): 1733 indu_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j 1734 ru_coarse = ceed.ElemRestriction(nelem, p_coarse, ncomp, nu_coarse, 1735 ncomp * nu_coarse, indu_coarse, 1736 cmode=libceed.USE_POINTER) 1737 1738 indu_fine = np.zeros(nelem * p_fine, dtype="int32") 1739 for i in range(nelem): 1740 for j in range(p_fine): 1741 indu_fine[p_fine * i + j] = i * (p_fine - 1) + j 1742 ru_fine = ceed.ElemRestriction(nelem, p_fine, ncomp, nu_fine, 1743 ncomp * nu_fine, indu_fine, 1744 cmode=libceed.USE_POINTER) 1745 strides = np.array([1, q, q], dtype="int32") 1746 1747 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 1748 1749 # Bases 1750 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 1751 bu_coarse = ceed.BasisTensorH1Lagrange(1, ncomp, p_coarse, q, libceed.GAUSS) 1752 bu_fine = ceed.BasisTensorH1Lagrange(1, ncomp, p_fine, q, libceed.GAUSS) 1753 1754 # QFunctions 1755 file_dir = os.path.dirname(os.path.abspath(__file__)) 1756 qfs = load_qfs_so() 1757 1758 qf_setup = ceed.QFunction(1, qfs.setup_mass, 1759 os.path.join(file_dir, "test-qfunctions.h:setup_mass")) 1760 qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT) 1761 qf_setup.add_input("dx", 1, libceed.EVAL_GRAD) 1762 qf_setup.add_output("rho", 1, libceed.EVAL_NONE) 1763 1764 qf_mass = ceed.QFunction(1, qfs.apply_mass_two, 1765 os.path.join(file_dir, "test-qfunctions.h:apply_mass_two")) 1766 qf_mass.add_input("rho", 1, libceed.EVAL_NONE) 1767 qf_mass.add_input("u", ncomp, libceed.EVAL_INTERP) 1768 qf_mass.add_output("v", ncomp, libceed.EVAL_INTERP) 1769 1770 # Operators 1771 op_setup = ceed.Operator(qf_setup) 1772 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 1773 libceed.VECTOR_NONE) 1774 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 1775 op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED, 1776 libceed.VECTOR_ACTIVE) 1777 1778 op_mass_fine = ceed.Operator(qf_mass) 1779 op_mass_fine.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata) 1780 op_mass_fine.set_field("u", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1781 op_mass_fine.set_field("v", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1782 1783 # Setup 1784 op_setup.apply(x, qdata) 1785 1786 # Create multigrid level 1787 p_mult_fine = ceed.Vector(ncomp * nu_fine) 1788 p_mult_fine.set_value(1.0) 1789 b_c_to_f = ceed.BasisTensorH1Lagrange( 1790 1, ncomp, p_coarse, p_fine, libceed.GAUSS_LOBATTO) 1791 interp_C_to_F = b_c_to_f.get_interp1d() 1792 [op_mass_coarse, op_prolong, op_restrict] = op_mass_fine.multigrid_create_tensor_h1(p_mult_fine, 1793 ru_coarse, bu_coarse, interp_C_to_F) 1794 1795 # Apply coarse mass matrix 1796 u_coarse.set_value(1.0) 1797 op_mass_coarse.apply(u_coarse, v_coarse) 1798 1799 # Check 1800 with v_coarse.array_read() as v_array: 1801 total = 0.0 1802 for i in range(nu_coarse * ncomp): 1803 total = total + v_array[i] 1804 assert abs(total - 2.0) < 1E-13 1805 1806 # Prolong coarse u 1807 op_prolong.apply(u_coarse, u_fine) 1808 1809 # Apply mass matrix 1810 op_mass_fine.apply(u_fine, v_fine) 1811 1812 # Check 1813 with v_fine.array_read() as v_array: 1814 total = 0.0 1815 for i in range(nu_fine * ncomp): 1816 total = total + v_array[i] 1817 assert abs(total - 2.0) < 1E-13 1818 1819 # Restrict state to coarse grid 1820 op_restrict.apply(v_fine, v_coarse) 1821 1822 # Check 1823 with v_coarse.array_read() as v_array: 1824 total = 0.0 1825 for i in range(nu_coarse * ncomp): 1826 total = total + v_array[i] 1827 assert abs(total - 2.0) < 1E-13 1828 1829# ------------------------------------------------------------------------------- 1830# Test creation, action, and destruction for mass matrix operator with 1831# multigrid level, non-tensor basis 1832# ------------------------------------------------------------------------------- 1833 1834 1835def test_553(ceed_resource): 1836 ceed = libceed.Ceed(ceed_resource) 1837 1838 nelem = 15 1839 p_coarse = 3 1840 p_fine = 5 1841 q = 8 1842 ncomp = 1 1843 nx = nelem + 1 1844 nu_coarse = nelem * (p_coarse - 1) + 1 1845 nu_fine = nelem * (p_fine - 1) + 1 1846 1847 # Vectors 1848 x = ceed.Vector(nx) 1849 x_array = np.zeros(nx) 1850 for i in range(nx): 1851 x_array[i] = i / (nx - 1.0) 1852 x.set_array(x_array, cmode=libceed.USE_POINTER) 1853 1854 qdata = ceed.Vector(nelem * q) 1855 u_coarse = ceed.Vector(ncomp * nu_coarse) 1856 v_coarse = ceed.Vector(ncomp * nu_coarse) 1857 u_fine = ceed.Vector(ncomp * nu_fine) 1858 v_fine = ceed.Vector(ncomp * nu_fine) 1859 1860 # Restrictions 1861 indx = np.zeros(nx * 2, dtype="int32") 1862 for i in range(nx): 1863 indx[2 * i + 0] = i 1864 indx[2 * i + 1] = i + 1 1865 rx = ceed.ElemRestriction(nelem, 2, 1, 1, nx, indx, 1866 cmode=libceed.USE_POINTER) 1867 1868 indu_coarse = np.zeros(nelem * p_coarse, dtype="int32") 1869 for i in range(nelem): 1870 for j in range(p_coarse): 1871 indu_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j 1872 ru_coarse = ceed.ElemRestriction(nelem, p_coarse, ncomp, nu_coarse, 1873 ncomp * nu_coarse, indu_coarse, 1874 cmode=libceed.USE_POINTER) 1875 1876 indu_fine = np.zeros(nelem * p_fine, dtype="int32") 1877 for i in range(nelem): 1878 for j in range(p_fine): 1879 indu_fine[p_fine * i + j] = i * (p_fine - 1) + j 1880 ru_fine = ceed.ElemRestriction(nelem, p_fine, ncomp, nu_fine, 1881 ncomp * nu_fine, indu_fine, 1882 cmode=libceed.USE_POINTER) 1883 strides = np.array([1, q, q], dtype="int32") 1884 1885 rui = ceed.StridedElemRestriction(nelem, q, 1, q * nelem, strides) 1886 1887 # Bases 1888 bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS) 1889 bu_coarse = ceed.BasisTensorH1Lagrange(1, ncomp, p_coarse, q, libceed.GAUSS) 1890 bu_fine = ceed.BasisTensorH1Lagrange(1, ncomp, p_fine, q, libceed.GAUSS) 1891 1892 # QFunctions 1893 file_dir = os.path.dirname(os.path.abspath(__file__)) 1894 qfs = load_qfs_so() 1895 1896 qf_setup = ceed.QFunctionByName("Mass1DBuild") 1897 qf_mass = ceed.QFunctionByName("MassApply") 1898 1899 # Operators 1900 op_setup = ceed.Operator(qf_setup) 1901 op_setup.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx, 1902 libceed.VECTOR_NONE) 1903 op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE) 1904 op_setup.set_field("qdata", rui, libceed.BASIS_COLLOCATED, 1905 libceed.VECTOR_ACTIVE) 1906 1907 op_mass_fine = ceed.Operator(qf_mass) 1908 op_mass_fine.set_field("qdata", rui, libceed.BASIS_COLLOCATED, qdata) 1909 op_mass_fine.set_field("u", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1910 op_mass_fine.set_field("v", ru_fine, bu_fine, libceed.VECTOR_ACTIVE) 1911 1912 # Setup 1913 op_setup.apply(x, qdata) 1914 1915 # Create multigrid level 1916 p_mult_fine = ceed.Vector(ncomp * nu_fine) 1917 p_mult_fine.set_value(1.0) 1918 b_c_to_f = ceed.BasisTensorH1Lagrange( 1919 1, ncomp, p_coarse, p_fine, libceed.GAUSS_LOBATTO) 1920 interp_C_to_F = b_c_to_f.get_interp1d() 1921 [op_mass_coarse, op_prolong, op_restrict] = op_mass_fine.multigrid_create_h1(p_mult_fine, 1922 ru_coarse, bu_coarse, interp_C_to_F) 1923 1924 # Apply coarse mass matrix 1925 u_coarse.set_value(1.0) 1926 op_mass_coarse.apply(u_coarse, v_coarse) 1927 1928 # Check 1929 with v_coarse.array_read() as v_array: 1930 total = 0.0 1931 for i in range(nu_coarse * ncomp): 1932 total = total + v_array[i] 1933 assert abs(total - 1.0) < 1E-13 1934 1935 # Prolong coarse u 1936 op_prolong.apply(u_coarse, u_fine) 1937 1938 # Apply mass matrix 1939 op_mass_fine.apply(u_fine, v_fine) 1940 1941 # Check 1942 with v_fine.array_read() as v_array: 1943 total = 0.0 1944 for i in range(nu_fine * ncomp): 1945 total = total + v_array[i] 1946 assert abs(total - 1.0) < 1E-13 1947 1948 # Restrict state to coarse grid 1949 op_restrict.apply(v_fine, v_coarse) 1950 1951 # Check 1952 with v_coarse.array_read() as v_array: 1953 total = 0.0 1954 for i in range(nu_coarse * ncomp): 1955 total = total + v_array[i] 1956 assert abs(total - 1.0) < 1E-13 1957 1958# ------------------------------------------------------------------------------- 1959