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