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 distutils.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.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_tet, 1228 libceed.VECTOR_NONE) 1229 op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE) 1230 op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, 1231 qdata_tet) 1232 1233 op_mass_tet = ceed.Operator(qf_mass_tet) 1234 op_mass_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED, qdata_tet) 1235 op_mass_tet.set_field("u", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1236 op_mass_tet.set_field("v", ru_tet, bu_tet, libceed.VECTOR_ACTIVE) 1237 1238 # ------------------------- Hex Elements ------------------------- 1239 1240 # Restrictions 1241 indx_hex = np.zeros(nelem_hex * p_hex * p_hex, dtype="int32") 1242 for i in range(nelem_hex): 1243 col = i % nx_hex 1244 row = i // nx_hex 1245 offset = (nx_tet * 2 + 1) * (ny_tet * 2) * (1 + row) + col * 2 1246 1247 for j in range(p_hex): 1248 for k in range(p_hex): 1249 indx_hex[p_hex * (p_hex * i + k) + j] = offset + \ 1250 k * (nx_hex * 2 + 1) + j 1251 1252 rx_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, dim, ndofs, 1253 dim * ndofs, indx_hex, 1254 cmode=libceed.USE_POINTER) 1255 1256 ru_hex = ceed.ElemRestriction(nelem_hex, p_hex * p_hex, 1, 1, ndofs, 1257 indx_hex, cmode=libceed.USE_POINTER) 1258 strides = np.array([1, q_hex * q_hex, q_hex * q_hex], dtype="int32") 1259 rui_hex = ceed.StridedElemRestriction(nelem_hex, q_hex * q_hex, 1, 1260 nqpts_hex, strides) 1261 1262 # Bases 1263 bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS) 1264 bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS) 1265 1266 # QFunctions 1267 qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d, 1268 os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d")) 1269 qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT) 1270 qf_setup_hex.add_input("dx", dim * dim, libceed.EVAL_GRAD) 1271 qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE) 1272 1273 qf_mass_hex = ceed.QFunction(1, qfs.apply_mass, 1274 os.path.join(file_dir, "test-qfunctions.h:apply_mass")) 1275 qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE) 1276 qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP) 1277 qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP) 1278 1279 # Operators 1280 op_setup_hex = ceed.Operator(qf_setup_tet) 1281 op_setup_hex.set_field("weights", libceed.ELEMRESTRICTION_NONE, bx_hex, 1282 libceed.VECTOR_NONE) 1283 op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE) 1284 op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, 1285 qdata_hex) 1286 1287 op_mass_hex = ceed.Operator(qf_mass_hex) 1288 op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex) 1289 op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1290 op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE) 1291 1292 # ------------------------- Composite Operators ------------------------- 1293 1294 # Setup 1295 op_setup = ceed.CompositeOperator() 1296 op_setup.add_sub(op_setup_tet) 1297 op_setup.add_sub(op_setup_hex) 1298 1299 # Apply mass matrix 1300 op_mass = ceed.CompositeOperator() 1301 op_mass.add_sub(op_mass_tet) 1302 op_mass.add_sub(op_mass_hex) 1303 1304 # View 1305 print(op_setup) 1306 print(op_mass) 1307 1308 stdout, stderr, ref_stdout = check.output(capsys) 1309 assert not stderr 1310 assert stdout == ref_stdout 1311 1312# ------------------------------------------------------------------------------- 1313# CeedOperatorApplyAdd for composite operator 1314# ------------------------------------------------------------------------------- 1315 1316 1317def test_524(ceed_resource): 1318 ceed = libceed.Ceed(ceed_resource) 1319 1320 nelem_tet, p_tet, q_tet = 6, 6, 4 1321 nelem_hex, p_hex, q_hex = 6, 3, 4 1322 nx, ny = 3, 3 1323 dim = 2 1324 nx_tet, ny_tet, nx_hex = 3, 1, 3 1325 ndofs = (nx * 2 + 1) * (ny * 2 + 1) 1326 nqpts_tet, nqpts_hex = nelem_tet * q_tet, nelem_hex * q_hex * q_hex 1327 1328 # Vectors 1329 x = ceed.Vector(dim * ndofs) 1330 x_array = np.zeros(dim * ndofs, dtype=ceed.scalar_type()) 1331 for i in range(ny * 2 + 1): 1332 for j in range(nx * 2 + 1): 1333 x_array[i + j * (ny * 2 + 1)] = i / (2 * ny) 1334 x_array[i + j * (ny * 2 + 1) + ndofs] = j / (2 * nx) 1335 x.set_array(x_array, cmode=libceed.USE_POINTER) 1336 1337 qdata_hex = ceed.Vector(nqpts_hex) 1338 qdata_tet = ceed.Vector(nqpts_tet) 1339 u = ceed.Vector(ndofs) 1340 v = ceed.Vector(ndofs) 1341 1342 # ------------------------- Tet Elements ------------------------- 1343 1344 # Restrictions 1345 indx_tet = np.zeros(nelem_tet * p_tet, dtype="int32") 1346 for i in range(nelem_tet // 2): 1347 col = i % nx 1348 row = i // nx 1349 offset = col * 2 + row * (nx * 2 + 1) * 2 1350 1351 indx_tet[i * 2 * p_tet + 0] = 2 + offset 1352 indx_tet[i * 2 * p_tet + 1] = 9 + offset 1353 indx_tet[i * 2 * p_tet + 2] = 16 + offset 1354 indx_tet[i * 2 * p_tet + 3] = 1 + offset 1355 indx_tet[i * 2 * p_tet + 4] = 8 + offset 1356 indx_tet[i * 2 * p_tet + 5] = 0 + offset 1357 1358 indx_tet[i * 2 * p_tet + 6] = 14 + offset 1359 indx_tet[i * 2 * p_tet + 7] = 7 + offset 1360 indx_tet[i * 2 * p_tet + 8] = 0 + offset 1361 indx_tet[i * 2 * p_tet + 9] = 15 + offset 1362 indx_tet[i * 2 * p_tet + 10] = 8 + offset 1363 indx_tet[i * 2 * p_tet + 11] = 16 + offset 1364 1365 rx_tet = ceed.ElemRestriction(nelem_tet, p_tet, dim, ndofs, dim * ndofs, 1366 indx_tet, cmode=libceed.USE_POINTER) 1367 1368 ru_tet = ceed.ElemRestriction(nelem_tet, p_tet, 1, 1, ndofs, indx_tet, 1369 cmode=libceed.USE_POINTER) 1370 strides = np.array([1, q_tet, q_tet], dtype="int32") 1371 rui_tet = ceed.StridedElemRestriction(nelem_tet, q_tet, 1, nqpts_tet, 1372 strides) 1373 1374 # Bases 1375 qref = np.empty(dim * q_tet, dtype=ceed.scalar_type()) 1376 qweight = np.empty(q_tet, dtype=ceed.scalar_type()) 1377 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 1378 libceed.lib.CEED_SCALAR_TYPE]) 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) < 10. * TOL 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) < 10. * TOL 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, dtype=ceed.scalar_type()) 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) < 10. * TOL 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) < 10. * TOL 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) < 10. * TOL 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, dtype=ceed.scalar_type()) 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_interp_1d() 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) < TOL 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) < TOL 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) < TOL 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, dtype=ceed.scalar_type()) 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_interp_1d() 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) < TOL 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) < TOL 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) < TOL 1957 1958# ------------------------------------------------------------------------------- 1959