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 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 QR factorization 64# ------------------------------------------------------------------------------- 65 66 67def test_301(ceed_resource): 68 ceed = libceed.Ceed(ceed_resource) 69 70 m = 4 71 n = 3 72 a = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], 73 dtype=ceed.scalar_type()) 74 qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], 75 dtype=ceed.scalar_type()) 76 tau = np.empty(3, dtype=ceed.scalar_type()) 77 78 qr, tau = ceed.qr_factorization(qr, tau, m, n) 79 np_qr, np_tau = np.linalg.qr(a.reshape(m, n), mode="raw") 80 81 for i in range(n): 82 assert abs(tau[i] - np_tau[i]) < TOL 83 84 qr = qr.reshape(m, n) 85 for i in range(m): 86 for j in range(n): 87 assert abs(qr[i, j] - np_qr[j, i]) < TOL 88 89# ------------------------------------------------------------------------------- 90# Test Symmetric Schur Decomposition 91# ------------------------------------------------------------------------------- 92 93 94def test_304(ceed_resource): 95 ceed = libceed.Ceed(ceed_resource) 96 97 A = np.array([0.2, 0.0745355993, -0.0745355993, 0.0333333333, 98 0.0745355993, 1., 0.1666666667, -0.0745355993, 99 -0.0745355993, 0.1666666667, 1., 0.0745355993, 100 0.0333333333, -0.0745355993, 0.0745355993, 0.2], 101 dtype=ceed.scalar_type()) 102 103 Q = A.copy() 104 105 lam = ceed.symmetric_schur_decomposition(Q, 4) 106 Q = Q.reshape(4, 4) 107 lam = np.diag(lam) 108 109 Q_lambda_Qt = np.matmul(np.matmul(Q, lam), Q.T) 110 assert np.max(Q_lambda_Qt - A.reshape(4, 4)) < TOL 111 112# ------------------------------------------------------------------------------- 113# Test Simultaneous Diagonalization 114# ------------------------------------------------------------------------------- 115 116 117def test_305(ceed_resource): 118 ceed = libceed.Ceed(ceed_resource) 119 120 M = np.array([0.2, 0.0745355993, -0.0745355993, 0.0333333333, 121 0.0745355993, 1., 0.1666666667, -0.0745355993, 122 -0.0745355993, 0.1666666667, 1., 0.0745355993, 123 0.0333333333, -0.0745355993, 0.0745355993, 0.2], 124 dtype=ceed.scalar_type()) 125 K = np.array([3.0333333333, -3.4148928136, 0.4982261470, -0.1166666667, 126 -3.4148928136, 5.8333333333, -2.9166666667, 0.4982261470, 127 0.4982261470, -2.9166666667, 5.8333333333, -3.4148928136, 128 -0.1166666667, 0.4982261470, -3.4148928136, 3.0333333333], 129 dtype=ceed.scalar_type()) 130 131 X, lam = ceed.simultaneous_diagonalization(K, M, 4) 132 X = X.reshape(4, 4) 133 np.diag(lam) 134 135 M = M.reshape(4, 4) 136 K = K.reshape(4, 4) 137 138 Xt_M_X = np.matmul(X.T, np.matmul(M, X)) 139 assert np.max(Xt_M_X - np.identity(4)) < TOL 140 141 Xt_K_X = np.matmul(X.T, np.matmul(K, X)) 142 assert np.max(Xt_K_X - lam) < TOL 143 144# ------------------------------------------------------------------------------- 145# Test GetNumNodes and GetNumQuadraturePoints for basis 146# ------------------------------------------------------------------------------- 147 148 149def test_306(ceed_resource): 150 ceed = libceed.Ceed(ceed_resource) 151 152 b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO) 153 154 p = b.get_num_nodes() 155 q = b.get_num_quadrature_points() 156 157 assert p == 64 158 assert q == 125 159 160# ------------------------------------------------------------------------------- 161# Test interpolation in multiple dimensions 162# ------------------------------------------------------------------------------- 163 164 165def test_313(ceed_resource): 166 ceed = libceed.Ceed(ceed_resource) 167 168 for dim in range(1, 4): 169 Q = 10 170 Qdim = Q**dim 171 Xdim = 2**dim 172 x = np.empty(Xdim * dim, dtype=ceed.scalar_type()) 173 uq = np.empty(Qdim, dtype=ceed.scalar_type()) 174 175 for d in range(dim): 176 for i in range(Xdim): 177 x[d * Xdim + i] = 1 if (i % 178 (2**(dim - d))) // (2**(dim - d - 1)) else -1 179 180 X = ceed.Vector(Xdim * dim) 181 X.set_array(x, cmode=libceed.USE_POINTER) 182 Xq = ceed.Vector(Qdim * dim) 183 Xq.set_value(0) 184 U = ceed.Vector(Qdim) 185 U.set_value(0) 186 Uq = ceed.Vector(Qdim) 187 188 bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO) 189 bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO) 190 191 bxl.apply(1, libceed.EVAL_INTERP, X, Xq) 192 193 with Xq.array_read() as xq: 194 for i in range(Qdim): 195 xx = np.empty(dim, dtype=ceed.scalar_type()) 196 for d in range(dim): 197 xx[d] = xq[d * Qdim + i] 198 uq[i] = eval(dim, xx) 199 200 Uq.set_array(uq, cmode=libceed.USE_POINTER) 201 202 # This operation is the identity because the quadrature is collocated 203 bul.T.apply(1, libceed.EVAL_INTERP, Uq, U) 204 205 bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS) 206 bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS) 207 208 bxg.apply(1, libceed.EVAL_INTERP, X, Xq) 209 bug.apply(1, libceed.EVAL_INTERP, U, Uq) 210 211 with Xq.array_read() as xq, Uq.array_read() as u: 212 for i in range(Qdim): 213 xx = np.empty(dim, dtype=ceed.scalar_type()) 214 for d in range(dim): 215 xx[d] = xq[d * Qdim + i] 216 fx = eval(dim, xx) 217 assert math.fabs(u[i] - fx) < 1E-4 218 219# ------------------------------------------------------------------------------- 220# Test grad in multiple dimensions 221# ------------------------------------------------------------------------------- 222 223 224def test_314(ceed_resource): 225 ceed = libceed.Ceed(ceed_resource) 226 227 for dim in range(1, 4): 228 P, Q = 8, 10 229 Pdim = P**dim 230 Qdim = Q**dim 231 Xdim = 2**dim 232 sum1 = sum2 = 0 233 x = np.empty(Xdim * dim, dtype=ceed.scalar_type()) 234 u = np.empty(Pdim, dtype=ceed.scalar_type()) 235 236 for d in range(dim): 237 for i in range(Xdim): 238 x[d * Xdim + i] = 1 if (i % 239 (2**(dim - d))) // (2**(dim - d - 1)) else -1 240 241 X = ceed.Vector(Xdim * dim) 242 X.set_array(x, cmode=libceed.USE_POINTER) 243 Xq = ceed.Vector(Pdim * dim) 244 Xq.set_value(0) 245 U = ceed.Vector(Pdim) 246 Uq = ceed.Vector(Qdim * dim) 247 Uq.set_value(0) 248 Ones = ceed.Vector(Qdim * dim) 249 Ones.set_value(1) 250 Gtposeones = ceed.Vector(Pdim) 251 Gtposeones.set_value(0) 252 253 # Get function values at quadrature points 254 bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO) 255 bxl.apply(1, libceed.EVAL_INTERP, X, Xq) 256 257 with Xq.array_read() as xq: 258 for i in range(Pdim): 259 xx = np.empty(dim, dtype=ceed.scalar_type()) 260 for d in range(dim): 261 xx[d] = xq[d * Pdim + i] 262 u[i] = eval(dim, xx) 263 264 U.set_array(u, cmode=libceed.USE_POINTER) 265 266 # Calculate G u at quadrature points, G' * 1 at dofs 267 bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS) 268 bug.apply(1, libceed.EVAL_GRAD, U, Uq) 269 bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones) 270 271 # Check if 1' * G * u = u' * (G' * 1) 272 with Gtposeones.array_read() as gtposeones, Uq.array_read() as uq: 273 for i in range(Pdim): 274 sum1 += gtposeones[i] * u[i] 275 for i in range(dim * Qdim): 276 sum2 += uq[i] 277 278 assert math.fabs(sum1 - sum2) < 10. * TOL 279 280# ------------------------------------------------------------------------------- 281# Test creation and destruction of a 2D Simplex non-tensor H1 basis 282# ------------------------------------------------------------------------------- 283 284 285def test_320(ceed_resource): 286 ceed = libceed.Ceed(ceed_resource) 287 288 P, Q, dim = 6, 4, 2 289 290 in_array = np.empty(P, dtype=ceed.scalar_type()) 291 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 292 qweight = np.empty(Q, dtype=ceed.scalar_type()) 293 294 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 295 libceed.lib.CEED_SCALAR_TYPE]) 296 297 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 298 299 print(b) 300 del b 301 302# ------------------------------------------------------------------------------- 303# Test integration with a 2D Simplex non-tensor H1 basis 304# ------------------------------------------------------------------------------- 305 306 307def test_322(ceed_resource): 308 ceed = libceed.Ceed(ceed_resource) 309 310 P, Q, dim = 6, 4, 2 311 312 xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 313 0., 0.5, 0.5, 1.], dtype=ceed.scalar_type()) 314 in_array = np.empty(P, dtype=ceed.scalar_type()) 315 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 316 qweight = np.empty(Q, dtype=ceed.scalar_type()) 317 318 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 319 libceed.lib.CEED_SCALAR_TYPE]) 320 321 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 322 323 # Interpolate function to quadrature points 324 for i in range(P): 325 in_array[i] = feval(xr[0 * P + i], xr[1 * P + i]) 326 327 in_vec = ceed.Vector(P) 328 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 329 out_vec = ceed.Vector(Q) 330 out_vec.set_value(0) 331 weights_vec = ceed.Vector(Q) 332 weights_vec.set_value(0) 333 334 b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec) 335 b.apply(1, libceed.EVAL_WEIGHT, libceed.VECTOR_NONE, weights_vec) 336 337 # Check values at quadrature points 338 with out_vec.array_read() as out_array, weights_vec.array_read() as weights_array: 339 sum = 0 340 for i in range(Q): 341 sum += out_array[i] * weights_array[i] 342 assert math.fabs(sum - 17. / 24.) < 10. * TOL 343 344# ------------------------------------------------------------------------------- 345# Test grad with a 2D Simplex non-tensor H1 basis 346# ------------------------------------------------------------------------------- 347 348 349def test_323(ceed_resource): 350 ceed = libceed.Ceed(ceed_resource) 351 352 P, Q, dim = 6, 4, 2 353 354 xq = np.array([0.2, 0.6, 1. / 3., 0.2, 0.2, 0.2, 355 1. / 3., 0.6], dtype=ceed.scalar_type()) 356 xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 357 0., 0.5, 0.5, 1.], dtype=ceed.scalar_type()) 358 in_array = np.empty(P, dtype=ceed.scalar_type()) 359 qref = np.empty(dim * Q, dtype=ceed.scalar_type()) 360 qweight = np.empty(Q, dtype=ceed.scalar_type()) 361 362 interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[ 363 libceed.lib.CEED_SCALAR_TYPE]) 364 365 b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight) 366 367 # Interpolate function to quadrature points 368 for i in range(P): 369 in_array[i] = feval(xr[0 * P + i], xr[1 * P + i]) 370 371 in_vec = ceed.Vector(P) 372 in_vec.set_array(in_array, cmode=libceed.USE_POINTER) 373 out_vec = ceed.Vector(Q * dim) 374 out_vec.set_value(0) 375 376 b.apply(1, libceed.EVAL_GRAD, in_vec, out_vec) 377 378 # Check values at quadrature points 379 with out_vec.array_read() as out_array: 380 for i in range(Q): 381 value = dfeval(xq[0 * Q + i], xq[1 * Q + i]) 382 assert math.fabs(out_array[0 * Q + i] - value) < 10. * TOL 383 384 value = dfeval(xq[1 * Q + i], xq[0 * Q + i]) 385 assert math.fabs(out_array[1 * Q + i] - value) < 10. * TOL 386 387# ------------------------------------------------------------------------------- 388