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