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