1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3# reserved. See files LICENSE and NOTICE for details. 4# 5# This file is part of CEED, a collection of benchmarks, miniapps, software 6# libraries and APIs for efficient high-order finite element and spectral 7# element discretizations for exascale applications. For more information and 8# source code availability see http://github.com/ceed. 9# 10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11# a collaborative effort of two U.S. Department of Energy organizations (Office 12# of Science and the National Nuclear Security Administration) responsible for 13# the planning and preparation of a capable exascale ecosystem, including 14# software, applications, hardware, advanced system engineering and early 15# testbed platforms, in support of the nation's exascale computing imperative. 16 17from _ceed_cffi import ffi, lib 18import tempfile 19import numpy as np 20from abc import ABC 21from .ceed_constants import TRANSPOSE, NOTRANSPOSE 22 23# ------------------------------------------------------------------------------ 24 25 26class Basis(ABC): 27 """Ceed Basis: finite element basis objects.""" 28 29 # Attributes 30 _ceed = None 31 _pointer = ffi.NULL 32 33 # Representation 34 def __repr__(self): 35 return "<CeedBasis instance at " + hex(id(self)) + ">" 36 37 # String conversion for print() to stdout 38 def __str__(self): 39 """View a Basis via print().""" 40 41 # libCEED call 42 with tempfile.NamedTemporaryFile() as key_file: 43 with open(key_file.name, 'r+') as stream_file: 44 stream = ffi.cast("FILE *", stream_file) 45 46 err_code = lib.CeedBasisView(self._pointer[0], stream) 47 self._ceed._check_error(err_code) 48 49 stream_file.seek(0) 50 out_string = stream_file.read() 51 52 return out_string 53 54 # Apply Basis 55 def apply(self, nelem, emode, u, v, tmode=NOTRANSPOSE): 56 """Apply basis evaluation from nodes to quadrature points or vice versa. 57 58 Args: 59 nelem: the number of elements to apply the basis evaluation to; 60 the backend will specify the ordering in a 61 BlockedElemRestriction 62 emode: basis evaluation mode 63 u: input vector 64 v: output vector 65 **tmode: CEED_NOTRANSPOSE to evaluate from nodes to quadrature 66 points, CEED_TRANSPOSE to apply the transpose, mapping 67 from quadrature points to nodes; default CEED_NOTRANSPOSE""" 68 69 # libCEED call 70 err_code = lib.CeedBasisApply(self._pointer[0], nelem, tmode, emode, 71 u._pointer[0], v._pointer[0]) 72 self._ceed._check_error(err_code) 73 74 # Transpose a Basis 75 @property 76 def T(self): 77 """Transpose a Basis.""" 78 79 return TransposeBasis(self) 80 81 # Transpose a Basis 82 @property 83 def transpose(self): 84 """Transpose a Basis.""" 85 86 return TransposeBasis(self) 87 88 # Get number of nodes 89 def get_num_nodes(self): 90 """Get total number of nodes (in dim dimensions) of a Basis. 91 92 Returns: 93 num_nodes: total number of nodes""" 94 95 # Setup argument 96 p_pointer = ffi.new("CeedInt *") 97 98 # libCEED call 99 err_code = lib.CeedBasisGetNumNodes(self._pointer[0], p_pointer) 100 self._ceed._check_error(err_code) 101 102 return p_pointer[0] 103 104 # Get number of quadrature points 105 def get_num_quadrature_points(self): 106 """Get total number of quadrature points (in dim dimensions) of a Basis. 107 108 Returns: 109 num_qpts: total number of quadrature points""" 110 111 # Setup argument 112 q_pointer = ffi.new("CeedInt *") 113 114 # libCEED call 115 err_code = lib.CeedBasisGetNumQuadraturePoints( 116 self._pointer[0], q_pointer) 117 self._ceed._check_error(err_code) 118 119 return q_pointer[0] 120 121 # Gauss quadrature 122 @staticmethod 123 def gauss_quadrature(q): 124 """Construct a Gauss-Legendre quadrature. 125 126 Args: 127 Q: number of quadrature points (integrates polynomials of 128 degree 2*Q-1 exactly) 129 130 Returns: 131 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 132 and array of length Q to hold the weights""" 133 134 # Setup arguments 135 qref1d = np.empty(q, dtype="float64") 136 qweight1d = np.empty(q, dtype="float64") 137 138 qref1d_pointer = ffi.new("CeedScalar *") 139 qref1d_pointer = ffi.cast( 140 "CeedScalar *", 141 qref1d.__array_interface__['data'][0]) 142 143 qweight1d_pointer = ffi.new("CeedScalar *") 144 qweight1d_pointer = ffi.cast( 145 "CeedScalar *", 146 qweight1d.__array_interface__['data'][0]) 147 148 # libCEED call 149 err_code = lib.CeedGaussQuadrature(q, qref1d_pointer, qweight1d_pointer) 150 self._ceed._check_error(err_code) 151 152 return qref1d, qweight1d 153 154 # Lobatto quadrature 155 @staticmethod 156 def lobatto_quadrature(q): 157 """Construct a Gauss-Legendre-Lobatto quadrature. 158 159 Args: 160 q: number of quadrature points (integrates polynomials of 161 degree 2*Q-3 exactly) 162 163 Returns: 164 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 165 and array of length Q to hold the weights""" 166 167 # Setup arguments 168 qref1d = np.empty(q, dtype="float64") 169 qref1d_pointer = ffi.new("CeedScalar *") 170 qref1d_pointer = ffi.cast( 171 "CeedScalar *", 172 qref1d.__array_interface__['data'][0]) 173 174 qweight1d = np.empty(q, dtype="float64") 175 qweight1d_pointer = ffi.new("CeedScalar *") 176 qweight1d_pointer = ffi.cast( 177 "CeedScalar *", 178 qweight1d.__array_interface__['data'][0]) 179 180 # libCEED call 181 err_code = lib.CeedLobattoQuadrature( 182 q, qref1d_pointer, qweight1d_pointer) 183 self._ceed._check_error(err_code) 184 185 return qref1d, qweight1d 186 187 # QR factorization 188 @staticmethod 189 def qr_factorization(ceed, mat, tau, m, n): 190 """Return QR Factorization of a matrix. 191 192 Args: 193 ceed: Ceed context currently in use 194 *mat: Numpy array holding the row-major matrix to be factorized in place 195 *tau: Numpy array to hold the vector of lengt m of scaling factors 196 m: number of rows 197 n: numbef of columns""" 198 199 # Setup arguments 200 mat_pointer = ffi.new("CeedScalar *") 201 mat_pointer = ffi.cast( 202 "CeedScalar *", 203 mat.__array_interface__['data'][0]) 204 205 tau_pointer = ffi.new("CeedScalar *") 206 tau_pointer = ffi.cast( 207 "CeedScalar *", 208 tau.__array_interface__['data'][0]) 209 210 # libCEED call 211 lib.CeedQRFactorization( 212 ceed._pointer[0], mat_pointer, tau_pointer, m, n) 213 214 return mat, tau 215 216 # Symmetric Schur decomposition 217 @staticmethod 218 def symmetric_schur_decomposition(ceed, mat, n): 219 """Return symmetric Schur decomposition of a symmetric matrix 220 via symmetric QR factorization. 221 222 Args: 223 ceed: Ceed context currently in use 224 *mat: Numpy array holding the row-major matrix to be factorized in place 225 n: number of rows/columns 226 227 Returns: 228 lbda: Numpy array of length n holding eigenvalues""" 229 230 # Setup arguments 231 mat_pointer = ffi.new("CeedScalar *") 232 mat_pointer = ffi.cast( 233 "CeedScalar *", 234 mat.__array_interface__['data'][0]) 235 236 lbda = np.empty(n, dtype="float64") 237 l_pointer = ffi.new("CeedScalar *") 238 l_pointer = ffi.cast( 239 "CeedScalar *", 240 lbda.__array_interface__['data'][0]) 241 242 # libCEED call 243 lib.CeedSymmetricSchurDecomposition( 244 ceed._pointer[0], mat_pointer, l_pointer, n) 245 246 return lbda 247 248 # Simultaneous Diagonalization 249 @staticmethod 250 def simultaneous_diagonalization(ceed, matA, matB, n): 251 """Return Simultaneous Diagonalization of two matrices. 252 253 Args: 254 ceed: Ceed context currently in use 255 *matA: Numpy array holding the row-major matrix to be factorized with 256 eigenvalues 257 *matB: Numpy array holding the row-major matrix to be factorized to identity 258 n: number of rows/columns 259 260 Returns: 261 (x, lbda): Numpy array holding the row-major orthogonal matrix and 262 Numpy array holding the vector of length n of generalized 263 eigenvalues""" 264 265 # Setup arguments 266 matA_pointer = ffi.new("CeedScalar *") 267 matA_pointer = ffi.cast( 268 "CeedScalar *", 269 matA.__array_interface__['data'][0]) 270 271 matB_pointer = ffi.new("CeedScalar *") 272 matB_pointer = ffi.cast( 273 "CeedScalar *", 274 matB.__array_interface__['data'][0]) 275 276 lbda = np.empty(n, dtype="float64") 277 l_pointer = ffi.new("CeedScalar *") 278 l_pointer = ffi.cast( 279 "CeedScalar *", 280 lbda.__array_interface__['data'][0]) 281 282 x = np.empty(n * n, dtype="float64") 283 x_pointer = ffi.new("CeedScalar *") 284 x_pointer = ffi.cast("CeedScalar *", x.__array_interface__['data'][0]) 285 286 # libCEED call 287 lib.CeedSimultaneousDiagonalization(ceed._pointer[0], matA_pointer, matB_pointer, 288 x_pointer, l_pointer, n) 289 290 return x, lbda 291 292 # Destructor 293 def __del__(self): 294 # libCEED call 295 err_code = lib.CeedBasisDestroy(self._pointer) 296 self._ceed._check_error(err_code) 297 298# ------------------------------------------------------------------------------ 299 300 301class BasisTensorH1(Basis): 302 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 303 H^1 discretizations.""" 304 305 # Constructor 306 def __init__(self, ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, 307 qref1d, qweight1d): 308 309 # Setup arguments 310 self._pointer = ffi.new("CeedBasis *") 311 312 self._ceed = ceed 313 314 interp1d_pointer = ffi.new("CeedScalar *") 315 interp1d_pointer = ffi.cast( 316 "CeedScalar *", 317 interp1d.__array_interface__['data'][0]) 318 319 grad1d_pointer = ffi.new("CeedScalar *") 320 grad1d_pointer = ffi.cast( 321 "CeedScalar *", 322 grad1d.__array_interface__['data'][0]) 323 324 qref1d_pointer = ffi.new("CeedScalar *") 325 qref1d_pointer = ffi.cast( 326 "CeedScalar *", 327 qref1d.__array_interface__['data'][0]) 328 329 qweight1d_pointer = ffi.new("CeedScalar *") 330 qweight1d_pointer = ffi.cast( 331 "CeedScalar *", 332 qweight1d.__array_interface__['data'][0]) 333 334 # libCEED call 335 err_code = lib.CeedBasisCreateTensorH1(self._ceed._pointer[0], dim, ncomp, 336 P1d, Q1d, interp1d_pointer, 337 grad1d_pointer, qref1d_pointer, 338 qweight1d_pointer, self._pointer) 339 self._ceed._check_error(err_code) 340 341 # 342 def get_interp1d(self): 343 """Return 1D interpolation matrix of a tensor product Basis. 344 345 Returns: 346 *array: Numpy array""" 347 348 # Compute the length of the array 349 nnodes_pointer = ffi.new("CeedInt *") 350 lib.CeedBasisGetNumNodes1D(self._pointer[0], nnodes_pointer) 351 nqpts_pointer = ffi.new("CeedInt *") 352 lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer) 353 length = nnodes_pointer[0] * nqpts_pointer[0] 354 355 # Setup the pointer's pointer 356 array_pointer = ffi.new("CeedScalar **") 357 358 # libCEED call 359 lib.CeedBasisGetInterp1D(self._pointer[0], array_pointer) 360 361 # Return array created from buffer 362 # Create buffer object from returned pointer 363 buff = ffi.buffer( 364 array_pointer[0], 365 ffi.sizeof("CeedScalar") * 366 length) 367 # return read only Numpy array 368 ret = np.frombuffer(buff, dtype="float64") 369 ret.flags['WRITEABLE'] = False 370 return ret 371 372 373# ------------------------------------------------------------------------------ 374 375 376class BasisTensorH1Lagrange(BasisTensorH1): 377 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange basis 378 objects for H^1 discretizations.""" 379 380 # Constructor 381 def __init__(self, ceed, dim, ncomp, P, Q, qmode): 382 383 # Setup arguments 384 self._pointer = ffi.new("CeedBasis *") 385 386 self._ceed = ceed 387 388 # libCEED call 389 err_code = lib.CeedBasisCreateTensorH1Lagrange(self._ceed._pointer[0], dim, 390 ncomp, P, Q, qmode, self._pointer) 391 self._ceed._check_error(err_code) 392 393# ------------------------------------------------------------------------------ 394 395 396class BasisH1(Basis): 397 """Ceed H1 Basis: finite element non tensor-product basis for H^1 discretizations.""" 398 399 # Constructor 400 def __init__(self, ceed, topo, ncomp, nnodes, 401 nqpts, interp, grad, qref, qweight): 402 403 # Setup arguments 404 self._pointer = ffi.new("CeedBasis *") 405 406 self._ceed = ceed 407 408 interp_pointer = ffi.new("CeedScalar *") 409 interp_pointer = ffi.cast( 410 "CeedScalar *", 411 interp.__array_interface__['data'][0]) 412 413 grad_pointer = ffi.new("CeedScalar *") 414 grad_pointer = ffi.cast( 415 "CeedScalar *", 416 grad.__array_interface__['data'][0]) 417 418 qref_pointer = ffi.new("CeedScalar *") 419 qref_pointer = ffi.cast( 420 "CeedScalar *", 421 qref.__array_interface__['data'][0]) 422 423 qweight_pointer = ffi.new("CeedScalar *") 424 qweight_pointer = ffi.cast( 425 "CeedScalar *", 426 qweight.__array_interface__['data'][0]) 427 428 # libCEED call 429 err_code = lib.CeedBasisCreateH1(self._ceed._pointer[0], topo, ncomp, 430 nnodes, nqpts, interp_pointer, 431 grad_pointer, qref_pointer, 432 qweight_pointer, self._pointer) 433 434# ------------------------------------------------------------------------------ 435 436 437class TransposeBasis(): 438 """Transpose Ceed Basis: transpose of finite element tensor-product basis objects.""" 439 440 # Attributes 441 _basis = None 442 443 # Constructor 444 def __init__(self, basis): 445 446 # Reference basis 447 self._basis = basis 448 449 # Representation 450 def __repr__(self): 451 return "<Transpose CeedBasis instance at " + hex(id(self)) + ">" 452 453 # Apply Transpose Basis 454 def apply(self, nelem, emode, u, v): 455 """Apply basis evaluation from quadrature points to nodes. 456 457 Args: 458 nelem: the number of elements to apply the basis evaluation to; 459 the backend will specify the ordering in a 460 Blocked ElemRestriction 461 **emode: basis evaluation mode 462 u: input vector 463 v: output vector""" 464 465 # libCEED call 466 self._basis.apply(nelem, emode, u, v, tmode=TRANSPOSE) 467 468# ------------------------------------------------------------------------------ 469