1# Copyright (c) 2017-2026, 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 Basis functionality 10 11import os 12import math 13import libceed 14import numpy as np 15import buildmats as bm 16import check 17 18TOL = libceed.EPSILON * 256 19 20# ------------------------------------------------------------------------------- 21# Utilities 22# ------------------------------------------------------------------------------- 23 24 25def eval(dim, x): 26 result, center = 1, 0.1 27 for d in range(dim): 28 result *= math.tanh(x[d] - center) 29 center += 0.1 30 return result 31 32 33def feval(x1, x2): 34 return x1 * x1 + x2 * x2 + x1 * x2 + 1 35 36 37def dfeval(x1, x2): 38 return 2 * x1 + x2 39 40# ------------------------------------------------------------------------------- 41# Test creation and distruction of a H1Lagrange basis 42# ------------------------------------------------------------------------------- 43 44 45def test_300(ceed_resource, capsys): 46 ceed = libceed.Ceed(ceed_resource) 47 48 b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS_LOBATTO) 49 print(b) 50 del b 51 52 b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS) 53 print(b) 54 del b 55 56 # Only run this test in double precision 57 if libceed.lib.CEED_SCALAR_TYPE == libceed.SCALAR_FP64: 58 stdout, stderr, ref_stdout = check.output(capsys) 59 assert not stderr 60 assert stdout == ref_stdout 61 62# ------------------------------------------------------------------------------- 63# Test GetNumNodes and GetNumQuadraturePoints for basis 64# ------------------------------------------------------------------------------- 65 66 67def test_306(ceed_resource): 68 ceed = libceed.Ceed(ceed_resource) 69 70 b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO) 71 72 p = b.get_num_nodes() 73 q = b.get_num_quadrature_points() 74 75 assert p == 64 76 assert q == 125 77 78# ------------------------------------------------------------------------------- 79# Test interpolation in multiple dimensions 80# ------------------------------------------------------------------------------- 81 82 83def test_313(ceed_resource): 84 ceed = libceed.Ceed(ceed_resource) 85 86 for dim in range(1, 4): 87 Q = 10 88 Qdim = Q**dim 89 Xdim = 2**dim 90 x = np.empty(Xdim * dim, dtype=ceed.scalar_type()) 91 uq = np.empty(Qdim, dtype=ceed.scalar_type()) 92 93 for d in range(dim): 94 for i in range(Xdim): 95 x[d * Xdim + i] = 1 if (i % 96 (2**(dim - d))) // (2**(dim - d - 1)) else -1 97 98 X = ceed.Vector(Xdim * dim) 99 X.set_array(x, cmode=libceed.USE_POINTER) 100 Xq = ceed.Vector(Qdim * dim) 101 Xq.set_value(0) 102 U = ceed.Vector(Qdim) 103 U.set_value(0) 104 Uq = ceed.Vector(Qdim) 105 106 bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO) 107 bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO) 108 109 bxl.apply(1, libceed.EVAL_INTERP, X, Xq) 110 111 with Xq.array_read() as xq: 112 for i in range(Qdim): 113 xx = np.empty(dim, dtype=ceed.scalar_type()) 114 for d in range(dim): 115 xx[d] = xq[d * Qdim + i] 116 uq[i] = eval(dim, xx) 117 118 Uq.set_array(uq, cmode=libceed.USE_POINTER) 119 120 # This operation is the identity because the quadrature is collocated 121 bul.T.apply(1, libceed.EVAL_INTERP, Uq, U) 122 123 bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS) 124 bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS) 125 126 bxg.apply(1, libceed.EVAL_INTERP, X, Xq) 127 bug.apply(1, libceed.EVAL_INTERP, U, Uq) 128 129 with Xq.array_read() as xq, Uq.array_read() as u: 130 for i in range(Qdim): 131 xx = np.empty(dim, dtype=ceed.scalar_type()) 132 for d in range(dim): 133 xx[d] = xq[d * Qdim + i] 134 fx = eval(dim, xx) 135 assert math.fabs(u[i] - fx) < 1E-4 136 137# ------------------------------------------------------------------------------- 138# Test grad in multiple dimensions 139# ------------------------------------------------------------------------------- 140 141 142def test_314(ceed_resource): 143 ceed = libceed.Ceed(ceed_resource) 144 145 for dim in range(1, 4): 146 P, Q = 8, 10 147 Pdim = P**dim 148 Qdim = Q**dim 149 Xdim = 2**dim 150 sum1 = sum2 = 0 151 x = np.empty(Xdim * dim, dtype=ceed.scalar_type()) 152 u = np.empty(Pdim, dtype=ceed.scalar_type()) 153 154 for d in range(dim): 155 for i in range(Xdim): 156 x[d * Xdim + i] = 1 if (i % 157 (2**(dim - d))) // (2**(dim - d - 1)) else -1 158 159 X = ceed.Vector(Xdim * dim) 160 X.set_array(x, cmode=libceed.USE_POINTER) 161 Xq = ceed.Vector(Pdim * dim) 162 Xq.set_value(0) 163 U = ceed.Vector(Pdim) 164 Uq = ceed.Vector(Qdim * dim) 165 Uq.set_value(0) 166 Ones = ceed.Vector(Qdim * dim) 167 Ones.set_value(1) 168 Gtposeones = ceed.Vector(Pdim) 169 Gtposeones.set_value(0) 170 171 # Get function values at quadrature points 172 bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO) 173 bxl.apply(1, libceed.EVAL_INTERP, X, Xq) 174 175 with Xq.array_read() as xq: 176 for i in range(Pdim): 177 xx = np.empty(dim, dtype=ceed.scalar_type()) 178 for d in range(dim): 179 xx[d] = xq[d * Pdim + i] 180 u[i] = eval(dim, xx) 181 182 U.set_array(u, cmode=libceed.USE_POINTER) 183 184 # Calculate G u at quadrature points, G' * 1 at dofs 185 bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS) 186 bug.apply(1, libceed.EVAL_GRAD, U, Uq) 187 bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones) 188 189 # Check if 1' * G * u = u' * (G' * 1) 190 with Gtposeones.array_read() as gtposeones, Uq.array_read() as uq: 191 for i in range(Pdim): 192 sum1 += gtposeones[i] * u[i] 193 for i in range(dim * Qdim): 194 sum2 += uq[i] 195 196 assert math.fabs(sum1 - sum2) < 10. * TOL 197 198# ------------------------------------------------------------------------------- 199# Test creation and destruction of a 2D Simplex non-tensor H1 basis 200# ------------------------------------------------------------------------------- 201 202 203def test_320(ceed_resource): 204 ceed = libceed.Ceed(ceed_resource) 205 206 P, Q, dim = 6, 4, 2 207 208 in_array = np.empty(P, dtype=ceed.scalar_type()) 209 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 210 qweight = np.empty(Q, dtype=ceed.scalar_type()) 211 212 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 213 libceed.lib.CEED_SCALAR_TYPE]) 214 215 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 216 217 print(b) 218 del b 219 220# ------------------------------------------------------------------------------- 221# Test integration with a 2D Simplex non-tensor H1 basis 222# ------------------------------------------------------------------------------- 223 224 225def test_322(ceed_resource): 226 ceed = libceed.Ceed(ceed_resource) 227 228 P, Q, dim = 6, 4, 2 229 230 xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 231 0., 0.5, 0.5, 1.], dtype=ceed.scalar_type()) 232 in_array = np.empty(P, dtype=ceed.scalar_type()) 233 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 234 qweight = np.empty(Q, dtype=ceed.scalar_type()) 235 236 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 237 libceed.lib.CEED_SCALAR_TYPE]) 238 239 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 240 241 # Interpolate function to quadrature points 242 for i in range(P): 243 in_array[i] = feval(xr[0 * P + i], xr[1 * P + i]) 244 245 in_vec = ceed.Vector(P) 246 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 247 out_vec = ceed.Vector(Q) 248 out_vec.set_value(0) 249 weights_vec = ceed.Vector(Q) 250 weights_vec.set_value(0) 251 252 b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) 253 b.apply(1, libceed.EVAL_WEIGHT, libceed.VECTOR_NONE, weights_vec) 254 255 # Check values at quadrature points 256 with out_vec.array_read() as out_array, weights_vec.array_read() as weights_array: 257 sum = 0 258 for i in range(Q): 259 sum += out_array[i] * weights_array[i] 260 assert math.fabs(sum - 17. / 24.) < 10. * TOL 261 262# ------------------------------------------------------------------------------- 263# Test grad with a 2D Simplex non-tensor H1 basis 264# ------------------------------------------------------------------------------- 265 266 267def test_323(ceed_resource): 268 ceed = libceed.Ceed(ceed_resource) 269 270 P, Q, dim = 6, 4, 2 271 272 xq = np.array([0.2, 0.6, 1. / 3., 0.2, 0.2, 0.2, 273 1. / 3., 0.6], dtype=ceed.scalar_type()) 274 xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 275 0., 0.5, 0.5, 1.], dtype=ceed.scalar_type()) 276 in_array = np.empty(P, dtype=ceed.scalar_type()) 277 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 278 qweight = np.empty(Q, dtype=ceed.scalar_type()) 279 280 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 281 libceed.lib.CEED_SCALAR_TYPE]) 282 283 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 284 285 # Interpolate function to quadrature points 286 for i in range(P): 287 in_array[i] = feval(xr[0 * P + i], xr[1 * P + i]) 288 289 in_vec = ceed.Vector(P) 290 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 291 out_vec = ceed.Vector(Q * dim) 292 out_vec.set_value(0) 293 294 b.apply(1, libceed.EVAL_GRAD, in_vec, out_vec) 295 296 # Check values at quadrature points 297 with out_vec.array_read() as out_array: 298 for i in range(Q): 299 value = dfeval(xq[0 * Q + i], xq[1 * Q + i]) 300 assert math.fabs(out_array[0 * Q + i] - value) < 10. * TOL 301 302 value = dfeval(xq[1 * Q + i], xq[0 * Q + i]) 303 assert math.fabs(out_array[1 * Q + i] - value) < 10. * TOL 304 305# ------------------------------------------------------------------------------- 306# Test interpolation with a 2D Quad non-tensor H(div) basis 307# ------------------------------------------------------------------------------- 308 309 310def test_330(ceed_resource): 311 ceed = libceed.Ceed(ceed_resource) 312 313 P, Q, dim = 4, 4, 2 314 315 in_array = np.ones(P, dtype=ceed.scalar_type()) 316 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 317 qweight = np.empty(Q, dtype=ceed.scalar_type()) 318 319 interp, div = bm.buildmatshdiv(qref, qweight, libceed.scalar_types[ 320 libceed.lib.CEED_SCALAR_TYPE]) 321 322 b = ceed.BasisHdiv(libceed.QUAD, 1, P, Q, interp, div, qref, qweight) 323 324 # Interpolate function to quadrature points 325 in_vec = ceed.Vector(P) 326 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 327 out_vec = ceed.Vector(Q * dim) 328 out_vec.set_value(0) 329 330 b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) 331 332 # Check values at quadrature points 333 with out_vec.array_read() as out_array: 334 for i in range(Q): 335 assert math.fabs(out_array[0 * Q + i] + 1.) < 10. * TOL 336 assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL 337 338# ------------------------------------------------------------------------------- 339# Test interpolation with a 2D Simplex non-tensor H(curl) basis 340# ------------------------------------------------------------------------------- 341 342 343def test_340(ceed_resource): 344 ceed = libceed.Ceed(ceed_resource) 345 346 P, Q, dim = 3, 4, 2 347 348 in_array = np.empty(P, dtype=ceed.scalar_type()) 349 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 350 qweight = np.empty(Q, dtype=ceed.scalar_type()) 351 352 interp, curl = bm.buildmatshcurl(qref, qweight, libceed.scalar_types[ 353 libceed.lib.CEED_SCALAR_TYPE]) 354 355 b = ceed.BasisHcurl(libceed.TRIANGLE, 1, P, Q, interp, curl, qref, qweight) 356 357 # Interpolate function to quadrature points 358 in_array[0] = 1. 359 in_array[1] = 2. 360 in_array[2] = 1. 361 362 in_vec = ceed.Vector(P) 363 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 364 out_vec = ceed.Vector(Q * dim) 365 out_vec.set_value(0) 366 367 b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) 368 369 # Check values at quadrature points 370 with out_vec.array_read() as out_array: 371 for i in range(Q): 372 assert math.fabs(out_array[0 * Q + i] - 1.) < 10. * TOL 373 374 # Interpolate function to quadrature points 375 in_array[0] = -1. 376 in_array[1] = 1. 377 in_array[2] = 2. 378 379 out_vec.set_value(0) 380 381 b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) 382 383 # Check values at quadrature points 384 with out_vec.array_read() as out_array: 385 for i in range(Q): 386 assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL 387 388# ------------------------------------------------------------------------------- 389