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